FEM -3
G. Puppo
Riassunto

Autovalori e aliasing
 Problema di convezione-diffusione agli
elementi finiti
 Problema di convezione-diffusione alle
differenze finite
 Stabilizzazione
Autovalori esatti e
approssimati
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 discreto
h=1/(n+1);
for l=1:n
landa_exa(l)=l^2*pi^2;
landa_disc(l)=(2-2*cos(l*pi*h));
end
landa_disc=landa_disc/h^2;
Ottengo questi risultati
N = 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
Aliasing
In 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
La griglia
di u(x)
“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 è:
Quindi,
su questi
questa
griglia, la
La griglia
“vede”
dati:
funzione u(x) = sin(10 π x)
è equivalente alla funzione nulla
u(x) = sin(11 π x)
Il grafico di questa funzione è:
La griglia quindi non
Questi
valorifraperò
sono
distingue
queste
due gli stessi
funzioni
che ottengo
su questa
con la
La griglia
“vede”griglia
questi valori:
funzione u(x) = -sin(9 π x):
La nostra griglia “vede” la funzione
verde, più “lenta”, al posto della
funzione blu
u(x) = sin(14 π x)
La funzione blu ha gli stessi valori sulla griglia della funzione
verde, u(x) = -sin(6 π x). Il sistema “vede” solo la funzione
verde.
Problema di convezione diffusione
Vogliamo risolvere il problema:
Qui β e’ il vettore velocita’, che considereremo costante.
Nel seguito considereremo sempre un carico uniforme, f(x,y)=1.
Metodo agli elementi finiti
Se uso elementi finiti, con elementi P1 su una triangolazione
uniforme, ottengo una matrice con questa struttura:
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:
dove:
I numeri di Peclet sono definiti dalle relazioni
e:
I blocchi B sono bidiagonali ed hanno la seguente struttura.
Sopra la diagonale principale i blocchi sono:
mentre sotto la diagonale principale abbiamo:
function 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 griglia
pex=betax*h/(2*nu) %numero di Peclet lungo x
pey=betay*h/(2*nu) %numero di Peclet lungo y
pex=pex/3; pey=pey/3;
% blocco diagonale
b1=(-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 diagonale
for i=1:n
inizio=(i-1)*n+1;
fine=i*n;
a(inizio:fine, inizio:fine)=g;
end
% costruisce le quattro diagonali lontane
b1=(-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 programma
La function che calcola la soluzione del problema di
convezione-diffusione ha la seguente struttura

Costruisce la matrice di rigidità.
 Costruisce il vettore di carico.
 Risolve il sistema lineare.
 Memorizza la soluzione su una matrice per
permetterne la visualizzazione
FEM e Differenze Finite
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:

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
function 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 centrali
if nargin < 4
metodo='FEM'
end
Scelta di default
Sceglie il metodo da applicare:
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
Risolve il sistema e memorizza la soluzione su un array bidimensionale
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
Aggiunge le condizioni di Dirichlet omogenee al bordo:
% 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
script script_2dcd
Infine, 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 centrali
n=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 ν aumenta la ripidità della soluzione nello strato limite
I numeri di
Peclet sono:
I dati sono:
n=19; nu=0.03; beta=[1,-1];
pex =
0.8333
pey =
-0.8333
Diminuendo ν ancora, la soluzione comincia ad oscillare:
Esercizi
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 ν e di h,
con β=(-1,0), calcolando la larghezza della regione in cui la
soluzione varia dal suo valore massimo a zero.
Differenze 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 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:
function 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 griglia
pex=betax*h/(2*nu) %numero di Peclet lungo x
pey=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 diagonale
for i=1:n
inizio=(i-1)*n+1;
fine=i*n;
a(inizio:fine, inizio:fine)=g;
end
% costruisce le due diagonali lontane
b1=(-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 Upwind
Quando υ è molto piccolo rispetto a β, mi aspetto che il comportamento 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 Upwind
La matrice alle differenze finite per il problema di convezionediffusione con metodo upwind ha la stessa struttura della
matrice alle differenze finite con differenze centrali
Il listato della function che calcola la matrice per il metodo upwind é:
function a=mat_cd2d_upw(n,nu,beta)
% A=MAT_CD2D_UPW(N,NU,[A,B]) calcola la matrice
% alle differenze finite con derivate upwind
% 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 griglia
pex=betax*h/(2*nu) %numero di Peclet lungo x
pey=betay*h/(2*nu) %numero di Peclet lungo y
bm=(betay-abs(betay))/2;
bp=(betay+abs(betay))/2;
am=(betax-abs(betax))/2;
ap=(betax+abs(betax))/2;
b1=(-1+bm*h/nu)*ones(n-1,1);
b2=(-1-bp*h/nu)*ones(n-1,1);
g=(4+(abs(betax)+abs(betay))*h/nu)*eye(n)+diag(b2,-1)+diag(b1,1);
% g contiene i blocchi sulla diagonale
for i=1:n
inizio=(i-1)*n+1;
fine=i*n;
a(inizio:fine, inizio:fine)=g;
end
% costruisce le due diagonali lontane
b1=(-1+am*h/nu)*ones(n^2-n,1);
b2=(-1-ap*h/nu)*ones(n^2-n,1);
a=a +diag(b1,n) +diag(b2,-n);
a=nu*a;
Considero il problema:
n=19; nu=0.01; beta=[-1,1]
Soluzione alle differenze centrali
Soluzione Upwind
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.
Metodo SUPG
Il metodo SUPG è derivato dal metodo Upwind.
Si inserisce una correzione stabilizzante nella matrice di rigidità,
che tiene conto della direzione del vento.
Il listato della function che calcola la matrice di rigidità per il
metodo FEM con correzione SUPG é:
function a=mat_cd2d_supg(n,nu,beta)
% A=MAT_CD2D(N,NU,[A,B]) calcola la matrice FEM
% di convezione-diffusione sul quadrato,
% con stabilizzazione SUPG
% 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 griglia
pex=betax*h/(2*nu) %numero di Peclet lungo x
pey=betay*h/(2*nu) %numero di Peclet lungo y
pex=pex/3; pey=pey/3;
Calcolo del coefficiente di stabilizzazione:
% Calcolo del coefficiente tau
pe= norm(beta)*h/(2*nu)
if pe <1
csi=pe;
else
csi=1;
end
tau=csi*h/(2*norm(beta)), tau=tau/nu;
Calcolo del blocco diagonale della matrice di rigidità
% Calcola le correzioni SUPG alle diagonali
supg0=tau*2*(betax^2+betax*betay+betay^2);
supg=-tau*betay*(betax+betay);
b1=(-1+pex+2*pey+supg)*ones(n-1,1);
b2=(-1-pex-2*pey+supg)*ones(n-1,1);
g=(4+supg0)*eye(n)+diag(b2,-1)+diag(b1,1);
% g contiene i blocchi sulla diagonale
for i=1:n
inizio=(i-1)*n+1;
fine=i*n;
a(inizio:fine, inizio:fine)=g;
end
Calcolo dei blocchi fuori della diagonale principale della matrice di
rigidità
% costruisce le quattro diagonali lontane
supg=-tau*betax*(betax+betay);
b1=(-1+2*pex+pey+supg)*ones(n^2-n,1);
b2=(-1-2*pex-pey+supg)*ones(n^2-n,1);
a=a +diag(b1,n) +diag(b2,-n);
supg=tau*betax*betay;
c=(pex-pey)*ones(n^2-n+1,1);
cs=supg*ones(n^2-n+1,1);
c(1:n:n^2-n+1,1)=0;
cs(1:n:n^2-n+1,1)=0;
a=a+diag(c+cs,n-1)+diag(-c+cs,-n+1);
a=nu*a;
Considero il problema:
Soluzione FEM
n=19; nu=0.01; beta=[-1,1]
Soluzione FEM stabilizzato
Se abbasso nu:
Soluzione FEM
n=19; nu=0.01; beta=[-1,1]
Soluzione SUPG
Questa volta, la stabilizzazione SUPG non é sufficiente