Heat Transfer

1. Notations and Units

Notation Quantity Unit

$\rho$

fluid density

$kg \cdot m^{-3}$

$C_p$

heat capacity at constant pressure

$J \cdot kg^{-1} \cdot K^{-1}$

$T$

temperature

$K$

$k$

thermal conductivity

$W \cdot m^{-1} \cdot K^{-1}$

$\boldsymbol{u}$

fluid velocity

$m \cdot s^{-1}$

$\beta$

coefficient of thermal expansion

$K^{-1}$

$\mu$

dynamic viscosity

$Pa \cdot s$

$\mathbf{g}$

gravitational acceleration

$m \cdot s^{-2}$

$\rho_0$

fluid density of air

$kg \cdot m^{-3}$

2. Equations

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 stem: [N_p = dim (\ mathcal {P}_h)]. Let $\{ \lambda_i \}_{i = 1, ..., N_u}$ is a base of stem: [\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.