FEM -2
Gabriella Puppo
Indice

Problema della trave elastica
 Visualizzazione della soluzione
 Problema della membrana elastica
Trave elastica
Vogliamo risolvere il problema:
Abbiamo già calcolato la matrice di rigidità di questo problema.
La struttura della matrice dipende dalla numerazione che scegliamo
per gli elementi della base dello spazio delle funzioni test
Matrice di rigidità
Costruisco la matrice di rigidità per una griglia con N nodi
interni numerando le funzioni di base in questo modo:
Φ(i), i=1,…,N
le funzioni di base che trasportano il valore
della soluzione nel nodo x(i)
Φ(N+i), i=1,…,N le funzioni di base che trasportano il valore
della derivata prima della soluzione nel
nodo x(i)
function mat_trave.m
function a=mat_trave(n)
% Crea la matrice di rigidita' della trave
% elastica, con N nodi interni
a=eye(2*n);
for i=1:n
a(i,i)=24;
end
for i=1:n-1
a(i,i+1)=-12;
a(i,i+n+1)=6;
end
Crea le righe
relative alle
funzioni di base per
i valori nodali della
soluzione…
for i=2:n
a(i,i-1)=-12;
a(i,i+n-1)=-6;
end
for i=n+1:2*n
a(i,i)=8;
end
for i=n+1:2*n-1
a(i,i+1)=2;
a(i,i-n+1)=-6;
end
for i=n+2:2*n
a(i,i-1)=2;
a(i,i-n-1)=6;
end
Crea le righe relative alle funzioni
di base per i valori nodali della
derivata prima della soluzione.
Per N=5 nodi interni, ottengo la seguente matrice
24
-12
0
0
0
0
6
0
0
0
-12
24
-12
0
0
-6
0
6
0
0
0
-12
24
-12
0
0
-6
0
6
0
0
0
-12
24
-12
0
0
-6
0
6
0
0
0
-12
24
0
0
0
-6
0
0
-6
0
0
0
8
2
0
0
0
6
0
-6
0
0
2
8
2
0
0
0
6
0
-6
0
0
2
8
2
0
0
0
6
0
-6
0
0
2
8
2
0
0
0
6
0
0
0
0
2
8
Numero di condizionamento
Il numero di condizionamento
del problema della trave è molto
più grande di quello del filo
elastico
Andamento
Il numero di condizionamento per il
problema della trave elastica cresce
come N 4. La curva a pallini è stata
ottenuta interpolando i dati sul
condizionamento con un polinomio
di grado 4, con il metodo dei minimi
quadrati
Soluzione del problema della
trave
La function che calcola la soluzione del problema della trave deve:

Costruire la matrice di rigidità
 Costruire il vettore di carico
 Aggiungere le condizioni al bordo
 Esportare la soluzione
function trave.m
Blocco di inizializzazione: calcolo della matrice di rigidità
function [uu,uux,x]=trave(f,n)
% Risolve il problema della trave elastica
% per condizioni al bordo omogenee di
% Dirichlet, con carico assegnato f
% e N punti interni
a=mat_trave(n);
h=1/(n+1); % spaziatura di griglia
Calcolo del vettore di carico: uso la regola di Simpson per stimare
gli integrali. Con questa regola di quadratura devo valutare f
anche sui nodi intermedi x i +h/2 (nodi sfalsati).
x=linspace(h,1-h,n); %nodi interni
xhalf=linspace(h/2,1-h/2,n+1); %nodi sfalsati
fx=feval(f,x);
%n componenti
fxhalf=feval(f,xhalf); %n+1 componenti
b=h^4/3*(fx+fxhalf(1:n)+fxhalf(2:n+1));
b(n+1:2*n)=h^4/12*(fxhalf(1:n)-fxhalf(2:n+1));
Soluzione ed esportazione dei risultati.
In output questa function fornisce u e u’ nei due vettori separati
uu e uux, con le rispettive condizioni al contorno.
u=a\b';
% Ora aggiungo le condizioni al contorno
uu(2:n+1)=u(1:n); uu(1)=0; uu(n+2)=0;
uux(2:n+1)=u(n+1:2*n); uux(1)=0; uux(n+2)=0;
% Calcola il vettore delle ascisse
for i=0:n+1; x(i+1) = i*h; end
Trave con carico uniforme
Rosso: N=10, curva blu N=20
Ho usato i seguenti comandi:
f=inline('-1+0*x');
[uu,uux,x]=trave(f,10);
subplot(1,2,1)
plot(x,uu,'r','Linewidth',2)
[uu,uux,x]=trave(f,20);
subplot(1,2,2)
plot(x,uu,'b','Linewidth',2)
Trave con carico puntiforme
Visualizzazione dei risultati
Nei grafici precedenti, ho utilizzato soltanto i valori nodali della
soluzione per tracciare l’andamento della deformazione della
trave.
In realtà ho a disposizione un’informazione molto più ricca.
Per sfruttare in modo più efficace la soluzione ottenuta, devo
calcolare l’interpolante cubico a tratti di Hermite che ho usato per
costruire il metodo agli elementi finiti della trave elastica
function visualizza_trave.m
function [uu,xx]=visualizza_trave(u,ux,x)
% [UU,XX]=VISUALIZZA_TRAVE(U,UX,X) calcola
% i valori UU dell'interpolante cubico a tratti
% di Hermite sulla griglia XX, utilizzando i
% valori U di una funzione e i valori della sua
% derivata UX sui nodi X.
% La griglia XX e' costruita infittendo X
%
n=length(u)-1; %Numero intervalli
h=1/n;
m=10; hh=h/m; % Raffinamento di griglia
xx(1)=x(1); uu(1)=u(1);
for i=1:n
for jj=1:m
ii=m*(i-1)+jj+1;
y=jj*hh/h;
xx(ii)=x(i)+jj*hh;
phi1=1-y^2+2*dx^2*(dx-1);
phi2=dx^2-2*dx^2*(dx-1);
phi3=dx-dx^2+dx^2*(dx-1);
phi4=dx^2*(dx-1);
uu(ii)=u(i)*phi1 +u(i+1)*phi2+ ...
ux(i)*phi3 +ux(i+1)*phi4;
end
end
Calcolo delle
funzioni di base
di elemento
Calcolo
dell’interpolante
sull’ elemento i,
nel punto y
Trave con carico puntiforme
Esercizi
Paragonare la deformazione della trave e del filo elastico sotto
l’azione dello stesso carico.
Studiare la velocita’ di convergenza del metodo FEM per una
trave soggetta al carico:
f(x)=2*exp(-x^2)*(x^2-6*x+1-7*x^4+14*x^3+2*x^6-4*x^5)
con soluzione esatta:
exa=x^2*(x-1)^2*exp(-x^2)
Membrana 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 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 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
b=ones(n^2-n,1);
a=a -diag(b,n) -diag(b,-n);
Metodo FEM in 2D
A 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
% uniforme
h = 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 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
Quindi la function membrana_unif deve contenere anche queste
istruzioni
% 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 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
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);
Per impostare il sistema lineare, devo costruire il vettore dei
termini noti, impostando il carico come vettore colonna, e
moltiplicando per h 2
% Trasforma il carico come vettore colonna
for i=1:n
inizio=(i-1)*n+1;
fine=i*n;
b(inizio:fine)=h^2*fxy(i,:);
end
u=a\b';
Per visualizzare la soluzione devo procedere come prima,
trasformando u nella matrice uu, definita sul quadrato, e
aggiungendo le condizioni al bordo.
% 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 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
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:
>> u=membrana(f,100);
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.
Esercizi
Calcolare 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