# Résolution d'un problème multiphysique couplé par décomposition en fonctions à variables séparées

Marianne Béringhier, Mikaël Gueguen, Jean-Claude Grandidier

► To cite this version:

Marianne Béringhier, Mikaël Gueguen, Jean-Claude Grandidier. Résolution d'un problème multi-physique couplé par décomposition en fonctions à variables séparées. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01408465

HAL Id: hal-01408465

https://hal.science/hal-01408465v1

Submitted on 5 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Résolution d'un problème multiphysique couplé par décomposition en fonctions à variables séparées.

M. Béringhier¹, M. Gueguen¹, J-C. Grandidier¹

¹Laboratoire de Mécanique et Physique des Matériaux (ENSMA – UMR CNRS 6617)
Téléport 2 – 1 av. Clément Ader – BP 40109 – F – 86961 Futuroscope Cédex
marianne.beringhier@lmpm.ensma.fr

**Résumé** – Les modèles de couplage multiphysique peuvent être résolus par la méthode des éléments finis. Cette méthode connaît certaines limites, principalement dans le cas de problèmes multiphysiques à temps caractéristiques fortement différents. Dans ce papier, l'application de la méthode de décomposition en fonctions à variables séparées, présentée dans [3] pour un problème de thermique 2D, est étudiée et discutée dans le cas d'un problème couplé thermomécanique 1D.

**Mots clés** – Problème multiphysique, séparation de variables, temps caractéristiques.

## 1 Introduction

La résolution des modèles de couplage multiphysique par la méthode des éléments finis nécessite un maillage adapté pour déterminer numériquement et précisément les états transitoires de température ou de diffusion d'espèces dans un solide. Dans les logiciels industriels classiques, afin de satisfaire la stabilité des schémas d'intégration temporelle utilisés, il faut ajuster la densité de maillage au pas de temps (condition CFL), or ce dernier est directement lié au temps caractéristique des processus physiques. Lorsque les problèmes multiphysiques couplés contiennent différents temps caractéristiques (par exemple dans le cas où l'on traite des équations de thermique, de diffusion et de mécanique avec une loi de comportement viscoélastique), il faut donc ajuster le pas de temps pour l'ensemble de ces phénomènes transitoires et par conséquence le maillage pour chacune de ces équations. A cette fin, une méthode consiste à découpler la résolution numérique des différentes équations de la physique. Néanmoins, quand les couplages deviennent très forts, c'est-à-dire par exemple quand les paramètres dépendent non linéairement des champs physiques calculés, la résolution découplée n'est plus raisonnable. Il faut alors résoudre simultanément les différentes équations. Toutefois, si l'on utilise un seul maillage avec des éléments possédant aux nœuds toutes les inconnues du problème physique, cela mène à des difficultés dans le choix de l'incrément de temps, à des temps de calcul importants ou voire à des impossibilités de convergence.

L'objectif de notre recherche s'inscrit dans cette problématique et porte sur la résolution des problèmes couplés complexes (par exemple [1] : thermique, schéma mécanistique, diffusion d'espèces et mécanique), et il est proposé dans ce papier d'utiliser une méthode alternative à la méthode des éléments finis, méthode initiée par Ammar et al. [2, 3], qui consiste à rechercher la solution sous la forme d'une somme de fonctions à variables séparées. Elle s'inscrit dans le cadre théorique des méthodes de Ritz et elle a été détaillée dans [3] pour la résolution d'un problème de thermique 2D. Nous présentons ci-dessous une extension de la méthode à un système couplé thermomécanique avec une loi de comportement viscoélastique en se limitant à un cas unidimensionnel.

La section 2 introduit les équations du problème ainsi que le formalisme général de la
méthode de résolution. Les particularités de la méthode sont explicitées dans le cas de la résolution d'un problème couplé. La section 3 illustre les résultats obtenus par la méthode dans le cas de la résolution d'un problème thermomécanique couplé. L'outil numérique basé sur la méthode est ensuite validé et discuté par comparaison aux résultats obtenus à l'aide de la méthode classique des éléments finis.

# 2 Présentation de la méthode de décomposition en fonctions à variables séparées dans le cas d'un problème couplé

## 2.1 Les équations du problème couplé

Dans le cadre de notre étude, nous considérons un problème couplé thermomécanique unidimensionnel en espace (variable notée x) avec uniquement un couplage de type thermodynamique (prise en compte de la dilatation thermique sur la mécanique et la thermique). La forme générique du problème couplé est la suivante :

$$ (\lambda+2 \mu) \frac{\partial^2 u}{\partial x^2} + (\lambda_v + 2 \mu_v) \frac{\partial^2}{\partial x^2} \left( \frac{\partial u}{\partial t} \right) - (3\lambda + 2\mu)\alpha_T \frac{\partial T}{\partial x} + f = 0 \quad (1) $$

$$ C \frac{\partial T}{\partial t} - \lambda_T \frac{\partial^2 T}{\partial x^2} + (3\lambda + 2\mu)\alpha_T T \frac{\partial}{\partial t} \left( \frac{\partial u}{\partial x} \right) - r = 0 \quad (2) $$

où dans l'équation (1) de la mécanique, $\lambda$, $\mu$, et $\lambda_v$, $\mu_v$ sont les coefficients de Lamé liés respectivement à l'élasticité et à la viscoélasticité, $\alpha_T$ le coefficient de dilatation, $f$ le chargement mécanique et, dans l'équation (2) de la thermique, $C$ est la chaleur spécifique, $\lambda_T$ la conductivité thermique et $r$ la source de chaleur. Tous ces coefficients sont considérés ici constants.

Le domaine est $\Omega = (0, L_x) \times (0, L_t)$ et on considère des conditions aux bords et initiales nulles pour le déplacement $u$ et la température $T$. Ce système présente deux temps caractéristiques : le rapport adimensionnalisé pour une épaisseur unité $C/\lambda_T$ pour la thermique et le rapport $(\lambda_v + 2 \mu_v)/(\lambda + 2 \mu)$ pour la mécanique.

## 2.2 Résolution des équations

### 2.2.1 La méthode de décomposition en fonctions à variables séparées

La méthode consiste à exprimer la solution comme une combinaison linéaire de fonctions à variables séparées, qui forment une famille que l'on enrichit par itérations successives, les coefficients de la combinaison étant ajustés à chaque itération [3]. Nous nous intéresserons ici uniquement aux spécificités liées aux termes de couplage. Pour plus de détails sur la méthode, le lecteur pourra se référer à [2, 3].

### 2.2.2 Les spécificités liées à la résolution d'un problème couplé

Le couple solution est recherché sous la forme :

$$ u(x,t) = \sum_{i=1}^{n} \alpha_i A_i(x) B_i(t), \quad T(x,t) = \sum_{i=1}^{n} \beta_i D_i(x) E_i(t). $$

Du point de vue de la discrétisation, on utilise des fonctions d'interpolation de degré 2 nulles aux bords du domaine pour l'espace et des fonctions d'interpolation d'ordre 1 nulles en $t=0$ pour le temps. Le nombre de points de discrétisation peut être différent pour la température et le déplacement. Pour les fonctions dépendantes du temps, les points peuvent être choisis en fonction des temps caractéristiques des phénomènes physiques.

A l'itération N, pour enrichir la famille, on cherche des fonctions $R(x)$, $P(t)$, $V(x)$ et $W(t)$ telles que
$$u(x,t) = \sum_{i=1}^{N} \alpha_i A_i(x) B_i(t) + R(x)P(t), \quad T(x,t) = \sum_{i=1}^{N} \beta_i D_i(x) E_i(t) + V(x)W(t).$$

En choisissant comme champ virtuel (fonctions tests) la fonction suivante :

$$u^*(x,t) = R^*(x)P(t) + R(x)P^*(t), \quad T^*(x,t) = V^*(x)W(t) + V(x)W^*(t),$$

nous obtenons la formulation intégrale (3)-(4) correspondant au problème (1)-(2) :

$$\int_{\Omega} \left[ - \left( (\lambda+2\mu)\frac{\partial u}{\partial x} + (\lambda_v+2\mu_v)\frac{\partial}{\partial x}\left(\frac{\partial u}{\partial t}\right) \right) \cdot \frac{\partial u^*(x,t)}{\partial x} - (3\lambda+2\mu)\alpha_T \frac{\partial T}{\partial x} u^* + \rho f u^* \right] d\Omega = 0 \quad (3)$$

$$\int_{\Omega} \left[ \lambda \frac{\partial T}{\partial x} \cdot \frac{\partial T^*}{\partial x} + \rho C \frac{\partial T}{\partial t} T^* + (3\lambda + 2\mu) \alpha_T \frac{\partial}{\partial t} \left( \frac{\partial u}{\partial x} \right) T \cdot T^* - r T^* \right] d\Omega = 0 \quad (4)$$

Le système non linéaire en R,P,V,W qui en résulte peut se mettre sous la forme générique suivante :

$$F_u(R, P, V, W) = 0$$

$$F_T(R, P, V, W) = 0$$

La décomposition en sous-systèmes linéaires proposée dans [3] pour la résolution de ce système non linéaire ne peut plus s'appliquer à cause des termes de couplage. Pour résoudre plus efficacement ce système non linéaire, nous l'avons approché par le système partiellement couplé suivant :

$$F_u(R, P, V=0, W=0) = 0$$

$$F_T(R=0, P=0, V, W) = 0$$

où $F_u$ est non linéaire en R, P et $F_T$ en V,W.

Cela correspond à un couplage faible car, pour déterminer R et P (respectivement V et W), la fonction de T (respectivement u) calculée à l'itération précédente est utilisée. Ce système qui reste non linéaire a été résolu à l'aide d'un algorithme de Newton modifié pour pallier le défaut d'unicité de la solution du système non linéaire.

Alors que l'étape précédente est partiellement découplée, on recherche les coefficients $\{\alpha_i, \beta_i\}_{i=1,N+1}$ comme solutions du système complètement couplé (3)-(4) avec :

$$u(x,t) = \sum_{i=1}^{N+1} \alpha_i A_i(x) B_i(t), \quad T(x,t) = \sum_{i=1}^{N+1} \beta_i D_i(x) E_i(t) \text{ et}$$

$$u^*(x,t) = \sum_{i=1}^{N+1} \alpha_i^* A_i(x) B_i(t), \quad T^*(x,t) = \sum_{i=1}^{N+1} \beta_i^* D_i(x) E_i(t).$$

En ce qui concerne l'influence du couplage sur la nature du système à résoudre pour calculer les coefficients $\{\alpha_i, \beta_i\}_{i=1,N+1}$, on remarque que c'est pour un $\alpha_T$ non nul dans (2) que le système (1)-(2) est non linéaire en u, T. Dans un tel cas, le système dont $\{\alpha_i, \beta_i\}_{i=1,N+1}$ est solution est également non linéaire. Sa résolution ne pose pas de problème particulier et se fait à l'aide d'un algorithme de Newton classique.

Avec la solution de u et T, le résidu est calculé pour le déplacement à partir de l'équation (1) et celui pour la température à partir de l'équation (2). On stoppe les itérations lorsque ces deux résidus sont inférieurs à une valeur donnée.
En comparaison de [3] ou d'un problème complètement découplé, le calcul du résidu nécessite de stocker un grand nombre de tableaux de grandes tailles. Dans le cas présent, ils sont à 4 entrées, liés aux calculs des intégrales des produits de 4 fonctions de forme alors que dans les cas découplés, seuls des tableaux à 2 entrées devaient être stockés. De même pour le calcul des coefficients et de l'enrichissement de la famille, nous devons stocker des tableaux à 3 entrées au lieu de 2 dans les cas découplés. Une gestion adéquate du stockage s'impose dans le cas couplé.

# 3 Application à un problème thermomécanique couplé

## 3.1 Présentation des simulations

Une barre unidimensionnelle de longueur 10 mm soumise à un chargement réparti uniformément $f = 1 \text{ N.mm}^{-1}$ et une source de chaleur homogène de $1 \text{ W.mm}^{-1}$ est étudiée. Les paramètres matériau utilisés pour la simulation sont donnés dans le tableau 1. La simulation est comparée à la solution obtenue avec la méthode des Eléments Finis, réalisée à l'aide du code de calcul ABAQUS Standard. Les couplages ainsi que la loi de comportement viscoélastique ont été implémentés en imbriquant un sous programme utilisateur HETVAL avec le sous programme UMAT. Ce développement s'est avéré nécessaire car ABAQUS n'utilise pas tout à fait la même loi viscoélastique et ne prend pas en compte tous les couplages du problème (1) et (2). Pour que les résultats soient comparables, les conditions aux limites ont été choisies pour avoir uniquement une déformation longitudinale dans la barre.

<table><caption>Tableau 1 – Paramètres matériau</caption><thead><tr><th>λ</th><th>μ</th><th>λ<sub>v</sub></th><th>μ<sub>v</sub></th><th>α<sub>T</sub></th><th>C</th><th>λ<sub>T</sub></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>0.1</td><td>1</td><td>1</td></tr></tbody></table>

Pour les deux simulations, le temps total est de 50 s. Dans le cas de la simulation avec la méthode de décomposition en fonctions à variables séparées (ML), 6 points sont utilisés en espace et 25 points en temps, soit un pas d'espace de 1 mm et un pas de temps de 2,04 s. Avec la méthode des Eléments Finis (EF), nous avons le même pas d'espace mais un pas de temps variable, beaucoup plus fin. Il a été demandé au logiciel une forte précision dans l'incrément de température autorisé pour chaque itération temporelle. Ce choix a été réalisé pour que le résultat de la simulation Eléments Finis serve de solution de référence, dans la mesure où pour le cas complètement couplé, nous ne disposons pas de solution analytique.

## 3.2 Résultats

Avant de présenter les résultats dans le cas complètement couplé, notons que la méthode de décomposition en fonctions à variables séparées a été validée et a donné des résultats satisfaisants dans le cas de problèmes unidimensionnels thermomécaniques découplés complètement voire partiellement couplés. La méthode s'est avérée tout à fait robuste pour résoudre l'équilibre mécanique avec un comportement mécanique viscoélastique, lorsque l'on sollicite la barre avec une force en échelon.

En ce qui concerne le problème complètement couplé, les résultats en terme de déplacement et de température sont représentés sur les graphes 1 et 2 respectivement. Sur ces 2 graphes, la figure b) détaille des sections particulières des surfaces représentées sur la figure a). La méthode de décomposition en fonctions à variables séparées est validée dans le cas complètement couplé. En effet, les différences obtenues entre la solution ML et la solution EF sont très faibles et principalement liées à la discrétisation utilisée dans la méthode ML. Effectivement, nous avons volontairement limité le nombre de points pour tester la robustesse de la ML, et la solution a été obtenue avec seulement 25 points, et cela donne des résultats tout à fait satisfaisants. Les différences sont plus importantes pour le déplacement, du fait que le problème de mécanique découplé possède un temps caractéristique de 8 s environ alors que le problème de thermique à
un temps caractéristique intrinsèque de 50 s environ. Il n'y a donc que 4 points pour décrire le transitoire dans le cas du problème de mécanique alors qu'il y en a 25 dans le cas du problème de thermique, malgré cela la méthode donne un résultat probant. Evidemment cet écart pourrait être diminué en utilisant un maillage différent en temps pour le déplacement et pour la thermique. En effet la possibilité de créer un maillage non structuré est une stratégie envisageable pour mieux décrire le transitoire dans le contexte des problèmes à différents temps caractéristiques.

<figure><img src="image_2.png" /><figcaption>Graphe 1 : Déplacement</figcaption></figure>

a : solution obtenue avec la méthode de décomposition en fonctions à variables séparées (ML).

b : comparaison de la solution ML avec la solution Abaqus (EF).

<figure><img src="image_3.png" /><figcaption>Graphe 2 : Température</figcaption></figure>

a : solution obtenue avec la méthode de décomposition en fonctions à variables séparées (ML).

b : comparaison de la solution ML avec la solution Abaqus (EF).
# 4 Discussion et conclusion

La méthode de décomposition en fonction à variables séparées a été validée dans le cas d'un problème couplé thermomécanique unidimensionnel en espace. Le temps de calcul des simulations n'a pas encore été étudié mais il est clair que ce dernier peut être optimisé en considérant des discrétisations adaptées aux temps caractéristiques des phénomènes.

Si l'on considère des couplages concernant les caractéristiques du matériau (par exemple avec un module d'Young dépendant de la température), la taille des tableaux à stocker augmente et une étude pour l'optimisation du stockage des tableaux mérite d'être menée. En comparaison, le passage à un espace tridimensionnel n'augmente pas la taille des tableaux mais seulement leur nombre et peut être envisagé sereinement.

Dans le cas de dépendance des propriétés matériau, si ces dernières sont données par des fonctions non linéaires, leur intégration dans la méthode ne devrait pas poser de problème, en revanche si elles ont données par table (ce qui est souvent le cas dans l'industrie), un pré traitement spécifique doit être envisagé.

L'application de la méthode à des systèmes fortement couplés semble être prometteuse et envisageable. Elle devrait permettre entre autre de réduire le temps de calcul des simulations (par rapport à la méthode des Eléments Finis) mais aussi de traiter des problèmes couplés avec des temps caractéristiques fortement différents qui conduisent parfois à des impasses dans la résolution des problèmes.

# Références

[1] G. Rambert, J-C. Grandidier, An approach to the coupled behaviour of polymers subjected to a thermomechanical loading in a gaseous environment. *European Journal of Mechanics A/Solids* 24 (1), 151-168, 2005.

[2] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A new family of solvers for some classes of multidimensionnal partial differential equations encountered in kinetic theory modeling of complex fluids. *J. Non-Newtonian Fluid Mech.* 139, 153-176, 2006.

[3] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A new family of solvers for some classes of multidimensionnal partial differential equations encountered in kinetic theory modelling fluids. Part II : Transient simulation using space-time separated representations. *J. Non-Newtonian Fluid Mech.* 144, 98-121, 2007.