Rapport Implementation de La MEF 2D

2016/2017

RAPPORT PROGRAMME MEF 2D (P1/P2/Q1/Q2) Réalisé par : 

MARTAB Rochdi



NAFAE Soukaina

Encadré par: M. A.RACHID

ENSAM Casablanca Implémentation de la MEF

Sommaire Introduction ................................................................... ........................................................................................ ..................... 2 Maillage : ............................................................ ............................................................................................. ................................. 3 Fonction chapeau : ......................................................... .............................................................................. ..................... 6 Quadrature .......................................................................................... .......................................................................................... 8 Jacobien : ........................................................... .......................................................................................... ............................... 10 MatElem2 .......................................................................................... .......................................................................................... 13 SMelem .............................................................. ............................................................................................. ............................... 15 Assemblage .................................................................... ....................................................................................... ................... 16 MAIN 2D ............................................................. ............................................................................................ ............................... 18 Conclusion ......................................................................................... ......................................................................................... 21

ENSAM Casablanca 2016/2017

Introduction La méthode des éléments- finis (MEF) est une méthode d’ approximation numérique de solutions de problèmes aux limites statiques ou dynamiques tels que – diffusion diffusion thermique – mécanique des milieux continus (solides et fluides) – électromagnétisme

mais en fait,

absolument tous les problèmes d’équations aux dérivées partielles (EDP) aux limites.

Il

s’agit, comme dans to utes les méthodes numériques, de trouver une approximation discrète. Pour faire bref, d’un problème différentiels aux limites linéaire, on trouve une formulation

vibrationnelle associe équivalente, dont on calcule une approximation de la solution en projetant sur un espace de dimension finie, ce qui revient a résoudre au final un système linéaire.

L’appellation éléments finis vient de la décomposition du domaine d’étude en

éléments : ils sont souvent représentés par un maillage. Historiquement, l’origine de la méthode peut se trouver dans les travaux de Fermat et

Bernoulli ´ (1743) avec le calcul des variations, puis il faut attendre le début du XX émet siècle avec les progrès en analyse avec la méthode de Galérien se basant sur des théorèmes de projection dans les espaces de Hilbert. En 1943 Robert Courant introduit le principe vibrationnel avec des fonctions de base a supports locaux ouvrant la voie à une division d’un domaine considéré en ” éléments”. Cependant ce n’est qu’avec le développement des

ordinateurs que ces travaux trouvent leurs applications avec les travaux pionniers de Zienckiewiz et Argyris qui définiront la méthode en 1960. Ce qui amène le succès de la méthode et sa puissance est l’apport du calcul matriciel, introduit par un ingénieur civil anonyme. La méthode connait alors un développement développement fulgurant accompagne par les progrès de l’informatique. La méthode des éléments-finis est une méthode puissante basée sur une théorie mathématique rigoureuse. Aujourd’hui, les éléments-

finis sont un outil majeur, incontournable en mécanique (fluides et solides, interactions, structures), et applicable dans de nombreux domaines impliquant des problèmes d’EDP aux limites comme par exemple en mathématiques financières ou l’électromagnétisme.

De

nombreux codes industriels (solveurs) existent et sont généralement couplés à un logiciel de CAO 1 ou Computer Aided Design (CAD) en Anglais. Citons Ansys, Abaqus, Robot, LSdyna, Feap, Code- Aster, Cast3M et bien d’autres.


On considère le problème suivant :

 ∗∆ ∗ ∆ = , , u(x,y)=0 ∀ ∈   Maillage : Pour le maillage, On doit comme dans le cas de 1D discrétiser notre domaine, il y a deux types de discrétisation, soit triangulaire ou quadrilatère. On commence tout d’abord par la table de coordonnes, o n la définit pour les deux types de

maillage. Tandis que pour la table de connectivite, ça diffère selon les cas, pour P1 la table se compose de trois colonnes, qui contiennent les positions des nœuds dans notre élément, cette table se remplit à l’aide d’une boucle et contient 2*nx*ny lignes (nx=npx-1, npx est le nombre de

point), Pour P2 la table contient 6 colonnes. Dans le cas de quadrilatère Q1 la table se compose de 4 colonnes et pour Q2 de 9 colonnes. Dans tous les cas le nombre de lignes est le même que P1.

function [X,T2,b,T] function [X,T2,b,T] = Maillage2d(x1,x2,y1,y2,npx,npy,elem,k) T2=0; T=0; if k==1 if k==1 X = zeros((npx)*(npy),2); zeros((npx)*(npy),2); % Nombre des ï¿½lï¿½ments dans chaque direction nx = npx-1; ny = npy-1; % hx = (x2-x1)/nx; hy = (y2-y1)/ny; xs = linspace(x1,x2,npx)'; unos = ones(npx,1); % % Coordonnï¿½es des noeuds % yys = linspace(y1,y2,npy); linspace(y1,y2,npy); for i=1:npy for i=1:npy ys = yys(i)*unos; posi = (i-1)*(npx)+1:i*(npx); X(posi,:)=[xs,ys]; end % P1 % Connectivit Connectivitï¿½ ï¿½ if(elem==1) if (elem==1) T2 = zeros(nx*ny,3); for i=1:ny for i=1:ny for j=1:nx for j=1:nx ielem = 2*((i-1)*nx+j)-1;


inode = (i-1)*(npx)+j; T2(ielem,:) = [inode inode+1 inode+(npx)]; T2(ielem+1,:) = [inode+1+npx inode+npx inode+1]; end end else if if(elem==0) (elem==0) T2 = zeros(nx*ny,4); zeros(nx*ny,4); ielem = 0; for i=1:ny for i=1:ny for j=1:nx for j=1:nx ielem = ielem+1; inode = (i-1)*(npx)+j; T2(ielem,:) T2(ielem,:) = [inode inode+1 inode+1+(npx) inode+1+(npx) inode+(npx)]; end end end end T=T2; %Noeuds sur le bord stockï¿½es dans "b" comme suit (bas,gauche,droite,haut) b=[1:npx,npx+1:npx:npx*npy,2*npx:npx:npx*npy,npx*npy-npx+2:npx*npy-1]; else if if k== k== 2 %%% P2 & Q2 npx2=2*npx-1;npy2=2*npy-1; X = zeros((npx*npy),2); zeros((npx*npy),2); % Nombre des ï¿½lï¿½ments dans chaque direction nx2 = npx2-1; ny2 = npy2-1; nx = npx-1; ny = npy-1; % hx = (x2-x1)/nx2; hy = (y2-y1)/ny2; xs = linspace(x1,x2,npx2)'; unos = ones(npx2,1); % % Coordonnï¿½es des noeuds % yys = linspace(y1,y2,npy2); linspace(y1,y2,npy2); for i=1:npy2 for i=1:npy2 ys = yys(i)*unos; posi = (i-1)*(npx2)+1:i*(npx2); X(posi,:)=[xs,ys]; end %%% P2: % Connectivit Connectivitï¿½ ï¿½ if(elem==1) if (elem==1) T2 = zeros(nx*ny,6); T = zeros(nx2*ny2,3); zeros(nx2*ny2,3); for i=1:ny2 for i=1:ny2 for j=1:nx2 for j=1:nx2 ielem = 2*((i-1)*nx2+j)-1; inode = (i-1)*(npx2)+j; T(ielem,:) = [inode inode+1 inode+(npx2) ]; T(ielem+1,:) = [inode+1+npx2 inode+npx2 inode+1];


end end for i=1:ny for i=1:ny for j=1:nx for j=1:nx ielem = 2*((i-1)*nx+j)-1; inode = 2*(i-1)*(npx2)+(j-1)*2+1; T2(ielem,:) T2(ielem,:) = [inode inode+2 inode+2*npx2 inode+1 inode+npx2+1 inode+(npx2)]; T2(ielem+1,:) = [inode+2*(npx2+1) inode+2*npx2 inode+2 inode+2*npx2 +1 inode+1+npx2 inode+npx2+2] ; %les connectivité de 6 noeuds par element end end else if if(elem==0) (elem==0) %%% Q2: T2 = zeros(nx*ny,9); zeros(nx*ny,9); T = zeros(nx2*ny2,4); ielem = 0; for i=1:ny2 for i=1:ny2 for j=1:nx2 for j=1:nx2 %ielem = 2*((i-1)*nx+j)-1; ielem = ielem+1; inode = (i-1)*(npx2)+j; T(ielem,:) = [inode inode+1 inode+1+(npx2) inode+(npx2)]; end end ielem = 0; for i=1:ny for i=1:ny for j=1:nx for j=1:nx ielem = ielem+1; inode = 2*(i-1)*(npx2)+(j-1)*2+1; T2(ielem,:)= [inode inode+2 inode+2*(npx2+ 1) inode+(2*npx 2) inode+1 inode+2+npx2 inode+1+2*npx2 inode+npx2 inode+1+npx2]; end end end end %Noeuds sur le bord stockï¿½es dans "b" comme suit (bas,gauche,droite,haut) b=[1:npx2,npx2+1:npx2:npx2*npy2,2*npx2:npx2:npx2*npy2,npx2*npy2npx2+2:npx2*npy2-1]; end end end


Fonction chapeau : Cette fonction nous calcule les fonctions de forme dans l’élément de référence pour les différents cas, le calcul des Phis est fait par mapel, on a procède par par la manière suivante pour les calculer : Cas P2 : On a :

(0,1) (0,1/2)

(1/2,1/2)

(0,0)

(1,0) (1/2,0)

On a Et

̂ =  ∗  2 +  ∗∗ 2 +  ∗ ∗  ∗ ∗ +  ∗ + ∗  +  ( , ) =  1 ≤  ≤ 6

Et par ces deux formules, on trouve les coefficients pour les 6 fonctions de forme. Cette fonction nous donne également les gradients des fonctions de formes, on les calcule aussi par Mapel en les dérivant.

function [phi,grad] function [phi,grad] = FoncChap(k,Xgauss,elem) s = Xgauss(:,1); t = Xgauss(:,2); Ngauss=size(Xgauss,1); if elem if elem == 1 if k if k == 1 %P1 phi = [1-(s+t),s,t]; [1-(s+t),s,t]; grad = [-ones(2,1) ,eye(2)]; % N3(0,1) % | \ % | \ % | \ % N6(0,1/2)|_ _ _ \ N5(1/2,1/2) % |\ /\ % |\ / \ % | \ / \ % |_ \/_ _ _ \


% N1(0,0) N4(1/2,0) (1,0)N2 else if if k k == 2 % P2 s1=0;t1=0;s2=1;t2=0;s3=0;t3=1; s4=0.5;t4=0;s5=0.5;t5=0.5;s6=0;t6=0.5; %%% Pour le calcule des phi A=[s1^2 t1^2 s1*t1 s1 t1 1; s2^2 t2^2 s2*t2 s2 t2 1; s3^2 t3^2 s3*t3 s3 t3 1; s4^2 t4^2 s4*t4 s4 t4 1; s5^2 t5^2 s5*t5 s5 t5 1; s6^2 t6^2 s6*t6 s6 t6 1]; coef = [s.^2 t.^2 s.*t s t ones(Ngauss,1)]; phi = coef*(A^-1); grad = [4*s+4*t-3 4*s-1 zeros(Ngauss,1) -8*s-4*t+4 4*t -4*t; 4*s+4*t-3 zeros(Ngauss,1) zeros(Ngauss,1) 4*t-1 -4*s 4*s -8*t-4*s+4]; end end else if if elem elem ==0 if k if k == 1 %Q1 phi = 0.25*[-s-t+ s.*t+1,s-t-s. *t+1,s+t+s.*t+ 1,-s+t-s.*t+ 1]; grad = [-1+t 1-t 1+t -1-t; -1+s -1-s 1+s 1-s]; else if if k k == 2 %Q2 % N4(-1,1) _ _ _N7 _ _ _N3=(1,1) % | | | % | | | % N8 |_ _ _N9 _ _ _|N6 % | | | % | | | % |_ _ _ | _ _ _| % N1=(-1,-1) N5 N2=(1,-1) phi = [0.25*(s.*t).*(1-s).*(1-t),-0.25*(s.*t).*(1+s).*(1t),0.25*(s.*t).*(1+s).*(1+t),-0.25*(s.*t).*(1-s).*(1+t),-0.5*(1+s).*(1s).*t.*(1-t),0.5*s.*(1+s).*(1-t).*(1+t),0.5*(1-s).*(1+s).*t.*(1+t),0.5*s.*(1-s).*(1-t).*(1+t),(1-s).*(1+s).*(1-t).*(1+t)]; % phi1 phi2 phi3 phi4 phi5 phi6 phi7 phi8 phi9 grad = [0.5*s.*(t.^2)-0.5*s.*t-0.25*(t.^2)+0.25*t 0.50*s.*(t.^2)0.5*s.*t+0.25*(t.^2)-0.25*t 0.50*(s.*(t.^2))+0.50*(s.*t)+0.25*(t.^2)+0.25*t 0.50*(s.*(t.^2))+0.50*(s.*t)-0.25*(t.^2)-0.25*t -s.*(t.^2)+s.*t -s.*(t.^2)0.5*(t.^2)+s+0.5 -s.*(t.^2)-s.*t -s.*(t.^2)+0.5*(t.^2)+s-0.5 2*s.*(t.^2)2*s; 0.50*(s.^2).*t-0.25*(s.^2)-0.50*s.*t+0.25*s 0.50*(s.^2).*t0.25*(s.^2)+0.50*s.*t-0.25*s 0.50*(s.^2).*t+0.25*(s.^2)+0.50*(s.*t)+.25*s 0.50*(s.^2).*t+0.25*(s.^2)-0.50*(s.*t)-0.25*s -(s.^2).*t+0.5*(s.^2)+t-0.5 (s.^2).*t-s.*t -(s.^2).*t-0.5*(s.^2)+t+0.5 -(s.^2).*t+s.*t 2*(s.^2).*t2*t]; %Calculés utilisant maple end end %%% end end end


Quadrature Cette fonction nous fournit les points et les poids de gauss selon le nombre de gauss choisi, Pour le cas de triangle on a le choix entre un nombre de gauss de 1, 3, 4, 7, 16 selon le degré de précision voulu, et pour le cas de rectangulaire, on a le choix entre 4 et 9.

function function [Xgauss,Wgauss]=Quadrature(elem, Ngauss) if elem if elem == 0 if Ngauss if Ngauss == 4 % degré 1 pos1 = 1/sqrt(3); Xgauss=[-pos1 -pos1 pos1 -pos1 pos1 pos1 -pos1 pos1]; Wgauss=[ 1 1 1 1]; elseif Ngauss elseif Ngauss == 9 % degré 2 pos1 = sqrt(3/5); Xgauss=[-pos1 -pos1 0 -pos1 pos1 -pos1 -pos1 0 0 0 pos1 0 -pos1 pos1 0 pos1 pos1 pos1]; pg1=5/9; pg2=8/9; pg3=pg1; Wgauss= [pg1*pg1 pg2*pg1 pg3*pg1 pg1*pg2 pg2*pg2 pg3*pg2 pg1*pg3 pg2*pg3 pg3*pg3]; else error(' error(' Cette quadrature n est pas programmée') programmée') end elseif elem elseif elem == 1 if Ngauss if Ngauss == 1 % degré 1 pos1 = 1/3; Xgauss=[pos1 pos1]; Wgauss = 1/2; elseif Ngauss elseif Ngauss == 3 % degré 2 pos1 = 1/2; Xgauss = [pos1 pos1 0 pos1 pos1 0]; pes1 = 1/6; Wgauss = [pes1 pes1 pes1]; elseif Ngauss elseif Ngauss == 4 % degré 3 Xgauss=[1/3 1/3 0.6 0.2 0.2 0.6 0.2 0.2]; Wgauss = [-27/96 25/96 25/96 25/96]; elseif Ngauss elseif Ngauss == 7 % degré 5 a = 0.101286507323456338800987361915123; b = 0.470142064105115089770441209513447; P1 = 0.0629695902724135762978419727500906;


P2 = 0.0661970763942530903688246939165759; Xgauss=[a a a 1-2*a 1-2*a a b b b 1-2*b 1-2*b b 1/3 1/3]; Wgauss = [P1, P1, P1, P2, P2, P2, 0.1125]; elseif Ngauss elseif Ngauss == 12 % degré 6 a = 0.0630890144915022283403316028708191; b = 0.249286745170910421291638553107019; c = 0.0531450498448169473532496716313981; d = 0.310352451033784405416607733956552; P1 = 0.0254224531851034084604684045534344; P2 = 0.0583931378631896830126448056927897; P3 = 0.0414255378091867875967767282102212; Xgauss = [a a a 1-2*a 1-2*a a b b b 1-2*b 1-2*b b c d c 1-c-d 1-c-d c 1-c-d d d 1-c-d d c ]; Wgauss = [P1, P1, P1, P2, P2, P2, P3, P3, P3, P3, P3, P3]; elseif Ngauss elseif Ngauss == 16 % degré 8 a = 0.170569307751760206622293501491464; b = 0.0505472283170309754584235505965989; c = 0.459292588292723156028815514494169; d = 0.728492392955404281241000379176061; e = 0.263112829634638113421785786284643; f = 1/3; P1 = 0.0516086852673591251408957751460645; P2 = 0.0162292488115990401554629641708902; P3 = 0.0475458171336423123969480521942921; P4 = 0.0136151570872174971324223450369544; P5 = 0.0721578038388935841255455552445323; Xgauss = [a a a 1-2*a 1-2*a a b b b 1-2*b 1-2*b b c c c 1-2*c 1-2*c c d e d 1-d-e 1-d-e d 1-d-e e e 1-d-e e d f f ];


Wgauss = [P1*ones(1,3), P2*ones(1,3), P3*ones(1,3), P4*ones(1,6), P5]; else error('Cette error('Cette quadrature n est pas programmée') programmée') end else error(' error(' Cette quadrature n est pas programmée') programmée') end

Jacobien : Cette fonction nous permet de calculer le jacobien, ce dernièr se calcule de la façon suivante :

J = grad grad∗∗     = [ ] Le gradient des phi selon le cas est déjà calculé dans la fonction fonchap. function [J]=jacobien(elem,k,S,x) function % Calcule de la matrice jacobienne % La matrice jacobienne, dont le déterminant nous intéresse, % est définie à partir des dérivées des polynômes de base. % _ _ _ _ % | dphi/ds | | X1 Y1 | % J = | | * | . . | % |_dphi/dt_| |_. . _| % % elem: - espace Q: elem = 0 % - espace P: elem = 1 % k: - degres d'espace d'approximation (Pk ou Qk) % s: - cordoonés du variables s = [Si Ti] par exemple: s = [0 0] % | xi yi | % x: - cordoonés du variables par exemple pour P1: x = [x y] = |xi+1 yi+1| % |xi+2 yi+2| s = S(1,1); t = S(1,2); if elem if elem == 1 if k if k == 1 % N1=(0,1) % | \ % | \ % | \ % | TR \ % |_ _ _\ % N1=(0,0) (1,0)=N2 % _ _ % | x2-x1 x3-x1 | % J = | | % |_y2-y1 y3-y1_| % J = [-1 1 0;-1 0 1]*x; else if if k k == 2 % N3(0,1) % | \ % | \ % | \


% N6(0,1/2)|_ _ _ \ N5(1/2,1/2) % |\ /\ % |\ / \ % | \ / \ % |_ \/_ _ _ \ % N1(0,0) N4(1/2,0) N2(1,0) % grad = [4*s+4*t-3 4*s-1 0 -8*s-4*t+4 4*t -4*t; % 4*s+4*t-3 0 4*t-1 -4*s 4*s -8*t4*s+4]; switch s switch s case 0 case 0 switch t switch t case 0 case 0 J = [-3 -1 0 4 0 0;-3 0 -1 0 0 4]*x; case 1 case 1 J = [1 -1 0 0 4 -4;1 0 3 0 0 -4]*x; case 0.5 case 0.5 J = [-1 -1 0 2 2 -2;-1 0 1 0 0 0]*x; end case 1 case 1 J = [1 3 0 -4 0 0;1 0 -1 -4 4 0]*x; case 0.5 case 0.5 switch t switch t case 0 case 0 J = [-1 1 0 0 0 0;-1 0 -1 -2 2 2]*x;; case 0.5 case 0.5 J = [1 1 0 -2 2 -2;1 0 1 -2 2 -2]*x; end % end end end % else if if elem elem == 0 if k if k == 1 % N4(-1,1) _ _ _ _ _ _ _N3=(1,1) % | | % | | % | QR | % | | % |_ _ _ _ _ _ _| % N1=(-1,-1) N2=(1,-1) switch s switch s case -1 case -1 switch t switch t case -1 case -1 J = ([-2 2 0 0;-2 0 0 2]*x)'; %Le jacobien en l’image d’un sommet de QR (-1,-1) case 1 case 1 J = ([0 0 2 -2;-2 0 0 2]*x)'; %Le jacobien en l’image d’un sommet de QR (-1, 1) end case 1 case 1 switch t switch t case -1 case -1 J = ([-2 2 0 0;0 -2 2 0]*x)'; %Le jacobien en l’image d’un sommet de QR (1,-1) case 1 case 1


J = ([0 0 2 -2;0 -2 2 0]*x)'; sommet de QR (1, 1) end

%Le jacobien en l’image d’un

end else if if k k == 2 % % % % % % % %

N4(-1,1) _ _ _N7 _ _ _N3=(1,1) | | | | | | N8 |_ _ _N9 _ _ _|N6 | | | | | | |_ _ _ | _ _ _| N1=(-1,-1) N5 N2=(1,-1)

switch s switch s case -1 case -1 switch t switch t case -1 case -1 J = ([-1.5 -0.5 0 0 2 0 0 0 0;-1.5 0 0 -0.5 0 0 0 2 0]*x);%Le jacobien en l’image d’un sommet de QR ( -1,-1) case 1 case 1 J = ([0 0 -0.5 -1.5 0 0 2 0 0;0.5 0 0 1.5 0 0 0 -2 0]*x);%Le jacobien en l’image d’un sommet de QR ( -1, 1) case 0 case 0 J = ([0 0 0 0 0 -0.5 0 -1.5 2;-0.5 0 0 0.5 0 0 0 0 0]*x);%Le jacobien en l’image d’un sommet de QR ( -1, 0) end case 1 case 1 switch t switch t case -1 case -1 J = ([0.5 1.5 0 0 -2 0 0 0 0;0 -1.5 -0.5 0 0 2 0 0 0]*x); %Le jacobien en l’image d’un sommet de QR (1, -1) case 1 case 1 J = ([0 0 1.5 0.5 0 0 -2 0 0;0 0.5 1.5 0 0 -2 0 0 0]*x); %Le jacobien en l’image d’un sommet de QR (1, 1)

case 0 case 0 J = ([0 0 0 0 0 1.5 0 0.5 -2;0 -0.5 0.5 0 0 0 0 0 0]*x); %Le jacobien en l’image d’un sommet de QR (1, 0 ) end case 0 case 0 switch t switch t case -1 case -1 J = ([-0.5 0.5 0 0 0 0 0 0 0;0 0 0 0 -1.5 0 -0.5 0 2]*x); %Le jacobien en l’image d’un sommet de QR (0, -1) case 1 case 1 J = ([0 0 0.5 -0.5 0 0 0 0 0;0 0 0 0 0.5 0 1.5 0 -2]*x); %Le jacobien en l’image d’un sommet de QR (0, 1)

case 0 case 0 J = ([0 0 0 0 0 0.5 0 -0.5 0;0 0 0 0 -0.5 0 0.5 0 0]*x); %Le jacobien en l’image d’un sommet de QR (0 , 0) end end end end end end end


MatElem2 Dans cette fonction On construit la matrice élémentaire qui s’écrit de la façon suivante

   = −′′ =    ∗ (− ∗ ) ∗ (− ∗  ) ∗   function function [Ke] = MatElem2(elem,k,X,T,ie,phi,grad,Xgauss,Wgauss) % Ngauss = size(Xgauss,1); if elem if elem == 0 if k if k == 1 Ngeom = 4; elseif k elseif k == 2 Ngeom=9; end elseif elem elseif elem == 1 if k if k == 1 Ngeom = 3; elseif k elseif k == 2 Ngeom = 6; elseif k elseif k == 3 Ngeom = 9; end end if k==1 if k==1 % Noeuds du triangle courant (sommet dans le cas oï¿½ Ngeom=3) x=X(T(ie,1:Ngeom),:); % matice du changement de variable

KKe=zeros(3*k);%initialiser KKe=zeros(3*k); %initialiser les tailles de Ke % % if elem if elem == 0 gradS=grad(1:Ngauss,:); % premiere partie du grad où les grad selon dS gradT=grad((Ngauss+1):(2*Ngauss),:); % deuxiéùe partie du grad où les grad selon dT SS=[-1 -1;1 -1;1 1;-1 1]; for i=1:Ngeom for i=1:Ngeom for j=1:Ngeom for j=1:Ngeom gh = 0; S=SS(j,:); J = jacobien(elem,k,S,x); for lmm=1:Ngauss for lmm=1:Ngauss Jint=[gradS(lmm,:);gradT(lmm,:)]*x; gh = gh + Wgauss(lmm)*((Jint\[gradS(lmm,i);gradT(lmm,i)])'*(Jint\[gradS(lmm,j);gra dT(lmm,j)]));


end KKe(i,j) = det(J)*gh; KKe(j,i)=KKe(i,j); end end Ke = KKe; %calcul de la matrice ï¿½lï¿½mentaire Ke else if if elem elem == 1 SS = [0 0;1 0;0 1]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% for i=1:Ngeom for i=1:Ngeom for j=i:Ngeom for j=i:Ngeom S=SS(j,:); J = jacobien(elem,k,S,x); gh = 0; for lmm=1:Ngauss for lmm=1:Ngauss gh = gh + Wgauss(lmm)*(J\grad(:,i))'*(J\grad(:,j)); % erreur corrigé (calcule d'integrale: faut multiplier % les Wgauss end KKe(i,j) = gh; KKe(j,i)=KKe(i,j); end end Ke=det(J)*KKe; %pas besoin d'utiliser l'air apres changement de variable (air = 1) end end % else if if k k == 2 x=X(T(ie,1:Ngeom),:); KKe=zeros(3*k);%initialiser KKe=zeros(3*k); %initialiser les tailles de Ke % %%% %%% P2 & Q2 gradS=grad(1:Ngauss,:); % premiere partie du grad où les grad selon dS gradT=grad((Ngauss+1):(2*Ngauss),:); % deuxiéùe partie du grad où les grad selon dT if elem if elem == 1 SS = [0 0;1 0;0 1;1/2 0;1/2 1/2;0 1/2]; else if if elem elem ==0 SS = [-1 -1;1 -1;1 1;-1 1;0 -1;1 0;0 1;-1 0;0 0]; end end for i=1:Ngeom for i=1:Ngeom for j=1:Ngeom for j=1:Ngeom S=SS(j,:); J = jacobien(elem,k,S,x); gh = 0; for lmm=1:Ngauss for lmm=1:Ngauss Jint=[gradS(lmm,:);gradT(lmm,:)]*x; gh = gh + Wgauss(lmm)*((Jint\[gradS(lmm,i);gradT(lmm,i)])'*(Jint\[gradS(lmm,j);gra dT(lmm,j)]));


end KKe(i,j) = det(J)*gh; KKe(j,i)=KKe(i,j); end end Ke = KKe; else error(' error(' ce type Pk n est pas encore programmï¿½') programmï¿½') end end end

SMelem On calcule ici le second membre de la façon suivante :

 ∗    ∗  ∗     =  + ∗  ∗   ∫−  ∗  = ∑  function [ Fe ] = SMelem(f,elem,k,X,T,ie,Ngauss) function [ %calcul du second membre élémentaire Fe % [Xgauss,Wgauss]=Quadrature(elem, Ngauss); %fonctions chapeaux [phi, grad] = FoncChap(k,Xgauss,elem); if elem if elem == 0 if k if k == 1 Ngeom = 4; elseif k elseif k == 2 Ngeom=9; end elseif elem elseif elem == 1 if k if k == 1 Ngeom = 3; elseif k elseif k == 2 Ngeom = 6; elseif k elseif k == 3 Ngeom = 9; end end % if k==1 if k==1 Fe=zeros(3*k,1); %initialise %initialiser r les tailles de Me et Fe % Noeuds du triangle courant (sommet dans le cas où Ngeom=3) x=X(T(ie,1:Ngeom),:); grad2

= [-1 1 1 -1; -1 -1 1 1]; if(elem if (elem == 1) J = (grad * x)'; else if if (elem (elem == 0) J = (grad2 * x)'; end end

%


%On calcule Xc=X1+J*S aux noeuds de Gauss dont le nombre est Ngauss Xc=zeros(Ngauss,2); for i=1:Ngauss for i=1:Ngauss Xc(i,:)=x(1,:)+(J*Xgauss(i,:)')'; end % %on calcule Fe1 Fe2 et Fe2 les composantes de Fe %if elem == 1 for i=1:Ngeom for i=1:Ngeom s=0; for j=1:Ngauss for j=1:Ngauss s=s+(Wgauss(j)*f(Xc(j,1),Xc(j,2))*phi(j,i)); end Fe(i)=det(J)*s; end else if if k==2 k==2 %%P2 Fe=zeros(3*k,1); %initialiser les tailles de Me et Fe % Noeuds du triangle courant (sommet dans le cas où Ngeom=3) x=X(T(ie,1:Ngeom),:); if(elem if (elem == 1) grad2 = [-3 -1 0 4 0 0; -3 0 -1 0 0 4]; else if if (elem (elem == 0) grad2 = [zeros(2,4),[0 0.5 0 -0.5 0;-0.5 0 0.5 0 0]]; end end J = (grad2 * x)'; % %On calcule Xc=X1+J*S aux noeuds de Gauss dont le nombre est Ngauss Xc=zeros(Ngauss,2); for i=1:Ngauss for i=1:Ngauss Xc(i,:)=x(1,:)+(J*Xgauss(i,:)')'; end % %on calcule Fe1 Fe2 et Fe2 les composantes de Fe for i=1:Ngeom for i=1:Ngeom s=0; for j=1:Ngauss for j=1:Ngauss s=s+(Wgauss(j)*f(Xc(j,1),Xc(j,2))*phi(j,i)); end Fe(i)=det(J)*s; end end end end

Assemblage Pour cette fonction, On assemble les matrices élémentaires pour construire la matrice globale


function [ K,F ] = Assemblage2D(elem,k,X,T,f,phi,grad,Xgauss,Wgauss,Ngauss) function [ % Assemblage des matrices élémentaires "Ke" dans la matrice globale K % Assemblage des seconds membres élémentaires "Fe" dans le second... % membre global F % On utilise la fonction "MatElem(k,X,T,ie,f,phi,grad,Xgauss,Wgauss)" % %Entrées: % elem: type d'élément (0 pour quadrilatère quadrilatères s et 1 pour triangles) % X: table des coordonnées % T: tbale des connectivités % Nn=size(X,1); %nombre des noeuds Nt=size(T,1); % nombre des éléments % %K=sparse(zeros(Nn,Nn));% initialisation de K comme matrice creuse K=zeros(Nn,Nn);% K=zeros(Nn,Nn); % initialisation de K F=zeros(Nn,1);% F=zeros(Nn,1);% initialisation de F comme matrice creuse % for ie for ie = 1:Nt Tie = T(ie,:); [Ke] = MatElem2(elem,k,X,T,ie,phi,grad,Xgauss,Wgauss); [Fe] = SMelem(f,elem,k,X,T,ie,Ngauss); K(Tie,Tie)=K(Tie,Tie)+Ke; F(Tie)=F(Tie)+Fe; clear Ke Fe Fe; ; end end 1- MEF2D

On impose dans cette fonction les conditions aux bords, et on résout le système K.U=F. function [Kb,Fb,U] = MEF2D(f,elem,k,X,T,b,Ngauss) function [Kb,Fb,U] % fonction traite l'equation alpha*laplace(u)=f sur Omega avec u|Omega=0 % moyennant la méthode des éléments finis de type Galerkin Pk %Entrées % f: fonction second membre % elem: type d'élément (0 pour quadrilatère quadrilatères s et 1 pour triangles) % k: type de la méthodes des éléments finis (Pk) % X: table des coordonnées % T: table de connectivité % b: noeuds sur le bords % Ngauss: Nombre des noeuds de Gauss dans l'élément courant %Sorties: % U: Vecteurs solutions aux noeuds du maillage. % % --------Quadrature de Gauss----------[Xgauss,Wgauss]=Quadrature(elem, Ngauss); % %---------Fonctions chapeaux-----------[phi,grad] = FoncChap(k,Xgauss,elem); % %-------- Assemblage-------------------[ K,F ] = Assemblage2D(elem,k,X,T,f,phi,grad,Xgauss,Wgauss,Ngauss); % %--------Conditions aux bords----------%On met "0" dans les lignes et les colonnes qui correspondent aux noeuds


%bords dans K et F K(b,:)=0; K(:,b)=0; F(b)=0; %On met "1" sur la diagonale de K correspodant aux aux noeuds bords %K(b,b)=speye(length(b %K(b,b)=spe ye(length(b),length(b) ),length(b)); ); % cas où K est sparse K(b,b)=eye(length(b),length(b)); % K et F après conditions aux bords Kb=K; Fb=F; % U=Kb\Fb; end

MAIN 2D Premièrement on choisit une fonction qui vérifie les conditions aux bords, et on construit le problème, puis on fait appel à les fonctions déjà élabore pour le résoudre, et on ajoute la solution exacte pour la comparaison. clear all all, , close all all, , home %-----------------------% Les entrées %-----------------------elem=input('Choisir elem=input('Choisir elment(0 quadrilatère quadrilatère, , 1 triangle) : '); k=input('Choisir k=input('Choisir k pour type Pk : '); Ngauss=input('Choisir Ngauss=input('Choisir le nombre nombre de de points points de Gauss : '); npx=input('Entrer npx=input('Entrer nombre de points suivant (Ox) : '); npy=input('Entrer npy=input('Entrer nombre de points suivant (Oy) : '); %-----------------------------% Les données du maillage x1=-2;x2=2; y1=-1;y2=1; % % Second membre f=@(x,y) -2*x^2-2*y^2+10; % Solution exacte ue=@(x,y) (x^2-4)*(y^2-1); %-------------------------------% maillage [X,T2,b,T] = Maillage2d(x1,x2,y1,y2,npx,npy,elem,k); % calcul de la solution [Kb,Fb,U] = MEF2D(f,elem,k,X,T2,b,Ngauss); %-----------------------------% Affichage et comparaison % Solution approchée figure(1); clf; trisurf(T,X(:,1),X(:,2),0*X(:,1),U,'edgecolor' trisurf(T,X(:,1),X(:,2),0*X(:,1),U, 'edgecolor', ,'k' 'k', ,'facecolor' 'facecolor', ,'interp' 'interp'); ); view(2),axis([x1 x2 y1 y2]),axis equal equal,colorbar; ,colorbar; % Solution exacte figure(2); clf; uee=(X(:,1).^2-4).*(X(:,2).^2-1); trisurf(T,X(:,1),X(:,2),0*X(:,1),uee,'edgecolor' trisurf(T,X(:,1),X(:,2),0*X(:,1),uee, 'edgecolor', ,'k' 'k', ,'facecolor' 'facecolor', ,'interp' 'interp'); ); view(2),axis([x1 x2 y1 y2]),axis equal equal,colorbar; ,colorbar; %-------------------------------% Calcul de l'erreur en norme L2


error_L2=norm(U-uee,2); err_rela=error_L2\norm(uee,2); fprintf('\t\t fprintf('\t\t %4d \t\t %20.16e \t\t %20.16e\n\n', %20.16e\n\n', npx*npy,error_L2,err_rela); fprintf('On fprintf('On a fini normalement\n'); normalement\n');

Résultats pour P1 :

Résultats pour P2 :


Résultats pour Q1 :

Résultats pour Q2 :


Conclusion Ce projet illustre bien l’importance des méthodes numériques pour la résolution des problèmes mathématique, et il nous a permis de concrétiser nos connaissances dans la résolution des problèmes en deux dimensions selon le type de maillage choisi. D’autre part, il nous a été très utile de travailler sur ce projet, afin de se familiariser davantage avec un outil important pour un élève ingénieur qu’est le Matlab et ainsi reconnaitre son utilité.


Rapport Implementation de La MEF 2D

Recommend Documents