Implémentation Des Éléments Finis en Matlab v5

EDP mixte : Implémentation des éléments finis en Matlab

Implémentation des éléments finis en Matlab 1. Introduction : Une courte exécution de Matlab pour les éléments finis P1-Q1 sur des triangles et des parallélogrammes est donnée pour la résolution numérique des problèmes elliptiques avec des conditions aux frontières mixtes sur des grilles non structurées. Selon la brièveté du programme et de la documentation donnée, n'importe quelle adaptation des exemples modèles simples, à des problèmes plus complexes, peut facilement être exécuté. Les exemples numériques prouvent la flexibilité de l'outil de Matlab. Les programmes proposés de Matlab utilisent la méthode d'élément fini pour calculer une solution numérique U qui rapproche la solution u du problème bidimensionnel de Laplace (P ) avec des conditions aux frontières mélangés (mixtes): 2. Le problème exact : Soit   2 un domaine de Lipschitz avec une frontière polygonale  . Sur un certain sous-ensemble fermé de la frontière D avec longueur positive, nous considérons des conditions de Dirichlet, alors que nous avons les conditions de Neumann sur la partie restante :    \ D Soient f  L () , uD  H 1 () et g  L2 (N ) . u  H ( ) Cherchons avec :  u  f Dans  (1) u  u D Sur D (2) u  g Sur  (3) (P) n N

2

1

N

D’après le théorème de Lax-Milgram, il existe toujours une solution faible de (1)-(3) ce qui donne la régularité intérieure (i.e., 2 u  H loc () ), et on a des frontières lisses et aussi un changement de conditions à la frontière. Les conditions non homogènes de Dirichlet (2) sont associées à la décomposition : donc v  0 sur D , i.e., ()  {   H 1 () /   0 surD }

v  u  uD

vH

1 D

Par : GRARI et KORIKACHE

1/30


Alors, la formulation faible (ou vibrationnelle) du problème (P) est de : Rechercher v  H 1D () tel que :

 v.wdx   f .wdx   g.wds   u






D

2/30



D

.wdx, w  H 1 () D

(4)


3. Discrétisation de Galerkin du problème : Pour l’implémentation, le problème (4) est discrétisé en utilisant la méthode standard de Galerkin ou H () et H D1 () sont remplacées par des sous-espaces de dimensions finis S et S D  S  H D1 , respectivement, Soit U D  S une approximation de uD sur D (On définit U D comme étant un interpolant (relatif au nœud) de uD sur D ). Donc la discrétisation de problème ( P ) est de chercher V  S D tel que ; 1

S

 V .Wdx   fWdx   gWds   U





N

D

.Wdx,

W  SD



(5) (1 ,........, N ) une

Soit

base de l’espace S de dimension fini, et soit (i ,........,i ) une base de S D ou I  {i1 ,........., i M }  {1,........, N} est un ensemble d'index de cardinale M  N  2 . 1

M

 V . dx   f dx   g ds   U j

j



j



N

D



. j dx, j  I

(6)

L’équation (5) est équivalente à : En outre soit, V   xk k k I

Donc, système

N

et U D  U k k . k 1

l’équation (6) donne le linéaire des équations :

Ax  b

(7)

Les coefficients de la matrice A  ( A b  (b )   sont définit par : j

j I

jk

) j , kI 2   M M et

le coté droit

M

A jk 

  j . k dx et b j 





f j dx 



N

 g j ds  U k   j .k dx

N

k 1

(8)



La matrice A est symétrique, définie positive, donc, le système (7) a exactement une solution x   M ce qui détermine la solution de Galerkin N

U  U D V   U j j   xk k j 1


3/30

kI


4. Représentation des données de la triangulation  : Supposant que le domaine  a une frontière polygonale  , nous pouvons recouvrir  par une triangulation régulière T formée de triangles et de quadrilatères, i.e.,    t où chaque t est un triangle fermé ou un quadrilatère fermé. tT

Figure 1. Exemple de maillage

coordinates.dat contient les coordonnées de chaque nœud de maillage donné. Chaque ligne est de la forme : node# x-coord y-coord Notre subdivision  est formée de triangles et de quadrilatères. Dans les deux cas, les nœuds sont numérotés dans le sens contraire des aiguilles d'une montre. Pour les triangles ; elements3.dat contient pour chacun triangle le nombre de nœuds (sommets), Chaque ligne est de la forme : element# node1 node2 node3. De même, les données pour les quadrilatères sont données dans elements4.dat. Ici, nous employons le format : element# node1 node2 node3 node4.


4/30


elements3.dat 1 2 3 13 2 3 4 13 3 4 5 15 4 5 6 15

1 2 3 4 5 6

elements4.dat 1 2 13 12 13 14 13 4 15 11 14 9 14 15 8 15 6 7

12 11 14 10 9 8

Les deux fichiers Neumann.dat et Dirichlet.dat contiennent dans chaque ligne les deux nombres de nœuds attacher au le bord de la frontière : Neumann edge# node1 node2 resp., Dirichlet edge # node1 node2.

Figue2 : fonction chapeau

neumann.dat 1 5 6 2 6 7 3 1 2 4 2 3

dirichlet.dat 1 3 4 2 4 5 3 7 8 4 8 9 5 9 10 6 10 11 7 11 12 8 12 1

Dans la figure2, deux fonctions chapeau typiques  définies pour chaque nœud ( x , y ) du maillage par : j

j

j

 j ( xk , yk )   jk ,


5/30

j , k  1,......., N .


Le sous espace S D  S est l’espace des splines engendré par tout les  pour tout ( x , y ) qui ne se sont pas sur D . D’autre part U D est définit comme étant un interpolant nodal de u D , dans S . Avec ces espaces S et S D et leurs bases correspondantes, les intégrales dans la relation (8) peuvent être calculé comme somme de tous les éléments et aussi somme de tous les bords de l’arc  , c-à-d., j

j

j

N

A jk 

   t T t

j

. k dx

b j    f j dx  t T t

(9) N

 g j ds  U k    jk dx



E  N E

k 1

(10)

tT t

5. la matrice de courbatures La matrice locale de courbatures est déterminée par les coordonnées des sommets de l'élément correspondant, elle est calculé par les fonctions stima3.m et stima4.m. Pour un élément de la triangulation T, soient ( x1 , y1 ) , ( x2 , y2 ) et ( x , y ) des sommets et 1 , 2 et  les fonctions de base correspondantes dans S ,i.e.,  ( x , y )   j , k  1,2,3. 3

3

j

k

k

jk

,

La réflexion d'un moment indique  1 x   j ( x, y )  det  1 x j 1  1 x (11) j2 

y 

 1 xj  y j 1  /det 1 x j 1  1 x y j  2  j2  

D’où


yj   y j 1  y j  2 

1  yj1  yj2  j(x, y)  2T  xj2  xj1  6/30

3


Avec

T

est donné par :  x x x  x  2 T  det 2 1 3 1   y2  y1 y3  y1 

L'entrée résultante de la matrice de courbatures est :

 yk1  yk2 

T   M jk   j (k ) dx  2  yj1  yj2, xj2  xj1   ( 2 T ) t x k2  xk1  t

Avec 3. Ceci pour M 

t 2

.GG

t

avec

 1  G   x1  y  1

1 x2 y2

1  x3  y3 

1

 0   1  0 

l’index modulo est écrit simultanément tous les index :

0  0 1 

Puisque nous obtenons les formules semblables pour trois dimensions, le programme en Matlab est simultanément pour d=2 et d=3 function M = stima3(vertices) d = size(vertices,2); G = [ones(1,d+1);vertices’] \ [zeros(1,d);eye(d)]; M = det([ones(1,d+1);vertices’]) * G * G’ / prod(1:d);

Pour un élément quadrilatéral T soient ( x1 , y1 ),........, ( x4 , y4 ) les sommets avec les fonctions chapeau correspondantes 1 ,......., 4 .puisque T est un parallélogramme, il y a un quadrillage


7/30


 x  2xx 1 4xx 1  x1   T(,)     y   2yy 1 4yy 1  y1 Pour les éléments [0,1] sur T. puis  fonctions de la forme 2

j

,

( x, y )   j ( T1 ( x, y )) avec

1 ( ,  ) : (1   )(1   ),

 2 ( ,  ) :  (1   ),

3 ( ,  ) :  ,

 4 ( ,  ) : (1   ) .

les

De la loi de substitution il suit pour les intégrales (9) : M

jk

M jk

   

T

 j ( x, y ). k ( x, y ) d ( x, y )

0 ,1 2

 (   T1 )( ( ,  ))(( k   T1 ))( ( ,  ))T det D T d ( ,  )

M jk  det  D T  

 0 ,1 2

 j ( ,  )(( D T ) T D T ) 1 ( k ( ,  ))T d ( ,  )

On résout ces intégrales à partir de la matrice locale de courbatures, pour un élément de quadrilatère on aura :


8/30


 2a  c  3b  (a  c) a  2c   3b  2(a  c)    2a  c  3b  2(a  c) a  2c 3b  2(a  c)  det  D T   M    3b  (a  c) a  2c 3b  2(a  c)  2a  c  6    a  2c 3b  (a  c)  2a  c 3b  2(a  c)  

Avec

 a b ((DT )T DT )1    b c

function M = stima4(vertices) D_Phi = [vertices(2,:)vertices(1,:); vertices(4,:) ... vertices(1,:)]’; B = inv(D_Phi’*D_Phi); C1 = [2,2;2,2]*B(1,1)+[3,0;0,3]*B(1,2)+[2,1;1,2]*B(2,2); C2 = [1,1;1,1]*B(1,1)+[3,0;0,3]*B(1,2)+[1,2;2,1]*B(2,2); M = det(D_Phi) * [C1 C2; C2 C1] / 6;

6. Assembler le côté droit de l’équation Les forces de volume sont employées pour assembler le côté droit. Utilisons la valeur de f au centre de gravité ( x , y ) de T l'intégrale T f  j dx en (10) est approximée par : S

S

1  x2  x1 x3  x1  T f jdx  kT det y2  y1 y3  y1  f (xS , yS ) Tel que kT  6 si T est un triangle et parallélogramme.

kT  4

si

T est un

% Volume Forces for j = 1:size(elements3,1) b(elements3(j,:)) = b(elements3(j,:)) + ... det([1 1 1; coordinates(elements3(j,:),:)’]) * ... f(sum(coordinates(elements3(j,:),:))/3)/6; end for j = 1:size(elements4,1) b(elements4(j,:)) = b(elements4(j,:)) + ... det([1 1 1; coordinates(elements4(j,1:3),:)’]) * ... f(sum(coordinates(elements4(j,:),:))/4)/4; end

Les valeurs de du problème.

f

sont données à partir de la fonction f.m qui dépend


9/30


La fonction est définit par les coordonnées des points qui se trouve dans Ω et elle renvoie la force de volume a ces endroits. Pour l’exemple numérique représenté sur le schéma 3 nous avons employé function VolumeForce = f(x); VolumeForce = ones(size(x,1),1);

De même, les conditions de Neumann contribuent au coté droit. En utilisant la valeur de g au centre ( xM , yM ) de E avec la longueur E , l’intégrale E g j ds dans (10) est approché par  g

j

ds 

E

E g ( x M , y M ). 2

% Neumann conditions for j = 1 : size(neumann,1) b(neumann(j,:))=b(neumann(j,:)) + ... norm(coordinates(neumann(j,1),:) ... coordinates(neumann(j,2),:)) * ... g(sum(coordinates(neumann(j,:),:))/2)/2; end

Ici, nous employons le fait que dans Matlab le taille d’une matrice vide est placé par zéro et qu'une boucle de 1 à 0 est totalement omis. De cette façon, la question de l'existence des données de frontière de Neumann doit être renoncée. Les valeurs de g sont donnés par la fonction g.m qui dépend encore du problème. La fonction est définit avec les coordonnées des points sur  et retours les efforts correspondants. Pour l'exemple numérique g.m était N

function Stress = g(x) Stress = zeros(size(x,1),1);

7. États d’incorporation de Dirichlet Avec une numérotation appropriée des nœuds, le système des équations linéaires résultant de la construction décrite dans la section précédente sans incorporer des états de Dirichlet peut être écrit comme suit :


 A1 A12 U b 

10/30


Avec U   M , U D   N M . Ici, U sont les valeurs aux nœuds libres qui sont à déterminer, U D sont les valeurs aux nœuds qui sont sur la frontière de Dirichlet ainsi sont connus a priori. Par conséquent, le premier bloc d'équations peut être récrit : A11 .U  b  A12 .U D

C'est la formulation de (6) avec U D  0 aux nœuds de non-Dirichlet. Dans le deuxième bloc d'équations dans (12) l'inconnu est bD mais puisqu'il n'a % Dirichlet conditions pas u = sparse(size(coordinates,1),1); u(unique(dirichlet)) = u_d(coordinates(unique(dirichlet),:)); d’intérêt, il b = b A * u; est omis dans le suivant.

Les valeurs u D aux nœuds sur D sont données par la fonction u.d.m qui dépend du problème. La fonction est appelée par les coordonnées aux points sur D et retourne les valeurs aux endroits correspondants. Pour l'exemple numérique u.d.m était : function DirichletBoundaryValue = u_d(x) DirichletBoundaryValue = zeros(size(x,1),1);

8. Calcul de la solution numérique Les lignes de (7) correspondant à la première M ligne de la forme (12) qui réduit le système des équations avec une matrice définie symétrique et positive de coefficient A11 . Il est obtenu du système original des équations en prenant les lignes et les colonnes et on les fait correspondant les nœuds libres du problème. La restriction peut être réalisée dans Matlab à travers l’indexation appropriée.


11/30


Le système des équations est résolu par l'opérateur binaire (installé dans Matlab) FreeNodes=setdiff(1:size(coordinates,1),unique(dirichlet)); qui donne u(FreeNodes)=A(FreeNodes,FreeNodes)\b(FreeNodes); l'inverse gauche d'une matrice.

Matlab se sert des propriétés d'une matrice définie positive, symétrique pour résoudre le système des équations efficacement. Une représentation graphique de la solution est donnée par la fonction show.m. function show(elements3,elements4,coordinates,u) trisurf(elements3,coordinates(:,1),coordinates(:,2),u’,... ’facecolor’,’interp’) hold on trisurf(elements4,coordinates(:,1),coordinates(:,2),u’,... ’facecolor’,’interp’) hold off view(10,40); title(’Solution of the Problem’)

Figure 3.solution du problème de Laplace


12/30


Dans Matlab trisurf(ELEMENTS,X,Y,U) est utilisée pour dessiner des triangulations pour les éléments de types égaux. Chaque ligne de la matrice ELEMENTS détermine un polygone où les x-, y-, et z- coordonnées de chaque sommet de ce polygone est donnée par la correspondance avec X, Y et U respectivement. La couleur des polygones est donnée par des valeurs de U. Les paramètres additionnels, 'facecolor', 'interp', mènent à une coloration interpolée. La figure 3 illustre la solution pour le maillage définie dans la section 4 et les fichiers de données f.m, g.m, et u_d.m donnés dans les sections 6 et 7 . La récapitulation sectionne 4-8, le programme principal, qui est       

Lignes 3-10: Chargement de la géométrie et initialisation du maillage. Lignes 11-19: Assemblée la matrice de rigidité dans deux boucles, d'abord aux éléments des triangulaires, puis aux éléments des quadrilatères. Lignes 20-30: Incorporation de la force de volume dans deux boucles, d'abord aux éléments des triangulaires, puis aux éléments des quadrilatères. Lignes 31-35: Incorporation de l'état de Neumann. Lignes 36-39: Incorporation de l'état de Dirichlet. Lignes 40-41: Solution du système linéaire réduit. Lignes 42-43: Représentation graphique de la solution numérique.

énuméré dans l'annexe A, est structuré comme suit (les lignes références sont selon la numérotation dans l'annexe A) : 9. L'équation de la chaleur Pour des méthodes numériques de l'équation de la chaleur, u

t

 u  f sur

   0, T 

avec un procédé implicite à temps d'Euler, nous avons devisé l'intervalle de temps [0, T] en N sous-intervalles de taille dt  T N qui mène à l'équation : (id  dt )u n  dtf n  u n1 ,

où

f n  f ( x, t n ) et u n

est l'approximation discrète de u au temps


13/30

(13)

tn  n


La forme faible de (13) est : u



n

vdx  dt  u n .vdx  dt (  f n vdx  



g



n

vdx ) 

u

n 1

vdx



Avec g  g ( x, t ) et les notations dans la section 2. Pour chaque étape, cette équation est résolue en utilisant les éléments finis qui mène au système linéaire suivant : n

n

( dtA  B )U n  dtb  BU n 1

La matrice de courbatures A et le côté droit b sont comme avant. La matrice de masse B (voir 8) est le résultat des limites  u n vdx , i.e., B jk 

  

T  T

j

k

dx.

Pour la triangulation, affinez par morceaux les éléments que nous obtenons :

 2 1 1  1  x2  x1 x3  x1       dx  det 1 2 1 T j k 24  y2  y1 y3  y1     1 1 2   L'annexe B montre le programme modifié de l'équation de la chaleur. L’exemple numérique a été basé sur le domaine dans la figure1, cette fois avec f  0 et u D  1 sur la frontière externe. La valeur sur le cercle (intérieur) est toujours u D  0 . Sur les frontière de Neumann, nous avons toujours f  f ( x, t ) , La figure 4 montre la solution pendant quatre fois différentes (T = 0:1, 0:2, 0:5 et T = 1). (T est la variable dans la ligne 10 du programme principal.) Le programme principal, listé dans l'annexe B, est structuré comme suit (les lignes référence ont la même numérotation dans l'annexe B) : n

n

 Lignes 3-11 : Chargement de la géométrie et initialisation du maillage  Lignes 12-16 : Assemblée la matrice de courbatures A dans une boucle en tous les éléments triangulaires.  Lignes 17-20 : Assemblée la matrice de masse B dans une boucle en tous les éléments triangulaires. Lignes 21-22 : Définir l'état initial du discret U.  Lignes 23-48 : Boucle (aux au-dessus) étapes de temps.  En particulier :  Ligne 25 : (Dégager) le vecteur du côté droit.  Lignes 26-31 : Incorporation de la force de volume à l'étape de temps n.  Lignes 32-37 : Incorporation de la condition de Neumann à l'étape de temps n.  Lignes 38-39 : Incorporation de la solution à l'étape précédente de temps n -1.  Lignes 40-43 : Incorporation de la condition de Dirichlet à l'étape de temps n.  Lignes 44-47 : Solution du système Par : GRARI et KORIKACHE 14/30linéaire réduit pour la solution à l'étape de temps n.  Lignes 49-50 : Représentation graphique de la solution numérique à l’étape temps final


L'annexe B montre le programme modifié de l'équation de la chaleur. L’exemple numérique a été basé sur le domaine dans la figure1, cette fois avec f  0 et u D  1 sur la frontière externe. La valeur sur le cercle (intérieur) est toujours u D  0 . Sur la frontière de Neumann, nous avons toujours g  0 .


15/30


10. Un problème non-linéaire Comme application simple du problème variationnel non convexe, nous considérons l'équation de Ginzburg-Landau u  u 3  u dans  , u  0 sur 

Pour 1 100

J (u , v ) :

 u.vdx   (u  u

, (15) c’est la formulation faible, i.e.,



3

(14)

)vdx  0



v  H 01 ()

On peut considérer également la condition nécessaire pour minimiser le problème variationnel 2 1 2   u   u 2  1  dx! 2 4  

min  

(16)

Nous visons à résoudre (15) avec la méthode de Newton-Raphson's. Commençant par un certain u 0 , dans chacun étape d'itération, nous calculons u n  u n1  H 01 () satisfaisant : DJ (u n , v; u n  u n 1 )  J (u n , v ),

Ou

DJ (u , v;  ) 

v  H 01 ()

 v.dx   (v  3vu



2

(17)

 ) dx.



(18)

Les intégrales dans J (U , V ) et DJ (U , V ; W ) sont de nouveau calculés comme somme de tous les éléments. Les intégrales locales résultantes peuvent être calculées analytiquement et sont implémentés en localj.m, respectivement localdj.m, comme donné dans l'annexe C. Le programme en Matlab a besoin encore de petites modifications, montrées dans l'annexe C. Essentiellement, on doit initialiser le programme (avec un vecteur de début qui remplit la condition de Par : GRARI et KORIKACHE

16/30


frontière de Dirichlet (lignes 9 et 10)), pour ajouter une boucle (lignes 12 et 45), pour mettre à jour la nouvelle approximation de newton (ligne 41), et pour fournir un critère d’arrêt en cas de convergence (lignes 42-44). On sait que les solutions ne sont pas uniques. En effet, pour tout minimum local u, -u est également un minimum et 0 résout aussi le problème. La fonction constante u  1 mène à l'énergie nulle, mais viole la continuité ou on a les conditions aux frontières. Par conséquent, on observe la frontière ou les couches internes qui séparent de grandes régions, où u est presque constant  1 . Dans le problème en dimension finie, les différentes valeurs initiales u 0 peut mener à différentes approximations numériques.

Figure 5. Solution de l’équation non-linéaire

La figure 5 montre deux solutions possibles trouvées pour deux différentes valeurs après environ 20-30 itérations. La figure du côté gauche est réalisée en des valeurs comme étant choisies dans le programme dans l'annexe C. Changer le rapport dans la ligne 9 dans l'annexe C à U = signe (coordonnées (:, 1)); montrer à la figure du côté droit.


17/30


Le programme principal, donné dans l'annexe C, est structuré comme suit (les lignes références sont selon la numérotation dans l'annexe  Lignes 3-7 : Chargement de la géométrie et initialisation du maillage.  Lignes 8-10 : Réglage du vecteur d’initialisation U pour le procédé de l'itération, incorporant la condition de Dirichlet sur la solution.  Lignes 11-45 : Boucle pour l'itération de Newton-Raphson. Il finit après un maximum de 50 itérations (dans la ligne 12) ou en cas de convergence (lignes 4244).  Lignes 13-18 : Assemblage de la matrice de la dérivé du fonctionnel J évalué à l'étape courante d'itération U.  Lignes 19-24 : Assemblage du vecteur du fonctionnel J évalué à la courante étape d'itération U.  Lignes 25-30 : Incorporation de la force de volume.  Lignes 31-35 : Incorporation de l'état de Neumann.  Lignes 36-38 : Incorporation des conditions homogènes de Dirichlet du vecteur de mise à jour W.  Lignes 39-40 : Solution du système linéaire réduit pour le vecteur de mise à jour W.  Ligne 41 : Mise à jour U.  Lignes 42-44 : Éclatement de la boucle si le vecteur de mise à jour W est suffisamment petit (sa norme étant plus petite que 10 10 ).  Lignes 46-47 : La représentation graphique de la finale itération.

C) :


18/30


11. Problèmes tridimensionnels : Avec quelques modifications, le programme de Matlab pour des problèmes linéaires en deux dimensions étudié dans les sections 5-8 peut être prolongé aux problèmes à trois dimensions. Tétraèdres sont utilisés en tant qu'éléments finis. Les fonctions de base sont correspondantes à celles définie en deux dimensions, par exemple, pour un élément de tétraèdre T soient ( x , y , z )( j  1,......,4) les sommets et  les fonctions de base correspondantes, c.-à-d., j

j

j

j

 j ( xk , y k , z k )   jk , j , k  1,.........,4.

Chacun des dossiers *.dat obtient une entrée additionnelle par ligne. Dans coordinates.dat, c’est le zéme-composant de chaque nœud P  ( x , y , z ) Une entrée typique dans elements3.dat se relit maintenant : j k l m n, j

j

j

j

Tel que k, l, m, n, sont les nombres de sommets P ,......., P du jéme élément. elements4.dat n'est pas utilisé pour des problèmes à trois dimensions. k

n

L'ordre des nœuds est organisé tels que le côté droit de  1   x 6 T  det  k y  k  z  k

1 xl yl

1 xm ym

zl

zm

1  xn  yn  

z n 

est positif, La numérotation des éléments définis dans neumann.dat et dirichlet.dat est fait avec le visionnement positif mathématique d'orientation de l'extérieur W sur la surface.


19/30


En utilisant le code de Matlab dans l'annexe A, l'annulation des lignes 5, 16-19 et 26-30 et substitution de 22-24, 33-34, 43 par les lignes suivantes b(elements3(j,:)) = b(elements3(j,:)) + ... donne un det([1,1,1,1;coordinates(elements3(j,:),:)’]) * ... outil f(sum(coordinates(elements3(j,:),:))/4) / 24; b(neumann(j,:)) = b(neumann(j,:)) + ... court et norm(cross(coordinates(neumann(j,3),:) ... flexible coordinates(neumann(j,1),:),coordinates(neumann(j,2),:) ... coordinates(neumann(j,1),:))) ... pour * g(sum(coordinates(neumann(j,:),:))/3)/6; showsurface([dirichlet;neumann],coordinates,full(u)); résoudre la grandeur scalaire, problèmes à trois dimensions linéaires : La représentation graphique pour des problèmes à trois dimensions peut être faite par raccourcis version de show.m de la section 8. function showsurface(surface,coordinates,u) trisurf(surface,coordinates(:,1),coordinates(:,2),... coordinates(:,3),u’, ’facecolor’,’interp’) axis off view(160,30)

La distribution de la température d'un piston simplifié est présentée sur le schéma 6. Calcul de la distribution de la température avec 3728 nœuds et 15111 éléments (y compris le rendement graphique) prend quelques minutes sur un poste de travail.

Figure 6. La distribution de la température d'un piston


20/30


Le programme principal, qui est énuméré dans l'annexe D, est       

Lignes 3-9 : Chargement de la géométrie et de l'initialisation de maille. Lignes 11-14 : Assemblée de la matrice de rigidité dans une boucle au-dessus de tous les tétraèdres. Lignes 16-20 : Incorporation de la force de volume dans une boucle au-dessus de tous les tétraèdres. Lignes 22-27 : Incorporation de l'état de Neumann. Lignes 29-31 : Incorporation de l'état de Dirichlet. Ligne 33 : Solution du système linéaire réduit. Ligne 35 : Représentation graphique de la solution numérique.

structuré comme suit (les lignes références sont selon la numérotation dans l'annexe D) :


21/30


Annexe A. Le code complet de Matlab pour le problème à deux dimensions de Laplace Le programme suivant peut être trouvé dans le paquet, sous le chemin acf/fem2d. Il s'appelle fem2d.m. Les autres dossiers sous ce chemin sont les fonctions fixes stima3.m, stima4.m, et show.m aussi bien que les fichiers de fonctions et de données cela décrivez la discrétisation et les données du problème, à savoir coordinates.dat, elements3.dat, elements4.dat, dirichlet.dat, neumann.dat, f.m, g.m, et u_d.m. Ces problèmes-qui décrivent des dossiers doivent être adaptés par l'utilisateur pour d'autres géométries, discrétisations, et/ou données.


22/30


1 % de FEM2D de méthode d'élément fini bidimensionnelle pour Laplacian. 2 % d'initialisation 3 charge coordinates.dat ; coordonnées (: , 1) = [] ; 4 eval (' charge elements3.dat ; elements3 (: , 1) = [] ; ',' elements3= [] ; ') ; 5 eval (' charge elements4.dat ; elements4 (: , 1) = [] ; ',' elements4= [] ; ') ; 6 eval (' charge neumann.dat ; neumann (: , 1) = [] ; ',' neumann= [] ; ') ; 7 charge dirichlet.dat ; dirichlet (: , 1) = [] ; 8 FreeNodes=setdiff (1 : taille (coordonnées, 1), unique (dirichlet)); 9 A = clairsemé (taille (coordonnées, 1), taille (coordonnées, 1)); 10 b = clairsemé (taille (coordonnées, 1), 1) ; 11 % d'Assemblée 12 pour j = 1 : taille (elements3,1) 13 A (elements3 (j :), elements3 (j :)) = A (elements3 (j :),… 14 elements3 (j :)) + stima3 (coordonnées (elements3 (j :) :)); extrémité 15 16 pour j = 1 : taille (elements4,1) 17 A (elements4 (j :), elements4 (j :)) = A (elements4 (j :),… 18 elements4 (j :)) + stima4 (coordonnées (elements4 (j :) :)); extrémité 19 20 % de forces de volume 21 pour j = 1 : taille (elements3,1) 22 b (elements3 (j :)) = b (elements3 (j :)) +… det 23 ([1.1.1 ; coordonnées (elements3 (j :) :)']) *… 24 f (somme (coordonnées (elements3 (j :) :))/3) /6 ; extrémité 25 26 pour j = 1 : taille (elements4,1) 27 b (elements4 (j :)) = b (elements4 (j :)) +… det 28 ([1.1.1 ; coordonnées (elements4 (j, 1 : 3) :)']) *… 29 f (somme (coordonnées (elements4 (j :) :))/4) /4 ; extrémité 30 31 % d'états de Neumann 32 pour j = 1 : taille (neumann, 1) 33 b (neumann (j :))=b (neumann (j :)) +… norme 34 (coordonnées (neumann (j, 1) :) coordonnées (neumann (j, 2) :))*… g (somme (coordonnées (neumann (j :) :))/2) /2 ; extrémité 35 36 % d'états de Dirichlet 37 u = clairsemé (taille (coordonnées, 1), 1) ; 38 u (uniques (dirichlet)) = u_d (coordonnées (uniques (dirichlet) :)); 39 b = b A * u ; 40 % de calcul de la solution 41 u (FreeNodes) = A) (de FreeNodes, de FreeNodes \ b (FreeNodes) ; 42 % de représentation de graphique exposition 43 (elements3, elements4, coordonnées, pleines (u));


23/30


B. Code de Matlab pour l'équation de la chaleur Le programme suivant peut être trouvé dans le paquet, sous le chemin acf/fem2d la chaleur. Il s'appelle fem2d heat.m. Les autres dossiers sous ce chemin sont les fixes fonctions stima3.m et show.m aussi bien que les fichiers de fonctions et de données qui décrivent la discrétisation et les données du problème, à savoir coordinates.dat, elements3.dat, dirichlet.dat,neumann.dat,f.m,g.m,et u_d.m.Ces problèmes-qui décrivent des dossiers doivent être adaptés par l'utilisateur pour d'autres géométries, discrétisations, et/ou données.


24/30


1 méthode d'élément fini de %FEM2D_HEAT pour l'équation bidimensionnelle de la chaleur. 2 %Initialisation 3 charge coordinates.dat ; coordonnées (: , 1) = [] ; 4 charge elements3.dat ; elements3 (: , 1) = [] ; 5 eval (' charge neumann.dat ; neumann (: , 1) = [] ; ',' neumann= [] ; ') ; 6 charge dirichlet.dat ; dirichlet (: , 1) = [] ; 7 FreeNodes=setdiff (1 : taille (coordonnées, 1), unique (dirichlet)); 8 A = clairsemé (taille (coordonnées, 1), taille (coordonnées, 1)); 9 B = clairsemé (taille (coordonnées, 1), taille (coordonnées, 1)); 10 T = 1 ; décollement = 0.01 ; N = T/dt ; 11 U = zéros (taille (coordonnées, 1), N+1) ; 12 % d'Assemblée 13 pour j = 1 : taille (elements3,1) 14 A (elements3 (j :), elements3 (j :)) = A (elements3 (j :),… 15 elements3 (j :)) + stima3 (coordonnées (elements3 (j :) :)); extrémité 16 17 pour j = 1 : taille (elements3,1) 18 B (elements3 (j :), elements3 (j :)) = B (elements3 (j :),… 19 elements3 (j :)) + det ([1.1.1 ; coordonnées (elements3 (j :) :)'])… * [2.1.1 ; 1.2.1 ; 1.1.2] /24 ; extrémité 20 21 % d'état initial 22 U (: , 1) = zéros (taille (coordonnées, 1), 1) ; 23 % d'étapes de temps 24 pour n = 2 : N+1 25 b = clairsemé (taille (coordonnées, 1), 1) ; 26 % de forces de volume 27 pour j = 1 : taille (elements3,1) 28 b (elements3 (j :)) = b (elements3 (j :)) +… det 29 ([1.1.1 ; coordonnées (elements3 (j :) :)']) *… dt*f 30 (somme (coordonnées (elements3 (j :) :))/3, n*dt) /6 ; extrémité 31 32 % d'états de Neumann 33 pour j = 1 : taille (neumann, 1) 34 b (neumann (j :)) = b (neumann (j :)) +… norme 35 (coordonnées (neumann (j, 1) :) coordonnées (neumann (j, 2) :))*… dt*g 36 (somme (coordonnées (neumann (j :) :))/2, n*dt) /2 ; extrémité 37 38 % de timestep précédent 39 b = b + B * U (: , n1) ; 40 % d'états de Dirichlet 41 u = clairsemé (taille (coordonnées, 1), 1) ; 42 u (uniques (dirichlet)) = u_d (coordonnées (uniques (dirichlet) :), n*dt) ; 43 b = b (décollement * A + B) * u ; 44 % de calcul de la solution 45 u (FreeNodes) = (dt*A (FreeNodes, FreeNodes) +… 46 B (FreeNodes, FreeNodes))\ b (FreeNodes) ; 47 U (: , n) = u ; extrémité 48 49 % de représentation de graphique exposition 50 (elements3, [], coordonnées, complètement (U (: , N+1)));


25/30


C. Code de Matlab pour le problème non-linéaire Le programme suivant peut être trouvé dans le paquet, sous le chemin acf/fem2d non-linéaire. Il s'appelle fem2d nonlinear.m. Les autres dossiers sous ce chemin sont la fonction fixe show.m aussi bien que les fichiers de fonctions et de données qui décrivent le fonctionnel J, son dérivé DJ, la discrétisation et les données du problème, à savoir, localj.m, localdj.m, coordinates.dat, elements3.dat,dirichlet.dat, f.m, et u_d.m. Ces problèmes-qui décrivent des dossiers doivent être adaptés par l'utilisateur pour d'autres problèmes, géométries, discrétisations, et/ou données non-linéaires (y compris ajouter probablement les dossiers appropriés neumann.dat et g.m).


26/30


1 % de FEM2D_NONLINEAR de méthode d'élément fini pour bidimensionnel % d'équation nonlinéaire. 2 % d'initialisation 3 charge coordinates.dat ; coordonnées (: , 1) = [] ; 4 charge elements3.dat ; elements3 (: , 1) = [] ; 5 eval (' charge neumann.dat ; neumann (: , 1) = [] ; ',' neumann= [] ; ') ; 6 charge dirichlet.dat ; dirichlet (: , 1) = [] ; 7 FreeNodes=setdiff (1 : taille (coordonnées, 1), unique (dirichlet)); 8 % de valeur initiale 9 U = ceux (taille (coordonnées, 1), 1) ; 10 U (uniques (dirichlet)) = u_d (coordonnées (uniques (dirichlet) :)); 11 % d'itération de NewtonRaphson 12 pour i=1 : 50 13 % d'Assemblée de DJ (U) 14 A = clairsemé (taille (coordonnées, 1), taille (coordonnées, 1)); 15 pour j = 1 : taille (elements3,1) 16 A (elements3 (j :), elements3 (j :)) = A (elements3 (j :),… elements3 (j :)) … 17 + localdj (coordonnées (elements3 (j :) :), U (elements3 (j :))); extrémité 18 19 % d'Assemblée de J (U) 20 b = clairsemé (taille (coordonnées, 1), 1) ; 21 pour j = 1 : taille (elements3,1) ; 22 b (elements3 (j :)) = b (elements3 (j :)) … 23 + localj (coordonnées (elements3 (j :) :), U (elements3 (j :))); 24 extrémités 25 % de forces de volume 26 pour j = 1 : taille (elements3,1) 27 b (elements3 (j :)) = b (elements3 (j :)) +… det 28 ([1 1 1 ; coordonnées (elements3 (j :) :)']) *… 29 f (somme (coordonnées (elements3 (j :) :))/3) /6 ; extrémité 30 31 % d'états de Neumann 32 pour j = 1 : taille (neumann, 1) 33 b (neumann (j :))=b (neumann (j :)) … norme (coordonnées (neumann (j, 1) :)… 34 coordonnées (neumann (j, 2) :)) *… *g (somme (coordonnées (neumann (j :) :))/2) /2 ; extrémité 35 36 % d'états de Dirichlet 37 W = zéros (taille (coordonnées, 1), 1) ; 38 W (uniques (dirichlet)) = 0 ; 39 % résolvant une étape de newton 40 W (FreeNodes) = A) (de FreeNodes, de FreeNodes \ b (FreeNodes) ; 41 U = U W ; 42 si norme (W) < 10ˆ(10) coupure 43 extrémité 44 extrémité 45 46 % de représentation de graphique exposition 47 (elements3, [], coordonnées, pleines (U)); fonction b = localj (sommets, U) Eps = 1/100 ;


27/30

EDP mixte : Implémentation des éléments finis en Matlab G = [ceux (1.3) ; ] de sommets' \ [zéros (1.2) ; oeil (2)]; Secteur = det ([ceux (1.3) ; ) de sommets']/2 ; b=Area* ((Eps*G*G' [2.1.1 ; 1.2.1 ; 1.1.2] /12) *U+ … [4*U (1) ˆ3+ U (2) ˆ3+U (3) ˆ3+3*U (1) ˆ2* (U (2) +U (3))+2*U (1)… * (U (2) ˆ2+U (3) ˆ2) *U de +U (2) (3) * (U (2) +U (3))1) *U du *U de +2*U ((2) (3) ; 4*U (2) ˆ3+ U (1) ˆ3+U (3) ˆ3+3*U (2) ˆ2* (U (1) +U (3))+2*U (2)… * (U (1) ˆ2+U (3) ˆ2) 1) *U de +U ((3) * (U (1) +U (3))1) *U du *U de +2*U ((2) (3) ; 4*U (3) ˆ3+ U (2) ˆ3+U (1) ˆ3+3*U (3) ˆ2* (U (2) +U (1))+2*U (3)… * (U (2) ˆ2+U (1) ˆ2) *U de +U (2) (1) * (U (2) +U (1))1) *U du *U de +2*U ((2) (3)]/60) ; fonction M = localdj (sommets, U) Eps = 1/100 ; G = [ceux (1.3) ; ] de sommets' \ [zéros (1.2) ; oeil (2)]; Secteur = det ([ceux (1.3) ; ) de sommets']/2 ; M = Area* (Eps*G*G' [2.1.1 ; 1.2.1 ; 1.1.2] /12 +… [12*U (1) ˆ2+2* *U (d'U (2) ˆ2+U (3) ˆ2+U (2) (3))+6*U (1) * (U (2) +U (3)),… *U de 3* (U (1) ˆ2+U (2) ˆ2) +U (3) ˆ2+4*U (1) (2) +2*U (3) * (U (1) +U (2)),… *U de 3* (U (1) ˆ2+U (3) ˆ2) +U (2) ˆ2+4*U (1) (3) +2*U (2) * (U (1) +U (3)); *U de 3* (U (1) ˆ2+U (2) ˆ2) +U (3) ˆ2+4*U (1) (2) +2*U (3) * (U (1) +U (2)),… 12*U (2) ˆ2+2* *U (d'U (1) ˆ2+U (3) ˆ2+U (1) (3))+6*U (2) * (U (1) +U (3)),… *U de 3* (U (2) ˆ2+U (3) ˆ2) +U (1) ˆ2+4*U (2) (3) +2*U (1) * (U (2) +U (3)); *U de 3* (U (1) ˆ2+U (3) ˆ2) +U (2) ˆ2+4*U (1) (3) +2*U (2) * (U (1) +U (3)),… *U de 3* (U (1) ˆ2+U (3) ˆ2) +U (2) ˆ2+4*U (1) (3) +2*U (2) * (U (1) +U (3)),… *U de 3* (U (2) ˆ2+U (3) ˆ2) +U (1) ˆ2+4*U (2) (3) +2*U (1) * (U (2) +U (3)),… 12*U (3) ˆ2+2* *U (d'U (1) ˆ2+U (2) ˆ2+U (1) (2))+6*U (3) * (U (1) +U (2))]/60) ;


28/30


D. Code de Matlab pour le problème à trois dimensions Le programme suivant peut être trouvé dans le paquet, sous le chemin acf/fem3d. Il s'appelle fem3d.m. Les autres dossiers sous ce chemin sont les fonctions fixes stima3.m, et showsurface.m aussi bien que les fichiers de fonctions et de données qui décrivent la discrétisation et les données du problème, à savoir coordinates.dat, elements3.dat, dirichlet.dat, neumann.dat, f.m, g.m, et u_d.m. Ces problèmes-qui décrivent des dossiers doivent être adaptés par l'utilisateur pour d'autres géométries, discrétisations, et/ou données. 1 % de FEM3D de méthode d'élément fini tridimensionnelle pour Laplacian. 2 % d'initialisation 3 charge coordinates.dat ; coordonnées (: , 1) = [] ; 4 charge elements3.dat ; elements3 (: , 1) = [] ; 5 eval (' charge neumann.dat ; neumann (: , 1) = [] ; ',' neumann= [] ; ') ; 6 charge dirichlet.dat ; dirichlet (: , 1) = [] ; 7 FreeNodes=setdiff (1 : taille (coordonnées, 1), unique (dirichlet)); 8 A = clairsemé (taille (coordonnées, 1), taille (coordonnées, 1)); 9 b = clairsemé (taille (coordonnées, 1), 1) ; 10 % d'Assemblée 11 pour j = 1 : taille (elements3,1) 12 A (elements3 (j :), elements3 (j :)) = A (elements3 (j :),… 13 elements3 (j :)) + stima3 (coordonnées (elements3 (j :) :)); extrémité 14 15 % de forces de volume 16 pour j = 1 : taille (elements3,1) 17 b (elements3 (j :)) = b (elements3 (j :)) +… det 18 ([1.1.1.1 ; coordonne (elements3 (j :) :)'])… 19 * f (somme (coordonnées (elements3 (j :) :))/4)/24 ; extrémité 20 21 % d'états de Neumann 22 pour j = 1 : taille (neumann, 1) 23 b (neumann (j :)) = b (neumann (j :)) +… 24 normes (croix (coordonnées (neumann (j, 3) :)… coordonnées (neumann (j, 1) :), … 25 coordonnées (neumann (j, 2) :) coordonnées (neumann (j, 1) :))) … 26 * g (somme (coordonnées (neumann (j :) :))/3) /6 ; extrémité 27 28 % d'états de Dirichlet 29 u = clairsemé (taille (coordonnées, 1), 1) ; 30 u (uniques (dirichlet)) = u_d (coordonnées (uniques (dirichlet) :)); 31 b = b A * u ; 32 % de calcul de la solution 33 u (FreeNodes) = A) (de FreeNodes, de FreeNodes \ b (FreeNodes) ; 34 % de représentation de graphique showsurface 35 ([dirichlet ; neumann], coordonnées, pleines (u));


29/30


Références : [1] S.C. Brenner et L.R. Scott, La théorie mathématique de méthodes d'élément fini, Textes dans Mathématiques appliquées, vol. 15 (Springer, New York, 1994). [2] P.G. Ciarlet, La méthode d'élément fini pour des problèmes elliptiques (la Nord-Hollande, Amsterdam, 1978). [3] L. Langemyr et autres., Le guide d'équation de l'utilisateur partiel de boîte à outils ( Math Works, Inc. 1995). [4] H.R. Schwarz, Der Finiten Elemente de Methode (Teubner, Stuttgart, 1991).


30/30

Implémentation Des Éléments Finis en Matlab v5

Recommend Documents