Équations de la thermique

Équation locale de la chaleur

L'équation de la chaleur, en tout point \(x\) d'un domaine matériel \(\Omega\), s'écrit :

(1)\[\frac{\partial H}{\partial t} + \vec{\nabla}.\vec{\phi} = q\]

avec :

\(H(T)\) l'enthalpie volumique (en J.m^-3)
\(\vec{\phi}(T,t)\) la densité de flux thermique (en W.m^-2)
\(q(T,t)\) la puissance volumique, terme source (en W.m^-3)
\(t\) le temps (en s)
\(T(x,t)\) la température (en K), principale inconnue du problème

L'enthalpie volumique \(H\) ne dépend que de la température, on peut donc écrire :

(2)\[\frac{\partial H}{\partial t} = \frac{\partial H}{\partial T} \frac{\partial T}{\partial t} = \rho c_p \frac{\partial T}{\partial t}\]

avec :

\(\rho\) la masse volumique (en kg.m^-3)
\(c_p\) la capacité thermique massique (en J.kg^-1.K^-1)

Le flux thermique \(\vec{\phi}\) est relié au gradient de température par la loi de Fourier :

\[\vec{\phi}=-\lambda \vec{\nabla} T\]

où \(\lambda\) est la conductivité thermique (en W.m^-1.K^-1)

Dans le cas où il n'y a pas de changement de phase et où les caractéristiques ne dépendent pas de la température, l'équation de la chaleur devient :

(3)\[\rho c_p \frac{\partial T}{\partial t} - \vec{\nabla}.\left(\lambda \vec{\nabla}T\right) = q\]

Conditions aux limites

Les conditions aux limites peuvent porter soit sur la température \(T\) :

(4)\[T=T_{\textrm{imp}} \quad \textsf{sur } \Gamma_T\]

soit sur le flux \(\vec{\phi}\) à travers un bord de normale \(\vec{n}\) :

(5)\[\begin{split}\begin{align} \vec{\phi}.\vec{n} & = \phi_{\textrm{imp}} & \textsf{sur } & \Gamma_{\phi} & \textsf{flux imposé} \\ \vec{\phi}.\vec{n} & = h(T_f - T) & \textsf{sur } & \Gamma_c & \textsf{convection} \\ \vec{\phi}.\vec{n} & = \varepsilon \sigma (T_{\infty}^4 - T^4) & \textsf{sur } & \Gamma_r & \textsf{rayonnement} \\ \end{align}\end{split}\]

avec :

\(h\) le coefficient d'échange (en W.m^-2.K^-1)
\(T_f\) la température du fluide en paroi (en K)
\(T_{\infty}\) la température à l'infini (en K)
\(\varepsilon\) l'émissivité de la surface
\(\sigma\) la constante de Stefan-Boltzmann (égale à en 5,67 10^-8 W.m^-2.K^-4)

Formulation éléments finis

Formulation faible

En multipliant chaque terme de l'équation (1) par un champ de température virtuel \(\tau\) non nul, puis en intégrant par partie sur \(\Omega\), on obtient une formulation faible de l'équation de la chaleur :

(6)\[\int_{\Omega} \left(\frac{\partial H}{\partial T}\frac{\partial T}{\partial t}\right)\tau d\Omega - \int_{\Omega} \lambda \vec{\nabla}T.\vec{\nabla}\tau d\Omega + \int_{\Gamma} \left(\lambda \vec{\nabla}T\right).\vec{n}\tau d\Gamma\]

L'intégrale sur le bord \(\Gamma\) peut se réduire aux bords portant les conditions aux limites (5) de flux imposé :

(7)\[\int_{\Gamma} \left(\lambda \vec{\nabla}T\right).\vec{n}\tau d\Gamma = - \int_{\Gamma_{\phi}}\phi_{\textrm{imp}}\tau d\Gamma_{\phi} - \int_{\Gamma_c}h(T_f-T)\tau d\Gamma_c - \int_{\Gamma_r}\varepsilon\sigma(T_{\infty}^4-T^4)\tau d\Gamma_r\]

Le terme du rayonnement peut être ramené à un terme de convection en faisant l'approximation suivante :

(8)\[T_{\infty}^4-T^4 = (T_{\infty}-T)(\tilde{T}^3 + \tilde{T}^2 T_{\infty} + \tilde{T} T_{\infty}^2 + T_{\infty}^3)\]

Cette approximation, permettant de linéariser le flux en fonction de \(T\), n'est valable qu'au voisinage de la température \(\tilde{T}\).

Discrétisation par éléments finis

En discrétisant l'espace \(\Omega\) par un maillage, les températures \(T\) et \(\tau\) sont discrétisées sur la base des fonctions d'interpolation \(\mathbfcal{N}\) :

\[T(x,t) = \mathbfcal{N}(x).T(t) = \sum_i \mathcal{N}_i(x) T_i(t)\]

\(T\) désignant maintenant le vecteur des températures aux noeuds et \(T_i\) la valeur de la température au noeud \(i\) du maillage. En injectant cette discrétisation dans la formulation faible (6), nous obtenons le système matriciel pour tout instant \(t\) :

(9)\[\mathbfcal{C}(T).\dot{T} + \mathbfcal{K}(T).T = Q\]

\(\mathbfcal{C}\) est la matrice de capacité :

(10)\[\mathbfcal{C}(T) = \int_{\Omega} \mathbfcal{N}^T.\left(\frac{\partial H}{\partial T}\frac{\partial T}{\partial t}\right).\mathbfcal{N} d\Omega\]

\(\mathbfcal{K}\) est la matrice de conductivité :

(11)\[\mathbfcal{K}(T) = \int_{\Omega} \nabla\mathbfcal{N}^T.\lambda(T).\nabla\mathbfcal{N} d\Omega + \int_{\Gamma_c} \mathbfcal{N}^T h(T) \mathbfcal{N} d\Gamma_c + \int_{\Gamma_r} \mathbfcal{N}^T \varepsilon\sigma \tilde{T}^3 \mathbfcal{N} d\Gamma_r\]

\(Q\) est le vecteur des puissances thermiques :

(12)\[Q = \int_{\Omega} \mathbfcal{N}^T q d\Omega + \int_{\Gamma_{\phi}} \mathbfcal{N}^T \phi_{\textrm{imp}} \mathbfcal{N} d\Gamma + \int_{\Gamma_c} \mathbfcal{N}^T hT_f d\Gamma_c + \int_{\Gamma_r} \mathbfcal{N}^T \varepsilon\sigma T_{\infty}^4 d\Gamma_r\]

Prise en compte des blocages/relations

La prise en compte des blocages ou des relations en température est faite de manière similaire à la mécanique. Les conditions de températures imposées (4) s'écrivent à l'aide d'une matrice de blocage \(\mathbfcal{A}\) :

(13)\[\mathbfcal{A}.T = T_{\textrm{imp}}\]

Cette égalité sera adjointe au système (9) et résolue à l'aide de multiplicateurs de Lagrange.

Opérateurs de Cast3M associés

Les termes des équations (9) et (13) sont calculés à l'aide des opérateurs suivants :

\(\mathbfcal{C}\) : CAPA pour la matrice de capacité
\(\mathbfcal{K}\) : COND pour les 2 premiers termes de conductivité et de convection
\(Q\) : SOUR (puissance volumique \(q\)), FLUX (flux imposé \(\phi_{\textrm{imp}}\)), CONV (convection)
\(\mathbfcal{A}\) : BLOQ (blocages), RELA (relations)
\(T_{\textrm{imp}}\) : DEPI (valeur des blocages ou relations)
La résolution du problème (9) nécessite la mise en oeuvre d'un schéma numérique d'intégration temporelle. Plusieurs méthodes sont proposées dans la procédure PASAPAS et décrites ci-après.
Les termes de rayonnement présents au premier membre \(\mathbfcal{K}\) (11) et second membre \(Q\) (12) sont calculés par la procédure PAS_RAYO