Considérations d'ordre numérique pour le problème d'hygro-mécanique

dans le cas de plaques de composite.

Exploration

La première série de tests porte sur M

la détermination du type de schéma d'intégration suffisant pour les calculs exigés.
L'"optimisation" du type d'interpolation nécessaire pour M et pour U.
L'effet du maillage sur la solution.

Les questions auxquels je tente de répondre sont les suivantes:

Quel est l'effet de l'interpolation sur M (passage linéaire à quadratique à Crouzeix-Raviart)?
Quel est l'effet d'un changement d'interpolation de M sur le résultat en U?
Sachant que l'interpolation en U doit être quadratique pour éviter le vérouillage quel est le "plus petit" quadratique que je peux employé?
On sait que le schéma d'intégration à un effet sur les résultats. Quel doit être le schéma d'intégration pour chacune des interpolations et pour les deux types d'éléments possibles (héxaèdre et prisme)?
Comment se comporte la solution relativement au maillage, essais d'adaptation et l'utilité de l'adapatation pour ce type de problème?

Le choix de l'interpolation à un effet direct sur la taille du système matricielle à résoudre et le schéma d'intégration à un effet direct sur le temps CPU requis pour l'assemblage du système. Les deux raisons majeures pour s'attarder à ces questions sont: (i) le problème en M est instationnaire et non-linéaire, on va le résoudrede manière itérative ce qui veut dire que l'on fera des réassemblage plus ou moins fréquent; (ii) on prévoit travailler sur des precessus d'optimisation et de sensibilité du système i.e. on a encore affaire à des processus itératifs où on résoud de nombreuse fois le système matricielle. Dans ces deux cas la taille du système global et le temps d'assemblage de celui-ci devient important.

Le cas test.

Comme cas tests j'ai choisi le problème de MDF avec coefficients constants. Il s'agit d'un cas où l'on retrouve des symétries ce qui permet de réduire la taille du problème à l'étude. On à un plaque de 560 mm x 560 mm et d'épaisseur 12 mm, composée de 3 couches: 2 mm pour la couche du dessus, 8 mm pour la couche interne et 2 mm pour la couche du dessous. Les données:

M(x,y,z,0) = 7
$q = - K \nabla M \cdot n = 3.2e-4 (M - 10)$ sur la surface du dessus
U_X = 0 sur (280,y,0)
U_Y = 0 sur (x,280,0)
U_Z = 0 sur (280,280,0), (280,336,0), (280,224,0), (224,280,0), (336,280,0)
K, tenseur diffusion effective kg/(m s %), (Id identité dans R³).
- 4.7e-10 Id pour les couches du dessus et dessous
- 5.04e-10 Id pour la couche du centre
ß, tenseur de dilatation m / (m %) (Id identité dans R³)
- 7.5e-4 Id pour les couches du dessus et dessous
- 4.5e-4 Id pour la couche du centre.
E (module d'Young GPa)
- 4.05 pour les couches du dessus et dessous
- 2.5 pour la couche du centre.
v (coefficient de Poisson) 0.3 partout

La symétrie me permet d'utiliser un quart de la plaque donc une plaque 280 mm x 280 mm x 12 mm. Les conditions aux bords mécanique deviennent alors

U_X = 0 sur (280,y,z)
U_Y = 0 sur (x,280,z)

Pour la composante en Z j'ai simplifié la condition (permet d'avoir des maillages plus grossier) et j'ai pris

U_Z = 0 sur (280,280,0).

Quand on voudra faire la comparaison avec les résultats déjà obtenus on prendra/construira un maillage (Maillage2, Maillage5 ou Maillage8) permettant de mettre la vrai condition au bord i.e.

U_Z = 0 sur (280,280,0), (280,224,0), (224,280,0).

J'ai construit 2 séries de maillages: une série de 9 maillages de héxaèdres, une série de 9 maillages de prismes. Comme point de départ j'ai trois grilles 2D composées de quadrangles: une grille 8x8, une grille 5x5 et une grille 3x3. Les maillages prismes sont obtenus en coupant chaque quadrangle de la "base" en deux. Pour chacunes de ces grilles j'extrude de trois façons différentes: une extrusion avec 2 éléments pour les couches externes et 4 éléments pour le milieu, une extrusion avec 4 éléments pour les couches externeset 8 éléments pour le milieu et une extrusion avec 6 éléments pour les couches externeset12 éléments pour le milieu. Ce qui donne 9 maillages en héxaèdre. On produit 9 maillages en prisme en utilisant le même procédé à partir des trois grilles de triangles.

Données sur les maillages/éléments pour les héxaèdres

Maillage1: 144 sommets et 72 éléments. Les éléments sont des "briques" de dimensions (mm) 93.33 x 93.33 x 1 et 93.33 x 93.33 x 2.
Maillage2: 324 sommets et 200 éléments. Les éléments sont des "briques" de dimensions (mm) 56 x 56 x 1 et 56 x 56 x 2.
Maillage3: 729 sommets et 512 éléments. Les éléments sont des "briques" de dimensions (mm) 35 x 35 x 1 et 35 x 35 x 2.
Maillage4: 272 sommets et 144 éléments. Les éléments sont des "briques" de dimensions (mm) 93.33 x 93.33 x 0.5 et 93.33 x 93.33 x 1.
Maillage5: 612 sommets et 400 éléments. Les éléments sont des "briques" de dimensions (mm) 56 x 56 x 0.5 et 56 x 56 x 1.
Maillage6: 1377 sommets et 1024 éléments. Les éléments sont des "briques" de dimensions (mm) 35 x 35 x 0.5 et 35 x 35 x 1.
Maillage7: 400 sommets et 216 éléments. Les éléments sont des "briques" de dimensions (mm) 93.33 x 93.33 x 0.33 et 93.33 x 93.33 x 0.66.
Maillage8: 900 sommets et 600 éléments. Les éléments sont des "briques" de dimensions (mm) 56 x 56 x 0.33 et 56 x 56 x 0.66.
Maillage9: 2025 sommets et 1336 éléments. Les éléments sont des "briques" de dimensions (mm) 35 x 35 x 0.33 et 35 x 35 x 0.66.

Pour l'équivalent du Maillage1 en prisme cliquer ici, pour l'équivalent du Maillage2 en prisme cliquer ici. Pour finir j'ai construit deux maillages (héxaèdres et prismes) très fin (10x10 avec 8 -16 - 8 couches) servant à calculer des solutions de référence.

boisqref: 3993 sommets et 3200 éléments. Les éléments sont des briques de dimensions (mm) 28 x 28 x 0.25 et 28 x 28 x 0.5.

Cliquer ici pour les étape de la création du maillage avec iMEF++

Cliquer ici pour avoir accès aux fichiers de données (.mail, enti, geom, CL, champs)

Le code boisHM.cc est le code de base à partir duquel j'ai produit les résultats présenté ici. Il s'agit d'un seul code permettant de calculer M et U à chaque pas de temps (U est calculé en post-traitement du calcul de M). Ce code est tout à fait général et peux être utilisé dans le cas où le phénomène d'hystérèse est présent. Cliquer ici pour avoir accès au source du programme boisHM.cc

Schéma d'intégration

Dans le cas de maillages composé de prismes et de héxaèdres le jacobien associé à la géométrie n'est pas constant. Le choix de la formule de quadrature pour approximer les différentes intégrales devient alors plus ou moins complexe suivant le problème à résoudre, le type d'élément (effet du jacobien) et de l'interpolation choisie (effet des fonctions de base). Certains cas sont cependant "classiques", entre autres dans le cas d'interpolation linéaire et quadratique avec des prismes ou des héxaèdres pour les problème thermique (M) ou élastique (U). Le cas de l'interpolation de Crouzeix-Raviart est moins usité et c'est surtout pour lui que j'ai fait ces tests.

Remarque. Le nombre de points d'intégration à un effet direct sur le temps de calcul et d'assemblage. Le temps pour l'assemblage des matrices et vecteurs peux se décomposer en deux parties:

calcul des matrices élémentaires et vecteurs élémentaires
assemblages des matrices/vecteurs global.

Le changement de schéma d'intégration n'a pas d'effet sur 2. mais à un effet sur 1. Par exemple, pour un maillage d'héxaèdres le passage d'un schéma 3x3x3 (27 points) à un schéma 4x4x4 (64 points) multipliera par un facteur > 2 (en fait 64/27) le temps pour le calcul de chaque matrices/vecteurs élémentaire. Comme le problème en M est instationnaire et non-linéaire on devra obligatoirement passer par une boucle (pour le résolution) dans laquelle on devra faire un réassemble plus ou moins fréquent du système.

Mise en oeuvre

J'ai calculé la vraie solution plutôt que de faire des comparaisons au niveau élémentaire. J'ai fait le calcul sur un pas de temps de longueur 1.
J'ai pris comme schéma d'interpolation de référence les schémas: 3 x 3 x 3 pour les héxaèdres et 7 x 3/12 x 3 pour les prismes.
Pour évaluer la qualité des schémas j'ai décidé de comparer les solutions à chaque DDL pour chaque schéma par rapport à la solution obtenue avec le schéma de référence. Si on note $F^r$ la solution obtenue avec le schéma de référence et $F$ la solution obtenue avec un autre schéma, les valeurs affichées sont:
Abs(F) = $\max\limits_{i \in DDL} (| F^r_i - F_i |)$ et Rel(F) = $\max\limits_{i \in DDL} \frac{| F^r_i - F_i |}{max(|F^r_i|,|F_i |)}$
J'ai utilisé les maillages no. 5 (héxaèdres et prismes) pour créé les tableaux mais les résultats sont indépendant, qualitativement, du maillage employés. Le même comportement numérique se retrouve en employant n'importe quelle des maillages proposés.
Le code pour ces calculs est relativement simple: on choisi un schéma de référence puis on compare avec la solution pour les schémas N x N x N pour N = 3,4,5,6,7,8 pour les héxadres et 6 x N, 7 x N, 12 x N et 16 x N pour N = 3,4,5,6,7,8 pour les prismes (l'explication de ces schémas est donnée plus bas).

Héxaèdre

On regarde 3 interpolations différentes pour M: linéaire, quadratique et Crouzeix-Raviart et deux interpolations pour U: quadratique et Crouzeix-Raviart. On va tester différents schémas d'intégrations en faisant le produit de schéma d'intégration 1D. On aura des schéma de la forme N x N x N avec N = 3,4,5,6,7,8. A noter que le schéma ordinairement employé est le 3 x 3 x 3 pour M et U quadratique., c'est pourquoi j'ai choisi ce schéma comme référence.

	4 x 4 x 4		5 x 5 x 5		6 x 6 x 6
M	Abs	Rel (%)	Abs	Rel (%)	Abs	Rel (%)
Linéaire	3.e-14	4.e-13	5.e-14	7.e-13	4.e-14	5.e-13
Quadratique	2.e-13	2.e-12	3.e-13	4.e-12	2.e-13	3.e-12
Crouzeix-Raviart	8.e-14	1.e-12	8.e-14	1.e-12	2.e-13	3.e-12

	4 x 4 x 4		5 x 5 x 5		6 x 6 x 6
U	Abs	Rel (%) ¹	Abs	Rel (%) ¹	Abs	Rel (%) ¹
Quadratique	3.e-11	6.e-3	3.e-11	5.e-3	4.e-11	7.e-3
Crouzeix-Raviart	3.e-11	4.e-3	3.e-11	4.e-3	6.e-11	9.e-3

¹ L'erreur relative, Rel, n'est grande qu'en apparence. En effet les DDL où l'erreur relative est la plus grande sont des DDL où les déplacements sont très petits.

Conclusion: Si pour M on utlise une interpolation linéaire, quadratique ou Crouzeix-Raviart et pour U une interpolation quadratique ou Crouzeix-Raviart on peut prendre le schéma 3 x 3 x 3. Dans le cas où les différents coefficients ne sont plus constants on utilisera si nécessaire le schéma 4 x 4 x 4.

Prisme

Même procédé que pour les héxaèdres. Cependant on devra distinguer le problème en M et le problème en U. Pour M la formule de quadrature de référence est obtenue en faisant le produit d'un schéma 2D à 7 points sur les triangles et d'un schéma 1D à 3 points. C'est le schéma 7 x 3. Pour U la formule de référence sera 7 x 3 (quad) ou 12 x 3 (Crou.-Rav.). A fin de comparaison j'ai aussi fait les calculs avec un schéma 7 x 3 pour U.

6 x 3

6 x 4

7 x 3

7 x 4

12 x 3

12 x 4

16 x 3

16 x 4

M

Abs

Rel (%)

Abs

Rel (%)

Abs

Rel (%)

Abs

Rel (%)²

Abs

Rel (%)²

Abs

Rel (%)²

Abs

Rel (%)²

Abs

Rel (%)²

Lin.

2.e-14

3.e-13

2.e-14

3.e-13

0.

0.

2.e-14

3.e-13

2.e-14

4.e-13

2.e-14

3.e-13

2.e-14

4.e-13

2.e-14

3.e-13

Quad.

1.e-13

2.e-12

1.e-13

2.e-12

0.

0.

1.e-13

2.e-12

1.e-13

2.e-12

1.e-13

2.e-12

2.e-13

3.e-12

2.e-13

3.e-12

Crou.-Rav.

3.e-8

100.

5.e-8

100.

0.

0.

3.e-13

200.

3.e-13

200.

3.e-13

200.

5.e-13

200.

5.e-13

200.

² Pour les schéma "plus fin" que 7 x 3 l'erreur relative dans le cas Crou.-Rav. est trompeuse puisqu'il s'agit de très petites valeurs produisant ces résultats. En limitant les calculs avec ces valeurs, par exemple en imposant $\max(|F^r_i|,|F_i |) > 10^{-8} \max\limits_{i \in DDL} |F^r_i|$ (exclus les valeurs trop petites par rapport à M), on obtient des résultats similaires à ceux du cas quadratique.

6 x 3

6 x 4

7 x 3

7 x 4

12 x 3

12 x 4

16 x 3

16 x 4

U

Abs

Rel (%)³

Abs

Rel (%)³

Abs

Rel (%)³

Abs

Rel (%)³

Abs

Rel (%)³

Abs

Rel (%)³

Abs

Rel (%)³

Abs

Rel (%)³

Quad. (7 x 3)

2.e-11

4.e-3

3.e-11

6.e-3

0.

0.

1.e-11

3.e-3

3.e-12

6.e-4

2.e-11

4.e-3

6.e-12

1.e-3

3.e-11

6.e-3

Crou.-Rav. (7 x 3)

1.e-6

199

1.e-6

199.

0.

0.

2.e-11

8.

4.e-9

199.

4.e-9

199.

4.e-9

199.

4.e-9

199.

Crou.-Rav. (12 x 3)

1.e-6

200.

1.e-6

200.

4.e-9

199.

4.e-9

199.

0.

0.

5.e-12

0.2

6.e-11

1.6

2.e-11

0.6

³ Pour les schéma "plus fin" que le schéma de référence (7 x 3 ou 12 x 3) l'erreur relative dans le cas Crou.-Rav. est trompeuse puisque, là encore, il s'agit de très petites valeurs produisant ces résultats. En modifiant le procédé appliqué à M, les valeurs en U sont petites: $\max(|F^r_i|,|F_i |) > 10^{-8}$ (cette fois on exclus les valeurs inférieure à 1e-8 m), on obtient, avec le schéma 12 x 3, des résultats comparable au quadratique pour l'erreur relative. Dans le cas du schéma 7 x 3 les résultats sont alors améliorés mais pas suffisamment; de plus l'erreur absolue nous permettent de déterminer le bon schéma dans le cas Crou.-Rav.

Conclusion: Si pour M on utilise une interpolation linéaire, quadratique ou Crouzeix-Raviart et pour U on choisi une interpolation quadratique on peut prendre le schéma 7 x 3. Si pour U on a choisi une interpolation Crouzeix-Raviart on doit prendre pour U le schéma 12 x 3. Dans le cas où les différents coefficients ne sont plus constants on utilisera si nécessaire un schéma plus fin.

Interpolation

L'idée ici est de comparer les différentes solutions obtenues suivant le couple d'interpolation choisies pour M et U. Tout d'abord une remarque importante concernant M: puisqu'on a choisi de prendre une bonne approximation de M (quad. ou Crou.-Rav.) le changement d'interpolation d'ordre 2 n'a que peu ou pas d'impact sur la solution. Cela a pour effet de "découpler" l'interpolation en M et en U i.e. si on choisi une interpolation quad ou Crou.-Rav. pour M la qualité de l'approximation en U ne dépend plus que du choix d'interpolation de U. Naturellement cette remarque ne s'applique pas si on considère une interpolation linéaire pour M.

Remarque. Le type d'interpolation à un effet direct sur la taille des matrices élémentaore et assemblée. Ce qui à un effet double sur la résolution

espace requis pour la résolution
temps CPU

L'effet de l'interpolation ce fait évidemment via les degrés de libertés. Par exemple, pour un maillage de 81 héxaèdres le passage d'une interpolation quadratique à une interpolation Crouzeix-Raviart multipliera par un facteur >1.5 (en fait 931/544) le nombre de DDL. Comme le problème en M est instationnaire et non-linéaire on devra obligatoirement passer par une boucle (pour le résolution) dans laquelle on devra faire un réassemblage plus ou moins fréquent du système.

Mise en oeuvre

J'ai calculé la vraie solution plutôt que de faire des comparaisons au niveau élémentaire. J'ai fait le calcul sur 10 pas de temps de longueur 1.
J'ai pris comme schéma d'interpolation de référence les schémas: 3 x 3 x 3 pour les héxaèdres et 12 x 3 pour les prismes.
Pour comparer les interpolations j'ai décidé de calculer la différence à chaque DDL pour chaque pas de temps ainsi que la norme . Si on note $F^r (j\Delta t)$ et $F (j \Delta t)$ la solution au pas de temps "j" pour les interpolations "0" et "1" les valeurs affichées sont:
Abs = $\max\limits_{j \Delta t} \max\limits_{i \in DDL} (| F^r_i - F_i |)$ et
Le code pour ces calculs est relativement simple: on calcul la solution pour les deux interpolations choisies puis en postraitement on fait les calculs décrit plus haut.

Comparaison Quadratique/Crouzeix-Raviart

	Maillage 1		Maillage 2		Maillage 4		Maillage 5
M	Abs		Abs		Abs		Abs
Héxaèdres	3.e-13	6.e-16	3.e-13	7.e-16	4.e-13	6.e-16	4.e-13	7.e-16
Prismes	3.e-13	7.e-16	3.e-13	7.e-16	3.e-13	7.e-16	4.e-13	7.e-16

Conclusion: Pour les maillages choisis le type d'interpolation d'ordre deux n'a pas d'effet sur la solution obtenue. C'est pourquoi j'ai choisi d'utiliser l'interpolation quadratique (20/15 DDL par élément) plutôt que Crouzeix-Raviart (27/21 DDL par élément).

De cette conclusion il découle que le choix d'interpolation sur M n'a pas d'effet sur le choix d'interpolation de U. Dans ce cas on a:

	Maillage 1		Maillage 2		Maillage 4		Maillage 5
U	Abs		Abs		Abs		Abs
Héxaèdres	6.e-8	5.e-10	4.e-8	2.e-10	1.e-8	5.e-11	8.e-9	3.e-11
Prismes	2.e-7	1.e-9	2.e-7	2.e-9	9.e-8	7.e-10	1.e-7	9.e-10

Conclusion: L'effet du type d'interpolation ne peux pas être considéré nul dans ce cas mais il y a lieu de se demander si il est nécessaire d'utiliser une interpolation contenant un si grand nombre de degrés de libertés (81/63 par éléments pour Crou.-Rav. au lieu de 60/45 pour quad.) pour un gain aussi petit.

Convergence et adaptation

On regarde ce qui se passe du point de vue convergence de la solution éléments finis et ce qui pourrait être fait relativement à l'adaptation du maillage. Premièrement il est clair que dans ce cas-ci M ne dépendra que de z (la variation sera nulle dans le plan XY puisque les matériaux sont homogènes). J'ai regardé l'effet du raffinement du maillage dans l'épaisseur et dans le plan séparément. J'ai travaillé avec les maillages de quadrangles, la solution de référence est obtenue en utilisant le maillage boisqref. J'ai pris 10 pas de temps de longueur 100 et

pour comparer les solutions. Cliquer pour agrandir.

Conclusion: La solution, M et U, ne dépend pas du raffienement du maillage de la base mais de la discrétisation dans l'épaisseur. Puisque M ne dépend que de z un raffinement dans le plan xy n'aura aucun effet sur M mais alors puisque U ne varient pas avec un raffinement dans le plan j'en conclus que la solution calculée est exacte i.e. U est quadratique.

Résultats et essais mpeg

Pour finir avec le cas constant voici mes premiers tests de mpeg pour visualiser la solution pour un intervalle de temps arbitraire. Ici j'ai choisi de calculer la solution pour 120000 sec. (à peu près 33 h. comme dans les articles). J'ai choisi de prendre 160 pas de temps de 750 sec. (12.5 min).

Les résultats semblent concordés (c'est à vérifier) sauf pour U où j'obtiens un déplacement maximal de 0.0125 m plutôt que 0.00189 m dans les articles.