Prima di procedere alla illustrazione rigorosa degli step che formano il FEM è utile avere bene in mente il principio dei lavori virtuali. SPOILER ALERT! Esso permette di riscrivere il problema passando DA una sua rappresentazione in “forma forte” – ovvero la sua forma differenziale – accompagnata, ovviamente, da tutti i vincoli necessari affinché sia possibile arrivare ad una unica soluzione, AD una sua rappresentazione in “forma debole”– una formulazione puramente integrale, in cui i vincoli vengono inglobati all’interno dell’integrale stesso.
Procediamo lentamente così da non lasciarci sfuggire nulla durante il percorso che ci condurrà ad implementare questo metodo nel codice. Il principio dei lavori virtuali, come suggerisce il nome, è un principio energetico e rientra nella categoria dei principi variazionali. Visto l’ampio uso che se ne fa in Fisica, ne esistono varie versioni, l’originaria delle quali, si narra, sia dovuta a Lagrange. Il principio tiene conto del fatto che esiste sempre e solo una unica configurazione vera per il sistema e che quest’ultimo evolva per minimizzare la energia associata ad ogni possibile stato di configurazione. Nel momento in cui il sistema si troverà ad essere in una sua configurazione di equilibrio, tutte le altre possibili configurazioni sono non reali o, appunto, virtuali.
Per lavoro virtuale si intende quel lavoro svolto da forze responsabili di spostamenti virtuali; uno spostamento è virtuale se rappresenta una deformazione infinitesima, reversibile e compatibile coi vincoli cui è soggetto il sistema.
Vediamo di essere un po’ più precisi introducendo una definizione più rigorosa riprendendo il sistema esempio della trave appesa al soffitto. Nella sua formulazione dovuta a Lagrange il principio enuncia:
Condizione Necessaria e Sufficiente perché un sistema materiale sia in equilibrio, è che sia nullo il lavoro virtuale delle forze attive associato a qualunque spostamento virtuale compatibile
Intuitivamente possiamo pensarla così; se il sistema è all’equilibrio, e quindi nella sua deformazione vera, e lo si perturba infinitesimamente allora le forze esterne applicate non svolgono lavoro e il sistema non cambia stato energetico. L’energia totale è in condizione di stazionarietà.
Prendiamo in considerazione il sistema meccanico composta da una barra elastica con entrambi gli estremi ben fissi. L’equazione del moto sarà la stessa che ci siamo ricavati la scorsa volta (vedi qui), a cambiare saranno solamente le condizioni al contorno:
$$\rho (x)\frac{\partial^2 u}{\partial t^2}(x,t) – \frac{\partial }{\partial x}(k(x)\frac{\partial u}{\partial x}(x,t))=f(x,t)$$
$$u(0,t)=0$$
$$\frac{\partial u}{\partial x}(l,t)=0$$
L’obiettivo da perseguire è il ri-esprimere l’equazione del sistema barra vincolata come il minimo di un certo potenziale, grazie all’aiuto del PLV. Un trucco prevede di moltiplicarla preliminarmente per $\frac{\partial u }{\partial t}$ e, dopodiché, integrare
$$\int_{0}^{l}\left [ \rho (x)\frac{\partial^2 u}{\partial t^2}(x,t)\frac{\partial u}{\partial t}(x,t) – \frac{\partial }{\partial x}(k(x)\frac{\partial u}{\partial x}(x,t))\frac{\partial u}{\partial t}\right ]dx=0$$
Utilizzo la integrazione per parti per il secondo addendo
$$\int_{0}^{l}\frac{\partial }{\partial x}(k(x)\frac{\partial u}{\partial x}(x,t))\frac{\partial u}{\partial t}dx=$$
$$ -k(x)\frac{\partial u}{\partial x}(x,t)\frac{\partial u}{\partial t}(x,t) + \int_{0}^{l} k(x)\frac{\partial u}{\partial x}(x,t)\frac{\partial^2 u}{\partial t \partial x}(x,t)$$
Il primo termine si annulla per i vincoli imposti al sistema, infatti, sia in x=0 che in x=l abbiamo posto nulla la funzione o la sua derivata. Successivamente, rielaboro sia il secondo termine che il primo addendo dell’integrale precedente come
$$\int_{0}^{l} k(x)\frac{\partial u}{\partial x}\frac{\partial^2 u}{\partial t \partial x}dx= \int_{0}^{l} k(x)\frac{\partial }{\partial t}\left [\frac{1}{2}(\frac{\partial u}{\partial x})^2\right ]dx$$
$$\int_{0}^{l}\rho (x)\frac{\partial^2 u}{\partial t^2}\frac{\partial u}{\partial t}dx = \int_{0}^{l}\rho (x)\frac{\partial }{\partial t}\left [ \frac{1}{2}(\frac{\partial u}{\partial t})^2 \right ]dx$$
Così da ottenere, complessivamente
$$\int_{0}^{l}\rho (x)\frac{\partial }{\partial t}\left [ \frac{1}{2}(\frac{\partial u}{\partial t})^2 \right ]dx + \int_{0}^{l} k(x)\frac{\partial }{\partial t}\left [ \frac{1}{2}(\frac{\partial u}{\partial x})^2\right ]dx=0$$
è possibile portare fuori dall’integrale la derivata temporale
$$\frac{\mathrm{d} }{\mathrm{d} t}\left [ \frac{1}{2}\int_{0}^{l}\rho (x)(\frac{\partial u}{\partial t})^2 dx + \frac{1}{2}\int_{0}^{l} k(x)(\frac{\partial u}{\partial x})^2dx\right ]=0$$
Questo deve esser vero $x \in (0,l)$ e ciò significa che è la quantità tra parentesi a doversi annullare. La somma dei due integrali è costante nel tempo. Qui scatta il campanello d’allarme. Conosco benissimo una quantità che rimane costante al variare di t: è la sua coniugata, l’energia. In particolare il termine proporzionale a $\rho(x)$ esprime l’energia cinetica del sistema meccanico (1/2*massa*velocità^2) mentre il termine proporzionale a $k(x)$ la sua energia potenziale. Questa rappresenta l’energia elastica interna dovuta alla deformazione della barra. Ora occorre prestare attenzione. Consideriamo la situazione in cui il sistema si trovi ad essere all’equilibrio: la derivata temporale non darà contributo, posso qundi lasciar perdere il contenuto cinetico e tenere solo quello potenziale.
$$\frac{1}{2}\int_{0}^{l} k(x)(\frac{\partial u}{\partial x})^2dx=0$$
L’energia potenziale mappa lo stato delle configurazioni ammissibili per il sistema. La configurazione all’equilibrio è la configurazione che minimizza l’apporto energetico, per cui se il sistema si trova ad essere all’equilibrio, allora qualsiasi altra deformazione dovrà restituire una energia potenziale maggiore. Quindi, qualsiasi altra energia potenziale che non si riferisce alla posizione di equilibrio del sistema, può essere vista come somma di due parti: quella all’equilibrio e una qualsiasi altra che dà luogo a un contributo positivo.
L’espressione completa per la energia potenziale totale della barra all’equilibrio è la seguente:
$$\varepsilon_{pot}(u)=\frac{1}{2}\int_{0}^{l} k(x)(\frac{\partial u}{\partial x})^2dx-\int_{0}^{l}f(x)u(x)dx + \varepsilon _{0}$$
- il primo termine $\frac{1}{2}\int_{0}^{l} k(x)(\frac{\partial u}{\partial x})^2dx$ rappresenta l’energia potenziale di tipo elastico che caratterizza il sistema. Essa dipende dalla dalla posizione relativa delle parti che costituiscono il sistema e viene chiamata energia potenziale interna
- il secondo termine, invece, $\int_{0}^{l}f(x)u(x)dx$ rappresenta l’energia potenziale di tipo gravitazionale o qualsiasi altro tipo che dipende solamente dalla posizione assoluta della barra. Per questo viene indicata come energia potenziale esterna
- l’ultimo contributo $\varepsilon _{0}$ è quella parte di energia potenziale esterna che corrisponde alla barra nella sua configurazione all’equilibrio
E’ possibile dimostrare quanto detto semplicemente andando a sostituire la funzione u(x,t) nella espressione integrale di cui sopra, la funzione w(x,t) intendendo con ciò, una qualsiasi funzione diversa da u(x,t) compatibile coi vincoli, che rappresenti la deformazione del sistema. Voglio poter dimostrare che
$$\epsilon_{pot}(w)>\epsilon_{pot}(u)$$
posso esprimere la nuova funzione come somma w(x,t)=u(x,t)+v(x,t) con la condizione che u(0)=w(0)=v(0)=0 e u(l)=w(l)=v(l)=0. Quindi
$$\epsilon_{pot}(u+v)>\epsilon_{pot}(u)$$
E’ anche possibile dimostrare che le due forme debole e forte sono del tutto equivalenti. Entrambe queste dimostrazioni sono ben sviluppate nel capitolo V del libro di Mark S. Gockenbach – Partial Differential Equations – 2a edizione che in generale sto seguendo.
Il risultato davvero importante quando si applica il PLV a questo genere di problemi e che ritengo giusto evidenziare ancora, sta nella nuova visione che ci da nella loro risoluzione. Uno stato all’equilibrio corrisponde allo stato in cui l’energia potenziale del sistema viene minimizzata.
La forma debole ottenuta applicando il PLV al problema differenziale coinvolge un numero infinito di equazioni visto che esistono infinite funzioni di test v(x,t) $\in$ V. Allo stesso modo, anche il problema di partenza espresso in forma forte comporta infinite equazioni poiché devono mantenersi valide per ogni x $\in$ [0,l]. In entrambi i casi è impossibile è impossibile arrivare alla soluzione direttamente utilizzando una procedura numerica che vanta comunque un numero finito di passi. Quale vantaggio ci porta il riformulare il problema in forma integrale, allora?
Semplice! La formulazione debole ammette l’uso del metodo di Galërkin, una macro-categoria che include il famoso metodo degli elementi finiti. Il medico di Galërkin permette di ridurre il numero infinito delle equazioni ad una collezione finita . Fantastico, no?
Galërkin lo vedremo nel prossimo articolo e il punto di partenza sarà la formulazione debole del problema al contorno che abbiamo ottenuto da una applicazione diretta del principio dei lavori virtuali. Ancora più chiaramente, posso riesprimere questa forma impiegando il prodotto interno su $L^2$ $$(f,v)=\int_{0}^{l}f(x)v(x)dx$$ e introducendo la seguente forma bilineare e simmetrica $$a(u,v)=\int_{0}^{l}k(x)\frac{\mathrm{d} u}{\mathrm{d} x}(x)\frac{\mathrm{d} v}{\mathrm{d} x}(x)dx$$Così il problema finale diventa trovare quella funzione u $\in$ V tale che
$$a(u,v)=(f,v)$$
per ogni v $\in$ V
Per proseguire con il prossimo articolo clicca qui