\begin{alignat}{2}
\rho \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}\right) - \nabla \cdot (\mu \nabla \mathbf{u}) + \nabla p &= - \rho_0 \beta(T-T_\text{ref}) \mathbf{g} & \quad\text{in} \; \Omega_F & \quad(2)\\
\nabla\cdot \mathbf{u} &= 0 & \quad\text{in} \; \Omega_F & \quad(3)\\
\mathbf{u} &= 0 & \quad\text{on} \; \partial \Omega_F & \quad\text{(Dirichlet)}
\end{alignat}
\(\quad\)We consider \(\mathbf{v} \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{v}\) and we integrate on \(\Omega_F\).
\[\rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} -\int_{\Omega_F} (\nabla \cdot (\mu \nabla \mathbf{u})) \cdot \mathbf{v} + \int_{\Omega_F} \nabla p \cdot \mathbf{v} = -\int_{\Omega_F} \rho_0 \beta(T-T_\text{ref}) \mathbf{g} \cdot \mathbf{v}\]
\(\quad\)We can also assume that \(\mu\) is constant. We have
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{v}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} -\mu \int_{\Omega_F} \Delta \mathbf{u} \cdot \mathbf{v} + \int_{\Omega_F} \nabla p \cdot \mathbf{v} = -\int_{\Omega_F} \rho_0 \beta(T-T_\text{ref}) \mathbf{g} \cdot \mathbf{v}\]
\(\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{v}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} +\mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{v} - \int_{\partial \Omega_F}\frac{\partial \mathbf{u}}{\partial n } \cdot \mathbf{v} + \int_{\Omega_F} \nabla p \cdot \mathbf{v} =- \int_{\Omega_F} \rho_0 \beta(T-T_\text{ref}) \mathbf{g} \cdot \mathbf{v}\]
\[ \rho \int_{\Omega_F}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\mathbf{v}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} +\mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{v} - \int_{\partial \Omega_F}\frac{\partial \mathbf{u}}{\partial n } \cdot \mathbf{v} - \int_{\Omega_F} p \cdot \nabla \mathbf{v} + \int_{\partial \Omega_F} \frac{\partial p}{\partial n} \cdot \mathbf{v} = - \int_{\Omega_F} \rho_0 \beta(T-T_\text{ref})\mathbf{g} \cdot \mathbf{v}\]
\(\quad\)As \(\mathbf{v}\) 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{v}} + \rho\int_{\Omega_F} (\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} +\mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{v} - \int_{\Omega_F} p \cdot \nabla \mathbf{v} =- \int_{\Omega_F} \rho_0 \beta(T-T_\text{ref})\mathbf{g} \cdot \mathbf{v}\]
\(\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{ and } 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{\mathbf{u}^{k+1}}{dt}\cdot \mathbf{\mathbf{v}} + \rho\int_{\Omega_F} (\mathbf{u}^{k+1} \cdot \nabla \mathbf{u}^{k+1}) \cdot \mathbf{v} +\mu \int_{\Omega_F} \nabla \mathbf{u}^{k+1} : \nabla \mathbf{v} - \int_{\Omega_F} p \cdot \nabla \mathbf{v} = \rho \int_{\Omega_F}\frac{\mathbf{u}^{k}}{dt}\cdot \mathbf{\mathbf{v}} - \int_{\Omega_F} \rho_0 \beta(T-T_\text{ref})\mathbf{g} \cdot \mathbf{v}\]
\(\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{\mathbf{u}^{k+1}}{dt}\cdot \mathbf{\mathbf{v}} + \rho\int_{\Omega_F} (\mathbf{u}^{k+1} \cdot \nabla \mathbf{u}^{k+1}) \cdot \mathbf{v} +\mu \int_{\Omega_F} \nabla \mathbf{u}^{k+1} : \nabla \mathbf{v} = \rho \int_{\Omega_F}\frac{\mathbf{u}^{k}}{dt}\cdot \mathbf{\mathbf{v}} - \int_{\Omega_F} \rho_0 \beta(T-T_\text{ref})\mathbf{g} \cdot \mathbf{v}\]
\(\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 (3). Let \(q \in L^2 (\Omega)\), we multiply (3) by \(q\):
\[\int_{\Omega_F} \nabla \mathbf{u} \cdot q = 0\]
Let \(\phi\in H^1(\Omega)\) be a test function associated to the temperature space. We multiply the equation (1) by \(\phi\) and integrate over \(\Omega\):
\[\begin{equation}
\int_{\Omega_H}\rho C_p\left(\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T\right)\phi - k\nabla^2 T\phi = \int_{\Omega_H}Q\phi.
\end{equation}\]
Doing another integration by parts, we obtain:
\[\begin{equation}
\int_{\Omega_H} \rho C_p\left(\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T \right)\phi + \int_{\Omega_H} k\nabla T\cdot\nabla\phi = \int_{\partial\Omega_H} k\nabla T\cdot\boldsymbol{n}\phi + \int_{\Omega_H}Q\phi.
\end{equation}\]
|
|
The right hand side term is updated according to the boundary conditions of the heat problem. In the following, we will consider homogeneous Neumann boundary conditions: \(k\nabla T\cdot\boldsymbol{n} = 0\) on \(\partial\Omega_H\).
|
As earlier, we use the Euler implicit method for the time term:
\[\frac{\partial T}{\partial t} (t^{k+1}) \approx \frac{T(t^{k+1}) - T(t^k)}{dt} =: \frac{T^{k+1}-T^k}{dt} \quad \forall t^k \in \mathbb{R^+} \text{ and } k \in \mathbb{N}\]
\[\begin{equation}
\int_{\Omega_H} \rho C_p \left(\frac{T^{k+1}}{dt} + \mathbf{u}\cdot\nabla T^{k+1}\right)\phi + \int_{\Omega_H} k\nabla T^{k+1}\cdot\nabla\phi = \int_{\Omega_H}Q\phi + \rho C_p\frac{T^k}{dt} \phi.
\end{equation}\]
\(\quad\)We then pose a trilinear form \(a \colon \mathcal{H}_0^1(\Omega_F)^3 \times \mathcal{H}_0^1(\Omega_F)^3 \times \to \mathcal{H}_0^1(\Omega_F)^3 \mathbb{R}\) defined by:
\[a(\mathbf{u}, \mathbf{u}, \mathbf{v})=\rho\int_{\Omega_F} ( \frac{\mathbf{u}}{dt} + \mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} + \mu \int_{\Omega_F} \nabla \mathbf{u} : \nabla \mathbf{v}\]
and two bilinear forms \(b \colon \mathcal{H}_0^1(\Omega_F)^3 \times L^2_0(\Omega_F) \to \mathbb{R}\) and \(d \colon \mathcal{H}^1(\Omega_F) \times \mathcal{H}_0^1(\Omega_F)^3 \times \to \mathbb{R}\) defined by:
\[\begin{align}
b(\mathbf{v}, p) &= -\int_{\Omega_F} p \cdot \nabla \mathbf{v}
\\
d(T, \mathbf{v}) &= -\int_{\Omega_F} \rho_0 \beta T \mathbf{g} \cdot \mathbf{v}
\end{align}\]
and a tri-linear form \(e \colon \mathcal{H}_0^1(\Omega_F)^3 \times \mathcal{H}^1(\Omega_F) \times \mathcal{H}^1(\Omega_F) \to \mathbb{R}\) and a bilinear form \(f \colon \mathcal{H}^1(\Omega_F) \times \mathcal{H}^1(\Omega_F) \to \mathbb{r}\) defined by:
\[\begin{align}
e(\mathbf{u}, T, \phi) &= \rho C_p\int_{\Omega_H} \left(\frac{T}{dt} + \mathbf{u}\cdot\nabla T\right)\phi
\\
f(T, \phi) &= k\int_{\Omega_H} \nabla T\cdot\nabla\phi
\end{align}\]
\(\quad\)The variational formulation is then written as follows: find \(\mathbf{u}^{k+1} \in \mathcal{H}_0^1(\Omega_F)^3\) and \(p \in L^2(\Omega_F)\) such as:
\[\begin{cases}
a(\mathbf{u}^{k+1},\mathbf{u}^{k+1}, \mathbf{v}) &+ b(\mathbf{v},p) &+ d(T, \mathbf{v}) &= \ell_1( \frac{\mathbf{u}^{k}}{dt} + \rho\beta T_\text{ref}\mathbf{g}, \mathbf{v}) \\
b(\mathbf{u}, q) & & &= 0 \\
e(\mathbf{u}^{k+1}, T, \phi) & &+ f(T, \phi) &= \ell_2(T^k + Q, \phi)
\end{cases}\]
\(\quad\)The problem for the Stokes equation \(a(\mathbf{u}, \mathbf{v}) + b(\mathbf{v}, p)\) is well settled if the 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{such that} \quad \sup_{\mathbf{v} \in \mathcal{H}_0^1(\Omega), \mathbf{v} \neq 0} \frac{b(\mathbf{v}, q)}{\lVert \mathbf{v} \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 that the pressure will be determined numerically during the simulations.
2.2. Discretization
Let \(\mathcal{V}_h\) is the discretized space of the velocity, \(\mathcal{P}_h\) is the discretized space of the pressure and \(\mathcal{T}_h\) is the discretized space of the temperature.
We set \(N_u = \mathrm{dim} (\mathcal {V}_h)\), \(N_p = \mathrm{dim} (\mathcal{P}_h)\) and \(N_T = \mathrm{dim}(\mathcal{T}_h)\). Let \(\{ \boldsymbol{\lambda}_i \}_{i = 1, ..., N_u}\), \(\{ \mu_j \}_{j = 1, ..., N_p}\) and \(\{ \theta_m \}_{ = 1, ..., N_T}\) be the bases of \(\mathcal{V}_h\), \(\mathcal{P}_h\) and \(\mathcal{T}_h\) respectively.
\[\begin{align}
\mathbf{u}_h &= \sum_{i=1}^{N_u} u_i \boldsymbol{\lambda}_i
\\
p_h &= \sum_{j=1}^{N_p} p_j \mu_j
\\
T_h &= \sum_{m=1}^{N_T} T_m \theta_m
\end{align}\]
\(\quad\) Returning everything into the variational formulation, we obtain the discrete formulation.
\[\begin{cases}
a\left(\sum_{i=1}^{N_u} u_i \boldsymbol{\lambda}_i, \mathbf{v}\right) &+ b\left(\mathbf{v}, \sum_{j=1}^{N_p} p_j \mu_j\right) &+ d\left(\sum_{m=1}^{N_T} T_m \theta_m, \mathbf{v}\right) &= \ell_1(\mathbf{v}) \\
b\left(\sum_{i=1}^{N_u} u_i \boldsymbol{\lambda}_i, q\right) & & &= 0 \\
e\left(\sum_{i=1}^{N_u} u_i \boldsymbol{\lambda}_i, \sum_{m=1}^{N_T} T_m \theta_m, \phi\right) & &+ f\left(\sum_{m=1}^{N_T} T_m \theta_m, \phi\right) &= \ell_2(\phi)
\end{cases}\]
\(\quad\)By putting \(\mathbf{v} = \boldsymbol{\lambda}_i, \; \forall i = 1, ..., N_u\) and \(q = \mu_j, \; \forall j = 1, ..., N_p\). We obtain :
\[\begin{cases}
\forall j \in \{1,...,N_u\}, \quad \sum_{i=1}^{N_u} u_i a(\boldsymbol{\lambda}_i, \boldsymbol{\lambda}_j) + \sum_{k=1}^{N_p} p_k b(\boldsymbol{\lambda}_j, \mu_k) = \ell_1(\boldsymbol{\lambda}_j)
\\
\forall k \in\{1,...,N_p\}, \quad \sum_{i=1}^{N_u} u_i b(\boldsymbol{\lambda}_i, \mu_k) = 0
\end{cases}\]
\(\quad\)By putting the following matrices:
\[\begin{equation}
U =(u_i)_{i=1,..,N_u}^T \quad p =(p_i)_{i=1,...,N_p}^T
\\
A = (a(\boldsymbol{\lambda}_i, \boldsymbol{\lambda}_j))_{1 \leq i,j \leq N_u} \quad B= (b(\boldsymbol{\lambda}_i, \mu_j))_{1 \leq i \leq N_u , 1 \leq j \leq N_p}
\\
F = (\ell_1(\boldsymbol{\lambda}_i))_{i=1,...,N_u}
\end{equation}\]
\(\quad\) We obtain the following linear system:
\[\begin{bmatrix}
A & B^T \\
B & 0 \\
\end{bmatrix}
\begin{bmatrix}
U \\
p \\
\end{bmatrix}
=
\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