The goal of this project is to simulate the fluid flow under natural convection: the heated fluid circulates towards the low temperature under the action of density and gravity differences. Thie phenomenon is important in the sense it models evacuation of heat, generated by friction forces for example, with a cooling fluid.
We shall put in place a simple convection problem in order to study the phenomenon without having to handle the difficulties of more complex domaines. We describe then some necessary transformations to the equations, then we define quantities of interest to be able to compare the simulations with different parameter values.
To study the convection, we use a model problem: it consists in a rectangular tank of height \(1\) and width \(W\), in which the fluid is enclosed :
![]() |
We wish to know the fluid velocity \(\mathbf{u}\), the fluid pressure \(p\) and fluid temperature \(\theta\).
We introduce the adimensionalized Navier-Stokes and heat equations parametrized by the Grashof and Prandtl numbers. These parameters allow to describe the various regimes of the fluid flow and heat transfer in the tank when varying them.
The adimensionalized steady incompressible Navier-Stokes equations reads:
\begin{equation*} \begin{split} \mathbf{u}\cdot\nabla \mathbf{u} +\nabla p -\frac{1}{\sqrt{\Gr}} \Delta \mathbf{u} &= \theta \mathbf{e}_2\\ \nabla \cdot \mathbf{u} & = 0\ \text{sur}\ \Omega\\ \mathbf{u} & = \mathbf{0}\ \text{sur}\ \partial \Omega \end{split} \end{equation*}
where \(\Gr\) is the Grashof number, \(\mathbf{u}\) the adimensionalized velocity and \(p\) adimensionalized pressure and \(\theta\) the adimensionalized temperature. The temperature is in fact the difference between the temperature in the tank and the temperature \(T_0\) on boundary \(\Gamma_1\).
The heat equation reads:
\begin{equation*} \begin{split} \mathbf{u} \cdot \nabla \theta -\frac{1}{\sqrt{\Gr}\Pr} \Delta \theta &= 0\\ \theta &= 0\ \text{sur}\ \Gamma_1\\ \frac{\partial \theta}{\partial n} &= 0\ \text{sur}\ \Gamma_{2,4}\\ \frac{\partial \theta}{\partial n} &= 1\ \text{sur}\ \Gamma_3 \end{split} \end{equation*}
where \(\Pr\) is the Prandtl number.
What are the effects of the Grashof and Prandtl numbers ? We remark that both terms with these parameters appear in front of the \(\Delta\) parameter, they thus act on the diffusive terms. If we increase the Grashof number or the Prandtl number the coefficients multiplying the diffusive terms decrease, and this the convection, that is to say the transport of the heat via the fluid, becomes dominant. This leads also to a more difficult and complex flows to simulate. The influence of the Grashof and Prandtl numbers are different but they generate similar difficulties and flow configurations. Thus we look only here at the influence of the Grashof number which shall vary in \([1, 1e7]\).
![]() |
\(h=0.01\) and \(\Pr=1\). |
We would like to compare the results of many simulations with respect to the Grashof defined in the previous section. We introduce two quantities which will allow us to observe the behavior of the flow and heat transfer.
We consider first the mean temperature on boundary \(\Gamma_3\)
\begin{equation*} T_3 = \int_{\Gamma_3} \theta \end{equation*}
This quantity should decrease with increasing Grashof because the fluid flows faster and will transport more heat which will cool down the heated boundary \(\Gamma_3\). We observe this behavior on this figure.
![]() |
\(h=0.02\) with \(\mathbb{P}_3\) Lagrange element for the velocity, \(\mathbb{P}_2\) Lagrange for the pressure and \(\mathbb{P}_1\) Lagrange for the temperature. |
Another quantity of interest is the flow rate through the middle of the tank. We define a segment \(\Gamma_f\) as being the vertical top semi-segment located at \(W/2\) with height \(1/2\). The flow rate, denoted \(\mathrm{D}_f\), reads
\begin{equation*} \mathrm{D}_f = \int_{\Gamma_f} \mathbf{u} \cdot \mathbf{e}_1 \end{equation*}
where \(\mathbf{e}_1=(1,0)\). Note that the flow rate can be negative or positive depending on the direction in which the fluid flows.
As a function of the Grashof, we shall see a increase in the flow rate. This is true for small Grashof, but starting at \(1e3\) the flow rate decreases. The fluid is contained in a boundary layer which is becoming smaller as the Grashof increases.
![]() |
\(h = 0.02\), \(\mathbb{P}_3\) for the velocity, \(\mathbb{P}_2\) for the pressure and \(\mathbb{P}_1\) for the temperature. |
section Implementation Implementation This application in implemented in "feel/doc/manual/heatns/convection.cpp"
. The implementation solve the full nonlinear problem using the nonlinear solver framework.
Consider the following problem,
\begin{equation*} \mbox{Stokes: }\left\{ \begin{array}{rcc} -\mu\Delta\mathbf{u} + \nabla p = \mathbf{f}\\ \nabla\cdot\mathbf{u} = 0\\ \mathbf{u}|_{\partial \Omega} = 0 \end{array} \right. \end{equation*}
where \(\Omega \subset \mathbb{R}^d\). There are no boundary condition on the pressure. This problem is ill-posed, indeed we only control the pressure through its gradient \(\nabla p\). Thus if \((\mathbf{u},p)\) is a solution, then \((\mathbf{u},p+c)\) is also a solution with \(c\) any constant. This comes from the way the problem is posed: the box is closed and it is not possible to determine the pressure inside. The remedy is to impose arbitrarily a constraint on the pressure, e.g. its mean value is zero. In other words, we add this new equation to the problem
\begin{equation} \int_\Omega p = 0 \end{equation}
In order to impose the previous condition, we introduce a new unknown, a Lagrange multiplier, \(\lambda \in \mathbb{R} \) and modify the incompressibility equation. Our problem reads now, find \((\mathbf{u},p,\lambda)\) such that
\begin{equation*} \mbox{Stokes 2: }\left\{ \begin{array}{rcl} -\mu\Delta\mathbf{u} + \nabla p &=& \mathbf{f}\\ \nabla\cdot\mathbf{u} + \lambda &=& 0\\ \mathbf{u}|_{\partial \Omega} &=& 0\\ \int_\Omega p &=& 0 \end{array} \right. \end{equation*}
The variational formulation now reads: find \((\mathbf{u}, p, \lambda) \in \mathbf{H}^1_0(\Omega) \times L^2_0(\Omega) \times \mathbb{R}\) such that for all \((\mathbf{v}, q, \eta) \in \mathbf{H}^1_0(\Omega) \times L^2_0(\Omega) \times \mathbb{R}\)
\begin{equation*} \label{notes:eq:20} \mbox{Stokes 3: }\left\{ \begin{array}{rcl} \int_\Omega \Big(\nabla \mathbf{u} \colon \nabla \mathbf{v} + \nabla \cdot \mathbf{v} p\Big) &=& \int_\Omega \mathbf{f} \cdot \mathbf{v}\\ \int_\Omega \Big(\nabla\cdot\mathbf{u} q + \lambda q\Big) &=& 0\\ \int_\Omega p \eta &=& 0 \end{array} \right. \end{equation*}
Summing up all three equations we get the following condensed formulation:
\begin{equation*} \label{notes:eq:19} \int_\Omega \nabla \mathbf{u} \colon \nabla \mathbf{v} + \nabla \cdot \mathbf{v} p + \nabla \cdot \mathbf{u} q + \lambda q + \eta p = \int_\Omega \mathbf{f} \cdot \mathbf{v} \end{equation*}
where \(\mathbf{H}^1_0(\Omega)= \Big\{ \mathbf{v} \in \mathbf{L}^2(\Omega), \nabla \mathbf{v} \in [L^2(\Omega)]^{d\times d},\ \mathbf{v} = 0\ \text{on}\ \partial \Omega \Big\}\), \(L^2_0(\Omega)= \Big\{ v \in L^2(\Omega),\ \int_\Omega v = 0\Big\}\), and \(\mathbf{L}^2(\Omega)= \Big\{ \mathbf{v} \in [L^2(\Omega)]^d\Big\}\) that is to say each component of a vector field of \(\mathbf{L}^2(\Omega)\) are in \(L^2(\Omega)\).
Consider the following steady incompressible Navier-Stokes equations, find \((\mathbf{u},p)\) such that
\begin{equation} \label{notes:eq:7} \begin{split} \underbrace{\rho \mathbf{u} \cdot \nabla \mathbf{u}}_{\text{convection}} - \underbrace{\nu \Delta \mathbf{u}}_{\text{diffusion}} + \nabla p &= \mathbf{f} \ \text{on}\ \Omega \\ \nabla \cdot \mathbf{u} &= 0 \\ \mathbf{u} &= \mathbf{0}\ \text{on}\ \partial \Omega \end{split} \end{equation}
where \(\rho\) is the density of the fluid, \(\nu\) is the dynamic viscosity of the fluid(la viscosité cinématique \(\eta = \nu/\rho\)) and \(\mathbf{f}\) is the external force density applied to the fluid, (e.g. \(\mathbf{f}=-\rho g \mathbf{e}_2\) with \(\mathbf{e}_2=(0,1)^T\) ). This equation system is nonlinear due to the \(\mathbf{u} \cdot \nabla \mathbf{u}\) convection term. A simple approach to solve \((2)\) is to use a fix point algorithm.
The fixpoint algorithm for NS reads as follows, find \((\mathbf{u}^{(k)},p^{(k)})\) such that
\begin{equation*} \label{notes:eq:13} \begin{split} \rho\mathbf{u}^{(k-1)} \cdot \nabla \mathbf{u}^{(k)} - \nu \Delta \mathbf{u}^{(k)} + \nabla p^{(k)} &= \mathbf{f} \ \text{on}\ \Omega \\ \nabla \cdot \mathbf{u}^{(k)} &= 0 \\ \mathbf{u}^{(k)} &= 0\ \text{on}\ \partial \Omega\\ (\mathbf{u}^{(0)},p^{(0)}) &= (\mathbf{0},0) \end{split} \end{equation*}
This system is now linear at each iteration \(k\) and we can write the variational formulation accordingly. A stopping criterium is for example that \(\|\mathbf{u}^{k}-\mathbf{u}^{(k-1)}\|+\|p^{k}-p^{(k-1)}\| < \epsilon\) where \(\epsilon\) is a given tolerance (e.g. \(1e-4\)) and \(\|\cdot\|\) is the \(L_2\) norm.
Here is the implementation using :
Recall that we have to solve two coupled problems :
\[ \mbox{Heat(\textbf{u}) }\left\{ \begin{array}{rcc} - \kappa\Delta T + \mathbf{u}\cdot\nabla T &=& 0 \\ T|_{\Gamma_1} &=& T_0 \\ \frac{\partial T}{\partial \mathbf{n}}|_{\Gamma_3} &=&1 \\ \frac{\partial T}{\partial \mathbf{n}}|_{\Gamma_2,\Gamma_4} &=& 0 \end{array} \right. \]
and
\[ \mbox{Stokes(T) : }\left\{ \begin{array}{rcc} -\nu\Delta\mathbf{u} + \frac{1}{\rho}\nabla p = \mathbf{F}\\ \nabla\cdot\mathbf{u} = 0\\ \mathbf{u}|_{\partial \Omega} = 0 \end{array} \right. \]
Where \(\mathbf{F}\) can be taken as \( \left( \begin{array}{c} 0 \\ \beta(T-T_0) \end{array} \right) \) for some \(\beta > 0\). \(\beta\) is called the dilatation coefficient.
Here is a simple algorithm fix point strategy in pseudo-code:
Another possiblity is to use a Newton method which allows us to solve the full nonlinear problem coupling velocity, pressure and temperature
\begin{equation*} \label{notes:eq:21} \text{Find}\ X\ \text{such that}\ F(X) = 0 \end{equation*}
the method is iterative and reads, find \(X^{(n+1)}\) such that
\begin{equation} \label{notes:eq:22} J_F(X^{(n)})( X^{(n+1)}-X^{(n)}) = - F (X^{(n)}) \end{equation}
starting with \(X^{(0)} = \mathbf{0}\) or some other initial value and where \(J_F\) is the jacobian matrix of \(F\) evaluated at \(X=((u_i)_i,(p_i)_i,(\theta_i)_i)^T\). For any \(\phi_k, \psi_l\) and \(\rho_m\) the test functions associated respectively to velocity, pressure and temperature, our full system reads, Find \(X=((u_i)_i,(p_i)_i,(\theta_i)_i)^T\) such that
\begin{equation*} \label{notes:eq:23} \begin{array}{rll} F_1((u_i)_i,(p_i)_i,(\theta_i)_i)&=\sum_{i,j} u_i u_j a(\phi_i,\phi_k,\phi_j) - \sum_i p_i b(\phi_k,\psi_i) + \sum_i \theta_i c(\rho_i, \phi_k)+\sum_i u_i d(\phi_i,\phi_k) &= 0\\ F_2((u_i)_i,(p_i)_i,(\theta_i)_i)&=\sum_i u_i b(\phi_i,\psi_l) &=0\\ F_3((u_i)_i,(p_i)_i,(\theta_i)_i)&=\sum_{i,j} u_i\theta_j e(\phi_i,\rho_j,\rho_m) + \sum_i \theta_i f(\rho_i,\rho_m)-g(\rho_m) &= 0 \end{array} \end{equation*}
where \(F=(F_1,F_2,F_3)^T\) and
\begin{equation*} \label{notes:eq:26} \begin{array}{rl} a(\mathbf{u},\mathbf{v},\beta) &= \int_\Omega \mathbf{v}^T ((\nabla \mathbf{u} )\beta)\\ b(\mathbf{v},p) &= \int_\Omega p (\nabla \cdot \mathbf{v}) - \int_{\partial \Omega} \mathbf{v}\cdot\mathbf{n} p\\ c(\theta,\mathbf{v})&= \int_\Omega \theta \mathbf{e}_2 \cdot \mathbf{v}\\ d(\mathbf{u},\mathbf{v}) &= \frac{1}{\sqrt{\mathrm{Gr}}} \Big(\int_\Omega \nabla \mathbf{u} \colon (\nabla \mathbf{v})^T - \int_{\partial \Omega} ((\nabla \mathbf{u}) \mathbf{n})\cdot \mathbf{v}\Big)\\ e(\mathbf{u},\theta,\chi) &= \int_\Omega (\mathbf{u}\cdot \nabla \theta) \chi \\ f(\theta,\chi) &=\frac{1}{\sqrt{\mathrm{Gr}}\mathrm{Pr}} \Big( \int_\Omega \nabla \theta \cdot \nabla \chi - \int_{\Gamma_1} (\nabla \theta \cdot \mathbf{n} ) \chi \Big)\\ g(\chi) &=\frac{1}{\sqrt{\mathrm{Gr}}\mathrm{Pr}} \int_{\Gamma_3} \chi \end{array} \end{equation*}
In order to apply the newton scheme, we need to compute the jacobian matrix \(J_F\) by deriving each equation with respect to each unknowns, ie \(u_i, p_i\) and \(\theta_i\). Consider the first equation
\begin{equation*} \label{notes:eq:30} \frac{\partial F_1}{\partial u_i} = \sum_j u_j a(\phi_i,\phi_k,\phi_j) + \sum_i u_i a(\phi_i,\phi_k,\phi_j) + d(\phi_i,\phi_k) \end{equation*}
\begin{equation*} \label{notes:eq:30} \frac{\partial F_1}{\partial p_i} = -b(\phi_k,\psi_l) \end{equation*}
\begin{equation*} \label{notes:eq:30} \frac{\partial F_1}{\partial \theta_i} = c(\rho_i,\rho_k) \end{equation*}
Consider the second equation, only the derivative with respect to \(u_i\) is non zero.
\begin{equation*} \label{notes:eq:31} \frac{\partial F_2}{\partial u_i} = b(\phi_i,\psi_l) \end{equation*}
Finally the third component
\begin{equation*} \label{notes:eq:33} \frac{\partial F_3}{\partial u_i} = \sum_j \theta_j e(\phi_i,\rho_j,\rho_m) \end{equation*}
\begin{equation*} \label{notes:eq:34} \frac{\partial F_3}{\partial p_i} = 0 \end{equation*}
\begin{equation*} \label{notes:eq:35} \frac{\partial F_3}{\partial \theta_i} = \sum_j u_j e(\phi_j,\rho_i,\rho_m) + f(\rho_i,\rho_m) \end{equation*}
\begin{equation*} \label{notes:eq:35} J_F = \begin{pmatrix} \frac{\partial F_1}{\partial u_i} & \frac{\partial F_1}{\partial p_i} & \frac{\partial F_1}{\partial \theta_i} \\ {\frac{\partial F_2}{\partial u_i}} & {\frac{\partial F_2}{\partial p_i}}(=0) & {\frac{\partial F_2}{\partial \theta_i}}(=0) \\ \frac{\partial F_3}{\partial u_i} & {\frac{\partial F_3}{\partial p_i}}(=0) & \frac{\partial F_3}{\partial \theta_i} \end{pmatrix} \end{equation*}
In order to implement \(J_F\) and solve \((3)\), \(J_F\) can be expressed as the matrix associated with the discretisation of\begin{equation*} \label{notes:eq:37} \begin{array}{rl} a(\mathbf{u},\mathbf{v},\beta_1) + a(\beta_1, \mathbf{v}, \mathbf{u})+d(\mathbf{u},\mathbf{v})-b(\mathbf{v},p)+c(\theta,\mathbf{v}) &= 0\\ b(\mathbf{u},q)&=0\\ e(\beta_1,\theta,\chi)+f(\theta,\chi)+e(\mathbf{u},\beta_2,\chi)&=0\\ \end{array} \end{equation*}
where \(\beta_1 = u^{(n)}\), \(\beta_2=\theta^{(n)}\) are known from the previous Newton iteration, indeed \(J_F\) is actually evaluated in \(X^{(n)}\).Now we use the Feel++ non linear framework in order to implement our Newton scheme. We need to define two new functions in our application
updateJacobian(X,J)
which takes as input X
\(=X^{(n)}\) and returns the matrix J=
\(J_F(X^{(n)})\) updateResidual(X,R)
which takes as input X
\(=X^{(n)}\) and returns the vector R=
\(F(X^{(n)})\)Here is a snippet of code that implements the nonlinear framework.
see bratu.cpp
or nonlinearpow.cpp
for example.