Appunti di Meccanica Meccanica Computazionale delle Strutture Strutture I Mauro Borri Brunetto A.A. 2006/07
Indice 1 Modellazione de dei pr problemi me meccanici
1.1 Introduzione . . . . . . . . . . . . . . . . . . . . . 1.2 Formulazione diretta . . . . . . . . . . . . . . . . 1.3 Formulazione variazionale . . . . . . . . . . . . 1.4 Sistemi dinamici . . . . . . . . . . . . . . . . . . 1.5 1.5 Verso erso la disc discre reti tizz zzaz azion ione: e: un prob proble lema ma mode modello llo Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. 1 . 5 . 8 . 12 . 14 . 20
2 I problemi della fisica-matematica
2.1 2.2 2.2 2.3 2.4
Introduzione . . . . . . . . . . . . . . . . . . . . . . . Clas Classi sific ficaz azio ione ne dell dellee equa equazio zioni ni del del seco second ndo o ordi ordine ne Procedimenti variaz iazion ionali . . . . . . . . . . . . . . . Regol egolee del del calc alcolo olo varia ariazi zion onaale . . . . . . . . . . . .
22
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
3 Soluzione ap a pprossimata di d i pr p roblemi al a l co c ontorno
3.1 Operatori differenziali . . . . . . . . . . . . . . . . . . . . . . 3.2 Concetti preliminari . . . . . . . . . . . . . . . . . . . . . . . 3.3 Metodi dei residu idui pesati . . . . . . . . . . . . . . . . . . . . . 3.3.1 Metodo di colloca ocazio zione . . . . . . . . . . . . . . . . . 3.3.2 .3.2 Collo olloccazio azione ne a sotto otto--dom domini ini . . . . . . . . . . . . . . 3.3.3 Metodo di Galerkin . . . . . . . . . . . . . . . . . . . . 3.3.4 Formulazioni deboli . . . . . . . . . . . . . . . . . . . 3.3.5 Form ormulazioni oni inverse . . . . . . . . . . . . . . . . . . . 3.4 Il metod metodo o degli degli elemen elementi ti finiti: finiti: il punto punto di vista vista matem matemati atico co 3.4.1 .4.1 Equiva uivale len nza dell dellee form formul ulaazion zionii . . . . . . . . . . . . . 3.4.2 Form ormulazione one di Galer lerkin . . . . . . . . . . . . . . . . 3.4.3 .4.3 Il metodo todo deg degli elem elemen entti finit finitii . . . . . . . . . . . . . 4 Il metodo degli elementi finiti
22 23 29 31 37
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
37 39 41 42 45 46 49 53 56 56 60 63 67
4.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 4.2 Disc Discre reti tizz zzaz azio ione ne delle delle equa equazio zioni ni di equi equilib libri rio o . . . . . . . . . . . . 67 II
Indice 1 Modellazione de dei pr problemi me meccanici
1.1 Introduzione . . . . . . . . . . . . . . . . . . . . . 1.2 Formulazione diretta . . . . . . . . . . . . . . . . 1.3 Formulazione variazionale . . . . . . . . . . . . 1.4 Sistemi dinamici . . . . . . . . . . . . . . . . . . 1.5 1.5 Verso erso la disc discre reti tizz zzaz azion ione: e: un prob proble lema ma mode modello llo Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. 1 . 5 . 8 . 12 . 14 . 20
2 I problemi della fisica-matematica
2.1 2.2 2.2 2.3 2.4
Introduzione . . . . . . . . . . . . . . . . . . . . . . . Clas Classi sific ficaz azio ione ne dell dellee equa equazio zioni ni del del seco second ndo o ordi ordine ne Procedimenti variaz iazion ionali . . . . . . . . . . . . . . . Regol egolee del del calc alcolo olo varia ariazi zion onaale . . . . . . . . . . . .
22
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
3 Soluzione ap a pprossimata di d i pr p roblemi al a l co c ontorno
3.1 Operatori differenziali . . . . . . . . . . . . . . . . . . . . . . 3.2 Concetti preliminari . . . . . . . . . . . . . . . . . . . . . . . 3.3 Metodi dei residu idui pesati . . . . . . . . . . . . . . . . . . . . . 3.3.1 Metodo di colloca ocazio zione . . . . . . . . . . . . . . . . . 3.3.2 .3.2 Collo olloccazio azione ne a sotto otto--dom domini ini . . . . . . . . . . . . . . 3.3.3 Metodo di Galerkin . . . . . . . . . . . . . . . . . . . . 3.3.4 Formulazioni deboli . . . . . . . . . . . . . . . . . . . 3.3.5 Form ormulazioni oni inverse . . . . . . . . . . . . . . . . . . . 3.4 Il metod metodo o degli degli elemen elementi ti finiti: finiti: il punto punto di vista vista matem matemati atico co 3.4.1 .4.1 Equiva uivale len nza dell dellee form formul ulaazion zionii . . . . . . . . . . . . . 3.4.2 Form ormulazione one di Galer lerkin . . . . . . . . . . . . . . . . 3.4.3 .4.3 Il metodo todo deg degli elem elemen entti finit finitii . . . . . . . . . . . . . 4 Il metodo degli elementi finiti
22 23 29 31 37
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
37 39 41 42 45 46 49 53 56 56 60 63 67
4.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 4.2 Disc Discre reti tizz zzaz azio ione ne delle delle equa equazio zioni ni di equi equilib libri rio o . . . . . . . . . . . . 67 II
III
INDICE
4.3 4.4 4.4 4.5 4.5 4.6 4.6
Equilibrio dinamico . . . . . . . . . . . . . . . Il meto metodo do dell dellee coor coordi dina nate te gene genera rali lizz zzat atee . . . Caratte atteri risstich tichee dell dellaa solu soluzi zion onee . . . . . . . . . Tratt attamen amento to loc locale ale deg degli elem elemen entti . . . . . . . 4.6. 4.6.11 Numer umeraz azio ione ne loca locale le dei dei grad gradii di liber libertà tà 4.6.2 .6.2 Siste istem mi di rif riferim erimen ento to loc locali ali . . . . . . . 4.7 4.7 Impos mposiz izio ione ne dell dellee cond condiz izio ioni ni al cont contor orno no . . . 4.8 4.8 Rela elazion zionii costi ostitu tuti tive ve elast lastic ich he . . . . . . . . . . 4.8.1 .8.1 Stato tato di tens ension ione monoa onoasssial sialee . . . . . 4.8.2 Stato di tensione one pian iana . . . . . . . . . . 4.8.3 .8.3 Stato tato di defo deform rmaz azio ione ne pian pianaa . . . . . . 4.8.4 .8.4 Cond ondizio izioni ni di assia ssials lsim imm metria tria . . . . . 4.8.5 Lastre inflesse . . . . . . . . . . . . . . . Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
5 Formulazione degli elementi finiti
98
5.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Elementi isoparametrici . . . . . . . . . . . . . . . . . . . . . . . 5.3 Elementi Elementi del continuo continuo n -dimensionale - dimensionale . . . . . . . . . . . . . . 5.3.1 5.3.1 Ma Matri trici ci H e B . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 .3.2 Integ ntegrrazio azion ne numer umeric icaa . . . . . . . . . . . . . . . . . . . . 5.4 Imple Implemen mentaz tazion ionee degli degli eleme elementi nti isopar isoparame ametri trici ci . . . . . . . . .
. . . . . .
. . . . . .
A Supporti matematici matematici
A.1 Spazi vettoriali . . . . . . A.1.1 Definizione . . . A.1.2 Sottospazi . . . . A.1.3 Prodotto interno A.1.4 Norme e metrica metrica
74 75 79 83 83 85 90 92 92 93 93 95 96 96 98 98 101 104 106 107 113
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
113 113 114 115 117
IV
INDICE
Capitolo 1 La modellazione matematica di problemi meccanici 1.1 Introduzione La Meccanica computazionale delle strutture si occupa della formulazione e della soluzione di problemi matematici che descrivono il comportamento meccanico di modelli dei sistemi fisici con i quali l’ingegnere si deve confrontare nella sua attività professionale. Sotto la dicitura generale di comportamento meccanico, si intende la risposta, in termini di tensioni e deformazioni, che una struttura fornisce a fronte delle azioni alla quale viene sottoposta. La struttura in esame è quindi il sistema fisico oggetto di indagine. È bene premettere subito che, dato un certo sistema fisico, magari convenientemente isolato dal resto dell’universo, non esiste in generale un solo modello matematico atto a descriverne il comportamento. Talvolta è difficile costruirne uno appena convincente, talaltra ne esistono diversi, che privilegiano aspetti differenti del problema. L’analisi fisica dovrebbe portare, sulla base dei principi fondamentali della meccanica, alla scelta di un modello convincente del sistema: tale modello conterrà certe relazioni tra grandezze note (i dati ) e grandezze da determinare (le incognite ), dipendenti dai parametri che descrivono alcune caratteristiche del sistema. Una volta costruito il modello, l’analisi si sposta su un piano totalmente matematico: da questo punto di vista si deve stabilire se il problema è posto in modo da non presentare contraddizioni, e se ammette una o più soluzioni. Stabilita l’esistenza di una soluzione si può poi cercare di determinarla esattamente (in forma chiusa ) o, quando ciò non sia possibile o conveniente, di ricercarne un’approssimazione . 1
2
1 MODELLAZIONE DEI PROBLEMI MECCANICI
Allo studio dei metodi di approssimazione per la soluzione di problemi discendenti dalla modellazione matematica dei sistemi fisici è dedicata la Meccanica Computazionale. Una prima classificazione dei modelli di sistemi meccanici li suddivide a seconda se l’incognita sia un campo , cioè una funzione che associa a ogni punto del dominio occupato dalla struttura un valore scalare (ad es. la pressione in un fluido) o vettoriale (ad es. la velocità o lo spostamento), oppure l’incognita sia un elenco finito di valori di una certa grandezza, calcolati in un certo numero di punti del dominio. Il primo tipo di modello è detto continuo , il secondo discreto . Nella Figura 1.1, a titolo di esempio, sono schematizzati due modelli, rispettivamente continuo e discreto, dello stesso sistema meccanico: la trave caricata assialmente. La colonna a sezione variabile raffigurata in Figura 1.1(a), è un esempio di schematizzazione di un sistema fisico: l’analista strutturale ha individuato tale sistema come rappresentativo di unasituazione reale, ne ha individuate la caratteristiche geometriche, le forze agenti e i vincoli. Stabilito poi di procedere a un’analisi elastica del sistema, tendente a determinare lo spostamento di alcuni punti della colonna, sono presentati due possibili modelli. Il primo (Figura 1.1(b)) è un modello continuo che generalizza il problema di Saint-Venant al caso di solidi a sezione variabile, i cui parametri sono le costanti elastiche E , ν, la funzione che descrive l’andamento della sezione A (z ), l’altezza h . I dati del problema sono la forza concentrata applicata all’estremo superiore e il peso proprio; l’incognita è la funzione u (z ) che esprime lo spostamento assiale dei punti della colonna. Il secondo modello è, invece, di tipo discreto: l’incognita è la coppia di valori {u 1 , u 2 } che rappresentano gli spostamenti di due punti del modello. Si tratta di un modello a parametri concentrati, dove i termini di rigidezza k 1 e k 2 sono scelti in modo da riprodurre la deformabilità della struttura di partenza. È evidente che, mentre la soluzione del modello continuo richiede la solu-
Figura 1.1: Un sistema meccanico reale (a); il suo modello continuo (b) e un modello discreto (c).
1.1 INTRODUZIONE
3
zione di un’equazione differenziale con opportune condizioni al contorno (in breve: un problema a valori al contorno ), quella del modello discreto sarà ottenuta attraverso un sistema di equazioni algebriche, i cui coefficienti saranno funzione dei parametri del sistema. Se si eccettuano i casi eccezionalmente semplici, la soluzione dei problemi continui risulta impossibile. D’altro canto, tranne che in situazioni particolari, non è ovvia la determinazione dei parametri concentrati da utilizzare in un modello discreto. Come vedremo, tutti i metodi della Meccanica computazionale sono rivolti alla costruzione di procedimenti coerenti di discretizzazione , cioè di tecniche che consentono di trasformare un modello continuo in un modello discreto equivalente, la cui soluzione possa essere resa il più possibile prossima a quella continua (Figura 1.2).
Figura 1.2: Procedimenti di soluzione della modellazione meccanica: modelli discreti (a sinistra); modelli continui (a destra); discretizzazione di modelli continui
Esempio 1.1
Si consideri una trave di materiale lineare elastico caratterizzato dal modulo di Young E , avente lunghezza L , con area della sezione retta variabile linearmente con la legge A (x ) = A 0 (1 + x /L ). Supponendo che l’unica sollecitazione agente sia lo sforzo normale costante N = F , dove F è una forza assegnata, si valuti l’allungamento della trave. Seguendo il procedimento suggerito dalla Figura 1.1, si può pensare alla struttura come se fosse formata da una serie di n aste di lunghezza L i = L /n , ciascuna caratterizzata da un’opportuna sezione A i , costante lungo l’asta i ,
4
1 MODELLAZIONE DEI PROBLEMI MECCANICI
la cui rigidezza sia pertanto k i = E A i /L i . Il collegamento in serie delle aste fornisce la rigidezza complessiva k secondo la relazione: 1 k
n
=
1.
i =1 k i
Per la valutazione dell’area A i di ciascun’asta, una scelta possibile (anche se L ¯i = n non l’unica) è calcolarla nel punto medio: x (i − 1/2), per cui
¯i ) = A 0 1 + A i = A (x
1 n
1 i −
2
.
L’allungamento cercato risulterà infine: F F L n 1 U = = . k E A 0 i =1 n + i − 12
1 La sommatoria α = n i =1 n +i − 1 si può quindi interpretare come un termine che, 2 al crescere del numero di suddivisioni n , tende al valore corretto per una variazione continua della sezione. Nella Figura 1.3 è riportato l’andamento di tale coefficiente. La soluzione esatta, che si può determinare risolvendo un problema continuo, è log2 = 0,693147 (l’asintoto orizzontale in figura); il valore di α calcolato per n = 100 è pari a 0,693144.
0.69
0.685
0.68
0.675
0.67
n
10
20
30
40
Figura 1.3: Rigidezza asintotica della trave a sezione variabile (Esempio 1.1).
1.2 FORMULAZIONE DIRETTA
5
1.2 Formulazione diretta Per iniziare la trattazione delle tecniche di modellazione matematica dei problemi di meccanica computazionale è conveniente affrontare in primo luogo la classe dei modelli discreti. Gran parte dei concetti meccanici qui utilizzati sono ben noti, e quindi ci si può concentrare sugli aspetti inerenti più strettamente la procedura di soluzione di problemi anche complessi, in vista della loro successiva applicazione all’interno delle tecniche di discretizzazione dei modelli continui che verranno affrontate in seguito. In un modello discreto di un problema fisico: – lo stato del sistema può essere descritto direttamente con precisione adeguata dai valori assunti da un numero finito di variabili di stato ; – le variabili di stato possono essere grandezze fisiche scalari (es: temperatura, pressione) o vettoriali (es: velocità, spostamento); – la soluzione numerica consiste in un elenco dei valori assunti dalle variabili di stato nella situazione definita dai dati del problema. Supponiamo di dover approntare un modello discreto atto a rappresentare la risposta meccanica di un dato sistema. I passi attraverso i quali si snoda l’analisi sono, in linea di massima, i seguenti: Schematizzazione il sistema meccanico viene rappresentato da un insieme di
elementi. Equilibrio si esprime l’equilibrio di ogni elemento in termini delle variabili di
stato. Assemblaggio si tiene conto delle mutue connessioni fra elementi per scrivere
un sistema di equazioni aventi come incognite le variabili di stato, i cui termini noti siano le azioni imposte. Calcolo della risposta si risolve il sistema di equazioni determinando i valori
delle variabili di stato. Nell’esempio successivo la procedura qui indicata viene dettagliatamente sviluppata in tutti i suoi passaggi. Esempio 1.2
L’analisi di un sistema meccanico ha portato alla schematizzazione illustrata nella Figura 1.4(a). Gli elementi elastici (molle) presentano un comportamento
6
1 MODELLAZIONE DEI PROBLEMI MECCANICI
lineare descritto nella Figura 1.4(b), caratterizzato dalla rigidezza k e , e = 1,...5. Le molle sono connesse a tre carrelli: elementi verticali rigidi, appoggiati al suolo mediante incastri scorrevoli. Fa eccezione la molla 1, che ha un estremo fisso. La configurazione del sistema è completamente individuata dalla terna (vettore) {U 1 , U 2 , U 3 }T che raccoglie gli spostamenti dei carrelli, considerati positivi verso destra. Ciascuno dei carrelli è soggetto a una forza, anch’essa positiva verso destra: tali forze si raccolgono nel vettore {R 1 , R 2 , R 3 }T . Indicando con F i (e ) la forza agente nella molla e in corrispondenza al carrello i , per l’equilibrio di ciascuna molla si avrà: k 1U 1 = F 1(1) ,
k k k 2
− 2
4
−k 4
(2) 1 1 = (2) 2 2 (4) 1 1 = U 3 F 3(4)
U F k , k U k F F k k U
−k 2
2
− 4
k 4
3
− 3
,
5
−k 5
(3) 1 = (3) 2 2 (5) 2 2 = U 3 F 3(5)
U F , k U F F k U
−k 3
3
− 5
k 5
1
.
L’assemblaggio si ottiene espandendo le equazioni matriciali di equilibrio alla dimensione propria del vettore delle incognite (in questo caso 3), aggiungendo zeri nei posti opportuni delle varie matrici e imponendo che ciascuno dei carrelli risulti in equilibrio con le forze applicate (notare che le forze F i (e ) sono considerate come applicate alle molle, quindi l’azione sui carrelli va cambiata di segno), perciò: F 1(1) + F 1(2) + F 1(3) + F 1(4) = R 1 ,
(1.1a)
F 2(2) + F 2(3) + F 2(5) = R 2 ,
(1.1b)
F 3(4) + F 3(5) = R 3 .
(1.1c)
Figura 1.4: Sistema discreto a tre gradi di libertà (a); legge di comportamento di una molla (b) (Esempio 1.2).
7
1.2 FORMULAZIONE DIRETTA
La relazione valida per la molla 1 si può scrivere quindi nella forma matriciale: k 1 0 0 U 1 F 1(1) 0 0 0 U 2 = 0 . 0 0 0 U 3 0
Analogamente si procede per le altre molle:
k k 0 k 0 2
− 2
4
−k 4
k 0 U k 0 U 0 0 U 0 k U 0 0 U
− 2
1
2
2
=
3
− 4
0
k 4
k k 0 U k k 0 U 0 0 0 U 0 0 0 U 0 k k U
=
− 3
1
3
2
=
3 1
5
F 3(4)
U 3
− 3
3
(4) 1
1 2
(2) 1 (2) 2
F F , F 0 0 ,
0
−k 5
− 5
k 5
2
=
U 3
(3) 1 (3) 2
F F , 00 F . (5) 2 F 3(5)
Sostituendo queste relazioni nell’eq.(1.1), si ottiene l’equazione di equilibrio globale della struttura.
k
1 + k 2 + k 3 + k 4 −(k 2 + k 3 ) −k 4
−(k 2 + k 3 )
−k 4
k 2 + k 3 + k 5 −k 5 k 4 + k 5 −k 5
U R U R , 1 2
1
=
U 3
2
R 3
ovvero, in forma matriciale, KU = R,
dove si è posto UT = [U 1 U 2 U 3 ], RT = [ R 1 R 2 R 3], e la matrice di rigidezza globale risulta:
K =
k
1 + k 2 + k 3 + k 4 −(k 2 + k 3 ) −k 4
−(k 2 + k 3 )
−k 4
−k 5 k 2 + k 3 + k 5 k 4 + k 5 −k 5
.
La soluzione, in termini del vettore spostamento, è data formalmente dall’espressione: U = K −1 R . Si noti che risulta det K = k 1 [k 4 k 5 + k 2 (k 4 + k 5 ) + k 3 (k 4 + k 5 )]. L’eventuale annullarsi del determinante, come è noto, comporta l’inesistenza della inversa K −1 e la perdita dell’unicità della soluzione: in tale caso l’equilibrio è possibile solo per determinate configurazioni di carico (v. l’Esercizio 1.2).
8
1 MODELLAZIONE DEI PROBLEMI MECCANICI
L’approccio alla soluzione qui illustrato è detto metodo degli spostamenti , in quanto l’incognita primaria del problema è lo spostamento di alcuni punti della struttura. Naturalmente, dopo aver determinato gli spostamenti, è immediato ricavare le sollecitazioni negli elementi strutturali. Esempio 1.3
Riprendendo la soluzione del sistema discusso nell’Esempio 1.2, supponendo ora di aver determinato il vettore degli spostamenti incogniti U = {U 1 U 2 U 3 }T , si può osservare che, ad esempio, per la molla 2, la sollecitazione è data dallo sforzo N 2 = −F 1(2) = k 2U 1 − k 2U 2
e, procedendo allo stesso modo per le altre molle si possono determinare tutti gli sforzi. Uno dei campi di applicazione più comuni del metodo degli spostamenti è la soluzione dei telai piani .1 Si tratta di strutture contenute e caricate in un piano, composte da travi generalmente ad asse rettilineo, connesse fra loro e al suolo mediante incastri o cerniere. In questo caso fra le incognite del problema figurano, oltre che gli spostamenti propriamente detti, anche le rotazioni dei nodi di incastro interno. Il procedimento di soluzione è del tutto simile a quanto descritto nell’esempio precedente, con l’avvertenza di considerare correttamente l’orientamento delle singole travi nella scrittura dell’equilibrio delle forze ai nodi.
1.3 Formulazione variazionale Come è noto, l’esistenza di un potenziale elastico comporta, per un sistema meccanico, la possibilità di ricercare la soluzione come configurazione associata a un minimo dell’energia potenziale totale del sistema. Più in generale, si dimostra che è possibile scrivere un’espressione, detta funzionale , che associa ad ogni configurazione del sistema un numero reale, il quale deve risultare stazionario (cioè presentare un minimo, un massimo o un punto di sella) in corrispondenza alla soluzione del problema. In base a questa proprietà è possibile, come vedremo, ricavare le equazioni di equilibrio del modello anche in casi in cui risulti difficile individuarle direttamente. In particolare, per un dato sistema, assegnato uno spazio vettoriale V (per definizioni e proprietà si veda l’Appendice A) i cui elementi (vettori) individui1 Si
veda ad es. la trattazione contenuta in A. Carpinteri, Scienza delle Costruzioni, Vol. 2, Pitagora Editrice, Bologna, 2a ed., 1993, cap. 14.
9
1.3 FORMULAZIONE VARIAZIONALE
no le possibili configurazioni del sistema medesimo, un funzionale è un’applicazione lineare che associa a ciascuna configurazione un numero reale Π
: V → R .
(1.2)
Nei casi discreti qui presentati, lo spazio V è a dimensione finita , cioè una configurazione qualsiasi è individuata da una ennupla di valori (le variabili di stato); più in generale, nei casi continui che vedremo in seguito, gli elementi di tale spazio possono essere funzioni definite su un certo dominio, ma gran parte dei concetti qui anticipati continueranno a valere. In ogni caso, le equazioni di equilibrio di un dato sistema si ottengono scrivendo le condizioni di estremo per il funzionale. Supponendo che la configurazione del sistema sia individuata da n grandezze U i , i = 1,..., n , con il linguaggio del calcolo delle variazioni si dirà che tali condizioni di estremo si raggiungono qualora risulti nulla la variazione prima del funzionale, cioè qualora: δΠ = 0 .
(1.3)
Con il simbolo δ Π, si intende un nuovo funzionale, che si ottiene quando alle variabili indipendenti U i vengano imposte certe variazioni arbitrarie δU i . Le variazioni arbitrarie δU i hanno come sola condizione quella di rispettare i vincoli, cioè devono essere nulle nel caso in cui la corrispondente variabile U i abbia un valore imposto. La relazione che lega la variazione del funzionale Π alle variazioni arbitrarie δU i è data da δΠ =
∂Π ∂Π δU 1 + · · · + δU n , ∂U 1 ∂U n
(1.4)
per cui la condizione di estremo (1.3), espressa in funzione delle variazioni arbitrarie delle variabili indipendenti si può scrivere ∂Π ∂Π δU 1 + · · · + δU n = 0 . ∂U 1 ∂U n
(1.5)
Il concetto fondamentale da ricordare è chel’eq.(1.4) deve valere comunque si scelgano le variazioni δU i , sempre nel rispetto dei vincoli. In base a questa idea, prendendo ad es. δU 1 = 1 e δU i = 0, i = 1 si ricava necessariamente ∂Π =0. ∂U 1
(1.6)
D’altra parte il ragionamento si può riproporre per qualsiasi delle U i , ottenendo le n equazioni ∂Π = 0, i = 1, n , ∂U i
(1.7)
10
1 MODELLAZIONE DEI PROBLEMI MECCANICI
che sono precisamente le equazioni di equilibrio cercate. Gli esempi seguenti applicano e mettono a fuoco questi concetti con riferimento a casi a uno o più gradi di libertà. Esempio 1.4
Consideriamo il sistema a un grado di libertà illustrato nella Figura 1.5 (a) sup ponendo di non conoscere l’equazione di equilibrio della molla soggetta all’allungamento U , ma anzi di voler determinarla a partire da una formulazione variazionale. Quando all’estremo della molla viene imposto lo spostamento U , il carico compie il lavoro V = PU , mentre la molla accumula l’energia di deformazione 2 U = kU /2; i grafici di queste funzioni sono riportati in Figura 1.5 (b). Per definizione, il potenziale totale è la differenza Π = U − V
.
In funzione di U , quindi (Figura 1.5 (c)): Π=
1 2 kU − PU . 2
Imponendo la condizione di stazionarietà (1.4), si ha: ∂Π δU = 0 ∂U
(1.8)
(kU − P )δU = 0 .
(1.9)
δΠ =
e quindi Poiché δU è arbitraria, affinché sia soddisfatta la condizione di estremo in ogni caso (anche per δU = 0) deve risultare kU − P = 0
che è l’equazione di equilibrio cercata. Si noti poi (v. Figura 1.5 (c)) che Π raggiunge il minimo in corrispondenza della configurazione equilibrata, cioè per U ∗ = P /k .
Esempio 1.5
Determiniamo Π per il sistema dell’Esempio 1.2, derivando le equazioni di equilibrio per via variazionale.
11
1.3 FORMULAZIONE VARIAZIONALE
Figura 1.5: Sistema a un grado di libertà (a); potenziale elastico U , potenziale dei carichi V (b); potenziale totale Π (c). (Esempio 1.4). Il procedimento è del tutto analogo a quello del caso monodimensionale dell’esempio precedente, con la sola differenza data dal carattere vettoriale delle variabili. Il potenziale dei carichi si può scrivere infatti T V = U R .
L’energia di deformazione è data da U =
1 T U KU 2
e quindi il potenziale totale, Π = U − V risulta Π=
1 T U KU − UT R . 2
La condizione di stazionarietà in forma vettoriale si può quindi scrivere: ∂ Π δΠ = δUT =0 ∂U
e le equazioni di equilibrio, stante l’arbitrarietà del vettore δU si possono ottenere esplicitando le condizioni ∂Π =0. ∂U
12
1 MODELLAZIONE DEI PROBLEMI MECCANICI
È molto interessante notare che l’energia elastica si può esprimere come somma dei contributi relativi a ciascuna molla, nel caso in esame U = U 1 + U 2 + U 3 + U 4 + U 5 , dove 1 2 U 1 = k 1U 1 2 1 2 U 2 = k 2 (U 2 − U 1 ) 2 1 2 U 3 = k 3 (U 2 − U 1 ) 2 1 2 U 4 = k 4 (U 3 − U 1 ) 2 1 2 U 5 = k 5 (U 3 − U 2 ) 2 Esplicitando l’espressione di Π in funzione delle U i , è immediato ottenere le equazioni di equilibrio (∂Π/∂U i ) = 0, i = 1,...5: 1 1 1 1 k 1U 12 + k 2 (U 2 − U 1 )2 + k 3 (U 2 − U 1 )2 + k 4 (U 3 − U 1 )2 + 2 2 2 2 1 k 5 (U 3 − U 2 )2 − R 1U 1 − R 2U 2 − R 3U 3 2 Ad esempio, la prima equazione è data da Π=
k 1U 1 − k 2 (U 2 − U 1 ) − k 3 (U 2 − U 1 ) − k 4 (U 3 − U 1 ) − R 1 = 0 ,
ovvero (k 1 + k 2 + k 3 + k 4 )U 1 − (k 2 + k 3 )U 2 − k 4U 3 − R 1 = 0 . Procedendo analogamente per le altre due variabili si perviene all’equazione matriciale di equilibrio del sistema KU = R.
1.4 Sistemi dinamici Nel contesto dei sistemi discreti, meritano un cenno i problemi dinamici. Essi si differenziano dai problemi statici per la presenza delle forze d’inerzia , legate, come è noto, all’accelerazione delle masse tramite la seconda legge di NEWTON: F = m a. In tutti quei casi in cui la presenza delle forze d’inerzia risulti rilevante (carichi variabili velocemente nel tempo, urti, eccitazione sismica ecc.) è necessario considerare fra le azioni anche il loro contributo. Seguendo il principio di D ’ALEMBERT, occorre considerare, oltre ai carichi statici , anche i termini inerziali proporzionali alla massa, come evidenziato nell’esempio seguente.
13
1.4 SISTEMI DINAMICI
Esempio 1.6
Consideriamo il sistema dell’Esempio 1.2 in condizioni dinamiche: i carrelli sono dotati di massa m 1 , m 2 , m 3 e i carichi sono applicati secondo una legge temporale assegnata (Figura 1.6).
Figura 1.6: Sistema dinamico a tre gradi di libertà (a); andamento dei carichi R i in funzione del tempo (b). (Esempio 1.6). Le equazioni di equilbrio (1.1) diventano, in questo caso equazioni differenziali ¨ 1 , F 1(1) + F 1(2) + F 1(3) + F 1(4) = R 1 (t ) − m 1U
(1.10a)
¨2 , F 2(2) + F 2(3) + F 2(5) = R 2 (t ) − m 2U
(1.10b)
¨ 3 , F 3(4) + F 3(5) = R 3 (t ) − m 3U
(1.10c)
¨ i si indica la derivata seconda, fatta rispetto al tempo, dello spodove con U stamento del carrello i e con la scrittura R i (t ) si rende esplicita la dipendenza dal tempo dei carichi assegnati. Di conseguenza, anche tutte le grandezze incognite ( U i e F i (e ) ) sono da intendersi funzioni del tempo t . Introducendo la matrice di massa
m 0
1
M=
0
0 m 2
0
0 0 , m 3
le equazioni (1.10) si possono scrivere in forma compatta come ¨ + KU = R(t ) MU ¨ 1 U ¨ 2 U ¨ 3 ]T . ¨ = [U dove si è introdotto anche il vettore delle accelerazioni U
14
1 MODELLAZIONE DEI PROBLEMI MECCANICI
1.5 Verso la discretizzazione: un problema modello Consideriamo, come esempio atto a introdurre l’argomento dei metodi di calcolo più usati nell’ambito della meccanica computazionale, una struttura tanto semplice quanto familiare: la trave ad asse rettilineo caricata assialmente. Obiettivo dei problemi che incontreremo nel seguito sarà la determinazione del cosiddetto campo di spostamento , in generale a valori vettoriali, ma in questo caso particolare, nel quale sono possibili solo spostamenti nella direzione dell’asse della trave, O x , la soluzione deve portare alla determinazione della funzione scalare u (x ). Supponiamo quindi di voler determinare il campo di spostamento per una trave di lunghezza unitaria, soggetta a un carico distribuito p (x ), al cui estremo sinistro (x = 0) è applicata la forza F , mentre all’estremo destro ( x = 1) è imposto lo spostamento U (Figura 1.7). Detto N (x ) lo sforzo normale agente nella sezione posta all’ascissa x , e indicata con p (x ) la densità di forza per unità di lunghezza, l’equazione indefinita d’equilibrio è data da dN (1.11) + p (x ) = 0 . dx Dalla teoria tecnica della trave, che generalizza i risultati relativi al solido di de Saint-Venant, ricordiamo che, detta A la sezione trasversale della trave, lo stato tensionale è completamente definito dalla tensione normale σxx (x ) =
N (x ) A
(1.12)
alla quale corrisponde la deformazione longitudinale xx (x ) =
σxx (x ) , E
dove E è il modulo di Young.
Figura 1.7: Schema del problema monodimensionale.
(1.13)
1.5 VERSO LA DISCRETIZZAZIONE: UN PROBLEMA MODELLO
15
La deformazione longitudinale, inoltre, può essere legata allo spostamento u tramite la definizione du . (1.14) xx (x ) = dx Esprimendo dapprima la tensione σ xx in funzione dello spostamento, mediante le equazioni (1.13) e (1.14), sostituendola poi nella (1.12), si ottiene: du . dx
N (x ) = E A
(1.15)
Sostituendo questa espressione nell’equazione di equilibrio (1.11) si ricava un’equazione differenziale per l’incognita u : d2 u p (x ) + =0. dx 2 E A
(1.16)
L’equazione differenziale (1.16) stabilisce la relazione che deve sussistere tra la funzione nota p (x ), cioè il carico applicato, e la funzione incognita u (x ). Questa relazione coinvolge la derivata seconda di u , quindi è un’equazione del secondo ordine. Si noti che è necessario che la funzione u (x ), ancora incognita, risulti derivabile con continuità almeno due volte, in modo che tale derivata esista; è altresì necessario che la funzione p (x ), che costituisce insieme ai parametri E , A i dati di questo problema, sia continua . Sono dunque escluse dal campo di validità dell’equazione (1.16) situazioni con dati irregolari , quali, ad esempio, quelle costituite da distribuzioni di carico p (x ) che presentino salti (discontinuità di prima specie finite).2 Osserviamo ora le condizioni agli estremi. Mentre all’estremo destro è richiesto allo spostamento di assumere un dato valore, cioè u (1) = U , la condizione imposta all’estremo sinistro è un po’ più complicata. Infatti, si richiede che la forza applicata sulla base del solido sia pari a un valore assegnato F .3 Occorre dunque esprimere tale condizione in termini legati all’incognita u (x ). Ciò non è difficile ricordando ancora l’equazione (1.15) e notando che N (0) = −F , quindi si ottiene: F du . (1.17) =− E A dx x =0
A questo punto possiamo cercare la soluzione imponendo che essa rispetti anche le condizioni agli estremi destro e sinistro: la funzione incognita u (x ) 2 Analogamente, una versione generalizzata del problema che ammettesse modulo di Young
variabile non potrebbe comprendere, per lo stesso motivo, ad es. situazioni con due materiali aventi moduli E 1 ed E 2 diversi fra loro. 3 Vedremo in seguito che le condizioni imposte allo spostamento si dicono condizioni essenziali , mentre quelle imposte alla sua derivata si chiamano condizioni naturali .
16
1 MODELLAZIONE DEI PROBLEMI MECCANICI
dovrà quindi soddisfare contemporaneamente l’equazione di campo e le condizioni al contorno : d2 u p (x ) + = 0, 0 < x < 1 , dx 2 E A F du (1.18) , =− E A dx x =0 u |x =1 = U .
Il procedimento di soluzione passa attraverso l’integrazione dell’equazione differenziale, che comporta la comparsa delle costanti d’integrazione . Queste ultime dovranno essere determinate in modo da rispettare tutte le condizioni al contorno imposte alla soluzione. Nel semplice caso in esame è possibile ottenere il risultato in forma chiusa per una vasta classe di forme di carico p (x ) (Esempio 1.7), ma in generale non sempre si può determinare analiticamente la funzione incognita u (x ) , anche se tale soluzione esiste. In questi casi si deve ricorrere a metodi di approssimazione di tipo numerico, che permettono di ottenere soluzioni che, pur non essendo quella esatta, si possono considerare vicine ad essa, a meno di un errore che può essere reso sufficientemente piccolo. Esempio 1.7
Un pilastro di materiale isotropo lineare elastico con modulo di Young E , avente altezza H e sezione costante A è incastrato alla base. Il carico consiste nel peso per unità di volume γ e in una forza verticale di compressione F applicata all’estremo superiore (Figura 1.8). Determinare il campo di spostamento u (x ). Per ricondurre il problema a quello descritto nel testo, occorre esprimere la forza assiale per unità di lunghezza p (x ) in funzione del peso volumico γ . Considerando un concio di lunghezza d x , deve risultare p (x )dx = −γ A dx , da cui p (x ) = −γ A . La condizione di vincolo impone spostamento nullo per x = 0; il carico ap plicato all’estremo superiore impone il valore della deformazione per x = H . Il problema differenziale da risolvere è quindi il seguente: Determinare la funzione u (x ) ∈ C 2 [0, H ] tale che: 2
d u dx du dx u 2
−
γ = 0, E =−
x =H
|x =0 = 0 .
0 < x < H ,
F , E A
1.5 VERSO LA DISCRETIZZAZIONE: UN PROBLEMA MODELLO
17
Figura 1.8: Schema del problema monodimensionale. Integrando due volte l’equazione differenziale si ottiene: u −
γ 2 x + C 1 x + C 2 = 0 . 2E
Dalla condizione essenziale u |x =0 si ottiene immediatamente C 2 = 0. Imponendo la condizione naturale all’estremo superiore, per x = H si ricava C 1 =
γ AH + F . E A
Sostituendo le costanti così determinate, si ottiene la funzione cercata: x F 1 x u (x ) = γH . −1 − E A 2 H
(1.19)
Si noti che, successivamente, è possibile ricavare lo stato di sollecitazione nel pilastro, ricordando l’eq.(1.15), per cui: N (x ) = −F − γ A ( H − x )
e, ad esempio, lo sforzo normale alla base si ricava ponendo x = 0, ottenendo N (0) = −F − γ AH ,
somma dei contributi del carico esterno e del peso proprio del pilastro.
18
1 MODELLAZIONE DEI PROBLEMI MECCANICI
Nel caso in cui non sia possibile la soluzione in forma chiusa del problema differenziale è necessario ricorrere a una procedura di tipo numerico tramite la quale le equazioni vengono discretizzate. A questo fine, si può cercare di approssimare gli operatori di derivata prima e seconda mediante incrementi finiti. Si considerino gli sviluppi in serie di T AYLOR della funzione spostamento u (x ), arrestati al secondo ordine, nell’intorno di un punto x i : u (x i +1 ) = u (x i ) + u (x i −1 ) = u (x i ) +
du 1 (x i +1 − x i ) + dx x =x i 2 du 1 (x i −1 − x i ) + dx x =x i 2
d2 u (x i +1 − x i )2 2 dx x =x i d2 u (x i −1 − x i )2 2 dx x =x i
(1.20) (1.21)
I punti x i −1 = x i − ∆x e x i +1 = x i + ∆x si scelgono entrambi distanti da x i della stessa quantità ∆x , pensata sufficientemente piccola. Sottraendo l’eq. (1.21) dall’eq. (1.20) e sommando le due equazioni si ottengono rispettivamente le espressioni d u u (x i +1 ) − u (x i −1 ) =2 ∆x , dx x =x i d2 u 2 u (x i +1 ) + u (x i −1 ) =2u (x i ) + ∆ x , 2 dx x =x i
(1.22) (1.23)
dalle quali si possono ricavare le espressioni approssimate delle derivate prima e seconda, dette differenze finite : u (x i +1 ) − u (x i −1 ) du , = dx x =x i 2∆x u (x i +1 ) − 2u (x i ) + u (x i −1 ) d2 u . = dx 2 x =x i ∆x 2
(1.24) (1.25)
Mediante le espressioni delle differenze finite, un’equazione differenziale viene discretizzata, scrivendo un insieme di equazioni algebriche la cui soluzione può essere ottenuta in modo da soddisfare anche le condizioni al contorno. Consideriamo a questo scopo il problema (1.18). L’equazione differenziale viene sostituita dalla opportuna espressione alle differenze finite, scritta in un generico punto x i del dominio: u (x i +1 ) − 2u (x i ) + u (x i −1 ) p (x i ) + =0. E A ∆x 2
(1.26)
Scrivendo l’espressione precedente per n + 2 punti, tenendo conto delle due condizioni al contorno, si possono ricavare n valori approssimati dello spostamento cercato, come illustrato nell’esempio seguente.
19
1.5 VERSO LA DISCRETIZZAZIONE: UN PROBLEMA MODELLO
Esempio 1.8
Risolvere con il metodo delle differenze finite il problema dell’Esempio 1.7. Consideriamo un’approssimazionedella soluzione mediante sei puntiequidistanti (con passo ∆x = 1/5) di ascissa compresa tra 0 e 1. A questi aggiungiamo due ulteriori punti esterni all’intervallo rispettivamente a x = −∆x e x = 1 + ∆x , i quali serviranno a imporre le condizioni al contorno. L’insieme di punti viene così a consistere di 8 punti, da x 1 a x 8 . Ricorrendo all’eq. (1.26), si possono scrivere le espressioni seguenti: −u 1 + 2u 2 − u 3 ∆ x 2
−u 2 + 2u 3 − u 4 ∆ x 2
−u 3 + 2u 4 − u 5 ∆ x 2
−u 4 + 2u 5 − u 6 ∆ x 2
−u 5 + 2u 6 − u 7 ∆ x 2
−u 6 + 2u 7 − u 8 ∆ x 2
p , E A p , = E A p , = E A p , = E A p , = E A p . = E A =
La condizione al contorno essenziale, per x = 0, viene imposta molto semplicemente ponendo: u 2 = 0 . La condizione naturale, che impone il valore della forza all’estremo superiore, si esprime mediante l’approssimazione della derivata prima: −u 6 + u 8
=−
F . E A
2∆x Tutte queste equazioni costituiscono un sistema lineare che viene convenientemente espresso nella forma matriciale seguente, dove si è posto ∆x = 1/5:
0 1 00 0 0 0 −
0
1 2 −1 0 0 0 0 0
0 −1 2 −1 0 0 0 0
0 0 −1 2 −1 0 0 0
0 0 0 −1 2 −1 0 0
0 0 0 0 −1 2 −1 −1
0 0 0 0 0 −1 2 0
0 0 0 0 0 0 −1 1
u u u u u u u 1 2 3 4 5 6 7
u 8
=
0 1 1 1 . 1 1 1 −10F /p
p 25E A
(1.27)
20
1 MODELLAZIONE DEI PROBLEMI MECCANICI
Nella Figura 1.9 sono diagrammati i valori degli spostamenti u i , i = 1,8 ottenuti dalla soluzione del sistema ponendo, per semplicità, p /E A = γ/E = −1 e F = 0. Nello stesso diagramma, con la linea continua, è anche riportata la funzione u (x ) ottenuta dalla soluzione analitica ricavata nell’esempio 1.7 (eq. (1.19)), che mostra la coincidenza fra i risultati. Non è difficile comprendere il motivo di una così sorprendente efficacia del metodo numerico nel caso particolare qui esaminato se si osserva che la soluzione esatta è un polinomio di secondo grado, per il quale le espressioni delle derivate (1.24), (1.25) sono esatte, anziché, come nel caso generale, approssimate.
Figura 1.9: Confronto tra la soluzione discreta (punti) e quella in forma chiusa (linea continua) per il problema dell’Esempio 1.8.
Esercizi 1.1. Determinare, per il modello illustrato nell’Esempio 1.2, la soluzione che si ottiene supponendo tutte le molle caratterizzate dalla rigidezza k . [U 1 = −(R 1 + R 2 + R 3 )/k , U 2 = −(5R 1 + 7R 2 + 6R 3 )/5k , U 3 = −(5R 1 + 6R 2 + 8R 3 )/5k ]. 1.2. Analizzare il modello illustrato nell’Esempio 1.2 nel caso in cui venga soppressa la molla 1 (si ponga ad esempio k 1 = 0). È possibile la soluzione? Se
sì, sotto quali condizioni? Discutere le diverse possibilità con riferimento alle caratteristiche della matrice di rigidezza globale. 1.3. Diagrammare l’andamento dello spostamento e della tensione normale in
un pilastro di conglomerato cementizio, sezione di 1000 cm2 e altezza di 10 m,
1.5 VERSO LA DISCRETIZZAZIONE: UN PROBLEMA MODELLO
21
caricato in testa da una forza di compressione pari a 100 kN. Si assuma, per il conglomerato, un modulo di Young pari a 25GPa e un peso di volume di 24kN/m3 . 1.4. Generalizzare il problema illustrato nel testo per studiare il caso di un pilastro a sezione A (x ) variabile con continuità. Ripetere l’esercizio 1.3 supponen-
do che la sezione quadrata abbia il lato variabile linearmente con l’altezza che vada dai 50cm, alla base, ai 20cm in sommità. 1.5. Adottando la formulazione ricavata nell’esercizio 1.4, ricavare la legge di variazione della sezione con l’altezza, A (x ), tale da indurre nella struttura una tensione σ xx costante per tutte le sezioni. Supponendo poi una sezione circolare di raggio r , diagrammare l’andamento r (x ) del raggio nel caso di carico esterno F = 0. 1.6. Risolvere il problema presentato nell’esercizio 1.3 mediante una discretiz-
zazione alle differenze finite, osservando l’influenza del numero di punti impiegati sulla soluzione, confrontando i risultati ottenuti con la soluzione analitica.
Capitolo 2 Formulazione dei problemi della fisica-matematica 2.1 Introduzione In questo capitolo si richiamano i principali concetti relativi alla definizione di alcuni classici problemi matematici, ai quali vengono spesso ricondotti i modelli di numerosi sistemi fisici di grande interesse. Ove possibile si farà riferimento a problemi meccanici, ma si tenga presente che, poiché la generalità dell’approccio è in larga misura indipendente dal particolare sistema fisico oggetto di studio, gran parte delle considerazioni svolte risulteranno applicabili a tutti i problemi trattati. Diversamente dal caso dei modelli discreti a cui si è fatto cenno nel Capitolo 1, per i quali è possibile un’immediata soluzione in termini di ennuple di valori rappresentanti la variabile incognita in alcuni punti, affronteremo in questo capitolo i cosiddetti modelli continui . In questo caso le variabili di stato costituiscono dei campi , cioè variano in genere con continuità nel dominio esaminato. Le variabili di stato possono essere scalari (ad es. la temperatura o la pressione di un fluido), o vettori (ad es. la velocità o lo spostamento). Nella ricerca di una formulazione adeguata di un certo modello, si devono costruire, di solito, equazioni differenziali alle derivate parziali che definiscono il bilancio di certe quantità su elementi infinitesimi, le quali forniscono le condizioni che la soluzione deve soddisfare in ogni punto della regione dello spazio in cui si definisce il problema. Normalmente, le condizioni di bilancio non sono sufficienti a consentire, da sole, la formulazione del problema. È pertanto necessario stabilire relazioni locali (valide in ciascun punto) fra i diversi campi in gioco, tramite le cosiddette relazioni costitutive , il cui classico prototipo è il legame fra tensioni e deformazioni, ad esempio di tipo elastico. 22
2.2 CLASSIFICAZIONE DELLE EQUAZIONI DEL SECONDO ORDINE
23
Come è noto, però, la scrittura delle equazioni differenziali di un certo modello non è sufficiente per ottenere una soluzione. Il problema deve essere completato dalle condizioni al contorno e, nel caso di problemi dipendenti dal tempo, dalle condizioni iniziali .
2.2 Classificazione delle equazioni del secondo ordine Al fine di introdurre una nomenclatura che consenta di delimitare accuratamente l’ambito dei problemi oggetto di studio, è utile procedere dapprima alla loro classificazione. Al di là della semplice distinzione formale, come vedremo, questa classificazione consente di distinguere molto nettamente fra classi di modelli atte a descrivere fenomeni fisici fra loro molto differenti. Una prima sostanziale divisione separa la classe delle equazioni lineari da quellegenericamente dette non lineari . Un’equazione alle derivate parziali nella variabile u si dice lineare se si può scrivere come combinazione lineare della funzione u e delle sue derivate, mediante coefficienti che siano indipendenti da u . Esempio 2.1
∂2 u ∂2 u ∂2 u ∂u x + y +2 + + − u + e =0 2 2 ∂x ∂x ∂ y ∂ y ∂ y
è un’equazione lineare a coefficienti costanti. 2 ∂2 u ∂2 u ∂u 2 ∂ u cos(x y ) 2 + 2x + + − u = 0 ∂x ∂x ∂ y ∂ y 2 ∂ y
è un’equazione lineare a coefficienti variabili. 2 ∂2 u ∂2 u 2 ∂ u + 2x + ∂x 2 ∂x ∂ y ∂ y 2
2
=0
è un’equazione non lineare. Le equazioni lineari godono di proprietà che ne rendono la soluzione più semplice rispetto alle equazioni non lineari. Per esse è spesso possibile determinare la soluzione analitica, in forma chiusa , per lo meno per casi particolari di dominio e di dati. In genere, le equazioni non lineari non possono essere risolte analiticamente e l’approccio numerico è, in tali casi, una scelta obbligata.
2 I PROBLEMI DELLA FISICA-MATEMATICA
24
Un’altra classificazione molto importante consente di distinguere le equazioni differenziali a seconda del tipo, evidenziandone nettamente le particolarità che le rendono adatte a modellare certe classi di sistemi fisici, che possono essere descritti da modelli matematici fra loro molto simili. Consideriamo a questo fine il caso particolare, per quanto molto importante, delle equazioni differenziali del secondo ordine. Ricordiamo che l’ordine di un’equazione differenziale è dato dal massimo grado di derivazione che compare nell’equazione medesima. Per semplicità limitiamo, inoltre, l’analisi al caso di un problema definito in una regione bidimensionale , nella quale sia definito un sistema di riferimento cartesiano Ox y . La più generale equazione differenziale del secondo ordine nel dominio bidimensionale, per un campo descritto dalla funzione u (x , y ) si può scrivere nella forma: ∂2 u ∂2 u ∂2 u ∂u ∂u A (x , y ) 2 + 2B (x , y ) , + C (x , y ) = Ψ x , y , u , ∂x ∂x ∂ y ∂ y 2 ∂x ∂ y
(2.1)
dove A , B , C , Ψ sono funzioni assegnate , cioè note, dei loro argomenti. Si consideri ora la funzione B 2 − AC , detta discriminante , in generale funzione delle coordinate (x , y ) : si danno tre casi, a seconda del segno di tale funzione. In particolare l’equazione differenziale (2.1) si dirà – di tipo ellittico se B 2 − AC < 0; – di tipo parabolico se B 2 − AC = 0; – di tipo iperbolico se B 2 − AC > 0. Ciascun tipo di equazione individua una classe di problemi fra loro affini in quanto al tipo di soluzione. Ciascuna classe può essere rappresentata da un problema prototipo , che descrive ben noti problemi fisici. Ad esempio, i prototipi di equazione per i diversi tipi possono essere i seguenti: – per il tipo ellittico: l’equazione di L AP LACE (es. filtrazione nei mezzi porosi, conduzione del calore in regime stazionario); – per il tipo parabolico: l’equazione della diffusione (es. consolidazione dei mezzi porosi saturi, propagazione del calore in regime non stazionario); – per il tipo ellittico: l’equazione delle onde (es. trasmissione del suono nei mezzi elastici, vibrazioni). L’analisi di alcuni esempi permette di considerare la formulazione di problemi differenziali per i diversi tipi di equazioni lineari del secondo ordine.
2.2 CLASSIFICAZIONE DELLE EQUAZIONI DEL SECONDO ORDINE
25
Esempio 2.2
La filtrazione nei mezzi porosi è governata da un’equazione differenziale di ti po ellittico. In questo esempio, dopo aver ottenuto l’equazione differenziale che governa il moto stazionario dell’acqua nel terreno permeabile, si completa il problema imponendo le opportune condizioni al contorno per la geometria illustrata nella Figura 2.1(a), che schematizza la situazione del moto di filtrazione al di sotto di una diga impermeabile in condizioni piane. La relazione tra flusso e gradiente idraulico è la legge di D ARCY con permeabilità isotropa k . L’equazione di conservazione della massa d’acqua che attraversa, in un dato tempo, un elemento infinitesimo di volume avente lati dx e d y (Figura2.1(b)), si scrive ( q y y − q y y +d y )dx − ( q x x − q x x +dx )d y = 0 , (2.2)
dove q = q x i + q y j è il vettore flusso, ovvero la portata in volume attraverso una superficie di area unitaria. Data la continuità di tutte le funzioni coinvolte, si esprimono i flussi uscenti attraverso i differenziali:
x +dx = q x x + (∂q x /∂x )dx , q y y +d y = q y y + (∂q y /∂ y )d y .
q x
(2.3a) (2.3b)
e si introduce la relazione costitutiva (legge di Darcy) che lega il flusso al potenziale φ: q x = − k (∂φ/∂x ) , q y = − k (∂φ/∂ y ) .
(2.4a) (2.4b)
Sostituendo le eq. (2.3) nelle (2.2) e tenendo conto delle eq. (2.4) si ottiene ∂2 φ ∂2 φ + =0, ∂x 2 ∂ y 2
(2.5)
che è la classica equazione di Laplace, esprimibile anche attraverso l’operatore di Laplace: ∂2 (·) ∂2 (·) 2 , ∇ (·) = + 2 2 ∂x
come
∂ y
2 ∇ φ=0 .
È immediato verificare l’ellitticità dell’equazione di Laplace osservandoche, per confronto con l’espressione canonica (2.1) si ha A = 1, B = 0, C = 1, quindi B 2 − AC = −1. Passando ora alla definizione della condizioni al contorno, il problema ne prevede di due tipi:
2 I PROBLEMI DELLA FISICA-MATEMATICA
26
Figura 2.1: Equazione di tipo ellittico (filtrazione): geometria del problema (a); conservazione della massa (b) (Esempio 2.2). – Essenziali, se definiscono il valore della funzione da determinare (in questo caso il carico totale φ) – Naturali, se definiscono il valore della derivata prima della funzione (il flusso q) Osservando lo schema della Figura 2.1(a), vediamo che si devono imporre condizioni al contorno essenziali all’interfaccia fra acqua e terreno, alla base dei due serbatoi di monte e di valle, ove è prescritto il carico idraulico, pari all’altezza del pelo libero:
φ y =L = h 1 , per − ∞ ≤ x ≤ −b /2 , φ y =L = h 2 , per b /2 ≤ x ≤ ∞ .
Lungo i bordi a contatto con zone impermeabili, cioè alla base della diga e sul piano del substrato roccioso, si devono imporre condizioni al contorno naturali; infatti, su tali linee, la componente del flusso perpendicolare al bordo deve essere nulla. Tale condizione si traduce, mediante la legge di Darcy, in una condizione sulla derivata prima del potenziale φ:
q y y =0 = −k (∂φ/∂ y ) y =0 = 0 ,
per − ∞ ≤ x ≤ ∞ ,
q y y =L = −k (∂φ/∂ y ) y =L = 0 ,
per − b /2 ≤ x ≤ b /2 .
2.2 CLASSIFICAZIONE DELLE EQUAZIONI DEL SECONDO ORDINE
27
Per completare la definizione delle condizioni al contorno, consideriamo il dominio di forma rettangolare allungata nella direzione dell’asse x , con i lati verticali molto lontani dalla zona di interesse. È un’ipotesi ragionevole considerare tali bordi impermeabili, ciò che equivale a consentire l’ingresso e l’uscita d’acqua nel sistema possa avvenire solo attraverso i due serbatoi di monte e di valle:
q x x =−∞ = −k (∂φ/∂x ) q x x =∞ = −k (∂φ/∂x )
x =−∞ = 0,
per 0 ≤ x ≤ L , per 0 ≤ x ≤ L .
x =∞ = 0,
Esempio 2.3
Uno strato di argilla satura di spessore 2 H (Figura 2.2(a)) viene sottoposto im provvisamente a un sovraccarico per unità di superficie pari a p ; ciò comporta una variazione nel tempo della pressione dell’acqua nei pori, descritta tramite la teoria della consolidazione monodimensionale di Terzaghi. L’equazione differenziale che lega la sovrappressione interstiziale u alla posizione z e all’istante t è ∂2 u ∂u c v 2 = . (2.6) ∂z
∂t
Il confronto con l’espressione canonica (2.1) permette di classificare l’eq. (2.6) fra le equazioni paraboliche. Infatti si ha A = c v , B = 0, C = 0, quindi B 2 − AC = 0. In questo caso abbiamo a che fare con un problema di evoluzione nel tem po: oltre alle condizioni da imporre sul bordo del dominio, si devono precisare le condizioni iniziali del campo incognito u . Le condizioni al contorno, agli estremi dello strato, sono di tipo essenziale: la sovrappressione deve essere nulla: u |z =0 = 0 , u |z =2H = 0 .
Le condizioni iniziali ipotizzano che, per tutti i punti appartenenti allo strato, all’istante iniziale t = 0 la sovrappressione assuma il valore del sovraccarico p : u (z , 0) = u 0 .
Come è noto, introducendo il tempo adimensionale T = c v t /H 2 , la soluzione del problema è m =∞ M z −M 2 T 2u 0 u (z , t ) = sen e , (2.7)
m =0
M
H
28
2 I PROBLEMI DELLA FISICA-MATEMATICA
Figura 2.2: Equazione di tipo parabolico (consolidazione): geometria del problema (a); andamento della sovrappressione nel tempo (b) (Esempio 2.3). dove si è posto M = Π2 (2m + 1). Nella Figura 2.2(b) sono riportate alcune delle isocrone della soluzione, cioè l’andamento della sovrappressione u per dati istanti di tempo. Il processo porta, al crescere del tempo, alla progressiva diminuzione della sovrappressione nello strato di argilla, ma la completa dissipazione richiede, teoricamente, un tempo infinito.
Esempio 2.4
Consideriamo una barra elastica, vincolata a un estremo, avente lunghezza L e sezione A , caratterizzata dal modulo di Young E edallamassavolumica ρ (Figura 2.3). Supponendo che all’estremo libero all’istante t = 0 venga applicato un carico R 0 mantenuto poi costante, è necessario analizzare il comportamento della struttura in campo dinamico. Detto u lo spostamento assiale, se si considera l’equilibrio di un concio infinitesimo, di lunghezza d x occorre considerarne l’inerzia, proporzionale alla massa dm e all’accelerazione, per cui: ∂2 u , − σx |x A + σx |x +dx A = dm ∂t 2
2.3 PROCEDIMENTI VARIAZIONALI
29
Figura 2.3: Equazione di tipo iperbolico (onde elastiche): geometria del problema e andamento del carico nel tempo (Esempio 2.4). ovvero, essendo σx |x +dx = σ x |x + (∂σx /∂x )dx e dm = ρ A dx : ∂σx ∂2 u =ρ 2 , ∂x ∂t
quindi, ricordando che σx = E (∂u /∂x ), si ottiene l’equazione differenziale ∂2 u ρ ∂2 u , = ∂x 2 E ∂t 2
(2.8)
tipica della propagazione delle onde; la grandezza c = E ρ è detta celerità ed esprime la velocità di propagazione delle onde lungo l’assedel solido. In questo caso il confronto con l’espressione canonica (2.1) permette di classificare l’eq. (2.8) fra le equazioni iperboliche. Infatti si ha A = 1, B = 0, C = −ρ /E , quindi B 2 − AC = ρ /E > 0. Le condizioni al contorno sono le seguenti:
– all’estremo vincolato la condizione essenziale u |x =0 = 0 ; – all’estremo caricato la condizione naturale E A (∂u /∂x )|x =L = R 0 . inoltre si devono imporre le condizioni iniziali di quiete nel dominio corrispondenti a spostamenti e velocità nulla in tutti i punti u |t =0 = 0 , per 0 < x < L ,
∂u ∂x
= 0 , per 0 < x < L . t =0
2.3 Procedimenti variazionali In alternativa ai metodi differenziali, nei quali si parte dalla scrittura di equazioni di bilancio valide per elementi infinitesimi del continuo, gli stessi problemi
2 I PROBLEMI DELLA FISICA-MATEMATICA
30
possono essere affrontati, talvolta convenientemente, attraverso formulazioni variazionali. Analogamente a quanto anticipato trattando i modelli discreti nel paragrafo 1.3, per seguire questa strada è necessario scrivere l’espressione di un funzionale Π (il potenziale) e quindi imporne la stazionarietà rispetto alle variabili di stato. Questa condotta risulta talvolta conveniente sia per la ricerca delle equazioni differenziali che governano un dato problema, sia per la definizione delle condizioni al contorno naturali. Come abbiamo già detto, un funzionale è un’applicazione lineare che associa ad ogni elemento di uno spazio vettoriale V un numero reale, cioè Π
: V → R .
(2.9)
Nel caso continuo, gli elementi dello spazio vettoriale V sono funzioni . Ad esempio, Π può essere l’area sottesa dal grafico di una funzione f (x ) continua su un intervallo [a , b ]: b
Π=
f (x )dx .
(2.10)
a
Un esempio familiare di funzionale definito in ambito meccanico è l’energia potenziale totale Π della trave elastica caricata trasversalmente (Figura 2.4): ad ogni deformata elastica (elemento di V ) è associato un numero reale che misura l’energia Π, differenza fra l’energia di deformazione elastica e il lavoro compiuto dal carico a seguito della deformazione. L’applicazione dei concetti variazionali ai modelli continui non è che una generalizzazione di quanto accennato a proposito della formulazione variazionale dei modelli discreti. In questo caso, rimanendo per semplicità nell’ambito delle funzioni di una variabile, la variazione prima di una funzione v (x ) ∈ V è, per definizione, la funzione δv (x ) = η(x ) , (2.11)
Figura 2.4: Il potenziale totale per una trave caricata come esempio di funzionale lineare Π : V → R.
2.4 REGOLE DEL CALCOLO VARIAZIONALE
31
dove ∈ R è un numero arbitrario, mentre η (x ) è una funzione regolare, nulla nei punti ove la funzione v (x ) soddisfa le condizioni essenziali. A titolo di esempio, nella Figura 2.5 sono illustrati i diversi termini per un caso ipotetico. È importante notare che δv è una nuova funzione della variabile x , che possiede le stesse doti di regolarità di v e può perciò essere derivata e integrata al pari di essa.
2.4 Regole del calcolo variazionale È utile riassumere qui alcune proprietà utili nella successiva trattazione di problemi con tecniche variazionali. Dato il carattere introduttivo, non si darà dimostrazione di quanto asserito, rimandando il lettore interessato ai numerosi testi disponibili. Assumiamo che una funzione F , per un certo valore di x dipenda dalla variabile di stato v e delle sue derivate spaziali: dv d2 v dp v F v , , 2 ,..., p . dx dx dx
Si definisce variazione prima di F la funzione ∂F ∂F dv ∂F dp v . + +···+ δF = δ δ ∂v ∂dv /dx dx ∂dp v /dx p dx p
(2.12)
La successione di operatore variazionale δ e operatore differenziale d, rispetto alle variabili indipendenti, si può invertire: dn η dn (δv ) dn v . = =δ dx n dx n dx n
Figura 2.5: Variazione di una funzione di una variabile.
2 I PROBLEMI DELLA FISICA-MATEMATICA
32
L’operatore δ gode delle stesse proprietà della derivata, ad esempio: δ(F + Q ) = δF + δQ , δ(FQ ) = (δF )Q + F (δQ ) , δ(F n ) = n (F n −1 )δF .
Inoltre, l’operatore δ passa sotto il segno di integrale (e viceversa): δ
F (x )dx =
δF (x )dx .
Nell’esempio seguente vedremo l’applicazione di queste tecniche a un problema meccanico, mostrando come, attraverso un procedimento variazionale si possa arrivare alla scrittura delle equazioni differenziali che lo governano e delle condizioni al contorno naturali. Esempio 2.5
Scrivere il potenziale totale per la struttura schematizzata in Figura 2.6 e ricavare, mediante la formulazione variazionale, l’equazione differenziale del problema (cioè l’equazione della linea elastica) e le condizioni al contorno naturali.
Figura 2.6: Formulazione variazionale per una trave (Esempio 2.5). Affrontando il problema in termini di spostamento, scegliamo, in prima ap prossimazione, di trascurare la deformazione assiale della trave, supponendola ininfluente rispetto agli effetti legati alla flessione. La variabile di stato sarà quindi la componente di spostamento trasversale all’asse w (x ). In modo analogo a quanto visto precedentemente, occorre dapprima formulare l’espressione del potenziale totale per il sistema, che in questo caso si può scrivere Π = U T + U M − V P ,
2.4 REGOLE DEL CALCOLO VARIAZIONALE
33
ove sono individuati rispettivamente l’energia di deformazione della trave, U T , l’energia di deformazione della molla, U M , l’energia potenziale del carico ap plicato V P . Calcoliamo ora i tre contributi, in funzione dello spostamento w . Energia di deformazione della trave. Considerando un tronco infinitesimo, avente lunghezza infinitesima d x , della trave: detti rispettivamente M e χ il
momento flettente e la curvatura, il teorema di C LAPEYRON permette di scrivere il corrispondente lavoro di deformazione 1 1 dL = M χdx , 2 quindi, ricordando che momento e curvatura sono legati dalla relazione M = E I χ e, inoltre, che nell’ipotesi di piccoli spostamenti e piccole rotazioni χ = (∂2 w /∂x 2 ), si ha 2
1 d2 w dL = E I dx . 2 dx 2
Infine, integrando lungo tutta la lunghezza della trave l’energia elastica risulterà 2 1 L d2 w E I dx . U T = 2 0 dx 2
Energia di deformazione della molla. Dipende in maniera quadratica dall’al-
lungamento secondo l’espressione (v. l’Esempio 1.4 del Capitolo 1): U M =
1 2 kw , 2 L
dove si è indicato con w L = w |x =L lo spostamento dell’estremo della trave. Energia potenziale del carico. La valutazione di questo termine è più compli-
cato, in quanto è necessario valutare lo spostamento assiale dell’estremo libero della trave u L = u |x =L : V P = Pu L .
Tale spostamento, avendo deciso di trascurare la deformazione assiale della trave, è legato alla rotazione dell’asse. Si consideri infatti un concio 1 Si veda ad es. A. Carpinteri, Scienza delle Costruzioni, Vol. 1, Pitagora Editrice, Bologna, 2 a
ed., 1993, par. 9.3.
2 I PROBLEMI DELLA FISICA-MATEMATICA
34
infinitesimo di lunghezza d x : a seguito di una rotazione φ, la sua proiezione sull’asse x sarà accorciata della quantità d s = (1 − cos φ)dx . Conviene poi introdurre lo sviluppo accorciato in serie di M C L AUR IN della funzione coseno: cos φ = 1 −
φ2
2 + o (φ )
2! per cui l’accorciamento totale dell’asse della trave, e quindi lo spostamento del punto di applicazione del carico P è dato da L
u L =
L φ2
0
ds =
2
0
dx .
Legando infine la piccola rotazione allo spostamento trasversale: φ = (∂w /∂x ), integrando sulla trave e moltiplicando per P si ottiene l’energia cercata: P L dw 2 V P = dx 2 0 dx
L’espressione del funzionale del problema, cioè l’energia totale, è quindi 1 Π= 2
L
d2 w 2 1 2 P E I x kw − d + dx 2 2 L 2
0
L
dw 2 dx . dx
0
Analogamente a quanto visto nel caso discreto, per ottenere le equazioni del problema si ricerca la stazionarietà del funzionale, imponendo l’annullarsi della sua variazione prima, δΠ = 0: 2
d2 w 1 2 P E I x kw − d + dx 2 2 L 2
1 δ 2
L
0
dw 2 dx dx
L
0
=0.
Sfruttando le proprietà dell’operatore δ acuisièfattocenno,sipuòscrivere: L
0
d2 w d2 δw E I 2 dx + kw L (δw )L − P dx dx 2
L dw
dδw dx = 0 . dx dx
0
(2.13)
Attraverso una integrazione per parti è possibile eliminare dall’eq. (2.13) le derivate della funzione δw di ordine più elevato (il secondo). Il primo termine richiede due integrazioni per parti successive: L
0
L
L d2 w d2 δw d2 w dδw dw 3 dδw E I 2 E I 3 dx = E I 2 dx = − dx dx 2 dx dx 0 0 d x dx L L L d2 w dδw dw 3 dw 4 E I 2 E I 4 δw dx , δw + − E I dx dx 0 d3 x d x 0 0
(2.14)
2.4 REGOLE DEL CALCOLO VARIAZIONALE
35
l’ultimo un’integrazione sola: P
L dw
dδw dw dx = P δw dx dx dx
0
L d2 w
L
− P
0
0
dx 2
δ w dx .
(2.15)
Sostituendo le identità (2.14) e (2.15) nell’espressione (2.13), esplicitando gli incrementi indicati e raccogliendo opportunamente, si ottiene: dw 4 d2 w E I 4 + P 2 δw dx d x dx 0 2 d w dδw d2 w dδw + E I − E I dx 2 dx x =L dx 2 dx x =0 dw 3 d2 w dw 3 d2 w w δ δw − E I + P − E I + P 3 x 2 x d3 x dx 2 d d x =L + kw L (δw )L = 0 . L
(2.16) x =0
Ricordiamo ora che le variazioni δ w e δ(dw /dx ) = dδw /dx devono rispettare le condizioni al contorno essenziali, quindi esse sono nulle all’estremo x = 0. Perciò risulta immediatamente: d2 w dδw E I 2 dx dx
e
=0
x =0
dw 3 d2 w E I 3 + P 2 δw =0. d x dx x =0 Sostituendo queste condizioni nell’eq. (2.16) si ricava
L
dw 4 d2 w d2 w dδw E I 4 + P 2 δw dx + E I 2 d x dx dx dx x =L dw 3 d2 w E I 3 + P 2 δw + kw L (δ w )L = 0 . d x dx x =L
0
−
(2.17)
Scegliendo oradelle variazioni particolari, è possibile ricavare sia l’equazione differenziale del problema, sia le condizioni al contorno naturali. Si prenda, a questo fine, una variazione δw nulla, insieme alla sua derivata prima dδw /dx , in x = L , ma diversa da zero negli altri punti del dominio. Con questa posizione si avrà d2 w dδw E I 2 =0, dx dx x =L dw 3 d2 w E I 3 + P 2 δw d x dx kw L (δw )L = 0 .
(2.18) =0,
(2.19)
x =L
(2.20)
2 I PROBLEMI DELLA FISICA-MATEMATICA
36
Sostituendo queste espressioni nell’eq. (2.17) si può vedere come essa si riduca a L dw 4 d2 w E I 4 + P 2 δw dx = 0 (2.21) d x dx 0 la quale, dovendo essere soddisfatta per variazioni δw = 0 impone
dw 4 d2 w E I 4 + P 2 = 0, d x dx
∀x ∈]0, L [ .
(2.22)
L’espressione (2.21) è l’equazione differenziale della linea elastica in presenza di carico assiale P , ed è l’equazione di campo del problema.
Capitolo 3 Soluzione approssimata di problemi al contorno 3.1 Operatori differenziali Indicando con L un generico operatore differenziale di tipo ellittico, la corrispondente equazione differenziale si può scrivere: L (u ) = b ,
in Ω
(3.1)
e si può interpretare nel modo seguente: l’operatore L trasforma una funzione u , definita in una regione Ω dello spazio, avente frontiera Γ, in una funzione b . In generale, la funzione sulla quale agisce l’operatore L può essere anche un vettore u. Esempio 3.1 Alcuni esempi di operatori differenziali, definiti su un intervallo di R (i primi due) o su una regione bidimensionale (il terzo):
d2 (•) , dx 2 d4 (•) d(•) L (•) = , + dx 4 dx ∂2 (•) ∂2 (•) L (•) = 2 + 2 . L (•) =
∂x
∂ y
(3.2) (3.3) (3.4)
Quando si applica l’operatore a una funzione u , essa va sostituita al posto del punto fra parentesi.
37
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
38
Data una qualsiasi funzione w , si può definire il prodotto interno :
〈L (u ), w 〉 =
L (u )w dΩ .
(3.5)
Ω
L’operatore si dice definito positivo se:
〈L (u ), u 〉 =
L (u )u dΩ > 0,
∀u = 0 .
Ω
(3.6)
Si può integrare per parti l’espressione del prodotto interno, ottenendone la forma trasposta:
Ω
L (u )w dΩ =
uL ∗ (w )dΩ +
Ω
∗
S (w )G (u ) G ∗ (w )S (u ) dΓ .
−
Γ
(3.7)
L’operatore L ∗ si dice aggiunto di L . Se L ∗ = L , quest’ultimo si dice autoaggiunto, e in tal caso anche S ∗ = S e G ∗ = G . I valori imposti a S (u ) sulla porzione Γs di frontiera sono le condizioni al contorno essenziali, mentre i valori imposti a G (u ) sulla porzione Γg di frontiera sono le condizioni al contorno naturali. Esempio 3.2 Data l’equazione differenziale omogenea
d2 u 2 − λ u = 0, 0 < x < 1 , dx 2
(3.8)
definire gli operatori L (u ), S (u ) e G (u ), verificando anche che L è autoaggiunto. Occorre scrivere l’espressione del prodotto interno, procedendo quindi alle successive integrazioni per parti:
〈L (u ), w 〉 =
1
0
d2 u 2 − λ u w dx . dx 2
(3.9)
Una prima integrazione per parti dà:
〈L (u ), w 〉 =
1
− 0
du dw 2 du λ uw dx + w − dx dx dx
1
.
(3.10)
0
Mediante una seconda integrazione per parti si ottiene:
〈L (u ), w 〉 =
1
0
d2 w 2 du λ − w u x w d + dx 2 dx
1
dw u dx
1
− 0
0
.
(3.11)
3.2 CONCETTI PRELIMINARI
39
Dal confronto fra l’espressione generica del prodotto interno in forma trasposta e quella qui ottenuta, si possono dedurre gli operatori:
d2 (•) ∗ L (•) = − λ2(•) ,
dx 2 S ∗ (•) = (•) , d(•) G ∗ (•) = , dx
(3.12) (3.13) (3.14)
e, come si può notare, l’operatore è autoaggiunto in quanto L ∗ (•) = L (•).
3.2 Concetti preliminari Si supponga di voler trovare la soluzione u 0 di un problema differenziale L (u 0 ) = b ,
in Ω , S (u 0 ) = s , su Γs , G (u 0 ) = g , su Γg .
(3.15)
Spesso la soluzione esatta, u 0 , è impossibile da determinare, in questo caso, anziché determinare tale soluzione, ci si accontenta di una soluzione u che, in qualche modo, approssimi quella esatta Si può scegliere una soluzione approssimata nella forma n
u = α0 +
αk φk ,
(3.16)
k =1
dove α k sono parametri da determinare e φk sono funzioni delle coordinate spaziali. Le funzioni φk sono tali da soddisfare certe condizioni di ammissibilità relative a condizioni al contorno e grado di continuità e le funzioni φk devono essere linearmente indipendenti, cioè α1 φ1 + α2 φ2 + · · · + αn φn = 0
solo nel caso in cui tutti i coefficienti siano nulli, cioè α1 = α2 = · · · = αn = 0 .
La successione delle funzioni φk deve essere completa , cioè, fissato un numero > 0 arbitrariamente piccolo, e una qualsiasi funzione ammissibile u 0 , è possibile trovare un numero di termini n superato il quale si abbia
Ω
(u − u 0 )2 ≤ .
(3.17)
40
3 SOLUZIONE SOLUZIONE APPROSSIMA APPROSSIMAT TA DI PROBLEMI PROBLEMI AL CONTORNO CONTORNO
La soluzione esatta del problema soddisfa le condizioni (3.15); nel caso, invece, in cui si introduca una soluzione approssimata u , che si suppone soddisfi esattamente tutte le condizioni al contorno, si può introdurre una funzione residuo o errore: R = L (u ) − b . (3.18)
Si dovranno dovranno determina determinare re i coefficient coefficientii incogniti incogniti αk cercand cercando o di rendere renderepicco piccolo, lo, in qualche senso medio sul dominio Ω, l’errore (3.18). Nel caso in cui la funzione u non non soddisfi tutte la condizioni al contorno, si dovranno considerare considerare anche gli errori sulle condizioni al contorno essenziali e naturali: R s s = S (u )
− s , R g g = G (u ) − g .
(3.19) (3.20)
Esempio Esempio 3.3 Si consideri l’equazione differenziale
d2 u − b = 0, 0 < x < 1 dx 2
(3.21)
dove b b = cost. e le condizioni al contorno sono omogenee, cioè u |x =0 = 0, La soluzione esatta è
1 2
u |x =1 = 0,
u 0 = b x (x 1)
−
Si consideri invece un’approssimazione sinusoidale che soddisfi le condizioni al contorno: u = α1 φ1 = α sen(πx ) La funzione residuo assume la forma:
d2u − b = α sen(πx ) R = dx 2
(3.22)
Per determinar determinare e l’unico parametro parametro incognito incognito α ,si può imporr imporre, e, ad esempi esempio, o, R = 0 per x x = 1/2: b 1 b = 0 = α = απ2 π2 2 La soluzione approssimata e il residuo corrispondente sono quindi:
−
−
⇒
− πb sen(πx ) R = b sen(πx ) − b
−
(3.23)
u =
(3.24)
3.3 METODI METODI DEI RESIDUI PESATI
41
Figura 3.1: Grafico soluzione esatta per x = 1/2.
Scegliendo, invece, di azzerare il residuo per x x = 1/4, si trova
b 2 u = sen(πx ) π2 R = b 2sen(πx ) b
−
−
(3.25)
(3.26)
3.3 3.3 Metod etodii dei dei res esid idui ui pe pesa sati ti In generale, scelta un’approssimazione u della della funzione incognita u 0 da determinare, esistono certe funzioni dette errori o residui , che misurano lo scostamento dalla soluzione esatta sia nel dominio, sia sul contorno: R (x ) R s s (x ) R g g (x )
in Ω in Γs in Γg
(3.27) (3.28) (3.29)
I metodi dei residui pesati cercano cercano di distribuire gli errori in modo da ottenere un errore medio nullo. Si assuma che l’approssimazione u soddisfi soddisfi esattamente le condizioni al contorno (R s s = R g g = 0), è ragionevole richiedere che, data una funzione peso
3 SOLUZIONE SOLUZIONE APPROSSIMA APPROSSIMAT TA DI PROBLEMI PROBLEMI AL CONTORNO CONTORNO
42
Figura 3.2: Grafico soluzione esatta per x = 1/4. w , definita in Ω, si abbia:
〈R , w 〉 =
R w dΩ = 0
(3.30)
Ω
I diversi metodi si differenziano nella scelta della funzione approssimante u ; ad esempio: n
u = α0 +
αk φk
(3.31)
k =1
o della funzione peso w ; ad esempio: n
w =
βk ψk
(3.32)
k =1
3.3. 3.3.1 1 Metod etodo o di coll colloc ocaz azio ione ne Si impone che che l’errore l’errore si annulli in un numero finito finito di punti. punti. Si può interinterpretare come metodo dei residui pesati, assumendo come funzione peso una combinazione lineare di funzioni delta di DIRAC 1 La delta di Dirac ∆ è una distribuzione , o funzione generalizzata , definita tramite le sue proprietà: (x ξ) = 0,
∆
−∞ +
−∞
∀x = ξ
(x − ξ)dx = 1
∆
1 P. A. M. Dirac, 1902-1984, premio Nobel per la fisica nel 1933.
(3.33) (3.34)
3.3 METODI DEI RESIDUI PESATI
43
Risulta, inoltre, per qualsiasi funzione continua f (x ): +
∞
f (x )∆(x ξ)dx = f (ξ)
−
−∞
L’impiego di delta di Dirac come funzione peso consente di dare un’importanza “infinita” al residuo nei punti ξi . Assumendo quindi come funzione peso una combinazione lineare del tipo: w = β1 ∆(x ξ1 ) + β2 ∆(x ξ2 ) + · · · + βk ∆(x ξk )
−
−
−
l’applicazione del metodo dei residui pesati porta al metodo di collocazione, cioè a imporre le condizioni in k punti discreti del dominio. Infatti, ponendo uguale a zero l’espressione del residuo pesato si ottiene
R w dΩ = 0 =
⇒ R (ξ1) = 0, R (ξ2) = 0, · · · R (ξk ) = 0
Ω
Esempio 3.4 Risolvere, mediante un’approssimazione polinomiale e il metodo di collocazione, il seguente problema a valori al contorno:
d2 u + u + x = 0, 0 < x < 1 dx 2 u |x =0 = 0 u |x =1 = 0
(3.35) (3.36) (3.37)
Il polinomio che soddisfa le condizioni al contorno è: x (1
− x )(α1 + α2 x + · · · )
Si consideri il caso di due parametri incogniti, α1 e α 2 , corrispondente a un polinomio di terzo grado: u = x (1
− x )(α1 + α2 x )
L’espressione del residuo diventa:
d2 u R = + u + x = dx 2
2 + x − x 2 α1 + 2 − 6x + x 2 − x 3 α2 + x
−
Volendo ad esempio annullare il residuo nei punti x 1 = 1/4 e x 2 = 1/2, la funzione peso di Dirac è:
1 4
1 2
− −
w = β1 ∆ x
+ β2 ∆ x
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
44
Con questa scelta si impone: 1
0
R w dΩ = 0 =
⇒
1 R 4
= 0,
1 2
R
=0
Le condizioni imposte si traducono nel sistema lineare:
29 16 7 4
− 35 64 7 8
1 α1 = 41 α2 2
(3.38)
Le cui soluzioni sono α 1 = 6/31 e α 2 = 40/217 . La soluzione approssimata cercata è quindi: 42 + 40x u = x (1 x )
−
217
L’andamento del residuo è: R =
−4 − 19x − 2x 2 − 40x 3 217
La figura 3.3 riporta la soluzione approssimata u a confronto con quella esatta ( u 0 = x + sen x / sen1 ), insieme con l’errore R
−
Figura 3.3: Metodo di collocazione.
3.3 METODI DEI RESIDUI PESATI
45
3.3.2 Collocazione a sotto-domini Si divide il dominio in sotto-regioni Ωi (eventualmente anche sovrapposte) e poi si impone l’annullamento della media del residuo in ciascuna di esse:
R dΩ = 0,per x
Ωi
∈ Ωi
Esempio 3.5 Si consideri ancora il problema:
d2 u + u + x = 0, 0 < x < 1 dx 2 u |x =0 = 0 u |x =1 = 0
(3.39) (3.40) (3.41)
con un’approssimazione polinomiale x (1
− x )(α1 + α2 x + · · · )
Studiamo dapprima l’approssimazione quadratica (solo α1 ), impiegando, quindi, un solo sotto-dominio coincidente con l’intero dominio Ω: u = x (1
− x )α1
L’espressione del residuo diventa:
d2 u R = + u + x = dx 2
2 + x − x 2 α1 + x
−
Dalla condizione di collocazione a sotto-domini: 1
0
⇒ R = 113
R dx = 0 =
Il risultato di prima approssimazione è dunque: u (1) =
3 x (1 − x ) 11
Cerchiamo ora un’approssimazione più ricca, mediante un polinomio di terzo grado; si ottiene il residuo:
d2 u R = + u + x = dx 2
2 + x − x 2 α1 + 2 − 6x + x 2 − x 3 α2 + x
−
46
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
In questo caso si devono determinare due parametri α 1 e α 2 , perciò sono necessarie due condizioni. Per imporre tali condizioni si considerano due sottodomini, ad es. Ω1 = (0,1/2) e Ω2 = (0,1), imponendo poi 1 2
0
0
1
R dx = 0
(3.42)
R dx = 0
(3.43)
La valutazione di questi integrali porta a un sistema di due equazioni lineari nelle incognite α1 e α2 che in forma matriciale è
− −
11 12 11 6
53 − 192 − 11 12
1 8 1 2
− −
α1 = α2
(3.44)
da cui risulta α1 = 97/517 e α2 = 24/141 . Il risultato di seconda approssimazione è dunque: 291 + 264x u (2) = x (1 x )
−
1151
Nella figura 3.4 sono riportati, a confronto con la soluzione esatta u 0 , i grafici delle due soluzioni e dei rispettivi residui. Si noti che il residuo R (1) ha media nulla su tutto il dominio, mentre R (2) ha media nulla anche su (0,1/2)
3.3.3 Metodo di Galerkin Dato il problema differenziale L (u 0 ) = b ,
in Ω
con condizioni al contorno omogenee, scelta una funzione approssimante
u =
αk φk
k
con il metodo di G AL ER KI N, il residuo R = L (u ) − b viene minimizzato utilizzando come funzione peso una funzione ottenuta come combinazione lineare delle stesse funzioni base:
w =
βk φk
k
La condizione
〈R , w 〉 =
Rw dΩ = 0
Ω
3.3 METODI DEI RESIDUI PESATI
47
Figura 3.4: Metodo di collocazione a sottodomini. fornisce, nel caso in cui L sia lineare, un sistema di equazioni lineari per determinare i coefficienti α k . La funzione peso di Galerkin può essere interpretata come variazione prima della funzione incognita, con βk = δαk : w = δu = δα1 φ1 + δα2 φ2 + · · · + δαn φn
La condizione sul residuo assume la forma:
Ω
Rw dΩ =
Ω
R β1 φ1 dΩ +
Ω
R β2 φ2 dΩ + · · · +
Ω
R βn φn dΩ = 0
e quindi, stante l’arbitrarietà dei coefficienti βk , scegliendo ad esempio β1 = 1 e βk = 0, k = 1 risulta: R φ1 dΩ = 0
Ω
oppure, assumendo β2 = 1 e βk = 0, k = 2:
Ω
R φ2 dΩ = 0
Questeequazioni costituiscono un sistema lineare perdeterminarei coefficienti αk dell’approssimazione cercata
48
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
Esempio 3.6 Risolvere, mediante un’approssimazione polinomiale di terzo grado e il metodo di Galerkin, il seguente problema a valori al contorno:
d2 u + u + x = 0, 0 < x < 1 dx 2 u |x =0 = 0 u |x =1 = 0
(3.45) (3.46) (3.47)
Il polinomio che approssima la soluzione e la funzione peso sono rispettivamente u = x (1 x )(α1 + α2 x ) = α1 x (1 x ) + α2 x 2 (1 x )
− − − w = x (1 − x )(β1 + β2x ) = β1 x (1 − x ) + β2 x 2 (1 − x )
Si scrivono quindi le equazioni: R φ1 dΩ = 0 R φ2 dΩ = 0 Ω
Ω
dove le funzioni base sono:
φ1 = x (1
− x ),
φ2 = x 2 (1
− x )
e il residuo è lo stesso ricavato negli esempi precedenti e risulta:
d2 u + u + x = R = dx 2
2 + x − x 2 α1 + 2 − 6x + x 2 − x 3 α2 + x
−
Dopo le sostituzioni e le integrazioni si ottiene il sistema simmetrico:
− 203 α1 − 121 (3.48) 13 α2 = − 105 − 201 da cui risulta α 1 = 71/369 e α 2 = 7/41. Infine, le espressioni della soluzione
− −
3 10 3 20
approssimata e del residuo sono: u = R =
71 + 63x x (1 − x ) 369
−16 − 62x − 8x 2 − 63x 3 369
Nel grafico della figura 3.5 sono riportate le soluzioni esatta u 0 e approssimata u , con l’andamento dell’errore R .
3.3 METODI DEI RESIDUI PESATI
49
Figura 3.5: Grafici della soluzione esatta e approssimata del problema descritto nell’esempio 3.6.
3.3.4 Formulazioni deboli Si consideri il problema seguente: trovare la funzione u tale che d2 u + u + x = 0, 0 < x < 1 dx 2 u |x =0 = 0 u |x =1 = 0
(3.49) (3.50) (3.51)
Se denotiamo con H 0 la classe delle funzioni u che soddisfano le condizioni al contorno essenziali (omogenee) sopra indicate, dotate di sufficiente regolarità , lo stesso problema si può porre nella forma integrale seguente: Determinare la funzione u ∈ H 0 tale che 1
0
d2 u + u + x w dx = 0, dx 2
∀w ∈ H 0
(3.52)
Naturalmente, perché abbia senso tale formulazione, la funzione incognita u , in quanto al requisito di regolarità, deve avere una derivata seconda almeno continua a tratti, mentre per w ciò non è necessario. Deve però risultare integrabile il loro prodotto. La funzione u può appartenere, in generale, a una classe di funzioni H ∗ diversa da H A partire dalla formulazione precedente è possibile, mediante integrazione per parti, arrivare a un’espressione ove gli operatori su w e u siano dello stesso
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
50
tipo: 1
0
du du − uw − xw dx = 0, dx dx
∀w ∈ H 01
(3.53)
dove H 01 indica la classe delle funzioni derivabili che soddisfano la condizioni al contorno del problema, che siano dotate della sufficiente regolarità: questa è una formulazione debole del problema. Il tema generale della regolarità delle funzioni esula dagli argomenti qui trattati e viene convenientemente trattata ricorrendo alle proprietà degli spazi di S OBOLEV ; per quanto riguarda la comprensione di quanto segue, basti notare che entrambe le funzioni u e w devono essere a derivata prima quadratointegrabile , in modo da garantire l’esistenza dell’integrale quando, per esempio, si assuma (Galerkin): u = α1 φ1 + α2 φ2 + · · · w = β1 φ1 + β2 φ2 + · · ·
(3.54) (3.55) (3.56)
Le funzioni base φk devono quindi avere derivata prima quadrato-integrabile. Per definizione, una funzione f (x ) si dice quadrato-integrabile, e si scrive f (x ) ∈ H 0 se
f 2 (x )dx <
Ω
∞
mentre f (x ) è a derivata prima quadrato-integrabile, cioè f (x ) ∈ H 1 se:
2
f (x ) +
Ω
d f (x ) dx
2
dx < ∞
Nella figura 3.6 sono riportati i grafici di funzioni con diversa regolarità nei confronti dell’integrazione. La funzione riportata nella parte sinistra della figura appartiene a H 0 , ed è meno regolare rispetto a quella di destra, che, invece, appartiene a H 1 (avendo la derivata, a sua volta, in H 0 ). Mediante una formulazione debole si possono affrontare anche problemi le cui condizioni al contorno siano di tipo naturale, le quali possono essere soddisfatte anche in modo approssimato, come dimostrato nell’esempio seguente. Esempio 3.7 Si consideri un problema in cui si pone una condizione alla derivata prima di u
3.3 METODI DEI RESIDUI PESATI
51
Figura 3.6: Classi di regolarità di funzioni: funzione quadrato-integrabile (a), funzione a derivata prima quadrato integrabile (b). all’estremo destro dell’intervallo:
d2 u + u + x = 0, 0 < x < 1 dx 2 u |x =0 = 0 du ¯ = q dx x =1 la cui soluzione esatta è
(3.57) (3.58)
(3.59)
1 + ¯q sen x cos1 Vediamo dapprima il caso in cui entrambe le condizioni suddette siano u = x +
−
verificate dalla funzione approssimante, scegliendo una forma u = α1 + α2 x + α3 x 2 Imponendo le due condizioni, la soluzione approssimata si può esprimere in funzione di un solo parametro α = α3 :
¯ 2αx + αx 2 u = qx
−
Assumendo un peso alla Galerkin, che soddisfi le condizioni essenziali omogenee si ha: w = 2βx + βx 2 L’espressione del residuo pesato è: 1
0
R w dx =
1
0
2α + ¯qx − 2αx + αx 2 + x
2x + x 2 βdx = 0
−
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
52
da cui, stante l’arbitrarietà di β, si ottiene
α=
q + 1) − 25 ( ¯ 48
La funzione approssimante risulta, infine:
¯ u = qx
− 25 q + 1)(x 2 − 2x ) ( ¯ 48
La formulazione debole del problema si ottiene a partire da un’espressione del residuo pesato che tenga conto anche dell’approssimazione delle condizioni al contorno naturali 1
0
d2 u + (u + x ) w dx = dx 2
du − ¯q w dx
(3.60) x =1
In questo modo, infatti, integrando per parti si ottiene la scrittura: 1
du dw − (u + x )w dx = q¯ w x =1 (3.61) x x d d 0 ¯ dettato dalla condizione al contorno. ove è possibile imporre il valore di q
Assumiamo ora una funzione approssimata che non soddisfa la condizione naturale in x = 1, e un peso di Galerkin: u = α1 x 2 + α2 x w = β1 x 2 + β2 x Derivando e sostituendo si ottiene: 1
0
(2α1 x + α2 ) 2β1 x + β2
¯ β1 + β2 α1 x 2 + α2 x + x β1 x 2 + β2 x dx = q
−
Con la consueta procedura di scelta arbitraria dei coefficienti βk si ottiene quindi il sistema: 17 3 α ¯ + 14 q 1 15 4 (3.62) 3 2 α = q 1 ¯ + 2 4 3 3
le cui soluzioni sono:
α1 = α2 =
60 − 139 (1 + ¯q )
1 q ) − 139 (137 + 276 ¯
(3.63) (3.64)
La condizione al contorno naturale è solo approssimata, infatti, come si può verificare, la derivata della soluzione approssimata nel punto x = 1 non ¯ , ma risulta coincide con q
q x =1 =
du 1 q ) (17 + 156 ¯ = dx x =1 139
3.3 METODI DEI RESIDUI PESATI
53
Figura 3.7: Grafico delle soluzioni: esatta, u 0 , approssimata con condizione naturale esatta, u 1 e approssimata, u 2 .
3.3.5 Formulazioni inverse Attraverso un’ulteriore integrazione per parti della formulazione integrale del problema in esame, è possibile, nel caso dei problemi del secondo ordine, ribaltare i ruoli della funzione approssimante e della funzione peso. Tale procedura porta così alla cosiddetta formulazione inversa . Consideriamo ancora il problema d2 u + u + x = 0, 0 < x < 1 dx 2 u |x =0 = 0 u |x =1 = 0
(3.65) (3.66) (3.67)
Supponendo di commettere errori nell’approssimazione di entrambi i tipi di condizioni, l’espressione dei residui si può scrivere: 1
0
d2 u + (u + x ) w dx − dx 2
du − ¯q dx
1
dw ¯) (u − u dx
+
0
1
=0
(3.68)
0
Integrando per parti due volte si ottiene la cosiddetta formulazione inversa: 1
0
d2 w + w u dx + dx 2
1
0
xw dx + [qw ]10
dw ¯ u dx
1
−
=0
(3.69)
0
Si può trarre partito da questa formulazione, se si riesce a determinare una funzione peso tale da rendere nullo l’integrando del primo integrale dell’equazione (3.69). Questa scelta porta al metodo di T REFFTZ, applicato nell’esempio seguente.
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
54
Esempio 3.8 Determinare, mediante il metodo di Trefftz, una soluzione approssimata della formulazione inversa 1
0
d2 w + w u dx + dx 2
1
0
xw dx + [qw ]10
1
dw ¯ u dx
−
=0
(3.70)
0
In via preliminare, si deve determinare la funzione peso che permette di annullare il primo integrale, cioè la soluzione del problema
d2 w + w = 0 dx 2 La soluzione dell’equazione differenziale fornisce w = β1 cos x + β2 sen x
dw = −β1 sen x + β2 cos x dx
(3.71)
(3.72)
Sostituendo tali espressioni nella forma inversa si ricava: 1
0
xw dx + [qw ]10 = 0
(3.73)
nella quale le incognite sono i valori delle derivate q 0 , q 1 nei punti x = 0 e x = 1, rispettivamente. Operando le sostituzioni si ottiene: 1
0
x β1 cos x + β2 sen x dx + q 1 β1 cos1 + β2 sen1
−
q 0 β1 = 0
(3.74)
Data l’arbitrarietà di β1 e β2 , si perviene alle equazioni: 1
0 0
1
x cos x dx = (q 1 cos1 + q 0 )
(3.75)
x sen x dx = q 1 sen1
(3.76)
−
−
che forniscono i risultati (esatti in questo caso): q 0 =
1 , sen1
q 1 =
cos1 −1 sen1
3.3 METODI DEI RESIDUI PESATI
55
La formulazione inversa consente anche l’impiego delle cosiddette soluzioni fondamentali che rappresentano la risposta del sistema a sollecitazioni singolari. Nel caso dei problemi elastici, ad esempio, tali soluzioni sono i campi di spostamento generati nel corpo da forze concentrate applicate sul contorno o all’interno del dominio. Qualora queste soluzioni siano note in forma chiusa, è possibile costruire una procedura di soluzione efficace, come descritto nell’esempio seguente. Esempio 3.9 Riprendiamo il problema in forma integrale inversa: 1
0
d2 w + w u dx + dx 2
1
0
xw dx + [qw ]10
dw ¯ u dx
1
−
=0
(3.77)
0
a cui associamo condizioni al contorno essenziali omogenee. Consideriamo il termine tra parentesi quadre del primo integrale: diremo soluzione fondamentale la funzione w ∗ soluzione dell’equazione:
d2 w ∗ ∗ + w = ∆(x − ξi ) 2 dx
(3.78)
Sfruttando le proprietà della delta di Dirac è possibile scrivere le condizioni necessarie a determinare l’approssimazione. Sostituendo l’eq. (3.78) nell’espressione precedente, oltre alle condizioni omogenee si ha 1
0
d2 w ∗ ∗ + w u dx = 2 dx
1
0
(x − ξ)u dx
∆
(3.79)
Ricordando poi la proprietà della delta di Dirac: 1
(x − ξi )u dx = u (ξi )
∆
0
(3.80)
si ricavano le equazioni per le incognite q 0 e q 1 : u (ξ) +
1
0
xw ∗ dx + [qw ∗ ]10 = 0
(3.81)
La soluzione fondamentale, cioè la soluzione dell’equazione : (3.78) si trova in letteratura, ed è
1 2
w ∗ = sen |x ξi |
−
sostituendola nell’espressione precedente, scegliendo ξ1 = 0 e ξ2 = 1 si ricavano nuovamente i valori esatti delle derivate: q 0 =
1 , sen1
q 1 =
cos1 −1 sen1
56
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
3.4 Il metodo degli elementi finiti: il punto di vista matematico A completamentodella carrellata sui metodi di soluzione approssimatadei problemi a valori al contorno, in questa sezione si vuole offrire uno spunto rivolto a inquadrare i capitoli successivi, di orientamento più spiccatamente meccanico, nel contesto delle idee e delle tecniche fin qui accennate. Si intende mostrare come le diverse formulazioni possibili per l’approssimazione dei problemi forniscano, sotto certe condizioni, la stessa soluzione. In particolare, osserveremo che il principio dei lavori virtuali non è altro che una conveniente e familiare via per porre i problemi meccanici in una forma particolarmente adatta alla soluzione approssimata e quindi tale approccio sarà, come vedremo, alla base delle tecniche numeriche di analisi strutturale che vanno sotto il nome di Metodo degli Elementi Finiti.
3.4.1 Equivalenza delle formulazioni Consideriamo come esempio il sistema meccanico costituito da una barra rettilinea, vincolata ad un estremo, caricata assialmente mediante un carico distribuito f B e un carico concentrato all’estremo libero R . Siano, inoltre, assegnati il modulo di Young E , l’area della sezione retta A , la legge di distribuzione del carico f B (x ) = ax . L’incognita del problema, ovvero la variabile di stato, è lo spostamento assiale u (x ). Vedremo ora come il problema possa essere analizzato in forma differenziale oppure integrale in forma debole (o variazionale) e come quest’ultimo approccio sia equivalente all’applicazione del principio degli spostamenti virtuali. Formulazione differenziale. si scrive l’equazione indefinita di equilibrio, vali-
Figura 3.8: Esempio affrontato con diversi metodi di soluzione del problema.
3.4 IL METODO DEGLI ELEMENTI FINITI: IL PUNTO DI VISTA MATEMATICO
57
da in ogni punto della trave, associata alle sue condizioni al contorno:
d2 u + ax = 0, 0 < x < L , dx 2 u |x =0 = 0 , du E A = R . dx x =H
E A
L’equazione differenziale viene, come è noto, determinata studiando l’equilibrio di un elemento infinitesimo di trave e vale in tutti i punti interni del dominio. Ad essa si aggiungono le condizioni al contorno, rispettivamente essenziale, per x = 0, e naturale, per x = L . La soluzione del problema fornisce la soluzione esatta u (x ) =
ax 3
1 E A
− 6
+ R +
aL 2
2
x
Si può notare che la funzione u è continua e derivabile due volte con continuità, come richiesto dall’equazione differenziale. Assegnato il carico f B = a x , u appartiene quindi allo spazio delle funzioni continue e derivabili due volte che soddisfano le condizioni al contorno. Formulazione variazionale. Il potenziale totale si ottiene come somma del ter-
mine relativo all’energia di deformazione: 1 U = 2
L
du 2 E A dx dx
0
con quello del potenziale del carico: L
V =
0
uax dx + Ru |x =L
Per cui il potenziale totale risulta: 1 Π = U − V = U = 2
L
du 2 E A dx dx
L
− 0
0
uax dx + Ru |x =L
Si impone ora la stazionarietà del potenziale totale sotto condizione di rispetto delle condizioni essenziali per la funzione e per la sua variazione. δΠ = 0 u |x =0 = 0 ,
δu |x =0 = 0 .
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
58
si noti che la variazione δ u è deve, in generale, soddisfare una versione omogenea delle condizioni essenziali (cioè δu = 0), in quanto laddove la funzione u ha un valore imposto, la sua variazione deve essere nulla. In questo caso, essendo omogenee le condizioni essenziali per la u , le condizioni della variazione δu coincidono con quelle della funzione. Integrando per parti l’equazione δΠ = 0 e imponendo le condizioni al contorno naturali ed essenziali, seguendo un procedimentosimile a quello illustrato dettagliatamente nell’Esempio XX, si ottiene la formulazione differenziale del problema illustrata al paragrafo precedente. Dopo le opportune integrazioni si perviene alla scrittura L
d2 u du E A 2 + ax δu dx + R − E A δu dx dx
0
=0
(3.82)
x =L
da cui, scegliendo le opportune variazioni arbitrarie, si ottengono l’equazione differenziale di campo e la condizione al contorno naturale. La formulazione variazionale, quindi, si può ricondurre a quella differenziale. Esse sono esattamente equivalenti se i dati sono sufficientemente regolari entrambe, quindi, hanno la stessa soluzione “esatta”. Principio dei lavori virtuali L’applicazione del PLV al caso in esame consente
di scrivere l’uguaglianza tra il lavoro virtuale interno, associato allo sforzo normale N (x ) e quello esterno compiuto dalle forze di volume f B e all’estremo libero R : L
0
N ˆdx =
L
0
ˆ dx + R u ˆ |x =L f B u
(3.83)
ˆ e ˆ sono indicati rispettivamente gli spostamenti virtuali e le dove con u deformazioni virtuali da essi derivate mediante le equazioni cinematiche. Tramite l’equazione costitutiva elastico-lineare si può esprimere lo sforzo normale in funzione della deformazione reale: N = E A = E A du /dx , ottenendo la scrittura: L
0
u du d ˆ dx = E A dx dx
L
0
ˆ dx + R u ˆ |x =L f B u
(3.84)
Integrando per parti l’integrale a primo membro, si ricava du ˆ E A u dx
L
L
− 0
0
d2 u ˆ dx = E A 2 u dx
L
0
ˆ dx + R u ˆ |x =L f B u
(3.85)
3.4 IL METODO DEGLI ELEMENTI FINITI: IL PUNTO DI VISTA MATEMATICO
59
e, raccogliendo e riordinando i termini, ponendo anche f B = ax e tenenˆ = 0, si ottiene do conto che, per x = 0 deve essere u L
d2 u du ˆ E A 2 + ax ˆ u dx + R − E A u dx dx
0
=0
(3.86)
x =L
È immediato riconoscere l’equivalenza tra questa espressione e la (3.82), ˆ (x ) come variazione interpretando il campo degli spostamenti virtuali u δu (x ) degli spostamenti reali u (x ). A commento di quanto trovato si possono svolgere le seguenti osservazioni: – le tre formulazioni descritte sono totalmente ed esattamente equivalenti. – La soluzione trovata è l’unica soluzione del problema. – Se si adotta un approccio variazionale si procede attraverso una formulazione debole che richiede solo la derivabilità al primo ordine della funzione incognita e della sua variazione, quindi è possibile ottenere soluzioni in una classe di funzioni più ampia di quella necessaria per la formulazione differenziale. – Il principio degli spostamenti virtuali (detto anche dei lavori virtuali) permette di tener conto di singolarità presenti nel problema, come carichi concentrati o discontinuità del materiale, dove la derivata prima di u è discontinua e la formulazione differenziale cade in difetto. – La soluzione di un problema meccanico mediante metodi di approssimazione viene riportata a una formulazione debole alla Galerkin – La scelta delle funzioni base come funzioni definite su sotto-domini porta al FEM. In modo più astratto e generale, detto V lo spazio vettoriale di funzioni nel quale ricercare la soluzione (dette trial functions ) , che sarà meglio precisato nel seguito, il problema variazionale si può porre nella forma: trovare u ∈ V , tale che
〉 ∀w ∈ V ˜
a (u , w ) = f , w ,
〈
(3.87)
dove: a (u , w ) è la forma bilineare del problema che, nel caso specifico qui in esame diventa: L du dw a (u , w ) = E A dx (3.88) dx dx 0
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
60
mentre 〈 f , w 〉 è la forma lineare del problema che assume nel caso particolare qui esaminato, l’espressione:
〈 f , w 〉 =
L
0
f B w dx + Rw |x =L
(3.89)
Si noti che l’equazione (3.87) deve valere per ogni funzione w (dette test func˜ , che in generale, ma tions ), appartenente a un opportuno spazio vettoriale V non necessariamente, può essere diverso da V e la cui scelta caratterizza i diversi metodi di soluzione approssimata. La scelta dello spazio V delle funzioni ammissibili deve essere fatta comprendendo in esso le funzioni che soddisfano le condizioni al contorno essenziali. Inoltre, tali funzioni devono possedere una regolarità sufficiente a garantire l’esistenza delle derivate coinvolte nella formulazione variazionale del problema (3.87). Altro requisito molto importante, che coinvolge anche la regolarità delle funzioni test, è che gli integrali definiti devono dare un valore finito.
3.4.2 Formulazione di Galerkin Prendiamo come esempio l’espressione in forma debole del problema variazionale seguente L
0
du dw − uw − xw dx = 0, dx dx
∀w ∈ H 01
dove H 01 indica la classe delle funzioni ammissibili con dominio Ω =]0,1[ che soddisfano le condizioni al contorno essenziali omogenee. Lo spazio H 01 è uno spazio vettoriale normato. È un esempio degli spazi di S OBOLEV (il cui studio esula i contenuti qui trattati, ma che risulta fondamentale nell’analisi funzionale). Nella notazione H 01, il pedice 0 indica che le funzioni sono nulle agli estremi, mentre l’apice 1 significa che il quadrato della derivata prima possiede un integrale finito se calcolato sul dominio Ω. Lo spazio H 01 ha dimensione infinita , vale a dire che non esiste un insieme finito di funzioni base la cui combinazione lineare permetta di esprimere una qualsiasi funzione dello spazio. Supponiamo però che sia possibile rappresentare ogni funzione w mediante una combinazione lineare di funzioni φ i , fra loro indipendenti, appartenenti all’insieme φi delle funzioni base di H 01 . w (x ) =
∞
βi φi (x ) .
i =1
In generale si suppone che una data funzione w si possa approssimare bene quanto si voglia prendendo un numero sufficientemente grande di funzioni
3.4 IL METODO DEGLI ELEMENTI FINITI: IL PUNTO DI VISTA MATEMATICO
61
base. Prendendo solo N termini della serie, si introduce un’approssimazione w N di w : w N (x ) =
N
βi φi (x ) .
i =1
Le N funzioni base {φ1 ... φN } definiscono un sottospazio N-dimensionale di H 01 , detto H 0(N ) . Per definizione di sottospazio di uno spazio lineare, ogni funzione in H 0(N ) è combinazione lineare delle N funzioni base. Ad esempio, la terna {φ1φ2 φ3 } è la base di H 0(3) e ogni funzione appartenente a H 0(3) si può esprimere nella forma w 3 (x ) = β1 φ1 (x ) + β2 φ2 (x ) + β3 φ3 (x ) .
Il metodo di G AL ER KI N ricerca le soluzioni approssimate al problema “accontentandosi” di muoversi solamente nel sottospazio H 0N , assumendo come approssimazione l’espressione u N (x ) =
N
αi φi (x ) ,
i =1
cioè usando per l’approssimazione la stessa base utilizzata per definire la funzione test, w . Ne deriva il problema debole approssimato L
0
du N dw N − u N w N − xw N dx = 0, dx dx
(N ) ∀w N ∈ . H 0
Poiché le φ i sono note, u N risulterà noto una volta determinati i parametri α i (gradi di libertà). Introducendo nell’espressione precedente le funzioni u N (x ) e w N (x ) si ottiene:
L
0
d dx
N
− −
αi φi (x )
i =1
N
d dx
N
βi φi (x )
i =1
N
x
αi φi (x )
i =1
βi φi (x )
dx = 0,
i =1
N
βi φi (x )
i =1
(3.90)
∀βi , i = 1,..., N .
Ponendo φ (x ) = d φ(x )/dx e riordinando, mettendo in evidenza i coefficienti arbitrari βi , si ottiene:
− N
1
N
βi
i =1
0
i =1
1
0
φi (x )φ j (x ) φi (x )φ j (x ) dx α j
−
x φi (x )dx = 0,
∀βi , i = 1,2,..., N .
(3.91)
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
62
Ponendo poi: 1
K i j = F i =
0
1
0
φi (x )φ j (x ) − φi (x )φ j (x ) dx ,
(3.92)
x φi (x )dx ,
(3.93)
l’espressione precedente assume la forma più compatta
N
i =1
N
βi
i =1
K i j α j F i = 0,
−
∀βi , i = 1,..., N
(3.94)
che rappresenta l’espressione discretizzata del problema variazionale è. A questo punto, si può osservare che l’espressione (3.94) deve essere comunque soddisfatta, per qualsiasi scelta dei i parametri βi che, ricordiamo, precisano le funzioni peso in H 0N . Seguendo questa osservazione, si vede che tale equazione permette di scrivere N condizioni per le incognite α j . Si scelga, ad esempio, la particolare ennupla di coefficienti: β1 = 1, β 2 = 0, β 3 = 0, ... , β N = 0 .
Sostituendo questi valori nella (3.94) si ottiene: N
K 1 j α j = F 1
i =1
Analogamente, usando ennuple di coefficienti tutti nulli tranne uno, si ricavano N equazioni, in sintesi: β1 = 1, β2 = 1,
βi = 0 per i = 1
=
i =1 N
βi = 0 per i = 2
=
···
βN = 1,
N
⇒ ⇒ ⇒
··· =
K 2 j α j = F 2 ,
i =1 N
βi = 0 per i = N
K 1 j α j = F 1 ,
K N j α j = F N .
i =1
Le N equazioni concorrono a formare infine il sistema lineare: N
K i j α j = F i ,
i = 1,..., N
i =1
I coefficienti K i j costituiscono la matrice di rigidezza del problema: K = [K i j ]
(3.95)
3.4 IL METODO DEGLI ELEMENTI FINITI: IL PUNTO DI VISTA MATEMATICO
63
Come si può facilmente osservare, tale matrice è simmetrica (scambiando i e j il risultato non cambia): 1
K i j = K j i =
0
φi (x )φ j (x ) − φi (x )φ j (x ) dx ,
e risulta anche invertibile, poiché le funzioni φi sono state scelte linearmente indipendenti e quindi risultano indipendenti anche le equazioni). I termini F i sono gli elementi del vettore dei carichi : F = [F i ]
La soluzione del sistema (3.95) fornisce i coefficienti α j che definiscono la soluzione approssimata u N (x ): N
α j =
(K −1 )i j F i
i =1
dove(K −1 )i j sono gli elementi dell’inversa di K (cioè della matrice K −1 tale che K −1 K = I, dove I è la matrice identità).
3.4.3 Il metodo degli elementi finiti Il metodo di Galerkin consente di calcolare i coefficienti della soluzione approssimata una volta che siano state fissate le funzioni base, ma non fornisce alcuna indicazione per la scelta di una base di funzioni appropriata. D’altro canto la bontà dell’approssimazione dipende fortemente dalle proprietà delle funzioni base e una scelta cattiva può portare a una matrice di rigidezza mal condizionata con conseguenti difficoltà nella soluzione del sistema lineare. Inoltre, particolarmente nel caso di problemi al contorno a due o tre dimensioni è difficile trovare funzioni base che rispettino le condizioni al contorno in domini con geometria complicata. Il Metodo degli Elementi Finiti (FEM) fornisce una tecnica generale per costruire funzioni base per l’applicazione del metodo di Galerkin alla soluzione approssimata di problemi al contorno. Le φi vengono definite a tratti su regioni del dominio dette elementi finiti e su ogni elemento si sceglie per le φ i un andamento molto semplice, ad esempio polinomiale. Si consideri, ad esempio, una suddivisione del dominio Ω =]0,1[ in quattro elementi Ωi , i = 1,2,3,4 (Figura 3.9). La lunghezza di ogni elemento sia h i = h , scelta costante per semplicità. Gli estremi degli elementi prendono il nome di nodi. Gli estremi degli elementi prendono il nome di nodi o punti nodali . L’insieme di nodi ed elementi si dice mesh di elementi finiti. Le funzioni base del metodo degli elementi finiti sono scelte in modo da soddisfare tre criteri:
64
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
Figura 3.9: Mesh di elementi finiti per un problema monodimensionale. 1. sono generate da funzioni semplici definite a tratti, cioè elemento per elemento; 2. sono abbastanza regolari da consentire l’integrabilità dei prodotti delle derivate e devono soddisfare la versione omogenea delle condizioni al contorno essenziali; 3. sono scelte in modo che i parametri α i ai che definiscono la soluzione approssimata u N (x ) siano proprio i valori di u N calcolati nei nodi. L’esempio più semplice si basa su funzioni base lineari a tratti con valore unitario in un nodo e nullo in tutti gli altri, i cui grafici sono illustrati nella Figura 3.10. Le funzioni φ1 , φ2 , φ3 prescelte soddisfano il primo criterio in quanto, all’interno di ciascun elemento, si tratta di polinomi di primo grado o di grado zero (costanti). L’espressione analitica delle funzioni base si ricava per interpolazione lineare nei tratti [i − 1, i ,] e [i , i + 1,].
φi (x ) =
0
per x x i −1 x x i −1 per x i −1 x x i h i x i +1 x per x i x x i +1 h i +1 0 per x x i +1
≤
−
−
≤ ≤ ≤ ≤ ≥
dove h i = x i − x i −1 è la lunghezza dell’elemento Ωi . Per quanto riguarda la classe di integrabilità delle funzioni base, si può verificare che esse appartengono allo spazio delle funzioni la cui derivata prima è quadrato-integrabile, e hanno valore nullo agli estremi del dominio, ovvero appartengono a H 01 . Come è ovvio, le φi si annullano, come richiesto, agli estremi dell’intervallo, cioè per x = 0 e x = 1, per costruzione. Più delicata è la verifica del rispetto delle condizioni di regolarità, cioè che le φi e le loro derivate φi siano quadratointegrabili.
3.4 IL METODO DEGLI ELEMENTI FINITI: IL PUNTO DI VISTA MATEMATICO
65
Figura 3.10: Funzioni base per la discretizzazione del problema monodimensionale. Per quanto riguarda le φ i come si può facilmente verificare applicando ad esempio la regola di Simpson: 1
x x i 1 2 φi (x ) dx = dx + h i x i −1 h i + h i +1 x i
x i
− − 0
2
=
<
3
x i +1
x i −1
− x
h i +1
2
dx
∞
Per quanto concerne, invece, le derivate φi , il punto delicato della questione è la presenza della discontinuità angolare in x = x i , dove la derivata prima non è definita (si veda la Figura 3.11). Le derivate prime delle funzioni base risultano:
φi (x ) =
0 1 h i
−1
h i +1
0
per x < x i −1
per x i −1 < x < x i per x i < x < x i +1
per x > x i +1
3 SOLUZIONE APPROSSIMATA DI PROBLEMI AL CONTORNO
66
Figura 3.11: Integrabilità delle derivate delle funzioni di interpolazione. Si verifica però facilmente che la funzione φ2i è integrabile: 1
x i
1
x i −1
h i
2
x i
1
x i −1
h i +1
2
0
2
φi (x ) dx = =
1
h i
+
1
h i +1
dx +
<
dx
∞
Passando infine al terzo criterio, verifichiamo che i parametri dell’approssimazione corrispondono ai valori nodali della variabile incognita, cioè α i = u N (x i ). A questo fine, basta notare che Si osservi che, posto u j = u N (x j ), si ha: u j = u N (x j ) =
N
αi φi (x j ) = α1 · 0 + · · · + α j · 1 + · · · + αN · 0 = α j
i =1
e quindi i coefficienti della combinazione lineare αi coincidono con i valori nodali u i della funzione incognita, e quindi quest’ultima si può esprimere come interpolazione fra i valori nodali: u N (x j ) =
N
i =1
αi φi (x j ) =
N
i =1
u i φi (x j )
Capitolo 4 Il metodo degli elementi finiti 4.1 Introduzione 4.2 Discretizzazione delle equazioni di equilibrio Consideriamo un corpo continuo che occupa il dominio Ω di uno spazio euclideo tridimensionale, ove è fissato un sistema di riferimento cartesiano OXY Z (Figura 4.1). La frontiera Γ del corpo è suddivisa in due parti: Γu , ove sono imposti gli spostamenti dei punti, e Γ f , sulla quale sono applicate forze di superficie note, f . La trattazione che segue, basata sull’applicazione del principio dei lavori virtuali, consente di introdurre, oltre a forze distribuite nel volume, b, anche l’azione di forze concentrate in alcuni punti del corpo, Fi . Obiettivo dell’analisi è descrivere il problema della ricerca della configurazione assunta dal corpo, definita attraverso il campo vettoriale dello spostamento, U( X , Y , Z ): U ( X , Y , Z ) U = V ( X , Y , Z ) . W ( X , Y , Z )
Come è noto, il campo delle deformazioni infinitesime può essere descritto tramite le componenti rispetto al sistema OXY Z raccolte nel vettore ( X , Y , Z ):
X X Y Y Z Z = . γ X Y γY Z γ Z X
La relazione differenziale tra il vettore spostamento e vettore deformazione si 67
4 IL METODO DEGLI ELEMENTI FINITI
68
Figura 4.1: Il corpo solido in equilibrio. può esplicitare introducendo l’operatore
∂=
per cui, in forma compatta:
∂ ∂ X
0
0
0
∂ ∂Y
0
0
0
∂ ∂Y
∂ ∂ X ∂ ∂ Z
0 ∂ ∂ Z
∂ ∂ Z ,
0 ∂ ∂Y ∂ ∂ X
0
( X , Y , Z ) = ∂U( X , Y , Z ).
(4.1)
Lo stato di tensione viene, analogamente, descritto dal vettore delle componenti di tensione σ( X , Y , Z ):
σ X X σY Y σ Z Z σ= . σ X Y σY Z σZ X
4.2 DISCRETIZZAZIONE DELLE EQUAZIONI DI EQUILIBRIO
69
Con queste notazioni, le equazioni indefinite di equilibrio si scrivono, in forma simbolica: T ∂ σ + b = 0, T N σ − f = 0,
in Ω,
(4.2)
su Γ f ,
(4.3)
dove si è introdotta la matrice contenente i coseni direttori del versore normale uscente da Γ f , n = [n X , n Y , n Z ]T :
N =
n X
0
0 0
n Y
0 0
0
n Z
n Y n X 0 0 n Z n Y n Z 0 n X
Per giungere a una formulazione debole del problema, si può applicare il principio dei lavori virtuali al sistema equilibrato di forze-tensioni. È noto, infatti, che, in tali condizioni, il lavoro virtuale delle forze applicate al corpo deve eguagliare il lavoro virtuale delle tensioni, per qualsiasi sistema di spostamentideformazioni, purchécompatibile con i vincoli e conle condizioni di congruen¯ ( X , Y , Z ) il campo za imposte dalla continuità del corpo. Indicati quindi con U di spostamenti virtuali e con ¯( X , Y , Z ) quello delle corrispondenti deformazioni virtuali, dovrà risultare:
Ω
¯ T σdΩ =
Ω
¯ T
U bdΩ +
Γ f
¯ T f dΓ + U
n
¯ i T Fi U
(4.4)
i
Si osservi che l’integrale relativo al lavoro virtuale delle forze di superficie si estende alla sola parte Γ f della frontiera del corpo, in quanto sulla restante parte Γu gli spostamenti virtuali, dovendo rispettare la condizione al contorno es¯ T f è senziale, sono nulli e, conseguentemente, su tale dominio l’integrando U pari a zero. Sfruttando la proprietà additiva degli integrali è possibile suddividere i domini di integrazione indicati nell’eq.(4.4) in regioni di forma poligonale,Ωk , k = 1,... m dette elementi , i cui vertici sono detti nodi o punti nodali. In ciascun elemento, inoltre, può essere conveniente introdurre un sistema di riferimento locale oxy z , rispetto al quale definire le diverse funzioni coinvolte. Ad esempio, con b (m ) sarà indicata la forza volumica agente sull’elemento m , espressa nel suo sistema di riferimento locale, mentre la funzione spostamento, sempre nel riferimento locale, sarà indicata con u(m ) .
4 IL METODO DEGLI ELEMENTI FINITI
70
Figura 4.2: Suddivisione del corpo in elementi. Supponendo, per semplicità, che la frontiera del corpo sia poligonale,1 e che le forze concentrate siano applicate nei nodi del reticolo, l’espressione del principio dei lavori virtuali diventa: n e
m n e m
Ω(m )
Ω(m )
(m ) T (m )
dΩ(m ) =
¯ (m ) T (m )
(m )
¯
u
σ
b
dΩ
n e
+
m
(m ) (m ) Γ ,...,Γ q 1
T (m )
¯ f u
(m )
dΓ
n
+
(4.5) i T i
¯ F u
i
dove l’integrale di superficie comprende le facce Γ(i m ) , i = 1,..., q dell’elemento m che appartengano eventualmente a Γ f . A questo punto si introduce il concetto più importante della formulazione a elementi finiti: si esprime la funzione incognita spostamento u (m ) , elemento per elemento, mediante un’approssimazione: scelte opportune funzioni di interpolazione H(m ) (x , y , z ) si scrive: ˆ u(m ) (x , y , z ) = H(m ) (x , y , z )U
(4.6)
ˆ raccoglie gli spostamenti nodali del sistema, cioè tutte le comdove il vettore U ponenti di spostamento dei nodi, espresse nel riferimento globale OXY Z : ˆ = [U 1V 1W 1 ... U N V N W N ]T = [U 1U 2 ... U 3N ]T U 1 In
caso contrario l’unione degli elementi finiti lascia fuori una parte del dominio Ω (v. Figura 4.2), che si potrebbe fare tendere a zero infittendo opportunamente la suddivisione.
4.2 DISCRETIZZAZIONE DELLE EQUAZIONI DI EQUILIBRIO
71
avendo indicato con N il numero di nodi della discretizzazione. Il campo di deformazione approssimato si ottiene per derivazione del campo di spostamento approssimato (m ) (x , y , z ) = ∂u(m ) (x , y , z ) per cui, sfruttando l’eq.(4.6) si ha:
(m )
ˆ = B(m ) (x , y , z )U ˆ (x , y , z ) = ∂H(m ) (x , y , z )U
(4.7)
dove si è introdotta la matrice B(m ) = ∂H(m ) che contiene le derivate delle funzioni di interpolazione dell’elemento m . Si noti che tali derivate sono fatte rispetto alle coordinate del sistema di riferimento globale OX Y Z , in quanto le componenti di deformazione contenute nel vettore (m ) sono espresse in tale riferimento. Adottando per gli spostamenti e le deformazioni virtuali le stesse interpolazioni impiegate per gli spostamenti reali, eq.(4.6), e per le deformazioni reali, eq.(4.7), sostituendo tali relazioni nell’eq.(4.5) si ottiene: n e
m n e m
Ω(m )
Ω(m )
¯ˆ T B(m ) T σ(m ) dΩ(m ) = U ¯ˆ T H(m ) T b(m ) dΩ(m ) + U
n e
m
(m ) (m ) 1 ,...,Γq
¯ˆ T H(m ) T f (m ) dΓ(m ) + ¯U ˆ T Fi U
(4.8)
Γ
L’eq.(4.8) contiene il campo delle tensioni negli elementi, σ(m ) , ancora da determinare; scopo di questa trattazione è istituire un problema in cui l’incognita sia il campo degli spostamenti o meglio, trattandosi di un procedimento di discretizzazione, il vettore degli spostamenti nodali . A questo scopo è necessario eliminare il riferimento alla tensione σ(m ) a favore di una presentazione del problema appunto in termini di spostamento. Perciò si deve assumere un legame costitutivo , cioè una legge che metta in relazione le tensioni con le corrispondenti deformazioni. Per rimanere nell’ambito di una trattazione lineare , ipotizziamo che tale legame costitutivo sia rappresentato dalla legge di Hooke. Assumiamo cioè che si possa scrivere: (m ) (m ) (m ) I (m ) σ (4.9) =C +σ dove si è introdotta la matrice costitutiva C(m ) dell’elemento, mentre σ I (m ) ha il significato di uno stato tensionale preesistente alla deformazione (m ) , potendo rappresentare ad esempio lo stato tensionale originario in un problema geotecnico o uno stato di autotensione nel caso di solidi in stato di coazione. Per essere utilizzata nell’eq.(4.8), anche l’eq.(4.9) deve essere espressa in funzione degli spostamenti nodali. Sostituendovi l’eq.(4.7) si ha (m )
σ
=C
(m ) (m ) ˆ
B
U + σI (m )
(4.10)
4 IL METODO DEGLI ELEMENTI FINITI
72
e quindi l’eq.(4.8) diventa n e
m n e m
Ω(m )
n e
¯ˆ T B(m ) T C(m ) B(m ) U ˆ dΩ(m ) = U
Ω(m )
m
n e
¯ˆ T H(m ) T f (m ) dΓ(m ) − U
(m ) (m ) 1 ,...,Γq
Γ
Ω(m )
m
¯ˆ T H(m ) T b(m ) dΩ(m ) + U
¯ˆ T B(m ) T σI (m ) dΩ(m ) + ¯U ˆ T Fi U
(4.11)
Mettendo in evidenza i vettori degli spostamenti nodali reali e virtuali che, essendo costanti, si possono portare fuori dagli integrali e dalle sommatorie e si può riscrivere l’eq.(4.11) nella forma ¯ˆ T U
n e
m
Ω(m )
(m ) T (m ) (m )
B
C
n e
m
B
(m )
dΩ
(m ) T (m )
(m ) (m ) 1 ,...,Γ q
H
f
n e
ˆU = U ˆ¯ T
m
(m )
dΓ
n e
−
Γ
m
Ω(m )
H(m ) T b(m ) dΩ(m ) +
(m ) T I (m )
Ω(m )
B
σ
(m )
dΩ
(4.12)
i
+F
L’equazione scalare (4.12) rappresenta il principio dei lavori virtuali per il sistema discretizzato. L’uguaglianza da essa indicata deve essere soddisfatta ˆ¯ di spostamenti virtuali dei nodi. Per comprendere la per qualsiasi vettore U portata di questa affermazione, conviene introdurre le definizioni seguenti: matrice di rigidezza
n e
K =
Ω(m )
m
B(m ) T C(m ) B(m ) dΩ(m )
(4.13)
H(m ) T b(m ) dΩ(m )
(4.14)
vettore delle forze di volume n e
RB =
Ω(m )
m
vettore delle forze di superficie n e
RS =
m
(m ) (m ) 1 ,...,Γq
H(m ) T f (m ) dΓ(m )
(4.15)
Γ
vettore delle tensioni iniziali n e
RI =
m
Ω(m )
B(m ) T σI (m ) dΩ(m )
(4.16)
vettore delle forze concentrate RC = Fi
(4.17)
4.2 DISCRETIZZAZIONE DELLE EQUAZIONI DI EQUILIBRIO
73
per cui l’eq.(4.12) si può riscrivere in forma compatta ¯ˆ T K U ˆ =U ˆ¯ T [RB + RS − RI + RC ] U
(4.18)
o anche, con la posizione R = RB + RS − RI + RC in maniera ancora più sintetica ¯ˆ T K U ˆ =U ˆ¯ T R U
(4.19)
L’equazione (4.19) deve essere soddisfatta per qualsiasi vettore di spostaˆ¯ . Sfruttando questa condizione, è possibile ottenere un sistemento virtuale U ˆ . Scegliendo ad esempio il particolare ma di equazioni lineari nell’incognita U vettore ¯ˆ = 1 0 ... 0 T U 1
è facile verificare che l’espressione (4.19) si riduce alla sola equazione N
¯ˆ = R K 1 j U 1 j
j =1
Scegliendo invece il vettore ¯ˆ = 0 1 ... 0 U 2 si ottiene
N
T
¯ˆ = R K 2 j U 2 j
j =1
ˆ¯ N = 0 0 . . . N T si ottengono Procedendo analogamente fino a scegliere U quindi N equazioni che devono essere soddisfatte contemporaneamente per rispettare la condizione (4.19). In definitiva, quindi, il procedimento porta alla scrittura del sistema lineare ˆ =R K U (4.20)
È importante osservare che i vettori dei carichi e la matrice di rigidezza globali si ottengono per somma delle matrici relative ai singoli elementi. Questo è possibile poiché si sta operando in modo da tener conto di tutti i gradi di libertà della struttura, cioè i vettori hanno dimensione n × 1, mentre la matrice di rigidezza ha dimensione n × n . Nel caso in cui ci si riferisca a matrici valutate con riferimento alla numerazione locale dei gradi di libertà, è necessario espandere le matrici stesse alle dimensioni opportune prima di poterle sommare fra loro.
4 IL METODO DEGLI ELEMENTI FINITI
74
4.3 Equilibrio dinamico Nel caso in cui siano presenti accelerazioni, conviene applicare il principio di D ’A L EMBERT, introducendo le forze d’inerzia tra i carichi applicati. Si approssima il campo delle accelerazioni con le stesse funzioni di interpolazione usate per gli spostamenti: ¨ ¨ (m ) (x , y , z ) = H(m ) (x , y , z )U u (4.21) Indicando con ρ (m ) la massa volumica, il vettore dei carichi di volume, diventa quindi: n e
RB =
m
Ω(m )
¨ ]dΩ(m ) = RB − MU ¨ H(m ) T [b(m ) − ρ (m ) H(m ) U
(4.22)
dove si è introdotta la matrice n e
M=
Ω(m )
m
ρ (m ) H(m ) T H(m ) dΩ(m )
(4.23)
detta matrice di massa della struttura. Con questa posizione, l’equazione di equilibrio dinamico della struttura diventa il sistema di equazioni differenziali : ¨ + KU = R MU
(4.24)
Alla stessa stregua delle forze d’inerzia si possono trattare altre forze di volume, come quelle generate dalla viscosità del materiale. Come è noto, le forze viscose sono proporzionali alla velocità delle particelle, quindi, introducendo un coefficiente di viscosità κ(m ) si può valutare la matrice di smorzamento n e
C=
Ω(m )
m
κ(m ) H(m ) T H(m ) dΩ(m )
(4.25)
per cui il vettore dei carichi di volume diventa
n e
RB =
m
Ω(m )
¨ − κ(m ) H(m ) U ˙ ]dΩ(m ) H(m ) T [b(m ) − ρ (m ) H(m ) U
¨ − CU ˙ = RB − MU
(4.26)
e quindi le equazioni differenziali che esprimono l’equilibrio dinamico della struttura assumono la forma: ¨ + CU ˙ + KU = R MU
(4.27)
4.4 IL METODO DELLE COORDINATE GENERALIZZATE
75
4.4 Il metodo delle coordinate generalizzate Uno degli aspetti più importanti nella formulazione a elementi finiti è la scelta delle funzioni di interpolazione, contenute, come abbiamo visto, nella matrice H(m ) di ogni elemento. Vedremo in seguito i requisiti che tali funzioni devono soddisfare e una tecnica che consente di generarle facilmente e in modo sistematico (v. sez. 5.2). Come primo approccio al problema, e per fissare le idee finora esposte mediante alcuni esempi, presentiamo qui il metodo delle coordinate generalizzate , che è stato impiegato parecchio nel primo periodo dell’applicazione del metodo degli elementi finiti. Sostanzialmente, questa tecnica consiste nello scegliere la forma delle funzioni di interpolazione, espresse in un sistema di riferimento opportuno, generalmente locale, per ogni elemento. Tali funzioni sono definite tramite espressioni polinomiali i cui coefficienti, detti appunto coordinate generalizzate , vanno determinati in modo che il campo di spostamenti che ne risulta sia coerente con i valori assunti dagli spostamenti nodali.
Figura 4.3: Discretizzazione di una struttura monodimensionale (Esempio 4.1).
4 IL METODO DEGLI ELEMENTI FINITI
76
Esempio 4.1
Esempio 4.2
Una lastra piana di forma quadrata (Figura 4.4) e caricata nel suo piano viene discretizzata con una suddivisione in quattro elementi uguali. Si vogliono valutare, per l’elemento 2 le espressioni della matrice di interpolazione, H(2) , della matrice delle derivate, B(2) e della matrice costitutiva, C(2). La lastra risulta sede di uno stato tensionale piano. Il primo passo consiste nel numerare tutti i gradi di libertà della struttura. Come illustrato nella Figura 4.4, le incognite sono le due componenti degli spostamenti nodali U i e V i , nelle direzioni degli assi X e Y , rispettivamente. Per comodità di esecuzione automatica delle procedure di calcolo, è conveniente ribattezzare tali componenti e raccoglierle in un vettore UT = [U 1 V 1 U 2 V 2 ... U 9 V 9 ] = [U 1 U 2 ... U 18 ]
L’elemento 2, con la numerazione globale e locale dei gradi di libertà è illustrato nella Figura 4.6: il vettore degli spostamenti nodali dell’elemento, nella numerazione locale è ˆ T = [u 1 u 2 u 3 u 4 | v 1 v 2 v 3 v 4 ] u L’approssimazione del campo degli spostamenti entro l’elemento si può scrivere u (x , y ) = α1 + α2 x + α3 y + α4 x y
v (x , y ) = β1 + β2 x + β3 y + β4 x y
(4.28) (4.29)
dove ciascuna delle componenti è considerata funzione delle coordinate locali x , y tramite quattro parametri. Tali parametri, rispettivamente α 1 ... α4 e β1 ... β4 sono detti coordinate generalizzate. In forma matriciale le precedenti espressioni si possono scrivere: α1 α2 α3 u (x , y ) 1 x y x y 0 0 0 0 α4 = v (x , y ) 0 0 0 0 1 x y x y β1 β2 β3 β4
(4.30)
4.4 IL METODO DELLE COORDINATE GENERALIZZATE
77
P
P 3
6
9
!
m c 4
$
2
5
V , Y
8
#
1
"
X,U
4 cm
4
7
2 cm
2 cm
Figura 4.4: Discretizzazione di una lastra piana (Esempio 4.2). 3
V 3 U 3
2
!
V 2 U 2
1
#
V 1 U 1
6
9
V 6 U 6
5
V 9 U 9
$
V 5 U 5
4
"
V 4 U 4
8
V 8 U 8
7
V 7 U 7
2 1
6
U
3 U 5
4
U
!
2 U 3
2
U
#
U
6U 11
0 1
U
1 U 1
$
5 U 9
8
U
8 1
"
4 U 7
U
9U 17
6 1
U
8U 15
4 1
U
7U 13
Figura 4.5: Numerazione globale dei gradi di libertà (Esempio 4.2). 3
V 3 U 3
6
V
6 U 6
1
2
v 1 u 1
v 2 u 2 v , y
x,u 2
V 2 U 2
5
V 5 U 5
3
v 3 u 3
4
v 4 u 4
Figura 4.6: Numerazione locale dei gradi di libertà per l’elemento 2. (Esempio 4.2).
4 IL METODO DEGLI ELEMENTI FINITI
78
Le funzioni contenute nella matrice Φ prendono il nome di funzioni base. Ponendo inoltre: 1 x y xy 0 0 0 0 Φ= (4.31) 0 0 0 0 1 x y xy e
α=
α1 α2 α3 α4 β1 β2 β3 β4
si scrive, in forma compatta:
T
(4.32) Il problemaconsiste oranel determinare i parametri α1 ... α4 e β1 ... β4 ,chiamati coordinate generalizzate, in modo che l’espressione (4.32) rappresenti un’ interpolazione fra i valori che lo spostamento assume nei nodi dell’elemento. A questo scopo, indichiamo con A una matrice che raccoglie i valori delle funzioni base calcolate nei nodi dell’elemento: 1 1 1 1 0 0 0 0 1 −1 1 −1 0 0 0 0 1 −1 −1 1 0 0 0 0 1 1 −1 −1 0 0 0 0 A = 0 0 0 0 1 1 1 1 0 0 0 0 1 −1 1 −1 0 0 0 0 1 −1 −1 1 0 0 0 0 1 1 −1 −1 u = Φα
Con tale posizione, possiamo quindi scrivere ˆ = A α u
(4.33)
La relazione di interpolazione si esprime in generale come: ˆ u = Hu
(4.34)
e quindi, dal confronto con la (4.33) si ricava H = Φ A −1
(4.35)
La determinazione delle funzioni di interpolazione, contenute nella matrice H si ottiene previa valutazione della matrice A −1 e fornisce l’espressione seguente: 1 (1 + x )(1 + y ) (1 − x )(1 + y ) (1 − x )(1 − y ) (1 + x )(1 − y ) H= 0 0 0 0 4 (4.36) 0 0 0 0 (1 + x )(1 + y ) (1 − x )(1 + y ) (1 − x )(1 − y ) (1 + x )(1 − y )
Nella Figura 4.7 sono riportati i grafici delle componenti di H, cioè le funzioni di interpolazione dell’elemento finito quadrilatero a quattro nodi.
4.5 CARATTERISTICHE DELLA SOLUZIONE
2
79
y 1
3 4
h 11
x
h 12
h 13
h 14
Figura 4.7: Grafici delle funzioni di interpolazione per l’elemento quadrilatero a quattro nodi (Esempio 4.2).
4.5 Caratteristiche della soluzione Dopo aver descritto il metodo di discretizzazione e la formulazione del problema che, come abbiamo visto porta alla determinazione del vettore degli spostamenti nodali, è opportuno riflettere sulla qualità del risultato ottenuto. Un primo aspetto da considerare è la verifica della soluzione dal punto di vista delle condizioni di equilibrio. Ci domandiamo pertanto se, un volta determinati gli spostamenti dei nodi, sussista l’equilibrio del corpo discretizzato. Consideriamo, in particolare, un punto appartenente a un dato elemento finito: l’equilibrio del suo intorno infinitesimo richiede, come abbiamo ricordato, che risultino soddisfatte le equazioni statiche. T ∂ σ+b = 0
(4.37)
ˆ Esprimiamo ora lo stato tensionale in un punto a partire dalla soluzione U ottenuta risolvendo il sistema (4.20). Scegliendo, per semplicità, una relazione costitutiva elastico-lineare, se il punto appartiene all’elemento m siha,nelcaso di tensioni iniziali nulle: (m ) (m ) (m ) ˆ B U σ (4.38) =C Sostituendo la (4.38) nell’eq. (4.37), si vede che l’equilibrio nel punto sussiste solo nel caso particolare in cui C(m )
∂
T (m )
B
ˆU = −b(m )
(4.39)
Dato che il primo membro discende dalla scelta delle funzioni di interpolazione, mentre il secondo membro corrisponde alle forze di volume assegnate, è evidente che le equazioni indefinite d’equilibrio non sono, in generale, soddisfatte dalla soluzione ottenuta con il metodo degli elementi finiti.
4 IL METODO DEGLI ELEMENTI FINITI
80
Esempio 4.3
Accertato che la soluzione a elementi finiti non soddisfa le equazioni differenziali di equilibrio punto per punto, possiamo però mostrare che tale soluzione è tale che 1. in ciascun nodo la somma delle forze nodali calcolate sugli elementi ivi connessi uguaglia il carico nodale applicato; 2. le forze nodali calcolate per un elemento, a partire dallo stato tensionale costituiscono un sistema equilibrato. Nelle proposizioni precedenti il termine di forza nodale va inteso come forza equivalente, nel senso del principio dei lavori virtuali, allo stato tensionale presente nell’elemento finito. Definiamo cioè, analogamente a quanto visto nella definizione dei carichi nodali (v. ad es. l’eq. (4.16)), il vettore delle forze nodali per l’elemento m : F
(m )
=
Ω(m )
B(m ) T σ(m ) dΩ(m )
(4.40)
Esprimendo la tensione nell’elemento m in funzione degli spostamenti nodali, (m ) = C(m ) B(m ) ˆ σ U, e quindi (m )
F
=
(m ) T (m ) (m )
Ω(m )
B
C
B
dΩ
(m )
ˆU
(4.41)
La proposizione 1 si dimostra facilmente sommando le forze nodali provenienti da tutti gli elementi, e ricordando la definizione della matrice di rigidezza (4.13). Si ottiene perciò: n e
m
(m )
F
n e
=
m
(m ) T (m ) (m )
Ω(m )
B
C
B
(m )
dΩ
ˆU
(4.42)
Osservando che il secondo membro di questa equazione corrisponde, in base ˆ , e quindi, secondo la (4.20) al vettore R dei carichi alla (4.13) al prodotto K U applicati ai nodi, risulta dimostrata la tesi, cioè n e
F(m ) = R
(4.43)
m
La dimostrazione della proposizione 2 risulta immediata tramite un’applicazione del principio dei lavori virtuali. Come è noto, un sistema di forze applicato a un corpo risulta equilibrato se il lavoro virtuale compiuto da tali forze
4.5 CARATTERISTICHE DELLA SOLUZIONE
81
è nullo per qualunque atto di moto rigido del corpo. In questo caso, usando le consuete notazioni, dobbiamo quindi dimostrare che ¯ˆ T F(m ) = 0, ∀u¯ˆ ∈ {insieme dei moti rigidi di m } u
(4.44)
Ricordando ancora la definizione (4.40) si può scrivere ¯ T (m )
ˆ F u
=
¯ T (m ) T
Ω(m )
ˆ B u
(m )
σ
dΩ
(m )
=
Ω(m )
(B(m ) u¯ˆ )T σ(m ) dΩ(m )
(4.45)
Osserviamo ora che le deformazioni virtuali dell’elemento m si esprimono in funzione degli spostamenti nodali come ¯(m ) = B(m ) u¯ˆ ; l’equazione precedente diventa perciò ¯ˆ T F(m ) = ¯ (m ) T σ(m ) dΩ(m ) u (4.46)
Ω(m )
Se si considerano i soli atti di moto virtuale rigido dell’elemento m , le deformazioni corrispondenti devono essere nulle in tutti i punti dell’elemento stesso, cioè ¯ (m ) T = 0, e quindi risulta nullo l’integrale a secondo membro, consentendo di scrivere infine ¯ˆ T F(m ) = 0 u
(4.47)
che dimostra la tesi. Esempio 4.4
Verificare l’equilibrio della soluzione dell’Esempio 4.1. La risposta della struttura, illustrata nella Figura 4.3, all’istante t = 1 s, è descritta dal sistema seguente: 15,4 240 −13,0 E
13,0 U 2 13,0 U 3
−
31,0 111,3
=
La soluzione fornisce gli spostamenti nodali: U 2 =
14230 E
cm, U 3 =
16285 E
cm
e, ovviamente, U 1 = 0cm. RIcordando l’eq. (4.40), le forze nodali equivalenti, per i due elementi, si calcolano come: (1)
F
=
Ω(1)
(1) T (1)
B
σ
(2)
(2)
dΩ , F
=
Ω(2)
B(2) T σ(2) dΩ(2)
4 IL METODO DEGLI ELEMENTI FINITI
82
Le matrici contenenti le derivate delle funzioni d’interpolazione sono: B(1) =
−
1 100
B(2) = 0
−
1 100 1 80
1 80
0
Le tensioni negli elementi si possono calcolare facilmente mediante le deformazioni, che sono state assunte costanti all’interno di un elemento: σ
U 2 (1) (1) 2 = C = [E ] = 142,30N/cm
(2)
U 3 − U 2 (2) (2) 2 = C = [E ] = 25,69N/cm
(1)
σ
100
80 Sostituendo i valori numerici, usando l’ascissa locale x come variabile di integrazione, si calcolano le forze nodali: F(1) =
F(2) =
1 100 − 100 1 100
0
80
0
142,30 142,30 · 1 · dx = 142,30 0 0 0 0 x 2 1 − 80 25,69 1 + 40 dx = −111,32 1 111,32 80
−
I risultati sono illustrati in Figura 4.8(a). In Figura 4.8(b) e 4.8(c), invece, sono riportate rispettivamente le forze nodali equivalenti ai carichi distribuiti b (1) e b (2) , e i carichi concentrati applicati nei nodi (v. Esempio 4.1). Si noti che, fra i carichi concentrati, all’estremo sinistro è stata indicata la reazione vincolare R 1b , ancora incognita. L’applicazione dell’eq. 4.43 ai due elementi fornisce i seguenti valori: 142,30 142,30 0
2
m =1
F(m ) =
−
+
0 −111,32 111,32
142,30 30,98 111,32
−
=
R 1b R 1b + 25,00 25,00 0,00 31,00 R = 25,00 + 6,00 + 0,00 = 0,00 11,33 100,00 111,33 Il confronto fra i due vettori mostra che, a meno di piccoli errori dovuti agli arrotondamenti, l’uguaglianza prevista teoricamente è soddisfatta. Inoltre, è possibile calcolare la reazione vincolare R 1b mediante la relazione R 1 = R 1a + R 1b , per cui: −142,30 = 25,00 + R 1b =⇒ R 1b = −167,30.
4.6 TRATTAMENTO LOCALE DEGLI ELEMENTI
-142,33
142,33
83
-111,37
111,37
6,00
11,33
0,00
100,00
(a) 25,00
25,00
(b) R1b
0,00
(c)
Figura 4.8: Calcolo delle forze nodali (Esempio 4.4): carichi equivalenti allo stato tensionale (a); carichi equivalenti alle forze distribuite (b); carichi concentrati e reazioni (c).
Esempio 4.5
Caso bidimensionale.
4.6 Trattamento locale degli elementi Nell’implementazione pratica del metodo degli elementi finiti, conviene costruire l’algoritmo che porta alla scrittura delle equazioni discretizzate d’equilibrio mediante operazioni standardizzate e ripetitive nel massimo grado possibile, per consentire all’elaboratore di dispiegare la sua massima efficacia. Un modo conveniente di procedere è operare con una procedura nella quale si prende in considerazione un elemento alla volta, con i suoi gradi di libertà, le sue funzioni di interpolazione e la sua risposta costitutiva. Per calcolare gli integrali necessari si può operare considerando: – una numerazione locale dei gradi di libertà – un sistema di riferimento locale. Nel seguito sono illustrati i concetti e la notazione necessari agli sviluppi della procedura di costruzione del sistema risolvente il problema discretizzato.
4.6.1 Numerazione locale dei gradi di libertà Consideriamo un dato elemento, relativamente al quale si vogliono valutare le matrici (di rigidezza, di massa, ecc.) e i vettori dei carichi equivalenti: ana-
4 IL METODO DEGLI ELEMENTI FINITI
84
logamente a quanto visto nel caso dei sistemi discreti, è sufficiente considerare dapprima solo i gradi di libertà dell’elemento di volta in volta esaminato, per poi assemblare le matrici globali sommando i contributi corrispondenti ai gradi di libertà globali coinvolti dall’elemento. L’espressione dell’approssimazione dello spostamento all’interno dell’elemento considerato si scrive: ˆ u(x , y , z ) = H(x , y , z )u
(4.48)
dove il vettore ˆu raccoglie gli spostamenti nodali dei soli nodi dell’elemento (si confronti con l’espressione (4.6)), espressi in un qualsiasi sistema di riferimento O x y z : ˆ = [u 1 v 1 w 1 ... u M v M w M ]T u avendo indicato con M il numero di nodi dell’elemento considerato. Il campo di deformazione approssimato si ottiene (si veda l’eq. (4.7) ) per derivazione del campo di spostamento (x , y , z ) = ∂u(x , y , z ) per cui, sfruttando l’eq.(4.48) si ha: ˆ = B(x , y , z )uˆ (x , y , z ) = ∂H(x , y , z )u (4.49) Per la valutazione di tutti gli integrali necessari si può quindi osservare che ciascuno degli integrali sotto il segno di sommatoria, nelle espressioni (4.13)(4.16) si può riscrivere omettendo l’esplicito riferimento all’indice m , ottenendo così le seguenti espressioni: matrice di rigidezza K =
BT CBdΩ
(4.50)
(4.51)
Ω
vettore delle forze di volume RB =
HT bdΩ
Ω
vettore delle forze di superficie RS =
Γ1
,...,Γq
HT f dΓ
(4.52)
vettore delle tensioni iniziali RI =
Ω
BT σI dΩ
(4.53)
4.6 TRATTAMENTO LOCALE DEGLI ELEMENTI
85
4.6.2 Sistemi di riferimento locali Gli integrali introdotti al paragrafo precedente si possono calcolare in un qualsiasi sistema di riferimento che risulti conveniente dal punto di vista pratico, ma, dovendo poi considerare una pluralità di elementi con riferimenti diversi, occorre l’accortezza di considerarne poi l’orientamento rispetto a un sistema di riferimento valido per l’intero corpo discretizzato, dettosistema di riferimento globale per distinguerlo dal sistema di riferimento locale valido per un dato elemento. Indichiamo con ˜u il vettore degli spostamenti nodali le cui componenti sono espresse rispetto al riferimento locale prescelto: ˜ (x , y , z )u˜ u=H
(4.54)
˜ (x , y , z ) deve essere valutata a partire da dove la matrice di interpolazione H quella consueta, H(x , y , z ). A questo fine, riprendiamo l’espressione analoga valida nel riferimento globale: ˆ u = H(x , y , z )u
(4.55)
e supponiamo di conoscere la matrice di rotazione T che lega le componenti di spostamento espresse nei due sistemi: ˜ = Tuˆ u
(4.56)
Sostituendo l’eq. (4.56) nella (4.54) e uguagliando quest’ultima alla (4.55) si trova l’espressione cercata: ˜ (x , y , z )T H(x , y , z ) = H
(4.57)
e, analogamente, si ricava la relazione fra le matrici delle derivate delle funzioni d’interpolazione: B(x , y , z ) = ˜B(x , y , z )T (4.58) Nella pratica, tuttavia, non conviene valutare esplicitamente le funzioni che ˜ e B˜ ; risulta infatti più efficiente trasformare gli incompaiono nelle matrici H tegrali (4.50)-(4.53) tenendo conto delle relazioni (4.57) e (4.58). Ad esempio, la matrice di rigidezza di un elemento si trasforma nel modo seguente: K =
BT CBdΩ
Ω
=
˜ T CBT ˜ dΩ TT B
Ω
T
=T
˜ T CB˜ TdΩ B
Ω
T ˜
= T KT
(4.59)
4 IL METODO DEGLI ELEMENTI FINITI
86
˜ = Ω ˜BT CB˜ dΩ relativa alle componenti di per cui si può valutare la matrice K spostamento ruotate e utilizzare poi l’espressione (4.59) per trasformarla nella consueta matrice K . In modo analogo, si dimostra che la matrice di massa si trasforma secondo la relazione: ˜ M = TT MT (4.60) I vettori dei carichi si valutano anch’essi mediante una procedura analoga a quella seguita per la matrice di rigidezza. Ad esempio, per le forze nodali equivalenti alle forze di volume la trasformazione si formula nel modo seguente:
RB = =
HT bdΩ
Ω
˜ T bdΩ TT H
Ω
T
=T
Ω
= T R˜B T
˜ T
H bdΩ
(4.61)
Esempio 4.6
Si consideri un elemento monodimensionale a due nodi (asta) (v. figura 4.9): le componenti di spostamento all’interno dell’elemento si possono esprimere fa˜ (x , y , z )u˜ ,dove u = [u (x ), v (x ))]T cilmente nel sistema di riferimento locale: u = H u 1 , ˜ v 1 , ˜ u 2 , ˜ v 2 ]T . e u˜ = [ ˜ Come è noto, nel caso bidimensionale, le componenti del vettore spostamento si trasformano, per una rotazione α del sistema di riferimento secondo la legge ˜1 ˆ1 u cos α sen α u (4.62) = ˜1 ˆ1 v v − sen α cos α Considerando entrambi i nodi dell’elemento si può quindi scrivere:
˜1 u ˜1 v = ˜2 u ˜2 v
cos α sen α − sen α cos α 0 0 0 0
per cui, la matrice di trasformazione cercata risulta
T=
cos α sen α − sen α cos α 0 0 0 0
0 0 0 0 cos α sen α − sen α cos α
0 0 0 0 cos α sen α − sen α cos α
ˆ1 u ˆ1 v ˆ2 u ˆ2 v
(4.63)
(4.64)
4.6 TRATTAMENTO LOCALE DEGLI ELEMENTI
87
Figura 4.9: Sistemi di riferimento globale e locale per un elemento monodimensionale (Esempio 4.6).
La tabella 4.1 raccoglie le diverse notazioni usate per indicare, nei diversi contesti, i vettori degli spostamenti nodali. Si noti che, per alleggerire un po’ la notazione, d’ora in poi il vettore delle componenti di spostamento nodale per ˆ. l’intera struttura si indicherà con U anziché con U u(m ) = H(m ) ˆ U ( ( m ) (o u = H m ) U) ˆ u = Hu ˜ u˜ u=H
u(m ) : spostamento entro l’elemento m ˆ , U: spostamenti nodali (num. globale) U ˆ : spostamenti nodi elemento (sistema globale) u ˜ : spostamenti nodi elemento (sistema locale) u
Tabella 4.1: Notazioni usate per i vettori spostamento.
Esempio 4.7
La struttura illustrata in figura 4.10 viene discretizzata considerandola composta di un elemento quadrilatero a quattro nodi in condizioni di tensione piana (A), di due aste a due nodi (B e C) e di una trave a due nodi (D), e vincolata come illustrato. Si determini la struttura della matrice di rigidezza a partire dalle matrici dei singoli elementi, tenendo conto delle loro mutue connessioni. È necessario dapprima individuare i gradi di libertà effettivi (cioè quelli non vincolati) del sistema, e numerarli per assegnare loro una posizione nel vettore U degli spostamenti globali della struttura. Tale numerazione è arbitraria, anche se influisce sulla srtuttura della matrice K . Si noti che l’estremo sinistro
88
4 IL METODO DEGLI ELEMENTI FINITI
Figura 4.10: Assemblaggio degli elementi (Esempio 4.7).
dell’elemento D ha tre gradi di libertà, comprendendo, oltre alle componenti di spostamento verticale e orizzontale, la rotazione (v. la figura 4.11).
Figura 4.11: Elencazione dei gradi di libertà (Esempio 4.7).
I vettori spostamento e le matrici di rigidezza degli elementi siano simboli-
4.6 TRATTAMENTO LOCALE DEGLI ELEMENTI
89
camente espressi, nelle rispettive numerazioni locali, nel modo seguente:
u 1 v 1 u 2 v ˆ A = 2 , u u 3 v 3 u 4 v 4
a 11 a 21 a 31 a 41 K A = a 51 a 61 a 71 a 81
a 12 a 22 a 32 a 42 a 52 a 62 a 72 a 82
a 13 a 23 a 33 a 43 a 53 a 63 a 73 a 83
a 14 a 24 a 34 a 44 a 54 a 64 a 74 a 84
a 15 a 25 a 35 a 45 a 55 a 65 a 75 a 85
a 16 a 26 a 36 a 46 a 56 a 66 a 76 a 86
a 17 a 27 a 37 a 47 a 57 a 67 a 77 a 87
a 18 a 28 a 38 a 48 a 58 a 68 a 78 a 88
(4.65)
u 1 v ˆ B = 1 , u u 2 v 2
b 11 b 21 K B = b 31 b 41
b 12 b 22 b 32 b 42
b 13 b 23 b 33 b 43
b 14 b 24 b 34 b 44
(4.66)
u 1 v ˆ C = 1 , u u 2 v 2
c 11 c 21 K C = c 31 c 41
c 12 c 22 c 32 c 42
c 13 c 23 c 33 c 43
c 14 c 24 c 34 c 44
(4.67)
d 13 d 23 d 33 d 43 d 53 d 63
d 14 d 24 d 34 d 44 d 55 d 66
d 15 d 25 d 35 d 45 d 55 d 65
u 1 v 1 θ ˆ D = 1 , u u 2 v 2 θ2
d 11 d 21 d 31 K D = d 41 d 51 d 61
d 12 d 22 d 32 d 42 d 52 d 62
d 16 d 26 d 36 d 46 d 56 d 66
(4.68)
Per l’assemblaggio dei termini di rigidezza è necessario istituire una corrispondenza tra la numerazione locale e quella globale dei gradi di libertà. È immediato verificare che sussistono le corrispondenze seguenti: elemento A:
componente di spost. u 1 v 1 u 2 v 2 u 3 v 3 u 4 v 4 ˆ 1 u ˆ2 u ˆ 3 u ˆ4 u ˆ5 u ˆ6 u ˆ7 u ˆ8 numerazione locale u numerazione globale U 2 U 3 − − − U 1 U 4 U 5 elemento B:
componente di spost. u 1 v 1 u 2 v 2 ˆ1 u ˆ2 u ˆ 3 u ˆ4 numerazione locale u numerazione globale U 6 U 7 U 4 U 5
(4.69)
(4.70)
4 IL METODO DEGLI ELEMENTI FINITI
90
elemento C:
componente di spost. u 1 v 1 u 2 v 2 ˆ1 u ˆ2 u ˆ 3 u ˆ4 numerazione locale u numerazione globale U 6 U 7 U 2 U 3
(4.71)
elemento D:
componente di spost. u 1 v 1 θ1 u 2 v 2 θ2 ˆ1 u ˆ2 u ˆ3 u ˆ4 u ˆ5 u ˆ6 numerazione locale u numerazione globale − − − U 6 U 7 U 8
(4.72)
Le matrici elementari vanno ora espanse alla dimensione della matrice di rigidezza globale, in questo caso 8 × 8, prima di procedere alla somma K = n e (m ) , e quindi, prendendo ad esempio l’elemento A, si avrà: m K
U 1 U 2 U 3 U 4 U= , U 5 U 6 U 7 U 8
a 66 a 16 a 26 a 76 ( A ) K = a 86
0 0 0 0
a 61 a 11 a 21 a 71 a 81
a 62 a 12 a 22 a 72 a 82
a 67 a 17 a 27 a 77 a 87
a 68 a 18 a 28 a 78 a 88
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
(4.73)
Procedendo allo stesso modo per gli altri elementi e sommando poi fra loro le matrici espanse, si ottiene la matrice di rigidezza dell’intera struttura: a 66 a 61 a 62 a 67 a 68 0 0 0 a 16 a 11+c 33 a 12+c 34 a 17 a 18 c 31 c 32 0 a 26 a 21+c 43 a 22+c 44 a 27 a 28 c 41 c 42 0 a 76 a 71 a 72 a 77+b 33 a 78+b 34 b 31 b 32 0 K = a 86 a 81 a 82 a 87+b 43 a 88+b 44 b 41 b 42 0 c 13 c 14 b 13 b 14 b 11+c 11+d 44 b 12+c 12+d 45 d 46 0 c 23 c 24 b 23 b 24 b 21+c 21+d 54 b 22+c 22+d 55 d 56 0 d 64 d 65 d 66 0 0 0 0 0
(4.74)
4.7 Imposizione delle condizioni al contorno L’ equazione di equilibrio dinamico per un sistema discretizzato libero da vincoli è ¨ + KU = R MU (4.75)
4.7 IMPOSIZIONE DELLE CONDIZIONI AL CONTORNO
91
Raccogliendo in U a gli spostamenti incogniti e in U b quelli noti (imposti), si può scrivere
¨ a U K aa K ab + ¨ Ub K ba K bb
Maa Mab Mba Mbb
Ua Ra = Ub Rb
(4.76)
Si ricava Ua risolvendo l’equazione Maa ¨ Ua + K aa Ua = Ra − Mab ¨ Ub − K ab Ub
(4.77)
Si vede quindi che, per determinare Ua basta assemblare le matrici di rigidezza e massa pertinenti ai soli gradi di libertà effettivamente incogniti. Occorre però “correggere” il termine noto nel caso di spostamenti imposti non nulli. Esso diventa: Ra − Mab ¨ Ub − K ab Ub (4.78) La soluzione del sistema (4.77) fornisce gli spostamenti incogniti U e le ac¨ , dopo di che si possono calcolare le forze nodali nei nodi vincocelerazioni U lati, mediante la seconda riga del sistema partizionato (4.76): Rb = Mba ¨ Ua + Mbb ¨ Ub + K ba Ua + K bb Ub
(4.79)
Tali forze nodali dovranno risultare pari a: b b b Rb = Rb B + RS + RI + RC + Rr
(4.80)
b dove R b B , R b S , R b I sono i carichi nodali equivalent, R C i carichi concentrati applicati nei nodi e R r sono le forze che nascono nei nodi vincolati allo scopo di mantenere l’equilibrio, vale a dire le reazioni vincolari , che risultano quindi:
b b b Rr = Rb − Rb B + RS + RI + RC
b b b Ua + Mbb ¨ Ub + K ba Ua + K bb Ub − Rb = Mba ¨ B + RS + RI + RC
(4.81)
Da quanto appena detto emerge il fatto che la soluzione primaria del problema, cioè la determinazione degli spostamenti dei nodi liberi, può essere ottenuta senza assemblare l’intera matrice di rigidezza della struttura, omettendo sin dall’inizio quei gradi di libertà, compresi nel vettore Ub , che a tutti gli effetti non sono da considerarsi incogniti, ma assegnati (ad esempio pari a zero, nel caso di vincoli perfetti).
4 IL METODO DEGLI ELEMENTI FINITI
92
4.8 Relazioni costitutive elastiche Finoraabbiamo affrontato il problema tensio-deformativo perun generico corpo tridimensionale, sotto le più generali ipotesi di sollecitazione e di vincoli imposti agli spostamenti. Non sempre, però, è necessario e utile affrontare il problema tramite questo modello strutturale. Mentre in linea di principio risulta sempre possibile l’approccio tridimensionale, ogniqualvolta sia possibile ridurre la complessità dello studio, tenendo conto di informazioni supplementari a disposizione dell’analista, si cerca di adottare gli schemi più semplici nella predisposizione del modello a elementi finiti. Se assumiamo, come di consueto, che il materiale sia elastico-lineare isotropo, descritto ad esempio dal modulo elastico longitudinale E e dal rapporto di Poisson ν, la relazione fra tensioni e deformazioni si può esprimere in forma matriciale come
1
ν σxx 1− νν σ y y E (1 − ν) 1−ν σzz = τx y (1 + ν)(1 − 2ν) 0 τ y z 0 τxz
0
ν 1−ν
1 ν 1−ν
ν 1− νν
1−ν 1
0
0
0
0
0
0
0
0
0 1−2ν 2(1−ν) 0
0
0
0
0
0
0 1−2ν 2(1−ν) 0
0
0
0
0 1−2ν 2(1−ν)
xx y y zz γx y γ y z γxz
(4.82) Vedremo ora brevemente alcune situazioni nelle quali la dimensionalità del problema si può ridurre. In alcuni casi ciò non comporta significative variazioni della procedura di analisi rispetto a quanto visto nella precedente trattazione; in altri casi, invece, è necessario il ricorso a modelli più complessi, che richiedono lo sviluppo preliminare di appropriate teorie.
4.8.1 Stato di tensione monoassiale Le condizioni di sollecitazione più semplici sono quelle associate a uno stato tensionale monoassiale . Questo tipo di situazione è tipico degli elementi monodimensionali che trasmettono solo lo sforzo normale, lungo l’asse, come le aste o bielle che costituiscono una travatura reticolare caricata nei nodi (v. figura 4.12). Con riferimento alla figura, detto x l’asse di un elemento, lo stato tensionale è quindi completamente caratterizzato dal valore σ xx , uniforme sulla sezione retta dell’elemento.
4.8 RELAZIONI COSTITUTIVE ELASTICHE
93
Figura 4.12: Esempio di struttura composta da elementi strutturali in condizioni di tensione monoassiale. Per questo elemento strutturale, la deformazione è completamente definita, tramite la componente assiale u dello spostamento, dalla dilatazione assiale xx . La relazione costitutiva, espressa in forma di prodotto matriciale, è dunque: σxx = E xx (4.83)
4.8.2 Stato di tensione piana Nel caso di lastre caricate nel loro piano, o di travi inflesse, lo stato tensionale è al più biassiale (figura 4.13). Se l’asse z è perpendicolare al piano della struttura e dei carichi, nel caso in esame risultano nulle le componenti di tensione σzz , τzx e τz y . Adottando la teoria tecnica della trave, la deformazione dell’elemento è descritta dalla curvatura κxx , duale del momento flettente M x x . La nota relazione momento-curvatura si scrive, in forma matriciale:
M x x = E I κxx
(4.84)
dove con I si indica il momento d’inerzia della sezione retta dell’elemento. Per quanto riguarda le lastre in stato piano di tensione, considerando le componenti di tensione non nulle e le corrispondenti componenti di deformazione, la relazione tensione-deformazione è data da:
1 ν ν 1
σxx E σ y y = (1 − ν)2 τx y 0 0
0 0
1−ν 2
xx y y γx y
(4.85)
4.8.3 Stato di deformazione piana Si consideri una struttura come quella schematizzata nella figura 4.14: si tratta di una diga, realizzata con un materiale isotropo, ottenuta dalla traslazione
94
4 IL METODO DEGLI ELEMENTI FINITI
Figura 4.13: Esempi di strutture in condizioni di tensione piana (travi e lastre). nella direzione z di una figura piana. Se la lunghezza nella direzione z risulta molto grande rispetto alle dimensioni trasversali della sezione, data la configurazione di carico assegnata (la spinta idrostatica e il peso proprio agente lungo l’asse y ), per simmetria, i punti appartenenti a una sezione non possono subire spostamenti fuori dal suo piano: si parla quindi di stato di deformazione piana . Un problema di questo tipo può dunque essere studiato considerando una sezione della struttura, opportunamente vincolata, per i punti della quale il vettore spostamento ha due sole componenti non nulle, u e v , essendo, come abbiamo detto w = 0 ovunque. Si noti che, in queste condizioni, anche le componenti di deformazione zz , γ zz e γ zx sono nulle. Essendo il materiale elastico-lineare e isotropo, è facile mostrare come il piano della sezione sia un piano principale delle tensioni sul quale, come è noto, non agiscono com-
Figura 4.14: Esempio di struttura in condizioni di deformazione piana.
4.8 RELAZIONI COSTITUTIVE ELASTICHE
95
ponenti di tensione tangenziale (quindi τ z x = τ z y = 0), ma risulta in generale diversa da zero la componente di tensione normale σz z . La relaz relazion ionee fra fra tensio tensioni ni e deform deformazi azioni oni è quind quindii espre espressa ssanel nel modo modo seguen seguen-te
1
σx x ν E (1 (1 − ν) σ y y = (1 + ν)(1 − 2ν) 1 − ν τx y
0
ν
1−ν 1 0
0 0 1 − 2ν 2(1 − ν)
x x y y γx y
(4.86)
E la componente di tensione perpendicolare è legata alle altre tramite la relazione
σz z = ν σx x + σ y y
4.8.4 4.8.4 Condiz Condizion ionii di assial assialsim simmet metria ria Un altr altro o caso caso in cui cui si può può ridur ridurre re la dime dimens nsion ional alit itàà del del prob proble lema ma è quel quello lo dett detto o assialsimmetrico . Si tratta tratta di strutture strutture con la forma del solido di rivoluzione rivoluzione intorno intorno a un asse. Natura Naturalment lmente, e, anche anche i carichi carichi devono rispettare rispettare la stessa stessa simmetria (v. figura 4.15). In questo caso, i punti appartenenti alla sezione che genera il solido di rivoluzione possono spostarsi solo nel piano della sezione e quindi, se z è è l’asse perpendicolare a tale piano, si avrà w = = 0. Per quanto riguarda lo stato tensionale, risultano nulle le componenti τ y z e
Figura 4.15: Stato assialsimmetrico. assialsimmetrico.
4 IL METODO METODO DEGLI ELEMENTI ELEMENTI FINITI FINITI
96
τx z , mentre le altre sono, in generale, diverse da zero. Risulta quindi:
ν
1
σx x ν E (1 − ν) σ y y 1−ν = τx y (1 + ν)(1 − 2ν) 0 σz z ν 1−ν
1−ν 1 0 ν
1−ν
0
ν
1− νν
0 1−ν 1 − 2ν 0 2(1 − ν) 0 1
x x y y γx y z z
(4.87)
4.8. 4.8.5 5 Lastr Lastre e infle infless sse e In quest questaa catego categoria ria ricado ricadono no elemen elementi ti con una dimens dimension ionee nettam nettament entee minor minoree rispet rispetto to alle alle altre altre due, due, carica caricate te trasv trasvers ersalm alment entee e, perciò perciò inflesse . Si parla parla usualusualmente di piastre , nel caso di lastre piane (figura 4.16), e di gusci , nel caso di lastre a semplice o doppia curvatura (figura 4.17). Per quanto riguarda le piastre, lo stato tensionale è piano, cioè la componente σ z z è nulla, ma, a differenza delle lastre caricate nel piano, tutte le altre componenti di tensione sono, in generale, diverse da zero. Anche in questo caso, analogamente a quello delle travi inflesse, i nflesse, le variabili duali atte a descrivere la relazione costitutiva sono i momenti flettenti e il momento torcente, torcente, e le curvature e la torsione.
1 ν ν 1
M x x E h 3 M y y y = 12(1 − ν)2 M x y 0 0
0 0
1−ν 2
κx x κ y y = κx y
(4.88)
Le curvature e la torsione si possono esprimere in funzione dello spostamento w (x , y ) tramite le relazioni: ∂2 w κx x = , ∂x 2
∂2 w κ y y = , ∂ y 2
∂2 w κx y = ∂x ∂ y
(4.89)
Esercizi 4.1. Ricavare la legge di trasformazione (4.60) della matrice di massa e quel-
la del vettore dei carichi di superficie relative a una rotazione del sistema di riferimento. 4.2. Ripete Ripetere re l’esem l’esempio pio 4.7, 4.7, sostit sostituen uendo do il carre carrello llo che che vincol vincolaa l’elem l’element ento o A con
una cerniera fissa.
4.8 RELAZIONI COSTITUTIVE ELASTICHE
Figura 4.16: Stato tensionale in una piastra.
Figura 4.17: Stato tensionale in un guscio.
97
Capitolo 5 Formulazione degli elementi finiti 5.1 Introduzione Abbiamo visto, nel capitolo precedente, che il metodo degli elementi finiti è una procedura di discretizzazione delle strutture, che permette di affrontare la soluzione delle equazioni differenziali che governano la deformazione dei corpi, valida per qualsiasi forma strutturale e per qualunque condizione di vincolo e di carico. In questo capitolo, soffermeremo l’attenzione su un aspetto fondamentale del metodo, che ne costituisce in certo qual modo il cuore: la descrizione delle leggi di approssimazione del campo di spostamento all’interno degli elementi. Tale scelta è condizionata, come vedremo, da vari fattori, ed è nostro scopo motivarla e renderla consapevole, anche se limiteremo l’analisi ad alcuni casi, per quanto di rilevante interesse applicativo. La scelta delle funzioni di interpolazione dell’elemento, associata alla definizione della legge costitutiva del materiale che lo costituisce, porta immediatamente alla determinazione della sua matrice di rigidezza. Vedremo quindi come costruire tale matrice, e come i suoi termini possano efficacemente essere valutati per via numerica, adatta all’impiego delle procedure di calcolo automatizzate.
5.2 Elementi isoparametrici L’obiettivo fondamentale della formulazione degli elementi finiti isoparametrici è ottenere la relazione fra spostamenti nodali e spostamenti interni direttamente mediante le funzioni d’interpolazione, scelte senza passare attraverso la tecnica delle coordinate generalizzate vista al capitolo precedente. 98
5.2 ELEMENTI ISOPARAMETRICI
99
L’idea di fondo della trattazione degli elementi isoparametrici è l’introduzione delle cosiddette coordinate naturali e l’applicazione dell’integrazione numerica su un elemento tipo, rappresentante di tutta una classe di elementi aventi geometria differente. A questo scopo, si introduce il concetto di interpolazione della geometria dell’elemento, che porta alla possibilità di istituire facilmente una corrispondenza biunivoca, detta mapping , fra i punti dell’elemento tipo (parent element ) e quelli dell’elemento reale a cui si fa di volta in volta riferimento (figura 5.1).
Figura 5.1: Corrispondenza (mapping) fra l’elemento tipo, con il sistema di coordinate naturali, e l’elemento generico. Il concetto di mapping e la sua efficacia si dimostra facilmente con un semplice esempio. Esempio 5.1 Si consideri un elemento di biella dotato di rigidezza assiale E A e lunghezza L : determinare la matrice di rigidezza.
Figura 5.2: Elemento isoparametrico a due nodi (Esempio 5.1). Si introduce il sistema di riferimento in coordinate naturali r ,
1 ≤ r ≤ 1,
−
5 FORMULAZIONE DEGLI ELEMENTI FINITI
100
legate alle coordinate globali dalla trasformazione lineare 1 1 X (r ) = (1 − r ) X 1 + (1 + r ) X 2 2 2
(5.1)
dove X 1 e X 2 sono le coordinate nodali dell’elemento nel sistema globale. La trasformazione di coordinate è un’interpolazione fra queste coordinate, che associa a ciascun valore di r un punto dell’elemento compreso fra i due nodi estremi. La relazione tra r e X sopra introdotta si può scrivere 2
X (r ) =
h (r ) X i
i
(5.2)
i =1
dove si sono introdotte le funzioni di forma 1 h 1 (r ) = (1 − r ), 2
1 h 2 (r ) = (1 + r ) 2
(5.3)
La formulazione isoparametrica deve il suo nome al fatto che gli spostamenti sono interpolati mediante le stesse funzioni usate nell’interpolazione delle coordinate, cioè le funzioni di interpolazione coincidono con le funzioni di forma, cioè si assume: 2
U (r ) =
h (r )U i
i
(5.4)
i =1
Ricordiamo che la deformazione nell’asta è data da = dU /d X , ma, dato che è nota la dipendenza di U da r , occorre esprimere tale deformazione tramite la regola di derivazione delle funzioni composte: =
dU dU dr = d X dr d X
(5.5)
Sfruttando le definizioni, è immediato ottenere le derivate: 2 dh (r ) dU i U i = dr i =1 dr
(5.6)
2 d X = dr i =1
(5.7)
dh (r ) X i
dr
i
Svolgendo i calcoli si ricava: dU U 2 − U 1 = dr 2 d X X 2 − X 1 L = = dr 2 2
(5.8) (5.9)
5.3 ELEMENTI DEL CONTINUO N -DIMENSIONALE
101
Ricordando la regola di derivazione della funzione inversa, dr /d X = 1/(d X /dr ), si ottiene l’espressione della deformazione:
=
U 2 − U 1 2 L 2
dU dr = d X d X
=
U 2 − U 1 L
(5.10)
Questa espressione si può mettere in forma matriciale per ritrovare la scrittura generale = BU nel modo seguente: 1 = BU = L
1 1 U 1
−
U 2
(5.11)
da cui ricaviamo appunto che B=
1 L
1 1
−
(5.12)
A questo punto è possibile valutare la matrice di rigidezza, ricordando l’espressione della matrice costitutiva appropriata, C = [E ]: X 2
K =
X 1
E A BT CB A dx = 2 L
1
1
−
1
−
1 1 J dr −
1
(5.13)
dove, essendo l’integrale valutato sul dominio normalizzato, è necessario introdurre il determinante jacobiano della trasformazione di coordinate: J =
d X L = dr 2
(5.14)
Completando i calcoli si ottiene, ovviamente, la ben nota matrice di rigidezza dell’asta: E A 1 −1 (5.15) L −1 1
5.3 Elementi del continuo n -dimensionale Soffermiamoora l’attenzione su un’ampia classe di elementi finiti, che si possono considerare porzioni di un corpo continuo. In generale l’analisi si condurrà nel contesto tridimensionale, associando a ciascun nodo i tre gradi di libertà u , v e w , cioè le componenti di spostamento nelle direzioni dei tre assi coordinati x , y , z . Nel caso in cui sia possibile ridurre la dimensionalità del problema, ad
5 FORMULAZIONE DEGLI ELEMENTI FINITI
102
esempio come nei casi di tensione o deformazione piana o di assialsimmetria, la trattazione resterà valida considerando solo gli assi coordinati necessari alla situazione particolare. Per un generico elemento tridimensionale dotatodi q nodi, l’interpolazione delle coordinate, cioè il mapping rispetto alle coordinate naturali r , s , t si scrive, in coordinate locali x , y , z : q
x =
q
h (r , s , t )x , i
i =1
i
y =
q
h (r , s , t ) y , i
i =1
i
z =
h (r , s , t )z i
i
(5.16)
i =1
dove x i , y i , z i , i = 1,...q sono le coordinate dei nodi dell’elemento considerato. Ricordiamo che le funzioni d’interpolazione h i sono definite nel sistema di coordinate naturali dell’elemento r , s , t , ciascuna delle quali varia tra −1 e 1 e, come anticipato, per elementi mono- o bidimensionali si usano solo le espressioni in r , o r e s , rispettivamente. Una funzione di interpolazione h i deve avere valore unitario nel nodo i e nullo negli altri nodi; sfruttando queste condizioni, le h i si potrebbero determinare in modo sistematico, una volta nota la disposizione dei nodi. Conviene però costruirle per successive addizioni e correzioni come nell’esempio seguente. Esempio 5.2 Costruire le funzioni di interpolazione per l’elemento asta a tre nodi illustrato nella figura 5.3.
Figura 5.3: Elemento monodimensionale a tre nodi (Esempio 5.2). Si consideri dapprima un’interpolazione lineare, valida per l’elemento a due nodi (figura 5.4). Sono le stesse utilizzate nell’esempio 5.1: 1 h 1 (r ) = (1 − r ), 2
1 h 2 (r ) = (1 + r ) 2
(5.17)
5.3 ELEMENTI DEL CONTINUO CONTINUO N -DIMENSIONALE
103
Figura 5.4: Funzioni di forma per l’elemento a due nodi (Esempio 5.2).
La prese presenza nza di un terzo terzo nodo nodo consen consente te un un’in ’inter terpol polazi azione one con un polinom polinomio io di seco second ndo o grad grado o. Si cost costru ruis isce ce dapp dappri rima ma la funz funzio ione ne di form forma a per per il no nodo do inte interrmedio, che sarà appunto una parabola con valore unitario in corrispondenza del nodo 3 e 3 e nullo negli altri due nodi (v. figura 5.5), cioè h 3 (r ) r ) = (1 − r )(1 r )(1 + r ) r )
(5.18)
Figura 5.5: Funzione di forma per il nodo centrale (Esempio 5.2). Le funzioni di forma lineari, valide per l’elemento a due nodi, non soddisfano la condizione di assumere valore unitario in un nodo e nullo negli altri: infatti presentano entrambe un valore pari a 1/2 nel 1/2 nel nodo centrale. Risulta pertanto necessario ripristinare ripristinare tale condizione. Le funzio funzioni ni h e h 2 si cost costru ruisc iscon ono o come come somma somma dell della a part parte e line linear are e e di un una a h 1 e h “correzione” quadratica, nulla nei nodi 1 e 1 e 2 e 2 e valore −1/2 nel 1/2 nel nodo 3 (figura 3 (figura 5.6). Esse assumono le espressioni 1 1 h 1 (r ) r ) = (1 − r ) r ) − (1 − r 2 ), 2 2
1 1 h 2 (r ) r ) = (1 + r ) r ) − (1 − r 2 ) 2 2
(5.19)
5 FORMULAZIONE FORMULAZIONE DEGLI ELEMENTI FINITI
104
Figura 5.6: Funzioni di forma per l’elemento a tre nodi (Esempio 5.2).
5.3.1 Matrici H e B Tornando a considerazioni generali valide per tutti gli elementi isoparametrici, gli spostamenti sono interpolati nello stesso modo della geometria q
u =
q
h (r , s , t )t )u , i i
i i
v =
i =1
q
h (r , s , t )t )v , i i
w =
i i
i =1
h (r , s , t )t )w i i
i i
(5.20)
i =1
dove u dove u i i , v i i , w i i , i = 1,... q sono sono le componenti dello spostamento dei nodi dell’elemento. Le funzioni h funzioni h i i (r , s , t ) t ) sono gli elementi della matrice delle funzioni di interpolazione (o di forma) H forma) H.. La matrice delle derivate, B, serve, come sappiamo, per calcolare la matrice di rigidezza dell’elemento. dell’elemento. Essa lega le deformazioni agli spostamenti nodali. I suoi suoi termin terminii si otteng ottengon ono o deriva derivando ndo rispet rispetto to alle alle coor coordin dinate ate local localii x , y , y , z le le funfunzioni di forma. Sorge a questo punto un problema: problema: tali funzioni sono costruite, come come abbiam abbiamo o visto visto,, nello nello spazio spazio delle delle coor coordin dinate ate natur naturali ali r , s , t , t , e non è quindi immediato ottenere, ottenere, per derivazione, i termini della matrice B matrice B.. In generale generale,, per ottenere ottenere le derivate derivate richieste richieste,, sembra sembra necessario necessario introdur introdurre re le regole di derivazione di funzioni composte, composte, compendiate negli operatori: ∂ ∂ x ∂ ∂ y
=
=
∂ ∂ z
=
∂ ∂r ∂r ∂ x ∂ ∂r ∂r ∂ y ∂ ∂r ∂r ∂ z
+
+
+
∂ ∂ s ∂ s ∂ x ∂ ∂s ∂ s ∂ y ∂ ∂ s ∂ s ∂ z
+
+
+
∂ ∂ t ∂ t ∂ x ∂ ∂ t ∂ t ∂ y ∂ ∂ t ∂ t ∂ z
(5.21) (5.22) (5.23)
In queste espressioni compaiono però le derivate delle funzioni inverse: r = r ( r (x , y , z ), ), s = s (x , y , z ), ), t = t ( t (x , y , z ), ), che non sono note e comunque difficili da ricavare.
5.3 ELEMENTI DEL CONTINUO CONTINUO N -DIMENSIONALE
105
Sono invece immediatamente valutabili le derivate rispetto alle coordinate coordinate naturali, mediante gli operatori: ∂ ∂x
∂ ∂r
=
∂ x ∂r ∂ ∂x
∂ ∂s
=
∂ x ∂ s ∂ ∂x
∂ ∂ t
=
∂ x ∂ t
+
+
+
∂ ∂ y
+
∂ y ∂r ∂ ∂ y
+
∂ y ∂ s ∂ ∂ y
+
∂ y ∂t
∂ ∂z ∂ z ∂r ∂ ∂z ∂ z ∂ s ∂ ∂z ∂ z ∂ t
(5.24) (5.25) (5.26)
ovvero, in forma matriciale: ∂
r s ∂ ∂ ∂ ∂
=
∂ t
∂
∂
∂ y
∂z
∂ ∂
∂r ∂ y
∂ ∂
∂ ∂
∂s ∂ y
∂ ∂
∂ ∂
∂t
∂ t
∂ t
∂z
x x r s x
x r z y s z ∂ ∂
(5.27)
Introducendo la matrice la matrice jacobiana , definita da
J=
x x r x s ∂
∂ y
∂ z
∂ ∂
∂r ∂ y
∂ ∂
∂ ∂
∂ s ∂ y
∂ ∂
∂ t
∂ t
∂t
∂
∂
r z s z
(5.28)
l’espressione l’espressione (5.27), con le definizioni ∂ ∂r
∂
=
∂ ∂ x
r
=
∂ s
∂
∂
T
t
∂
∂
∂
∂ x
∂ y
∂ z
(5.29)
T
(5.30)
si può quindi porre in forma simbolica: ∂ ∂r
=
J
∂ ∂ x
(5.31)
Nel caso in cui la relazione fra coordinate naturali e coordinate locali sia biunivoca, esiste l’inverso dello jacobiano, jacobiano, J, J, cosicché si può scrivere: ∂ ∂ x
=
J−1
∂ ∂r
(5.32)
5 FORMULAZIONE DEGLI ELEMENTI FINITI
106
Per quanto concerne la matrice B = x H, occorre fare riferimento, caso per caso, all’operatore cinematico x valido nella situazione particolare che si considera. Infatti, tale operatore contiene espressioni che sono funzioni delle derivate parziali calcolate rispetto al sistema locale O x y z , ciascuna delle quali deve essere sostituita dall’opportuna espressione ricavata dai termini dell’eq. (5.32). Esempio 5.3 Derivare l’operatore jacobiano J e la matrice B per l’elemento asta a tre nodi descritto nell’esempio 5.2. La matrice H contiene le funzioni di forma calcolate in precedenza: H=
1 2
(r
−
1 (r + 1)r 2
1)r
(1 − r 2 )
Ricordando che B = x H e che, nel caso monodimensionale, d d , a sua volta esprimibile come J−1 dr , si ha d x B = J−1
−
1 2
+
1 2
r
+
r
−
x coincide
con
2r
Il calcolo dello jacobiano è immediato, a partire dall’espressione del map ping isoparametrico: x 1 x (r ) = H x 2 x 3
dove le ascisse nodali assumono i valori x 1 = 0, x 2 = L , x 3 = L /2. Dall’espressione precedente si ottiene quindi x (r ) = L (r + 1)/2 e dunque J=
dx L dr
=
2
J−1 =
,
2 L
,
detJ =
L 2
L’espressione della matrice delle derivate risulta perciò: B=
2 L
−
1 2
+
r
1 2
+
r
−
1 2r
2r
=
L
−
1 2r + 1
−
4r
5.3.2 Integrazione numerica Nel paragrafo precedente abbiamo visto che è possibile, almeno formalmente, una volta scelte le funzioni di forma da inserire nella matrice H, determinare la matrice B delle derivate, mediante l’espressione (5.28). In essa compare la matrice jacobiana inversa (eq. (5.32)), e quindi sembrerebbe necessario determinare le derivate delle inverse delle funzioni di forma, che sono gli elementi di
5.4 INTEGRAZIONE NUMERICA
113
5.4 Integrazione numerica Nel paragrafo precedente abbiamo visto che è possibile, almeno formalmente, una volta scelte le funzioni di forma da inserire nella matrice H, determinare la matrice B delle derivate, mediante l’espressione (5.28). In essa compare la matrice jacobiana inversa (eq. (5.32)), e quindi sembrerebbe necessario determinare le derivate delle inverse delle funzioni di forma, che sono gli elementi di tali matrici. Vedremo ora che, ricorrendo alla tecnica dell’integrazione numerica, l’esplicitazione analitica di tali funzioni non è necessaria per la valutazione dei termini di rigidezza e di carico dei vari elementi. Prendiamo ad esempio il calcolo di K : detto Ω il dominio di integrazione, cioè la regione dello spazio occupata dall’elemento, la matrice di rigidezza è, per definizione: K =
BT CBdΩ
(5.33)
Ω
Ricordando le regole relative al cambio di variabile di integrazione, se si adottano le coordinate naturali r , s , t , occorre esprimere l’elemento differenziale dΩ = dx d y dz in funzione delle nuove coordinate. Si ha d Ω = detJdr ds dt = detJdV , dove dV rappresenta l’elemento di volume per l’elemento tipo, cioè una regione che è un cubo nello spazio delle coordinate naturali, come mostrato nella figura 5.11 per il caso bidimensionale. Con la sostituzione indicata, l’integrale diventa. K =
BT CBdet JdV
(5.34)
V
e analogamente si procede per la valutazione di tutti gli integrali necessari alla scrittura delle equazioni di equilibrio, introdotti al Capitolo 4. In generale, possiamo compendiare l’intera casistica degli integrali da valutare introducendo una funzione F (r , s , t ) per la quale si voglia calcolare il
Figura 5.11: Trasformazione dell’elemento di volume.
5 FORMULAZIONE DEGLI ELEMENTI FINITI
114
numero I =
1
F (r , s , t )dV =
V
1
1
F (r , s , t )dr ds dt
(5.35)
−1 −1 −1
Il valore di I può rappresentare un qualsiasi elemento dei vettori o delle matrici da valutare mediante integrazione. Gli integrali definiti del genere (5.35) si possono calcolare mediante diverse tecniche di approssimazione, che portano, alle cosiddette formule di quadratura . L’idea di formula di quadratura è semplice: si tenta di approssimare l’integrale con opportune sommatorie, del tipo: 1
I =
1
1
F (r , s , t )dr ds dt =
−1 −1 −1
α α α F (r , s , t ) i j k
i j
k
(5.36)
i , j ,k
cioè si determina la somma, pesata dai coefficienti αi , α j , αk , dei valori assunti dall’integrando F calcolati in certi punti del dominio normalizzato r i , s j , t k . Il numero di punti adottato, cioè il campo di valori di i , j e k dipende dalla precisione voluta. Senza entrare nella teoria matematica delle formule di quadratura, cui si fa cenno nel paragrafo 5.4.1, basti sapere che, nel nostro contesto, risultano particolarmente efficienti le formule di G AUS S . Tali formule richiedono la valutazione della funzione in punti del dominio individuati dall’annullamento di certi polinomi (di L EGENDRE , C HEBICHEV ecc.). Le posizioni di questi punti sono dette ascisse di G AUS S e i punti corrispondenti si chiamano punti di Gauss . Nella figura 5.12 sono rappresentati, ad esempio, i punti di Gauss per elementi bidimensionali quadrilateri. Nella tabella 5.5 sono riportate le ascisse dei punti di Gauss sull’intervallo ( 1,1), con i rispettivi pesi, per formule di quadratura di Gauss-Legendre da 1 a 6 punti.
−
Figura 5.12: Punti di Gauss per formule di quadratura di Gauss-Legendre con ordine compreso tra 2 e 4.
5.4 INTEGRAZIONE NUMERICA
115
n
r i
αi
1 2 3
0,000000000000000 ±0,577350269189626 ±0,774596669241483 0,000000000000000 ±0,861136311594053 ±0,339981043584856 ±0,906179845938664 ±0,538469310105683 0,000000000000000 ±0,932469514203152 ±0,661209386466265 ±0,238619186083197
2,000000000000000 1,000000000000000 0,555555555555556 0,888888888888889 0,347854845137454 0,652145154862546 0,236926885056189 0,478628670499366 0,568888888888889 0,171324492379170 0,360761573048139 0,467913934572691
4 5
6
Tabella 5.5: Ascisse e pesi di Gauss per n punti di integrazione.
5.4.1 Formule di quadratura Per il calcolo degli integrali definiti del tipo I =
b
f (x )dx
(5.37)
a
si può procedere per via numerica, approssimando la funzione f (x ) con un polinomio e integrando poi il polinomio stesso. Naturalmente, in generale, si commette un errore, in quanto l’integrale del polinomio approssimante non coincide necessariamente con quello della funzione f (x ). Per approssimare una funzione mediante un polinomio, può convenire costruire un’interpolazione mediante certipolinomi fondamentali , come, ad esempio, i polinomi di Lagrange. Scelti n punti x 1 , x 2 , . . . x n appartenenti all’intervallo di integrazione, tali polinomi, ovviamente di grado n 1, sono costruiti in modo da assumere valore unitario in un dato punto x i e nullo in tutti gli altri. La loro espressione generale è
−
l i (n −1) (x ) =
(x x 1 )(x x 2 ) · · · (x x i −1 )(x x i +1 ) · · · (x x n ) . (x i x 1 )(x i x 2 ) · · · (x i x i −1 )(x i x i +1 ) · · · (x i x n )
− −
− −
− −
− −
−
−
(5.38)
Esempio 5.4
Costruire i polinomi fondamentali di Lagrange per i punti x 1 = −1, x 2 = 0, x 3 = 1.
5 FORMULAZIONE DEGLI ELEMENTI FINITI
116
Essendo n = 3, il grado dei polinomi sarà n − 1 = 2. Applicando la formula (5.38) si ha: (x 0)(x 1) 1 = x (1 x ) ( 1 0)( 1 1) 2 (x + 1)(x 1) l 2(2)(x ) = x 2 ) = (1 (0 + 1)(0 1) (x + 1)(x 0) 1 l 3(2)(x ) = x (1 + x ) = (1 + 1)(1 0) 2 l 1(2)(x ) =
− − −− −− − − − − − − −
(5.39)
−
(5.40) (5.41)
Si sono ritrovati i ben noti polinomi usati come funzioni di forma per l’elemento finito a tre nodi. Un’approssimazione mediante un polinomio di grado n 1 della funzione f (x ) si può costruire tramite la combinazione lineare di n polinomi fondamentali, i cui coefficienti siano i valori assunti dalla funzione stessa nei punti x i :
−
n
(n 1)
(n 1) f (x i ) i
l −
P − (x ) =
i =1
(5.42)
Il polinomio P (n −1) (x ) prende il nome di polinomio interpolatore della funzione f (x ), in quanto i loro valori calcolati nei punti x i coincidono: P (n −1) (x i ) = f (x i ), i = 1,...n . Usando questa rappresentazione, l’integrale (5.37) diventa: b n
I =
a i =1
l i (n 1) f (x i )dx + E n ( f ) =
−
n
w f (x ) i
i + E n ( f )
(5.43)
i =1
dove si è messo in evidenza, indicandolo con E n ( f ), l’errore indotto dall’interpolazione su n punti e si sono introdotti i pesi w i =
b
a
l i (n −1)dx
(5.44)
I punti x i ove si valuta la funzione, e il loro numero, si scelgono in modo da rendere il più efficiente possibile il calcolo dell’integrale, cioè cercando di ridurre al minimo il numero delle della funzione f (x ) necessarie per ottenere la precisione richiesta. Per discutere questo argomento conviene introdurre il concetto di grado di precisione di una formula interpolatoria. Tutte le formule interpolatorie hanno grado di precisione almeno pari a n 1, cioè integrano esattamente polinomi con grado minore o uguale a n i . Ad esempio, tramite una formula su 3 punti (n = 3) si integrano esattamente polinomi con grado massimo pari a 2, indipendentemente dalla scelta dei punti x i . Particolari scelte dei punti x i portano a diverse famiglie di formule di quadratura, con diversi gradi di precisione.
−
−
5.4 INTEGRAZIONE NUMERICA
117
Formule di Newton-Cotes Se si scelgono nodi equidistanti che suddividono l’intervallo [ a , b ,] con passo −a costante h = b n −1 : x i = a + (i 1)h , i = 1,...n
−
si costruiscono le formule di N EWON -C OTES . Ad esempio, per n = 2, si ha h = b a , con x 1 = a e x 2 = b , si ottiene la formula dei trapezi : b a I = f (a ) + f (b ) . 2 Nel caso n = 3, invece, h = (b a )/2, x 1 = a , x 2 = (a + b )/2, x 3 = b si ottiene la formula di Cavalieri-Simpson :
−
−
−
I =
b a a + b f (a ) + 4 f 6 2
−
f (c ) . +
Si può dimostrare che, per questa famiglia di formule, il grado di precisione è n 1 se n è pari, mentre è n se n è dispari.
−
Formule di Gauss-Legendre Scegliendo opportunamente le posizioni dei punti ove si valuta la funzione, si ottengono formule con un grado di precisione più elevato. Queste formule appartengono alla famiglia delle formule gaussiane e hanno grado di precisione pari a 2n 1. In particolare, le formule di G AUS S -L EGENDRE scelgono i punti (detti punti di Gauss ) corrispondenti con le radici di certi polinomi, detti appunto polinomi di Legendre , definiti sull’intervallo [ 1,1]. Le posizioni r i in questo intervallo normalizzato prendonoqui il nome di ascisse di Gauss e i corrispondenti pesi αi sono detti pesi di Gauss e sono calcolati, come già illustrato per il caso generale, mediante l’integrazione dei rispettivi polinomi di Lagrange:
−
−
αi =
1
−1
l i (r )dr
La tabella 5.6 riporta alcuni dei polinomi suddetti con ascisse e pesi corrispondenti. Dato il loro grado di precisione più elevato, per l’integrazione sugli elementi finiti si ricorre alle formule di Gauss. Occorre però tenere presente che le ascisse di Gauss sono definite sull’intervallo normalizzato [ 1,1] e quindi si deve tener conto del mapping esistente fra questo e l’intervallo originario [ a , b ]:
−
I =
b
a
f (x )dx =
1
−1
f [x (r )]det J dr
(5.45)
5 FORMULAZIONE DEGLI ELEMENTI FINITI
118
Tabella 5.6: Ascisse e pesi di Gauss per m punti di integrazione. P (m ) (r )
m 1 2
r 1 (3r 2 2
r i
αi
0
2 1
− 1) − 13
1 3 1 3 − 3r ) − 35 2 (5r
3
1 5 9 8 9 5 9
0 3 5
dove det J = dx dr è lo iacobiano della trasformazione. Per cui si avrà: I =
m
1
g (r )dr =
−1
α g (r ) i
i + E m (g )
(5.46)
i =1
Si osservi che l’integrando g (r ), se polinomiale, può avere un ordine più elevato del polinomio f (x ) originale, essendo questo moltiplicato per lo jacobiano det J , in generale funzione di r . Esempio 5.5
Si supponga di dover integrare una funzione parabolica f (x ) = ax 2 + bx + c , su un elemento finito a tre nodi di lunghezza L il cui nodo intermedio sia posto a una distanza βL dall’estremo sinistro ( 0 < β < L ). Si impieghino le formule di quadratura di Newton-Cotes e di Gauss-Legendre, confrontando i risultati con 3 2 il valore esatto I = aL 3 + bL 2 + cL al variare del parametro β. Il mapping a 3 nodi, con x 1 = 0, x 2 = L , x 3 = βL si scrive r (1 r ) r (1 + r ) r 2 x (r ) = ·0+ · L + (1 r ) · βL = L (1 + r ) + β(1 r ) 2 2 2
−
−
−
e lo jacobiano risulta: dx 1 + 2r = L det J (r ) = dr 2
− 2r β
.
La funzione integranda, sull’intervallo normalizzato, g (r ) = f [x (r )]det J ,tenendo conto delle espressioni precedenti risulta essere il polinomio di quinto grado seguente:
1 2
g (r ) = aL
4
2
1 2β 2 β β r + aL β r + r + + bL 2 4 2 bL 1 + 2r r + β2 aL 2 + βbL + c 2βr L βaL 2 + 2 2
− −
2
4
2
1 aL −
3
−
−
5.5 IMPLEMENTAZIONE DEGLI ELEMENTI ISOPARAMETRICI
119
È interessante notare che, nel caso particolare in cui il nodo intermedio sia posto al centro dell’elemento, cioè per β = 1/2, determinante jacobiano è costante eparia L /2, e la funzione integranda si semplifica, diventando un polinomio di secondo grado:
g (r )
β= 12
L aL 2 2 L aL 2 bL r (aL b ) r = + + + + + c 2 4 2 4 2
Il grado di precisione delle formule di quadratura di Newton-Cotes, p N , e d i Gauss-Legendre, p G , a seconda del numero n di punti utilizzati è il seguente: n
1
2
3
4
5
6
p N 1 p G 1
1 3
3 5
3 7
5 9
5 11
Perciò, nel caso generale (polinomio di quinto grado), per ottenere il risultato esatto, risulta necessaria una formula di Newton-Cotes su 5 punti, mentre basta una formula di Gauss-Legendre con soli 3 punti. Si noti, inoltre, che nel caso dell’elemento a jacobiano costante, che si traduce in un integrando polinomiale di secondo grado, la quadratura risulta esatta usando 2 punti di Gauss. Queste considerazioni, svolte sul piano generale, si possono ripetere confrontando i risultati che si ottengono dall’applicazione delle formule di GaussLegendre con numero di punti crescente con quello esatto ottenuto analiticamente: n = 1: I = L βL (βaL + b ) + c . Risultato esatto solo nel caso particolare in cui a = b = 0, cioè f (x ) = c ; oppure per a = 0,b = 0 nel caso β = 1/2, cioè di
elemento non distorto.
L n = 2: I = 18 L 4aL 1 + 2β 2β2 + 9b + 18c . Risultato esatto se a = 0 , cioè f (x ) = b x + c , indipendentemente da β. Inoltre è pure esatto nel caso a = 0 se β = 1/2.
−
n = 3: I = L (2aL + 3b ) + 6c ]. Risultato esatto per tutti i valori dei parametri, 6 [L anche se l’elemento è distorto (infatti non dipende da β).
5.5 Implementazione degli elementi isoparametrici In questo paragrafo sono riportati, negli esempi 5.6 e 5.7, i listati, in linguaggio Fortran, di segmenti di un codice di calcolo 1 relativi a un elemento finito iso1
Tratto da: K.J. Bathe, Finite elements procedures, Prentice-Hall, 1996.
120
5 FORMULAZIONE DEGLI ELEMENTI FINITI
parametrico quadrilatero a quattro nodi, dedicati rispettivamente alla costruzione della matrice di rigidezza K e alla matrice B delle derivate delle funzioni di forma. Nell’interpretazione del codice Fortran, si tenga presente che le linee che iniziano con il carattere , dette commenti , non vengono eseguite dall’elaboratore, ma sono inserite dall’autore del programma a chiarimento e documentazione della procedura. Un carattere qualsiasi nella sesta colonna di una riga indica che essa è la prosecuzione dell’istruzione alla riga precedente. Le istruzioni vere e proprie iniziano dalla colonna 7 e sono eventualmente precedute da un numero, qualora sia necessario fare riferimento a tale istruzione. Esempio 5.6
La subroutine viene chiamata dal programma per ciascun elemento finito quadrilatero appartenente alla mesh, allo scopo di calcolarne il contributo alla matrice di rigidezza della struttura. Le variabili indicate tra parentesi alla linea 1 hanno il significato indicato nelle righe successive (da 2 a 27). La linea 28 afferma che tutte le variabili il cui nome inizia con una lettera compresa, in ordine alfabetico, tra A e H o tra O e Z sono numeri reali in doppia precisione, mentre le rimanenti (cioè quelle che iniziano con lettere tra I e N) sono implicitamente intese come numeri interi. La linea 35 definisce le dimensioni delle matrici e dei vettori. Ad esempio, la matrice , che contiene le coordinate dei nodi dell’elemento, ha 4 righe (una per nodo) e due colonne (rispettivamente la coordinate x e y di ciascun nodo). Nelle linee 39 -42 e 46 -49 si definiscono, nelle matrici e , rispettivamente le ascisse di Gauss e i pesi corrispondenti (v. tabella 5.5). A partire dalla linea 53 si definiscono le costanti elastiche e la matrice costitutiva C, qui indicata con , per i diversi casi (deformazione piana, assialsimmetria, tensione piana). Dalla linea 94 inizia il calcolo dei termini di rigidezza. La matrice di rigidezza K dell’elemento, qui indicata con viene dapprima inizializzata con valori nulli (linee 94-96). La parte tra le linee 99 e 123 è un doppio ciclo che considera il contributo per ogni punto di Gauss (r i , s i ) (sono in direzione r e in direzione s ). Alla linea 106 si chiama la subroutine per la costruzione della matrice K e il calcolo del determinante jacobiano detJ. Alla linea 111 si calcola il termine αi α j det J, costante rispetto alla sommatoria della formula di Gauss per il punto considerato e, dalla linea 112 alla 122, si esegue il doppio prodotto B T CB. I corrispondenti termini vengono assemblati nella matrice di rigidezza dopo averli moltiplicati per αi α j det J alla linea 121. Nelle linee 125-127 si completa, ricordandone la simmetria, la matrice di rigidezza dell’elemento.
5.5 IMPLEMENTAZIONE DEGLI ELEMENTI ISOPARAMETRICI
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
121
122
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
5 FORMULAZIONE DEGLI ELEMENTI FINITI
5.5 IMPLEMENTAZIONE DEGLI ELEMENTI ISOPARAMETRICI
123
125 126 127 128 129 130 131
Esempio 5.7
La subroutine viene utilizzata per il calcolo dei termini della matrice B delle derivate delle funzioni di forma. Essa viene chiamata alla linea 106 della subroutine (riportata nell’esempio 5.6) per ciascun punto di Gauss (r i , s i ). Si possono facilmente individuare le definizioni delle funzioni di forma del quadrilatero a quattro nodi nelle linee 13 -23. Nelle linee 29 -39 si calcolano le derivate rispetto a r e a s . Fra le linee 41 e 63 si valutano in successione la matrice jacobiana, J (qui chiamata ), il suo determinante (chiamato ) e l’inversa ( ). Nelle linee 67 -78 si calcolano gli elementi della matrice B , mentre la parte tra le linee 87 e 108 si riferisce al caso particolare di assialsimmetria per il quale è necessario calcolare la componente di deformazione circonferenziale. Da notare, infine, alla linea 53 , il controllo del valore del determinante jacobiano. Nel caso in cui esso non risulti positivo viene prodotto un messaggio d’errore e il calcolo si arresta dopo un salto alla linea 111. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
124
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
5 FORMULAZIONE DEGLI ELEMENTI FINITI
5.5 IMPLEMENTAZIONE DEGLI ELEMENTI ISOPARAMETRICI
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
125
Appendice A Supporti matematici A.1 Spazi vettoriali Uno spazio vettoriale U è un insieme di elementi, detti vettori , che possono essere moltiplicati per uno scalare e fra i quali è definita un’operazione di somma .
A.1.1 Definizione Formalmente, U è uno spazio vettoriale su guenti:
R se
sono definite le due leggi se-
Legge di composizione interna: U ×U → U che a ogni coppia u , v di elementi di U associa un elemento di U , indicato con u + v , detto somma di u e v . Legge di composizione esterna: R × U → U che ad ogni coppia costituita da uno scalare α ∈ R e da un elemento di U associa un elemento di U , indicato con αu , detto prodotto di α per u . Per cui sussistano le seguenti i proprietà : 1. u + v = v + u (proprietà commutativa della somma); 2. (u + v ) + w = u + (v + w ) (proprietà associativa della somma); 3. esiste un elemento di 0U ∈ U , detto vettore nullo , tale che u + 0U = u per ogni u ∈ U (esistenza dell’elemento neutro rispetto alla somma); 4. per ogni u ∈ U esiste un elemento, indicato con −u tale che u + (−u ) = 0U (esistenza dell’opposto); 5. α(u + v ) = αu + αv (proprietà distributiva del prodotto rispetto alla somma in U ); 115
A SUPPORTI MATEMATICI
116
6. (α + β)u = αu + βu (proprietà distributiva del prodotto rispetto alla somma in R); 7. α(βu ) = (αβ)u (proprietà associativa del prodotto); 8. 1 · u = u (esistenza dell’elemento neutro rispetto al prodotto). valide qualunque siano gli scalari α, β ∈ R e gli elementi u , v , w ∈ U . Dalla definizione assiomatica precedente si fanno discendere tutte le proprietà degli spazi vettoriali. Il concetto di spazio vettoriale (detto anche spazio lineare ) è uno dei più importanti della matematica. Esso permette di studiare insiemi di oggetti matematici, a prima vista fra loro diversissimi, ma dotati di una somiglianza strutturale profonda. Negli esempi seguenti si presentano alcuni spazi vettoriali utilizzati nel seguito. Esempio A.1
L’insieme Rn = R × R × · · · × R ( n volte) delle n-uple ordinate di scalari u ∈ Rn : u = (u 1 , u 2 , . . . ,u n ), con u 1 , u 2 , . . . ,u n ∈ R
(A.1)
con somma: u + v = (u 1 + v 1 , . . . ,u n + v n ) e prodotto αu = (αu 1 , . . . ,αu n ).
Esempio A.2
Le funzioni definite in un dominio Ω che posseggono derivate continue fino all’ordine m , con le operazioni di somma e prodotto definite in modo naturale, costituiscono uno spazio vettoriale, detto C m (Ω). Si verifica facilmente, infatti che la somma (αu + βv )( x ) = αu ( x ) + βv ( x ) appartiene ancora a C m (Ω), l’elemento neutro è la funzione identicamente nulla su Ω e l’opposta è −u ( x ) = −1 · u ( x ).
A.1.2 Sottospazi Non tutti i sottoinsiemi di elementi tratti da uno spazio vettoriale continuano a soddisfare gli assiomi di definizione. Quei sottoinsiemi i cui elementi costituiscono uno spazio vettoriale sono chiamati sottospazi . Formalmente, un sottoinsieme V di uno spazio vettoriale U su R , è un sottospazio di U se V è uno spazio vettoriale su R rispetto alle stesse leggi di composizione definite in U . Si dimostra facilmente che V ⊂ U è un sottospazio di U se e solo se u + v ∈ V , ∀u , v ∈ V αu ∈ V , ∀u ∈ V , ∀u ∈ R
(A.2)
(A.3)
A.1 SPAZI VETTORIALI
117
cioè è chiuso rispetto alle leggi di composizione definite in U . Esempio A.3
Si consideri lo spazio V 3 dei vettori liberi dello spazio tridimensionale. I vettori paralleli rispettivamente a una retta o a un piano sono sottospazi vettoriali. È immediato infatti verificare la sussistenza delle proprietà espresse dalle eq. (A.2) e (A.3), cioè, ad es. la somma di due vettori complanari è un vettore ancora appartenente allo stesso piano, così come continua ad appartenervi il vettore moltiplicato per uno scalare. Dati due sottospazi V , W ⊂ U , il sottoinsieme Z di U i cui elementi sono dati da z = v + w , al variare di v in V e di w in W è un sottospazio di U che si dice sottospazio somma di V e W e si indica: Z = V + W Se, in particolare, l’intersezione degli spazi coincide con il vettore nullo, cioè V ∩ W = {0U }, si dice che Z è somma diretta di V e W , e si scrive: Z = V ⊕ W
Si dimostra un utile teorema che afferma che la somma Z di due sottospazi V e W di U è diretta se e solo se ogni vettore z ∈ Z si scrive in un solo modo nella forma: z = v + w , v ∈ V , w ∈ W Nel caso in cui U = V ⊕ W , questi ultimi si dicono sottospazi supplementari . Esempio A.4
Si consideri nuovamente lo spazio V 3 dei vettori liberi dello spazio tridimensionale. I vettori paralleli rispettivamente a una retta e a un piano ad essa non parallelo sono sottospazi supplementari di V 3 .
A.1.3 Prodotto interno Sia U uno spazio vettoriale: si definisce prodotto interno 〈u , v 〉 di due vettori u , v ∈ U un’operazione U × U → R che soddisfa gli assiomi seguenti: 1. 〈u , v 〉 = 〈v , u 〉 (simmetria);
A SUPPORTI MATEMATICI
118
2. 〈αu + βv , w 〉 = α〈u , w 〉 + β〈v , w 〉 (linearità); 3. 〈u , u 〉 ≥ 0, con 〈u , u 〉 = 0 se e solo se u = 0 (definita positività). Uno spazio vettoriale dotato di prodotto interno è detto spazio euclideo . Esempio A.5
Date due funzioni u (x ) e v (x ) continue sull’intervallo Ω = (a , b ), considerate come elementi dello spazio C 0 (Ω), un valido prodotto interno è dato da: b
〈u , v 〉 =
u (x )v (x ) d x
(A.4)
a
L’assioma 1 è soddisfatto in quanto 〈v , u 〉 = a b v (x )u (x ) d x = 〈u , v 〉. La linearità (assioma 2) si verifica in quanto, essendo l’integrando non negativo in tutto l’intervallo:
b
〈αu + βv , w 〉 =
(αu (x ) + βv (x ))w (x ) d x
a
b
=α
b
u (x )w (x ) d x + β
a
v (x )w (x ) d x
(A.5)
a
= α〈u , w 〉 + β〈v , w 〉
Il terzo assioma è verificato in quanto b
〈u , u 〉 =
u 2 (x ) d x ≥ 0
(A.6)
a
e solamente nel caso u (x ) = 0 assume il valore nullo. Il prodotto interno qui introdotto nel caso di uno spazio vettoriale qualsiasi è una generalizzazione del concetto di prodotto scalare in V 3 , ed è quindi immediato generalizzare ad uno spazio vettoriale qualsiasi altre idee legate a tale prodotto. In tal caso, come è noto, due vettori u e v sono ortogonali se u · v = 0 (con · si è qui indicato l’usuale prodotto scalare in V 3 ). In maniera analoga è possibile introdurre il concetto di perpendicolarità in uno spazio vettoriale qualsivoglia, purché dotato di prodotto interno. Due vettori u , v ∈ U , sono fra loro ortogonali se 〈u , v 〉 = 0
Esempio A.6
(A.7)
A.1 SPAZI VETTORIALI
119
Siano u (x ) = αx e v (x ) = βx 2. Considerate come vettori dello spazio C ∞(Ω), Ω = (−a , a ), sono ortogonali, infatti: 〈u , v 〉 =
a
a
u (x )v (x ) d x =
−a
(αx )(βx 2 ) d x
−a
x 4 3 x d x = αβ = αβ 4 −a
a
(A.8)
a
=0
−a
A.1.4 Norme e metrica Detto u un qualsiasi vettore di uno spazio vettoriale U , si definisce norma di u , u , un’operazione che fornisce un valore reale in base ai seguenti assiomi: 1. u ≥ 0, con u = 0 se e solo se u = 0 (definita positività); 2. αu = |α| u (linearità); 3. u + v ≤ u + v (disuguaglianza triangolare). per ogni u , v ∈ U , α ∈ R. Gli assiomi qui introdotti si ispirano, come è facile verificare, al concetto di modulo di un vettore in V 3 e ne estendono la validità a spazi vettoriali qualsiasi. Esempio A.7
Si consideri lo spazio vettoriale Rn . L’espressione n
u=
1/α
α
|u k |
(A.9)
k =1
con 1 ≤ α < ∞, è una norma valida. Si noti che, scegliendo α = 2 si ottiene la norma euclidea del vettore u, la quale fornisce una generalizzazione del teorema di Pitagora in Rn . Uno spazio vettoriale dotato di norma è detto spazio normato . In generale, qualsiasi funzione che rispetti gli assiomi di definizione può essere usata come norma, ma nel caso particolare in cui lo spazio vettoriale sia dotato di prodotto interno è spesso conveniente definire una norma basandosi su tale prodotto. Negli spazi dotati di prodotto interno si può infatti definire una norma naturale, a partire dal prodotto interno nel modo seguente: 1/2
u = 〈u , u 〉
(A.10)