1 / 51

FEM -3

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.

nuru
Download Presentation

FEM -3

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. FEM -3 G. Puppo

  2. Riassunto • Autovalori e aliasing • Problema di convezione-diffusione agli elementi finiti • Problema di convezione-diffusione alle differenze finite • Stabilizzazione

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

  4. Ottengo questi risultati N = 19 N è il numero di nodi interni

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

  6. In due dimensioni l’andamento degli autovalori è simile N = 9 N è il numero di nodi interni per lato

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

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

  9. u(x) = sin(10 π x) Il grafico di questa funzione è: Quindi, su questa griglia, la funzioneu(x) = sin(10 π x)è equivalente alla funzione nulla La griglia “vede” questi dati:

  10. u(x) = sin(11 π x) Il grafico di questa funzione è: La griglia quindi non distingue fra queste due funzioni Questi valori però sono gli stessi che ottengo su questa griglia con la funzione u(x) = -sin(9 π x): La griglia “vede” questi valori: La nostra griglia “vede” la funzione verde, più “lenta”, al posto della funzione blu

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

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

  13. Metodo agli elementi finiti Se uso elementi finiti, con elementi P1 su una triangolazione uniforme, ottengo una matrice con questa struttura:

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

  15. I blocchi G sono tridiagonali ed hanno la seguente struttura: dove: I numeri di Peclet sono definiti dalle relazioni e:

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

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

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

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

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

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

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

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

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

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

  26. Ottengo questo grafico: Direzione del vento

  27. Cambiando la direzione del vento:

  28. Diminuendo ν aumenta la ripidità della soluzione nello strato limite I numeri di Peclet sono: pex = 0.8333 pey = -0.8333 I dati sono: n=19; nu=0.03; beta=[1,-1];

  29. Diminuendo ν ancora, la soluzione comincia ad oscillare:

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

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

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

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

  34. I blocchi lungo la diagonale principale sono ancora tridiagonali: Questa volta i coefficienti sono dati dalle seguenti formule:

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

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

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

  38. Anche con le differenze finite centrate ottengo oscillazioni spurie, se il numero di Peclet è troppo grande: n=19; nu=0.01; beta=[-1,1];

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

  40. Matrice Upwind La 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

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

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

  43. Considero il problema: n=19; nu=0.01; beta=[-1,1] n=19; nu=0.01; beta=[-1,1] Soluzione Upwind Soluzione alle differenze centrali

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

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

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

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

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

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

  50. Considero il problema: n=19; nu=0.01; beta=[-1,1] n=19; nu=0.01; beta=[-1,1] Soluzione FEM stabilizzato Soluzione FEM

More Related