Metodi FEM per problemi ellittici

Gabriella Puppo

Principali argomenti

• Elementi finiti P2

• Velocità di convergenza

• Problema di convezione/diffusione

Filo elastico con elementi P2

• 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 le3 funzioni della base.E’ facile vedere che si ottiene:

l 1 (x) = (2x - 1) (x - 1)l 2 (x) = - 4x (x - 1)l 3 (x) = x (2x - 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:

>> l1=(2*x-1)*(x-1)l1 =(2*x-1)*(x-1)

Per calcolare l’elemento di posto (1,1) della matrice di rigidità, devo:

Calcolare la derivata del polinomio:

Moltiplicare la derivata per se stessa e integrare fra 0 e 1:

>> dl1=diff(l1)dl1 =4*x-3

>> a(1,1)=int(dl1*dl1,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'*dlans =[ (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 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/(3h)

Ottengo così la matrice di rigidità globale.

Per N intervalli ho 2N -1 nodi interni. Nel caso N = 5:

a = 16 -8 0 0 0 0 0 0 0 -8 14 -8 1 0 0 0 0 0 0 -8 16 -8 0 0 0 0 0 0 1 -8 14 -8 1 0 0 0 0 0 0 -8 16 -8 0 0 0 0 0 0 1 -8 14 -8 1 0 0 0 0 0 0 -8 16 -8 0 0 0 0 0 0 1 -8 14 -8 0 0 0 0 0 0 0 -8 16

Devo infine moltiplicare ogni elemento per 1/(3h)

Esercizio

Scrivere una function che calcoli la matrice di rigidità globale, assegnando in input il numero N di intervalli, supponendo che il fattore 1/(3h) 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 P2, e condizioni di Dirichlet.% N e' il numero di intervalli % => A e' (2N-1)*(2N-1)% Il fattore 1/(3h) e’ stato portato a destra

continua ...

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

… continua

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 vaspezzato 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 P2.Inoltre, (x) 0 su un solo nodo della formula di quadratura.

Il vettore di carico quindi è:

Questa grandezza deve poi essere moltiplicata per 3h.

Programma FEM P2

Posso ( finalmente!) scrivere il programma che risolve il problema del filo elastico con condizioni di Dirichlet.Questo programma chiama la function penta che costruisce la matrice di rigidità per gli elementi P2.

function [x,u]=filo_p2(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 grigliam=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/2fx = feval(f,x);for i=1:2:m b(i) = 2*h/3 *fx(i); % nodi dispariendfor i = 2:2:m-1 b(i) = h/3 *fx(i); % nodi pariendb=3*h*b'; % ho portato 3h a destrau=a\b;% Ora aggiungo le condizioni al contornou(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_p2(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 P1 che gli elementi P2

Considero un problema di cui conosco la soluzione esatta, per esempio:

Valuto l’ errore su diverse griglie

exa=inline('x.*(x-1).*exp(x)');f=inline('-exp(x).*(x.^2+3*x)');n=10;fprintf(' n P1 P2 \n')for k=1:4 [x,u]=filo_dir(f,n); err_p1(k)=norm(exa(x)-u',inf); [x,u]=filo_p2(f,n); err_p2(k)=norm(exa(x)-u',inf); fprintf('%3.0f %12.5e %12.5e \n',n,... err_p1(k),err_p2(k)) n=n*2;end

Calcolo l’ errore usando griglie via via più fitte. Per semplicità uso la norma del massimo:

Risultati

N P1 P210 1.784725544444921e-03 2.082890535519071e-0520 4.910462530513526e-04 1.844384948512536e-0640 1.288307737212224e-04 1.630962657862395e-07

Ottengo questi risultati:

Qui N è il numero di intervalli. Il metodo agli elementi P2usa circa il doppio delle incognite del metodo P1, 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 P1 e con elementi P2.

Chiamo:CP1(N) = numero di condizionamento matrice P1 con N intervalliCP2(N) = numero di condizionamento matrice P2 con N intervalliDisegno il grafico di CP1 e CP2 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 P2, nel caso di un carico puntiforme.

Esercizio

Considerare una funzione di carico molto oscillante, tipof(x)= -1 + 10*sin(10*pi*x).Studiare la deformazione del filo, usando elementi P1 e elementi P2.Come posso sfruttare il fatto che la soluzione prodotta dagli elementi P2 è 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_cd1d(n,nu,beta)% A=MAT_CD1D(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 internia=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;enda=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_cd1d(n,nu,beta);h=1/(n+1); % spaziatura di griglia% calcola la posizione dei nodi interni:for i=1:n; x(i) = i*h; endb=h^2*feval(f,x); b=b'; % ho portato h a sinistrau=a\b;% Ora aggiungo le condizioni al contornou(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

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

Il numero di Peclet deve essere minore di 1 dappertutto o solo nello strato limite?

Esercizio

Modificare il problema di convezione diffusione aggiungendo una condizione di Neumann. Provare i due casi seguenti:

2) Condizione di Neumann in x=0, Dirichlet in x=1 con β >0;3) Condizione di Neumann in x=0, Dirichlet in x=1 con β <0.

Ottengo comunque uno strato limite per valori di β/ν elevati?