A high-order discontinuous Galerkin solver for the incompressible RANS equations coupled to the k− turbulence model

Accurate methods to solve the Reynolds-Averaged Navier-Stokes (RANS) equations coupled to turbulence models are still of great interest, as this is often the only computationally feasible approach to simulate complex turbulent ﬂows in large engineering applications. In this work, we present a novel discontinuous Galerkin (DG) solver for the RANS equations coupled to the k − (cid:2) model (in logarithmic form, to ensure positivity of the turbulence quantities). We investigate the possibility of modeling walls with a wall func- tion approach in combination with DG. The solver features an algebraic pressure correction scheme to solve the coupled RANS system, implicit backward differentiation formulae for time discretization, and adopts the Symmetric Interior Penalty method and the Lax-Friedrichs ﬂux to discretize diffusive and convective terms respectively. We pay special attention to the choice of polynomial order for any transported scalar quantity and show it has to be the same as the pressure order to avoid numerical instability. A manufactured solution is used to verify that the solution converges with the expected order of accuracy in space and time. We then simulate a stationary ﬂow over a backward-facing step and a Von Kármán vortex street in the wake of a square cylinder to validate our approach.


Introduction
Engineering applications often require the simulation of complex turbulent flows via accurate Computational Fluid Dynamics (CFD) methods. Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) constitute superior approaches in this regard, as they are able to resolve the stochastic fluctuations (though only the large-scale ones in case of LES) of any turbulent flow quantity [1] . However, nowadays they are still very computationally expensive and often unaffordable f or large engineering applications. For this reason, the Reynolds-Averaged Navier-Stokes (RANS) equations coupled to turbulence closure models often remain the preferred approach, even if it only allows for the modeling of the time-averaged flow quantities [2] . Accurate and efficient numerical methods to solve the RANS equations are therefore still of great relevance.
In this perspective, discontinuous Galerkin Finite Element methods (DG-FEM) are particularly attractive, due to their flexibility, high accuracy, and robustness. The characteristic feature of DG is that unknown quantities are approximated with polynomial ba-sis functions discontinuous across the mesh elements. This requires the use of numerical fluxes to deal with the discretization of face integrals, as in Finite Volume methods. However, thanks to the lack of continuity constraints, conservation laws are satisfied locally on each element, and the resulting algorithm stencil is compact, making the method suitable for efficient parallelization. As in continuous Galerkin FEM, the accuracy of the solution can be improved by increasing the order of the polynomial discretization. Moreover, DG methods can easily handle complex geometries, structured or unstructured meshes, and non-conformal local mesh and/or order refinement.
This class of Finite Element methods was originally developed in the early '70s to solve radiation transport problems [3] . However, it has become increasingly popular for CFD applications only in the past three decades, after the development of effective DG schemes for hyperbolic and elliptic problems. We refer to the reviews by Cockburn and Shu [4] and Arnold et al. [5] for a complete overview of these methods.
Nowadays, substantial experience has been gained in the DG-FEM discretization of the incompressible Navier-Stokes equations, and a variety of approaches can be found in literature. An early research effort in the field is the work of Cockburn et al. [6] , who proposed a locally conservative Local DG (LDG) method for the incompressible Oseen equations. The authors later extended the https://doi.org/10.1016/j.compfluid.2020.104710 0045-7930/© 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) approach to the full Navier-Stokes system [7] , employing a postprocessing operator to obtain an exactly divergence-free convective velocity field. Van der Vegt and Sudirham [8] presented an original DG scheme for the solution of the Oseen equations, where discontinuous basis functions are also used to discretize the time domain. The method, called space-time DG, is particularly suitable for simulations involving deforming or moving meshes. More recently, Rhebergen et al. [9] extended the approach to the full Navier-Stokes system. In [10] , the incompressibility constraint is relaxed by introducing an artificial compressibility term in the continuity equation. Thus, the numerical flux for the inviscid part of the Navier-Stokes equations can be computed by solving a Riemann problem just like for compressible flows. The coupled system of equations is solved by means of a Newton method with an implicit Euler time scheme. Higher-order Runge-Kutta schemes were instead employed by Bassi et al. [11] . Splitting methods based on pressure or velocity correction approaches have also been proposed in the context of DG solvers (e.g., [12][13][14][15] ). In particular, the solver presented in this paper is based on the work of Shahbazi et al. [16] , where a second-order accurate pressurecorrection scheme is applied at algebraic level. The Symmetric Interior Penalty (SIP) method is used to discretize the diffusive operator, and the Lax-Friedrichs flux is chosen for the hyperbolic term. Similar solvers, but based on a SIMPLE algorithm, were proposed by Klein et al. [17] and [18] .
Regarding turbulent flow simulations, DG solvers have been developed mostly for the compressible RANS equations coupled to either the Spalart-Allmaras (e.g., [19][20][21] ) or the k − ω models (e.g., [22][23][24][25] ). For the latter model, it is customary to solve for the logarithm of ω and to impose a realizability condition on it, to ensure the positivity of the solution and enhance the robustness of the solver. The shear-stress transport model was instead considered in [26] , where particular attention was given to the development of a stabilized continuous FE discretization of the eikonal equation for the computation of the distance to the nearest wall. Far fewer applications of DG schemes to incompressible turbulent flows can be found in literature. The work of Crivellini et al. [27] was the first. They extended the artificial compressibility flux method of [10] to deal with the set of incompressible RANS equations coupled to the SA model. Particular attention was paid to the treatment of negative values of the turbulent viscosity, thus increasing the robustness of the method. The approach was later tested on complex three-dimensional flows [28] and extended to the k − ω model [29] and to high-order Runge-Kutta time schemes [30,31] . More recently, Krank et al. [32] presented a DG solver for the incompressible RANS equations coupled to the SA model based on a semi-explicit velocity-correction splitting scheme augmented with an explicit step for the turbulence equation.
In this work, we extend a solver for laminar flows presented in [33] to handle turbulence through the set of incompressible RANS equations coupled to the k − model. We employ a pressure correction scheme to solve the RANS equations. Thus, contrary to most of previous literature, we do not rely on a free artificial compressibility parameter whose optimal value is problem-based. Moreover, our time scheme is fully implicit (as it relies on backward differentiation formulae) also for the turbulence equations. It constitutes the first step towards the development of a coupled CFD-neutronics solver for large multi-physics problems, such as transient scenarios in molten salt nuclear reactors. In these systems, the flow Reynolds number is of the order of 10 5 or higher, and resolving the steep gradients that characterize the velocity profile close to wall boundaries requires massive computational resources.
To the best of our knowledge, all previous literature on RANS DG solvers deals with wall-resolved turbulent flow simulations. An approach for efficient wall modeling in the DG context was pro-posed by Krank et al. [32] , who solved the RANS equations on coarse grids by enriching the polynomial function space for the velocity in the first layer of boundary elements. In this way, they could impose no-slip conditions at the walls and resolve the sharp solution gradients. However, the k − model is well-known to perform poorly in the vicinity of a wall. For this reason, in this work we investigate the possibility of using a more standard wall function approach [34] to bridge the gap between the viscous layer and the log layer.
The remainder of the paper is organized as follows. In Section 2 , we present the governing equations and the boundary conditions that close them. We then describe our spatial and temporal discretization scheme in Sections 3 and 4 , and we provide details on the solver in Section 5 . In Section 6 , we focus on the choice of polynomial order for scalar quantities, which differs slightly from what previously proposed in literature for similar DG solvers. In Section 7 , we verify our numerical scheme with a manufactured solution and we assess the soundness of our approach simulating the steady turbulent flow over a backward-facing step and a Von Kármán vortex street in the wake of a square cylinder, finally drawing some conclusions in Section 8 .

Governing equations
The RANS model equations for incompressible flows read (we omit dependencies for clarity) [2] ∂ u ∂t Here, u is the velocity and f represents a known momentum source. The total stress tensor is where ν and ν t are the molecular and the turbulent (or eddy) kinematic viscosities respectively. Finally, p = ˜ p ρ + 2 3 k is a pseudo kinematic pressure, where ˜ p is the fluid pressure, ρ the density, and k the turbulent kinetic energy.
The system is closed by adopting the k − turbulence model [34] , where indicates the dissipation rate of k . To ensure positivity of the turbulence quantities (which might not be preserved by the numerical scheme), the model equations are cast into a logarithmic form [35] , which is completely equivalent to the original model. Therefore, we solve for and the model equations read (omitting dependencies) where P k (u ) = ∇u : ∇u + (∇u ) T models the shear production of turbulent kinetic energy. We used the standard values of the model constants proposed by Launder and Spalding [34] and reported in Table 1 . The eddy viscosity can be then computed as Clearly, the latter expression ensures ν t is always positive. Moreover, potentially troublesome divisions by small values of (or by k in some terms of the turbulence equations) are avoided. Finally, q K and q E are extra source terms that, together with f , are used in the framework of the method of manufactured solutions to verify the numerical scheme (see Section 7.1 ). Normally, they are set to zero.

Boundary and initial conditions
The coupled system of the RANS (1a-1b ) and turbulence (4a-4b) equations can be solved only after imposition of proper boundary and initial conditions. The latter consist in imposing suitable values for the initial velocity, pressure, and turbulence fields, which we describe in detail for each test case reported in Section 7 . For the former, four types of boundary are considered in this work, and we describe them in the following.

Inflow boundary
On inflow boundaries, Dirichlet conditions are prescribed for the fluid velocity and the turbulence quantities:

Outflow boundary
On outflow boundaries we prescribe a zero-traction condition for the momentum equation and homogeneous Neumann conditions for the turbulence equations: [ −pI + τ ] · n = 0 , n · ∇K = 0 , n · ∇E = 0 .

(7)
Here, n is the outward unit normal vector at the boundary, and I represents the identity tensor.

Wall boundary
The standard k − model cannot be used to resolve turbulent flows up to wall boundaries, in the vicinity of which viscous effects are dominant and turbulence is damped. For this reason, we use the well-established wall function approach to describe the solution in these areas [34,36] .
The computational boundary is shifted by a distance δ w from the actual physical wall, so that the first degrees of freedom of the solution are placed in the region of fully turbulent flow [37] . The gap is bridged by analytical functions that are supposed to be universal for all wall-bounded flows. This leads to a great gain in computational efficiency, in particular for high-Reynolds flows, because mesh elements are not accumulated in the viscous layer to resolve the steep gradients of the solution.
In this work, we adopt the two-velocity scale wall function described by Ignat et al. [38] . Let u t be the velocity component tangential to the boundary and y the normal distance to the wall. They are made non-dimensional as follows: where u τ = √ τ w /ρ is the friction velocity ( τ w is the shear stress at the wall), whereas u k is a second velocity scale computed from k : where k w is the value of the turbulent kinetic energy at the boundary. The wall function is an analytical relation between u + and y + : where κ = 0 . 41 is the Von Kármán constant. Only smooth walls are considered in this work, so the roughness parameter is E = 9 . 0 . Eqs. (9) and (10) are used to prescribe a homogeneous Robinlike boundary condition for the momentum equation in the tangential direction. In fact, a linear relation exists between the shear stress and the tangential velocity component: where t is a unit vector that is tangential to the boundary and u − (u · n ) n = u t t . The boundary condition (11) is supplemented by a no-penetration condition in the normal direction: Regarding the turbulent quantities, homogeneous Neumann and Dirichlet boundary conditions are prescribed for k and respectively: The value of δ w must be sufficiently large for the wall function (10) to be valid, that is, the first degrees of freedom of the FEM solution must be in the logarithmic layer. This parameter can be either estimated as the distance from the wall of the centroid of the first boundary element or fixed to a predetermined value [39] . This work features standard and predictable flow configurations, so we adopted the second approach.

Symmetry boundary
On symmetry boundaries, we prescribe a homogeneous Neumann condition for both turbulence quantities: For the momentum equation, we prescribe vanishing normal velocity and tangential stress:

Spatial discretization
In this section, we describe in detail the spatial discretization of the governing equations with the discontinuous Galerkin Finite Element method. Some notation and definitions must be introduced first.
Let be the computational domain and its boundary. We indicate with D , N , W , and S the non-overlapping portions of boundary where Dirichlet (i.e., inflow), Neumann (i.e., outflow), wall, and symmetry conditions are imposed, such that The domain is meshed into a set of non-overlapping elements T h . The set of their interior faces is indicated with F i , whereas F D , F N , F W , and F S indicate the sets of Dirichlet, Neumann, wall, and symmetry boundary faces, which discretize the respective portions of boundary. Given an element T ∈ T h , the set of its faces is denoted with F T . The neighbors of face F are grouped into T F h . Any interior face F ∈ F i is equipped with a unit normal vector n F with an arbitrary but fixed direction, while for boundary faces n F coincides with n , which is the unit normal vector outward to .
Any scalar unknown is approximated on T h with a set of basis functions that are continuous within each element, but discontinuous at the element interfaces. We use a hierarchical set of orthogonal modal basis functions: the solution space within an element is the span of all polynomials up to an order P, with P ≥ 1 . For a generic unknown φ, we denote its FEM approximation by φ h and its polynomial order by P φ , so that the solution space can be written as where P P φ is the set of polynomials up to order P φ . The velocity FEM approximation lies in the corresponding vector space V d h,u , where d is the problem dimensionality. Despite not being a requirement of the numerical method, P φ is the same on each element and the meshes are conforming for all the test cases analyzed in this work.
for a point r on an interior face F ∈ F i . At the boundary, the jump and average operators coincide with the unique trace of φ h .

RANS Discretization
The semi-discrete weak formulation of System (1a-1b) subjected to the boundary conditions described in Section 2.1 is obtained by substituting ( u , p ) with their DG-FEM approximations and by subsequently integrating over the entire domain: and, letting f h be the Galerkin projection of All other terms are described more in detail in the following.

Convective term
The convective term of the momentum equation is discretized as where β is the convective field (i.e., the velocity), and H F is a numerical flux function defined on the internal face F . In this work, the Lax-Friedrichs flux is used [40] : where α F is evaluated point-wise along face F according to with = 2 for the momentum equation.
Note that usually ( n F · β D ) < 0 on any F ∈ F D , and therefore the corresponding term in Eq. (21) drops out. However, we included this term because it is relevant for the Taylor-vortex-like manufactured solution in Section 7.1 , where we impose Dirichlet conditions on the entire boundary, as it is usually done for this class of problems (e.g., [ 16,41 ]).

Diffusive term discretization
To discretize the diffusive term, we use a generalization of the Symmetric Interior Penalty (SIP) method [ 42 , p. 122]. Among the several methods available for elliptic problems [5] , we opted for the SIP because it is consistent and adjoint consistent, which guarantees optimal convergence rates for any order of the polynomial discretization. Moreover, it is characterized by compact stencil size, as the degrees of freedom on each element are coupled only with those on the first neighbors, with positive impact on memory requirements, computational cost, and parallelization efficiency.
The following bilinear operator and right-hand side term can be defined: and where L diff (w ) is a linear operator that coincides with the stress tensor τ ( Eq. (2) ) for the momentum equation, that is, L diff (u ) = τ . Note that the effective viscosity (i.e., ν + ν t ) is variable in space and so ∇ · τ = (ν + ν t ) ∇ 2 u , which is why the regular SIP cannot be used here.
The penalty parameter, η F , must be sufficiently large to ensure stability without impacting too negatively the condition number of the system. We follow the optimal value prescriptions provided by Hillewaert [43] and Drosson and Hillewaert [44] in case of anisotropic meshes and highly variable viscosity (as when turbulent flows are resolved): where card (F T ) indicates the number of faces of element T ; C P,T is a factor depending on the polynomial order of the solution and the mesh element type, for quadrilaterals and hexahedra L T is a length scale defined as where · leb represents the Lebesgue measure of a geometrical entity (i.e., length, area, or volume depending on the dimensionality), and f = 2 for F ∈ F i and f = 1 for boundary faces; finally, D is a scale depending on the diffusion parameter D (i.e., the effective viscosity for the momentum equation), which we evaluate pointwise along the face as The contribution of wall and symmetry faces must still be specified to complete the definition of the SIP bilinear form (25) . We extend what was proposed by Cliffe et al. [40] to arbitrarily oriented symmetry boundaries and to walls, where a non-zero shear stress is specified in the tangential direction, leading to a Robinlike boundary condition term (see, for example, [ 42 , p. 127]): where u k /u + is the proportionality constant between tangential stress and velocity component ( Eq. (11) ).

Continuity terms
The discretized continuity equation is constituted by the following discrete divergence operator and right-hand side term [7,16] :

Discretization of turbulence equations
The spatial discretization of the turbulence Eq. (4a-4b) , subjected to the boundary conditions described in Section 2.1 , proceeds along the same line of what was described for the momentum equation. The semi-discrete weak formulations are The temporal and right-hand side terms are straightforward extension to the case of scalar unknowns of those described for the momentum equation ( Eqs. (19) and (20) ). All terms on the righthand side of Eqs. (4a) and (4b) are treated explicitly in time and therefore added to the source integrals in l( The convective terms are discretized as for the momentum equation, though the α F parameter in the Lax-Friedrichs flux ( Eq. (24) ) is evaluated with = 1 . Diffusive terms are treated with the SIP method as in Eqs. (25) and (26) , where the linear operator These are also the diffusion coefficients used to evaluate the penalty parameter according to Eq. (27) . Finally, wall faces contribute to the SIP bilinear form of the E equation as other Dirichlet boundaries, whereas for the K equation their contribution is null, because homogeneous Neumann conditions are imposed on these boundaries. Likewise, for both equations symmetry faces have a null contribution to the SIP bilinear form.

Temporal discretization
Time derivatives in weak forms ( 18a-18b ) and ( 34a-34b ) are discretized implicitly using backward differentiation formulae (BDF) of order M [14,16,18] . For the generic unknown quantity φ and for a constant time step size t , it is ∂φ ∂t where n + 1 indicates the new time step and the BDF coefficients γ are reported in Table 2 . We use up to second-order schemes.
The discretized equations are solved in a segregated way, so that the discrete solution vectors at the new time step, u , p , K , E n +1 , are found with the following algorithm: 1. Obtain predictors for all unknowns at the new time step with an M th -order extrapolation from the previous time steps: where the weights ( ζ j ) are reported in Table 3 . A predictor for the eddy viscosity, ν * t , can be calculated from K * and E * according to Eq. (5) ; 2. Solve the RANS system (1) in a segregated way with a secondorder accurate algebraic splitting scheme described in detail in Section 4.1 . Use the predictors for the velocity, the turbulent kinetic energy 1 , and the eddy viscosity to linearize all terms; 3. Solve the K Eq. (4a) , using u n +1 h as convective field and ν * t ; 4. Solve the E Eq. (4b) , using u n +1 h as convective field and ν * t ; 5. If necessary, update all predictors with the new solutions and iterate steps 2 to 4.
When the BDF2 scheme is used, the time extrapolation of all quantities performed at the beginning of the time step guarantees second-order time convergence even when non-linearities are unresolved (as proven in Section 7.1.1 ). During particularly violent portions of transient (as it might happen at the start-up of a pseudo-transient towards steady-state conditions), iterations between RANS and turbulence equations (typically 2 or 3) are usually performed to better resolve non-linearities, thus avoiding too severe restrictions on the time step size.

Algebraic splitting scheme
The coupled momentum and continuity equations are solved in a segregated way with a second-order accurate pressure correction method [45] . Following [16] , we perform the splitting at algebraic level, which does not require the imposition of artificial pressure boundary conditions [46] . The fully discretized weak form of the RANS system (18) can be cast into the following linear form: where M is the mass matrix, N contains the implicit parts deriving from the discretization of the convective (21) and diffusive (25) terms, and D is the discrete divergence operator (32) . The right-hand side vector b u collects all the known terms in the momentum equation (i.e., boundary conditions, external forces, explicit terms from the discretization of the time derivative), while b p represents the fully-discrete right-hand side of the continuity Eq. (33) . As explained previously, predictors for the velocity field and the eddy viscosity at the new time step are used to linearize the convective and diffusive terms in N . The coupled System (38) is solved in a segregated way by replacing (and approximating) the left-hand-side matrix with its incomplete block LU factorization [46] , where I is the identity matrix, and by introducing δ p n +1 = p n +1 − p n , such that The discrete velocity and pressure fields at the new time step can therefore be calculated with the following predictor-corrector algorithm: 1. Find an approximate velocity field ˜ u n +1 by solving p n +1 = δ p n +1 + p n ; 3. Correct the velocity field, so that it satisfies the discrete continuity equation

Solution of linear systems
The numerical scheme described in Sections 3 and 4 has been implemented in the in-house parallel solver DGFlows . We employ METIS to partition the mesh [47] and the MPIbased software library PETSc [48] to assemble and solve all linear systems with iterative Krylov methods. The implicit convection treatment leads to non-symmetric matrices for the momentum and turbulence equations. For this reason, those linear systems are solved with the GMRES method combined with a block Jacobi preconditioner, with one block per process and an incomplete LU factorization on each block. On the other hand, the pressure-Poisson system is symmetric and positive definite. Hence, we solve it with the conjugate gradient method and an additive Schwarz preconditioner, with one block per process and an incomplete Cholesky decomposition on each block.
Solving the pressure Poisson equation is the most computationally expensive step of our algorithm. This is partially due to the fact that the pressure matrix in Eq. (42) corresponds to a Laplacian operator discretized with the local DG method, as explained by Shahbazi et al. [16] , and so it is characterized by larger stencil size than the SIP diffusive operator. However, note that the matrix is the same at each time step. Therefore, it is possible to assemble it and compute its preconditioner only once, thus partially alleviating the computational burden. On top of this, we initialize all Krylov solvers with the solution predictors described in Section 4 to speed convergence up.

Choice of the solution polynomial order
In this work, we used a mixed-order discretization for velocity and pressure (i.e., P p = P u − 1 ), which satisfies the inf-sup condition and guarantees optimal error convergence rates and stability in the low-t limit [33,49] .
The polynomial order of the turbulence quantities equals that of the pressure, which we believe is essential for a pressure-based DG solver. More precisely, the solution space of an arbitrary transported scalar quantity should be a subset of the solution space of the pressure. We have already mentioned this in [33] , but here we corroborate it with more rigorous theoretical and numerical arguments. This choice of polynomial order is in contrast with previous literature on mixed-order DG schemes. For example, Klein et al. [50] chose the same solution space for the temperature as for the components of the velocity field. We suspect that they found good results because their tests were done at a low Prandtl number of 0.7, whereas the problem with the solution spaces manifests itself when the convective term dominates, as shown in the following.
Our choice of polynomial order for any scalar solution stems from the fact that the continuity equation is weighed by the pressure basis functions, so that the numerical velocity satisfies the incompressibility constraint only in a weak sense up to order P p (i.e., D u h = 0 ). This means that the convective discretization can only be consistent up to an order P p .
To show this, consider the convection of a generic passive scalar quantity φ with the numerical velocity u h . For simplicity, assume that the domain is closed (i.e., n · u = 0 on ), so that substituting β ← u h and the continuous solution (i.e., w ← φ) into Eqs. (21) , (22) , and (23) gives the convective discretization for a test function v ∈ V h,φ . Here, we have used the fact that φ is continuous, so that φ = 0 and { φ} = φ, and so Compare this to the discrete continuity equation ( Eqs. (32) and (33) ), which, upon substituting v ← u h , integrating by parts, and using the fact that q u h = q { u h } + { q } u h on an interior face, can be rewritten as  which only vanishes if V h,φ is a subset of V h,p , that is, if the polynomial order of the pressure is at least as high as that of the scalar. If P φ > P p , the constant solution cannot be preserved in general due to numerical error introduced by the discontinuity across the elements of the velocity field. We corroborate this with a simple numerical example based on a standard lid-driven cavity problem. The domain is depicted in Fig. 1 . The Reynolds number is 40, based on the lid-velocity U and the characteristic length L , so the flow is laminar. We first found the steady-state velocity ( u h ) and pressure fields by solving the laminar Navier-Stokes equations. Then, we solved a transport equation for a passive scalar φ with the numerical velocity u h and a diffusion coefficient D : with homogeneous Neumann boundary conditions for φ. We set a constant initial condition φ 0 = φ(t = 0) = 1 , and took a single time step with the BDF1 scheme and U t/L = 1 , choosing a relative tolerance of 10 −12 for the GMRES and CG Krylov solvers, to minimize any numerical error in the solution of the linear systems. Clearly, the discrete solution φ h should remain unchanged after the time step. We therefore calculated the L 2 -norm error of φ h with respect to the expected exact solution φ 0 . To highlight the impact of the velocity divergence and continuity errors, we repeated the test for various spatial refinements. Following [52] , we used to estimate the mass conservation error in a strong sense. Table 4 collects the results for two structured meshes made of quadrilateral elements: (i) 50 by 50 uniform (indicated with M1 ); and (ii) 90 by 90 progressively finer towards the boundaries ( M2 ). On M2 , we studied the effect of varying the velocity polynomial order.
Results support our theoretical argument. When P φ > P p , the constant solution is not preserved, due to the numerical error introduced by discontinuity of the velocity field across the elements. In fact, as expected, the error decreases with spatial refinement, that is, with lower divergence and continuity errors. The problem is more pronounced in convection-dominated problems, as in case of turbulent flows. In fact, the biggest errors correspond to low values of the diffusion coefficient. Increasing D , as expected the error decreases, because the residual becomes progressively more dominated by the elliptic term. Asymptotically, the error is inversely proportional to D . On the other hand, the constant solution is always preserved (up to numerical precision) when P φ = P p .

Test cases
In this section, we verify and benchmark our numerical scheme with three test cases. First, we employ a manufactured solution for the RANS and K − E equations to verify the spatial and temporal convergence rates of the solver. Then, we simulate the stationary turbulent flow over a backward-facing step and finally a Von Kármán vortex street in the wake of a square cylinder. In both cases, we compare our results with those obtained from experiments and other numerical simulations reported in literature.
The meshes were generated with the open-source tool Gmsh [53] , which is also used to post-process the solution fields.

Manufactured solution
As first test-case, we employ the method of manufactured solutions to verify the correctness of the implementation of our spacetime numerical scheme. Starting from the well-know Taylor-Green vortex solution of the laminar Navier-Stokes equations (see, for example, [ 16 ]), we generalized it to the RANS equations and included an expression for the turbulence quantities. The analytical solution is This solution is defined on the domain (x, y ) ∈ [ −L, L ] 2 and t ≥ 0.
It does not satisfy the coupled RANS and K − E equations exactly, so the extra source terms f , q K and q E are non-zero. We computed Table 4 Lid-driven cavity test: L 2 -norm error in the scalar quantity φ with respect to the exact solution φ 0 = 1 , as a function of the diffusion coefficient ( D ) and for various meshes and polynomial orders, which affect the divergence ( E D ) and continuity ( E C ) errors of the velocity field. Note that φ 0 = 1 always lies in the solution space. The scalar quantity is only conserved when P φ = P p .  them by substituting (49) into Eqs. (1a-1b) and (4a-4b) symbolically with a Python script. We imposed Dirichlet conditions at all boundaries of the domain. Together with the initial conditions, they were evaluated from the analytical solution.
All relative errors reported in the following were evaluated for L = 1 , ν = 2 . 5 × 10 −4 , and at t = 4 , which corresponds to almost a threefold decay of the initial condition. We set a relative tolerance of 10 −12 for the GMRES and CG Krylov solvers. Moreover, all integrals were evaluated with a polynomial accuracy of 20, considerably higher than the usual 3 P u − 1 . We did this to minimize the quadrature error in the integration of the extra source terms, which are complicated non-polynomial functions.

Temporal convergence
To verify the temporal convergence of our numerical scheme, we solved Eqs. (1a-1b) and (4a-4b) on a structured mesh consisting of 100 by 100 quadrilateral elements adopting a polynomial order P u = 5 , in order to ensure that the time error was dominant. We carried out the study for both the first-order and second-order BDF schemes, progressively halving the time step size from 2 −6 to 2 −11 . Fig. 2 shows the relative errors in the L 2 -norm of the velocity, pressure, and turbulence quantities with respect to the analytical solution (49) as a function of the time step size. The dashed lines help verifying that the error convergence rate is the expected one for both BDF orders. A slight deviation from the second-order convergence trend can be noticed in the pressure and velocity plots at small t , due to the onset of the space-error saturation.
Given the availability of an exact expression for the turbulent viscosity, we could compute an L 2 error also for this quantity. As shown in Fig. 3 , this error also exhibits the correct convergence rates.

Spatial convergence
The spatial convergence was verified by solving Eqs. (1a-1b) and (4a-4b) on increasingly refined structured quadrilateral meshes, progressively halving the element size from 2 −1 to 2 −5 . The study was repeated for various polynomial orders To ensure the spatial error was dominant, we chose the BDF2 scheme and a time step size t = 2 . 5 × 10 −5 .    dashed lines). On the other hand, the velocity error exhibits a sub-optimal convergence rate, in particular for the second-order and third-order discretizations. This is due to the fact that the residual of the momentum equation is dominated by the error in the turbulent viscosity, which converges at the theoretical rate for the scalar quantities, as shown in Fig. 5 .
To corroborate this conclusion, we repeated the test, this time imposing the exact expression of the eddy viscosity (from Eq. (49) ) in the momentum equation, thus decoupling the RANS set from the turbulence equations. Results reported in Fig. 6 demonstrate that the theoretical convergence rate of the velocity error can be recovered in this case.

Flow over a backward-facing step
As second test case, we study the turbulent flow over a backward-facing step. The problem is a well-known benchmark for turbulence solvers and many variants have been investigated. Here, we consider the same configuration as in the experimental study by [54] 2 , depicted in Fig. 7 . The ratio of the inlet channel height ( H i ) to step height ( H ) is 2 and Re = 44737 , based on the reference free-stream velocity U 0 , the step height and the kinematic viscosity ν.
We imposed Dirichlet boundary conditions at the inlet, placed at a distance L i = 4 H upstream of the step. We used an analytical profile for the stream-wise velocity ( u x ) that matches the experimental data provided by Kim et al. [54] (reasonably) well, 2 Experimental measurements by Kim et al. [54] have been retrieved from the turbulent flow database described in [55] and available at https://turbmodels.larc. nasa.gov/bradshaw.html .   We computed the steady-state solution on three structured meshes, finer towards the wall boundaries, consisting of 12,420 (mesh indicated with M1 in the following), 29,700 ( M2 ), and 55,500 ( M3 ) quadrilaterals. A portion of the latter mesh is shown in Fig. 8 . To study the influence of the polynomial order, we performed calculations with both P u = 3 and P u = 4 . Fig. 9 shows the steady-state velocity, turbulent kinetic energy, and eddy viscosity fields obtained on the finest mesh for P u = 4 . They are in good qualitative agreement with those, for example, reported by Kuzmin et al. [37] . For a more quantitative assessment of our results, Table 5 compares the dimensionless recirculation length ( L R ) we obtained with the experimental measurement by Kim et al. [54] and other numerical predictions. Our results differ Table 5 Comparison of the dimensionless recirculation length ( L R ) obtained in the present work (for different meshes and polynomial orders) with experimental and numerical results (obtained with the k − model) reported in literature. The relative error is computed with respect to the experimental measurements by Kim et al. [54] .

Reference
L R / H Relative error Kim et al. [54] , experiment 7.0 ± 0.5 Zijlema [56] 5.90 −15.7% Thangam & Speziale [57] 6.25 −10.7% Muller et al. [58] 6.33 −9.6% Kuzmin et al. [37] 5.40 −22.9% Shih et al. [59] 6.35 −9.3% Ilinca et al. [35] 6.20 −11.4% Present, only a few percent at most from the most accurate predictions reported in literature, and always towards a better agreement with the experimental results. In all cases, the recirculation length is underestimated by around 6%, which is expected with the k − model. Differences in L R on M2 and M3 are limited to 0.5% at most, and the impact of a richer polynomial representation is limited to less than 1%, thus showing the mesh convergence of our calculations.
Finally, Fig. 10 shows the vertical profiles of the stream-wise velocity component at different locations downstream of the step, obtained with P u = 4 and for different meshes. Profiles on M2 and M3 are barely distinguishable. Results compare reasonably well with the experimental measurements of Kim et al. [54] up to x/H = 8 . 0 , then the discrepancies due to the limitations of the k − model become apparent, as already documented in [35] , for example.

Vortex-shedding in the wake of a square cylinder
The last simulation is of turbulent flow past a cylinder with square cross section, characterized by the appearance of a Von Kármán vortex street in the wake of the obstacle. Even if the flow features stochastic three-dimensional fluctuations, which can be resolved only with LES or DNS approaches (see, e.g., [ 60,61 ]), the problem is used as a two-dimensional test case for RANS solvers. The flow is in fact characterized by a mean periodicity that can be captured with unsteady RANS equations.   Symmetry conditions were imposed on the top and bottom boundaries (distant H = 30 D ), whereas an outlet was placed at L o = 40 D downstream of the obstacle to avoid any influence on the flow in the cylinder region. For the same reason, the inlet boundary was placed at L i = 20 D upstream of the cylinder. The following uniform Dirichlet condition for the velocity was imposed: Moreover, following [63] , we calculated the inlet values of K and E to match a turbulence intensity level of 2% and a viscosity ratio ν t /ν = 10 , which correspond to the experiments of Lyn et al. [62] .
Finally, wall functions were imposed on the edges of the cylinder, setting the distance of the computational domain from the physical wall to δ w = 0 . 015 .
We sampled the time variations in the drag and lift coefficients, defined as where F x and F y indicate the stream-wise and vertical components of the force acting on the cylinder. The mean periodicity

Table 6
Comparison of the force coefficients (in terms of average value, root mean square of the fluctuations, and frequency) obtained in the present work for three meshes and time step sizes with experimental and numerical results reported in literature.  Fig. 12 ) triangles, progressively more refined towards the obstacle. On the finest mesh then, we varied the time step size to t = 2 . 5 × 10 −3 and t = 5 . 0 × 10 −3 . Fig. 13 shows an example of the instantaneous velocity magnitude and eddy viscosity fields obtained for t = 1 . 25 × 10 −3 and P u = 3 on M3 . The vortex shedding is clearly visible. Table 6 reports the results we obtained in terms of Strouhal number, average C d , and root mean square values of the fluctuations of C d and C l , comparing them with experimental measurements, DNS and other RANS results. We ended our simulations at t end = 152 , corresponding to almost 2.5 flow-through periods which is sufficient to reach flow periodicity, and evaluated the quantities of interest over the final 8 oscillation periods.

Reference
Results on M3 and for t = 1 . 25 × 10 −3 are fully converged in time (they do not change even quadrupling the t ) and nearly mesh independent. The Strouhal number agrees with the experimental measurements and the DNS predictions, and it differs only 1.5% from the value reported by Bosch and Rodi [63] for a standard k − model. The discrepancy is much more relevant when comparing the other quantities of interest though. This might be due to the great sensitivity of the force fluctuations to the numerical details of the problem (e.g., position of inlet boundary, Dirichlet conditions imposed, blockage ratio) as documented, for example, by Bosch and Rodi [63] and Han et al. [66] . In any case, our values are inside the range of RANS results reported by Rodi et al. [64] and Voke [65] . The discrepancy with experimental and DNS results in terms of force coefficients is expected, due to the limitations of the k − model in predicting the strength of the shedding motion.

Conclusions
In this work, we have presented a novel high-order discontinuous Galerkin Finite Element (DG-FEM) solver for the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations coupled with the k − closure model. Positivity of the turbulence quantities is ensured by solving for their logarithm. To avoid expensive integration of the equations up to the wall boundary, which also leads to poor results in case of the standard k − model, we have investigated the possibility of employing a two-velocity scale wall function approach in combination with DG, thus bridging the gap between the viscous and the logarithmic layers.
Contrary to most of previous literature, our solver features a second-order accurate algebraic pressure correction scheme to solve the coupled RANS system. The approach does not rely then on an artificial compressibility parameter whose value must be tuned for each problem. Implicit backward differentiation formulae are adopted for time discretization of all the equations. Proper time extrapolation of the solution ensures global second-order accuracy even if no iterations are performed between the highly coupled RANS and k − equations. For space discretization, we have chosen the Symmetric Interior Penalty method to deal with elliptic terms and the Lax-Friedrichs flux for convection.
We have paid particular attention to the choice of polynomial order for the turbulence quantities. We have theoretically and numerically demonstrated that, to avoid any numerical instability, the polynomial order for the pressure and any transported scalar quantity must be the same when the velocity field is not point-wise divergence free. This contrasts with the choices made in previous literature.
We have verified the correct implementation of our space-time numerical scheme with a manufactured solution. However, it has shown that the convergence rate of the velocity error could be degraded by the slower convergence of the eddy-viscosity error (due to the lower polynomial order discretization of the turbulence quantities), in case the residual is dominated by the diffusive term. To prevent this, one could move to an equal-order discretization or, to avoid any instability problem, try to reduce the divergence and continuity errors of the velocity field with extra element-wise penalty terms in the momentum equation as proposed by Krank et al. [14] . Another possibility would be to apply a post-projection operator on the convective field to ensure that the incompressibility constraint is satisfied exactly [9,51] .
Simulations of turbulent flows over a backward-facing step and past a square cylinder, inducing vortex shedding, have shown the soundness of our approach. Though we are bound by the wellknown limitations of the k − model, our results are in general good agreement with those reported in literature. Compared to previous works that used the k − model, our discrepancies are always towards a better agreement with the experimental or DNS/LES results.

Disclaimer
The content of this paper does not reflect the official opinion of the European Union. Responsibility for the information and/or views expressed therein lies entirely with the authors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.