🧮 Logiciel de calcul thermique

But d'un logiciel de calcul thermique

En science et en ingénierie, la simulation thermique dynamique (STD) désigne le calcul de l'évolution temporelle de l'état thermique d'un système utilisant un modèle numérique approché de l'objet réel : par défaut, on obtient, à tous les instants choisis de la simulation, la température en un certain nombre de points des éléments le composant et qui évoluent selon les différentes lois régissant les échanges thermiques (convection, conduction, rayonnement, changement d'état).

Ici, nous allons écrire l'équation de l'évolution de l'énergie interne du système appliqué au volume d'air du bâtiment.

$$ \frac{dE}{dt} = m.{C_{p}}_{air}.\frac{T_{i}}{dt} = \sum \Phi_{Gains} - \sum \Phi_{Pertes} = \sum \Phi_{Internes} + \sum \Phi_{Solaire} - \left(\sum \Phi_{Transmission} + \sum \Phi_{Ventilation}\right) $$

Le but d'un logiciel de calcul thermique est de répérer les différents éléments pouvant créer des pertes ou des gains thermiques et de compiler ces données pour déterminer un modèle. Les parties suivantes regroupent les principales dissipations ou gain thermiques que l'on peut voir dans un bâtiment :

Pertes thermiques

Murs : Pertes par transmission

Les pertes par transmission reprennent toutes les pertes de chaleur à travers les éléments séparatifs entre le bâtiment et l'environnement extérieur, le sol et des structures voisines non chauffées.

$$ \Phi_{Transmission} = \sum_{m} U_{m}A_{m} \left( T_{i}-T_{e} \right) $$

Avec :

• $U_{m}$ : Coefficient de transmission thermique du mur "m" ($W/m^{2}.K$)
• $A_{m}$ : Surface de la paroi "m" ($m^{2}$)
• $T_{i}$ : Température intérieure ($Kelvin$)
• $T_{e}$ : Température extérieure ($Kelvin$)

Ventilation : Pertes par infiltration

Les pertes thermiques ne proviennent pas uniquement des surfaces isolées, mais aussi de la ventilation et des infiltrations d’air parasites.

$$ \Phi_{Ventilation} = \rho.C_{p}.\dot{V} \left( T_{i}-T_{s} \right) = 0.34nV \left( T_{i}-T_{s} \right) $$

Avec :

• $\rho$ : Masse volumique de l'air
• $C_{p}$ : Capacité calorifique de l'air
• $\dot{V}$ : Débit d'air dans le local
• $0.34$ : Approximation pour l'air : $\rho.C_{p}.\dot{V}/3600 = 0.34$)
• $n$ : Taux de renouvellement d'air ($Vol/h$) ($m^{2}$)
• $V$ : Volume du bâtiment
• $T_{i}$ : Température intérieure ($Kelvin$)
• $T_{S}$ : Température de soufflage ($Kelvin$)

Gains thermiques

Internes : Occupation, éclairage, équipements, chauffage

Les gains énergétiques internes peuvent venir de toutes les sources énergétiques. Cela peut être les équipements du local en fonctionnement dissipant de la chaleur, l'occupation du local ou encore la chaleur émise par l'éclairage.

Equipements/Chauffage : La dissipation peut être fixe s'il s'agit d'un équipement électrique ou variable s'il s'agit d'une tuyauterie

Eclairage : Des puissances surfaciques conventionnelles sont données dans la RE2020 (Logements : $1-2W/m^{2}$, Bureaux : $5-10W/m^{2}$)

Activité humaine : Des valeurs forfaitaires dépendantes de l'activité (Logements : $80-100W/pers$, Tertiaire : $70-120W/pers$)

Solaires : Rayonnement à travers le vitrage

Les gains solaires représente la puissance apportée par le rayonnement solaire au niveau des fenêtres.

$$ \Phi_{solaire} = \sum_{m} A_{m}.g_{m}.F_{m}.I_{m} $$

Avec :

• $A_{m}$ : Surface de la paroi "m" ($m^{2}$)
• $g_{m}$ : Facteur solaire pa la paroi "m" ($\%$)
• $F_{m}$ : Facteur d'ombrage de la paroi "m" ($\%$)
• $I_{m}$ : Irradiance solaire au niveau de cette paroi ($J/m^{2}$)

Solaires : Evolution de la température au niveau des murs.

Pour modéliser l'impact solaire au niveau des murs, il est nécessaire de calculer une température de surface extérieure du mur prenant en compte l'échauffement de la surface extérieure. La température de surface du mur est calculée de la façon suivante :

$$ T_{Surface~ext} = T_{ext} + \frac{\alpha .I}{h_{e}} $$

Avec :

• $T_{ext}$ : Température extérieure ($Kelvin$)
• $\alpha$ : Coefficient d'absorption du mur (peut évoluer avec la couleur, généralement pris entre 0.3 et 0.9)
• $g_{m}$ : Irradiance solaire ($W·m^{-2}$)
• $F_{m}$ : Coefficient d'échange extérieur ($W·m^{-2}·K^{-1}$)

Résolution matricielle du problème

Résolution mathématique

Lorsque l'on intègre toutes ces données dans la même équation, on a alors la formule suivante : $$ \underbrace{m.{C_{p}}_{air}.\frac{T_{i}}{dt}}_{\frac{dE}{dt}} = \sum \Phi_{Internes} + \underbrace{\sum_{m} A_{m}.g_{m}.F_{m}.I_{m}}_{\sum \Phi_{Solaire}} - \left( \underbrace{\sum_{m} U_{m}A_{m} \left( T_{i}-T_{e} \right)}_{\sum \Phi_{Transmission}} + \underbrace{\rho.C_{p}.\dot{V} \left( T_{i}-T_{e} \right)}_{\sum \Phi_{Ventilation}} \right) $$

Le but d'un logiciel thermique dynamique est de simuler l'évolution de la température de l'air en fonction du temps. L'ordinateur va alors calculer l'évolution de température estimée pour un petit écart de temps. Pour chaque local, on va donc discrétiser la formule selon la température du local :

$$ \begin{matrix} \frac{T_{i}}{dt} = \frac{ T_{i}.\left(\sum_{m} U_{m}.A_{m} + \rho.C_{p}.\dot{V} \right) - \left(\sum_{m} U_{m}.A_{m}.T_{e} + \rho.C_{p}.\dot{V}.T_{s} - \sum \Phi_{Internes} - \sum_{m} A_{m}.g_{m}.F_{m}.I_{m} \right) }{m.{C_{p}}_{air}} \\ \frac{T_{i}}{dt} = \frac{T^{n+1}_{i} - T^{n}_{i}}{\Delta T} \\ T^{n+1}_{i} = T^{n}_{i} + \frac{\Delta t}{m.{C_{p}}_{air}} \left[ \left( \sum_{m} U_{m}.A_{m}.T_{e} + \rho.C_{p}.\dot{V}.T_{s} - \sum \Phi_{Internes} - \sum_{m} A_{m}.g_{m}.F_{m}.I_{m} \right) - T^{n}_{i}.\left( \sum_{m} U_{m}.A_{m} + \rho.C_{p}.\dot{V} \right) \right] \end{matrix} $$

On peut donc réduire pour chaque local à l'équation suivante :

$$ T^{n+1}_{i} = T^{n}_{i}.\left( 1 - \mu_{\beta} \right) + \mu_{\alpha} $$

En sachant que pour chaque local, les valeurs sont les suivantes : $$ \begin{matrix} \mu_{\alpha} = \frac{\Delta t.\left(\sum_{m} U_{m}.A_{m}.T_{e} + \rho.C_{p}.\dot{V}.T_{s} - \sum \Phi_{Internes} - \sum_{m} A_{m}.g_{m}.F_{m}.I_{m}\right)}{m.{C_{p}}_{air}} \\ \mu_{\beta} = \frac{\Delta t.\left(\sum_{m} U_{m}.A_{m} + \rho.C_{p}.\dot{V}\right)}{m.{C_{p}}_{air}} \end{matrix} $$

Cette simplification marche très bien si l'on considère un bâtiment dont les locaux communiquent sans obstacles limitant le mélange de l'air. Si la modélisation doit prendre en compte les différences de températures entre locaux et leurs intéractions, il est nécessaire de modifier légèrement l'équation. On considèrera un local adjacent "$j$".

$$ \begin{matrix} \mu_{\alpha} = \frac{\Delta t.\left(\rho.C_{p}.\dot{V}.T_{s} - \sum \Phi_{Internes} - \sum_{m} A_{m}.g_{m}.F_{m}.I_{m}\right)}{m.{C_{p}}_{air}} \\ \mu_{\beta} = \frac{\Delta t.\rho.C_{p}.\dot{V}}{m.{C_{p}}_{air}} \\ \mu_{\gamma} = \frac{\Delta t.U_{m}.A_{m}}{m.{C_{p}}_{air}} \end{matrix} $$

Ce qui donne l'équation suivante :

$$ T^{n+1}_{i} = T^{n}_{i}.\left( 1 - \mu_{\beta} - \sum_{m} \mu_{\gamma} \right) + \mu_{\alpha} + \sum_{m_{int}} \mu_{\beta}.T_{j} + + \sum_{m_{ext}} \mu_{\gamma}.T_{ext} $$

Placement matricielle

Exemple global

Pour cette partie, nous allons supposer un bâtiment très simple avec 3 locaux.

Pour optimiser les calculs sur l'ordinateur, nous allons utiliser des matrices et des vecteurs. $$ \left\{ T^{n+1} \right\} = \left\{ T^{n} \right\}.\left[ M \right] + T_{e}.\left\{ Ext \right\} + \left\{ A \right\} $$

Avec :

• $\left\{ T^{n+1} \right\}$, $\left[ T^{n} \right]$ : Vecteur des températures des différents noeuds du système. Ce vecteur est à la fois la donnée d'entrée de l'équation représentant les températures au temps $n$, $n+1$, $n+2$...
• $\left[ M \right]$ matrice de calcul des éléments affectés par la températures des locaux ($T^{n}_{i}$)
• $\left\{ A \right\}$ vecteur des puissances associés aux locaux
• $\left\{ Ext \right\}$ vecteur des intéractions extérieures

Pour commencer, nous allons remplir la matrice de calcul :

$$ \left[ M \right] = \left[ \begin{array}{c|ccc} & \text{A} & \text{B} & \text{C} \\ \hline \text{A} & \underbrace{ \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right)}_{\text{Local A}} & \underbrace{{\mu_{\gamma}}_{B}}_{\text{Flux paroi A vs B}} & \underbrace{{\mu_{\gamma}}_{C}}_{\text{Flux paroi A vs C}} \\ \text{B} & {\mu_{\gamma}}_{A} & \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) & {\mu_{\gamma}}_{C} \\ \text{C} & {\mu_{\gamma}}_{A} & {\mu_{\gamma}}_{B} & \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) \end{array} \right] $$

La diagonale correspond à l'intéraction du local avec lui-même. Chaque croisement entre deux locaux doivent comporter le $\mu_{\gamma}$ correspondant à la paroi entre ces deux locaux. Les coefficients d'échanges convectifs peuvent être inclus au niveau du coefficient d'échange thermique de la paroi.

Ensuite, le vecteur correspondant aux échanges avec l'extérieur. Ce vecteur permet d'associer les coefficients d'échange à travers la paroi à la tyempérature extérieure pour affecter ces échanges à la température du local :

$$ \left\{ Ext \right\} = \left\{ \begin{array}{c} {\mu_{\gamma}}_{A} \\ {\mu_{\gamma}}_{B} \\ {\mu_{\gamma}}_{C} \\ \end{array} \right\} $$

Le vecteur des puissances des puissances va lui apporter ou retirer de l'énergie du système en fonction des équipements présents dans le local en fonction de son volume :

$$ \left\{ A \right\} = \left\{ \begin{array}{c} {\mu_{\alpha}}_{A} \\ {\mu_{\alpha}}_{B} \\ {\mu_{\alpha}}_{C} \\ \end{array} \right\} $$

Pour finir, le vecteur des températures va être créé. Chaque élément de ce vecteur $\left\{ T^{n} \right\}$ va correspondre à la température d'un local. Lorsque ce vecteur est utilisé dans la formule une fois, les évolutions de valeurs correspondent à l'évolution de température dans le local associé durant la période $\Delta T$. A $t=0$, les températures sont souvent initialisées à la même valeur. On parle de température d'initialisation $T_{i}$.

$$ \left\{ T^{n} \right\} = \left\{ \begin{array}{c} {T_{i}}_{A} \\ {T_{i}}_{B} \\ {T_{i}}_{C} \\ \end{array} \right\} $$

Cela nous donne l'équation matricielle suivante :

$$ \left\{ \begin{array}{c} {T^{n+1}}_{A} \\ {T^{n+1}}_{B} \\ {T^{n+1}}_{C} \\ \end{array} \right\} = \left\{ \begin{array}{c} {T^{n}}_{A} \\ {T^{n}}_{B} \\ {T^{n}}_{C} \\ \end{array} \right\} . \left[ \begin{matrix} \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) & {\mu_{\gamma}}_{B} & {\mu_{\gamma}}_{C} \\ {\mu_{\gamma}}_{A} & \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) & {\mu_{\gamma}}_{C} \\ {\mu_{\gamma}}_{A} & {\mu_{\gamma}}_{B} & \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) \end{matrix} \right] + T_{e}.\left\{ \begin{array}{c} {\mu_{\gamma}}_{A} \\ {\mu_{\gamma}}_{B} \\ {\mu_{\gamma}}_{C} \\ \end{array} \right\} + \left\{ \begin{array}{c} {\mu_{\alpha}}_{A} \\ {\mu_{\alpha}}_{B} \\ {\mu_{\alpha}}_{C} \\ \end{array} \right\} $$

Exemple spécifique

Pour cette partie et pour mieux comprendre le placement des éléments, nous allons considérer des locaux agencés de la manière suivante.

• Local A : Le local A est le seul local en contact avec l'extérieur. Ce local est en contact avec seulement le local B.
• Local B : Le local B est en contact avec le local A et le local C.
• Local C : Le local C est en contact uniquement avec le local B.

Dans ce cas, l'équation matricielle doit être la suivante :

$$ \left\{ \begin{array}{c} {T^{n+1}}_{A} \\ {T^{n+1}}_{B} \\ {T^{n+1}}_{C} \\ \end{array} \right\} = \left\{ \begin{array}{c} {T^{n}}_{A} \\ {T^{n}}_{B} \\ {T^{n}}_{C} \\ \end{array} \right\} . \underbrace{ \left[ \begin{matrix} \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) & {\mu_{\gamma}}_{B} & 0 \\ {\mu_{\gamma}}_{A} & \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) & {\mu_{\gamma}}_{C} \\ 0 & {\mu_{\gamma}}_{B} & \left(1 - \mu_{\beta} - \sum_{m} \mu_{\gamma}\right) \end{matrix} \right] }_{\text{Les locaux A et C ne sont pas en contact}} + T_{e}. \underbrace{ \left\{ \begin{array}{c} {\mu_{\gamma}}_{A} \\ 0 \\ 0 \\ \end{array} \right\} }_{\text{Seul A en contact exterieur}} + \left\{ \begin{array}{c} {\mu_{\alpha}}_{A} \\ {\mu_{\alpha}}_{B} \\ {\mu_{\alpha}}_{C} \\ \end{array} \right\} $$

Modéliser les équipements

Puissance fixe

La modélisation des équipements va passer par la description du vecteur des puissances $\left\{ A \right\} =\left\{\begin{array}{c}{\mu_{\alpha}}_{A} \\{\mu_{\alpha}}_{B} \\{\mu_{\alpha}}_{C} \\\end{array}\right\}$. Ce vecteur doit estimer l'impact de réchauffement ou de refroidissement des équipements sur la température du local.

Pour chaque local, on a donc l'équation suivante : $\partial T = \frac{\Phi}{m.{C_{p}}_{air}}$

Puissance variable

En cas de puissance variable, les équipements vont s'activer et se désactiver en fonction de seuils de température. Ces seuils de température vont soit activer l'équipement soit le désactiver. On parle de seuil d'enclenchement la valeur de température à laquelle l'équipement s'active. On appelle seuil de déclenchement la valeur de température à laquelle l'équipement se désactive.

On peut voir 2 cas différents ici, soit la température d'enclenchement est supérieure à la température de déclenchement (clim), soit la température d'enclenchement est inférieure à la température de déclenchement (chauffage).

Le chauffage peut être gérer en respectant la formule suivante : $\text{Statut} = \left( \text{Statut} \lor T < T_{encl} \right) \land T < T_{decl}$.

La climatisation peut être gérer en respectant la formule suivante : $\text{Statut} = \left( \text{Statut} \lor T > T_{encl} \right) \land T > T_{decl}$.

Chauffage			Climatisation
T°	Montée en température	Descente de température	T°	Montée en température	Descente de température
$T > T_{decl}$ $T_{decl}$ $T_{decl} > T > T_{encl}$ $T_{encl}$ $T < T_{encl}$	$0$ $0$ $1$ $1$ $1$	$0$ $0$ $0$ $1$ $1$	$T > T_{encl}$ $T_{encl}$ $T_{encl} > T > T_{decl}$ $T_{decl}$ $T < T_{decl}$	$1$ $1$ $0$ $0$ $0$	$1$ $1$ $1$ $0$ $0$

On peut remarquer que le statut lors du chauffage est l'inverse de celui de la climatisation. C'est à dire que : $\text{Statut}_{Chauffage} = \neg \text{Statut}_{Clim}$.

Si l'on veut faire une seule equation qui va permettre de gérer les seuils quelque soit leurs valeurs, l'équation suivante va être utilisée :

$$ \text{Statut}_{Equipement} = \underbrace{\left( T_{decl} > T_{encl} \right)}_{\text{1 si chauffage, 0 si clim}} $$