Page concernant la formulation et la notation utilisée pour l'implémentation dans MEF++ de codes pour la résolution E.F. des problèmes hygro-mécanique (HM) et thermo-hygro-mécanique (THM)

Cette page se veut une référence pour la compréhension des codes que j'ai construit et pour "fixé" la notation utilisé. Naturellement la notation ne suit peut être pas les conventions établies mais elle sera respectée dans tout les codes.

EDP


Les trois variables d'états sont T pour la température,  M pour l'humidité et U pour les déplacements. D'abord les 3 équations régissant le phénomène:

\[  C_T \frac{\partial T}{\partial t} - \nabla \cdot ( K_T \nabla T) = 0  \]      (1)
\[  C_M \frac{\partial M}{\partial t} - \nabla \cdot ( K_M \nabla M + K_{TM} \nabla T) = 0  \]      (2)
\[  \rho \frac{\partial^2 U}{\partial t^2} - \nabla \cdot (E(\epsilon - \beta_M \Delta M - \beta_T \Delta T)) = 0  \]      (3)




\[  C_T = C_T(\rho,T,M)  \]      \[  C_M = C_M(\rho,T,M)  \]

\[  K_T = K_T(\rho,T,M)  \]      \[  K_M = K_M(\rho,T,M)  \]      \[  K_{TM} = K_{TM}(\rho,T,M)  \]

\[  E = E(\rho,T,M)  \]      \[  \epsilon_{ij} = \frac{1}{2}(U_{i,j} + U_{j,i})  \]      \[  \beta_T = \beta_T(\rho,T,M)  \]      \[  \beta_M = \beta_M(\rho,T,M)  \]


Un problème HM correspond à un problème où l'on suppose que la température est constante. Dans ce cas le gradient de la température est nulle tout comme les variations de température. On se retrouve avec le système (2)-(3) à deux variables à résoudre. D'un point de vue numérique la présence de la température n'a pas d'impact quand à l'approche utilisée pour résoudre. On résoud à chaque pas de temps pour le couple (T,M) puis on résoud le problème d'élasticité. Le seul effet étant l'augmentation du temps CPU, l'équation en T est du même type que l'équation en M on peut donc dire que l'on double le temps de calcul en passant de la résolution de M seule à la résolution de T et M.

D'un point de vue pratique on peut utiliser le même code pour les deux types de problèmes puisque L'équation (1) correspond dans MEF++ à l'ajout de terme de formulation pour le problème associé à (2) tout comme les termes en T dans (2) et (3).

L'équation (3) peut être résolue de deux manières:
  1. dynamique on ne fait aucune modification à l'équation (les effets d'inertie ne sont pas négligé) dans ce cas \[  \Delta M = (M - M_0)  \]     \[  \Delta T = (T - T_0)  \]
  2. quasi-statique (ou incrémental) l'équation ne dépend plus explicitement du temps, à chaque pas de temps l'incrément des déplacements satisfait l'équation (l'indice p dénote la valeur au pas de temps précédent):
\[  \nabla \cdot (E(\epsilon - \beta_M \Delta M - \beta_T \Delta T)) = 0  \]
  \[  \Delta M = (M - M_p)  \]     \[  \Delta T = (T - T_p)  \]

Conditions initiales et aux bords

panneau
Conditions initiales:

 \[  T = T_0  \]      \[  M = M_0  \]      \[  U = 0  \]      \[  \sigma = 0  \]      \[  \sigma = E(\epsilon - \beta_M \Delta M - \beta_T \Delta T)  \]

Deux approches: sans contact ou avec contact

Jusqu'ici le comportement mécanique a toujours été décrit sans tenir compte du contact ce qui fait que l'on avait en fait une plaque supportée uniquement en son centre, sur le point \[  \Gamma_0  \] plus précisément. Avec l'addition de terme de formulation et d'un solveur permettant de tenir compte du contact on peux maintenant formuler exactement le problème. Je garde quand même les deux formulations jusqu'a ce que je sois sur d'avoir une résolution du problème de contact stable et validé.

Conditions aux bords:
(HM sans contact)     \[  q_M = h(M - M_{\infty})  \]      sur \[  \Gamma_b  \]    \[  U = 0  \]      sur \[  \Gamma_0  \]           \[  U_x = 0  \]      sur \[  \Gamma_{sx}  \]      \[  U_y = 0  \]     sur \[  \Gamma_{sy}  \]        
(HM avec contact)     \[  q_M = h(M - M_{\infty})  \]      sur \[  \Gamma_b  \]      \[  U_z \ge 0  \]      sur \[  \Gamma_b  \]        \[  U_x = 0  \]      sur \[  \Gamma_{sx}  \]      \[  U_y = 0  \]     sur \[  \Gamma_{sy}  \]    

(THM)    
\[  t \in [0,4\ min[  \]

 \[  \sigma_n = P  \]      \[  T = T_P  \]    sur \[  \Gamma_h  \] et \[  \Gamma_b  \]        \[  U_x = 0  \]      sur \[  \Gamma_{sx}  \]      \[  U_y = 0  \]     sur \[  \Gamma_{sy}  \]    

\[  t \in [4\ min, \infty  \]     

(sans contact)   \[  U = 0  \]      sur \[  \Gamma_0  \]           \[  U_x = 0  \]      sur \[  \Gamma_{sx}  \]      \[  U_y = 0  \]     sur \[  \Gamma_{sy}  \]     
(avec contact)  \[  U_z \ge 0  \]        sur \[  \Gamma_b  \]        \[  U_x = 0  \]      sur \[  \Gamma_{sx}  \]      \[  U_y = 0  \]     sur \[  \Gamma_{sy}  \]   


Condition au bord supplémentaire pour l'approche quasi-statique avec contact

Dans le cas avec contact si on utilise une approche quasi-statique le problème n'est pas bien posé puisqu'il n'y a pas unicité de la solution (la matrice a un pivot nul ou presque). En effet il nous manque une condition éliminant les mouvements de corps rigide dans la direction vertical. On devra donc ajouter une condition au bord fixant un DDL de la "composante z" des déplacements. Si l'on utilise une approche totalement dynamique le problème est bien posé car la matrice masse a pour effet de rendre le système bien conditionné.

Dans le cas du problème HM on peut ajouter le blocage de tout les DDL sur
\[  \Gamma_0  \] . Cette condition a cependant un impact important sur le comportement du panneau. Naturellement dans le cas du THM il n'est pas question d'imposer une telle condition puisque l'on veut avant tout modéliser un processus de fabrication existant. On veut alors imposer une condition qui n'aura pas d'effet sur la solution.

Je propose d'ajoute comme condition que le panneau est toujours en contact avec le plan xy ce qui se traduit par:

\[  \forall t \in [4\ min, \infty\ \ \ \exists X \in \Gamma_b\ tel\ que\ U_z(X) = 0  \]

On approxime cette condition en travaillant sur les DDL de U . En utilisant l'hypothèse que les déplacements verticaux sont positifs ou nul on peut alors traduire la condition de la façon suivante:

\[  \min\limits_{DDL} U_z = 0 \leftrightarrow \prod\limits_{DDL} U_z = 0  \]

Cette dernière condition étant traitée comme une condition de contact plutôt qu'un condition au bord. C'est-à-dire que cette condition ne sera pas utilisé lors de l'assemblage de la matrice mais plutôt lors de la résolution du problème de contact, ce sera une contrainte (au sens de l'optimisation) imposée à la solution. On favorise la deuxième forme de la condition sur les DDL car elle peut s'exprimer sous la forme d'un produit scalaire et s'ajoute de façon relativement simple au code déjà en place dans MEF++ pour les problèmes de contacts unilatéraux.