The zero-equation turbulence model

We obtain a non-linear system of three coupled equations

This system is supplemented with boundary conditions and initial values for the different fields of interest. These closure conditions will be detailed for each model and test case studied but we can already enumerate some usual boundary conditions and notations.

We consider the problem on a domain \(\Omega \in \mathbb{R}^d\), \(d = 2\), \(3\), with a sufficiently smooth boundary, \(\partial \Omega = \Gamma\). This boundary is partitioned (in two different ways) to settle the different boundary conditions

where the exponents \(F\) and \(T\) stand for fluid or temperature and the indices \(D\), \(N\) and \(R\) stand for Dirichlet, Neumann and Robin conditions. On these subsets of \(\Gamma\), the boundary conditions can be described as

where \(u_D\), \(T_D\), \(\phi_{\Gamma}\) and \(\kappa_{\Gamma}\) are at least \(\mathcal{L}^2\) functions on \(\Gamma\).

In order to write the variational formulation, we now introduce definite function spaces \(\forall t\).
We have the solution function spaces associated respectively with the speed and pressure:

as well as the space of the fluid velocity test functions:

Then, we multiply equation \((2.1)\) (respectively \((2.2)\)) with a test function \(v \in W\) (respectively \(q \in Q\)). We then integrate each of the equations on the domain \(\Omega\), and we perform an integration by part of the term containing the stress tensor \(\sigma\), which gives us:

We now introduce the solution space assiociated with the temperature:

as well as the space of the temperature test functions:

Then, we multiply equation \((2.3)\) with a test function \(r \in Y\). We then integrate each of the equations on the domain \(\Omega\), and we perform an integration by part of the term containing the stress tensor \(\sigma\), which gives us:

Finally, we must take into account the boundary conditions. Thus, we obtain the incompressible Navier-Stokes variational formulation for the problem we have considered:
Find \((u, p, T) \in V \times Q \times X\) such that \(\forall (v, q, r) \in W \times Q \times Y\) we have:

We will start by introducing some notations associated with the geometrical discretization. Let \(\mathcal{T}_{\delta}\) be the mesh resulting from the discretization of the field \(\Omega\) and thus forming the computational field of reference \(\Omega_{\delta}\). We also set \(\Gamma_{D,\delta}^F, \Gamma_{N,\delta}^F, \Gamma_{D,\delta}^T, \Gamma_{N,\delta}^T\) and \(\Gamma_{R,\delta}^T\) the discrete borders associated respectively with \(\Gamma_D^F, \Gamma_N^F, \Gamma_D^T, \Gamma_N^T\) and \(\Gamma_R^T\).
We sample the time dimension in \(N_T\) time steps, noted \(\{t_n\}_{1\leqslant n\leqslant N_T}\). Let \(t_i\) and \(t_f\) be the initial and final times respectively. We then have \(t_i\leqslant t_n\leqslant t_f,\) \(\forall n\in \{0,\dots,N_T\}\) and \(t_a<t_b\) if and only if \(a<b\). We note \(u_{\delta}^n\), \(p_{\delta}^n\) and \(T_{\delta}^n\) the discrete speed, pressure and temperature fields at time \(t_n\).

We start by describing the approximation spaces that we will use to write the discrete variational formulation. Let:

be the discrete spaces associated respectively with the velocity, the pressure and the temperature. The polynomial degrees \(M\) and \(N\) cannot be chosen freely, because the discrete problem obtained from the continuous variational formulation will not necessarily be well posed according to this choice. This discrete problem is of the saddle point type, the approximation spaces must then check a compatibility condition to ensure the existence and uniqueness of the problem. We can find for example the mathematical details in [10]. In the literature, this condition is often referred to as the Babuška – Brezzi condition or the Ladyshenskaya – Babuška – Brezzi condition. This condition states that the discrete problem is well posed if and only if the spaces \(V_{\delta}\) and \(Q_{\delta}\) are such that there exists a constant \(\beta_{\delta}\) such that:

If this constraint is not checked, the numerical solution can reveal instabilities. This is the case when \(M = N\). Solutions exist to solve this problem. They consist in adding stabilization terms to the equations. Several examples of stabilization exist, we can cite for example the PSPG [11], SUPG [12], GLS [13] and CIP [14] methods. We choose to use one of the more common finite element choices, the Taylor-Hood elements [9, 14]. The polynomial orders are then connected by \(N = M - 1\).

Since we have non-homogeneous Dirichlet conditions, we also need to introduce the discrete spaces of test functions associated respectively with the speed and temperature:

The discrete variational formulation results in the following problem:
Find \((u_{\delta}, p_{\delta}, T_{\delta}) \in V_{\delta} \times Q_\delta \times X_\delta\) such that \(\forall (v, q, r) \in W_\delta \times Q_\delta \times Y_\delta\) we have:

We will now approach the derivatives in time using implicit schemes called backward differentiation formulation \((\mathrm{BDF})\). These schemes can be written for an arbitrary order \(q\) and they are thus denoted by \(\mathrm{BDF_q}\). We will present \(\mathrm{BDF_q}\) diagrams up to order 4 in the following. Let \(\Delta t\) be the time step assumed to be constant over time, we have \(t_0=t_i\) and \(t_n = t_0 + n\Delta t\).

The following formula describes the approximation of the \(q\)-order time derivative of the speed using the \(\beta_j\) coefficients described in Table 1:

This expression is made up of two terms. The first contains the unknown \(u^{n+1}_\delta\). The second shows the solutions of the previous time steps, that is to say that they are known.
We define:

where \(q^F\) and \(q^T\) are diagram orders for speed and temperature, respectively.
We finish the temporal discretization by adding the initial conditions necessary to be able to correctly construct the temporal derivatives. If we use a \(\mathrm{BDF_1}\) diagram, only one initial condition is necessary \(u^0_\delta = u_0\) in \(\Omega_\delta\), where \(u_0=u(x, t_0)\). We can then calculate the solutions at the following time step \(t_n\) for \(n > 1\). When we use a \(q\) diagram of order greater than 1, we have two possibilities. Either we use \(q\) initial conditions, which represents the mathematically cleanest solution, but this information is not always accessible in practice. Either we gradually increase the order of the diagram, i.e. for the first time step we use a \(\mathrm{BDF_1}\) scheme which was initialized with a single initial condition \(u^0_\delta\), then for the second time step we use \(\mathrm{BDF_2}\) until reaching the desired order. However, this method can have the disadvantage of creating discontinuities at this start if the simulated phenomena are important.

Finally, the discrete variational formulation of our standard problem is given by the following problem:
Find \((u_{\delta}, p_{\delta}, T_{\delta}) \in V_{\delta} \times Q_\delta \times X_\delta\) such that \(\forall (v, q, r) \in W_\delta \times Q_\delta \times Y_\delta\) we have:

In order to simplify and speed up the calculations, we will use an extrapolation [16] of the velocity and separate the equations.
Let \(u^*\) be the extrapolation of \(u_\delta^{n+1}\), the problem becomes:

where \(l\) is the distance to the nearest closure and \(U^*\) is the extrapolated local mean velocity.