💨 Convection thermique

De Navier-Stokes à Lorenz

Pour déterminer le système d'équation de Lorenz, nous devons partir du système d'équation de Navier-Stokes:

$$ \left\lbrace \begin{array}{l} \frac{\partial \rho}{\partial t}+\nabla(\rho \vec{v})=0 \\ \frac{\partial v}{\partial t}+\left( v.\nabla \right)v = F - \frac{\nabla p}{\rho} + \frac{\mu}{\rho} \nabla^{2}v\\ \frac{\partial T}{\partial t} + \vec{v}.\nabla T = \lambda.\Delta T \end{array} \right. $$

Pour simplifier le problème, nous considérons une cellule de convection en deux dimensions de hauteur h, d'une température minimale au sommet T₀ et d'une température maximale au point le plus bas T₀+δT.

$$ T_{1} = T_{0}+\delta T $$

Les conditions aux limites sont les suivantes:

$$ \left\lbrace \begin{array}{l} T \left( x, y, z=0, t \right)= T_{0} +\delta T\\ T \left( x, y, z=h, t \right)= T_{0} \end{array} \right. $$

Nous allons supposer que notre problème est indépendant selon y. Nous allons également postuler que les coefficients ne varient pas selon la température. Nous pouvons alors modifier l'équation de la conservation de la masse:

$$ \frac{\partial v_{x}}{\partial x}+\frac{\partial v_{z}}{\partial z}=0 $$

Nous introduisons la fonction courant : ψ tel que $\psi$ tel que $\vec{v} = \vec{rot} (\psi)$. Ce qui nous donne:

$$ \left\lbrace \begin{array}{l} v_{x} = \frac{\partial \psi}{\partial y}-\frac{\partial \psi}{\partial z}=-\frac{\partial \psi}{\partial z}\\ v_{y} = \frac{\partial \psi}{\partial z}-\frac{\partial \psi}{\partial x}=0\\ v_{z} = \frac{\partial \psi}{\partial x}-\frac{\partial \psi}{\partial y}=\frac{\partial \psi}{\partial x}\\ \end{array} \right. $$

En tenant compte du système précédent et de celui de Navier Stokes, nous pouvons poser les équations suivantes:

$$ \left\lbrace \begin{array}{l} %equation de conservation de la masse \frac{\partial v_{x}}{\partial x}+\frac{\partial v_{z}}{\partial z}=0\\ %bilan de mouvement \left\lbrace \begin{array}{l} \frac{\partial v_x}{\partial t}+ v_x \frac{\partial v_x}{\partial x}+ v_z \frac{\partial v_x}{\partial z} = - F_x - \frac{1}{\rho} \frac{\partial P}{\partial x} + \frac{\mu}{\rho}\nabla^{2}v_x\\ \frac{\partial v_z}{\partial t}+ v_x \frac{\partial v_z}{\partial x}+ v_z \frac{\partial v_z}{\partial z} = - F_z - \frac{1}{\rho} \frac{\partial P}{\partial z} + \frac{\mu}{\rho}\nabla^{2}v_z\\ \end{array} \right. \\ %Energie \frac{\partial T}{\partial t} + \left\lbrace \begin{array}{c} -\frac{\partial \psi}{\partial z}\\ 0\\ \frac{\partial \psi}{\partial x}\\ \end{array} \right\rbrace .\nabla T = \lambda.\Delta T \end{array} \right. $$

Pour simplifier les calculs, nous posons la fonction θ(x,z,t) tel que:

$$ T(x,z,t) = T_{0} + \delta T - \frac{\delta T}{h}z + \theta(x,z,t) $$

Dans l'équation de la quantité de mouvement, nous avons F qui représente les forces massiques s'exerçant dans le fluide. Dans notre cas, ces forces sont représentées par la force de gravitation. Nous pouvons donc transformer F en vecteur :

$$ F = \left\lbrace \begin{array}{c} 0\\ 0\\ \rho\left(T\right).g \end{array} \right\rbrace $$

Pour éliminer la pression de notre équation du mouvement, nous allons faire la transformation suivante: $\frac{\partial}{\partial z}$(équation du mouvement selon x)$-\frac{\partial}{\partial x}$(équation du mouvement selon z).

Cela va nous donner l'équation suivante:

$$ \begin{split} \frac{\partial}{\partial t} \left( \frac{\partial v_x}{\partial z}-\frac{\partial v_z}{\partial x} \right) +\frac{\partial v_x}{\partial z}.\frac{\partial^2 v_x}{\partial x \partial z}+\frac{\partial v_z}{\partial z}.\frac{\partial^2 v_x}{\partial^2 z}-\frac{\partial v_x}{\partial x}.\frac{\partial^2 v_z}{\partial^2 x}-\frac{\partial v_z}{\partial x}.\frac{\partial^2 v_z}{\partial x \partial z}\\ =\frac{\mu}{\rho}\left( \frac{\partial \nabla^2 v_x}{\partial z} - \frac{\partial \nabla^2 v_z}{\partial x} \right) + \frac{\partial}{\partial x}\left( \frac{\rho\left(T\right)}{\rho} \right)g \end{split} $$

En prenant en compte les équations précédentes, nous pouvons ré-écrire le système de Navier Stokes:

$$ \left\lbrace \begin{array}{l} %equation de conservation de la masse \frac{\partial v_{x}}{\partial x}+\frac{\partial v_{z}}{\partial z}=0\\ %Quantité de mouvement \frac{\partial}{\partial t}\nabla^2 \psi = -\frac{\partial \left( \psi, \nabla^2 \psi \right)}{\partial \left( x, z \right)} + \nu \nabla^4 \psi + \frac{\partial}{\partial x}\left( \frac{\rho\left(T\right)}{\rho} \right)g\\ %Energie \frac{\partial \theta}{\partial t} = -\frac{\partial \left( \psi, \theta \right)}{\partial \left( x, z \right)} + \frac{\delta T}{h}\frac{\partial \psi}{\partial x} + K \nabla^2 \theta \end{array} \right. $$

Nous devons considérer le ρ(T) variable. En effet, le phénomène de convection est généré par les changement de densité de notre fluide qui se fait suivant la température de ce dernier. Nous allons ici utiliser l'approximation de Boussinesq-Oberheck. Cette approximation nous donne la formule suivante:

$$ \alpha \equiv -\frac{1}{\rho}\left( \frac{\partial \rho^{-1}}{\partial T} \right)_{P=constante} $$

Pour un ΔT faible, nous pouvons alors avoir la fonction suivante:

$$ \rho \left( T \right) = \rho_0 \left( 1 - \alpha \left( T-T_0 \right) \right) $$

Le système de Navier Stokes devient:

$$ \left\lbrace \begin{array}{l} %equation de conservation de la masse \frac{\partial v_{x}}{\partial x}+\frac{\partial v_{z}}{\partial z}=0\\ %Quantité de mouvement \frac{\partial}{\partial t}\nabla^2 \psi = -\frac{\partial \left( \psi, \nabla^2 \psi \right)}{\partial \left( x, z \right)} + \nu \nabla^4 \psi + g.\alpha \frac{\partial \theta}{\partial x}\\ %Energie \frac{\partial \theta}{\partial t} = -\frac{\partial \left( \psi, \theta \right)}{\partial \left( x, z \right)} + \frac{\delta T}{h}\frac{\partial \psi}{\partial x} + K \nabla^2 \theta \end{array} \right. $$

Nous pouvons alors faire des hypothèses sur certains points de notre problème:

• Lorsque z=0 et z=h, nous pouvons observer que θ est nul. Nous pouvons alors poser: θ(x, 0, t) = θ(x, h, t) = 0
• Nous savons également que notre courant est nul au sommet et à la base de notre cellule de convection. Nous pouvons alors poser: ψ(x, 0, t) = ψ(x, h, t) = 0
• Enfin, nous allons également poser: ∇² ψ(x, 0, t) = ∇² ψ(x, h, t) = 0

Avec cela, nous allons pouvoir déterminer notre système de Lorenz. Pour cela, nous allons développer nos fonctions θ^*(x^*, z^*, t^*) et ψ^*(x^*, z^*, t^*) en séries de Fourier. Pour cela, nous allons poser:

$$ \begin{array}{c c c} x = h.x^{*} & z = h.z^{*} & t = \frac{h^2}{K}t^{*} \end{array} $$

Les séries de Fourier sont un outil très utilisé dans l'étude des fonctions périodiques. Notre cellule de convection répond donc à ce critère. Nous allons donc développer les fonctions ψ^* et θ^* avec une longueur d'onde égale à l dans la direction de x et 2h dans la direction de z. Nous avons donc les développements suivant:

$$ \begin{array}{l} \psi^{*}(x^{*}, z^{*}, t^{*}) = \sum\limits_{m=-\infty}^{\infty} \sum\limits_{n=-\infty}^{\infty} \psi_{mn}(m, n, t^{*}).e^{2\pi i\left( \frac{mh}{l}x^{*}+\frac{n}{2}z^{*} \right)}\\ \theta^{*}(x^{*}, z^{*}, t^{*}) = \sum\limits_{m=-\infty}^{\infty} \sum\limits_{n=-\infty}^{\infty} \theta_{mn}(m, n, t^{*}).e^{2\pi i\left( \frac{mh}{l}x^{*}+\frac{n}{2}z^{*} \right)} \end{array} $$

Pour le développement des séries, Lorenz a choisi seulement 3 amplitudes: X(t), Y(t), Z(t). Il a donc retrouvé les formules suivantes:

$$ \begin{array}{l} \psi(x, z, t) = \frac{K(1+a^2)\sqrt{2}}{a}X(t) \sin\left( \frac{\pi a}{h}x \right) \sin\left( \frac{\pi}{h}z \right)\\ \theta(x, z, t) = \frac{\delta T. Ra_c}{\pi. Ra}\left[ \sqrt{2}Y(t)\cos\left( \frac{\pi a}{h}x \right) \sin\left( \frac{\pi}{h}z \right)-Z(t)\sin\left( \frac{2\pi}{h}z \right)\right] \end{array} $$

Où Ra est le nombre de Rayleigh et Ra_c est le nombre de Rayleigh critique. On a donc:

$$ \begin{array}{c c} Ra = \frac{\alpha g h^3 \delta T}{K \nu}& Ra_c = \frac{\pi^2 (1+a^2)^3}{a^2} \end{array} $$

En injectant ce résultat dans le système de Navier Stokes avant développement de Fourier, nous obtenons le système de Lorenz:

$$ \left\lbrace \begin{array}{l} \overset{\cdot}{x} = Pr (y - x)\\ \overset{\cdot}{y} = Rx - y - xz\\ \overset{\cdot}{z} = xy - \beta z \end{array} \right. $$

Où:

$$ \begin{array}{c c c} R = \frac{Ra}{Ra_c}&\beta = \frac{4}{()1+a^2} \end{array} $$

Méthode de Runge Kutta

Pour utiliser le système d'équation de Lorenz dans notre programme, nous avons recours aux dérivations de Runge Kutta d'ordre 4. Les formules sont faciles à programmer et les résultats sont généralement très proches de la fonction réelle. Cette méthode tire les avantages des méthodes de Taylor en gardant une grande simplicité d'exécution.

La formule de Runge Kutta à l'ordre 4 est la suivante pour une équation de la forme $y\left(t\right)$.

$$ \begin{array}{c} k_{1} = h \times f \left( t_{n}, y(t_{n}) \right)\\ k_{2} = h \times f \left( t_{n} + \frac{1}{2}h, y(t_{n}) +\frac{1}{2}k_{1} \right)\\ k_{3} = h \times f \left( t_{n} + \frac{1}{2}h, y(t_{n}) +\frac{1}{2}k_{2} \right)\\ k_{4} = h \times f \left( t_{n} + h, y(t_{n}) +k_{3} \right)\\ \\ y(t_{n+1}) = y(t_{n}) + \frac{1}{6}\left( k_{1} + 2.k_{2} + 2.k_{3} + k_{4} \right) + O(h^{4}) \end{array} $$

Où $O(h^{4})$ est l'erreur produite. Nous devons réfléchir au choix du pas de temps "h" nous permettant d'avoir un bon compromis entre la précision du résultat et le temps de calcul de l'ordinateur.

Dans notre cas, nous avons le système de Lorenz suivant:

$$ \begin{array}{l l l} \frac{dx}{dt} &= Pr \left( y-x \right) &= f(x,y)\\ \frac{dy}{dt} &= r.x - y -x.z &= g(x,y,z)\\ \frac{dz}{dt} &= x.y - b.z &= k(x,y,z) \end{array} $$

Nous allons donc avoir une approximation de notre système avec les suites suivantes:

$$ \begin{array}{l} X_{n+1} = X_n + \left( a_x + 2.b_x + 2.c_x + d_x \right)\\ Y_{n+1} = Y_n + \left( a_y + 2.b_y + 2.c_y + d_y \right)\\ Z_{n+1} = Z_n + \left( a_z + 2.b_z + 2.c_z + d_z \right) \end{array} $$

Avec:

$$ \begin{array}{l} a_x = h.f\left( X_n, Y_n \right)\\ a_y = h.g\left( X_n, Y_n, Z_n \right)\\ a_z = h.k\left( X_n, Y_n, Z_n \right)\\ \\ b_x = h.f\left( X_n+\frac{a_x}{2}, Y_n+\frac{a_x}{2} \right)\\ b_y = h.f\left( X_n+\frac{a_y}{2}, Y_n+\frac{a_y}{2}, Z_n+\frac{a_y}{2} \right)\\ b_z = h.f\left( X_n+\frac{a_z}{2}, Y_n+\frac{a_z}{2}, Z_n+\frac{a_z}{2} \right)\\ \\ c_x = h.f\left( X_n+\frac{b_x}{2}, Y_n+\frac{b_x}{2} \right)\\ c_y = h.f\left( X_n+\frac{b_y}{2}, Y_n+\frac{b_y}{2}, Z_n+\frac{b_y}{2} \right)\\ c_z = h.f\left( X_n+\frac{b_z}{2}, Y_n+\frac{b_z}{2}, Z_n+\frac{b_z}{2} \right)\\ \\ d_x = h.f\left( X_n+c_x, Y_n+c_x \right)\\ d_y = h.f\left( X_n+c_y, Y_n+c_y, Z_n+c_y \right)\\ d_z = h.f\left( X_n+c_z, Y_n+c_z, Z_n+c_z \right)\\ \end{array} $$

Nous pouvons ici nous rendre compte du nombre d'opérations nécessaires à la modélisation de notre attracteur étrange. La modélisation de notre phénomène est donc confié à un ordinateur.