Convective heat equation
\[\rho C_p \left( \frac{\partial T}{\partial t} + \boldsymbol{u} \cdot \nabla T \right) - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega_H\]
which is completed with boundary conditions and initial value
\[\text{at } t=0, \quad T(x,0) = T_0(x)\]
Equation of air movement (Navier-Stokes)
\[\begin{cases}
\rho (\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}) -\nabla \cdot (\mu \nabla \mathbf{u}) + \nabla P = - \rho_0 \beta(T-T_{ref}) \mathbf{g} & dans \; \Omega_F \quad (1) \\
\nabla \mathbf{u} = 0 & dans \; \Omega_F \quad (2) \\
\mathbf{u}=0 & sur \; \partial \Omega_F \quad \text{(boundary of Dirichlet)}
\end{cases}\]
\(\quad\)The equation (1) is the momentum equation inherited from Newton’s law and (2) is the mass conservation equation for incompressible flows.
\(\quad\)We consider \(\mathbf{\phi} \in \mathcal{H}_0^1(\Omega)^d\) a test function with compact support in the Sobolev space in dimension d. We multiply our equation by \(\mathbf{\phi}\) and we integrate on \(\Omega_F\).
\[\rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{\phi} -\int_{\Omega_F} (\nabla \cdot (\mu \nabla \mathbf{u})) \cdot \mathbf{\phi} + \int_{\Omega_F} \nabla P \cdot \mathbf{\phi} = -\int_{\Omega_F} \rho_0 \beta(T-T_{ref}) \mathbf{g} \cdot \mathbf{\phi}\]
\(\quad\)We can also assume that \(\mu\) is constant. We have
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{\phi} -\mu \int_{\Omega_F} \Delta \mathbf{u} \cdot \mathbf{\phi} + \int_{\Omega_F} \nabla P \cdot \mathbf{\phi} = -\int_{\Omega_F} \rho_0 \beta(T-T_{ref}) \mathbf{g} \cdot \mathbf{\phi}\]
\(\quad\)Using successively the formulas of Green on the term in \(\Delta \mathbf{u}\), then the term in pressure, we obtain:
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{\phi} +\mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{\phi} - \int_{\partial \Omega_F}\frac{\partial \mathbf{u}}{\partial n } \cdot \mathbf{\phi} + \int_{\Omega_F} \nabla P \cdot \mathbf{\phi} =- \int_{\Omega_F} \rho_0 \beta(T-T_{ref}) \mathbf{g} \cdot \mathbf{\phi}\]
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{\phi} +\mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{\phi} - \int_{\partial \Omega_F}\frac{\partial \mathbf{u}}{\partial n } \cdot \mathbf{\phi} - \int_{\Omega_F} P \cdot \nabla \mathbf{\phi} + \int_{\partial \Omega_F} \frac{\partial P}{\partial n} \cdot \mathbf{\phi} = - \int_{\Omega_F} \rho_0 \beta(T-T_{ref})\mathbf{g} \cdot \mathbf{\phi}\]
\(\quad\)As \(\phi\) is compact support, the terms on the edges vanish.
We will then obtain:
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{\phi} +\mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{\phi} - \int_{\Omega_F} P \cdot \nabla \mathbf{\phi} =- \int_{\Omega_F} \rho_0 \beta(T-T_{ref})\mathbf{g} \cdot \mathbf{\phi} \quad (3)\]
\(\quad\)Using the implicit Euler method for the time term:
\[\frac{\partial \mathbf{u}}{\partial t} (t^{ k+1}) \approx \frac{ \mathbf{u} (t^{ k+1}) - \mathbf{u}(t^k)}{ dt} \quad \forall t^k \in \mathbb{ R^+} \text{ et } k \in \mathbb{N}\]
\(\quad\)Denoting \(\mathbf{u}^k = \mathbf{u}(t^k)\), we write the formula in \(t^{ k+1}\), we obtain:
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}^{k+1}}{dt}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u}^{k+1} \cdot \nabla \mathbf{u}^{k+1}) \cdot \mathbf{\phi} +\mu \int_{\Omega_F} \nabla \mathbf{u}^{k+1} : \nabla \mathbf{\phi} - \int_{\Omega_F} P \cdot \nabla \mathbf{\phi} = \rho \int_{\Omega_F}\frac{\partial \mathbf{u}^{k}}{dt}\cdot \mathbf{\mathbf{\phi}} - \int_{\Omega_F} \rho_0 \beta(T-T_{ref})\mathbf{g} \cdot \mathbf{\phi}\]
\(\quad\)If you restrict the space of the test functions to the next space \(\mathcal{V}(\Omega)=\{(v \in \mathcal{H}_0^1(\Omega))^3 | \nabla \cdot v=0 \}\), we obtain the following weak wording:
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}^{k+1}}{dt}\cdot \mathbf{\mathbf{\phi}} + \rho\int_{\Omega_F} (\mathbf{u}^{k+1} \cdot \nabla \mathbf{u}^{k+1}) \cdot \mathbf{\phi} +\mu \int_{\Omega_F} \nabla \mathbf{u}^{k+1} : \nabla \mathbf{\phi} = \rho \int_{\Omega_F}\frac{\partial \mathbf{u}^{k}}{dt}\cdot \mathbf{\mathbf{\phi}} - \int_{\Omega_F} \rho_0 \beta(T-T_{ref})\mathbf{g} \cdot \mathbf{\phi}\]
\(\quad\)As the space \(\mathcal{V}\) is difficult to build, so we will use the formulation (3). We now look at what happens with equation (2). Let \(q \in L^2 (\Omega)\), we multiply (2) by \(q\):
\[\int_{\Omega_F} \nabla \mathbf{u} \cdot q = 0\]
\(\quad\)We then pose 2 bilinear forms \(a : \mathcal{H}_0^1(\Omega_F)^3 \times \mathcal{H}_0^1(\Omega_F)^3 \rightarrow \mathbb{R}\) and \(b : \mathcal{H}_0^1(\Omega_F)^3 \times L^2(\Omega_F) \rightarrow \mathbb{R}\):
\[\begin{align}
a(u,\phi)= \rho \int_{\Omega_F}\frac{\partial u}{dt}\cdot \mathbf{\mathbf{\phi}} + \int_{\Omega_F} \nabla u : \nabla \mathbf{\phi}
\\
b(\phi,p) =- \int_{\Omega_F} p \cdot \nabla \mathbf{\phi}
\end{align}\]
\(\quad\)And a trilinear form \(c: \mathcal{H}_0^1(\Omega_F)^3 \times \mathcal{H}_0^1(\Omega_F)^3 \times \mathcal{H}_0^1(\Omega_F)^3 \rightarrow \mathbb{R}\) :
\[c(u,u,\phi)=\int_{\Omega_F} (u \cdot \nabla u) \cdot \mathbf{\phi}\]
\(\quad\)The variational formulation is then written as follows:
\[\begin{cases}
c(\mathbf{u}^{k+1},\mathbf{u}^{k+1},\mathbf{\phi}) + a(\mathbf{u}^{k+1},\mathbf{\phi}) + b(\mathbf{\phi},P)= l( \mathbf{u}^{k}{dt} - \rho_0 \beta(T-T_{ref}) \mathbf{g},\mathbf{\phi}) \\
b(\mathbf{u},q)=0
\end{cases}\]
\(\quad\)The problem for the Stokes equation \(a(u,\phi) + b(\phi,p)\) is well settled, if the stem form: [a] is coercive on \(\mathcal {H}_0^1 (\Omega)\) and the form b satisfies the condition 'inf-sup', that is to say:
\[\exists \beta >0 \; \text{tel que} \quad \sup_{\phi \in \mathcal{H}_0^1(\Omega), \phi \neq 0} \frac{b(\phi,q)}{\lVert \phi \rVert_{\mathcal{H}^1}}\geq \beta \lVert q \rVert_{L^2} \quad \forall q \in L^2(\Omega)\]
\(\quad\) We must have at least these verified hypotheses to have a solution for the Navier-Stokes equation. In our case, the pressure is then defined to a constant, to have a single pressure we can take it in the space \(L^2(\Omega)\) to average zero. But nothing assures that it’s the right space, so the pressure will be determined numerically during the simulations.
\(\quad\) We first consider the Stokes problem for discretization.
Let \(\mathcal{V}_h\) is the discretized space of the velocity and \(\mathcal{P}_h\) is the discretized space of the pressure.
We pose \(N_u = dim (\mathcal {V}_h)\) and \(N_p = dim (\mathcal{P}_h)\). Let \(\{ \lambda_i \}_{i = 1, ..., N_u}\) is a base of \(\mathcal{V}_h\) and \(\{ \mu_i \}_{ i = 1, ..., N_p}\) is a base of \(\mathcal{P}_h\).
\[\begin{align}
u_h = \sum_{i=1}^{N_u} u_i \lambda_i
\\
p_h = \sum_{j=1}^{N_p} p_j \mu_j
\end{align}\]
\(\quad\) Returning everything into the variational formulation, we obtain the formulation discrete.
\[\begin{cases}
a(\sum_{i=1}^{N_u} u_i \lambda_i , \phi) + b(\phi , \sum_{j=1}^{N_p} p_j \mu_j ) = l(\phi)
\\
b(\sum_{i=1}^{N_u} u_i \lambda_i , q) = 0
\end{cases}\]
\(\quad\)By putting \(\phi = \lambda_i, \; \ forall i = 1, ..., N_u\) and \(q = \mu_j, \; \forall j = 1, ..., N_p\). We obtain :
\[\begin{cases}
j=1,...,N_u, \quad \sum_{i=1}^{N_u} u_i a(\lambda_i, \lambda_j) + \sum_{k=1}^{N_p} p_k b(\lambda_j, \mu_k) = l(\phi_j)
\\
k=1,...,N_p, \quad \sum_{i=1}^{N_u} u_i b(\lambda_i, \mu_k) = 0
\end{cases}\]
\(\quad\)By putting the following matrices:
\[\begin{align}
U =(u_i)_{i=1,..,N_u}^T \quad P =(p_i)_{i=1,...,N_p}^T
\\
A = (a(\lambda_i, \lambda_j))_{1 \leq i,j \leq N_u} \quad B= (b(\lambda_i, \mu_j))_{1 \leq i \leq N_u , 1 \leq j \leq N_p}
\\
F = (l(\phi_i))_{i=1,...,N_u}
\end{align}\]
\(\quad\) We obtain the following linear system:
\[\begin{pmatrix}
A & B^T \\
B & 0 \\
\end{pmatrix}
\begin{pmatrix}
U \\
P \\
\end{pmatrix}
=
\begin{pmatrix}
F \\
0 \\
\end{pmatrix}\]
\(\quad\) Since the Navier-Stokes problem is not linear, we will use the Newton’s or Picard’s method for solving.
.pdf