r5.02.04 Thermique non linéaire en repère mobile#

Résumé

On présente la formulation et l’algorithme du problème de convection-diffusion en thermique non linéaire stationnaire introduit au sein de la commande THER_NON_LINE_MO [U4.33.04].

Le but est de résoudre l’équation de la chaleur dans un référentiel mobile lié à un chargement et se déplaçant dans une direction et à une vitesse données.

Les non linéarités du problèmes proviennent aussi bien des caractéristiques du matériau qui dépendent de la température, que des conditions aux limites de type rayonnement.

Les problèmes de ce type peuvent être traités avec des modèles utilisant des éléments finis de structure plans, axisymétriques et tridimensionnels.

Table des matières

Conditions aux limites. Problème de référence à résoudre#

On se reportera, par exemple, à [R5.02.01] pour plus d’informations sur les conditions aux limites thermiques de type Dirichlet, Neumann et Fourier linéaire, et à [R5.02.02] pour les conditions aux limites de type flux normal non linéaire (Fourier non linéaire).

En formulation enthalpique, le problème de thermique stationnaire consiste donc à résoudre dans un domaine \(O\) de frontière sur \(\partial \Omega\) .

\(V.\text{grad}u(T)-div(k(T)\text{grad}T)=Q\) dans \(\Omega\) éq 2-1

avec \(k(T)\frac{\partial T}{\partial n}=\gamma ({T}_{\text{ext}}-T)\) sur \({\partial}_{1}\Omega\) éq 2-2

\(k(T)\frac{\partial T}{\partial n}={q}_{0}\) sur \({\partial}_{2}\Omega\) éq 2-3

\(k(T)\frac{\partial T}{\partial n}=\alpha (T)\) sur \({\partial}_{3}\Omega\) éq 2-4

\(T={T}_{0}\) sur \({\partial}_{4}\Omega\) éq 2-5

où:

  • \({T}_{0}\) : est la température imposée sur \({\partial}_{4}\Omega\) ;

  • \({q}_{0}\) : est le flux normal imposé sur \({\partial}_{2}\Omega\) ;

  • \(\gamma\) : est le coefficient d’échange thermique ;

  • \({T}_{\text{ext}}\) : est la température extérieure ;

  • \(\alpha (T)\) : est le flux normal de type Fourier non linéaire (rayonnement).

Les équations [éq 2-2], [éq 2-5] représentent les conditions aux limites de types, respectivement : Fourier linéaire, Neumann, Fourier non linéaire et Dirichlet.

Le problème de référence [éq 2-1], [éq 2-5] est fortement non linéaire en raison des non linéarités sur \(k(T),u(T)\) (changement de phase) et \(\alpha (T)\) (rayonnement).

Formulation variationnelle du problème#

Soit \(\Omega\) un ouvert de \({R}^{3}\) , de frontière \(\partial \Omega ={\partial}_{1}\Omega \cup {\partial}_{2}\Omega \cup {\partial}_{3}\Omega \cup {\partial}_{4}\Omega\) telle que,

pour \(i\mathrm{\ne }j\) et \(i,j=1,\mathrm{...},4\) , on a: \({\partial}_{i}\Omega \cap {\partial}_{j}\Omega =\varnothing\) .

Soit encore \(\psi\) une fonction suffisamment régulière qui s’annule sur \({\partial}_{4}\Omega\) : \(\psi \in V=\left\lbrace y\text{régulière}\text{et}\psi {\mid}_{{\partial}_{4}\Omega }=0\right\rbrace\) .

Multiplions par \(\psi\) les deux membres de l’équation [éq 2-1], puis intégrons sur \(\Omega\) . Une intégration par parties donne alors:

\(\begin{array}{ccc}\underset{\Omega}{\int}\mathrm{Qyd}\Omega & =\underset{\Omega}{\int}V.\text{grad}u(T)\mathrm{yd}\Omega -& \underset{\Omega}{\int}div(k(T)\text{grad}T)\mathrm{yd}\Omega \\ & =\underset{\Omega}{\int}V.\text{grad}u(T)\mathrm{yd}\Omega +& \underset{\Omega}{\int}k(T)\text{grad}T.\text{grad}\mathrm{yd}\Omega -\underset{\partial \Omega -{\partial}_{4}\Omega }{\int}(k(T)\frac{\partial T}{\partial n}y)\mathrm{dG}\end{array}\) éq 3-1

puisque \(\psi\) est nul sur \({\partial}_{4}\Omega\) .

D’où, en tenant compte des conditions aux limites [éq 2-2], [éq 2-3] et [éq 2-4], la formulation variationnelle du problème de référence qui est donnée par l’équation suivante:

\(\forall \psi \in V\)

\(\begin{array}{}\begin{array}{cc}\underset{\Omega}{\int}k(T)\text{grad}T.\text{grad}\mathrm{yd}\Omega & +\underset{\Omega}{\int}V.\text{grad}u(T)\mathrm{yd}\Omega +\underset{{\partial}_{1}\Omega }{\int}\mathrm{gTydG}-\underset{{\partial}_{3}\Omega }{\int}\alpha (T)\mathrm{ydG}\end{array}\\ =\underset{\Omega}{\int}\mathrm{Qyd}\Omega +\underset{{\partial}_{1}\Omega }{\int}{\gamma T}_{\text{ext}}\mathrm{ydG}+\underset{{\partial}_{2}\Omega }{\int}{q}_{0}\mathrm{yd}\Omega ,\end{array}\) éq 3-2

Traitements des non linéarités#

En vue de la résolution numérique du problème non linéaire que nous considérons, il est nécessaire de traiter toutes les non linéarités.

Dans notre cas, citons la forte non linéarité liée à la fonction enthalpie \(u(T)\) qui prend en compte le changement de phase solide-liquide, ainsi que la non linéarité liée à la présence éventuelle d’une condition aux limites de flux normal non linéaire (rayonnement).

Rappelons que dans le cas classique des problèmes de thermique transitoire non linéaires sans convection, soit \(V=0\) , plusieurs méthodes de résolution sont proposées dans la littérature. Il existe aussi bien des méthodes utilisant des formulations enthalpiques que des méthodes utilisant des formulations en température, toutes ayant pour but de traiter au mieux la non linéarité liée à l’enthalpie (changement de phase).

Nous renvoyons le lecteur à la référence [bib5] pour un résumé des principales méthodes rencontrées dans la littérature. Toutefois, notons qu’en raison de la difficulté liée à la présence du terme de transport \(V.\mathit{grad}u(T)\) dans le problème, aucune de ces méthodes ne sera employée dans la suite.

Comme dans tout processus itératif, le but du schéma numérique en vue est de trouver un champ de température \({T}^{n+1}\) à l’itération \(n+1\) , à partir du champ de température \({T}^{n}\) , solution de l’itération précédente.

Traitement de la non linéarité liée à l’enthalpie#

Afin de traiter cette non linéarité, la stratégie employée dans cette étude a été inspirée d’une technique de résolution des problèmes de frontières libres [bib3], qui, à l’origine a été proposée dans [bib4].

Considérons la fonction enthalpie \(u(T)\) comme étant donnée sous une forme réciproque: Température fonction de l’enthalpie (inverse de la fonction \(u(T)\) ). En d’autre termes on aura à traiter la relation Température-enthalpie suivante :

\(T=\tau (u)\) éq 4.1-1

La raison de ce choix sera plus claire dans ce qui suit. En effet nous aurons à traiter un problème à deux champs : un champ de température et un champ enthalpie. La discrétisation de la fonction inverse [éq 4.1-1] permet d’incrémenter le champ enthalpie en fonction du champ actuel de température (et non l’inverse) comme suit :

Le développement au premier ordre de la fonction \(\tau (u)\) est le suivant,

\({T}^{n+1}=\tau ({u}^{n})+\tau '({u}^{n})({u}^{n+1}-{u}^{n}),\) éq 4.1-2

où \(\tau '\) est la dérivée de la fonction définie par [éq 4.1-1] par rapport à son argument.

Afin de prendre en compte cette non linéarité, et à partir de [éq 4.1-2], on remplace \({u}^{n+1}\) par une approximation en fonction du champ de température inconnu \({T}^{n+1}\) de la façon suivante:

\({u}^{n+1}-{u}^{n}=\omega ({T}^{n+1}-\tau ({u}^{n})),\) éq 4.1-3

où \(\omega\) est un paramètre de relaxation, constant sur tout le domaine et durant toute le processus itératif, représentant le terme \(\frac{1}{\tau '({u}^{n})}\) .

En raison de la non-convexité de la fonction \(\tau (u)\) , ce paramètre de relaxation doit nécessairement vérifier la condition suivante [bib2], [bib3]:

\(\omega \le \frac{1}{\underset{n}{\max}\tau '({u}^{n})}\) éq 4.1-4

Dans la pratique on prend \(\omega =\frac{1}{\underset{n}{\max}\tau '({u}^{n})}\) .

En prenant en compte l’approximation [éq 4.1-3], la discrétisation du deuxième terme de l’équation [éq3-2] est exprimée de la façon suivante:

\(\begin{array}{cc}\underset{\Omega}{\int}V.\text{grad}{u}^{n+1}\psi d\Omega & =\underset{\Omega}{\int}V.\text{grad}{u}^{n}\psi d\Omega +\underset{\Omega}{\int}\omega V.\text{grad}{T}^{n+1}\psi d\Omega -\underset{\Omega}{\int}\omega V.\text{grad}\tau ({u}^{n})\psi d\Omega ,\end{array}\) éq 4.1-5

Traitement des non linéarités liées à la condition de Fourier non linéaire et à la conductivité thermique#

La non linéarité liée à la condition de flux normal non linéaire est traitée en considérant le développement au premier ordre de la fonction (supposée suffisamment régulière) \(\alpha (T)\) qui est donné par:

\(\alpha ({T}^{n+1})=\alpha ({T}^{n})+\alpha '({T}^{n})({T}^{n+1}-{T}^{n}),\) éq 4.2-1

où \((.)'\) désigne la dérivée de la fonction \((.)\) par rapport à son argument.

Il est apparu nécessaire de décider d’une stratégie de discretisation du terme \(k(T)\mathit{grad}T\) dans l’équation [éq 3-2] afin de pouvoir traiter cette non linéarité pour le problème stationnaire que nous considérons. Pour cela, nous avons adopté l’approximation suivante:

\(\begin{array}{cc}k({T}^{n+1})\text{grad}{T}^{n+1}& =k({T}^{n})\text{grad}{T}^{n+1}-\left[k({T}^{n})-k({T}^{n-1})\right]\text{grad}{T}^{n}\end{array}\) éq 4.2-2

Cette discrétisation est en fait une simplification du développement au premier ordre du terme \(k(T)\text{grad}T\) . Elle se trouve être efficace notamment en raison de la faible non linéarité de la fonction \(k(T)\) dans la pratique.

Remarque:

Notons également que l’approximation purement explicite suivante:

\(k({T}^{n+1})\text{grad}{T}^{n+1}\approx k({T}^{n})\text{grad}{T}^{n+1},\)

donne également des résultats satisfaisants. Cette observation a été vérifiée à partir de plusieurs expériences numériques.

Algorithme implanté dans Code_Aster#

Le schéma numérique employé pour la résolution du problème de référence [éq 2-1], [éq 2-5] est déduit de la formulation variationnelle [éq 3-2] et du traitement des différentes non linéarités, [éq 4.1-5], [éq 4.2-1], [éq 4.2-2], discutées dans la section précédente.

L’algorithme de résolution est constitué par l’enchaînement de deux opérations successives à chaque itération du calcul.

Connaissant les champs solutions à l’itération \(n\) : \({T}^{n}\) aux nœuds et \({u}^{n}\) aux points de Gauss, on cherche les solutions \({T}^{n+1}\text{et}{u}^{n+1}\) à l’itération \(n+1\) comme suit:

\(\forall \psi \in V,\)

\(\begin{array}{}\begin{array}{cc}\underset{\Omega}{\int}k({T}^{n})\text{grad}{T}^{n+1}.\text{grad}\psi d\Omega & +\underset{\Omega}{\int}\omega V.\text{grad}{T}^{n+1}\psi d\Omega \\ & +\underset{{\partial}_{1}\Omega }{\int}\gamma {T}^{n+1}\psi d\Gamma -\underset{{\partial}_{3}\Omega }{\int}\alpha '({T}^{n}){T}^{n+1}\psi d\Gamma \end{array}\\ =\underset{\Omega}{\int}Q\psi d\Omega +\underset{{\partial}_{1}\Omega }{\int}\gamma {T}_{\text{ext}}\psi d\Gamma +\underset{{\partial}_{2}\Omega }{\int}{q}_{0}\psi d\Omega \\ +\underset{{\partial}_{3}\Omega }{\int}(\alpha ({T}^{n})-\alpha '({T}^{n}){T}^{n})\psi d\Gamma +\underset{\Omega}{\int}\left[k({T}^{n})-k({T}^{n-1})\right]\text{grad}{T}^{n}.\text{grad}\psi d\Omega \\ +\underset{\Omega}{\int}\omega V.\text{grad}t({u}^{n})\psi d\Omega -\underset{\Omega}{\int}V.\text{grad}{u}^{n}\psi d\Omega ,\end{array}\) éq 5-1

\({u}^{n+1}={u}^{n}+\omega ({T}^{n+1}-\tau ({u}^{n}))\) éq 5-2

A chaque itération, un problème linéaire de convection-diffusion est résolu pour obtenir le champ aux nœuds \({T}^{n+1}\) [éq 5-1], et ensuite une simple correction locale est effectuée pour obtenir le champ aux points de Gauss \({u}^{n+1}\) [éq 5-2].

Le critère d’arrêt adopté dans Code_Aster fait intervenir en même temps les deux champs solutions : le champ de température, et le champ enthalpie.

L’algorithme continue les itérations tant que au moins une des variations relatives des itérés est supérieure à la tolérance correspondante donnée par l’utilisateur:

\(\begin{array}{}\frac{{(\sum_{i=1,...,\text{nddl}}{({T}_{i}^{n+1}-{T}_{i}^{n})}^{2})}^{1/2}}{{(\sum_{i=1,...,\text{nddl}}{({T}_{i}^{n+1})}^{2})}^{1/2}}>\text{tole 1}\\ \\ \frac{{(\sum_{i=1,...,\text{npg}}{({u}_{i}^{n+1}-{T}_{i}^{n})}^{2})}^{1/2}}{{(\sum_{i=1,...,\text{npg}}{({u}_{i}^{n+1})}^{2})}^{1/2}}>\text{tole 2}\end{array}\)

où \(\mathit{nddl}\) est le nombre total des degrés de libertés aux nœuds, et \(\mathit{npg}\) est le nombre total des points de Gauss.

\(\mathit{tole}1\) est renseigné sous le mot clé crit_temp_rela du mot clé facteur convergence de l’opérateur ther_non_line_mo.

\(\mathit{tole}2\) est renseigné sous le mot clé crit_enth_rela du mot clé facteur convergence de l’opérateur ther_non_line_mo.

Principales options de calcul dans Code_Aster#

On présente ci-dessous les principales options de Code_Aster spécifiques au déroulement de l’algorithme [éq 5-1], [éq 5-2] ci-dessus. En revanche, nous ne mentionnerons pas les options non spécifiques de Code_Aster et qui sont utilisées dans le calcul :

  • Conditions aux limites :

Fourier Linéaire

RIGI_THER_COET_R

RIGI_THER_COET_F \(\underset{{\partial}_{1}\Omega }{\int}\gamma {T}^{n+1}\psi d\Gamma\)

Fourier Non Linéaire

RIGI_THER_FLUTNL \(\underset{{\partial}_{3}\Omega }{\int}\alpha '({T}^{n}){T}^{n+1}\psi d\Gamma\)

CHAR_THER_FLUTNL \(\underset{{\partial}_{3}\Omega }{\int}(\alpha ({T}^{n})-\alpha '({T}^{n}){T}^{n})\psi d\Gamma\)

  • Matrices élémentaires et second membre :

RIGI_THER_TRANS

\(\underset{\Omega}{\int}k({T}^{n})\text{grad}{T}^{n+1}.\text{grad}\psi d\Omega\)

RIGI_THER_CONV_T

\(\underset{\Omega}{\int}\omega V.\text{grad}{T}^{n+1}\psi d\Omega\)

CHAR_THER_TNL

\(\begin{array}{}\underset{\Omega}{\int}\left[k({T}^{n})-k({T}^{n-1})\right]\text{grad}{T}^{n}.\text{grad}\psi d\Omega \\ +\underset{\Omega}{\int}\omega \text{V}.\text{grad}\tau ({u}^{n})\psi d\Omega -\underset{\Omega}{\int}V.\text{grad}{u}^{n}\psi d\Omega \end{array}\)

Bibliographie#

  1. Duvaut, G., Mécanique des milieux continus, Masson, 1990.

  2. Nochetto, R. H., Numerical methods for free boundary problems, Free Boundary Problems: Theory ans Applications, éds. K. H. Hofmann et J.Spreckels, (Pitman, Boston, 1988).

  3. Paolini, M., Sacchi, G. et Verdi, C., Finite Element approximations of singular parabolic problems, Int. J. Numer., Meth. Engng., 26 , pp.1989-2007, 1988.

  4. Friedman, A., Variational Principles and Free Boundary Problems, Pure and Applied Mathematics (Wiley-Interscience, New york, 1982).

  5. Tamma, K. K. et Namburu, R. R., Recent Advences, Trends and New Perspectives via enthalpy-based finite element formulations for applications to solidification problems, Int. J.Numer., Meth. Engng., 30 , pp.803-820, 1990.

Description des versions du document#

Version Aster

Auteur(s) Organisme(s)

Description des modifications

04/01/00

F. WAECKEL, B. NEDJAR (EDF/IMA/MMN, ENPC)

Texte initial