Recall of the fundamental equations of fluid dynamics

This section recalls the fundamental equations of fluid dynamics. It introduces the main notations and the fundamental partial differential equations without simplifying physical hypotheses. The derivations can be found in the classical references. Several mathematically equivalent forms exist in the literature: local form, local form in vector notation, local conservative form, local conservative form in vector notation, macroscopic forms, etc. Here, the local conservative form with vector notations is chosen.

Conservation of mass

\[\frac{\partial\rho}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\mathbf{u})=0 \]

where \(\rho\equiv\rho(\mathbf{x},\,t)\) is the density, with \(\mathbf{x}\) the position and \(t\) the time, and \(\mathbf{u}\equiv\mathbf{u}(\mathbf{x},\,t)\) the velocity. Other equivalent forms of this equation can be obtained by applying the vector identity \(\boldsymbol{\nabla}\cdot(\rho\mathbf{u})=\rho\boldsymbol{\nabla}\cdot\mathbf{u}+\mathbf{u}\cdot\boldsymbol{\nabla}\rho\) and by introducing the material derivative \(d\rho/dt=\partial\rho/\partial t+\mathbf{u}\cdot\boldsymbol{\nabla}\rho\).

Conservation of momentum

The equation of conservation of momentum expresses the fundamental principle of dynamics, which states that the variation of momentum inside a control volume is equal to the sum of all the external forces applied to it. The forces acting on the elementary volume can be separated into volume forces and surface forces. The latter are expressed as a stress vector acting on a surface, and this stress vector is itself expressed as the dot product of a stress tensor \(\mathbf{T}\) (with components \(T_{ij}\)) and the unit normal vector to the surface \(\mathbf{n}\). The total stress is decomposed into two parts: \(T_{ij}=-p\delta_{ij}+\tau_{ij}\). The first is the stress tensor associated with the pressure, \(-p\delta_{ij}\), where \(p\equiv p(\mathbf{x},\,t)\) is the pressure and \(\delta_{ij}\) is the Kronecker symbol, which is one if \(i=j\) and zero otherwise. The second, denoted \(\tau_{ij}\), is associated with the viscous stresses. The pressure acts isotropically and its value depends on the thermodynamic state of the fluid. The viscous stresses are related to the state of deformation of the fluid.

As for the conservation of mass, several mathematically equivalent forms can be deduced. Its local conservative form is:

\[\frac{\partial(\rho\mathbf{u})}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\mathbf{u}\mathbf{u})=-\boldsymbol{\nabla}p+\boldsymbol{\nabla}\cdot\boldsymbol{\tau}+\mathbf{F}_{v} \]

In this momentum equation, the left-hand side represents the amount of acceleration per unit volume. The terms of the right-hand side represent respectively (i) the forces associated with the pressure per unit volume, (ii) the viscous stresses per unit volume and (iii) the external force per unit volume. When only gravity is considered, it is expressed as \(\mathbf{F}_{v}=\rho\mathbf{g}\).

Form of the viscous stress tensor

The viscous stress tensor \(\boldsymbol{\tau}\) is generally expressed as a function of the deformation rates in the flow. The definitions of the strain-rate and rotation-rate tensors, from which the viscous stress tensor will be expressed, are recalled below.

Recall of the strain-rate tensor

The velocity increment of two fluid particles positioned respectively at \(\mathbf{r}\) and \(\mathbf{r}+d\mathbf{r}\), with velocities \(\mathbf{u}\) and \(\mathbf{u}+d\mathbf{u}\), is expressed as \(du_{i}=\sum_{j=1}^{3}(\partial u_{i}/\partial x_{j})dx_{j}\) to first order with respect to the components \(dx_{j}\) (for \(j=1,2,3\)). In this expression, the quantities \(G_{ij}=\partial u_{i}/\partial x_{j}\) are the elements of a rank-two tensor, the deformation-rate tensor of the fluid (or velocity gradients). In three dimensions, it is written as a \(3\times3\) matrix that can be decomposed into a symmetric part and an antisymmetric part:

\[G_{ij}=\frac{\partial u_{i}}{\partial x_{j}}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)+\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right) \]

The first term is the strain-rate tensor:

\[S_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right) \]

and it is symmetric ( \(S_{ij}=S_{ji}\)). The second term is the rotation-rate tensor:

\[\Omega_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right) \]

and this tensor is antisymmetric ( \(\Omega_{ij}=-\Omega_{ji}\)).

Form of the viscous stress tensor for a Newtonian fluid

When the fluids are Newtonian, the stress–strain relation is linear and isotropic. The general relation is written:

\[\tau_{ij}=\eta\left(2S_{ij}-\frac{2}{3}S_{kk}\delta_{ij}\right)+\zeta S_{kk}\delta_{ij} \]

which introduces two viscosities, the dynamic viscosity \(\eta\equiv\eta(\mathbf{x},\,t)\) and the volume viscosity \(\zeta\) (or second viscosity). The first term corresponds to a deformation without change of volume, while the second term corresponds to an isotropic dilatation. In most applications, the volume viscosity is not taken into account ( \(\zeta=0\)) and the stress tensor is written \(\tau_{ij}=2\eta S_{ij}-(2/3)\eta S_{kk}\delta_{ij}\), that is, using the strain-rate tensor relation:

\[\tau_{ij}=\eta\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)-\frac{2}{3}\eta\left(\frac{\partial u_{k}}{\partial x_{k}}\right)\delta_{ij} \]

or, in vector notation:

\[\boldsymbol{\tau}=\eta(\boldsymbol{\nabla}\mathbf{u}+\boldsymbol{\nabla}^{T}\mathbf{u})-\frac{2}{3}\eta(\boldsymbol{\nabla}\cdot\mathbf{u})\mathbf{I} \]

where \(\mathbf{I}\) is the unit diagonal matrix.

Newtonian fluid of constant viscosity

When the viscosity of the fluid \(\eta\) is constant (i.e. \(\eta(\mathbf{x},\,t)=\eta_{0}=\text{Cte}\)), the viscous stress term \(\boldsymbol{\nabla}\cdot\boldsymbol{\tau}\) in the momentum equation is expressed as:

\[\begin{aligned} \frac{\partial\tau_{ij}}{\partial x_{j}} & = \eta_{0}\left[\frac{\partial^{2}u_{i}}{\partial x_{j}\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left(\frac{\partial u_{j}}{\partial x_{i}}\right)\right]-\frac{2}{3}\eta_{0}\frac{\partial}{\partial x_{j}}(\boldsymbol{\nabla}\cdot\mathbf{u})\delta_{ij}\\ & = \eta_{0}\left[\frac{\partial^{2}u_{i}}{\partial x_{j}\partial x_{j}}+\frac{\partial}{\partial x_{i}}(\boldsymbol{\nabla}\cdot\mathbf{u})\right]-\frac{2}{3}\eta_{0}\frac{\partial}{\partial x_{i}}(\boldsymbol{\nabla}\cdot\mathbf{u})\\ & = \eta_{0}\boldsymbol{\nabla}^{2}u_{i}+\frac{\eta_{0}}{3}\frac{\partial}{\partial x_{i}}(\boldsymbol{\nabla}\cdot\mathbf{u}) \end{aligned} \]

That is, in vector notation:

\[\boldsymbol{\nabla}\cdot\boldsymbol{\tau}=\eta_{0}\boldsymbol{\nabla}^{2}\mathbf{u}+\frac{\eta_{0}}{3}\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{u}) \]

Summary

General case for a Newtonian fluid

The conservation equations of mass and momentum are written in the general form:

\[\begin{aligned} \frac{\partial\rho}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\mathbf{u}) & = 0\\ \frac{\partial(\rho\mathbf{u})}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\mathbf{u}\mathbf{u}) & = -\boldsymbol{\nabla}p+\boldsymbol{\nabla}\cdot\boldsymbol{\tau}+\rho\mathbf{F}_{v} \end{aligned} \]

where \(\rho\equiv\rho(\mathbf{x},\,t)\) and \(\mathbf{u}\equiv\mathbf{u}(\mathbf{x},\,t)\) are the unknowns of the system of equations. The equation of state for the pressure and the Newtonian-fluid hypothesis for the stress tensor allow the system to be closed. The equation of state for the pressure and the "low Mach" type models, which separate the pressure into a thermodynamic pressure and a hydrodynamic pressure, are presented in the following chapter. For a Newtonian fluid, the stress tensor \(\boldsymbol{\tau}\) is related to the deformation tensor \(\mathbf{D}\) by the relation:

\[\boldsymbol{\tau}=\eta(\boldsymbol{\nabla}\mathbf{u}+\boldsymbol{\nabla}^{T}\mathbf{u})-\frac{2}{3}\eta(\boldsymbol{\nabla}\cdot\mathbf{u})\mathbf{I} \]

where the dynamic viscosity is denoted \(\eta\equiv\eta(\mathbf{x},\,t)\) and the volume viscosity \(\zeta\) has been neglected.

Particular case of a constant viscosity

When the viscosity is considered constant (i.e. \(\eta(\mathbf{x},\,t)=\eta_{0}=\text{Cte}\)), the divergence term of the stress tensor \(\boldsymbol{\nabla}\cdot\boldsymbol{\tau}\) simplifies and the system of equations becomes:

\[\begin{aligned} \frac{\partial\rho}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\mathbf{u}) & = 0\\ \frac{\partial(\rho\mathbf{u})}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\mathbf{u}\mathbf{u}) & = -\boldsymbol{\nabla}p+\eta_{0}\boldsymbol{\nabla}^{2}\mathbf{u}+\frac{\eta_{0}}{3}\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{u})+\rho\mathbf{F}_{v} \end{aligned} \]

Particular case of incompressible flows with variable viscosity

When the fluid is considered incompressible (i.e. \(\rho(\mathbf{x},\,t)=\rho_{0}=\text{Cte}\)), the conservation of mass becomes \(\boldsymbol{\nabla}\cdot\mathbf{u}=0\) (since \(\partial\rho_{0}/\partial t=0\) and \(\boldsymbol{\nabla}\rho_{0}=\mathbf{0}\)), the nonlinear term is written \(\boldsymbol{\nabla}\cdot(\rho_{0}\mathbf{u}\mathbf{u})=\rho_{0}\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{u}\) and the viscous stress tensor also simplifies to \(\boldsymbol{\tau}=\eta(\boldsymbol{\nabla}\mathbf{u}+\boldsymbol{\nabla}^{T}\mathbf{u})\).

\[\begin{aligned} \boldsymbol{\nabla}\cdot\mathbf{u} & =0,\\ \rho_{0}\frac{\partial\mathbf{u}}{\partial t}+\rho_{0}\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{u} & =-\boldsymbol{\nabla}p+\boldsymbol{\nabla}\cdot\left[\eta(\boldsymbol{\nabla}\mathbf{u}+\boldsymbol{\nabla}^{T}\mathbf{u})\right]+\rho_{0}\mathbf{F}_{v} \end{aligned} \]

In this formulation, the dynamic viscosity \(\eta\) remains inside the bracketed term because it may depend on the position, as in turbulence models.

Particular case of incompressible flows with constant viscosity

Finally, when the fluid is considered incompressible and of constant dynamic viscosity \(\eta_{0}\), the viscous stress tensor simplifies once again to \(\boldsymbol{\nabla}\cdot\boldsymbol{\tau}=\eta_{0}\boldsymbol{\nabla}^{2}\mathbf{u}\) and the system of equations becomes:

\[\begin{aligned} \boldsymbol{\nabla}\cdot\mathbf{u} & =0,\\ \rho_{0}\frac{\partial\mathbf{u}}{\partial t}+\rho_{0}\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{u} & =-\boldsymbol{\nabla}p+\eta_{0}\boldsymbol{\nabla}^{2}\mathbf{u}+\rho_{0}\mathbf{F}_{v} \end{aligned} \]

Conservation of energy

Several forms of the energy balance equation are possible, depending on whether one considers the conservation of total energy, internal energy, total enthalpy or enthalpy. Formulations can also be deduced for the temperature or even the entropy. In the following, the presentation is restricted to the conservation of internal energy, which will be written in an equivalent form on the temperature equation.

The balance equation for the internal energy \(e\) is written [14] (p. 126):

\[\frac{\partial\rho e}{\partial t}+\boldsymbol{\nabla}\cdot(\rho e\mathbf{u})=-\boldsymbol{\nabla}\cdot\mathbf{q}-p\boldsymbol{\nabla}\cdot\mathbf{u}+\boldsymbol{\tau}:\boldsymbol{\nabla}\mathbf{u} \]

In this equation, the left-hand side represents the rate of variation of the internal energy per unit volume. On the right-hand side, the first term represents the heat flux per unit volume; the second term represents the power of the pressure forces per unit volume; and the last term represents the power of the viscous forces per unit volume. This last term is the viscous dissipation function, which is always positive or zero. Thus the viscous forces always lead to an increase of the internal energy of the fluid, and therefore of its temperature.

The conservation of internal energy can be reformulated into an equation on the temperature \(T\). This equation takes two different forms depending on whether one uses the specific heat at constant volume \(C_{v}\) or at constant pressure \(C_{p}\). Formulated in \(C_{v}\), it is written:

\[\frac{\partial(\rho C_{v}T)}{\partial t}+\boldsymbol{\nabla}\cdot(\rho C_{v}T\mathbf{u})=-\boldsymbol{\nabla}\cdot\mathbf{q}-T\left(\frac{\partial p}{\partial T}\right)_{\rho}\boldsymbol{\nabla}\cdot\mathbf{u}+\boldsymbol{\tau}:\boldsymbol{\nabla}\mathbf{u}+\rho T\frac{dC_{v}}{dt} \]

Formulated in \(C_{p}\), it is written:

\[\frac{\partial(\rho C_{p}T)}{\partial t}+\boldsymbol{\nabla}\cdot(\rho C_{p}T\mathbf{u})=-\boldsymbol{\nabla}\cdot\mathbf{q}-\left(\frac{\partial\ln\rho}{\partial\ln T}\right)\frac{dp}{dt}+\boldsymbol{\tau}:\boldsymbol{\nabla}\mathbf{u}+\rho T\frac{dC_{p}}{dt} \]

Note: Remark: when the gas is considered ideal, i.e. \(\rho=p_{th}/(RT)\) where \(R\) is the ideal gas constant and \(p_{th}\) is the thermodynamic pressure, the coefficient \((\partial\ln\rho/\partial\ln T)\) is equal to \(-1\). Assuming that the heat flux is defined by Fourier's law \(\mathbf{q}=-\lambda\boldsymbol{\nabla}T\) where \(\lambda\) is the thermal conductivity, the \(C_{p}\) equation becomes:
\[\frac{\partial(\rho C_{p}T)}{\partial t}+\boldsymbol{\nabla}\cdot(\rho C_{p}T\mathbf{u})=\boldsymbol{\nabla}\cdot(\lambda\boldsymbol{\nabla}T)+\frac{dp_{th}}{dt}+\boldsymbol{\tau}:\boldsymbol{\nabla}\mathbf{u}+\rho T\frac{dC_{p}}{dt} \]

The introduction of the thermodynamic pressure \(p_{th}\) will be useful for the "low Mach" models.

Boussinesq approximation in the incompressible case

For incompressible flows, when the density is assumed constant \(\rho=\rho_{0}=\text{Cte}\), the density is neither a function of the temperature nor of the composition of the fluid (for miscible mixtures). In this case, the buoyancy effects are only taken into account through the gravitational forces. This simplification is known as the "Boussinesq approximation" and is valid by considering that the density variation \(\Delta\rho\ll\rho_{0}\). In this case, the force term in the incompressible equations is written:

\[\mathbf{F}_{v}=-\mathbf{g}\beta_{T}(T-T_{0}) \]

where \(\beta_{T}\) is the coefficient of thermal expansion and \(T_{0}\) a reference temperature. In this relation, the negative sign indicates that if the temperature difference is positive \(\Delta T=T-T_{0}>0\) (i.e. near the hot wall in natural convection), then the force is directed in the opposite direction to gravity \(\mathbf{g}\).

Numerical implementation

The numerical methods implemented in TrioCFD for the incompressible model defined by the equations above (time scheme, VEF space scheme, and the P1NC/P0 finite-element-volume matrices) are described in Discretization in TrioCFD.