FEM -2

1. FEM -2 Gabriella Puppo

2. Indice • Problema della trave elastica • Visualizzazione della soluzione • Problema della membrana elastica

3. 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

4. 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)

5. function mat_trave.m Crea le righe relative alle funzioni di base per i valori nodali della soluzione… 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

6. 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.

7. 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

8. Numero di condizionamento Il numero di condizionamento del problema della trave è molto più grande di quello del filo elastico

9. 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

10. Soluzione del problema della trave • Costruire la matrice di rigidità • Costruire il vettore di carico • Aggiungere le condizioni al bordo • Esportare la soluzione La function che calcola la soluzione del problema della trave deve:

11. 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

12. 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));

13. 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

14. Trave con carico uniforme Rosso: N=10, curva blu N=20 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) Ho usato i seguenti comandi:

15. Trave con carico puntiforme

16. 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

17. 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

18. 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

19. Trave con carico puntiforme

20. 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: con soluzione esatta: f(x)=2*exp(-x^2)*(x^2-6*x+1-7*x^4+14*x^3+2*x^6-4*x^5) exa=x^2*(x-1)^2*exp(-x^2)

21. 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à.

22. 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.

23. 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;

24. 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);

25. 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

26. 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...

27. 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

28. 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

29. 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

30. >> uu=membrana_unif(9); >> mesh(uu) Per eseguire questo programma e visualizzare i risultati, devo dare i seguenti comandi: Ottengo la figura seguente:

31. Matlab stampa il grafico in funzione degli indici della matrice

32. 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

33. >> f=inline('x.^2-y.^2'); >> x=linspace(-1,1,21); >> y=linspace(-1:1,21); Vorrei ora poter definire un carico più generale. Devo costruire una funzione f(x,y), che mi dia i valori di f sul quadrato 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

34. >> 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); Per ottenere la funzione sul quadrato devo dare i seguenti comandi L’ istruzione meshgrid crea la matrice xx e la matrice yy che contengono le ascisse e le ordinate di ogni punto della griglia

35. 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);

36. 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';

37. 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

38. Ottengo questa figura, per 19 punti interni

39. 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.

40. 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