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 le
3 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:
>> dl1=diff(l1)
dl1 =
4*x-3
Moltiplicare la derivata per se stessa e integrare fra 0 e 1:
>> 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'*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
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 0
0 0 0
0 0 0
0 0 0
-8 1 0
16 -8 0
-8 14 -8
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 ...
… 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 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 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 3h a destra
u=a\b;
% 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_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
Calcolo l’ errore
usando griglie via 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
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
Risultati
Ottengo questi risultati:
N
10
20
40
P1
1.784725544444921e-03
4.910462530513526e-04
1.288307737212224e-04
P2
2.082890535519071e-05
1.844384948512536e-06
1.630962657862395e-07
Qui N è il numero di intervalli. Il metodo agli elementi P2
usa 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 intervalli
CP2(N) = numero di condizionamento matrice P2 con N intervalli
Disegno 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, tipo
f(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 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_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; end
b=h^2*feval(f,x); b=b'; % ho portato h a sinistra
u=a\b;
% 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?