A projection hybrid high order finite volume/finite element method for incompressible turbulent flows

A projection hybrid high order finite volume/finite element method for incompressible turbulent flows S. Bustoa, J.L. Ferr´ına,b, E. F. Toroc, M.E. V´azquez-Cend´ona,b,∗ aDepartamento de Matem´atica Aplicada, Universidade de Santiago de Compostela, Facultad de Matem´aticas. ES-15782 Santiago de Compostela, Spain bITMATI, Campus Sur. ES-15782 Santiago de Compostela, Spain cLaboratory of Applied Mathematics, DICAM, University of Trento. IT-38100 Trento, Italy Abstract In this paper the projection hybrid FV/FE method presented in [1] is extended to account for species transport equations. Furthermore, turbulent regimes are also considered thanks to the k−εmodel. Regarding the transport diffusion stage new schemes of high order of accuracy are developed. The CVC Kolgan-type scheme and ADER methodology are extended to 3D. The latter is modified in order to profit from the dual mesh employed by the projection algorithm and the derivatives involved in the diffusion term are discretized using a Galerkin approach. The accuracy and stability analysis of the new method are carried out for the advection-diffusion-reaction equation. Within the projection stage the pressure correction is computed by a piecewise linear finite element method. Numerical results are presented, aimed at verifying the formal order of accuracy of the scheme and to assess the performance of the method on several realistic test problems. Keywords: incompressible flows, k−εmodel, species transport, projection method, finite volume method, LADER, finite element method. ∗Corresponding author Email addresses: [email protected] (S. Busto), [email protected] (J.L. Ferr´ın), [email protected] (E. F. Toro), [email protected] (M.E. V´azquez-Cend´on) Preprint submitted to Journal of Computational Physics February 26, 2024 2 1. Introduction1 Finite volume methods combined with approximate Riemann solvers have2 been successfully developed for different kinds of flows in the 1980’s (see, [2]3 and the references therein). Focusing on the incompressible case, pressure4 results in a Lagrange multiplier that adapts itself to ensure that the velocity5 satisfies the incompressibility condition. In order to handle this situation,6 the typical explicit stage of finite volume methods has to be complemented7 with the so-called projection stage where a pressure correction is computed8 in order to get a divergence-free velocity. Many papers exist in the literature9 devoted to introduce and analyse projection finite volume methods for in-10 compressible Navier-Stokes equations (see, for instance, [3] or [4]). In order11 to get stability, staggered grids have been used to discretize velocity and pres-12 sure. While this can be done straightforwardly in the context of structured13 meshes, the adaptation to unstructured meshes is more challenging (see [5],14 [6], [7], [8], [9], [10]).15 On the other hand, projection methods have been also used in combina-16 tion with finite element discretizations (see [11]). Within this approach, the17 divergence-free condition for the velocity is replaced by an equation prescrib-18 ing the divergence of the linear momentum density which is a conservative19 variable.20 The scope of this paper is to extend the hybrid FV/FE projection method21 introduced in [1] for both laminar and turbulent flows considering also trans-22 port of species. Furthermore, new methods to increase the accuracy of the23 methodology are developed.24 Starting from a 3D tetrahedral finite element mesh of the computational25 domain, the equation of the transport-diffusion stage is discretized by a fi-26 nite volume method associated with a dual finite volume mesh where the27 nodes of the volumes are the barycentre of the faces of the initial tetrahedra.28 These volumes, which allow us for an easy implementation of flux boundary29 conditions, have already been used, among others, for the 2D shallow water30 equation (see [5]), for solving conservative and non conservative systems in31 2D and 3D (see [12] and [13]) and for DG schemes employed to solve com-32 pressible Navier-Stokes equations (see [14]). For time discretization we use33 the explicit Euler scheme. The convective term is upwinded using the Ru-34 sanov scheme (see [2] and [15]). Concerning the projection stage, the pressure35 correction is computed by continuous piecewise linear finite elements asso-36 ciated with the initial tetrahedral mesh. The use of the above “staggered”37 3 meshes together with a simple specific way of passing the information from38 the transport-diffusion stage to the projection one and vice versa leads to39 a stable scheme. The former is done by redefining the conservative variable40 (i.e. the momentum density) constant per tetrahedron. Conversely, the finite41 element pressure correction is redefined to constant values on the faces of the42 finite volumes and then used in the transport-diffusion stage.43 The coupling of Navier-Stokes equations and the turbulence model intro-44 duces turbulent viscosity which is typically computed by solving an additional45 pair of advection-diffusion-reaction equations, that is equations for the tur-46 bulent kinetic energy and the dissipation rate. One issue here is the time47 dependency of the viscous terms. This requires the use of methods that are48 at least second-order accurate in space and time for all terms involved (see49 [16] and [17]).50 For advection equations, several approaches for constructing high-order51 methods have been put forward. A classical example is the Lax-Wendroff52 scheme (see [18] and [19]). This scheme is linear in the sense of Godunov53 and thus oscillatory, according to Godunov’s theorem, [20]. A major step54 forward in this direction was the work of Kolgan [21], who introduced, for55 the first time, a numerical scheme that circumvents Godunov’s theorem, via56 the construction of a non-linear scheme using non-linear reconstructions (see57 [22]). Following these works, a new Kolgan-type method has been presented58 for the shallow water equations in [23]. In what follows, we will refer to this59 scheme as the CVC Kolgan-type scheme.60 In the present paper, the CVC Kolgan-type scheme is analysed and im-61 plemented at the transport-diffusion stage for the convective terms of the62 considered conservation laws: momentum conservation, transport equations63 and k-εmodel. The obtained scheme, second order in space and first order64 in time is combined with a Galerkin approach of the gradients involved in65 the viscous term. An alternative option will be the decomposition of the66 diffusion term into its orthogonal and non-orthogonal parts as introduced in67 [1].68 More advanced non-linear methods for advection dominated problems69 have appeared in the literature since the introduction of Kolgan scheme.70 Some of them are: Total Variation Diminishing Methods (TVD), Flux Lim-71 iter Methods, MUSCL-Hancock, semi-discrete ENO or WENO (see, for in-72 stance, [24], [25], [26], [24], [27], [28], [29] and [30]). Comprehensive reviews73 are found in [2] and [31], for example. Focusing on high order in time and74 space methodologies, we highlight the ADER approach, first put forward in75 4 [32]. It is also a fully discrete approach that relies on non-linear reconstruc-76 tions and the solution of the generalised Riemann problem, to any order of77 accuracy. The resulting schemes are arbitrarily accurate in both space and78 time in the sense that they have no theoretical accuracy barrier. An intro-79 duction to ADER schemes is found in Chapters 19 and 20 of [2]. Further80 developments and applications are found, for example, in [33], [34], [35], [36],81 [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49].82 In [50] an extension of ADER methodology to solve the advection-diffusion-83 reaction equation, admitting space and time dependent diffusion coefficients84 was introduced. The present work includes a modification of this scheme, the85 Local ADER method (LADER), which profits from the dual mesh. More-86 over, an ENO-based reconstruction is considered in order to prevent spurious87 oscillations.88 To asses the performance of the methodology different manufactured so-89 lution tests are introduced and the numerical results obtained with the devel-90 oped computer code are shown. Furthermore, several classical test problems91 from fluid mechanics are presented and some results are compared with ex-92 perimental data (see [17] and [1]).93 The paper is organized as follows. In Section 2 the mathematical model94 for incompressible flows is recalled. Then, the RANS k−εmodel for the95 turbulence and the species transport equations are described. In Section96 3 the numerical discretization is detailed. Special attention is paid to the97 description of the finite volume algorithm. Aiming to achieve a high order98 scheme, two different methodologies for the flux terms are developed: the99 CVC Kolgan-type method, second order in space and first order in time,100 and the LADER methodology, second order in space and time. The needed101 modifications on the approximation of remaining terms of the equations to102 achieve a high order scheme are also presented and a Galerkin approach103 to compute the diffusion terms is introduced. Finally, in Section 4 some104 of the numerical results obtained with the developed code are shown. On105 the one hand, the order of convergence of the method is analysed using the106 method of manufactured solutions. On the other hand, several classical test107 problems are analysed. The appendix includes the theoretical analysis of108 the LADER numerical method applied to the one-dimensional advection-109 diffusion-reaction equation.110 5 2. Governing equations111 In this section, the system of equations to be solved is introduced. The112 model for incompressible newtonian fluids, recalled in [1], is extended con-113 sidering a turbulent regime and taking into account the transport of species.114 2.1. Mass conservation and momentum equations115 The incompressible Navier-Stokes equations reduce to the mass conserva-116 tion equation and the momentum equation. Hence, the system of equations117 written in conservative variables reads118 divwu= 0,(1) ∂wu ∂t + divFwu(wu) + ∇π−div(τ) = fu.(2) Standard notation is used:119  ρis the density (kg/m3),120  p=π+πis the pressure (N/m2),121 –πis the mean pressure,122 –πis the pressure perturbation,123  u= (u1, u2, u3)tis the velocity vector (m/s),124  wu:= ρuis the vector of the conservative variables related to the125 velocity (kg/s m2),126  Fwuis the flux tensor:127 Fwu i(wu) = 1 ρwiwu=uiwu, i = 1,2,3,  τis the viscous part of the Cauchy stress tensor,128  fuis a generic source term used for the manufactured test problems.129 6 2.2. Turbulence model130 Special care for the viscous part of Cauchy stress tensor, τ, is required for131 turbulent regimes. In order to avoid the high computational cost of a direct132 simulation of the turbulence, the k−εstandard model is used (see [51] and133 [52]). The Reynolds-averaged viscous stress tensor is given by134 τ=τu+τR.(3) Denoting µthe laminar viscosity (kg/(m s)), the averaged stress tensor, τu,135 reads136 τu=µ∇u+∇uT.(4) The fluctuation, τR, called the Reynolds tensor, is given by137 τR=µt∇u+∇uT−2 3ρkI.(5) To obtain the turbulent viscosity,138 µt=ρCµ k2 ε,(6) two new variables are introduced: the turbulent kinetic energy, k(J/kg), and139 the energy dissipation rate, ε(J/(kg s)). They are computed from a new pair140 of partial differential equations, namely,141 ∂wk ∂t + divFwk(wk,u)−div µ+µt σk∇wk ρ+wε=Gk+fk,(7) ∂wε ∂t + divFwε(wε,u)−div µ+µt σε∇wε ρ+C2ε w2 ε wk =C1ε wε wk Gk+fε,(8) where142  wk(J), wε(J/s) are the conservative variables corresponding to kand143 ε, that is144 wk:= ρk, wε:= ρε,  Fwk,Fwεare the fluxes related to the turbulence variables,145 Fwk i(wk,u) = uiwk,Fwε i(wε,u) = uiwε, 7  Gkis the term of kinetic energy production, due to the mean velocity146 gradients, of the Reynolds stress tensor,147 Gk=µt 2"3 X i=1 3 X j=1 ∂ui ∂xj +∂uj ∂xi#2 ,(9)  fk,fεare the source terms related to manufactured solutions for ana-148 lytical tests; they have zero value in physical problems,149  σk= 1.0, σε= 1.3 are the turbulent Prandtl numbers,150  Cµ= 0.09, C1ε= 1.44, C2ε= 1.92 are the closure coefficients of the151 model whose values were taken from the literature.152 2.3. Species transport153 The equations of transport of species are also included in the system to154 be solved:155 ∂wy ∂t + divFwy(wy,u)−div ρD+µt Sct∇1 ρwy=fy,(10) with156  y= (y1,...,yNe)Tthe mass fraction vector of the species to be consid-157 ered. yicorresponds to species iand Neis the total number of species158 to be considered,159  wy:= ρythe conservative variable vector related to the mass fraction160 vector,161  Fwythe flux,162 Fwy i(wy,u) = uiwy  Dthe mass diffusivity coefficient (m2/s),163  Sct= 0.7 the turbulent Schmidt number,164  fythe source term for manufactured test problems.165 8 2.4. Complete system166 The complete system of equations to be solved is167 divwu= 0,(11) ∂wu ∂t + divFwu(wu) + ∇π−div(τ) = fu,(12) ∂wk ∂t + divFwk(wk,u)−div µ+µt σk∇wk ρ+wε=Gk+fk,(13) ∂wε ∂t +divFwε(wε,u)−divµ+µt σε∇wε ρ+C2ε w2 ε wk =C1ε wε wk Gk+fε,(14) ∂wy ∂t + divFwy(wy,u)−div ρD+µt Sct∇1 ρwy=fy.(15) Moreover, the vector of the conservative variables is w= (wu,b w)T,with b w168 the vector of the conservative variables related with turbulence, and species,169 i.e., b w= (wk, wε,wy)T. The flux tensor of the complete system has three170 components:171 F= (F1|F2|F3)(3+2+Ne)×3,Fi(w) = wi ρw, i = 1,2,3.(16) 3. Numerical discretization172 The numerical discretization of the complete system is performed by ex-173 tending the projection method first put forward in [1]. The developed nu-174 merical method solves, at each time step, equations (12)-(15) with a finite175 volume method (FVM) and, so, an approximation of wis obtained. The176 next step is applying projection to system (11)-(12). The pressure correction177 is provided by a piecewise linear finite element method (FEM). In the post-178 projection step, an approximation of wuverifying the divergence condition,179 (11), is obtained. Furthermore, the production terms of the turbulence equa-180 tions are also computed in this step to account for the corrected velocities.181 The reaction terms are treated via a semi-implicit method.182 We start by considering a two-stage in time discretization scheme: in or-183 der to get the solution at time tn+1, we use the previously obtained approx-184 imations Wnof the conservative variables w(x, y, z, tn), Unof the velocity185 u(x, y, z, tn) and πnof the pressure perturbation π(x, y, z, tn), and compute186 Wn+1 and πn+1 from the following system of equations:187 1 ∆tf Wn+1 u−Wn u+ divFwu(Wn u) + ∇πn−div(τn) = fn u,(17) 9 1 ∆tWn+1 u−f Wn+1 u+∇(πn+1 −πn) = 0,(18) divWn+1 u= 0,(19) 1 ∆tf Wn+1 k−Wn k+ divFwk(Wn k,Un)−div µ+µt σk∇Wn k ρ= 0,(20) 1 ∆tWn+1 k−f Wn+1 k+Wn ε−Gn+1 k=fn k,(21) 1 ∆tf Wn+1 ε−Wn ε+ divFwε(Wn ε,Un)−div µ+µn t σε∇Wn ε ρ= 0,(22) 1 ∆tWn+1 ε−f Wn+1 ε+C2ε Wn+1 εWn ε Wn k −C1ε Wn ε Wn k Gn+1 k=fn ε,(23) 1 ∆tf Wn+1 y−Wn y+divFwyWn y,Un−div ρD+µn t Sct∇1 ρWn y = 0,(24) 1 ∆tWn+1 y−f Wn+1 y=fn y.(25) Concerning the discretization of mass conservation and momentum equa-188 tions, by adding equations (17)-(18), we easily see that the scheme is actually189 implicit for the pressure term. However, the writing above shows that the190 pressure and the velocity can be solved in three uncoupled stages. The first of191 them corresponds to equation (17) and will be called the transport-diffusion192 stage; it is explicit and allows us to compute the intermediate approximation193 of the conservative variables f Wn+1 u(we notice that, in general, this approx-194 imation does not satisfy the divergence condition (19)). The second one, to195 be called the projection stage, is implicit; it consists of solving the coupled196 equations (18) and (19) with a finite element method to obtain the pressure197 correction δn+1 := πn+1 −πn. The last one is the post-projection stage; the198 intermediate approximation for the velocity conservative variables is updated199 with the pressure correction providing the final approximations Wn+1 uand200 πn+1 (see [1] for further details).201 As a novelty in this paper, for the remaining conservation laws the ap-202 proximation of the conservative variables is obtained in two steps. At the203 transport-diffusion stage we compute an approximation of the conservative204 variables, f Wn+1 k,f Wn+1 εand f Wn+1 y, taking into account the corresponding205 flux and diffusion terms. Let us remark that at this stage the update of the206 approximations involves all the neighbouring nodes of the finite volume Ci.207 16 Step 2. Computation of boundary extrapolated values at the barycenter of328 the faces, Nij:329 Wi Nij =pi ij(Nij) = Wi+ (Nij −Ni) (∇W)i ij ,(41) Wj Nij =pj ij(Nij)Wj+ (Nij −Nj) (∇W)j ij .(42) Step 3. Computation of the flux terms with second order of accuracy us-330 ing the mid-point rule. Taylor series expansion in time and Cauchy -331 Kovalevskaya procedure are applied to locally approximate the conser-332 vative variables at time ∆t 2. This methodology accounts for the con-333 tribution of the advection and diffusion terms to the time evolution of334 the flux term. The resulting evolved variables read335 Wi Nij =Wi Nij −∆t 2Lij Z(Wi Nij ,ηij) + Z(Wj Nij ,ηij) +∆t 2L2 ij αiNij (∇W)i ij ηij +αjNij (∇W)j ij ηij,(43) Wj Nij =Wj Nij −∆t 2Lij Z(Wi Nij ,ηij) + Z(Wj Nij ,ηij) +∆t 2L2 ij αiNij (∇W)i ij ηij +αjNij (∇W)j ij ηij.(44) We have denoted Lij = min n|Ci| S(Ci),|Cj| S(Cj)owith S(Ci) the area of the336 surface of cell Ci. Two different options will be considered in the scheme337 concerning the evolved variables. The first one corresponds to the338 previous definition of the evolved variables. Meanwhile, the second339 one neglects the evolution of the diffusion term.340 Step 4. Computation of the numerical flux considering Rusanov scheme.341 ϕWn i Nij , Wn j Nij ,ηij=1 2ZWn i Nij ,ηij+ZWn j Nij ,ηij −1 2αRS Wn i Nij ,Wn j Nij ,ηijWn j Nij −Wn i Nij .(45) 3.4. Viscous terms342 We come next to describe the computation of the integrals involving the343 viscous terms. First, applying Gauss’ theorem we relate the volume integral344 17 of the diffusion term with a surface integral over the boundary, Γi. Next,345 this integral is split into the integrals over the cell interfaces Γij. Thus, the346 viscous term of the momentum conservation equation reads347 ZCi divτndV = X Nj∈KiZΓij τn˜ ηij dS =X Nj∈KiZΓij (µ+µn t)∇Un+ (∇Un)T−2 3ρKnIe ηij dS.(46) Two different approaches can be considered in order to compute the above348 integral.349 On the one hand, decomposition with semi-implicit and explicit dis-350 cretizations can be applied when using the CVC Kolgan-type scheme. This351 methodology splits the diffusion flux into its orthogonal and non-orthogonal352 parts and relax the stability condition on the time step size (see [1] for further353 details).354 On the other hand, the dual mesh ease the use of Galerkin approach to355 compute the derivatives involved in (46). We introduce a numerical diffusion356 function φusuch that357 ZΓij (µ+µn t)∇Un˜ ηijdS ≈φuUn i,Un j, µn t, i, µn t, j,ηij(47) and we consider358 φuUn i,Un j, µn t, i, µn t, j,ηij=µ+µn t, ij(∇Un)Tij ηij,(48) with359 µn t, ij =1 2µn t, i +µn t, j.(49) Since we know the value of the turbulent kinetic energy at the nodes of360 the finite volumes, we approximate the turbulent kinetic energy term as the361 average of the values obtained at the two nodes related to the face362 −ZΓij 2 3Wn k˜ ηijdS = −1 3Wn k, i +Wn k, jηij.(50) 18 Finally, the viscous terms for the remaining equations are obtained equally363 to the gradient term of the momentum equation:364 ZCi 1 ρdiv Dn∇c WndV = 1 ρX Nj∈KiZΓij Dn∇c Wn˜ ηij dS,(51) where365 Dn= Dn k0 0 0Dn ε0 0 0 Dn y  =       µ+µn t σk 0 0 0µ+µn t σε 0 0 0 ρD+µn t Sct        .(52) Thus, we can introduce the diffusion flux function, φb w, verifying366 ZΓij Dn∇c Wn˜ ηijdS ≈φb wc Wn i,c Wn j, µn t, i, µn t, j,ηij, φb wc Wn i,c Wn j, µn t, i, µn t, j,ηij=Dn ij ∇d WnTij ηij (53) The above methodologies are used to approximate the viscous terms when367 choosing a first order method or the CVC Kolgan-type scheme to compute the368 advection term. Nevertheless, the LADER methodology requires a special369 treatment.370 3.4.1. LADER approach: the viscous terms371 As was already introduced in Section 3.3.2, to apply LADER and to372 obtain a second order in space and time scheme, instead of computing the373 diffusion flux functions, φuand φc W, with the value of the variables at the374 previous time step, Un,Kn,Enand Yn, its is necessary to use some evolved375 values, Un,Kn,Enand Yn.376 It is important to remark that the former evolved variables do not match377 the already computed ones for the flux term (see Appendix A for a detailed378 analysis of the scalar advection-diffusion-reaction equation). Taylor series379 expansion in time and Cauchy-Kovalevskaya procedure are applied neglect-380 ing the advection term so that a second order in space and time scheme is381 attained:382 Un=Un+∆t 2div (µ+µn t)∇Un−2 3ρKnI,(54) 19 Kn=Kn+∆t 2µ+µn t σk∇Kn,(55) En=En+∆t 2µ+µn t σε∇En,(56) Yn=Yn+∆t 2ρD+µn t Sct∇Yn.(57) In what follows, we describe the computation of the evolved velocities at383 an arbitrary node Ni:384 1. The gradients of the original variables are computed at each auxiliary385 tetrahedra of the FE mesh, Tij (see, on the 2D representation in Figure386 2, the triangle with green contour). The value of the gradient at each387 node, Ni, is obtained as the average of the values on the two tetrahedra388 containing the node, TijL (green filled triangle in Figure 2) and Tij.389 Taking into account the viscosity coefficients and the turbulent kinetic390 energy term, we introduce the auxiliary variable:391 Vn i:= µ+µn t, i1 2(∇Un)TijL + (∇Un)Tij −2 3ρKn iI. (58) 2. The divergence is computed as the average of the divergences of Vn 392 obtained on the auxiliary tetrahedra:393 Un i=Un i+∆t 4tr (∇Vn)TijL + (∇Vn)Tij .(59) 3. The diffusion function φuis evaluated on the evolved variables:394 φuUn i,Un j, µn t, i, µn t, j,ηij=µ+µn t, ij∇UnTij ηij.(60) The remaining evolved variables are similarly obtained. Hence the related395 diffusion function reads396 φb wc Wn i,c Wn j, µn t, i, µn t, j,ηij=Dn ij ∇c WnTij ηij.(61) 3.5. Pressure term397 For the integral of the pressure gradient we follow [1]. We split the bound-398 ary Γiinto the cell interfaces Γij using Gauss’ theorem and we compute the399 pressure as the arithmetic mean of its values at the three vertices of face Γij 400 and the barycentre of the tetrahedra to which the face belongs. Then, the401 corresponding approximation of the integral is given by402 ZΓij πne ηijdS ≈5 12 (πn(V1) + πn(V2)) + 1 12 (πn(V3) + πn(V4))ηij.(62) 20 3.6. Projection stage403 Within the projection stage, the pressure is computed using a standard404 finite element method. The incremental projection method presented in [11]405 is adapted to solve (18)-(19) obtaining the following weak problem:406 Find δn+1 ∈V0:= z∈H1(Ω) : RΩz= 0verifying407 ZΩ ∇δn+1 · ∇zdV = 1 ∆tZΩf Wn+1 · ∇zdV −1 ∆tZ∂Ω Gn+1zdS ∀z∈V0, (63) where δn+1 := πn+1 −πn(see [1] for further details).408 3.7. Post-projection stage409 Once the pressure is computed, we can update Wn+1 uwith ∇δn+1 i, that410 is,411 Wn+1 u, i =f Wn+1 u, i + ∆t∇δn+1 i.(64) The previous computation of the updated velocities allows for an implicit412 approach of the production term Gkon the turbulence equations. Meanwhile,413 for the dissipative terms a semi-implicit scheme is used:414 Wn+1 k, i −f Wn+1 k, i ∆t+Wn ε, i −Gk, i(Un+1) = fn k, i,(65) Wn+1 ε, i −f Wn+1 ε, i ∆t+C2ε Wn ε, i Wn k, i Wn+1 ε, i −C1ε Wn ε, i Wn k, i Gk, i(Un+1) = fn ε, i (66) where the derivatives involved in the production term, Gk, i(Un+1), are com-415 puted as the averaged of the auxiliary tetrahedra related to the node Ni.416 Finally, the source terms fb ware pointwise evaluated.417 3.8. Boundary conditions418 The boundary conditions were defined following [1]:419  Dirichlet boundary conditions for inviscid fluids: the normal component420 of the conservative variable is set at the boundary nodes.421  Dirichlet boundary conditions for viscous fluids: the value of the con-422 servative variable is imposed at the boundary nodes.423 21  Neumann boundary conditions: the definition of f Wn+1 takes into ac-424 count the inflow/outflow boundary condition with no need for any ad-425 ditional treatment.426 Moreover, in the manufactured tests designed to analyse the order of ac-427 curacy of the numerical discretizations, it is a usual practice to impose the428 values of the exact solution at the boundary nodes. This practice avoids that429 the accuracy of the method can be affected by the treatment of the bound-430 ary conditions. From the mathematical point of view, it is like considering431 Dirichlet boundary conditions.432 4. Numerical results433 In this section, we present the results obtained for several test problems.434 In order to define the time step, two different options are implemented in the435 code. On the one hand, we can simply introduce a fixed time step. On the436 other hand, we can provide the CFL from which the code will compute the437 time step at each time iteration. The latest option is the one chosen to run438 the test cases presented in this paper. Therefore, to determine the time step439 at each time iteration, we compute a local value for the time step at each440 cell Ci,441 ∆tCi=CFL L2 i 2|Ui| Li+ max µ+µt, i, ρD+µt, i Sct(67) with Li:= |Ci| S(Ci). Finally, as global time step at each time iteration, ∆t, we442 choose the minimum time steps obtained at each cell.443 Remark 1. The above definition of ∆tCiis valid if the transport of species444 equation is solved, otherwise its value is given by445 ∆tCi=CFL L2 i 2|Ui| Li+µ+µt, i .(68) 4.1. Manufactured test 1. Laminar flow446 The first test to be posed was obtained using the method of the manufac-447 tured solutions (MMS). We consider the domain Ω = [0,1]3and we assume448 the flow being defined by449 ρ= 1,(69) 22 π(x, y, z, t) = cos(πt(x+y+z)),(70) u(x, y, z, t) = sin(πyt) cos(πzt),−cos(πz3t),exp(−2πxt2)T,(71) with µ= 10−2. The related source terms are included in Appendix B.450 To perform the error and order of accuracy analysis we employ the three451 uniform meshes with different cell sizes presented in Table 1. We have de- Mesh NElements Vertices Nodes Vm h(m3)VM h(m3) M14 384 125 864 6.51E−04 1.30E−03 M28 3072 729 6528 8.14E−05 1.63E−04 M316 24576 4913 50688 1.02E−05 2.03E−05 Table 1: Manufactured test 1. Laminar flow. Mesh features. 452 noted N+1 the number of points along the edges, h= 1/N,Vm hthe minimum453 volume of the finite volumes and VM hthe maximum volume of the finite vol-454 umes.455 Four different methods are used to solve the problem: the first order456 method presented in [1], CVC method with an orthogonal decomposition457 of the diffusion term (CVC-orth), CVC method combined with a Galerkin458 approach for the diffusion term (CVC-G) and LADER. The errors and orders,459 depicted in Table 2, were computed as follows:460 E(π)Mi=∥π−πMi∥l2(L2(Ω)) E(wu)Mi=∥wu−wuMi∥l2(L2(Ω)3),(72) 461 oπMi/Mj=log E(π)Mi/E(π)Mj log hMi/hMj, owuMi/Mj=log E(wu)Mi/E(wu)Mj log hMi/hMj. (73) We can observe that CVC-G method provides an order of convergence close462 to two. This is in accordance with the theoretical order of this scheme, first463 order in time and second order in space, and the high time-dependency of464 the solution. Whereas, with LADER the expected second order of accuracy465 is achieved.466 23 Method Variable EM1EM2EM3oM1/M2oM2/M3 Order 1 π1.24E−01 5.70E−02 2.96E−02 1.12 0.95 wu6.40E−02 3.32E−02 1.78E−02 0.95 0.90 CVC-orth. π6.30E−02 1.91E−02 8.84E−03 1.72 1.11 wu5.51E−02 2.06E−02 8.95E−03 1.42 1.20 CVC-G π5.98E−02 1.58E−02 4.58E−03 1.92 1.78 wu5.41E−02 1.88E−02 6.52E−03 1.52 1.53 LADER π4.10E−02 8.74E−03 2.03E−03 2.23 2.11 wu2.61E−02 5.76E−03 1.24E−03 2.18 2.22 Table 2: Manufactured test1. Laminar flow. Observed errors and convergence rates. CFL = 1. 4.2. Manufactured test 2. Turbulent flow with species transport467 The second academic test to be considered is a modification of Test 1 to468 account for the turbulence and species transport equations. Let us define the469 flow as470 ρ= 1,(74) π(x, y, z, t) = cos(πt(x+y+z)),(75) u(x, y, z, t) = sin(πyt) cos(πzt),−cos(πz3t),exp(−2πxt2)T,(76) k(x, y, z, t) = sin(πxt)+2,(77) ε(x, y, z, t) = exp(−πzt)+1,(78) y(x, y, z, t) = sin(πxt)+2.(79) with parameters µ= 10−2,D= 10−3. For the exact solution to verify the471 equations, taught expressions of the source terms have to be taken. They472 have been included in Appendix B.473 We consider the meshes already defined in Table 1 and a CFL = 10474 (the reason why this large value of CFL is admitted was studied in [50]).475 Dirichlet boundary conditions are set for all the equations on the boundary.476 The computed errors are presented in Table 3. The results obtained for CVC-477 orth confirm that using only second order in space for computing the flux478 terms and neglecting the non orthogonal component will not capture properly479 the turbulence. Second order in space must also be used to approximate the480 diffusion terms and the whole flux should be computed. Furthermore, a481 second order in time scheme improves the results and order attained.482 24 Method Variable EM1EM2EM3oM1/M2oM2/M3 Order 1 π6.97E−01 5.52E−01 4.93E−01 0.34 0.16 wu4.40E−02 2.80E−02 2.18E−02 0.65 0.36 wk3.85E−02 2.45E−02 2.09E−02 0.65 0.23 wε1.53E−02 8.47E−03 6.40E−03 0.85 0.40 wy2.93E−02 2.04E−02 1.76E−02 0.52 0.21 CVC-orth. π6.23E−01 5.13E−01 4.75E−01 0.28 0.11 wu4.08E−02 2.65E−02 2.10E−02 0.62 0.33 wk3.18E−02 2.13E−02 1.95E−02 0.58 0.13 wε1.51E−02 8.17E−03 6.15E−03 0.89 0.41 wy2.51E−02 1.85E−02 1.68E−02 0.44 0.14 CVC-G π2.70E−01 7.60E−02 2.09E−02 1.83 1.86 wu1.50E−02 5.18E−03 1.49E−03 1.54 1.80 wk1.54E−02 3.24E−03 8.22E−04 2.25 1.98 wε1.06E−02 2.39E−03 6.35E−04 2.15 1.91 wy7.27E−03 1.89E−03 4.86E−04 1.94 1.96 LADER π2.68E−01 7.61E−02 2.10E−02 1.82 1.86 wu1.51E−02 5.17E−03 1.50E−03 1.55 1.79 wk1.37E−02 2.51E−03 5.89E−04 2.45 2.09 wε9.87E−03 1.80E−03 4.09E−04 2.46 2.14 wy7.25E−03 1.60E−03 3.79E−04 2.18 2.08 Table 3: Manufactured test 2. Turbulent flow. Observed errors and convergence rates. CFL = 10. 25 4.3. Test 3. Gaussian sphere483 The next problem to be analysed is the Gaussian sphere test introduced484 in [55] and [17]. We consider a normal distribution function in the domain485 Ω = [−0.9,0.9] ×[−0.9,0.9] ×[−0.3,0.3] with standard deviation 0.08 and486 mean 0.25. The density is one, the velocity vector is defined as u(x, y, z, t) =487 (−y, x, 0)Tand we assume that the diffusion matrix is given by D=µ.488 Hence, the solution of the problem is given by489 y(x, y, z, t) = σ0 σ(t)3 exp −r 2σ(t)2(80) with490 r(x, y, z, t) = (¯x+ 0.25)2+ ¯y2+z2, σ(t) = qσ2 0+ 2tD,(81) ¯x=xcos(t) + ysin(t),¯y=−xsin(t) + ycos(t).(82) The flow definition is completed setting the source terms491 fu(x, y, z, t) = (−x, −y, 0)T, fy(x, y, z, t)=0,(83) and considering Dirichlet boundary conditions.492 In order to analyse the accuracy in time and space, five structured meshes493 were generated. The properties of these meshes can be seen in Table 4, where494 hdenotes the size of the cubes used to generate the tetrahedra of the finite495 element mesh.496 Mesh Finite elements Vertices Nodes h M111664 2527 24408 0.1 M218522 3872 38514 0.0857 M354000 10571 111000 0.06 M493312 17797 190944 0.05 M5182250 33856 256711 0.04 Table 4: Test 3. Gaussian sphere. Mesh features. Table 5 shows the results obtained for the test considering µ= 10−3. On497 the other hand, in Table 6 the errors and orders of accuracy for µ= 10−2 498 are presented. In both test cases we have assumed a final time tend = 2π499 32 Figure 11: Test 4. Flow around a cylinder. Velocity magnitude on z= 0.205. Figure 12: Test 4. Flow around a cylinder. Pressure over the cylinder and streamlines. 5. Summary and conclusions538 In this paper a projection hybrid high order FV/FE method for incom-539 pressible flows has been presented. Navier-Stokes equations have been cou-540 pled with the k−εmodel in order to simulate turbulent flows. The system541 to be solved was enlarged with respect to [1] considering species transport.542 High order of accuracy has shown necessary for the proper computation of543 turbulent effects. Two different methodologies to achieve second order were544 presented. Firstly, CVC Kolgan provided a second order in space and first or-545 der in time scheme. To attain second order in space and time a new method546 was proposed: LADER. The corresponding accuracy and stability analy-547 sis were presented for the advection-diffusion-reaction equation. Godunov’s548 theorem was circumvented thanks to an ENO-based approach. The compu-549 tation of the gradients involved on diffusion terms was done via Galerkin.550 33 The method was applied to manufactured test problems in order to asses the551 accuracy. Furthermore, different benchmarks were considered and the results552 obtained were successfully confronted with experimental data.553 Acknowledgements554 This work was financially supported by Spanish MICINN projects MTM2008-555 02483, CGL2011-28499-C03-01 and MTM2013-43745-R; by the Spanish MECD556 under grant FPU13/00279; by the Xunta de Galicia Conseller´ıa de Cultura557 Educaci´on e Ordenaci´on Universitaria under grant Axudas de apoio ´a etapa558 predoutoral do Plan I2C; by Xunta de Galicia and FEDER under research559 project GRC2013-014 and by Fundaci´on Barri´e under grant Becas de pos-560 grado en el extranjero.561 References562 References563 [1] A. Berm´udez, J. L. Ferr´ın, L. Saavedra, M. E. V´azquez-Cend´on, A564 projection hybrid finite volume/element method for low-Mach number565 flows, J. Comp. Phys. 271 (2014) 360–378. doi:10.1016/j.jcp.2013.09.029.566 URL http://www.sciencedirect.com/science/article/pii/S0021999113006438567 [2] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics:568 A practical introduction, Springer, 2009. doi:10.1007/b79761.569 [3] J. B. Bell, A. S. Day, A. S. Almgren, M. J. Lijewski, C. A. Rendleman,570 A parallel adaptive projection method for low Mach number flows, Int.571 J. Numer. Methods Fluids 40 (2002) 209–216. doi:10.1002/fld.310.572 URL http://onlinelibrary.wiley.com/doi/10.1002/fld.310/full573 [4] R. B. Pember, A. S. Almgren, J. B. Bell, P. Colella, M. Howell, M. Lai,574 A high order projection method for the simulation of unsteady turbulent575 non premixed combustion in an industrial burner, in: Proceedings of the576 8th International Symposium on Transport Phenomena in Combustion,577 1995, pp. 16–20.578 [5] A. Berm´udez, A. Dervieux, J. A. Desideri, M. E. V´azquez-Cend´on, Up-579 wind schemes for the two-dimensional shallow water equations with vari-580 able depth using unstructured meshes, Comput. Methods Appl. Mech.581 Eng. 155 (1) (1998) 49–72. doi:10.1016/S0045-7825(97)85625-3.582 34 [6] T. Gallou¨et, L. Gastaldo, R. Herbin, J.-C. Latch´e, An uncondition-583 ally stable pressure correction scheme for the compressible barotropic584 Navier-Stokes equations, ESAIM: Mathematical Modelling and Numer-585 ical Analysis 42 (2) (2008) 303–331. doi:10.1051/m2an:2008005.586 URL https://doi.org/10.1051/m2an:2008005587 [7] W. Gao, H. Li, R. Liu, An unstructured finite volume projection588 method for pulsatile flows through an asymmetric stenosis, J. Eng.589 Math. 72 (1) (2012) 125–140. doi:10.1007/s10665-011-9469-1.590 URL https://link.springer.com/article/10.1007/s10665-011-9469-1591 [8] S. Perron, S. Boivin, J.-M. H´erard, A finite volume method592 to solve the 3D Navier-Stokes equations on unstructured col-593 located meshes, Computers & fluids 33 (10) (2004) 1305–1333.594 doi:https://doi.org/10.1016/j.compfluid.2003.10.006.595 URL http://www.sciencedirect.com/science/article/pii/S004579300400057X596 [9] S. Tu, S. Aliabadi, Development of a hybrid finite volume/element597 solver for incompressible flows, Int. J. Numer. Methods Fluids 55 (2)598 (2007) 177–203. doi:10.1002/fld.1454.599 URL http://onlinelibrary.wiley.com/doi/10.1002/fld.1454/full600 [10] S. Tu, S. Aliabadi, R. Patel, M. Watts, An implementation of the601 Spalart-Allmaras DES model in an implicit unstructured hybrid finite602 volume/element solver for incompressible turbulent flow, Int. J. Numer.603 Methods Fluids 59 (9) (2009) 1051–1062. doi:10.1002/fld.1864.604 URL http://onlinelibrary.wiley.com/doi/10.1002/fld.1864/abstract605 [11] J. L. Guermond, P. Minev, J. Shen, An overview of projection methods606 for incompressible flows, Comput. Methods Appl. Mech. Eng. 195 (2006)607 6011–6045. doi:10.1016/j.cma.2005.10.010.608 [12] E. F. Toro, A. Hidalgo, M. Dumbser, FORCE schemes on unstructured609 meshes I: Conservative hyperbolic systems, J. Comput. Phys. 228 (9)610 (2009) 3368 – 3389. doi:https://doi.org/10.1016/j.jcp.2009.01.025.611 URL http://www.sciencedirect.com/science/article/pii/S0021999109000412612 [13] M. Dumbser, A. Hidalgo, M. Castro, C. Par´es, E. F. Toro, FORCE613 schemes on unstructured meshes II: Non-conservative hyperbolic614 systems, Comput. Methods Appl. Mech. Eng. 199 (2010) 625–647.615 35 doi:https://doi.org/10.1016/j.cma.2009.10.016.616 URL http://www.sciencedirect.com/science/article/pii/S0045782509003612617 [14] M. Tavelli, M. Dumbser, A pressure-based semi-implicit space-618 time discontinuous Galerkin method on staggered unstructured619 meshes for the solution of the compressible Navier-Stokes equations620 at all Mach numbers, J. Comput. Phys. 341 (2017) 341 – 376.621 doi:https://doi.org/10.1016/j.jcp.2017.03.030.622 URL http://www.sciencedirect.com/science/article/pii/S0021999117302255623 [15] M. E. V´azquez-Cend´on, Solving Hyperbolic Equations with Finite Vol-624 ume Methods, Springer, 2015. doi:10.1007/978-3-319-14784-0.625 [16] L. Cea, J. R. French, M. E. V´azquez-Cend´on, Numerical modelling of626 tidal flows in complex estuaries including turbulence: An unstructured627 finite volume solver and experimental validation, Int. J. Numer. Meth.628 Engng. 67 (2006) 1909–1932. doi:10.1002/nme.1702.629 URL http://onlinelibrary.wiley.com/doi/10.1002/nme.1702/full630 [17] L. Saavedra, Modelizaci´on matem´atica y resoluci´on num´erica de prob-631 lemas de combusti´on de carb´on pulverizado, Ph.D. thesis, Departa-632 mento de Matem´atica Aplicada Universidade de Santiago de Compostela633 (2011).634 [18] P. Lax, B. Wendroff, Systems of conservation laws, Commun. Pur. Appl.635 Math. 13 (2) (1960) 217–237. doi:10.1002/cpa.3160130205.636 [19] P. D. Lax, Hyperbolic systems of conservation laws II, Commun. Pur.637 Appl. Math. 10 (4) (1957) 537–566. doi:10.1002/cpa.3160100406.638 [20] S. K. Godunov, A finite difference method for the computation of dis-639 continuous solutions of the equations of fluid dynamics, Mat. Sb. 47640 (1959) 357–393.641 [21] V. Kolgan, Application of the principle of minimizing the derivative642 to the construction of finite-difference schemes for computing discon-643 tinuous solutions of gas dynamics, J. Comput. Phys. 230 (7) (2011)644 2384–2390. doi:10.1016/j.jcp.2010.12.033.645 URL http://www.sciencedirect.com/science/article/pii/S0021999110007072646 36 [22] C. Berthon, Why the MUSCL-Hancock scheme is L1-stable, Numer.647 Math. 104 (2006) 27–46. doi:10.1007/s00211-006-0007-4.648 [23] L. Cea, M. E. V´azquez-Cend´on, Analysis of a new Kolgan-type scheme649 motivated by the shallow water equations, Appl. Num. Math. 62 (4)650 (2012) 489–506. doi:10.1016/j.apnum.2011.06.002.651 [24] B. van Leer, On the relation between the upwind-differencing schemes652 of Godunov, Engquist-Osher and Roe, SIAM J. Sci. Stat. Comp. 5 (1)653 (1984) 1–20. doi:10.1137/0905001.654 [25] B. van Leer, Towards the ultimate conservative difference scheme, J.655 Comp. Phys. 135 (2) (1997) 229–248, doi:10.1006/jcph.1997.5704.656 [26] P. K. Sweby, High resolution schemes using flux limiters for hyper-657 bolic conservation laws, SIAM J. Num. Anal. 21 (5) (1984) 995–1011.658 doi:10.1137/0721062.659 [27] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly660 high order accurate essentially non-oscillatory schemes, III, in: Up-661 wind and High-Resolution Schemes, Springer, 1987, pp. 218–290.662 doi:10.1006/jcph.1996.5632.663 [28] C. E. Castro, E. F. Toro, Solvers for the high-order Riemann problem664 for hyperbolic balance laws, J. Comp. Phys. 227 (4) (2008) 2481–2513.665 doi:10.1016/j.jcp.2007.11.013.666 [29] C.-W. Shu, S. Osher, Efficient implementation of essentially non-667 oscillatory shock-capturing schemes, J. Comp. Phys. 77 (2) (1988) 439–668 471. doi:10.1016/0021-9991(88)90177-5.669 [30] X.-D. Liu, S. Osher, T. Chan, Weighted essentially non-670 oscillatory schemes, J. Comp. Phys. 115 (1) (1994) 200–212.671 doi:10.1006/jcph.1994.1187.672 URL http://www.sciencedirect.com/science/article/pii/S0021999184711879673 [31] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cam-674 bridge Texts in Applied Mathematics, 2002.675 [32] E. F. Toro, R. C. Millington, L. A. M. Nejad, Godunov meth-676 ods, Springer, 2001, Ch. Towards very high order Godunov schemes.677 doi:10.1007/978-1-4615-0663-8.678 37 [33] E. F. Toro, V. A. Titarev, ADER schemes for scalar non-linear hyper-679 bolic conservation laws with source terms in three-space dimensions,680 J. of Comp. Phys. 202 (1) (2005) 196–215. doi:10.1016/j.jcp.2004.06.014.681 URL http://www.sciencedirect.com/science/article/pii/S0021999104002736682 [34] M. Dumbser, C. D. Munz, ADER discontinuous Galerkin schemes683 for aeroacoustics, CR Acad. Sci. II B 333 (9) (2005) 683–687.684 doi:10.1016/j.crme.2005.07.008.685 URL http://www.sciencedirect.com/science/article/pii/S1631072105001191686 [35] V. Titarev, Derivative Riemann problem and ADER schemes, Ph.D.687 thesis, Universit`a degli studi di Trento (2005).688 [36] E. F. Toro, V. A. Titarev, TVD fluxes for the high-order ADER schemes,689 J. Sci. Comp. 24 (3) (2005) 285–309. doi:10.1007/s10915-004-4790-8.690 [37] V. A. Titarev, E. F. Toro, ADER schemes for three-dimensional691 non-linear hyperbolic systems, J. Comp. Phys. 204 (2) (2005) 715–736.692 doi:10.1016/j.jcp.2004.10.028.693 URL http://www.sciencedirect.com/science/article/pii/S0021999104004358694 [38] E. F. Toro, M. Dumbser, V. A. Titarev, M. K¨aser, The derivative Rie-695 mann problem: the basis for high order ADER schemes, in: ECCOMAS696 CFD 2006: Proceedings of the European Conference on Computational697 Fluid Dynamics, Egmond aan Zee, The Netherlands, September 5-8,698 2006, Delft University of Technology; European Community on Compu-699 tational Methods in Applied Sciences (ECCOMAS), 2006.700 [39] E. F. Toro, V. A. Titarev, Derivative Riemann solvers for systems of701 conservation laws and ADER methods, J. Comp. Phys. 212 (1) (2006)702 150–165. doi:10.1016/j.jcp.2005.06.018.703 URL http://www.sciencedirect.com/science/article/pii/S0021999105003141704 [40] V. A. Titarev, E. F. Toro, ADER schemes for hyperbolic conservation705 laws with reactive terms, in: ECCOMAS CFD 2006: Proceedings of the706 European Conference on Computational Fluid Dynamics, Egmond aan707 Zee, The Netherlands, September 5-8, 2006, Delft University of Tech-708 nology; European Community on Computational Methods in Applied709 Sciences (ECCOMAS), 2006.710 38 [41] Y. Takakura, Direct-expansion forms of ADER schemes for conserva-711 tion laws and their verification, J. Comp. Phys. 219 (2) (2006) 855–878.712 doi:10.1016/j.jcp.2006.05.013.713 [42] V. A. Titarev, E. F. Toro, Analysis of ADER and ADER-WAF schemes,714 IMA J Num. Anal. 27 (3) (2007) 616–630. doi:10.1093/imanum/drl033.715 [43] H. Zahran, Central ADER schemes for hyperbolic conser-716 vation laws, J. Math. Anal. Appl. 346 (1) (2008) 120–140.717 doi:10.1016/j.jmaa.2008.05.032.718 URL http://www.sciencedirect.com/science/article/pii/S0022247X08005192719 [44] M. Dumbser, D. S. Balsara, E. F. Toro, C.-D. Munz, A unified frame-720 work for the construction of one-step finite volume and discontinuous721 Galerkin schemes on unstructured meshes, J. Comput. Phys. 227 (18)722 (2008) 8209–8253. doi:https://doi.org/10.1016/j.jcp.2008.05.025.723 URL http://www.sciencedirect.com/science/article/pii/S0021999108002829724 [45] M. Dumbser, C. Enaux, E. F. Toro, Finite volume schemes725 of very high order of accuracy for stiff hyperbolic bal-726 ance laws, J. Comput. Phys. 227 (8) (2008) 3971 – 4001.727 doi:https://doi.org/10.1016/j.jcp.2007.12.005.728 URL http://www.sciencedirect.com/science/article/pii/S0021999107005578729 [46] M. Dumbser, Arbitrary high order PNPM schemes on unstructured730 meshes for the compressible Navier-Stokes equations, Comput. Fluids731 39 (1) (2010) 60–76. doi:10.1016/j.compfluid.2009.07.003.732 URL http://www.sciencedirect.com/science/article/pii/S0045793009001030733 [47] A. Hidalgo, M. Dumbser, ADER schemes for nonlinear systems of stiff734 advection–diffusion–reaction equations, J. Sci. Comput. 48 (1-3) (2011)735 173–189. doi:10.1007/s10915-010-9426-6.736 [48] W. Boscheri, M. Dumbser, A direct arbitrary-lagrangian-eulerian737 ADER-WENO finite volume scheme on unstructured tetrahedral738 meshes for conservative and non-conservative hyperbolic systems in 3D,739 J. Comput. Phys. 275 (2014) 484–523. doi:10.1016/j.jcp.2014.06.059.740 URL http://www.sciencedirect.com/science/article/pii/S002199911400477X741 39 [49] G. I. Montecinos, E. F. Toro, Reformulations for general advection–742 diffusion–reaction equations and locally implicit ADER schemes, J.743 Comput. Phys. 275 (2014) 415–442. doi:10.1016/j.jcp.2014.06.018.744 [50] S. Busto, E. F. Toro, M. E. V´azquez-Cend´on, Design and analisis of745 ADER–type schemes for model advection–diffusion–reaction equations,746 J. Comp. Phys. 327 (2016) 553–575. doi:10.1016/j.jcp.2016.09.043.747 [51] A. Berm´udez, Continuum thermomechanics, Vol. 43 of Progress in748 Mathematical Physics, Birkh¨auser Verlag, Basel, 2005.749 [52] T. Chac´on, R. Lewandowski, Mathematical and numerical foundations750 of turbulence models and applications, Modeling and simulation in sci-751 ence, engineering and technology, Birkhauser, 2014. doi:10.1007/978-1-752 4939-0455-6.753 [53] V. V. Rusanov, The calculation of the interaction of non-stationary754 shock waves and obstacles, USSR Computational Mathematics and755 Mathematical Physics 1 (1962) 304–320.756 [54] L. Cea, M. E. V´azquez-Cend´on, Unstructured finite volume dis-757 cretization of two-dimensional depth-averaged shallow water equations758 with porosity, Int. J. Numer. Methods Fluids 63 (8) (2010) 903–930.759 doi:10.1002/fld.2107.760 URL http://onlinelibrary.wiley.com/doi/10.1002/fld.2107/full761 [55] R. Bermejo, L. Saavedra, Modified Lagrange-Galerkin methods of first762 and second order in time for convection-diffusion problems, Numer.763 Math. 120 (2012) 601–638. doi:10.1007/s00211-011-0418-8.764 [56] M. Sch¨afer, S. Turek, F. Durst, E. Krause, R. Rannacher, Benchmark765 Computations of Laminar Flow Around a Cylinder, Vieweg+Teubner766 Verlag, Wiesbaden, 1996, pp. 547–566. doi:10.1007/978-3-322-89849-4-767 39.768 URL http://dx.doi.org/10.1007/978-3-322-89849-4-39769 [57] J. Volker, Higher order finite element methods and multigrid solvers in a770 benchmark problem for the 3D Navier-Stokes equations, Int. J. Numer.771 Methods Fluids 40 (6) (2002) 775–798. doi:10.1002/fld.377.772 URL http://dx.doi.org/10.1002/fld.377773 40 Appendix A. LADER774 ADER methodology was successfully extended in [50] to solve advection-775 diffusion-reaction equations. The developed method, ADER-ADRE, is of776 second order in space and time and stability can be ensured by determining777 the time step taking into account the advection, diffusion and reaction coef-778 ficients. Despite this method is easily programmed for the one dimensional779 code, the computation of the fluxes focusing on a particular finite volume780 at each time is not suitable for the mesh structure we have in the three-781 dimensional case. In this case, we would like to take profit from the loop782 on the faces of the finite volumes and reduce the computational cost. To do783 that, the LADER method, which preserves the second order and the stability784 of the ADER-ADRE method, was developed.785 Let us consider the advection-diffusion-reaction equation786 ∂tq(x, t) + λ∂xq(x, t) = ∂x(α(x, t)∂xq(x, t)) + βq(x, t) (A.1) where787  q(x, t) is the conservative variable,788  x, t are the spatial and temporal independent variables,789  λis the characteristic speed,790  α(x, t) is the diffusion coefficient, a prescribed function,791  βis the coefficient of the reaction term.792 Then, LADER method is divided into the following steps:793 Step 1. Polynomial reconstruction. We consider a reconstruction of the794 data in terms of piecewise first-degree polynomials of the form795 pi(x) =      pi L(x) = qn i+ ∆n i L(x−xi), x ∈xi−1 2, xii, pi R(x) = qn i+ ∆n i R(x−xi), x ∈hxi, xi+1 2, (A.2) where ∆n i L, ∆n i R denote the approximations of the spatial derivatives796 of q(x, t) at time tnrelated to two auxiliary elements of volume Ci=797 hxi−1 2, xi+1 2i:798 Ti−1i L = [xi, xi+1], Tii+1 R= [xi−1, xi] (A.3) (see Figure A.13).799 41 qi−2Rqi−1Lqi−1Lqi+1 Lqi+2 Lqi+3 L qi−1Rqi R qi+1 Rqi+2 R i−1i i + 1 i+ 2 xi− 3 2xi− 1 2xi+1 2xi+3 2xi+5 2 Ti−1i R Ti−1i L =Ti+1 i+2 RTi+1 i+2 L Ti i+1 RTi i+1 L Figure A.13: Mesh and nomenclature. Step 2. Solution of the generalized Riemann problem (GRP). To construct800 the numerical flux at xi+1 2the following generalizations of the Classical801 Riemann Problem are made. On the one hand, the initial condition is802 assumed to be a piecewise first-degree polynomial. On the other hand,803 the partial differential equation accounts for the diffusion and reaction804 terms. That leads to the problem805    ∂tq(x, t) + λ∂xq(x, t) = ∂x(α∂xq) (x, t) + βq (x, t), q(x, 0) = pi R(x), x < 0, pi+1 L(x), x > 0. (A.4) Step 3. Diffusion and reaction terms. These terms are computed by approx-806 imating the integrals by the mid-point rule in both space and time.807 The solution of the GRP at the interface xi+1 2, set in Step 2, is expressed808 as a Taylor series expansion in time, namely,809 qi+1 2=q(0,0+) + τ∂tq(0,0+).(A.5) The first term of the above equation is computed as the solution of the810 classical Rieman problem811    ∂tq(x, t) + λ∂xq(x, t)=0, q(x, 0) = qi R, x < 0, qi+1 L, x > 0, (A.6) 48 fu2(x, y, z, t) = πz3sin(πtz3)−πt sin(πt(x+y+z)) −6πtzµ sin(πtz3) −9π2t2z4µcos(πtz3)+3πtz2exp(−2πt2x) sin(πtz3),(B.2) fu3(x, y, z, t) = −πt sin(πt(x+y+z)) −4π2t4µexp(−2πt2x) −4πtx exp(−2πt2x) −2πt2sin(πty) exp(−2πt2x) cos(πtz).(B.3)  Manufactured test 2. Turbulent flow with species transport.862 fu1(x, y, z, t) = (2πt cos(πtx))/3−πt sin(πt(x+y+z)) + πy cos(πty) cos(πtz) −πz sin(πty) sin(πtz)+2π2t2µsin(πty) cos(πtz) −πt cos(πtz3) cos(πty) cos(πtz)−πt sin(πty) sin(πtz) exp(−2πt2x) +(2π2Cµt2sin(πty) cos(πtz)(sin(πtx) + 2)2)/(exp(−πtz) + 1) +(π2Cµt2sin(πty) sin(πtz) exp(−πtz) (sin(πtx) + 2)2)/(exp(−πtz) + 1)2,(B.4) fu2(x, y, z, t) = πz3sin(πtz3)−πt sin(πt(x+y+z)) −6πtzµ sin(πtz3) 9π2t2z4µcos(πtz3)+3πtz2exp(−2πt2x) sin(πtz3) −(9π2Cµt2z4cos(πtz3)(sin(πtx) + 2)2)/(exp(−πtz) + 1) −(6πCµtz sin(πtz3)(sin(πtx) + 2)2)/(exp(−πtz) + 1) −(3π2Cµt2z2exp(−πtz) sin(πtz3) (sin(πtx) + 2)2)/(exp(−πtz) + 1)2,(B.5) fu3(x, y, z, t) = (4π2Cµt3exp(−2πt2x) cos(πtx)(sin(πtx) + 2))/(exp(−πtz) + 1) −4π2t4µexp(−2πt2x)−4πtx exp(−2πt2x) −2πt2sin(πty) exp(−2πt2x) cos(πtz) −(4π2Cµt4exp(−2πt2x)(sin(πtx) + 2)2)/(exp(−πtz) + 1) −πt sin(πt(x+y+z)),(B.6) fk(x, y, z, t) = exp(−πtz) + πx cos(πtx) −(Cµ(sin(πtx) + 2)2(2πt2exp(−2πt2x) +πt sin(πty) sin(πtz))2)/(exp(−πtz) + 1) +π2t2µsin(πtx) + πt sin(πty) cos(πtx) cos(πtz) −(9π2Cµt2z4sin(πtz3)2(sin(πtx) + 2)2)/(exp(−πtz) + 1) −(π2Cµt2cos(πty)2cos(πtz)2(sin(πtx) + 2)2)/(exp(−πtz) + 1) −(2π2Cµt2cos(πtx)2(sin(πtx) + 2))/(σk(exp(−πtz) + 1)) +(π2Cµt2sin(πtx)(sin(πtx) + 2)2)/(σk(exp(−πtz) + 1)) + 1,(B.7) fε(x, y, z, t)=(C2ε(exp(−πtz) + 1)2)/(sin(πtx) + 2) −πz exp(−πtz) −π2t2µexp(−πtz) 49 −(C1ε(exp(−πtz) + 1)((Cµ(sin(πtx) + 2)2(2πt2exp(−2πt2x) +πt sin(πty) sin(πtz))2)/(exp(−πtz) + 1) +(9π2Cµt2z4sin(πtz3)2(sin(πtx) + 2)2)/(exp(−πtz) + 1) +(π2Cµt2cos(πty)2cos(πtz)2(sin(πtx) +2)2)/(exp(−πtz) + 1)))/(sin(πtx) + 2) −πt exp(−πtz) exp(−2πt2x) −(π2Cµt2exp(−πtz)(sin(πtx) + 2)2)/(σε(exp(−πtz) + 1)) +(π2Cµt2exp(−2πtz)(sin(πtx) + 2)2)/(σε(exp(−πtz) + 1)2),(B.8) fy(x, y, z, t) = πx cos(πtx) + π2t2sin(πtx)(D+ (Cµ(sin(πtx) +2)2)/(Sct(exp(−πtz) + 1))) +πt sin(πty) cos(πtx) cos(πtz) −(2π2Cµt2cos(πtx)2(sin(πtx) + 2))/(Sct(exp(−πtz) + 1)).(B.9) 863