r5.02.02 Thermique non linéaire#

Résumé

L’opérateur THER_NON_LINE [U4.54.02] permet de résoudre les problèmes de thermique transitoire dans les solides en présence de non-linéarités des propriétés des matériaux (capacité calorifique et conductivité), ou des conditions aux limites (échange thermique de type rayonnement). On présente ici la formulation et l’algorithme employé, ce dernier étant proche de celui lié à l’opérateur STAT_NON_LINE [R5.03.01]. Les différentes options de calcul nécessaires ont été présentées dans les éléments de structure plans, axisymétriques et tridimensionnels [U3.22.01], [U3.23.01] et [U3.24.01].

Conditions aux limites, chargement et condition initiale#

On se reportera à [R5.02.01] pour les conditions aux limites thermiques et les chargements conduisant à des équations linéaires en température ainsi que pour la condition initiale.

Flux normal non-linéaire#

Ce sont des conditions de type Neumann, définissant le flux entrant dans le domaine.

\(-q(x,t).n=g(x,T)\) sur la frontière :math:`Gamma ` éq 2.1-1

où \(g(x,T)\) est une fonction de la température et éventuellement de la variable d’espace \(x\) et/ou de temps \(t\) et \(n\) désigne la normale extérieure à la frontière \(\Gamma\) , \(q\) est le vecteur flux de chaleur (dirigé suivant les températures décroissantes).

Cette expression permet d’introduire par exemple des conditions du type échange avec un coefficient d’échange convectif dépendant de la température:

\(-q(x,t).n=g(x,T)=h(x,T)({T}_{\text{ext}}(x,t)-T)\) éq 2.1-2

Flux normal non-linéaire - condition de type rayonnement à l’infini#

Un cas particulier des conditions aux limites précédentes est le rayonnement à l’infini de corps gris qui se traduit par un cas particulier de fonction \(g(x,T)\) :

\(-q(x,t).n=\sigma \epsilon \left[(T(x)+273.15{)}^{4}-({T}_{\infty}+273.15{)}^{4}\right]\) éq 2.2-1

Les caractéristiques à définir lors de la définition de ce chargement sont l’émissivité , la constante de Stefan-Boltzmann \(\sigma =5,73{.10}^{-8}\mathrm{usi}\) et la température à l’infini.

\(T(r)\) et \({T}_{\infty}\) sont alors exprimées en degrés Celsius. \(–273.15°C\) est la température du zéro absolu.

Formulation variationnelle du problème#

Nous nous bornerons ici à présenter le problème avec uniquement les conditions aux limites de température imposée [R5.02.01 §2.1], de flux normal imposé [R5.02.01 §2.3], d’échange [R5.02.01§2.4], de flux non linéaire [§2.1] et de rayonnement [§2.2].

Soit \(\Omega\) un ouvert de \({R}^{3}\) , de frontière \(\Gamma ={\Gamma}_{1}\cup {\Gamma}_{2}\cup {\Gamma}_{3}\cup {\Gamma}_{4}\cup {\Gamma}_{5}\) .

On doit résoudre l’équation [éq 1.1-4] en \(T\) sur \(\Omega \times ]0,t[\) avec les conditions aux limites:

\(\lbrace \begin{array}{ccccc}T={T}^{d}(r,t)& & & \text{sur}& {\Gamma}_{1}\\ \lambda (T)\frac{\partial T}{\partial n}& =& f(r,t)& \text{sur}& {\Gamma}_{2}\\ \lambda (T)\frac{\partial T}{\partial n}& =& h(r,t)({T}_{\text{ext}}(r,t)-T)& \text{sur}& {\Gamma}_{3}\\ \lambda (T)\frac{\partial T}{\partial n}& =& g(r,T)& \text{sur}& {\Gamma}_{4}\\ \lambda (T)\frac{\partial T}{\partial n}& =& \sigma \epsilon \left[(T+273.15{)}^{4}-({T}_{\infty}+273.15{)}^{4}\right]& \text{sur}& {\Gamma}_{5}\end{array}\) éq 3-1

et avec, éventuellement, des conditions initiales \(T(t=0)\) . Si ces dernières ne sont pas précisées, on résoud au préalable le problème stationnaire, c’est à dire l’équation [éq 1.3-1] sans le terme d’évolution temporelle.

Soit \(v\) une fonction suffisamment régulière s’annulant sur \({\Gamma}_{1}\) , en remarquant:

\(\begin{array}{}\frac{d}{\text{dt}}(\underset{\Omega}{\int}\beta (T).\mathrm{v.d}\Omega )=\underset{\Omega}{\int}\dot{\beta}(T).\mathrm{v.d}\Omega \\ \begin{array}{cc}\underset{\Omega}{\int}\lambda (T)\nabla T.\nabla \mathrm{v.d}\Omega =& -\underset{\Omega}{\int}div(\lambda (T)\nabla T).\mathrm{v.d}\Omega +\underset{\Gamma}{\int}\lambda (T)\frac{\partial T}{\partial n}.\mathrm{v.d}\Gamma \\ & \end{array}\end{array}\) éq 3-2

la formulation faible de l’équation de la chaleur peut alors s’écrire :

\(\frac{d}{\text{dt}}(\underset{\Omega}{\int}\beta (T).\mathrm{v.d}\Omega )+\underset{\Omega}{\int}\lambda (T)\nabla T.\nabla \mathrm{vd}\Omega -\underset{\Gamma}{\int}\lambda (T)\frac{\partial T}{\partial n}.\mathrm{v.d}\Gamma =\underset{\Omega}{\int}{r}_{\text{vol}}.\mathrm{v.d}\Omega\) éq 3-3

On en déduit la formulation variationnelle du problème:

\(\begin{array}{c}\underset{\Omega}{\int}\frac{d\beta (T)}{\text{dt}}.\mathit{v.d}\Omega +\underset{\Omega}{\int}\lambda (T)\nabla \mathit{T.}\nabla \mathit{v.d}\Omega +\underset{{\Gamma}_{3}}{\int}\mathit{hT.v.d}{\Gamma}_{3}=\\ \underset{\Omega}{\int}{r}_{\text{vol}}.\mathit{v.d}\Omega +\underset{{\Gamma}_{2}}{\int}\mathit{f.v.d}{\Gamma}_{2}+\underset{{\Gamma}_{3}}{\int}{\mathit{hT}}_{\text{ext}}.\mathit{vd}{\Gamma}_{3}+\\ \underset{{\Gamma}_{4}}{\int}\mathit{g.v.d}{\Gamma}_{4}+\underset{{\Gamma}_{5}}{\int}\sigma \epsilon .\left[(T+273.15{)}^{4}-({T}_{\infty}+273.15{)}^{4}\right].\mathit{v.d}{\Gamma}_{5}\end{array}\) éq 3-4

Discrétisation en temps de l’équation différentielle#

Introduction de la thêta-méthode#

Une façon classique de discrétiser une équation différentielle du premier ordre est la :math:`Theta ` -méthode. Considérons l’équation différentielle suivante :

\(\lbrace \begin{array}{c}\dot{y}(t)=\phi (t,y(t))\\ y(0)={y}_{0}\end{array}\) éq 4.1-1

La \(\Theta\) ‑méthode consiste à discrétiser l’équation [éq 4.1-1] par un schéma aux différences finies

\(\frac{1}{\Delta t}({y}_{n+1}-{y}_{n})=\theta .\varphi ({t}_{n+1},{y}_{n+1})+(1-\theta ).\varphi ({t}_{n},{y}_{n})\) éq 4.1-2

où \({y}_{n+1}\) est une approximation de \(y({t}_{n+1})\) , \({y}_{n}\) étant supposée connue et \(\theta\) est le paramètre de la méthode, \(\theta \in \left[0,1\right]\) .

Remarque :

si \(\theta =0\) le schéma est dit explicite,

si \(\theta \ne 0\) le schéma est dit implicite.

Application à l’équation de la chaleur#

Utilisons la \(\Theta\) -méthode dans la formulation variationnelle de l’équation de la chaleur, où l’on a posé:

\(\begin{array}{cccc}{T}^{+}=T(r,t+\Delta t)& {T}^{-}=T(r,t)& {h}^{+}=h(r,t+\Delta t)& {h}^{-}=h(r,t)\\ {f}^{+}=f(r,t+\Delta t)& {f}^{-}=f(r,t)& {T}_{\text{ext}}^{+}={T}_{\text{ext}}(r,t+\Delta t)& {T}_{\text{ext}}^{-}={T}_{\text{ext}}(r,t)\\ {r}_{\text{vol}}^{+}={r}_{\text{vol}}(r,t+\Delta t)& {r}_{\text{vol}}^{-}={r}_{\text{vol}}(r,t)& {T}^{d+}={T}^{d}(r,t+\Delta t)& {T}^{d-}={T}^{d}(r,t)\\ {g}^{+}=g(r,{T}^{+})& {g}^{-}=g(r,{T}^{-})& & \end{array}\)

où \({T}^{d}(r,t)\) représente la température imposée sur la frontière du domaine, en fonction du temps et de l’espace.

Introduisons les espaces de fonctions suivants:

\(\begin{array}{c}{V}_{{t}^{+}}=\left\lbrace v\in {H}^{1}(\Omega ){v}_{\text{/}{\Gamma}_{1}}={T}^{d+}\right\rbrace \\ {V}_{{t}^{-}}=\left\lbrace v\in {H}^{1}(\Omega ){v}_{\text{/}{\Gamma}_{1}}={T}^{d-}\right\rbrace \\ {V}_{0}=\left\lbrace v\in {H}^{1}(\Omega ){v}_{\text{/}{\Gamma}_{1}}=0\right\rbrace \end{array}\)

Le champ \({T}^{-}\in {V}_{{t}^{-}}\) étant supposé connu, on cherche \({T}^{+}\in {V}_{{t}^{+}}\) tel que \(\forall v\in {V}_{0}\) :

\(\begin{array}{cc}\underset{\Omega}{\int}\frac{\beta ({T}^{+})-\beta ({T}^{-})}{\Delta t}\mathit{v.d}\Omega & +\underset{\Omega}{\int}(\theta \lambda ({T}^{+})\nabla {T}^{+}.\nabla v+(1-\theta )\lambda ({T}^{-})\nabla {T}^{-}.\nabla v)d\Omega \\ -\underset{{\Gamma}_{2}}{\int}(\theta {f}^{+}+(1-\theta ){f}^{-})\mathit{v.d}{\Gamma}_{2}& -\underset{{\Gamma}_{3}}{\int}(\theta {h}^{+}{T}_{\text{ext}}^{+}+(1-\theta ){h}^{-}{T}_{\text{ext}}^{-}-\theta {h}^{+}{T}^{+}-(1-\theta ){h}^{-}{T}^{-})\mathit{v.d}{\Gamma}_{3}\\ -\underset{{\Gamma}_{4}}{\int}(\theta {g}^{+}+(1-\theta ){g}^{-})\mathit{v.d}{\Gamma}_{4}& =\\ \underset{\Omega}{\int}(\theta {r}_{\text{vol}}^{+}+(1-\theta ){r}_{\text{vol}}^{-})\mathit{v.d}\Omega & +\underset{\Omega}{\int}(\theta {r}_{\text{v}}({T}^{+})+(1-\theta ){r}_{\text{v}}({T}^{-}))\mathit{v.d}\Omega \end{array}\) éq 4.2-1

Pour ne pas alourdir excessivement l’écriture et dans la mesure où le procédé est identique aux autres termes, on n’a pas fait figurer le terme de rayonnement dans ces équations (intégrale sur \({\Gamma}_{5}\) ).

En posant :

\(\begin{array}{c}({\text{hT}}_{\text{ext}}{)}^{\theta}=\theta {h}^{+}{T}_{\text{ext}}^{+}+(1-\theta ){h}^{-}{T}_{\text{ext}}^{-}\\ {f}^{\theta}=\theta {f}^{+}+(1-\theta ){f}^{-}\\ {r}^{\theta}=\theta {r}_{\text{vol}}^{+}+(1-\theta ){r}_{\text{vol}}^{-}\end{array}\) on obtient finalement :

\(\begin{array}{c}\underset{\Omega}{\int}\frac{\beta ({T}^{+})}{\Delta t}\mathit{v.d}\Omega +\theta \underset{\Omega}{\int}\lambda ({T}^{+})\nabla {T}^{+}.\nabla \mathit{v.}d\Omega +\theta \underset{{\Gamma}_{3}}{\int}{h}^{+}{T}^{+}\mathit{v.d}{\Gamma}_{3}\\ -\theta \underset{{\Gamma}_{4}}{\int}g({T}^{+}).\mathit{v.d}{\Gamma}_{4}-\theta \underset{\Omega}{\int}{r}_{v}({T}^{+}).\mathit{v.d}\Omega ={L}_{1}(v,{T}^{-})\\ \forall v\in {V}_{0}\end{array}\) éq 4.2-2

où on a posé :

\(\begin{array}{c}{L}_{1}(v,{T}^{-})=\underset{\Omega}{\int}\frac{\beta ({T}^{-})}{\Delta t}\mathit{v.d}\Omega -\underset{\Omega}{\int}(1-\theta )\lambda ({T}^{-})\nabla {T}^{-}.\nabla v.d\Omega +\underset{{\Gamma}_{2}}{\int}{f}^{\theta}\mathit{v.d}{\Gamma}_{2}\\ +\underset{{\Gamma}_{3}}{\int}(({\text{hT}}_{\text{ext}}{)}^{\theta}-(1-\theta ){h}^{-}{T}^{-})\mathit{v.d}{\Gamma}_{3}+\underset{\Omega}{\int}{r}^{\theta}\mathit{v.d}\Omega \\ +(1-\theta )\underset{{\Gamma}_{4}}{\int}g({T}^{-})\mathit{v.d}{\Gamma}_{4}+(1-\theta )\underset{\Omega}{\int}{r}_{v}({T}^{-})\mathit{v.d}\Omega \end{array}\) éq 4.2-3

A un instant de calcul donné, ce terme est connu. En effet, seule la température à l’instant précédent, \({T}^{-}\) , ainsi que les valeurs à l’instant courant de fonction connues du temps, interviennent.

Dans le cas où la répartition de température dans le système à l’instant initial n’est pas fournie, on résout le problème stationnaire. Les termes d’évolution disparaissent, \(\Theta =1\) ; le champ de température à l’instant initial est donné par:

\(\begin{array}{c}\underset{\Omega}{\int}\lambda ({T}^{t=0})\nabla {T}^{t=0}.\nabla \mathit{v.d}\Omega +\underset{{\Gamma}_{3}}{\int}{h}^{t=0}{T}^{t=0}\mathit{v.d}{\Gamma}_{3}-\underset{{\Gamma}_{4}}{\int}g({T}^{t=0})\mathit{v.d}{\Gamma}_{4}\\ =\underset{{\Gamma}_{2}}{\int}{f}^{t=0}\mathit{v.d}{\Gamma}_{2}+\underset{{\Gamma}_{3}}{\int}{h}^{t=0}{T}_{\text{ext}}^{t=0}\mathit{v.d}{\Gamma}_{3}+\underset{\Omega}{\int}{r}^{t=0}\mathit{v.d}\Omega \\ \forall v\in {V}_{0}\end{array}\) éq 4.2-4

Le problème s’écrit finalement sous la forme condensée:

\(\lbrace \begin{array}{c}\mathrm{Soit}{T}^{-}\in {V}_{{t}^{-}}\mathrm{connu},\mathrm{trouver}{T}^{+}\in {V}_{{t}^{+}}\mathrm{tel}\mathrm{que}\\ \forall v\in {V}_{0}:a(v,{T}^{+})={L}_{1}(v,{T}^{-})\end{array}\) éq 4.2-5

Discrétisation spatiale et adaptation de l’algorithme de Newton au problème#

Le principe de la méthode de Newton est très détaillé dans [R5.03.01], on n’exposera ici que les adaptations propres à l’algorithme de thermique non linéaire.

Discrétisation spatiale#

Soit \({P}_{h}\) un découpage de l’espace \(\Omega ` , désignons par :math:`N\) le nombre de nœuds du maillage, \({p}_{i}\) la fonction de forme associé au nœud \(i\) . On désigne par \(J\) l’ensemble des nœuds appartenant à la frontière \({G}_{1}\) .

Soient :

\(\begin{array}{c}{V}_{{t}^{+}}^{h}=\left\lbrace v=\sum_{i=1,N}{v}_{i}{p}_{i}(x);{v}_{j}={T}^{d}({x}_{j},{t}^{+})\text{}j\in J\right\rbrace \\ {V}_{{t}^{-}}^{h}=\left\lbrace v=\sum_{i=1,N}{v}_{i}{p}_{i}(x);{v}_{j}={T}^{d}({x}_{j},{t}^{-})\text{}j\in J\right\rbrace \\ {V}_{0}^{h}=\left\lbrace v=\sum_{i=1,N}{v}_{i}{p}_{i}(x);{v}_{j}=0\text{}j\in J\right\rbrace \end{array}\) éq 5.1-1

Le problème [éq 4.2-5] peut être remplacé par le problème discrétisé à nombre d’inconnues fini suivant:

Soit \({T}^{-}\in {V}_{t}^{h}\) connu, trouver \({T}^{+}\in {V}_{{t}^{+}}^{h}\) tel que

\({v}_{h}\in {V}_{0}^{h}a({v}_{h},{T}^{+})={L}_{1}({v}_{h},{T}^{-})\)

éq 5.1-2

qu’on peut aussi écrire, avec le même formalisme que STAT_NON_LINE [R5.03.01], sous forme vectorielle:

\({v}^{T}R({T}^{+},{t}^{+})={v}^{T}L({T}^{-},{t}^{+})\)	\(\forall v\) tel que \(Bv=0\)
\({\mathit{BT}}^{+}={T}^{d}({t}^{+})\)

éq 5.1-3

où l’opérateur \(B\) exprime la condition aux limites de température imposée \({T}^{+}\in {V}_{{t}^{+}}^{h}\) . Il est défini par:

\((\text{Bv}{)}_{j}=\lbrace \begin{array}{c}0\mathrm{si}j\notin J\\ {v}_{j}\mathrm{si}j\in J\end{array}\) éq 5.1-4

Le cas où l’application \(R\) est linéaire est traité par la commande THER_LINEAIRE [R5.02.01].

La dualisation des conditions aux limites, détaillée dans [R3.03.01], conduit au problème non linéaire en \({T}^{+}\) :

\(\lbrace \begin{array}{c}R({T}^{+},{t}^{+})+{B}^{T}{\lambda}^{+}=L({T}^{-},{t}^{+})\\ {\mathit{BT}}^{+}={T}^{d}({t}^{+})\end{array}\) éq 5.1-5

Les inconnues sont le couple \(({T}^{+},{\lambda}^{+})\) , où \({\lambda}^{+}\) représente les «multiplicateurs de Lagrange» des conditions aux limites de Dirichlet.

Résoudre le système [éq 5.1-5] revient à annuler en \(({T}_{i}^{+},{\lambda}_{i}^{+})\) le vecteur \(F({T}^{+},{\lambda}^{+})\) , appelé résidu, défini par:

\(F({T}^{+},{\lambda}^{+})=(\begin{array}{c}L({T}^{-},{t}^{+})-R({T}^{+},{t}^{+})-{B}^{T}{\lambda}^{+}\\ {T}^{d}({t}^{+})-{\mathit{BT}}^{+}\end{array})\) éq 5.1-6

La méthode de Newton consiste à construire une suite de vecteurs \({\left\lbrace {x}^{n}\right\rbrace }_{n}\) convergeant vers la solution de \(F(x)=0\) à l’aide de l’application linéaire tangente de \(F\) .

Calcul stationnaire#

Le problème variationnel est celui de l’équation [éq 4.2-4]. A noter: en calcul stationnaire, l’enthalpie n’intervient pas dans l’application \(R\) .

On introduit la matrice de l’application linéaire tangente de la fonction \(R({T}^{n})\) :

\({K}^{n}=\frac{\partial R}{\partial T}{\mid}_{{T}_{n}}\)

Celle de la fonction \(F({T}^{n},{\lambda}^{n})\) est alors:

\(\left[\begin{array}{cc}{K}^{n}& {B}^{T}\\ B& 0\end{array}\right]\)

Dans le cas du calcul stationnaire, on doit itérer à partir d’une valeur d’initialisation du champ de température qui est donnée par le mot-clef ETAT_INIT. La première itération du calcul, dite itération de prédiction, consiste à résoudre le système suivant:

\(\left[\begin{array}{cc}K({T}_{0})& {B}^{T}\\ B& 0\end{array}\right]\left[\begin{array}{c}{T}_{1}-{T}_{0}\\ {\lambda}_{1}-{\lambda}_{0}\end{array}\right]=\left[\begin{array}{c}L-R({T}_{0})-{B}^{T}{\lambda}_{0}\\ {T}^{d}-{\mathrm{BT}}_{0}\end{array}\right]\) éq 5.2-1

Comme on peut le voir dans l’équation du problème stationnaire [éq 4.2-4], la température n’apparaît pas au second membre: on écrit \(L\) et non \(L({T}_{0})\) .

Si le problème est linéaire, \(R({T}_{0})=K({T}_{0}){T}_{0}=K.{T}_{0}\) . Tous les termes en \({T}_{0}\) disparaissent par simplification. La solution est obtenue en une itération par inversion d’un système identique à celui décrit dans [R5.02.01 §6].

Les itérations suivantes sont des itérations de Newton, avec réactualisation ou non de la matrice tangente \(K\) .

\(\left[\begin{array}{cc}K({T}_{(i)})& {B}^{T}\\ B& 0\end{array}\right]\left[\begin{array}{c}{T}_{i+1}-{T}_{i}\\ {\lambda}_{i+1}-{\lambda}_{i}\end{array}\right]=\left[\begin{array}{c}L-R({T}_{i})-{B}^{T}{\lambda}_{i}\\ 0\end{array}\right]\) éq 5.2-2

Pour l’itération de prédiction, l’écriture du sous-système inférieur de l’équation [éq 5.2-1], après simplification, nous assure que \({\mathrm{BT}}_{1}={T}^{d}\) . Le premier itéré et tous les suivants vérifient donc les conditions de Dirichlet.

Les parenthèses autour de l’indice d’itération dans l’expression \(K({T}_{(i)})\) signifient qu’on peut réactualiser ou non la matrice tangente au fil des itérations.

Remarque:

La température d’initialisation \({T}_{0}\) n’a d’influence que pour un calcul stationnaire non linéaire. En étant de l’ordre de grandeur des températures attendues, elle permettrait de «partir» de moins loin de la solution qu’un champ nul partout; et ainsi diminuerait le nombre d’itérations.

Calcul transitoire#

Pour la première itération du pas de temps, dite itération de prédiction, on «fait comme si» le problème décrit dans [éq 5.1-5] était linéaire. Cette formulation doit permettre d’obtenir directement la solution pour un problème de thermique linéaire. Mais ici, la situation est un peu différente du calcul stationnaire à cause de la formulation en enthalpie. La linéarisation de [éq 5.1-5] donne:

\(\lbrace \begin{array}{c}R({T}^{-},{t}^{+})+K({T}^{-},{t}^{+})({T}^{+}-{T}^{-})+{B}^{T}{\lambda}^{+}=L({T}^{-},{t}^{+})\\ {\mathit{BT}}^{+}={T}^{d}({t}^{+})\end{array}\) éq 5.3-1

Ce qui revient à résoudre, pour le problème présenté sous forme matricielle:

\(\left[\begin{array}{cc}K({T}^{-})& {B}^{T}\\ B& 0\end{array}\right]\left[\begin{array}{c}{T}_{1}^{+}\\ {\lambda}_{1}^{+}\end{array}\right]=\left[\begin{array}{c}L({T}^{-},{t}^{+})+K({T}^{-}){T}^{-}-R({T}^{-})\\ {T}^{d}({t}^{+})\end{array}\right]\) éq 5.3-2

La fonction enthalpie est connue à une constante d’intégration près qui apparaît dans la relation liant \(R({T}^{-})\) à \(K({T}^{-}){T}^{-}\) . Cette même constante se retrouve dans l’expression de \(L({T}^{-}{\text{,t}}^{+})\) . On peut alors l’éliminer en aboutissant au système d’équations suivant:

\(\left[\begin{array}{cc}K({T}^{-})& {B}^{T}\\ B& 0\end{array}\right]\left[\begin{array}{c}{T}_{1}^{+}\\ {\lambda}_{1}^{+}\end{array}\right]=\left[\begin{array}{c}\tilde{L}({T}^{-},{t}^{+})\\ {T}^{d}({t}^{+})\end{array}\right]\) éq 5.3-3

où \(\tilde{L}({T}^{-}{\text{,t}}^{+})\) est le second membre calculé avec la capacité calorifique et non l’enthalpie (option CHAR_THER_EVOLNI [§6.2]).

Enfin, comme pour le cas stationnaire vu au chapitre précédent, les itérations suivantes sont des itérations de Newton:

\(\left[\begin{array}{cc}K({T}_{(i)},{t}^{+})& {B}^{T}\\ B& 0\end{array}\right]\left[\begin{array}{c}{T}_{i+1}^{+}-{T}_{i}^{+}\\ {\lambda}_{i+1}^{+}-{\lambda}_{i}^{+}\end{array}\right]=\left[\begin{array}{c}L({T}^{-},{t}^{+})-R({T}_{i},{t}^{+})-{B}^{T}{\lambda}_{i}\\ 0\end{array}\right]\) éq 5.3-4

Cette fois-ci, par contre, \(L({T}^{-}{\text{,t}}^{+})\) est calculé avec l’enthalpie et non la capacité calorifique pour être cohérent avec \(R({T}_{i}^{+})\) .

Convergence#

Puisque le temps intervient dans l’expression de la matrice tangente, et également le pas de temps, on préfère actualiser systématiquement celle-ci au début de chaque pas pour ne pas trop dégrader la vitesse de convergence. En revanche, liberté est laissée à l’utilisateur de contrôler sa fréquence de calcul au cours d’un pas de temps.

A chaque itération, on peut effectuer la recherche d’un pas de progression optimum vers la solution par quelques itérations (2 ou 3) de recherche linéaire. Cette méthode est décrite en détail dans [R5.03.01].

Le calcul est réputé convergé quand le vecteur résidu est nul [éq 5.1-6]:

\(F({T}_{{i}^{}}^{+},{\lambda}_{{i}^{}}^{+},{t}^{+})=(\begin{array}{c}L({T}^{-},{t}^{+})-R({T}_{{i}^{}}^{+},{t}^{+})-{B}^{T}{\lambda}_{{i}^{}}^{+}\\ {T}^{d}({t}^{+})-{\mathit{BT}}_{{i}^{}}^{+}\end{array})\) éq 5.4-1

La partie inférieure du vecteur est toujours nulle (conditions de Dirichlet). On vérifie donc:

\(\frac{{\mathrm{\parallel }L({T}^{-},{t}^{+})-R({T}_{{i}^{}}^{+},{t}^{+})-{B}^{T}{\lambda}_{{i}^{}}^{+}\mathrm{\parallel }}_{2}}{{\mathrm{\parallel }L({T}^{-},{t}^{+})-{B}^{T}{\lambda}_{{i}^{}}^{+}\mathrm{\parallel }}_{2}}\le \epsilon\) éq 5.4-2

L’utilisateur a également la possibilité d’arrêter les itérations sur un critère absolu:

\({\mathrm{\parallel }L({T}^{-},{t}^{+})-R({T}_{{i}^{}}^{+},{t}^{+})-{B}^{T}{\lambda}_{{i}^{}}^{+}\mathrm{\parallel }}_{\infty}\le \epsilon\) éq 5.4-3

Principales options de thermique non-linéaire calculées dans Code_Aster#

Conditions aux limites et chargements#

On se reportera à [R5.02.01] pour les conditions aux limites et les chargements linéaires.

Flux non linéaire	CHAR_THER_FLUNL	\(\underset{{\Gamma}_{4}}{\int}(1-\theta )g({T}^{-})\mathit{v.d}{\Gamma}_{4}\)
Rayonnement	CHAR_THER_RAYO_R CHAR_THER_RAYO_F	\(\underset{{\Gamma}_{4}}{\int}\sigma \epsilon \left[(T+273.15{)}^{4}-(1-\theta )({T}^{-}+273.15{)}^{4}\right].\mathit{v.d}{\Gamma}_{4}\)
Source non-linéaire	CHAR_THER_SOURNL	\(\underset{\Omega}{\int}(1-\theta ){r}_{v}({T}^{-})\mathit{v.d}\Omega\)

Calcul des matrices élémentaires et terme transitoire#

Conductivité	RIGI_THER_TANG	:math:`underset{Omega}{int}theta lambda ({T}^{+})nabla v.nabla v.dOmega `
Inertie thermique	MASS_THER_TANG	:math:`underset{Omega}{int}frac{rho text{Cp}}{Delta t}v.v.dOmega `
Rayonnement	MTAN_THER_RAYO_R MTAN_THER_RAYO_F	\(\underset{{\Gamma}_{4}}{\int}\theta \mathrm{.4.}\sigma .\epsilon ({T}^{+}+273.15{)}^{3}\mathrm{v.v.d}{\Gamma}_{4}\)
Coefficient d’échange	MTAN_THER_COEF_R MTAN_THER_COEF_F	\(\underset{{\Gamma}_{4}}{\int}\theta \mathrm{h.v.vd}{\Gamma}_{4}\)
Flux non linéaire	MTAN_THER_FLUXNL	\(-\underset{{\Gamma}_{4}}{\int}\theta \frac{\text{dg}}{\text{dT}}({T}^{+})v.\mathit{v.d}{\Gamma}_{4}\)
Source non linéaire	MTAN_THER_SOURNL	\(-\underset{\Omega}{\int}\theta \frac{\text{d}{r}_{v}}{\text{dT}}({T}^{+})v.\mathit{v.d}\Omega\)
Terme transitoire	CHAR_THER_EVOLNI	\(\underset{\Omega}{\int}\frac{\beta ({T}^{-})}{\Delta t}.\mathrm{v.d}\Omega -\underset{\Omega}{\int}(1-\theta )\lambda ({T}^{-})\nabla {T}^{-}.\nabla \mathrm{v.d}\Omega\)
		\(\underset{\Omega}{\int}\frac{\rho {\text{CpT}}^{-}}{\Delta t}.\mathrm{v.d}\Omega -\underset{\Omega}{\int}(1-\theta )\lambda ({T}^{-})\nabla {T}^{-}.\nabla \mathrm{v.d}\Omega\)

Calcul du résidu#

	RAPH_THER	:math:`underset{Omega}{int}theta lambda ({T}^{i})nabla {T}^{i}.nabla v.dOmega `
	MASS_THER_RESI	:math:`underset{Omega}{int}frac{1}{Delta t}beta ({T}^{i})v.dOmega `
Rayonnement	RESI_THER_RAYO_R RESI_THER_RAYO_F	\(\underset{{\Gamma}_{4}}{\int}\theta \sigma \epsilon ({T}^{i}+273.15{)}^{4}\mathrm{v.}d{\Gamma}_{4}\)
Coefficient d’échange	RESI_THER_COEF_R RESI_THER_COEF_F	\(\underset{{\Gamma}_{3}}{\int}(\theta {h}^{+}{T}^{i})\mathrm{v.}d{\Gamma}_{3}\)
Flux non linéaire	RESI_THER_FLUXNL	\(-\underset{{\Gamma}_{3}}{\int}\theta g({T}^{i})\mathit{v.d}{\Gamma}_{3}\)
Source non linéaire	RESI_THER_SOURNL	\(-\underset{\Omega}{\int}\theta {r}_{v}({T}^{i})\mathit{v.d}\Omega\)

Bibliographie#

SALENCON. Mécanique des milieux continus. Ellipses. 1988.
RUUP, PENIGUEL. SYRTHES - Conduction et rayonnement, Manuel théorique de la version 3.1. HE-41/98/048/A