Metodi FEM per problemi ellittici Gabriella Puppo Principali

Metodi FEM per problemi ellittici Gabriella Puppo

Principali argomenti • Elementi finiti P 2 • Velocità di convergenza • Problema di convezione/diffusione

Filo elastico con elementi P 2 • • Costruisco gli elementi della base Calcolo la matrice di rigidità di elemento Calcolo la matrice di rigidità globale Costruisco una function che calcola il carico e risolve il sistema • Disegno la soluzione, come polinomio quadratico a tratti

Costruzione delle funzioni di base di elemento Considero l’intervallo [0, 1] e costruisco su questo intervallo le 3 funzioni della base. E’ facile vedere che si ottiene: l 1 (x) = (2 x - 1) (x - 1) l 2 (x) = - 4 x (x - 1) l 3 (x) = x (2 x - 1)

Costruendo i grafici delle tre funzioni, si ottiene:

Costruzione della matrice di rigidità di elemento La costruzione della matrice di base di elemento è un’operazione noiosissima, nella quale inoltre è facile commettere errori di calcolo. La situazione migliora tantissimo se si ha a disposizione un pacchetto per il calcolo simbolico. Poiché la Student edition di Matlab ha un pacchetto di calcolo simbolico, nel seguito calcoliamo la matrice di rigidità, facendo fare i conti a Matlab.

Fare i conti con il calcolo simbolico L’utilizzo del pacchetto di calcolo simbolico è molto semplice. Innanzitutto, devo dire che x è una variabile simbolica: >> syms x A questo punto, ogni variabile alla quale venga assegnata una espressione contenente x diventa a sua volta simbolica: >> l 1=(2*x-1)*(x-1) l 1 = (2*x-1)*(x-1)

Per calcolare l’elemento di posto (1, 1) della matrice di rigidità, devo: Calcolare la derivata del polinomio: >> dl 1=diff(l 1) dl 1 = 4*x-3 Moltiplicare la derivata per se stessa e integrare fra 0 e 1: >> a(1, 1)=int(dl 1*dl 1, 0, 1) a= 7/3

In realtà è più efficiente procedere in modo vettoriale: Calcolo un vettore con i tre polinomi: >> l=[(2*x-1)*(x-1), -4*x*(x-1), x*(2*x-1)] Ne calcolo la derivata: >> dl=diff(l) dl = [ 4*x-3, -8*x+4, 4*x-1]

Devo ora moltiplicare le derivate in tutti i modi possibili, cioè devo calcolare dl’*dl. Matlab però considera x una variabile complessa: devo specificare che x è reale: >> syms x real A questo punto posso calcolare: >> dl'*dl ans = [ (4*x-3)^2, (4*x-3)*(-8*x+4), (4*x-3)*(4*x-1)] [ (4*x-3)*(-8*x+4), (-8*x+4)^2, (-8*x+4)*(4*x-1)] [ (4*x-3)*(4*x-1), (-8*x+4)*(4*x-1), (4*x-1)^2]

Infine, calcolo l’integrale fra 0 e 1 di quello che ho ottenuto: >> a=int(ans, 0, 1) a= [ 7/3, -8/3, 1/3] [ -8/3, 16/3, -8/3] [ 1/3, -8/3, 7/3] Questa è la matrice di rigidità per un elemento di lunghezza 1. Per riscalare su un elemento di lunghezza h devo dividere per h. Notare che le frazioni sono rimaste indicate: questo risultato è esatto

Costruisco ora la matrice di rigidità globale La prima riga della matrice di rigidità è relativa alla funzione di base a “bolla”, centrata sul primo nodo interno: [16 -8 0 0 0. . . La seconda riga della matrice di rigidità è relativa alla funzione di base a “tenda”, centrata sul secondo nodo interno: [ -8 14 -8 1 0 0 0. . . Tutti gli elementi devono poi essere moltiplicati per 1/(3 h)

Ottengo così la matrice di rigidità globale. Per N intervalli ho 2 N -1 nodi interni. Nel caso N = 5: a= 16 -8 0 0 0 0 -8 14 -8 1 0 0 0 -8 16 -8 0 0 0 1 -8 14 -8 1 0 0 0 -8 16 -8 0 0 0 1 -8 14 -8 0 0 0 0 -8 16 Devo infine moltiplicare ogni elemento per 1/(3 h)

Esercizio Scrivere una function che calcoli la matrice di rigidità globale, assegnando in input il numero N di intervalli, supponendo che il fattore 1/(3 h) sia portato a destra del segno =.

Ottengo questo listato: function a=penta(n) % A=PENTA(N) costruisce la matrice di rigidita' per % il problema del filo elastico, trattato con % elementi P 2, e condizioni di Dirichlet. % N e' il numero di intervalli % => A e' (2 N-1)*(2 N-1) % Il fattore 1/(3 h) e’ stato portato a destra continua. . .

… continua m=2*n-1 % dimensioni della matrice a=16*eye(m); for i=2: 2: m-1; a(i, i)=14; end b=ones(m-1, 1); a=a+diag(-8*b, -1)+diag(-8*b, 1); % inserisce i +1 sulla diagonale +2 for i=2: 2: m-3 a(i, i+2)=1; end % inserisce i +1 sulla diagonale -2 for i=4: 2: m-1 a(i, i-2)=1; end

Calcolo del vettore di carico Il vettore b è dato da L’indice i denota tutti i nodi. Dobbiamo distinguere i nodi pari dai nodi dispari, perché le funzioni di base hanno forma e supporto diversi.

Nodi pari e nodi dispari Sui nodi dispari, la funzione della base è una funzione “a bolla”. Il suo supporto è l’intervallo [x i - h/2, x i + h/2]. Sui nodi pari, la funzione della base è una funzione “a tenda”. Il suo supporto è l’intervallo [x i - h, x i + h]. Inoltre, poiché queste funzioni hanno un punto angoloso in x i , l’integrale va spezzato in due contributi, da [x i - h, x i ] a [x i , x i + h].

Calcolo gli integrali usando la formula di Simpson Questa formula ha un errore O(b-a) 5, ed ha quindi precisione sufficiente per gli elementi P 2. Inoltre, (x) 0 su un solo nodo della formula di quadratura.

Il vettore di carico quindi è: Questa grandezza deve poi essere moltiplicata per 3 h.

Programma FEM P 2 Posso ( finalmente!) scrivere il programma che risolve il problema del filo elastico condizioni di Dirichlet. Questo programma chiama la function penta che costruisce la matrice di rigidità per gli elementi P 2. function [x, u]=filo_p 2(f, n) % Risolve il problema del filo elastico su (0, 1) % per condizioni al bordo omogenee di % Dirichlet in x=0 e in x=1, % con carico definito dalla funzione % F e N intervalli

a=penta(n); h=1/n; % spaziatura di griglia m=2*n-1 % numero dei nodi interni % calcola la posizione dei nodi interni: for i=1: m; x(i) = i*h/2; end % la spaziatura e' h/2 fx = feval(f, x); for i=1: 2: m b(i) = 2*h/3 *fx(i); % nodi dispari end for i = 2: 2: m-1 b(i) = h/3 *fx(i); % nodi pari end b=3*h*b'; % ho portato 3 h a destra u=ab; % Ora aggiungo le condizioni al contorno u(2: m+1)=u; u(1)=0; u(m+2)=0; x(2: m+1)=x; x(1)=0; x(m+2)=1;

Esempio Risolvo il problema del filo elastico con carico uniforme >>f=inline('-1+0*x'); >> [x, u]=filo_p 2(f, 10); >> plot(x, u, 'Linewidth', 2)

Velocità di convergenza Valutiamo la velocità di convergenza per il problema del filo elastico, utilizzando sia gli elementi P 1 che gli elementi P 2 Considero un problema di cui conosco la soluzione esatta, per esempio:

Valuto l’ errore su diverse griglie Calcolo l’ errore usando griglie via più fitte. Per semplicità uso la norma del massimo: exa=inline('x. *(x-1). *exp(x)'); f=inline('-exp(x). *(x. ^2+3*x)'); n=10; fprintf(' n P 1 P 2 n') for k=1: 4 [x, u]=filo_dir(f, n); err_p 1(k)=norm(exa(x)-u', inf); [x, u]=filo_p 2(f, n); err_p 2(k)=norm(exa(x)-u', inf); fprintf('%3. 0 f %12. 5 e n', n, . . . err_p 1(k), err_p 2(k)) n=n*2; end

Risultati Ottengo questi risultati: N 10 20 40 P 1 1. 784725544444921 e-03 4. 910462530513526 e-04 1. 288307737212224 e-04 P 2 2. 082890535519071 e-05 1. 844384948512536 e-06 1. 630962657862395 e-07 Qui N è il numero di intervalli. Il metodo agli elementi P 2 usa circa il doppio delle incognite del metodo P 1, ma l’ errore è molto più piccolo

Esercizio Costruire un programma che ricostruisca l’andamento del numero di condizionamento della matrice di rigidita’ del filo elastico al variare di N, con elementi P 1 e con elementi P 2. Chiamo: CP 1(N) = numero di condizionamento matrice P 1 con N intervalli CP 2(N) = numero di condizionamento matrice P 2 con N intervalli Disegno il grafico di CP 1 e CP 2 in funzione di N, e provo ad interpolare i dati con un polinomio di grado assegnato nel senso dei minimi quadrati (cioe’ con polyfit)

Esercizio Studiare il comportamento del problema del filo elastico con elementi P 2, nel caso di un carico puntiforme.

Esercizio Considerare una funzione di carico molto oscillante, tipo f(x)= -1 + 10*sin(10*pi*x). Studiare la deformazione del filo, usando elementi P 1 e elementi P 2. Come posso sfruttare il fatto che la soluzione prodotta dagli elementi P 2 è quadratica a tratti per migliorare il grafico?

Problema di convezione-diffusione Vogliamo risolvere questo problema al contorno: Che è equivalente alla formulazione variazionale:

Matrice di rigidità La matrice di rigidità di questo problema è una matrice tridiagonale, con: N. B. Il fattore h che moltiplica la matrice di rigidità è stato portato a destra dell’ uguale e quindi moltiplica il vettore di carico.

Questo è il listato di una function che costruisce la matrice di convezione - diffusione. function a=mat_cd 1 d(n, nu, beta) % A=MAT_CD 1 D(N, NU, BETA) Crea la matrice N per N tridiagonale del % filo elastico per il problema di convezione-diffusione, % con campo di velocita' BETA sull'intervallo (0, 1) % e diffusione NU % N e' il numero di nodi interni a=2*eye(n); h=1/(n+1); peclet = h*beta/2/nu; for i=1: n-1 a(i, i+1)=-1 +peclet; a(i+1, i)=-1 -peclet; end a=a*nu;

Il listato del programma principale è molto simile a quelli già visti: function [x, u]=c_d(f, n, nu, beta) % [X, U]=C_D(F, N, NU, BETA) % Risolve il problema del filo elastico su (0, 1) % per condizioni al bordo omogenee di % Dirichlet, con carico definito dalla funzione % F e N punti interni % NU e' il coefficiente di diffusione % BETA e' il campo di velocita' a=mat_cd 1 d(n, nu, beta); h=1/(n+1); % spaziatura di griglia % calcola la posizione dei nodi interni: for i=1: n; x(i) = i*h; end b=h^2*feval(f, x); b=b'; % ho portato h a sinistra u=ab; % Ora aggiungo le condizioni al contorno u(2: n+1)=u; u(1)=0; u(n+2)=0; x(2: n+1)=x; x(1)=0; x(n+2)=1;

Carico uniforme, =1, N=100

Carico puntiforme, =1, N=100

Negli esempi precedenti, il numero di Peclet era sempre minore di 1: Pe = h / 2 Se aumentiamo il numero di Peclet, la soluzione diventa oscillante

Convezione-diffusione, =100, N=100

Esercizio Studiare il comportamento delle oscillazioni spurie per il problema di convezione diffusione rispetto a variazioni nel numero di Peclet. Come si comporta l’ampiezza delle oscillazioni rispetto a Pe?

Esercizio Risolvere il problema: e verificare che anche in questo caso il comportamento oscillatorio della soluzione e’ legato al numero di Peclet.

Esercizio Il numero di Peclet deve essere minore di 1 dappertutto o solo nello strato limite? 1) Scrivere una function che calcoli la matrice di convezione diffusione per un vettore x, che contenga le ascisse di nodi di griglia assegnati (cioe’ non uniformi). 2) Modificare di conseguenza anche il vettore di carico. 3) Costruire una griglia ad-hoc che sia piu’ fitta all’interno dello strato limite e risolvere il problema di convezione-diffusione su quella griglia

Esercizio Modificare il problema di convezione diffusione aggiungendo una condizione di Neumann. Provare i due casi seguenti: 1) Condizione di Neumann in x=0, Dirichlet in x=1 con β >0; 2) Condizione di Neumann in x=0, Dirichlet in x=1 con β <0. Ottengo comunque uno strato limite per valori di β/ν elevati?