Metodi FEM in 2DMetodi FEM in 2D

G. Puppo

RiassuntoRiassunto

Autovalori e aliasingMembrana elasticaProblema di convezione-diffusione agli

elementi finitiProblema di convezione-diffusione alle

differenze finiteStabilizzazione

Autovalori esatti e Autovalori esatti e approssimatiapprossimati

Per studiare l’andamento degli autovalori esatti del problema del filo elastico e degli autovalori approssimati, uso questa function:

function [landa_exa,landa_disc]=autovalori_1d(n)% Calcola i primi autovalori esatti del problema del % filo elastico, e i corrispondenti % autovalori del problema discretoh=1/(n+1);for l=1:n landa_exa(l)=l^2*pi^2; landa_disc(l)=(2-2*cos(l*pi*h));endlanda_disc=landa_disc/h^2;

Ottengo questi risultatiN = 19

N è il numero di nodi interni

N = 49

N. B.: Solo i primi autovalori sono approssimati con precisione

In due dimensioni l’andamento degli autovalori è simile

N = 9

N è il numero di nodi interni per lato

AliasingAliasingIn questa sezione vorrei illustrare la rappresentazione delle autofunzioni del problema del filo elastico su una griglia.Considereremo il caso di una funzione che può essere risolta su una determinata griglia e il caso di una autofunzione troppo oscillante per essere individuata correttamente sulla griglia assegnata.

Su una griglia con 9 nodi interni, quindi con 9 autovettori linearmente indipendenti, consideriamo le seguenti funzioni:u(x) = sin(7 x)u(x) = sin(10 x)u(x) = sin(11 x)

u(x) = sin(7 x)

Il grafico di u(x) è:La griglia “vede” questi dati:

Questi punti individuano l’autovettore sin( 7 π x j ), j = 1,…,9

u(x) = sin(10 x)

Il grafico di questa funzione è:

La griglia “vede” questi dati:Quindi, su questa griglia, la funzione u(x) = sin(10 x)è equivalente alla funzione nulla

u(x) = sin(11 x)

Il grafico di questa funzione è:

La griglia “vede” questi valori:

Questi valori però sono gli stessi che ottengo su questa griglia con la funzione u(x) = sin(9 π x):

La griglia quindi non distingue fra queste due funzioni

La nostra griglia non distingue la funzione verde dalla funzione blu

u(x) = sin(14 x)

La funzione blu ha gli stessi valori sulla griglia della funzione verde. Il sistema non distingue le due funzioni.

Membrana elasticaMembrana elastica

Vogliamo risolvere il problema:

Abbiamo già calcolato la matrice di rigidità per questo problema.Dobbiamo creare una function che abbia in input il numero N dei nodi interni ad ogni lato, e in output la matrice di rigidità.

Costruzione della matriceCostruzione della matrice

Costruisco la matrice di rigidità come matrice tridiagonale a blocchi

Se N è il numero di punti interni di ogni lato, allora la matrice A h è una matrice N 2 per N 2, e i singoli blocchi sono N per N.

I blocchi I sono matrici identità N per N.I blocchi G sono tridiagonali ed hanno la seguente struttura:

Inizio dai blocchi sulla diagonale principale. Siccome ogni blocco occupa N righe, l’i-esimo blocco comincia alla riga:

inizio=(i-1)*n+1;

e finisce alla riga

fine=i*n;

function a=laplaciano(n)% A=LAPLACIANO(N) calcola la matrice del laplaciano% sul quadrato, con N nodi interni su ogni lato.b=ones(n-1,1);g=4*eye(n)-diag(b,-1)-diag(b,1);% g contiene i blocchi sulla diagonalefor i=1:n inizio=(i-1)*n+1; fine=i*n; a(inizio:fine, inizio:fine)=g;end% costruisce le due diagonali lontaneb=ones(n^2-n,1);a=a -diag(b,n) -diag(b,-n);

Metodo FEM in 2DMetodo FEM in 2DA questo punto, posso creare una function che risolva il problema agli elementi finiti della membrana elastica. Questa function deve:

Creare la matrice;Creare il vettore di carico (per ora considero

un carico uniforme)Risolvere il sistema lineare

Ottengo questo programma

function uquad=membrana_unif(n)% UQUAD=membrana(N) trova la soluzione del% problema della membrana elastica, sul % quadrato unitario, con N% nodi interni per lato, con un carico % uniformeh = 1/(n+1); a=laplaciano(n);b= -h^2*ones(n^2,1);u=a\b;

…continua...

Attenzione!

Il vettore soluzione u è un unico vettore colonna con N2 componenti.

Per visualizzare la soluzione, devo riordinare i valori contenuti nel vettore u distribuendoli su un array bidimensionale

Visualizzare la soluzioneVisualizzare la soluzione

Trasformare il vettore soluzione u in un array bidimensionale che dia la soluzione in ogni punto del quadrato

Aggiungere le condizioni al contorno.Fare un grafico con il comando mesh

% scrive u come array bidimensionale, a partire da (1,1)for i=1:n inizio=(i-1)*n+1; fine=i*n; uu(:,i)=u(inizio:fine);end% aggiunge le condizioni al bordouquad(2:n+1,2:n+1)=uu; % aggiunge la cornicefor i=1:n+2 uquad(n+2,i)=0; uquad(i,n+2)=0;end

Quindi la function membrana_unif deve contenere anche queste istruzioni

Per eseguire questo programma e visualizzare i risultati, devo dare i seguenti comandi:

>> uu=membrana_unif(9);>> mesh(uu)

Ottengo la figura seguente:

Matlab stampa il grafico in funzione degli indici della matrice

Per avere le scale corrette sugli assi devo preparare dei vettori che contengano le coordinate dei nodi:

>> x=linspace(0,1,21);>> y=linspace(0,1,21);>> mesh(x,y,u)

N.B.: ci sono 21 nodi per lato, perché devo considerare anche i nodi sui bordi

Vorrei ora poter definire un carico più generale. Devo costruire una funzione f(x,y), che mi dia i valori di f sul quadrato

>> f=inline('x.^2-y.^2');>> x=linspace(-1,1,21);>> y=linspace(-1:1,21);

Se pero’ora calcolo f(x,y) ottengo un vettore, perche’ Matlab valuta f(x(i),y(i)), per i che va da 1 a length(x): quindi otterrei un vettore con 21 componenti

Per ottenere la funzione sul quadrato devo dare i seguenti comandi

>> f=inline('x.^2-y.^2');>> x=linspace(-1,1,21);>> y=linspace(-1:1,21);>> [xx,yy]=meshgrid(x,y);>> fxy=f(xx,yy);

L’ istruzione meshgrid crea la matrice xx e la matrice yy che contengono le ascisse e le ordinate di ogni punto della griglia

Posso ora dare la funzione che calcola la deformazione della membrana, con un carico assegnato

function uquad=membrana(f,n)% UQUAD=membrana(N) trova la soluzione del% problema della membrana elastica, sul % quadrato unitario, con N% nodi interni per lato, con un carico % assegnato tramite la funzione f(x,y)a=laplaciano(n); h=1/(n+1);% Calcola il carico sui nodi interni del quadrato:x=linspace(0+h,1-h,n); y=linspace(0+h,1-h,n); [xx,yy]=meshgrid(x,y);fxy=feval(f,xx,yy);

% Trasforma il carico come vettore colonnafor i=1:n inizio=(i-1)*n+1; fine=i*n; b(inizio:fine)=h^2*fxy(i,:);endu=a\b';

Per impostare il sistema lineare, devo costruire il vettore dei termini noti, impostando il carico come vettore colonna, e moltiplicando per h 2

% scrive u come array bidimensionale, a partire da (2,2)for i=1:n inizio=(i-1)*n+1; fine=i*n; uu(:,i)=u(inizio:fine);end% aggiunge le condizioni al bordouquad(2:n+1,2:n+1)=uu; % aggiunge la cornicefor i=1:n+2 uquad(n+2,i)=0; uquad(i,n+2)=0;end

Per visualizzare la soluzione devo procedere come prima, trasformando u nella matrice uu, definita sul quadrato, e aggiungendo le condizioni al bordo.

Ottengo questa figura, per 19 punti interni

Se ora desidero ottenere più dettagli, devo aumentare il numero di nodi su ogni lato del quadrato.Se però dò questo comando:

pianto il computer (o almeno il mio PC va in crisi…)Cercheremo prima di capire che cosa succede. Poi, nel capitolo sui metodi iterativi, vedremo come è possibile migliorare le prestazioni del nostro programma.

>> u=membrana(f,100);

EserciziEserciziCalcolare il tempo di esecuzione del programma della membrana per N=10, N=20, N=30, N=40, N=50. Con che velocità aumenta il tempo di esecuzione, rispetto ad N? (Usare polyfit per stimare l’andamento del tempo di calcolo)

Modificare il programma membrana.m, in modo da assegnare condizioni al bordo di Neumann omogenee su uno dei lati del quadrato

Problema di convezione -Problema di convezione -diffusionediffusione

Vogliamo risolvere il problema:

Nel seguito considereremo sempre un carico uniforme, f(x,y)=1.

Qui β e’ il vettore velocita’, che considereremo costante. Chiamo (a,b) le componenti di β.

Metodo agli elementi finitiMetodo agli elementi finitiSe uso elementi finiti, con elementi P1 su una triangolazione uniforme, ottengo una matrice con questa struttura:

Matrice di rigiditàMatrice di rigidità

Costruisco la matrice di rigidità come matrice tridiagonale a blocchi

Se N è il numero di punti interni di ogni lato, allora la matrice A h è una matrice N 2 per N 2, e i singoli blocchi sono N per N.

I blocchi G sono tridiagonali ed hanno la seguente struttura:

I numeri di Peclet sono definiti dalle relazioni

e:

dove:

I blocchi B sono bidiagonali ed hanno la seguente struttura.Sopra la diagonale principale i blocchi sono:

mentre sotto la diagonale principale abbiamo:

functionfunction mat_cd2d

function a=mat_cd2d(n,nu,beta)% A=MAT_CD2D(N,NU,[A,B]) calcola la matrice FEM % di convezione-diffusione sul quadrato, % con N nodi interni su ogni lato.% Nu e' la diffusione, BETA e' il vettore velocita'.betax=beta(1); betay=beta(2);h=1/(n+1); %ampiezza di grigliapex=betax*h/(2*nu) %numero di Peclet lungo xpey=betay*h/(2*nu) %numero di Peclet lungo ypex=pex/3; pey=pey/3;

% blocco diagonaleb1=(-1+pex+2*pey)*ones(n-1,1);b2=(-1-pex-2*pey)*ones(n-1,1);g=4*eye(n)+diag(b2,-1)+diag(b1,1);% g contiene i blocchi sulla diagonalefor i=1:n inizio=(i-1)*n+1; fine=i*n; a(inizio:fine, inizio:fine)=g;end% costruisce le quattro diagonali lontaneb1=(-1+2*pex+pey)*ones(n^2-n,1);b2=(-1-2*pex-pey)*ones(n^2-n,1);a=a +diag(b1,n) +diag(b2,-n);c=(pex-pey)*ones(n^2-n+1,1); c(1:n:n^2-n+1,1)=0;a=a+diag(c,n-1)-diag(c,-n+1);a=nu*a;

Struttura del programmaStruttura del programma

Costruisce la matrice di rigidità.Costruisce il vettore di carico.Risolve il sistema lineare.Memorizza la soluzione su una matrice per

permetterne la visualizzazione

La function che calcola la soluzione del problema di convezione-diffusione ha la seguente struttura

FEM e Differenze FiniteFEM e Differenze Finite

FEM : utilizza elementi P1 agli elementi finiti. DIF : utilizza differenze finite centrate UPW : è basato su differenze finite upwind SUP : costruisce il metodo agli elementi finiti

SUPG stabilizzato

Nel caso di carico uniforme e griglia uniforme, il metodo FEM e il metodo alle differenze finite si distinguono solo per la struttura diversa della matrice di rigidità. Uso quindi lo stesso programma per provare metodi diversi. In particolare:

functionfunction membrana_cd membrana_cd

function uquad=membrana_cd(n,nu,beta,metodo)% UQUAD=membrana(N) trova la soluzione del% problema di convezione diffusione, sul % quadrato unitario, con N nodi interni per % lato, con un carico uniforme.% METODO='FEM' (default) utilizza la matrice % degli elementi finiti% METODO='DIF' utilizza la matrice alle% differenze finite centraliif nargin < 4 metodo='FEM'end Scelta di default

h = 1/(n+1);

if metodo=='FEM'

%Metodo FEM

a=mat_cd2d(n,nu,beta);

elseif metodo =='DIF'

%Metodo alle differenze finite

a=mat_cd2d_fd(n,nu,beta);

end

Sceglie il metodo da applicare:

b=h^2*ones(n^2,1);

u=a\b;

% scrive u come array bidimensionale, a partire da (1,1)

for i=1:n

inizio=(i-1)*n+1;

fine=i*n;

uu(:,i)=u(inizio:fine);

end

Risolve il sistema e memorizza la soluzione su un array bidimensionale

% aggiunge le condizioni al bordo

uquad(2:n+1,2:n+1)=uu; % aggiunge la cornice

for i=1:n+2

uquad(n+2,i)=0;

uquad(i,n+2)=0;

end

Aggiunge le condizioni di Dirichlet omogenee al bordo:

scriptscript script_2dcd script_2dcdInfine, uso questo script per lanciare il programma:

% Calcola la soluzione del problema di convezione% diffusione e disegna la soluzione% METODO='FEM' (default) utilizza la matrice % degli elementi finiti% METODO='DIF' utilizza la matrice alle% differenze finite centralin=19; nu=0.05; beta=[-1,1];metodo='FEM'uquad=membrana_cd(n,nu,beta,metodo);x=linspace(0,1,n+2);y=linspace(0,1,n+2);mesh(x,y,uquad)

Ottengo questo grafico:

Direzione del vento

Cambiando la direzione del vento:

Diminuendo nu aumenta la ripidità della soluzione nello strato limite

I dati sono: n=19; nu=0.03; beta=[1,-1];

I numeri di Peclet sono:

pex = 0.8333pey = -0.8333

Diminuendo nu ancora, la soluzione comincia ad oscillare:

EserciziEsercizi

Studiare il comportamento della soluzione in funzione del numero di Peclet. La comparsa delle oscillazioni spurie dipende dalla direzione del vento? Posso stabilire in base ad un unico parametro Pe = Pe (Pex, Pey) se ci saranno oscillazioni spurie?

Studiare lo spessore dello strato limite in funzione di nu e di h, con β=(-1,0), calcolando la larghezza della regione in cui la soluzione varia dal suo valore massimo a zero.

Differenze FiniteDifferenze Finite

I comandi:

n=19; nu=0.05; beta=[-1,1];metodo=‘DIF'

lanciano la soluzione del problema di convezione-diffusione con il metodo delle differenze finite centrate. Qualitativamente, ottengo la stessa soluzione che avevo calcolato con gli stessi dati ed il metodo FEM.

La matrice di convezione-diffusione alle differenze finite ha una struttura diversa dalla matrice di rigidità degli elementi finiti:

Infatti i blocchi fuori dalla diagonale principale sono diagonali.

Matrice alle differenze finiteMatrice alle differenze finite

La matrice alle differenze finite ha la stessa struttura tridiagonale a blocchi della matrice di rigidita’ agli elementi finiti:

I blocchi lungo la diagonale principale sono ancora tridiagonali:

Questa volta i coefficienti sono dati dalle seguenti formule:

I blocchi sopra e sotto la diagonale principale sono diagonali.I blocchi sopra la diagonale principale sono:

I blocchi sotto la diagonale principale sono:

functionfunction mat_cd2d_fd mat_cd2d_fd

La function che crea la matrice per il metodo alle differenze finite è:

function a=mat_cd2d_fd(n,nu,beta)% A=MAT_CD2D_FD(N,NU,[A,B]) calcola la matrice % alle differenze finite% di convezione-diffusione sul quadrato, % con N nodi interni su ogni lato.% Nu e' la diffusione beta e' il vettore velocita'.betax=beta(1); betay=beta(2);h=1/(n+1); %ampiezza di grigliapex=betax*h/(2*nu) %numero di Peclet lungo xpey=betay*h/(2*nu) %numero di Peclet lungo y

b1=(-1+pey)*ones(n-1,1); b2=(-1-pey)*ones(n-1,1); g=4*eye(n)+diag(b2,-1)+diag(b1,1);% g contiene i blocchi sulla diagonalefor i=1:n inizio=(i-1)*n+1; fine=i*n; a(inizio:fine, inizio:fine)=g;end% costruisce le due diagonali lontaneb1=(-1+pex)*ones(n^2-n,1); b2=(-1-pex)*ones(n^2-n,1); a=a +diag(b1,n) +diag(b2,-n);a=nu*a;

Anche con le differenze finite centrate ottengo oscillazioni spurie, se il numero di Peclet è troppo grande:

n=19; nu=0.01; beta=[-1,1];

Metodo UpwindMetodo Upwind

Quando è molto piccolo rispetto a β, mi aspetto che il compor-tamento della soluzione del problema di convezione diffusione si avvicinerà a quella del problema iperbolico con =0.

I metodi alle differenze finite per problemi iperbolici utilizzano discretizzazioni upwind per la derivata prima. Infatti un metodo basato sulla discretizzazione centrale risulta instabile.

Mi aspetto quindi di poter migliorare la soluzione scegliendo una discretizzazione upwind della derivata prima anche per il problema di convezione diffusione.

Matrice UpwindMatrice UpwindLa matrice alle differenze finite per il problema di convezione-diffusione con metodo upwind ha la stessa struttura della matrice alle differenze finite con differenze centrali

I blocchi sulla diagonale principale sono dati da:

Notare che rispetto alle differenze centrali è aumentato il peso della diagonale principale. Qui (a,b) sono le componenti del vettore .

Inoltre:

I blocchi sulla sopradiagonale e sulla sottodiagonale sono dati rispettivamente da:

Considero il problema: n=19; nu=0.01; beta=[-1,1]

Soluzione alle differenze centrali Soluzione Upwind

n=19; nu=0.01; beta=[-1,1]

Con il metodo upwind, ottengo una soluzione stabile anche se non c’è nessun punto di griglia all’interno dello strato limite.

n=19; nu=0.001; beta=[-1,1];

I due numeri di Peclet valgono circa 25.