Staggered Semi-Implicit Hybrid Finite Volume/Finite Element Schemes for Turbulent and Non-Newtonian Flows

Abstract

This paper presents a new family of semi-implicit hybrid finite volume/finite element schemes on edge-based staggered meshes for the numerical solution of the incompressible Reynolds-Averaged Navier–Stokes (RANS) equations in combination with the $k-\epsilon$ turbulence model. The rheology for calculating the laminar viscosity coefficient under consideration in this work is the one of a non-Newtonian Herschel–Bulkley (power-law) fluid with yield stress, which includes the Bingham fluid and classical Newtonian fluids as special cases. For the spatial discretization, we use edge-based staggered unstructured simplex meshes, as well as staggered non-uniform Cartesian grids. In order to get a simple and computationally efficient algorithm, we apply an operator splitting technique, where the hyperbolic convective terms of the RANS equations are discretized explicitly at the aid of a Godunov-type finite volume scheme, while the viscous parabolic terms, the elliptic pressure terms and the stiff algebraic source terms of the $k-\epsilon$ model are discretized implicitly. For the discretization of the elliptic pressure Poisson equation, we use classical conforming ${\mathcal{P}}^1$ and ${\mathcal{Q}}^1$ finite elements on triangles and rectangles, respectively. The implicit discretization of the viscous terms is mandatory for non-Newtonian fluids, since the apparent viscosity can tend to infinity for fluids with yield stress and certain power-law fluids. It is carried out with ${\mathcal{P}}^1$ finite elements on triangular simplex meshes and with finite volumes on rectangles. For Cartesian grids and more general orthogonal unstructured meshes, we can prove that our new scheme can preserve the positivity of k and $\epsilon$. This is achieved via a special implicit discretization of the stiff algebraic relaxation source terms, using a suitable combination of the discrete evolution equations for the logarithms of k and $\epsilon$. The method is applied to some classical academic benchmark problems for non-Newtonian and turbulent flows in two space dimensions, comparing the obtained numerical results with available exact or numerical reference solutions. In all cases, an excellent agreement is observed.

1. Introduction

The study of the incompressible Navier–Stokes equations has started many decades ago, motivated by the necessity of modelling numerous industrial processes and natural phenomena. However, more recently, with the growing development of high-performance supercomputers, numerical methods have become essential for its solution. Even if nowadays we account for a large number of well-established approaches, including different kinds of finite differences (FD), finite volumes (FV), and finite element (FE) methods, the numerical solution of incompressible flows is still a very active field of research.

According to the time discretization approach, we can divide most numerical methods for time-dependent PDEs among explicit and implicit approaches. Explicit methods directly compute the solution at a new time step from the previous one, but the allowed time step is limited by a CFL-type condition that is necessary for the stability of the scheme. On the other hand, fully implicit schemes allow in principle the use of arbitrarily large time steps, only bounded by accuracy and resolution requirements for the problem consideration, but require the solution of complex nonlinear systems. Trying to keep the best of each approach, semi-implicit schemes have been developed, showing an excellent performance for the solution of many PDE systems using either FD [1,2,3,4,5,6], FV [6,7,8,9,10,11,12,13] or continuous and discontinuous Galerkin FE methods [14,15,16,17,18,19].

Focusing on the numerical solution of the incompressible Navier–Stokes equations, a classical approach consists of applying an operator splitting strategy decoupling the computation of momentum and pressure variables [1,20,21,22]. Consequently, we can end up with a hyperbolic dominated subsystem for the velocity unknown and a Poisson problem for the pressure. If we now take into account that explicit FV methods are well suited for solving hyperbolic equations and implicit FE methods are well known to be suitable and efficient for the solution of elliptic Poisson-type problems, we will end up with a hybrid FV/FE methodology that solves each subproblem with the most suitable scheme. The operator splitting between pressure-independent purely convective terms and pressure terms can also be extended to the compressible case, see [23,24,25,26,27,28], keeping a reasonable time step restriction depending only on the velocity field, rather than on the sound speed.

Let us note that, despite turbulence playing an important role in many industrial applications, such as in the automotive and aerospace industry, most of the aforementioned schemes have been proposed in the context of laminar flows. Similarly, the Newtonian behavior of the fluid is a general assumption when addressing incompressible flows. However, there are many natural processes in which non-Newtonian fluids need to be considered, like blood and lava flows, or when studying the behavior of different materials, like toothpaste or honey. A common challenge between these two kinds of flows is the increasing complexity and weight of viscous terms with respect to the convective terms that usually dominate the incompressible Navier–Stokes equations. A careful discretization of the governing PDE system is therefore needed, which is able to deal with both low and high Reynolds number flows. In this paper, we aim at extending the hybrid FV/FE methodology presented in the series of papers [27,28,29,30,31,32], to the case of turbulent flows and non-Newtonian power-law fluids with yield stress.

To deal with turbulent flows, we consider the Reynolds averaged Navier–Stokes (RANS) system coupled with the $k-\epsilon$ turbulence model [33], a classical approach providing good results while keeping a low computational cost. Accordingly, the RANS equations are enlarged with two extra advection-diffusion-reaction equations for the turbulent kinetic energy k and the energy dissipation rate $\epsilon$, coupled with the momentum conservation equation through the viscous stress tensor. Let us note that one of the major difficulties of the new pair of PDEs is the numerical treatment of the algebraic source terms that depend on the flow variables and may become highly stiff. To overcome this issue, we will define an ODE subsystem whose solution will be used to update the advection-diffusion equations and that perfectly adjusts to the operator splitting strategy used in this paper. Following the seminal ideas in [34,35,36,37], instead of working directly with the turbulence variables k and $\epsilon$, we will make use of an ODE subsystem based on its logarithmic counterparts so that the positivity preserving property of k and $\epsilon$ is essentially verified by construction. For a certain class of methods, we are even able to prove strict positivity preservation in each substep of the scheme. This is crucial for the algorithm to provide physically consistent solutions [38,39,40], and is often also called realizable turbulence model.

From the rheology point of view, we focus on the so-called Herschel–Bulkley fluids [41]. The corresponding model is a generalization of the power-law model, which includes the yield stress and a non-linearity of the shear stress, and that have Newtonian and Bingham [42] fluids as special sub-cases. We then need to deal with the discontinuity in the rheological behavior of the flow arising from the presence of the yield stress. In the literature, there are two main ways of overcoming this issue: the first one is a variational reformulation using multipliers and combining the constitutive law with the momentum equation [43,44,45], and the second one is the regularization of the constitutive law [46,47,48,49]. In this paper, we choose the second option, i.e., we regularize the constitutive model taking into account the regularization proposed in [49]. Further advances in the simulation of non-Newtonian fluids in the framework of Stokes and Navier–Stokes equations can be found in [48,50] and of the GPR model for continuum mechanics in [51,52].

As already mentioned, one of the most significant numerical difficulties that we encounter when addressing turbulence and non-Newtonian fluids is the time step restriction that may result in an explicit discretization of the viscous terms. Therefore, differently from what has been done in former hybrid FV/FE methods [53], we will apply operator splitting again to define two subsystems. The first one contains the advection terms and will be solved explicitly. Meanwhile, the viscous subsystem will be treated implicitly. Consequently, the CFL condition of the explicit part of the scheme does not depend on the diffusion terms and allows larger time steps than a completely explicit discretization. In fact, this combination of implicit/explicit schemes to solve advection-diffusion equations is a classical approach used to take advantage of the strengths of both FV and FE methods in mixed schemes. The usual procedure for combining them is to consider a finite volume scheme for the approximation of convective terms while diffusion is handled by a continuous or discontinuous Galerkin finite element approach. Many examples can be found in the literature, including both collocated [40,54,55] and staggered grids [56,57,58]. Furthermore, theoretical analysis for this kind of combined FV-FE methods, applied to general advection-diffusion-reaction equations or the particular case of Navier–Stokes equations, is also available, see for instance [59,60,61,62].

Another important point that we have not raised yet is the spatial discretization of the 2D computational domain. We will consider two grid arrangements: unstructured simplex meshes and Cartesian grids for the finite element method. In both cases, edge-based dual grids will be generated for the finite volume scheme. Apart from the location of the discrete unknowns, one of the main differences between the algorithms proposed for the two kinds of meshes is how diffusion is addressed. If the domain is discretized by means of unstructured grids, in this paper we will use implicit continuous finite elements, while implicit finite volumes are chosen for Cartesian grids.

To attain a second order finite volume scheme on unstructured grids, a local ADER methodology is used. Furthermore, in this context, we take advantage of the dual mesh structure, and a Galerkin approach is employed for the computation of the gradients needed in the reconstruction and half-in-time evolution steps [63,64]. Further details on the classical ADER approach are presented in [65,66,67,68,69], and the modern ADER-DG methodology can be found, for instance, in [70,71]. For a more comprehensive review of the latest advances in this field, we refer to [72,73]. Regarding the Cartesian-based algorithm, the classical MUSCL-Hancock method is employed [74,75].

The rest of the paper is organized as follows. In Section 2, we present the governing equations and introduce the proposed flux splitting. Section 3 is devoted to the numerical method. We first concentrate on the case of unstructured grids recalling the computation of convection terms using explicit FV. Then, we introduce the implicit algorithm for viscous terms and define the weak problem associated with the pressure subsystem. The interpolations used to pass data between the staggered grids are provided, and the identity interpolation property is proven. The second part of the section focuses on the novel hybrid FV/FE algorithm on staggered non-uniform Cartesian grids. Again, all stages are detailed, highlighting the differences with respect to the scheme on unstructured meshes. Furthermore, the positivity preserving discretization of the source terms of the $k-\epsilon$ model using an ODE system based on the time evolution of the logarithms of k and $\epsilon$ is described, and the positivity preserving property on orthogonal grids is rigorously demonstrated for each substep of the scheme. The validation of the proposed methods is carried out in Section 4. Several numerical test cases for both, non-Newtonian fluids and simple turbulent flows are presented, and the obtained numerical results are compared against available analytical and numerical reference solutions. Finally, in Section 5, the main conclusions and an outlook are drafted.

2. Governing Partial Differential Equations

In this paper, we consider the incompressible Reynolds-Averaged Navier–Stokes (RANS) equations with the k − $\epsilon$ turbulence model that consist of the mass and momentum conservation laws and the evolution equations for the turbulence quantities k and $\epsilon$. The mathematical model for incompressible flows, written in terms of conservative variables, reads

$$\begin{array}{ccc}\hfill \nabla ·\mathbf{v}& =& 0,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}+\nabla ·\left(\rho \mathbf{v}\otimes \mathbf{v}\right)+\nabla p-\nabla ·\mathit{\tau }& =& 0,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{\partial \rho k}{\partial t}+\nabla ·\left(\mathbf{v}\rho k\right)-\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right]& =& {\mu }_t{G}_k-\rho \epsilon ,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \frac{\partial \rho \epsilon }{\partial t}+\nabla ·\left(\mathbf{v}\rho \epsilon \right)-\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]& =& C_{1\epsilon }\frac{\epsilon }{k}{\mu }_t{G}_k-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k},\hfill \end{array}$$

(4)

where $\rho =const.$ is the constant fluid density, $\mathbf{v}=(u,v,w)$ is the velocity vector, p is the pressure, k is the turbulent kinetic energy, $\epsilon$ the turbulent dissipation rate and the shear stress tensor $\mathit{\tau }$ is given by

$$\mathit{\tau }=\left(\mu +{\mu }_t\right)\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right)-\frac{2}{3}\rho k\mathbf{I}:=\sigma -\frac{2}{3}\rho k\mathbf{I},\mathrm{with}{\mu }_t=C_{\mu }\rho \frac{k^2}{\epsilon },$$

(5)

where $\mu$ denotes the laminar viscosity and ${\mu }_t$ the turbulent viscosity. The production term $G_k$ in (3) depends on the velocity gradient and reads

$$G_k=\frac{1}{2}{∥\nabla \mathbf{v}+\nabla {\mathbf{v}}^T∥}^2\ge 0.$$

Throughout this paper, we will make use of the following notation for the vector of conservative variables: $\mathbf{w}={\left(\rho {\mathbf{v}}^T,\rho k,\rho \epsilon \right)}^T={\left({\mathbf{w}}_{\mathbf{v}}^T,w_k,w_{\epsilon }\right)}^T$, which allows rewriting the evolutionary part of the PDE system, i.e., (2)–(4), in more compact form as

$$\frac{\partial \mathbf{w}}{\partial t}+\nabla ·{\mathcal{F}}^c+\nabla ·{\mathcal{F}}^p+\nabla ·{\mathcal{F}}^v=\mathcal{S}\left(\mathbf{w}\right).$$

(6)

Here, the convective flux tensor ${\mathcal{F}}^c=\left({\mathbf{f}}^c,{\mathbf{g}}^c,{\mathbf{h}}^c\right)={\left(\mathbf{v}\otimes {\mathbf{w}}_{\mathbf{v}}+\frac{2}{3}{w}_k\mathbf{I},\mathbf{v}{w}_k,\mathbf{v}{w}_{\epsilon }\right)}^T$ contains the hyperbolic part of the PDE system (2)–(4), while the parabolic terms are in the viscous flux tensor ${\mathcal{F}}^v=\left({\mathbf{f}}^v,{\mathbf{g}}^v,{\mathbf{h}}^v\right)={\left(-\sigma ,-\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k,-\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right)}^T$. Furthermore, ${\mathcal{F}}^p={\left(p\mathbf{I},0,0\right)}^T$ is the flux tensor due to the pressure, and $\mathcal{S}\left(\mathbf{w}\right)$ is the algebraic source term corresponding to the right-hand side of (2)–(4). While the definition of the stress tensor (5) is universal in this paper, the rheology for calculating the viscosity depends on the nature of the flow. More precisely, we consider either laminar non-Newtonian flows or turbulent Newtonian flows. In principle, also a combination of both is straightforward, but later we want to compare with existing results from classical test problems that are available in the literature for each flow rheology separately, and therefore this topic is left to future research. Below, we describe both the laminar non-Newtonian rheology, as well as the turbulent Newtonian one. Throughout this paper the international system of units (SI) is employed. The governing equations are presented for the most general case in three space dimensions, and the unstructured hybrid FV/FE method has also been implemented in 2D and 3D, see [27,28,30,31]. However, the numerical test problems shown later are restricted to the two-dimensional case, i.e., assuming planar problems with $\partial /\partial z=0$.

2.1. Laminar Non-Newtonian Fluids

In order to describe a laminar flow, it is sufficient to set the coefficients $C_{\mu }=C_{1\epsilon }=C_{2\epsilon }=0$ in the model (1)–(4). Furthermore, one needs to set $k=\epsilon =0$.

For non-Newtonian fluids, the viscosity coefficient $\mu$ depends on the rate of shear tensor

$$\dot{\gamma }=\dot{\epsilon }-\frac{1}{3}\mathrm{tr}\left(\dot{\epsilon }\right)\mathbf{I},with\dot{\epsilon }=\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right),$$

where $\frac{1}{2}\dot{\epsilon }$ is the strain rate tensor. The rheology for the viscosity $\mu$ considered in this paper is the one of a Herschel–Bulkley fluid [76] with yield stress ${\sigma }_y$ and reads as follows:

$$\mu =\kappa {\dot{\gamma }}^{n-1}+{\sigma }_y{\dot{\gamma }}^{-1},\mathrm{with}\dot{\gamma }=\parallel \dot{\gamma }\parallel =\sqrt{\frac{1}{2}{\dot{\gamma }}_{ij}{\dot{\gamma }}_{ij}}.$$

(7)

Here, $\kappa$ is the so-called consistency index, ${\sigma }_y$ is the yield stress of the fluid, and n is the power-law index which represents the deviation from Newtonian behavior. The Bingham model is a particular case of the Herschel–Bulkley (HB) model (7), where $n=1$ and the power-law model is a particular case where ${\sigma }_y=0$. For $\dot{\gamma }\to 0$ the viscosity $\mu \to \infty$ in the cases where $n<1$ or ${\sigma }_y>0$, and this fact leads into numerical difficulties (see [77]), hence making an implicit discretization of the viscous terms mandatory. To avoid infinite viscosities, in this paper we consider the widely used regularization of Papanastasiou, proposed in [49,78]:

$$\mu =\kappa {\dot{\gamma }}^{n-1}+{\sigma }_y{\dot{\gamma }}^{-1}\left(1-exp(-m\dot{\gamma })\right),$$

with the regularization constant set to $m=1000$. Moreover, we can define the generalized Bingham number $Bi$ for the HB model as

$$Bi=\frac{{\sigma }_y}{\kappa }{\left(\frac{h}{V}\right)}^n,$$

(8)

where V is a characteristic velocity, h is the characteristic length (usually the channel width) and the generalized Reynolds number can be written as (see [79,80])

$$Re=\frac{\rho V^{2-n}{h}^n}{\kappa }.$$

(9)

2.2. Turbulent Newtonian Flows

To model turbulent flow regimes at the aid of the RANS equations with the $k-\epsilon$ model, i.e., with Equations (1)–(4), one chooses in general a constant molecular viscosity $\mu$, while the turbulent viscosity is given by the usual relation

$${\mu }_t=C_{\mu }\frac{k^2}{\epsilon },$$

see also (5). The parameters for the k − $\epsilon$ model are typically chosen as $C_{\mu }=0.09$, $C_{1\epsilon }=1.44$ and $C_{2\epsilon }=1.92$, which are semi-empirical closure constants, and ${\sigma }_k=1.0,$${\sigma }_{\epsilon }=1.3$ that are the turbulent Prandtl numbers.

A critical point when simulating turbulent flows is the near-wall region within turbulent boundary layers. From the numerical point of view, capturing a turbulent boundary layer correctly is very demanding and requires the use of very fine grids. Moreover, the standard $k-\epsilon$ model is no longer valid in the viscous sublayer very close to the wall. Nevertheless, several approaches have been developed in the literature, aiming at overcoming this issue, like wall laws and low Reynolds number turbulence models.

This work will focus on a realizable discretization of the $k-\epsilon$ model that assures the positivity of k and $\epsilon$. Furthermore, we will make use of the low Reynolds number $k-\epsilon$ model of Lam and Bremhorst [81] that can be used to solve the entire turbulent boundary layer, from the outer region down to the wall, capturing both the viscous sublayer and the logarithmic layer.

The Lam-Bremhorst model suggests modifying the coefficients $C_{\mu }$, $C_{1\epsilon }$ and $C_{2\epsilon }$ of the standard $k-\epsilon$ model based on two characteristic Reynolds numbers $R_y$ and $R_t$ that contain the distance to the wall as well as a characteristic size of the turbulent eddies, respectively. They read

$$R_y=\frac{\rho \sqrt{k}y}{\mu },\mathrm{and}{R}_t=\frac{\rho k^2}{\mu \epsilon }.$$

(10)

The standard model coefficients $C_{\mu }$, $C_{1\epsilon }$, and $C_{2\epsilon }$ are then modified in terms of $R_y$ and $R_t$ and three functions $f_{\mu }$, $f_1$ and $f_2$ as follows:

$$C_{\mu }={\overline{C}}_{\mu }{f}_{\mu },C_{1\epsilon }={\overline{C}}_{1\epsilon }{f}_1,C_{2\epsilon }={\overline{C}}_{2\epsilon }{f}_2,$$

(11)

where

$$\begin{array}{ccc}\hfill f_{\mu }& =& {\left(1-exp(-0.0165R_y)\right)}^2\left(1+20.5/R_t\right),\hfill \\ \hfill f_1& =& \left(1+{\left(\frac{0.05}{f_{\mu }}\right)}^3\right),\hfill \\ \hfill f_2& =& max\left(1.01/{\overline{C}}_{2\epsilon },\left(1-exp(-R_t^2)\right)\right)\hfill \end{array}$$

(12)

with ${\overline{C}}_{\mu }=0.09$, ${\overline{C}}_{1\epsilon }=1.44$ and ${\overline{C}}_{2\epsilon }=1.92$ the constants of the original $k-\epsilon$ model. Compared to the original Lam-Bremhorst model in this work the function $f_2$ has been slightly modified in order to guarantee the realizability condition $C_{2\epsilon }>1$, which will be required later.

For the Lam-Bremhorst model, the following boundary conditions are generally imposed at wall boundaries ${\mathsf{\Gamma }}_w$, see [82]:

$$k=0,\nabla \epsilon ·\mathbf{n}=0,\forall \mathbf{x}\in {\mathsf{\Gamma }}_w,$$

i.e., homogeneous Dirichlet boundary conditions for the turbulent kinetic energy k and homogeneous Neumann boundary conditions for the turbulent dissipation rate $\epsilon$. These boundary conditions are usually straightforward to implement in a numerical code.

2.3. Flux Splitting

As usual in the context of semi-implicit schemes [3,8,18,24,83,84,85,86], to decouple the evolution system (6) into convection, diffusion, pressure terms and algebraic source terms, the above system (1)–(4) is now split into four subsystems:

Convective subsystem

The convective (hyperbolic) subsystem,

$$\frac{\partial \mathbf{w}}{\partial t}+\nabla ·{\mathcal{F}}^c=0,$$

(13)

which in this paper will be discretized explicitly, reads

$$\begin{array}{ccc}\hfill \frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}+\nabla ·\left(\rho \mathbf{v}\otimes \mathbf{v}+\frac{2}{3}\rho k\mathbf{I}\right)& =& 0,\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill \frac{\partial \rho k}{\partial t}+\nabla ·\left(\mathbf{v}\rho k\right)& =& 0,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \frac{\partial \rho \epsilon }{\partial t}+\nabla ·\left(\mathbf{v}\rho \epsilon \right)& =& 0.\hfill \end{array}$$

(16)

Viscous subsystem

The viscous (parabolic) subsystem,

$$\frac{\partial \mathbf{w}}{\partial t}+\nabla ·{\mathcal{F}}^v=0,$$

(17)

which will be discretized implicitly, reads

$$\begin{array}{ccc}\hfill \frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}-\nabla ·\left[\left(\mu +{\mu }_t\right)\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right)\right]& =& 0,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \frac{\partial \rho k}{\partial t}-\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right]& =& 0,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill \frac{\partial \rho \epsilon }{\partial t}-\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]& =& 0.\hfill \end{array}$$

(20)

Pressure subsystem

The elliptic pressure subsystem

$$\begin{array}{ccc}\hfill \nabla ·\mathbf{v}& =& 0,\hfill \\ \hfill \frac{\partial \mathbf{w}}{\partial t}+\nabla ·{\mathcal{F}}^p& =& 0,\hfill \end{array}$$

that must be discretized implicitly reads

$$\begin{array}{ccc}\hfill \nabla ·\mathbf{v}& =& 0,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill \frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}+\nabla p& =& 0.\hfill \end{array}$$

(22)

ODE subsystem including the source terms of the turbulence model

Last but not least, the system of ODE

$$\frac{\partial \mathbf{w}}{\partial t}=\mathcal{S}\left(\mathbf{w}\right),$$

containing the potentially stiff algebraic source terms of the turbulence model and which requires an implicit positivity preserving time discretization, reads

$$\begin{array}{ccc}\hfill \frac{\partial \rho k}{\partial t}& =& {\mu }_t{G}_k-\rho \epsilon ,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \frac{\partial \rho \epsilon }{\partial t}& =& C_{1\epsilon }\frac{\epsilon }{k}{\mu }_t{G}_k-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}.\hfill \end{array}$$

(24)

3. The Hybrid Finite Volume/Finite Element Method

To solve the incompressible RANS equations in combination with the $k-\epsilon$ turbulence model, we extend the family of hybrid finite volume/finite element methods described in [27,28,29,30,31,32]. This methodology relies on a specific combination of explicit and implicit FV and FE methods to solve the subsystems obtained from the flux splitting introduced in the previous section. The main stages can be summarized as:

• Transport stage. The convective subsystem (13) is solved using an explicit FV method, and an intermediate approximation of the conservative variables vector is obtained.
• Viscous stage. The viscous subsystem (17) is discretized at the aid of implicit FE or FV methods.
• Interpolation stage. The intermediate states of the conservative variables are interpolated from the dual mesh to the primal one. See Section 3.1.1 and Section 3.2.1 for further details about how the staggered grids are defined.
• Projection stage. The new discrete pressure field is obtained with an implicit FE method by solving the corresponding Poisson problem that is obtained from the saddle-point problem (21) and (22).
• Post-projection stage. The intermediate approximation of the conservative variables is updated considering the solution computed on the projection stage and the contribution from the source terms of the $k-\epsilon$ equations.

In what follows, we describe the semi-implicit hybrid FV/FE scheme on two different kinds of computational grids: staggered unstructured simplex meshes and staggered Cartesian grids. For each of them, a description of the mesh arrangement is given, the notation is introduced, and the numerical discretization of the four subsystems is provided. Moreover, the positivity-preserving property of the discretization of the source terms in the $k-\epsilon$ model is proven. Throughout this section, capital letters refer to the discrete variables, and the super-index ${}^{*}$ denotes intermediate approximations, that is, if the state vector is defined as $\mathbf{w}={({\mathbf{w}}_{\mathbf{v}},w_k,w_{\epsilon })}^T=(\rho u,\rho v,\rho w,\rho k,\rho \epsilon )$ then, ${\mathbf{W}}^n$ is the discrete approximation of $\mathbf{w}(\mathbf{x},t^n)={(\rho u(\mathbf{x},t^n),\rho v(\mathbf{x},t^n),\rho w(\mathbf{x},t^n),\rho k(\mathbf{x},t^n),\rho \epsilon (\mathbf{x},t^n))}^T$, and ${\mathbf{W}}^{*}$ and ${\mathbf{W}}^{**}$ are the intermediate approximations. For the sake of simplicity, the hybrid FV/FE method is in the following only presented in two space dimensions, but it has already been implemented as well for the unstructured three-dimensional case, see [27,28,30,31].

3.1. Unstructured Simplex Meshes

We first focus on extending the hybrid FV/FE schemes on unstructured simplex meshes to solve turbulent flows and to deal with non-Newtonian fluids. To this end, an implicit discretization of the viscous terms is proposed, and a novel numerical treatment of the k and $\epsilon$ equations is described.

3.1.1. Staggered Unstructured Grids

Let $\mathcal{T}$ denote an unstructured triangular partition of the spatial computational domain $\mathsf{\Omega }\in {\mathbb{R}}^2$, that is, $\overline{\mathsf{\Omega }}={\displaystyle \bigcup _{T_k\in \mathcal{T}}}{T}_k$, where $T_k,i=1,\dots ,n_{\mathrm{elem}}$ are the elements of the primal mesh. From this mesh, we built a dual staggered mesh, such that the barycenters of the edges of the triangles of the primal mesh will be taken as the nodes $N_i,i=1,\dots ,n_{\mathrm{vol}}$ of the elements of the dual mesh, $C_i$, which are called cells (see the details in Figure 1). There are two different kinds of dual volumes: the interior cells, built by merging the two triangles formed by connecting the barycenters B and $B^{\prime }$ of two primal elements, that share an interior edge, with this edge (gray elements in the central and right-hand sketches of Figure 1), and the boundary cells, built connecting a boundary edge with the barycenter of the triangle to which the edge belongs (white elements in the central and right-hand sketches of Figure 1). Moreover, the following mesh-related notation needs to be defined:

• ${\mathcal{K}}_i$ is a subset of nodes formed by the barycenters $N_i$ of the dual elements sharing an edge with the cell $C_i$, i.e., ${\mathcal{K}}_i$ is the set of neighbors of $C_i$.
• $\left|C_i\right|$ is the area of $C_i$, ${\displaystyle {\mathsf{\Gamma }}_i=\bigcup _{N_j\in {\mathcal{K}}_i}{\mathsf{\Gamma }}_{ij}}$ its boundary, and $\tilde{\mathbf{n}}{}_i$ its outward unit normal.
• ${\mathsf{\Gamma }}_{ij}$ is the shared edge between the dual elements $C_i$ and $C_j$, $N_{ij}$ its barycenter, and $\tilde{\mathbf{n}}{}_{ij}$ its outward unit normal vector. Moreover, $\mathbf{n}{}_{ij}:={\tilde{\mathbf{n}}}_{ij}\parallel \mathbf{n}{}_{ij}\parallel$, where $\parallel \mathbf{n}{}_{ij}\parallel =|{\mathsf{\Gamma }}_{ij}|$ represents the length of ${\mathsf{\Gamma }}_{ij}$.

3.1.2. Explicit Discretization of the Convective Terms

In the first stage, the convective subsystem (14)–(16) is discretized explicitly with a finite volume method. Integrating

$${\mathbf{W}}^{*}={\mathbf{W}}^n-\Delta t\nabla ·{\mathcal{F}}^c\left({\mathbf{W}}^n\right),$$

(25)

on each control volume $C_i$ and applying the Gauss theorem, we get the intermediate approximation of the conservative variables, including the convective terms,

$${\mathbf{W}}_i^{*}={\mathbf{W}}_i^n-\frac{\Delta t}{|C_i|}{\int }_{{\mathsf{\Gamma }}_i}{\mathcal{F}}^c\left({\mathbf{W}}^n\right)\tilde{\mathbf{n}}{}_i\mathrm{dS}={\mathbf{W}}_i^n-\frac{\Delta t}{|C_i|}\sum _{C_j\in {\mathcal{K}}_i}∥\mathbf{n}{}_{ij}∥\varphi ({\overline{\mathbf{W}}}_i^n,{\overline{\mathbf{W}}}_j^n,\tilde{\mathbf{n}}{}_{ij}),$$

where ${\mathcal{K}}_i$ is the set of neighbors of the control volume $C_i$. To approximate the flux term with second order of accuracy in space and time, the Rusanov scheme [87], combined with a local ADER methodology [30], has been considered. Accordingly, the integral over the cell boundary ${\mathsf{\Gamma }}_i$ is split into the sum of the integrals over the cell edges ${\mathsf{\Gamma }}_{ij}$, and the integral on ${\mathsf{\Gamma }}_{ij}$ is approximated using a numerical flux function $\varphi$, which reads

$$\varphi ({\overline{\mathbf{W}}}_i^n,{\overline{\mathbf{W}}}_j^n,\tilde{\mathbf{n}}{}_{ij})=\frac{1}{2}\left({\mathcal{F}}^c\left({\overline{\mathbf{W}}}_i^n\right)+{\mathcal{F}}^c\left({\overline{\mathbf{W}}}_j^n\right)\right)·\tilde{\mathbf{n}}{}_{ij}-\frac{1}{2}{\tilde{\alpha }}_{ij}^n\left({\overline{\mathbf{W}}}_j^n-{\overline{\mathbf{W}}}_i^n\right),$$

where the Rusanov coefficient ${\tilde{\alpha }}_{ij}^n$ is given by

$${\tilde{\alpha }}_{ij}^n={\tilde{\alpha }}_{ij}^n({\overline{\mathbf{W}}}_i^n,{\overline{\mathbf{W}}}_j^n,{\tilde{\mathbf{n}}}_{ij}):={\alpha }_{ij}^n({\overline{\mathbf{W}}}_i^n,{\overline{\mathbf{W}}}_j^n,\tilde{\mathbf{n}}{}_{ij})+c_{\alpha },$$

with ${\alpha }_{ij}^n$ the maximum signal speed on the edge, $c_{\alpha }\in {\mathbb{R}}_0^{+}$ an optional artificial viscosity coefficient that can be used to improve the stability properties of the scheme, and ${\overline{\mathbf{W}}}_{i,j}^n$ are the half-in-time evolved reconstructed variables on the left and the right of the element interface, respectively, which are obtained following the local ADER approach [30].

3.1.3. Implicit Discretization of the Viscous Terms

The viscous subsystem is discretized implicitly with a finite element method. After semi-discretization in time the viscous subsystem reads

$${\mathbf{W}}^{**}={\mathbf{W}}^{*}-\Delta t\nabla ·{\mathcal{F}}^v\left({\mathbf{W}}^{**}\right),$$

(26)

where ${\mathbf{W}}^{*}$ is the intermediate solution obtained from (25), which involves only convective terms, and ${\mathbf{W}}^{**}$ is the intermediate approximation including also viscous terms. Let us denote $\mathsf{\Gamma }:=\partial \mathsf{\Omega }$ the boundary of the computational domain, $\mathsf{\Omega }$, and $\mathbf{n}$ the outward-pointing normal vector. Multiplying Equation (26) by a test function $z\in V_0$, $V_0:=\left\{z\in H^1\left(\mathsf{\Omega }\right):{\int }_{\mathsf{\Omega }}z \mathrm{dV}=0\right\}$, and integrating in $\mathsf{\Omega }$, the weak formulation associated with the former system (26), after integration by parts, reads:

Weak problem 1.

Find ${\mathbf{W}}^{**}\in V_0$ verifying

$${\int }_{\mathsf{\Omega }}{\mathbf{W}}^{**}·z \mathrm{dV}-\Delta t{\int }_{\mathsf{\Omega }}{\mathcal{F}}^v\left({\mathbf{W}}^{**}\right)·\nabla z \mathrm{dS}={\int }_{\mathsf{\Omega }}{\mathbf{W}}^{*}·z \mathrm{dV}-\Delta t{\int }_{\mathsf{\Gamma }}{\mathcal{F}}^v\left({\mathbf{W}}^{n+1}\right)·\mathbf{n}z \mathrm{dS},$$

for all $z\in V_0$.

The former weak problem is discretized using ${\mathcal{P}}^1$ finite elements. To avoid nonlinearities, we assume a known constant viscosity coefficient per element, computed from the solution obtained at the previous time step. Hence, the resulting system can be efficiently solved using a matrix-free conjugate gradient algorithm. Let us also note that, since the value of the intermediate approximation ${\mathbf{W}}^{*}$ is computed on the dual mesh, to solve (26), we need to interpolate ${\mathbf{W}}^{*}$ from the dual mesh nodes to the primal mesh vertices. First, the local interpolation polynomials at each primal element passing through the barycenters of the edges are computed and evaluated on each vertex. Then, the final value on each vertex is computed as a weighted average of the values obtained at that vertex in each primal element. Further details on the properties of the employed interpolations are provided in Section 3.1.6.

3.1.4. Finite Element Discretization of the Pressure System

Once the intermediate value for the velocity field ${\mathbf{V}}^{**}$, including the convective and viscous terms, has been computed, the pressure and the final velocity fields are obtained using a projection method.

Projection stage

The semi-discrete pressure system results

$$\begin{array}{ccc}\hfill \rho {\mathbf{V}}^{n+1}& =& \rho {\mathbf{V}}^{**}-\Delta t\nabla P^{n+1},\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill \nabla ·{\mathbf{V}}^{n+1}& =& 0.\hfill \end{array}$$

(28)

Introducing (27) in Equation (28) multiplied by the constant density $\rho$, yields

$$-\nabla ·\nabla P^{n+1}=-\frac{1}{\Delta t}\nabla ·\rho {\mathbf{V}}^{**},$$

which can be seen as a Poisson-type problem to be solved using ${\mathcal{P}}^1$ finite elements. Multiplication by a test function $z\in V_0$ and integration over $\mathsf{\Omega }$, after integration by parts leads to the following weak formulation:

Weak problem 2.

Find $P^{n+1}\in V_0$ verifying

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}\nabla P^{n+1}·\nabla z \mathrm{dV}=\frac{1}{\Delta t}{\int }_{\mathsf{\Omega }}\rho {\mathbf{V}}^{**}·\nabla z \mathrm{dV}-\frac{1}{\Delta t}{\int }_{\mathsf{\Gamma }}\rho {\mathbf{V}}^{n+1}·\mathbf{n}z \mathrm{dS},\end{array}$$

(29)

for all $z\in V_0$.

The algebraic system resulting from the discretization of (29) is again solved with a matrix-free conjugate gradient method.

Post-projection stage

Once the new discrete solution of the pressure $P^{n+1}$ has been obtained, the vector of conservative variables can be updated, taking into account that according to (27)

$$\begin{array}{ccc}\hfill \rho {\mathbf{V}}^{n+1}& =& \rho {\mathbf{V}}^{**}-\Delta t\nabla P^{n+1}.\hfill \end{array}$$

Let us note that the degrees of freedom of the intermediate variables ${\mathbf{W}}^{**}$ and the new pressure $P^{n+1}$ are defined in the vertices of the primal mesh, while the conservative variables ${\mathbf{W}}^{n+1}$ are needed in the dual control volumes of the finite volume scheme to be able to run the next time step. Therefore, we first interpolate the contribution of viscous terms to the dual nodes and then add the contribution of the pressure gradient by using a weighted average of the values obtained at the finite elements related to the dual cell.

3.1.5. Algebraic Source Terms of the $k-\epsilon$ Model

The contribution of source terms in the k and $\epsilon$ equations is computed solving the ODE system (23) and (24). The methodology employed is the same for both the unstructured and Cartesian-based meshes algorithms and can be found in Section 3.3, where also the positivity-preserving property is demonstrated.

3.1.6. Interpolation between Staggered Grids

When dealing with staggered grids, special attention should be paid to the interpolation technique used to minimize artificial diffusion effects that might even ruin the overall accuracy of the scheme [55]. In the hybrid FV/FE algorithm for unstructured grids, we have seen the necessity of passing data between the vertices of the primal elements and the nodes of the dual cells. The assumption of the dual nodes to be located not in the barycenter of dual cells but in the barycenter of the edges of the primal grid allow us to define a special interpolation which has the property that the interpolation between the vertices of the primal mesh and the nodes of the dual grid is exactly reversible. In other words, the data interpolation from the primal mesh to the dual grid and back to the primal mesh is the identity operator.

To show this property, we start focusing on a triangle $T_k$ of the primal mesh with vertices of coordinates ${\mathbf{V}}_1,{\mathbf{V}}_2,{\mathbf{V}}_3\in {\mathbb{R}}^2$, and edge midpoints ${\mathbf{N}}_I=\frac{1}{2}\left({\mathbf{V}}_1+{\mathbf{V}}_2\right),$${\mathbf{N}}_{II}=\frac{1}{2}\left({\mathbf{V}}_2+{\mathbf{V}}_3\right),$${\mathbf{N}}_{III}=\frac{1}{2}\left({\mathbf{V}}_3+{\mathbf{V}}_1\right)$, see Figure 2. Denoting by ${\varphi }_1=1-\xi -\eta$, ${\varphi }_2=\xi$ and ${\varphi }_3=\eta$ the basis functions on the reference triangle of the ${\mathcal{P}}^1$ finite elements and being q an arbitrary variable with $q_i=q\left({\mathbf{V}}_i\right)$ the values at the vertices of the primal element, then there exists a unique interpolation polynomial of degree one, $P_1\left(\mathbf{x}\right)$, satisfying the interpolation property

$$P_1\left({\mathbf{V}}_i\right)=q_i,i\in \left\{1,2,3\right\},$$

which reads

$$P_1\left(\mathbf{x}\right)=P_1\left(\xi ,\eta \right)=q_1\left(1-\xi -\eta \right)+q_2\xi +q_3\eta ,\mathbf{x}={\mathbf{V}}_1+\xi \left({\mathbf{V}}_2-{\mathbf{V}}_1\right)+\eta \left({\mathbf{V}}_3-{\mathbf{V}}_1\right).$$

We, therefore, get the following interpolation from the vertices into the edge midpoints:

$$\begin{array}{c}{P}_1\left({\mathbf{N}}_1\right)=P\left(\frac{1}{2},0\right)=q_1\left(1-\frac{1}{2}\right)+q_2\frac{1}{2}=\frac{1}{2}\left(q_1+q_2\right):=q_I,\\ P_1\left({\mathbf{N}}_2\right)=P\left(\frac{1}{2},\frac{1}{2}\right)=q_1\left(1-\frac{1}{2}-\frac{1}{2}\right)+q_2\frac{1}{2}+q_3\frac{1}{2}=\frac{1}{2}\left(q_2+q_3\right):=q_{II},\\ P_1\left({\mathbf{N}}_3\right)=P\left(0,\frac{1}{2}\right)=q_1\left(1-0-\frac{1}{2}\right)+q_3\frac{1}{2}=\frac{1}{2}\left(q_1+q_3\right):=q_{III},\end{array}$$

corresponding to the proposed interpolation from FE vertices to FV nodes.

We now consider the backward interpolation from dual nodes to primal vertices. Given a primal element $T_k$ the basis functions associated with the barycenters of its edges read ${\psi }_I=1-2\eta$, ${\psi }_{II}=-1+2\xi +2\eta$, ${\psi }_{III}=1-2\xi$, which are the basis functions of linear Crouzeix-Raviart elements. Again, there exists a unique first order interpolation polynomial,

$$P_{k,1}\left(\mathbf{x}\right)=P_{k,1}\left(\xi ,\eta \right)=q_I\left(1-2\eta \right)+q_{II}\left(-1+2\xi +2\eta \right)+q_{III}\left(1-2\xi \right),$$

passing through ${\mathbf{N}}_I$, ${\mathbf{N}}_{II}$ and ${\mathbf{N}}_{III}$. Evaluating it in the vertices of $T_k$ we get

$$\begin{array}{c}{P}_{k,1}\left({\mathbf{V}}_1\right)=P\left(0,0\right)=q_I-q_{II}+q_{III}=q_1,\\ P_{k,1}\left({\mathbf{V}}_2\right)=P\left(1,0\right)=q_I+q_{II}-q_{III}=q_2,\\ P_{k,1}\left({\mathbf{V}}_3\right)=P\left(0,1\right)=-q_I+q_{II}+q_{III}=q_3.\end{array}$$

Hence, the original values in the vertices are identically reproduced for any primal element $T_k$ selected.

Finally, denoting $S_i$ the set of elements containing a vertex ${\mathbf{V}}_i$, and choosing any set of weights ${\left\{{\omega }_j\right\}}_{j\in S_i}$ such that ${\sum }_{{\mathsf{\Omega }}_j\in S_i}{\omega }_j=1$, the value of q at each vertex of the primal mesh can be obtained as the weighted average

$$q_i=\sum _{{\mathsf{\Omega }}_j\in S_i}{\omega }_j{P}_{j,1}\left({\mathbf{V}}_i\right).$$

3.2. Cartesian Grids

To the best knowledge of the authors, this is the first time that a semi-implicit hybrid finite volume/finite element scheme is proposed on staggered Cartesian meshes, since previous work on hybrid FV/FE schemes was focused on unstructured simplex meshes in two and three space dimensions, see [27,28,30,31,32]. Here, we propose to use a staggered Cartesian mesh similar to the one used in [8,9,12,25,88]. Overall, we have three overlapping edge-based staggered meshes.

3.2.1. Staggered Cartesian Grid Configuration and Notation

In this section, the barycenters of the non-uniform Cartesian primal mesh are denoted by ${\mathbf{x}}_{i,j}=(x_i,y_j)$, and the corresponding cell size is denoted by $\Delta x_i=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and $\Delta y_j=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$; the barycenters of the staggered edge-based dual cells in $x-$direction for the velocity component u are denoted by ${\mathbf{x}}_{i+\frac{1}{2},j}=(x_{i+\frac{1}{2}},y_j)$ with respective cell size $\Delta x_{i+\frac{1}{2}}=\frac{1}{2}(\Delta x_i+\Delta x_{i+1})=x_{i+1}-x_i$ and $\Delta y_j$; the barycenters of the staggered edge-based dual cells in $y-$direction for the velocity component u are denoted by ${\mathbf{x}}_{i,j+\frac{1}{2}}=(x_i,y_{j+\frac{1}{2}})$ with respective cell size $\Delta x_i$ and $\Delta y_{j+\frac{1}{2}}=\frac{1}{2}(\Delta y_j+\Delta y_{j+1})=y_{j+1}-y_j$, see Figure 3 for clarity. The elements of the primal mesh are ${\mathsf{\Omega }}_{i,j}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times [y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$; those of the edge-based staggered dual mesh in $x-$direction are denoted by ${\mathsf{\Omega }}_{i+\frac{1}{2},j}=[x_i,x_{i+1}]\times [y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$, and those of the edge-based staggered mesh in $y-$direction are defined as ${\mathsf{\Omega }}_{i,j+\frac{1}{2}}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times [y_j,y_{j+1}]$.

The pressure field is discretized via a conforming ${\mathcal{Q}}^1$ finite element method with the pressure defined in the vertices of the primal control volumes ${\mathsf{\Omega }}_{i,j}$, i.e., the pressure degrees of freedom at time $t^n$ are denoted by $P_{i+\frac{1}{2},j+\frac{1}{2}}^n$. The velocity components are defined in the edges of the primal control volumes as $u_{i+\frac{1}{2},j}^n$ and $v_{i,j+\frac{1}{2}}^n$, respectively. The turbulence quantities k and $\epsilon$ are defined in the barycenters of the primal mesh, i.e., $k_{i,j}^n$ and ${\epsilon }_{i,j}^n$.

The interpolation of the velocity field from one mesh to another is achieved via the following simple interpolation rules whenever needed:

$$u_{i,j}^n=\frac{1}{2}\left(u_{i+\frac{1}{2},j}^n+u_{i-\frac{1}{2},j}^n\right),v_{i,j}^n=\frac{1}{2}\left(v_{i,j+\frac{1}{2}}^n+v_{i,j-\frac{1}{2}}^n\right),$$

and

$$u_{i+\frac{1}{2},j}^n=\frac{\Delta x_i{u}_{i,j}^n+\Delta x_{i+1}{u}_{i+1,j}^n}{2\Delta x_{i+\frac{1}{2}}},v_{i,j+\frac{1}{2}}^n=\frac{\Delta y_j{v}_{i,j}^n+\Delta y_{j+1}{v}_{i,j+1}^n}{2\Delta y_{j+\frac{1}{2}}}.$$

3.2.2. Explicit Discretization of the Convective Terms

If ${\mathbf{W}}_{i,j}^n=(\rho u_{i,j}^n,\rho v_{i,j}^n,\rho k_{i,j}^n,\rho {\epsilon }_{i,j}^n)$ is the discrete state vector, the convective subsystem (14)–(16) is discretized using an explicit Godunov-type finite volume scheme on the main grid as follows:

$${\mathbf{W}}_{i,j}^{*}={\mathbf{W}}_{i,j}^n-\frac{\Delta t}{\Delta x_i}\left({\mathbf{f}}_{i+\frac{1}{2},j}^c-{\mathbf{f}}_{i-\frac{1}{2},j}^c\right)-\frac{\Delta t}{\Delta y_j}\left({\mathbf{g}}_{i,j+\frac{1}{2}}^c-{\mathbf{g}}_{i,j-\frac{1}{2}}^c\right),$$

with the numerical flux functions for a first order upwind scheme in x- and y-direction defined as:

$${\mathbf{f}}_{i+\frac{1}{2},j}^c=\frac{1}{2}{u}_{i+\frac{1}{2},j}^n\left({\mathbf{W}}_{i+1,j}^n+{\mathbf{W}}_{i,j}^n\right)-\frac{1}{2}\left|u_{i+\frac{1}{2},j}^n\right|\left({\mathbf{W}}_{i+1,j}^n-{\mathbf{W}}_{i,j}^n\right)$$

and

$${\mathbf{g}}_{i,j+\frac{1}{2}}^c=\frac{1}{2}{v}_{i,j+\frac{1}{2}}^n\left({\mathbf{W}}_{i,j+1}^n+{\mathbf{W}}_{i,j}^n\right)-\frac{1}{2}\left|v_{i,j+\frac{1}{2}}^n\right|\left({\mathbf{W}}_{i,j+1}^n-{\mathbf{W}}_{i,j}^n\right).$$

We emphasize that here the velocities $u_{i+\frac{1}{2},j}^n$ and $v_{i,j+\frac{1}{2}}^n$ are the normal velocity components on the edges resulting from the projection stage, i.e., from the solution of the elliptic pressure Poisson equation at the previous time step. In order to reach second order of accuracy, the boundary-extrapolated and time-evolved values of a second order TVD MUSCL-Hancock scheme can be used instead of the cell averages. The second order fluxes read

$${\mathbf{f}}_{i+\frac{1}{2},j}^c=\frac{1}{2}{u}_{i+\frac{1}{2},j}^n\left({\mathbf{W}}_{i+\frac{1}{2},j}^{+}+{\mathbf{W}}_{i+\frac{1}{2},j}^{-}\right)-\frac{1}{2}\left|u_{i+\frac{1}{2},j}^n\right|\left({\mathbf{W}}_{i+\frac{1}{2},j}^{+}-{\mathbf{W}}_{i+\frac{1}{2},j}^{-}\right)$$

and

$${\mathbf{g}}_{i,j+\frac{1}{2}}^c=\frac{1}{2}{v}_{i,j+\frac{1}{2}}^n\left({\mathbf{W}}_{i,j+\frac{1}{2}}^{+}+{\mathbf{W}}_{i,j+\frac{1}{2}}^{-}\right)-\frac{1}{2}\left|v_{i,j+\frac{1}{2}}^n\right|\left({\mathbf{W}}_{i,j+\frac{1}{2}}^{+}-{\mathbf{W}}_{i,j+\frac{1}{2}}^{-}\right).$$

with the boundary extrapolated values calculated as

$${\mathbf{W}}_{i\pm \frac{1}{2},j}^{\mp }={\mathbf{W}}_{i,j}^n\pm \frac{1}{2}\Delta x_i{\partial }_x{\mathbf{W}}_{i,j}^n+\frac{1}{2}\Delta t{\partial }_t{\mathbf{W}}_{i,j}^n,$$

$${\mathbf{W}}_{i,j\pm \frac{1}{2}}^{\mp }={\mathbf{W}}_{i,j}^n\pm \frac{1}{2}\Delta y_j{\partial }_y{\mathbf{W}}_{i,j}^n+\frac{1}{2}\Delta t{\partial }_t{\mathbf{W}}_{i,j}^n.$$

The limited slopes in space are given by

$${\partial }_x{\mathbf{W}}_{i,j}^n=\mathrm{minmod}\left(\frac{{\mathbf{W}}_{i+1,j}^n-{\mathbf{W}}_{i,j}^n}{\Delta x_{i+\frac{1}{2}}},\frac{{\mathbf{W}}_{i,j}^n-{\mathbf{W}}_{i-1,j}^n}{\Delta x_{i-\frac{1}{2}}}\right),$$

$${\partial }_y{\mathbf{W}}_{i,j}^n=\mathrm{minmod}\left(\frac{{\mathbf{W}}_{i,j+1}^n-{\mathbf{W}}_{i,j}^n}{\Delta y_{j+\frac{1}{2}}},\frac{{\mathbf{W}}_{i,j}^n-{\mathbf{W}}_{i,j-1}^n}{\Delta y_{j-\frac{1}{2}}}\right),$$

with the classical minmod slope limiter, see [75], while the approximation of the temporal derivative is

$$\begin{array}{ccc}\hfill {\partial }_t{\mathbf{W}}_{i,j}^n& =& -\frac{u_{i+\frac{1}{2},j}^n\left({\mathbf{W}}_{i,j}^n+\frac{1}{2}\Delta x_i{\partial }_x{\mathbf{W}}_{i,j}^n\right)-u_{i-\frac{1}{2},j}^n\left({\mathbf{W}}_{i,j}^n-\frac{1}{2}\Delta x_i{\partial }_x{\mathbf{W}}_{i,j}^n\right)}{\Delta x_i}\hfill \\ & & -\frac{v_{i,j+\frac{1}{2}}^n\left({\mathbf{W}}_{i,j}^n+\frac{1}{2}\Delta y_j{\partial }_y{\mathbf{W}}_{i,j}^n\right)-v_{i,j-\frac{1}{2}}^n\left({\mathbf{W}}_{i,j}^n-\frac{1}{2}\Delta y_j{\partial }_y{\mathbf{W}}_{i,j}^n\right)}{\Delta x_i}.\hfill \end{array}$$

For more details on the MUSCL-Hancock method, see, e.g., the well-known textbook of Toro [75].

3.2.3. Implicit Discretization of the Viscous Terms

Considering the viscous subsystem (18)–(20) and applying a finite volume scheme, we get

$$\Delta x_i\Delta y_j{\mathbf{W}}_{i,j}^{**}=\Delta x_i\Delta y_j{\mathbf{W}}_{i,j}^{*}-\Delta t\Delta y_j\left({\mathbf{f}}_{i+\frac{1}{2},j}^v-{\mathbf{f}}_{i-\frac{1}{2},j}^v\right)-\Delta t\Delta x_i\left({\mathbf{g}}_{i,j+\frac{1}{2}}^v-{\mathbf{g}}_{i,j-\frac{1}{2}}^v\right).$$

(30)

The appropriate discretization of the viscous fluxes is non-trivial since the correct discretization of the stress tensor of the Navier–Stokes equations with variable viscosity coefficient requires cross derivatives due to the term $\nabla {\mathbf{v}}^T$, the discretization of which needs corner gradients and thus leads to a nine-point stencil, while a provable positivity preserving discretization of the parabolic terms in the k and $\epsilon$ equations is best achieved at the aid of a classical two-point flux, hence leading to a five-point stencil. In other words, we split the viscous flux vectors as follows:

$${\mathbf{f}}_{i+\frac{1}{2},j}^v=\left(\begin{array}{c}{\mathbf{f}}_{i+\frac{1}{2},j}^{v,\rho \mathbf{v}}\\ {\mathbf{f}}_{i+\frac{1}{2},j}^{v,\rho k}\\ {\mathbf{f}}_{i+\frac{1}{2},j}^{v,\rho \epsilon }\end{array}\right),\mathrm{and}{\mathbf{g}}_{i,j+\frac{1}{2}}^v=\left(\begin{array}{c}{\mathbf{g}}_{i,j+\frac{1}{2}}^{v,\rho \mathbf{v}}\\ {\mathbf{g}}_{i,j+\frac{1}{2}}^{v,\rho k}\\ {\mathbf{g}}_{i,j+\frac{1}{2}}^{v,\rho \epsilon }\end{array}\right).$$

For the momentum equation, the viscous fluxes are defined via the trapezoidal rule and discretizations of the viscous stress tensor in the corners as

$${\mathbf{f}}_{i+\frac{1}{2},j}^{v,\rho \mathbf{v}}=\frac{1}{2}\left({\sigma }_{i+\frac{1}{2},j+\frac{1}{2}}+{\sigma }_{i+\frac{1}{2},j-\frac{1}{2}}\right)·{\mathbf{e}}_x,{\mathbf{g}}_{i,j+\frac{1}{2}}^{v,\rho \mathbf{v}}=\frac{1}{2}\left({\sigma }_{i+\frac{1}{2},j+\frac{1}{2}}+{\sigma }_{i-\frac{1}{2},j+\frac{1}{2}}\right)·{\mathbf{e}}_y,$$

with the unit normal vectors ${\mathbf{e}}_x={(1,0)}^T$ and ${\mathbf{e}}_y={(0,1)}^T$ in the $x-$ and $y-$direction, respectively, and the stress tensor is calculated in a semi-implicit manner

$${\sigma }_{i+\frac{1}{2},j+\frac{1}{2}}=\left({\mu }_{i+\frac{1}{2},j+\frac{1}{2}}^n+{\mu }_t\left({\mathbf{W}}_{i+\frac{1}{2},j+\frac{1}{2}}^n\right)\right)\left(\nabla {\mathbf{V}}_{i+\frac{1}{2},j+\frac{1}{2}}+\nabla {\mathbf{V}}_{i+\frac{1}{2},j+\frac{1}{2}}^T\right),$$

based on an explicit discretization of the nonlinear viscosity coefficient and the implicit corner gradients of the velocity field

$$\nabla {\mathbf{V}}_{i+\frac{1}{2},j+\frac{1}{2}}=\frac{1}{2}\left(\begin{array}{c}\frac{\left(u_{i+1,j+1}^{**}+u_{i+1,j}^{**}\right)-\left(u_{i,j+1}^{**}+u_{i,j}^{**}\right)}{\Delta x_{i+\frac{1}{2}}}\\ \frac{\left(v_{i+1,j+1}^{**}+v_{i,j+1}^{**}\right)-\left(v_{i+1,j}^{**}+v_{i,j}^{**}\right)}{\Delta y_{j+\frac{1}{2}}}\end{array}\right).$$

For the discretization of the parabolic terms in the k and $\epsilon$ equations, we do not need any cross derivatives, hence a classical two-point flux based on the mid-point rule is enough:

$${\mathbf{f}}_{i+\frac{1}{2},j}^{v,\rho q}=\frac{q_{i+1,j}^{**}-q_{i,j}^{**}}{\Delta x_{i+\frac{1}{2}}},{\mathbf{g}}_{i,j+\frac{1}{2}}^{v,\rho q}=\frac{q_{i,j+1}^{**}-q_{i,j}^{**}}{\Delta y_{j+\frac{1}{2}}},\mathrm{with}q\in \left\{k,\epsilon \right\}.$$

This method has the additional advantage that it is provably positivity preserving, see Section 3.4 later, which is a very important feature for the discretization of the k − $\epsilon$ model. In practice, the resulting linear algebraic system for the unknown velocity field given by (30) can be conveniently solved at the aid of a matrix-free conjugate gradient method, since the system is symmetric and positive definite.

3.2.4. Finite Element Discretization of the Pressure System

Once the preliminary velocity field ${\mathbf{V}}^{**}$ has been computed, which includes the nonlinear convective and the viscous terms, the final velocity field can be obtained via a projection method, as already detailed previously for the case of an unstructured mesh.

Projection stage

The pressure Poisson equation resulting from (21) and (22) reads

$$-\nabla ·\nabla P^{n+1}=-\frac{1}{\Delta t}\nabla ·\rho {\mathbf{V}}^{**},$$

and can be conveniently solved at the aid of a classical conforming ${\mathcal{Q}}^1$ finite element method. The associated weak problem of the projection stage then results:

Weak problem 3.

Find the pressure $P^{n+1}\in V_0$ that satisfies

$${\int }_{\mathsf{\Omega }}\nabla P^{n+1}·\nabla z \mathrm{dV}=\frac{1}{\Delta t}{\int }_{\mathsf{\Omega }}\rho {\mathbf{V}}^{**}·\nabla z \mathrm{dV}-\frac{1}{\Delta t}{\int }_{\mathsf{\Gamma }}\rho {\mathbf{V}}^{n+1}·\mathbf{n}z \mathrm{dS},$$

(31)

for all $z\in V_0$.

Post-projection stage

Once the new pressure field $P_{i,j}^{n+1}$ is known from (31), the velocity components can be updated as

$$\begin{array}{ccc}\hfill \rho u_{i+\frac{1}{2},j}^{n+1}& =& \rho u_{i+\frac{1}{2},j}^{**}-\frac{\Delta t}{2\Delta x_{i+\frac{1}{2}}}\left(\Delta x_i{\partial }_x{P}_{i,j}^{n+1}+\Delta x_{i+1}{\partial }_x{P}_{i+1,j}^{n+1}\right),\hfill \\ \hfill \rho v_{i,j+\frac{1}{2}}^{n+1}& =& \rho v_{i,j+\frac{1}{2}}^{**}-\frac{\Delta t}{2\Delta y_{j+\frac{1}{2}}}\left(\Delta y_j{\partial }_y{P}_{i,j}^{n+1}+\Delta y_{j+1}{\partial }_y{P}_{i,j+1}^{n+1}\right),\hfill \end{array}$$

with the average pressure gradients in the primal cell ${\mathsf{\Omega }}_{i,j}$ resulting from the ${\mathcal{Q}}^1$ discretization as

$${\partial }_x{P}_{i,j}^{n+1}=\frac{1}{2}\left(\begin{array}{c}\frac{\left(P_{i+\frac{1}{2},j+\frac{1}{2}}^{n+1}+P_{i+\frac{1}{2},j-\frac{1}{2}}^{n+1}\right)-\left(P_{i-\frac{1}{2},j+\frac{1}{2}}^{n+1}+P_{i-\frac{1}{2},j-\frac{1}{2}}^{n+1}\right)}{\Delta x_i}\\ \frac{\left(P_{i+\frac{1}{2},j+\frac{1}{2}}^{n+1}+P_{i-\frac{1}{2},j+\frac{1}{2}}^{n+1}\right)-\left(P_{i+\frac{1}{2},j-\frac{1}{2}}^{n+1}+P_{i-\frac{1}{2},j-\frac{1}{2}}^{n+1}\right)}{\Delta y_j}\end{array}\right).$$

3.3. Positivity-Preserving Discretization of the Source Terms of the $k-\epsilon$ Turbulence Model

We first concentrate on the positivity-preserving discretization of the—potentially stiff—algebraic source terms in the $k-\epsilon$ model, which is identical for both, the unstructured and the Cartesian hybrid FV/FE method presented in the previous sections. For $\rho =1$ the ODE subsystem associated with the $k-\epsilon$ model reads

$$\begin{array}{ccc}\hfill \frac{\partial k}{\partial t}& =& C_{\mu }\frac{k^2}{\epsilon }{G}_k-\epsilon ,\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill \frac{\partial \epsilon }{\partial t}& =& C_{1\epsilon }{C}_{\mu }k{G}_k-C_{2\epsilon }\frac{{\epsilon }^2}{k},\hfill \end{array}$$

(33)

with the non-negative production term

$$G_k=\frac{1}{2}{∥\nabla \mathbf{v}+\nabla {\mathbf{v}}^T∥}^2\ge 0.$$

Since throughout this paper we are making use of operator splitting, the initial condition for the above ODE system will be denoted by $k^{**}>0$ and ${\epsilon }^{**}>0$, which have been obtained by a successive explicit discretization of the nonlinear convective terms and an implicit discretization of the dissipative (parabolic) terms, as shown in the previous sections, and which we assume to be positive (for a proof of this statement, see the next section). From (32) and (33), one can easily obtain the evolution equations for the logarithms of k and $\epsilon$, since ${\partial }_{t}lnk=\left({\partial }_{t}k\right)/k$ and ${\partial }_{t}ln\epsilon =\left({\partial }_t\epsilon \right)/\epsilon$. In the following, we will use the abbreviations $\alpha :=lnk$, $\beta :=ln\epsilon$ and $\delta :=\alpha -\beta$. The evolution Equations (32) and (33) for the logarithms of k and $\epsilon$ become

$$\frac{\partial lnk}{\partial t}=C_{\mu }\frac{k}{\epsilon }{G}_k-\frac{\epsilon }{k},\frac{\partial ln\epsilon }{\partial t}=C_{1\epsilon }{C}_{\mu }\frac{k}{\epsilon }{G}_k-C_{2\epsilon }\frac{\epsilon }{k},$$

which now depend only on the ratio $k/\epsilon$ and thus in terms of $\alpha$ and $\beta$ read

$$\begin{array}{ccc}\hfill \frac{\partial \alpha }{\partial t}& =& C_{\mu }{G}_k{e}^{\alpha -\beta }-e^{\beta -\alpha },\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \frac{\partial \beta }{\partial t}& =& C_{1\epsilon }{C}_{\mu }{G}_k{e}^{\alpha -\beta }-C_{2\epsilon }{e}^{\beta -\alpha }.\hfill \end{array}$$

(35)

Subtracting (35) from (34) yields the following scalar ODE for the only unknown $\delta$:

$$\frac{\partial \delta }{\partial t}=C_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{\delta }-\left(1-C_{2\epsilon }\right)e^{-\delta }.$$

(36)

Equation (36) can be integrated exactly using separation of variables. For the initial condition $\delta \left(0\right)={\delta }_0$, one obtains the solution

$$\delta \left(t\right)=ln\left[-\frac{\sqrt{ab}}{a}tanh\left(t\sqrt{ab}-{tanh}^{-1}\left(\frac{aexp\left({\delta }_0\right)}{\sqrt{ab}}\right)\right)\right],$$

with $a=C_{\mu }\left(1-C_{1\epsilon }\right)$ and $b=\left(1-C_{2\epsilon }\right)$. However, in the following, we do not make further use of this solution but proceed with a numerical discretization of (34)–(36), instead. Discretization of (36) via the implicit Euler scheme yields

$${\delta }^{n+1}={\delta }^{**}+\Delta t{C}_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{{\delta }^{n+1}}-\Delta t\left(1-C_{2\epsilon }\right)e^{-{\delta }^{n+1}},$$

or, equivalently, leads to the following nonlinear scalar algebraic equation:

$$g\left({\delta }^{n+1}\right)={\delta }^{n+1}-{\delta }^{**}-\Delta t{C}_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{{\delta }^{n+1}}+\Delta t\left(1-C_{2\epsilon }\right)e^{-{\delta }^{n+1}}=0.$$

(37)

In the following, we prove that for $C_{1\epsilon }>1$ and $C_{2\epsilon }>1$, which is the case for the standard $k-\epsilon$ model, Equation (37) has exactly one real root, i.e., we prove the existence and uniqueness of the discrete solution ${\delta }^{n+1}$ that satisfies (37). Obviously, the function $g\left({\delta }^{n+1}\right)$ is continuous, hence a change in the sign of g would guarantee the existence of at least one real root via the intermediate value theorem (Bolzano’s theorem). In the following, we show that g actually changes sign in $\mathbb{R}$. For sufficiently large ${\delta }^{n+1}$, i.e., in the limit ${\delta }^{n+1}\to +\infty$ we have ${\delta }^{n+1}-{\delta }^{**}-\Delta t{C}_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{{\delta }^{n+1}}\to +\infty$, since $G_k\ge 0$ and $1-C_{1\epsilon }<0$ due to the assumption $C_{1\epsilon }>1$, while the term $\Delta t\left(1-C_{2\epsilon }\right)e^{-{\delta }^{n+1}}\to 0$ for ${\delta }^{n+1}\to +\infty$. On the contrary, for ${\delta }^{n+1}\to -\infty$ we have $-\Delta t{C}_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{{\delta }^{n+1}}\to 0$ and ${\delta }^{n+1}-{\delta }^{**}+\Delta t\left(1-C_{2\epsilon }\right)e^{-{\delta }^{n+1}}\to -\infty$, since $C_{2\epsilon }>1$. Hence, $g\left({\delta }^{n+1}\right)\to +\infty$ for ${\delta }^{n+1}\to +\infty$ and $g\left({\delta }^{n+1}\right)\to -\infty$ for ${\delta }^{n+1}\to -\infty$ and therefore the function $g\left({\delta }^{n+1}\right)$ changes sign in $\mathbb{R}$ and thus must have at least one real root. In order to prove the uniqueness of the solution, we make use of the monotonicity of the function $g\left({\delta }^{n+1}\right)$. Indeed, the first derivative of g reads

$$g^{\prime }\left({\delta }^{n+1}\right)=1-\Delta t{C}_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{{\delta }^{n+1}}-\Delta t\left(1-C_{2\epsilon }\right)e^{-{\delta }^{n+1}}>0\forall {\delta }^{n+1}\in \mathbb{R},$$

hence g is a monotonically increasing function, and thus there exists exactly one real root that satisfies $g\left({\delta }^{n+1}\right)=0$.

In practice, the root can be conveniently found via the bisection algorithm using a sufficiently large starting interval. For all numerical experiments carried out in this paper, we found that the starting interval ${\delta }^{n+1}\in [-100,+100]$ was sufficient. We emphasize that $\alpha$ and $\beta$ represent the logarithms of k and $\epsilon$, and that $\delta$ is their difference, hence the previously mentioned starting interval for the bisection method can actually be considered as very generous for practical applications. Once the root ${\delta }^{n+1}$ of (37) has been found, the corresponding values of ${\alpha }^{n+1}$ and ${\beta }^{n+1}$ can be easily obtained from an implicit Euler discretization of (34) and (35) as follows:

$$\begin{array}{ccc}\hfill {\alpha }^{n+1}& =& {\alpha }^{**}+\Delta t{C}_{\mu }{G}_k{e}^{{\delta }^{n+1}}-\Delta t{e}^{-{\delta }^{n+1}},\hfill \\ \hfill {\beta }^{n+1}& =& {\beta }^{**}+\Delta t{C}_{1\epsilon }{C}_{\mu }{G}_k{e}^{{\delta }^{n+1}}-\Delta t{C}_{2\epsilon }{e}^{-{\delta }^{n+1}}.\hfill \end{array}$$

It is, of course, trivial to see that $k^{n+1}=exp\left({\alpha }^{n+1}\right)>0$ and ${\epsilon }^{n+1}=exp\left({\beta }^{n+1}\right)>0$ are always positive.

3.4. Positivity-Preserving Discretization of the $k-\epsilon$ Model on Orthogonal Unstructured Meshes

In the following, we prove that the proposed discretization of the k − $\epsilon$ model on Cartesian grids is positivity preserving. To make the proof as general as possible, in this section, we will assume a general orthogonal unstructured mesh, instead of a simple Cartesian mesh that covers the computational domain $\mathsf{\Omega }$. The orthogonal mesh is composed of non-overlapping polygons ${\mathsf{\Omega }}_i\in \mathsf{\Omega }$ with centroids ${\mathbf{x}}_i$ and the orthogonality condition $({\mathbf{x}}_j-{\mathbf{x}}_i)\perp \partial {\mathsf{\Omega }}_{ij}$ with $\partial {\mathsf{\Omega }}_{ij}={\mathsf{\Omega }}_i\cap {\mathsf{\Omega }}_j$ the shared edge between elements ${\mathsf{\Omega }}_i$ and ${\mathsf{\Omega }}_j$. It is obvious that a non-uniform Cartesian mesh is just a particular case of an orthogonal unstructured mesh.

Theorem 1.

Positivity-preserving property. Assume $\rho =1=const.$ and a divergence-free velocity field $\mathbf{v}$ that satisfies the discrete divergence-free property

$$\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}{u}_{ij}^n|\partial {\mathsf{\Omega }}_{ij}|=0,\mathrm{with}{u}_{ij}^n={\mathbf{v}}_{ij}^n·{\mathbf{n}}_{ij},$$

(38)

where $u_{ij}^n$ is the normal velocity component on the element boundary, where $\partial {\mathsf{\Omega }}_{ij}={\mathsf{\Omega }}_i\cap {\mathsf{\Omega }}_j$ is the common edge between elements ${\mathsf{\Omega }}_i$ and ${\mathsf{\Omega }}_j$, $|\partial {\mathsf{\Omega }}_{ij}|$ is its length, ${\mathbf{n}}_{ij}$ is the outward-pointing unit normal vector pointing from element ${\mathsf{\Omega }}_i$ to its neighbor element ${\mathsf{\Omega }}_j$ and ${\mathcal{N}}_i$ is the set of neighbors of ${\mathsf{\Omega }}_i$. Then, the semi-implicit finite volume scheme for quantity $q\in \left\{k,\epsilon \right\}$ with ${\mu }_{ij}^n=\frac{1}{2}\left({\mu }_i^n+{\mu }_j^n\right)>0$ and ${\mu }_{t,ij}^n=\frac{1}{2}\left({\mu }_{t,i}^n+{\mu }_{t,j}^n\right)>0$ given by

$$\begin{array}{ccc}\hfill q_i^{*}& =& q_i^n-\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|\left(\frac{1}{2}{u}_{ij}^n\left(q_j^n+q_i^n\right)-\frac{1}{2}\left|u_{ij}^n\right|\left(q_j^n-q_i^n\right)\right),\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill |{\mathsf{\Omega }}_i|q_i^{**}& =& |{\mathsf{\Omega }}_i|q_i^{*}+\Delta t\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|\left({\mu }_{ij}^n+{\mu }_{t,ij}^n\right)\frac{q_j^{**}-q_i^{**}}{|{\mathbf{x}}_j-{\mathbf{x}}_i|},\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill g\left({\delta }_i^{n+1}\right)& =& {\delta }_i^{n+1}-{\delta }_i^{**}-\Delta t{C}_{\mu }\left(1-C_{1\epsilon }\right)G_k{e}^{{\delta }_i^{n+1}}+\Delta t\left(1-C_{2\epsilon }\right)e^{-{\delta }_i^{n+1}}=0,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill {\alpha }_i^{n+1}& =& {\alpha }_i^{**}+\Delta t{C}_{\mu }{G}_k{e}^{{\delta }_i^{n+1}}-\Delta t{e}^{-{\delta }_i^{n+1}},\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill {\beta }_i^{n+1}& =& {\beta }_i^{**}+\Delta t{C}_{1\epsilon }{C}_{\mu }{G}_k{e}^{{\delta }_i^{n+1}}-\Delta t{C}_{2\epsilon }{e}^{-{\delta }_i^{n+1}},\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill k_i^{n+1}& =& exp{\alpha }_i^{n+1}>0,\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill {\epsilon }_i^{n+1}& =& exp{\beta }_i^{n+1}>0,\hfill \end{array}$$

(45)

is positivity-preserving in each fractional step of the scheme, in the sense that if $k_i^n>0$, ${\epsilon }_i^n>0$ for all ${\mathsf{\Omega }}_i\in \mathsf{\Omega }$, then $k_i^{*}>0$, ${\epsilon }_i^{*}>0$, $k_i^{**}>0$, ${\epsilon }_i^{**}>0$ and $k_i^{n+1}>0$, ${\epsilon }_i^{n+1}>0$ for all ${\mathsf{\Omega }}_i\in \mathsf{\Omega }$, under the CFL-type condition

$$-\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|u_{ij}^{-}<1,u_{ij}^{-}=\frac{1}{2}\left(u_{ij}^n-\left|u_{ij}^n\right|\right).$$

(46)

Proof. 

Multiplying the divergence-free property (38) with the factor $\frac{\Delta t}{|{\mathsf{\Omega }}_i|}{q}_i^n$ and subtracting it from (39) leads to

$$q_i^{*}=q_i^n-\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|\left(\frac{1}{2}{u}_{ij}^n\left(q_j^n-q_i^n\right)-\frac{1}{2}\left|u_{ij}^n\right|\left(q_j^n-q_i^n\right)\right),$$

which can be rewritten more compactly as

$$q_i^{*}=q_i^n-\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|u_{ij}^{-}\left(q_j^n-q_i^n\right),$$

using the abbreviation $u_{ij}^{-}=\frac{1}{2}\left(u_{ij}^n-\left|u_{ij}^n\right|\right)\le 0$. Rearranging terms yields the convex combination

$$q_i^{*}=q_i^n\left(1+\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|u_{ij}^{-}\right)+\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|\left(-u_{ij}^{-}\right)q_j^n>0,$$

with $q\in \left\{k,\epsilon \right\}$ since $k_i^n>0$, ${\epsilon }_i^n>0$, $-u_{ij}^{-}\ge 0$ and $\left(1+\frac{\Delta t}{|{\mathsf{\Omega }}_i|}\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|u_{ij}^{-}\right)>0$ due to (46). The implicit diffusion step (40) can be rewritten as

$$|{\mathsf{\Omega }}_i|q_i^{**}-\Delta t\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}|\partial {\mathsf{\Omega }}_{ij}|\left({\mu }_{ij}^n+{\mu }_{t,ij}^n\right)\frac{q_j^{**}-q_i^{**}}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}=\left|{\mathsf{\Omega }}_i\right|q_i^{*},$$

or, more compactly as

$$\sum _j{A}_{ij}^n{q}_j^{**}=\left|{\mathsf{\Omega }}_i\right|q_i^{*},$$

with the diagonal elements of matrix $\mathbf{A}=A_{ij}^n$ given by

$$A_{ii}^n=\left|{\mathsf{\Omega }}_i\right|+\Delta t\sum _{{\mathsf{\Omega }}_j\in {\mathcal{N}}_i}\left({\mu }_{ij}^n+{\mu }_{t,ij}^n\right)\frac{|\partial {\mathsf{\Omega }}_{ij}|}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}>0,$$

and the non-zero off-diagonal elements

$$A_{ij}^n=-\Delta t\left({\mu }_{ij}^n+{\mu }_{t,ij}^n\right)\frac{|\partial {\mathsf{\Omega }}_{ij}|}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}=A_{ji}^n<0,\forall {\mathsf{\Omega }}_j\in {\mathcal{N}}_i.$$

It is obvious that matrix $A_{ij}^n$ is symmetric and positive definite because $A_{ij}^n=A_{ji}^n$ and the matrix is diagonally dominant

$$\left|A_{ii}^n\right|>\sum _{j\ne i}\left|A_{ij}^n\right|.$$

Since all off-diagonal elements are non-positive, matrix $A_{ij}^n$ is a Stieltjes matrix, and thus an M-matrix, whose inverse ${\mathbf{A}}^{-1}$ has only non-negative entries. This guarantees that $q_i^{**}>0$ if $q_i^{*}>0$, with $q\in \left\{k,\epsilon \right\}$. From the results shown in the previous section, we know that (41) has a unique solution ${\delta }_i^{n+1}$ from which ${\alpha }_i^{n+1}$ and ${\beta }_i^{n+1}$ can be computed from (42) and (43). Thanks to (44) and (45), it is obvious that $k_i^{n+1}=exp{\alpha }_i^{n+1}>0$ and ${\epsilon }_i^{n+1}=exp{\beta }_i^{n+1}>0$, which completes the proof. □

4. Numerical Results

To assess the proposed hybrid FV/FE methodology, we run several test cases with both the Cartesian and the unstructured mesh-based algorithms. We first focus on the laminar case, presenting numerical results for both Newtonian and non-Newtonian fluids on unstructured grids and comparing them with available analytical and numerical reference solutions. Secondly, we study the behavior of the methodology for solving turbulent flows, including, e.g., the quite demanding turbulent boundary layer benchmark, in which the entire boundary layer is resolved, including the viscous sublayer. If not stated otherwise, we assume $\rho =1$ in all test cases. As already stated previously, for the sake of simplicity in this section we only consider the two-dimensional planar case, hence assuming $\partial /\partial z=0$ and setting the velocity component in the z-direction to $w=0$.

4.1. Non-Newtonian Flows

Validation of the methodology on unstructured meshes in the framework of laminar Newtonian and non-Newtonian flows is carried out by considering different rheologies for the Couette and Hagen–Poiseuille flows, and for the lid-driven cavity benchmark, similar to what was done in [52] for an alternative mathematical model for the description of non-Newtonian fluids. Furthermore, the flow of a non-Newtonian fluid around a circular cylinder is computed, in order to show the capability of the hybrid FV/FE method to deal also with more complex geometries thanks to the use of unstructured meshes.

4.1.1. Couette Flow

The first set of test cases analysed consist of modified versions of the Couette flow benchmark in $\mathsf{\Omega }={[0,1]}^2$. Apart from the classical problem for Newtonian fluids, which corresponds to $n=1$, $\kappa =1$ and ${\sigma }_y=0$, we also consider the cases $n\in \left\{0.5,1,1.5\right\}$ with yield stress ${\sigma }_y\in \left\{0,0.5\right\}$. Periodic boundary conditions have been set in the $x-$direction, while no-slip walls are defined at $y=0$ and $y=1$. While the bottom is supposed to be steady, we assume that the upper boundary, $y=1$, moves horizontally with velocity $u_C\in {\mathbb{R}}^{+}$, in particular, ten non-equidistant velocities between $0.001$ and 2 are considered. The initial conditions are set to $p(\mathbf{x},t)=0$, $\mathbf{v}(\mathbf{x},0)=\mathbf{0}$, and the final time for all simulations is $t=2$. The exact stationary solution of this simple flow problem is $p\left(\mathbf{x}\right)=0$ and $\mathbf{v}\left(\mathbf{x}\right)={(y{u}_C,0)}^T$ for $t\to \infty$. In Figure 4, we show the comparison of the numerical results obtained on a mesh of 2906 simplex elements and the analytical solutions computed with the Herschel–Bulkley (HB) model. In the left figure, the Couette flow of a fluid without yield stress is plotted, with $\kappa =1$ and power-law indexes $n\in \left\{0.5,1,1.5\right\}$, while in the right plot, we report the results for a fluid with yield stress ${\sigma }_y=0.5$, $\kappa =1$ and $n\in \left\{0.5,1,1.5\right\}$. An excellent agreement for the power-law rheology of the Herschel–Bulkley model is observed in all cases.

4.1.2. Hagen–Poiseuille Flow

As second test, a Hagen–Poiseuille flow of a Herschel–Bulkley (HB) fluid is studied. The computational domain under consideration is $\mathsf{\Omega }={[0,1]}^2$ with periodic boundary conditions in $x-$direction and homogeneous viscous wall boundary conditions on the surfaces $y=0$ and $y=1$. The analytical solution of this flow [48], is given by

$$u\left(y\right)=\left\{\begin{array}{cc}\frac{1}{m}f\left({\left(\frac{y_0}{h}\right)}^m-{\left(\frac{y-y_0}{h}\right)}^m\right),\hfill & \mathrm{if}y\le y_0,\hfill \\ \frac{1}{m}f{\left(\frac{y_0}{h}\right)}^m,\hfill & \mathrm{if}{y}_0\le y\le h-y_0,\hfill \\ \frac{1}{m}f\left({\left(\frac{y_0}{h}\right)}^m-{\left(\frac{y-(h-y_0)}{h}\right)}^m\right),\hfill & \mathrm{if}y>h-y_0,\hfill \end{array}\right.$$

(47)

where $f=\frac{h}{V}{\left(\frac{\Delta ph}{\kappa }\right)}^{\frac{1}{n}}$, $m=1+\frac{1}{n}$, $y_0=h\left(\frac{1}{2}-\frac{\mathrm{Bi}}{f^n}\right)$ with h the channel width, $V=1$ the reference velocity, and $\mathrm{Bi}$ the Bingham number given in Equation (8). Since the computational domain has length $h=1$, the pressure drop, $\Delta p$, is set to 1. To easily impose this condition, we follow [12] and add a source term of the form ${\mathbf{S}}_{\mathbf{v}}={\left(1,0\right)}^T$ in the right-hand side of momentum conservation Equation (2). The initial condition is set to $p=0$, $\mathbf{v}=0$. We carry out a transient simulation that automatically stops when the steady state solution is reached with a tolerance of ${10}^{-10}$. In Figure 5, we compare the numerical solutions computed by using the hybrid FV/FE method and the analytical solutions given by (47) for different rheologies. The left subplot corresponds to a fluid without yield stress, ${\sigma }_y=0$, and power-law indexes $n\in \left\{0.5,1,1.5\right\}$. In the central image, the power-law index has been set to $n=1$, and the yield stress has been chosen to get Bingham numbers from 0 to $0.5$, that is, the yield stress is set to ${\sigma }_y\in \left\{0,0.1,0.2,0.3,0.4,0.5\right\}$. The right subplot depicts the solution for $n=0.5$ and ${\sigma }_y\in \left\{0,0.1,0.2,0.3,0.4\right\}$, which implies that the Bingham number goes from 0 to $0.4$. To run all these simulations, the consistency parameter was set to $\kappa =1$, and a primal simplex mesh of 46496 elements has been used. As it can be observed, an excellent agreement between the computed and the analytical solution is achieved.

4.1.3. Lid-Driven Cavity

We now consider two setups of the classical lid-driven cavity benchmark for the incompressible Navier–Stokes equations: a Herschel–Bulkley fluid with yield stress, ${\sigma }_y\ne 0$, and a power-law fluid without yield stress, ${\sigma }_y=0$. In both cases, the computational domain is $\mathsf{\Omega }={[0,1]}^2$ and the initial conditions are $p(\mathbf{x},0)=0$ and $\mathbf{v}(\mathbf{x},0)=0$. We consider no-slip boundary conditions on all the boundaries except on $y=1$, where the lid velocity is imposed equal to one.

In Figure 6, the results of the simulations with ${\sigma }_y=0$ are shown. The consistency parameter is set to $\kappa ={10}^{-2}$ so using the expression (9), we get a Reynolds number of Re $=100$. To discretize the computational domain, a mesh with 146,826 triangular elements is used. The values for the power-law indexes are $n=0.5$ (top row), $n=1.0$ (central row) and $n=1.5$ (bottom row). In the left column, the velocity norm alongside with the streamlines is displayed, and in the central and right columns, we report the 1D cuts of the velocity fields $u(0.5,y)$ and $v(x,0.5)$ and the reference solution given in [89].

In Figure 7, the results of the simulations with non-zero yield stress and power-law index $n=1$ are shown. The consistency parameter is set to $\kappa ={10}^{-2}$, so the Reynolds number is Re $=100$. To discretize the computational domain, a mesh of 589612 triangular elements is employed. The values of the yield stress are ${\sigma }_y={10}^{-2}$ (top row), ${\sigma }_y={10}^{-1}$ (central row) and ${\sigma }_y=1$ (bottom row). Hence, we are considering $\mathrm{Bi}=1$, $\mathrm{Bi}=10$ and $\mathrm{Bi}=100$, respectively. In the left column, the ratio $\sigma /{\sigma }_y$ is shown, and in the central and right columns, we provide the comparison of the 1D cuts of the velocity fields $u(0.5,y)$ and $v(x,0.5)$ against a reference solution taken again from [89].

4.1.4. Flow around a Cylinder

As the last test case for non-Newtonian fluids, we study the behavior of an incompressible flow around a cylinder. The computational domain under consideration is $\mathsf{\Omega }=[-20,80]\times [-20,20]$ with an embedded circular cylinder of diameter $d=2$ centered at the origin. We first consider a fluid without yield stress with a power-law index $n=1.5$, and then a fluid with yield stress ${\sigma }_y={10}^{-3}$ also with $n=1.5$. The initial velocity and pressure are set to $u(\mathbf{x},0)=0.2$ and $p(\mathbf{x},0)=\frac{1}{\gamma }$, respectively, with $\gamma =1.4$, and $\rho =1$ is the density. The corresponding Reynolds number is Re $\approx 1265$. The first plot in Figure 8 shows the von Kármán vortex street obtained at time $t=250$ for a fluid without yield stress. Meanwhile, the second plot in Figure 8 corresponds to the fluid with ${\sigma }_y={10}^{-3}$. To provide a quantitative reference for this numerical experiment, we compute the Strouhal number, that is the dimensionless frequency of vortex shedding,

$$\mathrm{St}=\frac{fd}{u},$$

where f is the frequency of vortex production, d is the cylinder diameter and $u=0.2$ is the inflow velocity. In both simulations, with and without yield stress, the Strouhal number obtained at time $t=250$ is $\mathrm{St}\approx 0.13$. Figure 9 shows the time series of the velocity $v({\mathbf{x}}_p,t)$ at ${\mathbf{x}}_p=(10,0)$. The simulations have been run for a primal mesh made of 182,165 triangular elements on 1024 CPU cores of the SuperMUC-NG supercomputer.

4.2. Turbulent Flows

In this section, we present several numerical tests aiming at assessing the behavior of the proposed methodology when solving turbulent flows. First, we perform a convergence study on a manufactured test. Next, we address three benchmarks from [23], originally proposed to compute the constants arising in the turbulence model, but that also results in simple test cases with available analytical reference solutions. Once the methodology is validated, a 2D planar mixing layer benchmark is addressed. Finally, to check the capability of the methodology to solve also the flow in the proximity of a wall properly, we present numerical results obtained for the turbulent boundary layer flow over a flat plate using the low Reynolds modification of Lam-Bremhorst [81], see (10)–(12).

4.2.1. Manufactured Solution Test Case

The order of convergence of the proposed methodology for the $k-\epsilon$ model is analyzed using a manufactured solution test with a prescribed analytical solution

$$\begin{array}{c}p=\frac{1}{4}\left(cos\left(2x\right)+cos\left(2y\right)\right),\mathbf{v}=\left(\begin{array}{c}\phantom{-}sin\left(x\right)cos\left(y\right)\\ -cos\left(x\right)sin\left(y\right)\end{array}\right),\\ k={10}^{-3}\left(sin\left(x\right)cos\left(y\right)+2\right),\epsilon ={10}^{-3}\left(cos(-y)sin(-x)+2\right).\end{array}$$

(48)

Since the former expressions do not verify the PDEs (2)–(4), we add a corresponding vector of non-stiff algebraic source terms $\tilde{\mathbf{S}}(\mathbf{x},t)$ to the right-hand side of the Equations (2)–(4) so that (48) becomes a solution of the system. To ensure that the errors obtained are not spoiled due to the approximation of these analytical sources, a fourth order quadrature rule is employed.

We consider the computational domain $\mathsf{\Omega }={\left[0,2\pi \right]}^2$ and assume periodic boundary conditions everywhere. Two different sets of meshes are employed, according to the code used to run the simulations. The mapped triangular primal meshes described in

Table 1. Description of the meshes and time steps used to perform the convergence analysis in triangular primal grids.

Mesh	Elements	Vertices	Dual Elements	$\mathbf{\Delta }\mathit{t}$
$M_1$	512	289	800	$4.0·{10}^{-2}$
$M_2$	2048	1089	3136	$1.0·{10}^{-2}$
$M_3$	8192	4225	12,416	$2.5·{10}^{-3}$
$M_4$	32,768	16,641	49,408	$6.25·{10}^{-4}$
$M_5$	131,072	66,049	202,687	$1.5625·{10}^{-4}$
$M_6$	524,288	263,169	787,456	$3.90625·{10}^{-5}$
$M_7$	2,097,152	1,050,625	3,171,158	$9.765625·{10}^{-6}$

were used for the hybrid FV/FE algorithm for triangular grids. Meanwhile, the quadrilateral grids needed for the Cartesian-based algorithm were defined taking a uniform Cartesian mesh with $N_x=N_y\in \left\{20,40,80,160\right\}$ elements in the x and y direction, respectively. The $L^2$ errors and convergence rates obtained at time $t=0.1$ are reported in

Table 2. Spatial $L^2$ error norms and convergence rates for the MMS test obtained at time $t=1$ using the hybrid FV/FE method for Cartesian grids.

${\mathit{N}}_{\mathit{x}}={\mathit{N}}_{\mathit{y}}$	${\mathit{L}}_{\mathbf{\Omega }}^2\left(\mathit{k}\right)$	$\mathcal{O}\left(\mathit{k}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^2\left(\mathit{\epsilon }\right)$	$\mathcal{O}\left(\mathit{\epsilon }\right)$
20	$1.5130·{10}^{-4}$		$5.1100·{10}^{-4}$
40	$3.5703·{10}^{-5}$	$2.08$	$1.3307·{10}^{-4}$	$1.94$
80	$8.9135·{10}^{-6}$	$2.00$	$3.3596·{10}^{-5}$	$1.99$
160	$2.2305·{10}^{-6}$	$2.00$	$8.4151·{10}^{-6}$	$2.00$

and

Table 3. Spatial $L^2$ error norms and convergence rates for the MMS test obtained at time $t=1$ using the hybrid FV/FE method for unstructured grids.

Mesh	${\mathit{L}}_{\mathbf{\Omega }}^2\left(\mathit{k}\right)$	$\mathcal{O}\left(\mathit{k}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^2\left(\mathit{\epsilon }\right)$	$\mathcal{O}\left(\mathit{\epsilon }\right)$
M1	$3.7640·{10}^{-4}$		$3.6172·{10}^{-4}$
M2	$6.6398·{10}^{-5}$	$2.50$	$9.5886·{10}^{-5}$	$1.92$
M3	$1.1782·{10}^{-5}$	$2.49$	$2.5060·{10}^{-5}$	$1.94$
M4	$2.1802·{10}^{-6}$	$2.43$	$6.3398·{10}^{-6}$	$1.98$
M5	$4.8186·{10}^{-7}$	$2.18$	$1.5959·{10}^{-7}$	$1.99$
M6	$1.1770·{10}^{-7}$	$2.03$	$4.0073·{10}^{-7}$	$1.99$
M7	$2.9464·{10}^{-8}$	$2.00$	$1.0045·{10}^{-7}$	$2.00$

. The expected second order of accuracy is achieved with both, the Cartesian and the unstructured code.

4.2.2. Isotropic Turbulence Decay

The second test case analyses the decay of homogeneous turbulence. Following [82], we set $\mathbf{v}=0$, hence $\nabla \mathbf{v}=0$, so that the $k-\epsilon$ model reduces to the ODE system

$${\partial }_{t}k=-\epsilon ,{\partial }_t\epsilon =-C_{2\epsilon }\frac{{\epsilon }^2}{k},$$

(49)

which is compatible with a polynomial decay solution and can be solved using classical ODE solvers obtaining a reference solution to test the proposed hybrid methodologies. To this aim, we consider the computational domain $\mathsf{\Omega }={\left[-0.5,0.5\right]}^2$ with periodic boundary conditions everywhere and a final simulation time of $t=20$. The initial conditions read

$$p\left(\mathbf{x},0\right)=0,\mathbf{v}\left(\mathbf{x},0\right)=\mathbf{0},k\left(\mathbf{x},0\right)=\epsilon \left(\mathbf{x},0\right)={10}^{-4},$$

and the laminar viscosity is set to $\mu ={10}^{-5}$. The Cartesian algorithm is run using a mesh of $N_x=N_y=50$ divisions along each axis, while the mesh for the unstructured-based code is made of 226 simplex elements. Figure 10 shows an excellent agreement for the time evolution of the kinetic energy and the dissipation rate obtained with both hybrid schemes and the reference solution computed from (49) applying a fourth order Runge–Kutta scheme with fixed time step $\Delta t={10}^{-3}$.

4.2.3. Turbulent Couette Flow

Considering the setup of a Couette flow for the velocity field,

$$\mathbf{v}\left(\mathbf{x},0\right)={\left(Cy,0,0\right)}^T,C\in {\mathbb{R}}^{+},p=0,$$

and constant initial conditions for the turbulence variables, Equations (3) and (4) reduce to the simple nonlinear algebraic system

$$C_{\mu }\frac{k^2}{\epsilon }{C}^2=\epsilon ,C_{1\epsilon }{C}_{\mu }k{C}^2=C_{2\epsilon }\frac{{\epsilon }^2}{k}.$$

(50)

Consequently, the turbulent kinetic energy, the dissipation rate and the velocity field will remain constant in time. To check that this property has been correctly conveyed to the proposed numerical schemes, we define a computational domain $\mathsf{\Omega }={\left[0,1\right]}^2$, set the shear velocity and the laminar viscosity to $C=0.1$, $\mu ={10}^{-5}$, respectively, compute the relations between the model parameters and the turbulence variables to verify (50), and fix their values according to [82] as:

$$k=\frac{1}{C\sqrt{C_{\mu }}}\epsilon ,\epsilon =1.5·{10}^{-6},C_{1\epsilon }=C_{2\epsilon }=1.92,C_{\mu }=0.09.$$

Simulations for single, double, and quadruple machine precision are run for each grid arrangement: an $N_x\times N_y=50\times 50$ Cartesian grid and an unstructured mesh made of 5616 triangular elements. The $L^2-$errors obtained at time $t=10$ are presented in

Table 4. Spatial $L^2-$error norms obtained for the turbulent Couette flow at time $t=10$ using the hybrid FV/FE method for Cartesian grids and for unstructured simplex meshes.

Hybrid FV/FE Scheme on Cartesian Grid
	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{u}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{v}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{k}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{\epsilon }\right)$
Single precision	$1.2205·{10}^{-8}$	$3.0451·{10}^{-9}$	$1.8644·{10}^{-11}$	$3.6947·{10}^{-13}$
Double precision	$2.4610·{10}^{-18}$	$2.1039·{10}^{-18}$	$4.0658·{10}^{-20}$	$4.2352·{10}^{-22}$
Quadruple precision	$2.4439·{10}^{-35}$	$5.1321·{10}^{-36}$	$3.5278·{10}^{-38}$	$2.0443·{10}^{-40}$
Hybrid FV/FE Scheme on Unstructured Triangular Mesh
	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{u}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{v}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{k}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^{\mathbf{2}}\left(\mathbf{\epsilon }\right)$
Single precision	$2.2135·{10}^{-6}$	$7.0541·{10}^{-8}$	$1.7449·{10}^{-10}$	$3.3065·{10}^{-10}$
Double precision	$1.9774·{10}^{-17}$	$4.4820·{10}^{-19}$	$4.0658·{10}^{-20}$	$4.2352·{10}^{-22}$
Quadruple precision	$5.3593·{10}^{-36}$	$9.7438·{10}^{-36}$	$7.5226·{10}^{-35}$	$2.2563·{10}^{-36}$

. Machine precision is obtained for all cases, demonstrating the good behavior of the implicit approach proposed to solve the production terms of the $k-\epsilon$ equations.

4.2.4. Logarithmic Velocity Profile

As already mentioned in Section 2, the boundary layer in the vicinity of a wall can be divided into the viscous, buffer, and logarithmic layers. While the use of wall laws may avoid the necessity of simulating the flow within the inner layer, the log layer, where the first mesh point is expected to be located, still needs to be correctly captured. Hence, considering future extensions of the hybrid methodology to real-world applications, either using a low Reynolds number model or wall laws, needs for a good performance of the proposed methodology when dealing with the logarithmic profile. Therefore, even if solving such complex flows is not the objective of this paper, to assess the proposed methodology, we run a logarithmic profile test case in the computational domain $\mathsf{\Omega }=\left[0,1\right]\times \left[4·{10}^{-5},4·{10}^{-3}\right]$, assuming a flow parallel to the horizontal axis, with periodic conditions in $x-$direction, and Dirichlet boundary conditions in the bottom and top boundaries. The velocity and length scales are then defined as

$$u^{*}=\sqrt{\frac{\mu }{\rho }\frac{\partial u}{\partial y}},y^{*}=\frac{\mu }{\rho u^{*}},\mathrm{so}{u}^{+}=\frac{u}{u^{*}},y^{+}=\frac{y}{y^{*}}.$$

(51)

From experimental results, we know that for $20\le y^{+}\le 100$ we are in the log layer and

$$u^{+}=\frac{1}{\chi }{y}^{+}+5.5,\chi =0.41.$$

(52)

Moreover, the mean flow is supposed to be stationary and $G_k={\mu }_t{\left({\partial }_{y}u\right)}^2$ so Equations (3) and (4) for $\rho =1$ become

$$\begin{array}{c}-\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right]-G_k+\rho \epsilon =0,\end{array}$$

(53)

$$\begin{array}{c}-\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]-C_{1\epsilon }\frac{\epsilon }{k}{G}_k+C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}=0.\end{array}$$

(54)

Neglecting the laminar viscosity, i.e., setting $\mu =0$, we get that

$$k=\frac{u^{{*}^2}}{\sqrt{C_{\mu }}},\epsilon =\frac{u^{{*}^3}}{\chi y},{\left({\partial }_{y}u\right)}^2=\frac{u^{{*}^2}}{{\chi }^2{y}^2},$$

(55)

is the solution of (53) and (54) given that

$${\sigma }_{\epsilon }=\frac{{\chi }^2}{\sqrt{C_{\mu }}\left(C_{2\epsilon }-C_{1\epsilon }\right)},$$

see [82] for further details. Let us note that, from (55), the turbulent viscosity, ${\mu }_t$, is a linear function in y.

Gathering (51), (52), and (55), we observe that the initial condition for this test case reads

$$\begin{array}{c}u\left(\mathbf{x},0\right)=u^{*}\left(\frac{1}{\chi }ln\frac{y}{y^{*}}+5.5\right),v\left(\mathbf{x},0\right)=0,k\left(\mathbf{x},0\right)=\frac{u^{{*}^2}}{\sqrt{C_{\mu }}},\epsilon \left(\mathbf{x},0\right)=\frac{u^{{*}^3}}{\chi y}.\end{array}$$

The final setup is completed defining

$$\begin{array}{c}{u}^{*}={10}^{-2},y^{*}={10}^{-6},\\ C_{\mu }=0.09,C_{1\epsilon }=1.44,C_{2\epsilon }=1.92,{\sigma }_k=1.0,{\sigma }_{\epsilon }=1.16736111111111111.\end{array}$$

The performance of both the Cartesian and the unstructured hybrid FV/FE schemes is studied. To build the Cartesian grid, we take advantage of the possibility of including a bias on the grid to get a boundary layer mesh. In particular, we define $N_x=10$ divisions along the x-axis and $N_y=100$ in $y-$direction with a refinement factor of $1.05$ towards the bottom boundary. Meanwhile, the triangular grid has a total number of 40,000 uniform elements. The computed velocity, mean dissipation rate, and turbulent viscosity profiles at $t=10$ are depicted in Figure 11. We observe a good agreement between the numerical solutions obtained with the hybrid schemes and the exact solution. Indeed, the $L_{\mathsf{\Omega }}^2-$error norms obtained for this test case at the final time with the Cartesian hybrid FV/FE scheme were $2.8713·{10}^{-6}$ for the velocity component u, $1.0098·{10}^{-8}$ for k, and $3.5296·{10}^{-6}$ for $\epsilon$. Meanwhile, with the unstructured hybrid FV/FE scheme we got $L_{\mathsf{\Omega }}^2\left({\mathbf{v}}_1\right)=1.6335·{10}^{-5}$, $L_{\mathsf{\Omega }}^2\left(k\right)=2.7570·{10}^{-7}$ and $L_{\mathsf{\Omega }}^2\left(\epsilon \right)=1.2719·{10}^{-5}$. As expected, more significant errors are obtained in the latter simulation due to the relatively coarse mesh used in the vicinity of the bottom boundary compared to the height of the quadrilateral elements of the non-uniform Cartesian grid located in the same region.

4.2.5. Turbulent Planar Shear Layer

The fifth test considered is the turbulent planar shear layer benchmark, also known as the turbulent mixing layer. It consists of two adjacent layers of fluid, both of them traveling in the direction parallel to its interface but with different velocities leading to the generation of a stationary shear layer in between. We define the computational domain $\mathsf{\Omega }=\left[0,1\right]\times \left[-0.25,0.25\right]$ and consider an initial condition of the form

$$\begin{array}{c}u(\mathbf{x},0)=\frac{1}{2}\left(u_R+u_L\right)+\frac{1}{2}\left(u_R-u_L\right)\mathrm{erf}\left(\frac{y}{\zeta }\right),\zeta ={10}^{-2},u_L=1,u_R=2,\\ v(\mathbf{x},0)=0,p(\mathbf{x},0)=0,k(\mathbf{x},0)={10}^{-3},\epsilon (\mathbf{x},0)=2.5·{10}^{-3},\end{array}$$

where an error function has been introduced to smooth the discontinuity of the initial data. Since this test does not have a known analytical solution, a mesh convergence study has been carried on. To this end, we have considered the Cartesian hybrid FV/FE scheme and three different grids with $\left(N_x,N_y\right)\in \left\{\left(50,200\right),\left(100,500\right),\left(200,1000\right)\right\}$ divisions in the horizontal and vertical directions, respectively. The horizontal velocity, dissipation rate, and turbulent kinetic energy profiles obtained at $x=0.8$ are reported in Figure 12. The results correspond to $t=10$, when the stationary solution has already been reached. We notice that mesh convergence is achieved. Moreover, we also observe a good agreement with the solution obtained with the unstructured hybrid FV/FE scheme run on a primal grid made of ${10}^5$ simplex elements. To illustrate the 2D flow pattern, the contour plot of the turbulent viscosity field is provided in Figure 13 for the hybrid FV/FE scheme on both, the Cartesian grid and on the unstructured simplex mesh.

4.2.6. Turbulent Boundary Layer over a Flat Plate

The turbulent boundary layer over a flat plate is an essential test case for the validation of numerical methods for turbulence models in near-wall regions. Since all the phenomena to be analyzed take place within the region $0\le y^{+}\le {10}^3$, the main issue when addressing this problem is the high mesh resolution needed in the vicinity of the wall. In fact, we know that the viscous sublayer arrives up to around five wall units, where the flow starts changing its behavior as viscous effects lose importance with respect to the inertial effects, which dominate the logarithmic layer starting around $y^{+}=10$. To be able to design a fine enough mesh with a good aspect ratio, we run this test case using the Cartesian hybrid FV/FE scheme, which has also the advantage of having the elements parallel to the flow direction. It is indeed common practice to use quadrilateral grids in the vicinity of the wall to resolve the boundary layer properly, while unstructured meshes can be used outside the boundary layer. Accordingly, to discretize the computational domain $\mathsf{\Omega }=\left[0,2\right]\times \left[0,1\right]$, we consider 100 divisions along the horizontal axis and 120 in the vertical direction, with a refinement factor of 1.1 towards the wall located at $y=0$. At this boundary, homogeneous Dirichlet boundary conditions are imposed for the velocity field, the turbulent kinetic energy is set to a minimum value $k={10}^{-14}$, and Neumann boundary conditions are defined for the pressure and energy dissipation rate. A horizontal inflow velocity of $\mathbf{v}=(1,0)$ is imposed at the left boundary of the domain, where k and $\epsilon$ are set to their reference values to be computed from the turbulence intensity, $I_t$, and the turbulent Reynolds number, $R{e}_t$, as

$$k_0=\frac{3}{2}{I}_t^2,{\epsilon }_0=\frac{k_0^2}{\mu R{e}_t},$$

where we have taken $I_t={10}^{-2}$, $R{e}_t={10}^3$, $\mu =1/Re$ with $Re=11.1·{10}^6$ the global Reynolds number of the flow with respect to a reference velocity of unity and a unitary reference length scale. Pressure outflow boundary conditions are employed at the top and right boundaries, that is, we set the pressure, $p=0$, while Neumann conditions are considered for the remaining unknowns. Finally, the initial condition for the horizontal velocity component u corresponds to the solution of the laminar Blasius boundary layer, see, e.g., [90], which is computed using a fourth order Runge–Kutta scheme in combination with the shooting method. The vertical velocity and the pressure are initialized to zero and $k(\mathbf{x},0)=k_0$, $\epsilon (\mathbf{x},0)={\epsilon }_0$. Let us note that to run this test case, the low Reynolds $k-\epsilon$ turbulence model, (10)–(12), is used instead of its standard version so that we can compute the solution right to the wall [81,91].

In Figure 14, we show the usual profile of $u^{+}$ as a function of $y^{+}$ obtained at $x=1$ and $t=10$. To provide a direct comparison with available reference solutions for the test case, also the law of the wall is plotted, together with law of the wake profile proposed by Coles in [92]. A pretty good agreement is observed between the hybrid FV/FE solution and the wall law given by

$$u^{+}=\left\{\begin{array}{cc}{y}^{+},\hfill & y^{+}<11.44531911,\hfill \\ \frac{1}{\chi }ln\left(y^{+}\right)+5.5\hfill & y^{+}\ge 11.44531911.\hfill \end{array}\right.$$

For the definitions of $u^{+}$ and $y^{+}$, see Equation (51), and the von Kármán constant is chosen as usual as $\chi =0.41$. Moreover, the solution is also very close to the law of the wake of Coles, which, apart from the profile above the log law region, also provides the behavior of the flow in the intermediate buffer region. Finally, our numerical results obtained with the new hybrid FV/FE scheme are also compared against a numerical reference solution obtained using the semi-implicit finite volume scheme proposed in [9,12]. An excellent agreement is observed between both numerical results.

5. Conclusions

In this paper, we have introduced a new family of semi-implicit hybrid finite volume/finite element methods on edge-based staggered simplex meshes and on edge-based staggered Cartesian grids for the numerical solution of the Reynolds-Averaged Navier–Stokes (RANS) equations based on the $k-\epsilon$ turbulence model and including non-Newtonian power-law fluids with yield stress (Herschel–Bulkley fluids). The coupled governing PDE system, which contains hyperbolic, parabolic, and elliptic terms, as well as stiff algebraic source terms, is then discretized at the aid of operator splitting, in order to divide the complete problem into a set of simpler subproblems, each of which is then discretized with the most appropriate method.

For the discretization of the hyperbolic convection terms, we use classical Godunov-type finite volume schemes of order one or two. The viscous terms require an implicit time discretization, since the viscosity coefficient of power-law fluids with $n<1$ and of all non-Newtonian fluids with yield stress ${\sigma }_y>0$ tends to infinity when the shear rate tends to zero. In general unstructured meshes, we use implicit continuous ${\mathcal{P}}^1$ finite elements for the discretization of the viscous terms, while implicit finite volume schemes are used for the Cartesian case. The elliptic pressure Poisson equation is solved at the aid of continuous ${\mathcal{P}}^1$ and ${\mathcal{Q}}^1$ finite elements on simplex meshes and on Cartesian grids, respectively. For unstructured simplex meshes, this paper shows a special interpolation between the discrete solution of the finite volume scheme on the edge-based dual mesh and the vertices of the finite element mesh, which is exactly reversible, i.e., the interpolation of the data from the vertices of the FE mesh to the FV mesh and back to the vertices of the FE mesh gives exactly the identity operator, hence no spurious numerical errors and numerical dissipation are introduced by the chosen interpolation between the FV and the FE mesh.

Last but not least, a new and provably positivity-preserving discretization of the stiff algebraic source terms of the $k-\epsilon$ model has been proposed, based on the implicit discretization of a linear combination of the evolution equations of the logarithms of k and $\epsilon$. In addition to the positivity, which is trivial to show for the chosen discretization, it is also possible to prove the existence and uniqueness of the discrete solution of the ODE subsystem. Last but not least, in this paper we have also rigorously proven the positivity preserving property of a full semi-implicit finite volume scheme applied to the $k-\epsilon$ model on orthogonal unstructured meshes, which contains nonuniform Cartesian grids as a special case.

The proposed method has been applied to several academic benchmark problems for non-Newtonian flows and for simple planar turbulent flows in two space dimensions and was compared with existing exact or numerical reference solutions. In all cases, a very good agreement between our numerical results and the reference solutions has been obtained. In the future we plan to extend our method also to more realistic and more complex three-dimensional turbulent flows including chemical reactions, making use of the massively parallel implementation of the hybrid FV/FE algorithm recently documented in [31]. However, the simulation of complex turbulent flows in 3D would require the implementation of suitable wall functions in order to reduce the computational effort of a 3D simulation and also needs a careful validation against available experimental results, which is out of scope of the present paper.

Further work will concern the application of the new method to non-Newtonian flows in the human cardiovascular system, coupled with the one-dimensional network models used in [93,94,95,96]. We will also consider the extension of our provably positivity-preserving scheme to other turbulence models, such as the $k-\omega$ model and to more complex full Reynolds stress models, in particular to the ones recently introduced in [97,98,99,100] for shallow water turbulence.