r5.02.01 Algorithme de thermique linéaire transitoire

Résumé:

On présente l’algorithme de thermique transitoire linéaire implanté au sein de la commande THER_LINEAIRE [U4.33.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 [U1.22.01], [U1.23.01] et [U1.24.01].

Conditions aux limites, chargement et condition initiale

On décrit ici uniquement les conditions aux limites thermiques conduisant à des équations linéaires en température, ce qui exclut les conditions de type rayonnement.

Températures imposées

Les conditions de type Dirichlet sont traitées habituellement par dualisation (cf. [R3.03.01]), mais elles peuvent aussi être éliminées dans certains cas (charges cinématiques).

\[T(r,t) = {T}_{1}(r,t) \text{ sur } {\Gamma}_{1}\]

où \({T}_{1}(r,t)\) est une fonction de la variable d’espace et/ou du temps.

Relations linéaires

Ce sont des conditions de type Dirichlet, permettant de définir une relation linéaire entre les valeurs de la température :

• entre deux ou plusieurs noeuds : avec une équation de la forme

\[\sum_{i=1}^{n}{\alpha}_{i}{T}_{i}(r,t) = \beta (t)\]

• entre des couples de noeuds : avec une équation de la forme

\[\sum_{i=1}^{{n}_{1}}{\alpha}_{{1i}}{T}_{i/{\Gamma}_{12}}(r,t) + \sum_{i=1}^{{n}_{2}}{\alpha}_{{2i}}{T}_{i/{\Gamma}_{21}}(r,t) = \beta (t)\]

où \({\Gamma}_{12}\) et \({\Gamma}_{21}\) sont deux sous-parties de la frontière auxquelles on lie deux à deux les valeurs de la température. Ce type de condition aux limites permet de définir des conditions de périodicité.

Flux normal imposé

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

\[-\vector{q}(r,t) \cdot \vector{n} = f(r,t)\text{ sur }{\Gamma}_{2}\]

où \(f(r,t)\) est une fonction de la variable d’espace et/ou de temps et :math:` vector{n}` désigne la normale à la frontière \({\Gamma}_{2}\) .

Échange

Ce sont des conditions de type Neumann modélisant les transferts convectifs sur les bords du domaine.

\[-\vector{q}(r,t) \cdot \vector{n} = h(r,t)({T}_{\text{ext}}(r,t)-T(r,t))\text{ sur }{\Gamma}_{3}\]

où \({T}_{\text{ext}}(r,t)\) est une fonction de la variable d’espace et/ou de temps représentant la température du milieu extérieur, et \(h(r,t)\) est une fonction de la variable d’espace et/ou de temps représentant le coefficient d’échange convectif sur la frontière \({\Gamma}_{3}\) .

Échange paroi

Ce sont des conditions de type Neumann, mettant en jeu deux sous parties de la frontière en vis à vis. Ce type de condition aux limites modélise une résistance thermique d’interface.

\[\lambda \frac{\partial {T}_{1}}{\partial \vector{n}_{1}} = h(r,t)({T}_{2}(r,t)-{T}_{1}(r,t))\text{ sur }{\Gamma}_{12}\]

avec \(\vector{n}_{1}\) normale extérieure à \({\Gamma}_{12}\)

\[\lambda \frac{\partial {T}_{2}}{\partial \vector{n}_{2}} = h(r,t)({T}_{1}(r,t)-{T}_{2}(r,t))\text{ sur }{\Gamma}_{21}\]

avec \(\vector{n}_{2}\) normale extérieure à \({\Gamma}_{21}\) (et donc \(\vector{n}_{1}=-\vector{n}_{2}\) en général).

Source volumique

C’est le terme \(s(r,t)\) fonction de la variable d’espace et/ou de temps.

Condition initiale

C’est l’expression du champ de température à l’instant initial \(t=0\) :

\[T(r,0) = {T}_{0}(r)\]

où \({T}_{0}(r)\) est une fonction de la variable d’espace.

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 [§2.1], de flux normal imposé [§2.3] ou d’échange [§2.4]. Les conditions aux limites d’échange par paroi [§2.5] sont traitées au [§4] et celles avec relations linéaires [§2.2] se ramènent sans difficulté à celles du [§2.1].

Soit \(\Omega\) un ouvert de \({\setR}^{3}\) , de frontière \(\Gamma ={\Gamma}_{1}\cup {\Gamma}_{2}\cup {\Gamma}_{3}\) avec \({\Gamma}_{1}\) la portion de la frontière sur laquelle on impose une température donnée, \({\Gamma}_{2}\) la portion de la frontière sur laquelle on impose un flux normal donné et \({\Gamma}_{3}\) la portion de la frontière sur laquelle on impose une condition d’échange.

La formulation faible de l’équation de la chaleur est :

\[{\int}_{\Omega} { \rho {C}_{p}\frac{\partial T}{\partial t} \, v } \, d\Omega + {\int}_{\Omega} { \lambda \vector{\nabla} T \cdot \vector{\nabla} v} \, d\Omega - {\int}_{\Gamma} { \lambda \frac{\partial T}{\partial \vector{n}} v } \, d\Omega = {\int}_{\Gamma} { s \, v } \, d\Omega\]

où \(v\) est une fonction suffisamment régulière s’annulant uniformément sur \({\Gamma}_{1}\) . Avec les conditions aux limites suivantes :

\[\begin{split}\left \lbrace \begin{array}{lll} T={T}_{1}(r,t)& \text{ sur }{\Gamma}_{1}\\ \lambda \frac{\partial T}{\partial \vector{n}}=f(r,t)& \text{ sur }{\Gamma}_{2}\\ \lambda \frac{\partial T}{\partial \vector{n}}=h(r,t)({T}_{\text{ext}}(r,t)-T)& \text{ sur }{\Gamma}_{3} \end{array} \right .\end{split}\]

La formulation variationnelle du problème est :

\[{\int}_{\Omega} \rho {C}_{p}\frac{\partial T}{\partial t} v \, d\Omega + {\int}_{\Omega} \lambda \vector{\nabla} T \cdot \vector{\nabla} v \, d\Omega + {\int}_{\Gamma_3} h \, T \, v \, d\Gamma = {\int}_{\Omega} s \, v \, d\Omega + {\int}_{\Gamma_2} f \, v \, d\Gamma + {\int}_{\Gamma_3} h \, {T}_{\text{ext}} \, v \, d\Gamma\]

Formulation variationnelle du problème avec condition d’échange entre deux parois

On considère le problème « simplifié » où n’apparaît plus de terme source et où les conditions aux limites sont uniquement du type température imposée et échange paroi.

Soit \(O\) un ouvert de \({ℝ}^{3}\) , de frontière \(\Gamma ={\Gamma}_{1}\cup {\Gamma}_{12}\cup {\Gamma}_{21}\) .

Les conditions aux limites sont dans ce cas :

\[\begin{split}\left \lbrace \begin{array}{ccc}T={T}_{1}(r,t)& & \text{ sur }{\Gamma}_{1}\\ \lambda \frac{\partial {T}_{1}}{\partial {n}_{1}}=h(r,t)({T}_{2}(r,t)-{T}_{1}(r,t))& & \text{ sur }{\Gamma}_{12}\\ \lambda \frac{\partial {T}_{2}}{\partial {n}_{2}}=h(r,t)({T}_{1}(r,t)-{T}_{2}(r,t))& & \text{ sur }{\Gamma}_{21} \end{array} \right .\end{split}\]

En substituant dans la formulation faible de l’équation de la chaleur, on obtient :

\[\underset{\Omega}{\int}\rho {C}_{p}\frac{\partial T}{\partial t}.vd\Omega +\underset{\Omega}{\int}\lambda \nabla T.\nabla vd\Omega +\underset{{\Gamma}_{12}}{\int}h({T}_{/{\Gamma}_{12}}-{T}_{/{\Gamma}_{21}}).vd{\Gamma}_{12}+\underset{{\Gamma}_{21}}{\int}h({T}_{/{\Gamma}_{21}}-{T}_{/{\Gamma}_{12}}).vd{\Gamma}_{21}=0`\]

où \(v\) s’annule uniformément sur \({\Gamma}_{1}\) .

Ce type de conditions aux limites fait apparaître des termes nouveaux mettant en relation des degrés de liberté situés sur les deux frontières en relation.

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

Une façon classique de discrétiser une équation différentielle du premier ordre consiste à utiliser une \(\theta\) -méthode. Considérons l’équation différentielle suivante :

\[\begin{split}\left \lbrace \begin{array}{c}\frac{\partial}{\partial t}y(t)=\varphi (t,y(t))\\ y(0)={y}_{0} \end{array} \right .\end{split}\]

La \(\theta -\) méthode consiste à discrétiser l’équation 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})\]

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 \ge 0\) le schéma est dit implicite.

Supposons \(y\) suffisamment régulière (au moins 3 fois différentiable), par un développement de Taylor au point \({t}_{n}\) on obtient :

\[y({t}_{n+1})-y({t}_{n})=\Delta ty'({t}_{n})+\frac{\Delta {t}^{2}}{2}y\text{''}({t}_{n})+O(\Delta {t}^{2})\]

et

\[\begin{split}\begin{array}{cc}\theta \varphi ({t}_{n+1},y({t}_{n+1}))+(1-\theta )\varphi ({t}_{n},y({t}_{n}))& =\theta y'({t}_{n+1})+(1-\theta )y'({t}_{n})\\ & =y'({t}_{n+1})+\theta (y'({t}_{n+1})-y'({t}_{n}))\\ & =y'({t}_{n})+\theta \Delta ty\text{''}({t}_{n})+O(\Delta {t}^{2})\end{array}\end{split}\]

La solution vérifie donc approximativement :

\[\frac{1}{\Delta t}(y({t}_{n+1})-y({t}_{n}))=\theta \varphi ({t}_{n+1},y({t}_{n+1}))+(1-\theta )\varphi ({t}_{n},y({t}_{n}))+(\frac{1}{2}-\theta )\Delta ty\text{''}({t}_{n})+O(\Delta {t}^{2})\]

Le schéma est d’ordre 1 en temps si \(\theta \ne \frac{1}{2}\) , et d’ordre 2 si \(\theta =\frac{1}{2}\) (schéma de Crank‑Nicolson).

Considérons l’équation différentielle suivante :

\[\begin{split}\left \lbrace \begin{array}{ccc}y'=-\lambda y& t\ge 0& \lambda \in R\\ y(0)={y}_{0}& & \end{array} \right .\end{split}\]

En utilisant la \(\theta\) -méthode dans cette équation différentielle on obtient :

\[{y}_{n+1}=\frac{1-(1-\theta )\lambda \Delta t}{1+\theta \lambda \Delta t}{y}_{n}\text{}0\le n\le N-1\]

Soit encore :

\[{y}_{n+1}={r}^{n}(\lambda \Delta t){y}_{0}\text{avec}r(x)=\frac{1-(1-\theta )x}{1+\theta x}\]

La solution approchée \({y}_{n}\) doit être bornée (la solution exacte du problème initial l’étant), ce qui impose la condition suivante :

\(\mid r(\lambda \Delta t)\mid \le 1\)

En étudiant les variations de la fonction \(r(x)\) , on constate facilement que :

• si \(\theta \ge\) la condition est vérifiée quel que soit \(\Delta t\) , le schéma est inconditionnellement stable;

• si \(\theta <\) la condition n’est vérifiée que si \(\Delta t\le \frac{2}{\lambda (1-2\theta )}\) , le schéma est conditionnellement stable.

Dans la commande THER_LINEAIRE [U4.33.01], le paramètre \(\theta\) est une donnée pouvant être fournie par l’utilisateur, la valeur par défaut est fixée à 0.57. Cette valeur a la réputation d’être préférable à la valeur de Crank‑Nicolson (0,5) et « optimale » pour les interpolations quadratiques, mais nous n’avons pas retrouvé trace des justifications.

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

\[\begin{split}\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)\\ {s}^{+}=s(r,t+\Delta t)& {s}^{-}=s(r,t)& {T}_{1}^{+}={T}_{1}(r,t+\Delta t)& {T}_{1}^{-}={T}_{1}(r,t) \end{array}\end{split}\]

Introduisons les espaces de fonctions suivants :

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

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

\[\begin{split}\begin{array}{}\underset{\Omega}{\int}\rho {C}_{p}\frac{{T}^{+}-{T}^{-}}{\Delta t}vd\Omega +\underset{\Omega}{\int}(\theta \lambda \nabla {T}^{+}.\nabla v+(1-\theta )\lambda \nabla {T}^{-}.\nabla v)d\Omega \\ -\underset{{\Gamma}_{2}}{\int}(\theta {f}^{+}+(1-\theta ){f}^{-})vd{\Gamma}_{2}-\underset{{\Gamma}_{3}}{\int}(\theta {h}^{+}{T}_{\text{ext}}^{+}+(1-\theta ){h}^{-}{T}_{\text{ext}}^{-}-\theta {h}^{+}{T}^{+}-(1-\theta ){h}^{-}{T}^{-})vd{\Gamma}_{3}\\ =\underset{\Omega}{\int}(\theta {s}^{+}+(1-\theta ){s}^{-})vd\Omega \\ \forall v\in {V}_{0}\end{array}\end{split}\]

En posant :

\[\begin{split}\begin{array}{}({\text{hT}}_{\text{ext}}{)}^{\theta}=\theta {h}^{+}{T}_{\text{ext}}^{+}+(1-\theta ){h}^{-}{T}_{\text{ext}}^{-}\\ {f}^{\theta}=\theta {f}^{+}+(1-\theta ){f}^{-}\end{array}\end{split}\]

on obtient finalement :

\[\begin{split}\begin{array}{}\underset{\Omega}{\int}\frac{\rho {C}_{p}}{\Delta t}{T}^{+}vd\Omega +\underset{\Omega}{\int}\theta \lambda \nabla {T}^{+}.\nabla vd\Omega +\underset{{\Gamma}_{3}}{\int}\theta {h}^{+}{T}^{+}vd{\Gamma}_{3}\\ =\underset{\Omega}{\int}\frac{\rho {C}_{p}}{\Delta t}{T}^{-}vd\Omega -\underset{\Omega}{\int}(1-\theta )\lambda \nabla {T}^{-}.\nabla vd\Omega +\underset{{\Gamma}_{2}}{\int}{f}^{\theta}vd{\Gamma}_{2}\\ +\underset{{\Gamma}_{3}}{\int}(({\text{hT}}_{\text{ext}}{)}^{\theta}-(1-\theta ){h}^{-}{T}^{-})vd{\Gamma}_{3}+\underset{\Omega}{\int}(\theta {s}^{+}+(1-\theta ){s}^{-})vd\Omega \\ \forall v\in {V}_{0}\end{array}\end{split}\]

Discrétisation spatiale

Soit \({P}_{h}\) un découpage de l’espace \(\Omega\) , désignons par \(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 \({\Gamma}_{1}\) .

Soient :

\[\begin{split}\begin{array}{}{V}_{{t}^{+}}^{h}=\lbrace v=\sum_{i=1,N}{v}_{i}{p}_{i}(x)\text{};\text{}{v}_{j}={T}_{1}({x}_{j},{t}^{+})\text{}j\in J\rbrace \\ {V}_{{t}^{-}}^{h}=\lbrace v=\sum_{i=1,N}{v}_{i}{p}_{i}(x)\text{};\text{}{v}_{j}={T}_{1}({x}_{j},{t}^{-})\text{}j\in J\rbrace \\ {V}_{0}^{h}=\lbrace v=\sum_{i=1,N}{v}_{i}{p}_{i}(x)\text{};\text{}{v}_{j}=0\text{}j\in J\rbrace \end{array}\end{split}\]

Posons :

\[\begin{split}\begin{array}{}{K}_{ij}{T}_{i}=\underset{{\Omega}_{h}}{\int}\frac{\rho {C}_{p}}{\Delta t}{T}_{i}{p}_{i}{p}_{j}d{\Omega}_{h}+\underset{{\Omega}_{h}}{\int}\theta \lambda {T}_{i}\nabla {p}_{i}.\nabla {p}_{j}d{\Omega}_{h}+\underset{{\Gamma}_{\mathrm{h3}}}{\int}\theta {h}^{+}{T}_{i}{p}_{i}d{\Gamma}_{\mathrm{h3}}\\ \\ \begin{array}{cc}{L}_{j}=& \underset{{\Omega}_{h}}{\int}\frac{\rho {C}_{p}}{\Delta t}{T}^{-}{p}_{j}d{\Omega}_{h}-\underset{{\Omega}_{h}}{\int}(1-\theta )\lambda \nabla {T}^{-}.\nabla {p}_{j}d{\Omega}_{h}+\underset{{\Gamma}_{\mathrm{h2}}}{\int}{f}^{\theta}{p}_{j}d{\Gamma}_{\mathrm{h2}}\\ & +\underset{{\Gamma}_{\mathrm{h3}}}{\int}(({\text{hT}}_{\text{ext}}{)}^{\theta}-(1-\theta ){h}^{-}{T}^{-}){p}_{j}d{\Gamma}_{\mathrm{h3}}+\underset{{\Omega}_{h}}{\int}(\theta {s}^{+}+(1-\theta ){s}^{-}){p}_{j}d{\Omega}_{h}\end{array}\\ \end{array}\end{split}\]

En dualisant les conditions aux limites en température imposée ([R3.03.01]), on fait apparaître l’opérateur \(B\) défini par :

\[\begin{split}(\text{Bv}{)}_{j}= \left \lbrace \begin{array}{ccc}0& \text{si}& j\notin J\\ {v}_{j}& \text{si}& j\notin J\end{array} \right .\end{split}\]

On obtient finalement le système suivant :

\[\begin{split}\left \lbrace \begin{array}{ccccccc}\sum_{i=1}^{N}{K}_{ij}{T}_{i}& +& {({}^{t}\text{}B\lambda )}_{j}& =& {L}_{j}& & \forall j\\ & & {(\mathrm{BT})}_{j}& =& {T}_{1}({x}_{j},t)& & j\in J\end{array} \right .\end{split}\]

Implémentation dans Code_Aster

Equations introduites

La commande THER_LINEAIRE [U4.33.01] permet de traiter l’équation dans le cas transitoire telle qu’elle est décrite ci-dessus, mais elle permet aussi de résoudre le problème stationnaire qui se réduit à l’équation suivante :

\(-div(\lambda \nabla T)=s\text{ dans }\Omega\)

et les conditions aux limites suivantes :

\[\begin{split}\left \lbrace \begin{array}{ccc}T={T}_{1}(r,{t}_{s})& & \text{sur}{\Gamma}_{1}\\ \lambda \frac{\partial T}{\partial n}=q(r,{t}_{s})& & \text{sur}{\Gamma}_{2}\\ \lambda \frac{\partial T}{\partial n}=h(r,t)({T}_{\text{ext}}(r,{t}_{s})-T)& & \text{sur}{\Gamma}_{3}\end{array} \right .\end{split}\]

\({t}_{s}\) étant l’instant pris pour évaluer les conditions aux limites de l’équation.

Dans le cas transitoire, il est nécessaire de fournir un état initial, cet état initial (champ de température) peut être choisi parmi les suivants :

• un champ qui peut être uniforme ou quelconque créé par la commande CREA_CHAMP,

• un champ résultat d’un problème stationnaire décrit par les équations ci-dessus, l’instant de calcul est pris au premier instant défini dans la liste de réels décrivant la discrétisation temporelle définie par l’utilisateur,

• un champ extrait du résultat d’un problème transitoire.

La discrétisation en temps (valeur de \(\Delta t\) ) doit être fournie sous la forme d’une ou plusieurs listes d’instants. Ces listes sont créées par l’utilisateur par la commande DEFI_LIST_REEL [U4.21.04].

Un transitoire thermique peut être calculé en effectuant plusieurs appels à la commande THER_LINEAIRE [U4.33.01] en enrichissant le même concept de type evol_ther en fournissant à partir du deuxième appel l’instant initial de reprise du calcul (pour obtenir \({T}^{-}\) ) et éventuellement l’instant final.

Les champs de températures issus d’un calcul contiennent à la fois la valeur aux nœuds du maillage et aux nœuds de Lagrange. Lors d’une reprise du calcul, il est possible de faire varier le type des conditions aux limites, le champ utilisé pour initier le nouveau calcul en interne est alors réduit aux seuls nœuds du maillage. Le concept résultat de type evol_ther contiendra alors des champs aux nœuds s’appuyant sur des numérotations différentes. Les opérateurs de Code_Aster interpolent ensuite uniquement aux nœuds du maillage lorsque la numérotation diffère.

Principales options thermiques calculées dans Code_Aster

Conditions aux limites et chargements

TEMP_IMPO	DDLI_R DDLI_F	\(\underset{{\Gamma}_{1}}{\int}{T}^{+}{\Phi}^{\ast }d{\Gamma}_{1}+\underset{{\Gamma}_{1}}{\int}{\Phi}^{\theta}{\mathrm{vdG}}_{1}\)
	DDLI_R DDLI_F	\(\underset{{\Gamma}_{1}}{\int}{\Phi}^{\ast }{T}_{1}d{\Gamma}_{1}\)
FLUX_REP	CHAR_THER_FLUN_R CHAR_THER_FLUN_F	\(\underset{{\Gamma}_{2}}{\int}{q}^{\theta}{\mathrm{vdG}}_{2}\)
ECHANGE	CHAR_THER_COEF_R CHAR_THER_COEF_F	\(\underset{{\Gamma}_{3}}{\int}\theta {h}^{+}{T}^{+}vd{\Gamma}_{3}\)
	CHAR_THER_TEXT_R CHAR_THER_TEXT_F	\(\underset{{\Gamma}_{3}}{\int}(({\text{hT}}_{\text{ext}}{)}^{\theta}-(1-\theta ){h}^{-}{T}^{-})vd{\Gamma}_{3}\)
ECHANGE_PAROI	RIGI_THER_PARO_R RIGI_THER_PARO_F	\(\underset{{\Gamma}_{12}}{\int}\theta {h}^{+}({T}_{/{\Gamma}_{12}}^{+}-{T}_{/{\Gamma}_{21}}^{+}){v}_{1}d{\Gamma}_{12}\)
	CHAR_THER_PARO_R CHAR_THER_PARO_F	\(\underset{{\Gamma}_{12}}{\int}(1-\theta ){h}^{-}({T}_{/{\Gamma}_{21}}^{-}-{T}_{/{\Gamma}_{12}}^{-}){v}_{1}d{\Gamma}_{12}\)
SOURCE	CHAR_THER_SOUR_R CHAR_THER_SOUR_F	\(\underset{\Omega}{\int}(\theta {s}^{+}+(1-\theta ){s}^{-})\mathrm{vd}\Omega\)

Calcul des matrices élémentaires et terme transitoire

	RIGI_THER	\(\underset{\Omega}{\int}\theta \lambda \nabla {T}^{+}.\nabla \mathrm{vd}\Omega\)
	MASS_THER	\(\underset{\Omega}{\int}\frac{\rho {C}_{p}}{\Delta t}{T}^{+}\mathrm{vd}\Omega\)
	CHAR_THER_EVOL	\(\underset{\Omega}{\int}\frac{\rho {C}_{p}}{\Delta t}{T}^{-}\mathrm{vd}\Omega -\underset{\Omega}{\int}(1-\theta )\lambda \nabla {T}^{-}.\nabla \mathrm{vd}\Omega\)