A semi-implicit hybrid finite volume/finite element scheme for all Mach number flows on staggered unstructured meshes

A semi-implicit hybrid nite volume / nite element scheme for all Mach number ows on staggered unstructured meshes S. Busto a, ∗ , L. Río-Martín a,b , M.E. Vázquez-Cendón b , M. Dumbser a a Laboratory of Applied Mathematics, DICAM, University of Trento, via Mesiano 77, 38123 Trento, Italy b Department of Applied Mathematics, University of Santiago de Compostela, 15782 Santiago de Compostela, Spain Abstract In this paper a new hybrid semi-implicit nite volume / nite element (FV/FE) scheme is presented for the numerical solution of the compressible Euler and Navier-Stokes equations at all Mach numbers on unstructured staggered meshes in two and three space dimensions. The chosen grid arrangement consists of a primal simplex mesh composed of triangles or tetrahedra, and an edge-based / face-based staggered dual mesh. The governing equations are discretized in conservation form. The nonlinear convective terms of the equations, as well as the viscous stress tensor and the heat ux, are discretized on the dual mesh at the aid of an explicit local ADER nite volume scheme, while the implicit pressure terms are discretized at the aid of a continuous P1 nite element method on the nodes of the primal mesh. In the zero Mach number limit, the new scheme automatically reduces to the hybrid FV/FE approach forwarded in [1] for the incompressible Navier-Stokes equations. As such, the method is asymptotically consistent with the incompressible limit of the governing equations and can therefore be applied to ows at all Mach numbers. Due to the chosen semi-implicit discretization, the CFL restriction on the time step is only based on the magnitude of the ow velocity and not on the sound speed, hence the method is computationally ecient at low Mach numbers. In the chosen discretization, the only unknown is the scalar pressure eld at the new time step. Furthermore, the resulting pressure system is symmetric and positive denite and can therefore be very eciently solved with a matrix-free conjugate gradient method. In order to assess the capabilities of the new scheme, we show computational results for a large set of benchmark problems that range from the quasi incompressible low Mach number regime to compressible ows with shock waves. Keywords: all Mach number ow solver, pressure-based projection method, nite element method, nite volume scheme, semi-implicit scheme on unstructured staggered meshes, ADER methodology 1. Introduction 1 Since their rst formulation more than 200 years ago, the Euler and Navier-Stokes equations describing 2 the ow of inviscid and viscous uids have always been a big challenge, both from the theoretical as well 3 as from the numerical point of view. The Euler equations can be directly derived from rst principles by 4 considering the conservation of mass, momentum, and total energy. Their extension to the Navier-Stokes 5 equations is then achieved at the aid of appropriate assumptions for the viscous stress tensor and the 6 heat ux. In the most general case, the uid is assumed to be compressible , but dierent ow regimes 7 can be identied at the aid of the dimensionless Mach number M=∥v∥/c , where v and c are the uid 8 velocity and the sound speed, respectively. For M→0 the behaviour of the uid becomes the one of an 9 incompressible medium, with the well-known condition ∇·v= 0 , which states that the velocity eld must 10 ∗ Corresponding author Email addresses: saray.bustoul[email protected] (S. Busto), laura.delr[email protected] (L. Río-Martín), elena.vazqu[email protected] (M.E. Vázquez-Cendón), michael.dumbs[email protected] (M. Dumbser) Preprint submitted to Applied Mathematics and Computation December 13, 2024 2 become divergence-free when the ow becomes incompressible in the limit M→0 . This asymptotic limit was 11 rigorously studied for the rst time by Klainerman and Majda in [2, 3]. The asymptotic analysis shows that 12 in the incompressible limit and without compression from the boundary, the pressure can be decomposed 13 into two dierent contributions: a spatially constant part of the pressure satisfying the equation of state 14 and some uctuations of the pressure around that constant, governed by the well-known elliptic pressure 15 Poisson equation. This change of behaviour for M→0 is important since the original governing equations 16 are hyperbolic-parabolic they can even exhibit shock waves for high Mach numbers. Because of this changing 17 behaviour of the equations according to the Mach number, it is notoriously dicult to construct suitable 18 numerical schemes which can be simultaneously applied to compressible high Mach number ows with shock 19 waves and also to incompressible or nearly incompressible low Mach number ows. 20 Typically, the incompressible Euler and Navier-Stokes equations are solved via semi-implicit pressure21 based schemes of the nite dierence type on staggered grids, see e.g. [4, 5, 6, 7, 8, 9, 10, 11, 12], or at 22 the aid of continuous nite elements [13, 14, 15, 16, 17, 18, 19]. Instead, for the simulation of compressible 23 ows at higher Mach numbers and with shock waves, explicit density-based nite volume schemes of the 24 Godunov-type on collocated grids are usually more popular, see [20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. 25 A rst attempt to generalize semi-implicit methods to the more general case of compressible ows was 26 made by Casulli and Greenspan in [30], but the proposed scheme was not conservative and therefore could 27 not be used for the treatment of high Mach number ows with shock waves. Semi-implicit schemes that 28 explicitly make use of the low Mach number asymptotics of the governing partial dierential equations 29 can be found in [31, 32, 33, 34], while the rst conservative staggered semi-implicit pressure-based scheme 30 for compressible ows was introduced by Park and Munz in [35]. The scheme [35] can be considered as 31 one of the rst all Mach number ow solvers ever proposed in the literature. The particular splitting of 32 explicit convective terms and implicit pressure terms used in [35] was later studied in more detail in [36] in 33 order to construct a novel ux-vector splitting method. Since the pioneering work of Park and Munz, the 34 development of all Mach number ow solvers, i.e., of numerical schemes that work at the same time for high 35 Mach number ows with shock waves and in the incompressible limit of the equations, has become a very 36 active research eld with many relevant contributions, see e.g. [37, 38, 39, 40, 41, 42, 43, 44, 45, 46] and 37 references therein. For special low Mach number corrections to explicit density-based nite volume schemes, 38 the reader is referred to [47, 48]. 39 On unstructured simplex meshes, classical continuous nite element methods can nowadays be considered 40 as standard for the numerical solution of the incompressible Navier-Stokes equations. Instead, the construc41 tion of discontinuous Galerkin nite element schemes for the solution of the compressible and incompressible 42 Navier-Stokes equations on unstructured meshes is still the topic of ongoing research. For an overview of 43 high order DG schemes for the compressible and incompressible Navier-Stokes equations, see for example 44 [49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], but this list does not pretend to be complete. Concerning high 45 order semi-implicit discontinuous Galerkin methods on collocated grids we refer to [61, 62, 63], while a new 46 family of semi-implicit staggered discontinuous Galerkin schemes for the discretization of the incompressible 47 and compressible Navier-Stokes equations was recently forwarded in [64, 65, 66, 67, 68, 69]. 48 To round-up this brief literature review, we would also like to point the reader to a very recent and 49 completely dierent approach for the solution of the Navier-Stokes equations, which consists in embedding 50 the Navier-Stokes equations in a more general rst order hyperbolic system with sti relaxation source terms 51 that is able to describe continuum mechanics as a whole, from nonlinear elastic solids over visco-plastic solids 52 to Newtonian and non-Newtonian uids, and from which for small enough relaxation times the Navier-Stokes 53 equations are retrieved in the limit of a much more general model that contains continuum mechanics as a 54 whole, see [70, 71, 72, 73, 74]. This universal model is based on the pioneering work of Godunov and Romenski 55 on symmetric hyperbolic and thermodynamically compatible systems and on nonlinear hyperelasticity, see 56 e.g. [75, 76, 77, 78, 79] and references therein. For an alternative hyperbolic relaxation approach for the 57 discretization of the Navier-Stokes equations, the reader is referred to [80]. 58 The numerical methods previously discussed were all of a specic type, say nite volume, nite dierence, 59 or nite element schemes. More recently in a series of papers a new class of hybrid nite volume / continuous 60 nite element methods on staggered unstructured meshes in 2D and 3D has been proposed in [81, 1, 82, 83] 61 for the solution of the incompressible Navier-Stokes equations and for the low Mach number limit of the 62 3 weakly compressible equations. In these hybrid schemes, the nonlinear convective part of the equations was 63 solved at the aid of an explicit nite volume scheme on an edge-based staggered dual mesh, see [84, 1] for a 64 more detailed analysis, and the pressure equation was solved with a continuous nite element method on the 65 primal grid. The advantage of this hybrid approach is that for each part of the governing PDE system the 66 most appropriate numerical method could be used, since it is well-known that explicit nite volume methods 67 are more suitable for the discretization of nonlinear hyperbolic PDE systems, while the clear strength of 68 continuous nite element methods lies in the discretization of elliptic problems. 69 It is therefore the aim of the present paper to provide a novel pressure-based semi-implicit hybrid nite 70 element / nite volume method on staggered unstructured meshes that can solve the compressible Euler 71 and Navier-Stokes equations in a wide range of Mach numbers, which is a very substantial generalization 72 compared to the incompressible and weakly-compressible ow solvers presented in [81, 1, 83]. Following 73 the seminal ideas outlined in [35, 36, 39], the nonlinear convective part of the equations will be discretized 74 via an explicit nite volume scheme, while the resulting pressure equation, which is more complex than 75 the simple pressure Poisson equation that typically results from the discretization of the incompressible 76 Navier-Stokes equations, is discretized on the primal mesh at the aid of classical continuous nite elements. 77 The semi-implicit discretization allows choosing a time step that is not limited by the sound speed, but only 78 by the velocity magnitude. The new hybrid scheme of this paper is designed to work simultaneously for 79 incompressible and low Mach number ows, as well as for compressible ows including shock waves. For 80 M→0 the scheme reduces to the hybrid FV/FE method for the incompressible Navier-Stokes equations 81 forwarded in [1]. As such, the proposed method is an asymptotic preserving (AP) all Mach number ow 82 solver. 83 The rest of the paper is organized as follows: in Section 2 we rst introduce the governing equations 84 considered in this paper; next, in Section 3 we present their discretization via the new semi-implicit hybrid 85 nite-volume / nite-element scheme on staggered meshes. In Section 4 we present numerical results for a 86 wide range of Mach numbers, from almost incompressible ows to supersonic ows with shock waves. The 87 conclusions and an outlook to further work are given in Section 5. 88 2. Governing partial dierential equations 89 Let us denote by ρ the density, u= (u, v, w) is the velocity vector and E is the specic total energy then, 90 the compressible Navier-Stokes equations given in conservative form read 91 ∂ρ ∂t +∇·(ρu)=0, (1) ∂ρu ∂t +∇·(ρu⊗u) + ∇p−∇·τ=ρg, (2) ∂ρE ∂t +∇·[u(ρE +p)] −∇·(τu) + ∇·q=ρg·u, (3) where g is the gravity vector, τ is the tensor of the viscous stresses, 92 τ=µ∇u+∇uT−2 3µ(∇·u)I, (4) and q denotes the heat ux, 93 q=−λ∇θ. (5) Here, θ is the temperature and λ denotes the thermal conductivity. In this paper we use the simple ideal 94 gas equation of state (EOS) to close the system: 95 p=ρRθ, (6) 4 where R=cp−cv is the specic gas constant, with cp being the heat capacity at constant pressure, while 96 cv denotes the heat capacity at constant volume. Accounting for (6), the relation between the total energy, 97 the kinetic energy k and the specic internal energy e reads 98 ρE =ρe +ρk =1 γ−1p+1 2ρ|u|2 (7) with γ=cp cv being the ratio of specic heats. Introducing the enthalpy, 99 h=e+p ρ=γ γ−1 p ρ (8) we get 100 ∂ρE ∂t +∇·(ρku) + ∇·(ρhu)−∇·(τu) + ∇·q=ρg·u. (9) The former system, (1)-(3), can be rewritten in terms of the conservative variables, w= (wρ,wu, wE)T= 101 (ρ, ρu, ρE)T as 102 ∂wρ ∂t +∇·(wu) = 0, (10) ∂wu ∂t +∇·1 ρwu⊗wu+∇p−∇·τ=ρg, (11) ∂wE ∂t +∇·1 ρwu(wE+p)−∇·1 ρτwu+∇·q=g·wu. (12) 3. Numerical method 103 Discretization of system (10)-(12) is performed extending the hybrid nite volume / nite element method 104 proposed in [81, 85, 1, 86]. We start by considering a semi-discrete scheme where only time discretization is 105 applied leading to 106 1 ∆tWn+1 ρ−Wn ρ+∇·(Wn u)=0, (13) 1 ∆tWn+1 u−Wn u+∇·1 ρnWn u⊗Wn u+∇Pn+1 −∇·Tn=ρng, (14) 1 ∆tWn+1 E−Wn E+∇·(KnWn u) + ∇·Hn+1Wn+1 u−∇·1 ρnTnWn u+∇·Qn=g·Wn u, (15) with Wn, Pn approximations of the solution, w(x, tn), p (x, tn) , at time tn∈R+ and x∈Rd the spa107 tial coordinate. We now introduce the following notation for an intermediate approximation of the linear 108 momentum 109 f W u=Wn u−∆t∇·1 ρnWn u⊗Wn u−∇·Tn−ρng (16) and dene 110 f f Wu:= f W u−∆t∇Pn, δPn+1 := Pn+1 −Pn (17) so 111 f W u=f f Wu+ ∆t∇Pn, (18) 5 Wn+1 u=f W u−∆t∇Pn+1 =f f Wu−∆t∇δPn+1. (19) In such a way, we have derived a pressure-correction formulation in which the computation of the pressure 112 and the linear momentum are decoupled". Similarly, we can dene an intermediate auxiliary variable for 113 the computation of the total energy, 114 f WE=Wn E−∆t∇·(KnWn u)−∇·1 ρnTnWn u+∇·Qn−g·Wn u. (20) Later, Wn+1 E would be recovered from 115 Wn+1 E=f WE−∆t∇·Hn+1Wn+1 u. (21) On the other hand, equation (15) can be rewritten in terms of the pressure and the kinetic energy by using 116 relation (7) and the ideal gas equation of state as follows: 117 Pn+1 γ−1+ (ρK)n+1 =Pn γ−1−Pn γ−1+f WE−∆t∇·Hn+1Wn+1 u (22) Substitution of (19) yields 118 Pn+1 γ−1−Pn γ−1=−(ρK)n+1 +f WE−Pn γ−1−∆t∇·Hn+1 f f Wu+ ∆t2∇·Hn+1 ∇δPn+1. (23) Hence, 119 1 γ−1δPn+1 −∆t2∇·Hn+1 ∇δPn+1=f WE−Pn γ−1−(ρK)n+1 −∆t∇·Hn+1 f f Wu. (24) A Picard procedure is applied to deal with the crossed tn+1 terms, i.e., Pn+1 in equation (19), (ρK)n+1 = 120 1 2ρn+1 Wn+1 u2 and Hn+1 in (24), and Wn+1 u in (21), since we do not want to solve a highly nonlinear 121 system. The nal system of equations to be discretized in space reads 122 Wn+1 ρ=Wn ρ−∆t∇·(Wn u), (25) f f Wu=Wn u−∆t∇·1 ρnWn u⊗Wn u+∇Pn−∇·Tn−ρng, (26) f WE=Wn E−∆t∇·(KnWn u)−∇·1 ρnTnWn u+∇·Qn−g·Wn u, (27) 1 γ−1δPn+1,k+1 −∆t2∇·Hn+1,k ∇δPn+1,k+1=f WE−Pn γ−1−(ρK)n+1,k −∆t∇·Hn+1,k f f Wu, (28) Wn+1,k+1 u=f f Wu−∆t∇δPn+1,k+1, (29) Pn+1,k+1 =Pn+δPn+1,k+1, (30) Wn+1 E=f WE−∆t∇·Hn+1,k+1Wn+1,k+1 u, (31) with k the Picard iteration index, k= 1, . . . , NPic . 123 Let us remark that the method is by construction asymptotic preserving in the low Mach number limit. 124 In the limit M→0 , we have c2→ ∞ with c2=H γ−1 . According to [2, 3], the pressure, and for constant 125 density also the enthalpy, tend to a constant. In this limit, we can now divide equation (28) by the enthalpy 126 and neglecting terms of the order 1/c2 we obtain the following equation 127 ∇2δPn+1 =1 ∆t∇·f f Wu, (32) 6 which, together with the momentum equation (26), corresponds to the pressure correction system obtained 128 for the incompressible Navier-Stokes equations in [1]. 129 For general equations of state, relation (24) needs to be replaced by 130 WE(Pn+1, W n+1 ρ)−∆t2∇·Hn+1 ∇Pn+1=f WE−(ρK)n+1 −∆t∇·Hn+1 f f Wu, (33) where the density at the new time tn+1 is readily available from eqn. (25), and thus the only unknown remains 131 the scalar pressure eld Pn+1 at the new time, see also [39]. When appropriate mass-lumping is used within 132 a nite-element discretization of (33), the resulting mildly nonlinear pressure system can be solved very 133 eciently at the aid of the (nested) Newton-type methods of Brugnano and Casulli [87, 88, 89] and Casulli 134 and Zanolli [90, 91], and for which convergence has been rigorously proven. For nite elements without mass 135 lumping, i.e., for non-diagonal mass matrices, the theorems which are the basis of the convergence proofs 136 in the aforementioned works of Casulli et al. do not directly apply and still need to be generalized to more 137 general non-diagonal but symmetric positive denite mass matrices. In the rest of this paper, we therefore 138 assume that the simple ideal gas equation of state holds. 139 Spatial discretization is done by choosing a numerical method adapted to the nature of each equation: 140 nite volumes are applied to approximate the transport-diusion equations, whereas the Poisson problem 141 is solved using continuous nite elements. The use of staggered grids avoids the checker-board phenomena, 142 which are typical for many numerical methods on collocated grids. The use of unstructured meshes increases 143 the applicability of the methodology with respect to Cartesian grids since the meshing of complex domains 144 becomes straightforward. The overall algorithm can be divided into four main stages: 145  Transport-diusion stage. The equations (25), (26), and (27) are solved using explicit nite volumes in 146 the dual mesh. To attain second order in space and time, a local ADER method combined with an ENO 147 reconstruction is considered. Within this stage we get the new density, ρn+1 , and the intermediate 148 approximations of the momentum, f f Wu , as well as the total energy density, f WE at each cell of the dual 149 mesh. 150  Pre-projection stage. The intermediate states for the total energy density and for the linear momentum 151 are transferred from the dual to the primal grid. Next, the auxiliary variables that will be needed within 152 the next stage, as the enthalpy, Hn+1,k , and the kinetic energy density, (ρK)n+1,k , are also computed. 153 Let us note that the intermediate variables are calculated only once per time step, meanwhile the 154 auxiliary variables are updated at each Picard iteration. 155  Projection stage. A P1 nite element scheme is employed in order to solve the pressure equation (28) 156 implicitly. The resulting δPn+1,k+1 is computed on the vertexes of the primal simplex mesh. 157  Post-projection stage. The pressure correction at the new time tn+1 is substituted in (29) and (30) to 158 update the linear momentum Wn+1,k+1 u and the pressure Pn+1,k+1 . Once the Picard iterations have 159 nished, the total energy, Wn+1 E , is recovered following (31). 160 In what follows, we will further detail each stage of the algorithm. 161 3.1. Staggered unstructured mesh 162 In this paper we make use of two overlapping unstructured staggered meshes to discretize the domain 163 Ω . For the sake of simplicity, we focus here on the 2D case introducing the main notation needed. Further 164 details on the construction of three-dimensional face-type staggered meshes can be found in [81, 65, 1, 69, 86]. 165 Let us consider a triangular primal mesh {Tk, k = 1, . . . , nel} with vertex {Vj, j = 1, . . . , nver} (Figure 1 166 left). We now dene the two triangles with basis one interior edge of a primal element and opposite vertex 167 the barycenters, B , B′ , of the two primal elements sharing this face. The dual element, Ci , is then built by 168 merging these two triangles (Figure 1 center). Similarly, a dual boundary element is a triangle which has 169 as basis a primal boundary edge and as opposite vertex the barycenter of the primal element. Let us note 170 that on the dual mesh the nodes {Ni, i = 1, . . . , nnod} are associated with the edges/faces of the simplex 171 elements on the primal mesh. The remaining notation related to the mesh is as follows: 172 7  Ki is the set of neighbouring nodes of a node Ni consisting of the barycenters of the dual cells sharing 173 a face with Ci . 174  Γi is the boundary of a cell Ci and e ηi its outward unit normal. 175  Γij is the edge between cells Ci and Cj . Nij is the barycenter of Γij (Figure 1 right). Note that 176 Γi=[ Nj∈Ki Γij . 177  |Ci| is the area of Ci . 178  e ηij is the outward unit normal vector to Γij . We dene ηij := e ηij||ηij|| , where, ||ηij|| represents the 179 length of Γij . 180 Sketches of the 2D and 3D staggered meshes are depicted in Figures 1 and 2, respectively. Tk Tl Tm V1V2 V3 V4 V5 Tk Tl Tm V1V2 V3 V4 V5 B B0 B00 Tk Tl Tm V1V2 V3 V4 V5 Figure 1: Construction of face-type dual elements from a 2D primal mesh. Left: primal mesh made of elements Tk , Tl , Tm and vertex Vn, n = 1,...,5 . Center: dual interior cells Ci , Cj (shadowed in grey), white triangles correspond to boundary cells. Right: boundary face, Γij , between the dual elements Ci , Cj . 181 V2 V3 V1 V4 V′ 4 B B′ Ni V2 V3 V1 V4 B Ni Figure 2: Example of an interior nite volume (left) and of a boundary nite volume (right) which constitute the staggered dual mesh in three space dimensions. 3.2. Transport-diusion stage 182 To solve the transport dominated equations, we use a nite volume method on the dual grid. As 183 result, we will obtain the value of the averaged new density, Wn+1 ρ=ρn+1 , on each dual cell, as well as 184 8 intermediate approximations for the cell averaged linear momentum, f f Wu , and total energy density, f WE . We 185 start integrating equations (25)-(27) on each dual cell Ci and applying Gauss theorem which yields 186 ρn+1 i=ρn i−∆t |Ci|ZΓiFW ρ(Wn)e ηidS, (34) f f Wu, i =Wn u, i −∆t |Ci|ZΓiFWu(Wn)e ηidS +ZCi∇PndV −ZΓi Tne ηidS −ZCi ρngdV , (35) f WE,i =Wn E,i −∆t |Ci|ZΓiFW E(Wn)·e ηidS −ZΓi1 ρnTnWn u·e ηidS +ZΓi Qn·e ηidS −ZCi g·Wn udV . (36) where 187 FW ρ(Wn) := Wn u,FWu(Wn) := 1 ρnWn u⊗Wn u,FW E(Wn) := KnWn u (37) are the convective uxes of mass, momentum, and total energy, respectively. 188 3.2.1. Convective numerical uxes 189 The global normal ux through the boundary of a dual cell is denoted by 190 Z(Wn,˜ηi) := F(Wn)˜ηi,with F(Wn) = FW ρ(Wn),FWu(Wn),FW E(Wn)T. (38) Let us recall that the ux contribution in the energy equation accounts only for the kinetic energy density 191 contribution, ρk , instead of the total energy density, ρE . Within the ux computation the value of Kn is 192 recovered from the linear momentum and density elds, Kn=1 2ρn|Wn u|2 . The integral of the ux term on 193 Γi can be split into the sum of the integral on the cell faces, Γij , 194 ZΓiF(Wn)e ηidS = X Nj∈KiZΓij Z(Wn,˜ηij) dS, (39) and approximated using an upwind scheme to get a stable discretization. In particular, we consider a 195 modied Rusanov ux function, [92, 86], 196 ϕWn i,Wn j,ηij=ϕρWn i,Wn j,ηij,ϕuWn i,Wn j,ηij, ϕEWn i,Wn j,ηijT =1 2(Z(Wn i,ηij) + Z(Wn j,ηij)) −1 2αn RS, ij Wn j−Wn i (40) with 197 W:= (Wρ,Wu, ρK)T (41) the modied conservative variables vector, 198 αn RS, ij =αRS(Wn i,Wn j,ηij) := max Un i·ηij,Un j·ηij+cαηij (42) the maximum signal speed on the edge and cα∈R+ 0 an articial viscosity coecient that may be activated 199 on particular tests to increase the stability properties of the nal scheme when large variations of the density 200 and energy elds are encountered in the presence of small velocities. Substitution in (34)-(36) gives 201 ρn+1 i=ρn i−∆t |Ci|X Nj∈Ki ϕρWn i,Wn j,ηij, (43) 9 f f Wu, i =Wn u, i −∆t |Ci| X Nj∈Ki ϕuWn i,Wn j,ηij+ZCi∇PndV −ZΓi Tne ηidS −ZCi ρngdV  , (44) f WE,i =Wn E,i −∆t |Ci| X Nj∈Ki ϕEWn i,Wn j,ηij−ZΓi1 ρnTnWn u·e ηidS +ZΓi Qn·e ηidS −ZCi g·Wn udV  . (45) The scheme proposed above would result in a rst order scheme in space and time. To increase the order 202 of accuracy attaining second order in space and time, a local ADER methodology (LADER) is employed, see 203 [1, 93, 86]. The reader is referred to [94, 29, 95, 84] for further details on the original ADER methodology 204 and to [96, 97, 72] for an alternative variant of ADER schemes that allow to avoid the cumbersome Cauchy205 Kovalevskaya procedure thanks to the use of a general space-time nite element predictor. In what follows, 206 we briey recall the main steps to be performed in the LADER algorithm: 207 Step 1. Piecewise polynomial reconstruction in the neighbourhood of each boundary edge of the cell. Con208 sidering an scalar conservative variable, W , and one of the cell boundaries, Γij , the related reconstruc209 tion polynomials read 210 Pi ij(N) = Wi+ (N−Ni) (∇W)i ij , Pj ij(N) = Wj+ (N−Nj) (∇W)j ij . (46) To circumvent Godunov's theorem and to develop a second order scheme avoiding spurious oscil211 lations, we introduce a non linearity via the use of a nonlinear Essentially Non-Oscillatory (ENO) 212 reconstruction. Accordingly, the gradients are computed as 213 (∇W)i ij =   (∇W)TijL , if (∇W)TijL ·(Nij −Ni)≤(∇W)Tij ·(Nij −Ni), (∇W)Tij ,otherwise; 214 (∇W)j ij =   (∇W)TijR , if (∇W)TijR ·(Nij −Nj)≤(∇W)Tij ·(Nij −Nj), (∇W)Tij ,otherwise, with Tij , TijL and TijR the centered and upwind primal elements to the face Γij where the gradients 215 are computed using a Galerkin approach (Crouzeix-Raviart nite elements). An alternative to the 216 ENO-based reconstruction is the use of classical slope limiters like the Barth and Jespersen limiter 217 [98], or the minmod limiter of Roe [99, 29]. Also, a posteriori limiting strategies like the MOOD 218 approach, [100], could be used and will be part of future research. 219 Step 2. Calculation of the necessary boundaryextrapolated data in xNij of each edge/face of the FV mesh, 220 Wi Nij =pi ij(Nij) = Wi+ (Nij −Ni) (∇W)i ij , (47) Wj Nij =pj ij(Nij) = Wj+ (Nij −Nj) (∇W)j ij . (48) Step 3. Use of the mid-point rule to get a second order of accuracy approximation in time. A tempo221 ral Taylor series expansion in combination with the Cauchy-Kovalevskaya procedure, based on the 222 mass, momentum, and energy equations (10), (11), (28), are employed in order to approximate the 223 conservative variables at the time tn+∆t 2 : 224 Wi Nij =ρi Nij +W∗ Nij ,Wj Nij =ρj Nij +W∗ Nij (49) where 225 W∗ ρ Nij := −∆t 2Lij Zρ(Wi Nij ,ηij) + Zρ(Wj Nij ,ηij), (50) 16 4.2. Riemann problems 368 In this section, we analyse the performance of the proposed methodology for the compressible Euler 369 equations in presence of medium to strong shocks. We consider a two-dimensional computational domain 370 with x∈[−0.5,0.5] and a variable width depending on the number of cells in the horizontal direction so 371 that the nal elements have a good aspect ratio and a small number of layers in the y -direction to reduce 372 the computational cost of the simulation. The initial condition is dened as 373 ρ0(x) = ρLif x≤xc, ρRif x>xc;u0 1(x) = uLif x≤xc, uRif x>xc;u0 2(x)=0 p0(x) = pLif x≤xc, pRif x>xc; (82) where ρL , ρR , pL , pR , uL , uR , xc are summarized in Table 3 for the diverse tests, selected among those 374 presented in [29, 66]. The nal time of each simulation, as well as the number of mesh divisions along 375 the x -axis ( Nx ), have also been reported in Table 3. The characteristic mesh spacing is therefore equal 376 to h= 1/Nx . All tests have been run with the rst order and LADER schemes using Dirichlet boundary 377 conditions in the x -direction and periodic boundary conditions in the y -direction. The rst test analysed, Test ρLρRuLuRpLpRxctend Nx RP1 1 0.125 0 0 1 0.1 0 0.2 200 RP2 1 1 −1 1 0.4 0.4 0 0.15 300 RP3 0.445 0.5 0.698 0 3.528 0.571 0 0.14 200 RP4 5.99924 5.99242 19.5975 −6.19633 460.894 46.095 −0.2 0.035 200 RP5 1 1 −19.59745 −19.59745 1000.0 0.01 0.3 0.01 300 RP6 1 1 2 −2 0.1 0.1 0 0.8 200 Table 3: Riemann problems. Initial condition, initial position of the discontinuity, xc , nal time, tend , and number of mesh cells on x -direction, Nx , for each Riemann problem. 378 RP1, is the classical Sod problem presented for the rst time in [104]. Figure 3 shows a good agreement 379 between the numerical and the exact solution for the shock, the contact, and the rarefaction waves. Figure 3: Riemann problem 1 (Sod). 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.2 using the rst order method ( CFLc= 3.35 , cα= 1 , M≈0.93 ). 380 RP2 corresponds to a double rarefaction problem. Overall the shape of the exact solution is captured 381 even if a ner mesh would be useful to better approximate the constant contact discontinuity between the 382 two rarefactions, Figures 5-6. 383 17 Figure 4: Riemann problem 1 (Sod). 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.2 using the LADER-ENO method ( CFLc= 3.35 , cα= 1 , M≈0.93 ). Figure 5: Riemann problem 2 (Double rarefaction). 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.15 using the rst order method ( CFLc= 0.4 , cα= 2 , M≈1.37 ). Figure 6: Riemann problem 2 (Double rarefaction). 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.15 using the LADER-ENO method ( CFLc= 0.4 , cα= 2 , M≈1.37 ). 18 The third test, RP3, corresponds to the Lax shock tube and is used to assess the ability of the method to 384 deal with simple waves. The obtained results, presented in Figures 7-8, match pretty well the exact reference 385 solution. Figure 7: Riemann problem 3 (Lax). 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.14 using the rst order method on mesh M1 ( CFLc= 2.78 , M≈0.94 ). Figure 8: Riemann problem 3 (Lax). 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.14 using the LADER-ENO method on mesh M1 ( CFLc= 2.78 , M≈0.94 ). 386 The fourth Riemann problem, RP4, presents three strong discontinuities travelling to the right originated 387 from two shock colliding waves. Figures 9-10 show the solution obtained with the rst and second order 388 schemes. Note that the highly restrictive Barth and Jespersen limiter has been employed jointly with an 389 articial viscosity coecient, cα= 5 , to keep the stability of the scheme. 390 RP5 is a severe test dened as a modication of the left half of the blast problem introduced in [105]. 391 It accounts for a left rarefaction wave, a right-travelling shock wave and a stationary contact discontinuity 392 generated by an initial large pressure jump of order 105 and a small velocity variation. The second order 393 scheme has been run using two dierent limiter strategies. We observe that the minmod limiter, Figure 13, 394 is more severe on damping the oscillation appearing after the rarefaction wave on the velocity eld than the 395 ENO-based reconstruction, Figure 12, which captures better the right shock. The results obtained with the 396 rst order scheme are reported in Figure 11. The robustness of the developed methodology and its capability 397 to deal with slowly moving contact discontinuities, for really high Mach numbers, are clearly proven. 398 The last Riemann problem considered, RP6, is characterised by two shock waves travelling in opposite 399 directions. An excellent agreement with the exact solution is observed in Figures 14-15. 400 19 Figure 9: Riemann problem 4. 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.035 using the rst order method ( CFLc= 0.42 , cα= 5 , M≈1.97 ). Figure 10: Riemann problem 4. 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.035 using LADER-BJ method with reconstruction thorough primitive variables instead the conservative ones ( CFLc= 0.42 , cα= 5 , M≈1.97 ). Figure 11: Riemann problem 5. 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.01 using the rst order scheme ( CFLc= 2.7·10−3 , cα= 2 , M≈956.42 ). 20 Figure 12: Riemann problem 5.1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.01 using the LADER-ENO method ( CFLc= 2.7·10−3 , cα= 2 , M≈513.68 ). Figure 13: Riemann problem 5. 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.01 using the LADER-minmod method ( CFLc= 2.7·10−3 , cα= 2 , M≈691.58 ). Figure 14: Riemann problem 6. 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.8 using the rst order method ( CFLc= 9.35 , cα= 2 , M≈7.75 ). 21 Figure 15: Riemann problem 6. 1D cut through the numerical results along the line y= 0 for ρ , u and p at tend = 0.8 using the rst order method ( CFLc= 9.35 , cα= 2 , M≈7.29 ). 4.3. 2D circular explosion 401 The circular explosion problem presented here is based on an initial radial solution given by the Sod 402 shock tube 403 ρ0(x) = 1 if r≤0.5, 0.125 if r > 0.5,u0(x)=0, p0(x) = 1 if r≤0.5, 0.1 if r > 0.5, (83) see [29, 106, 71]. We consider the computational domain Ω=[−1,1] ×[−1,1] and periodic boundary 404 conditions everywhere. The simulation is run until time tend = 0.25 on a primal triangular mesh of 85344 405 elements. To get a reference solution, a one-dimensional PDE in the radial direction obtained from the 406 compressible Euler equations when using convenient geometrical source terms, [29], is solved using a second 407 order TVD scheme on a very ne mesh made of 10000 elements. The results obtained with the rst order 408 scheme and the LADER-ENO methodology, Figures 16-17, present a good agreement with the reference 409 solution. Figure 18 allows for a direct comparison of the solution obtained with both schemes along a 1D 410 cut. The second order LADER method with ENO reconstruction provides a better approximation of the 411 solution compared to the rst order scheme, as expected. 412 4.4. 3D spherical explosion 413 In this section, we study the behaviour of the method for a 3D spherical explosion benchmark based 414 on the Sod problem. The computational domain is dened to be the sphere of unit radius centered at the 415 origin. Initial conditions read 416 ρ0(x) = (1 if r≤1 2, 0.125 if r > 1 2,p0(x) = (1 if r≤1 2, 0.1 if r > 1 2,u0(x) = 0, (84) with r=px2+y2+z2 . Dirichlet boundary conditions are imposed and the domain is covered by 2280182 417 tetrahedra. 418 The solution obtained using the LADER-ENO scheme with CFL = 1 , cα= 3 , up to tend = 0.25 is 419 depicted in Figure 19. As reference solution we solve again the 1D code for Euler equations introduced in 420 Section 4.3 updated with appropriate source terms to account for three-dimensional eects. The agreement 421 observed for the 1D cuts of density, velocity magnitude and pressure prove the capability of the method to 422 handle three-dimensional problems. 423 22 Figure 16: Circular explosion. The left top image corresponds to the 3D plot of the obtained ρ at the nal time whereas the 1D plots containing the reference solution (black continuous line), a 1D cut (blue squared line) and the scatter plot (red dots), correspond to the ρ , p , and |u| elds obtained using the rst order scheme ( cα= 1 ). 23 Figure 17: Circular explosion. The left top image corresponds to the 3D plot of the obtained ρ at the nal time whereas the 1D plots containing the reference solution (black continuous line), a 1D cut (blue squared line) and the scatter plot (red dots), correspond to the ρ , p , and |u| elds obtained using the LADER-ENO scheme ( cα= 1 ). Figure 18: Circular explosion. Comparison between the numerical solution obtained with the rst order scheme (blue squares) and the LADER-ENO approximation (green circles). 24 Figure 19: LADER-ENO solution for the 3D spherical explosion test at nal time. Surface mesh on the boundary and density contours on the interior surfaces obtained after taking away the rst quadrant and 1D plots for density, velocity magnitude and pressure elds: 1D cut on x∈[0,1] , y =z= 0 (blue squares), scatter plot (red dots), reference solution (black line). 25 4.5. First problem of Stokes 424 To further analyse the behaviour of the developed method in the incompressible limit, we now consider 425 the rst problem of Stokes, [107]. The initial condition, dened in Ω=[−0.5,0.5] ×[−0.5,0.5] , reads 426 ρ0(x) = 1, p0(x) = 1 γ, u0 1(x)=0, u0 2(x) = −0.1 if y≤0, 0.1 if y > 0 (85) In the incompressible limit, this test case has an exact analytical solution for u2 given by 427 u2(x, t) = 1 10erf x 2√µt. (86) To complete the physical set up, we dene γ=cp= 1.4 , λ= 0 , leading to M≈10−1 . Regarding boundary 428 conditions, we impose the exact values for density and velocity in the x -direction, while on the top and 429 bottom boundaries, we set periodic boundary conditions in y -direction. Meanwhile, the exact values for 430 density and velocity are employed in the remaining boundaries. Finally, three dierent simulations are run 431 attending to the value for the viscosity coecient: µ= 10−2 , µ= 10−3 , and µ= 10−4 . The simulations are 432 run on a triangular primal mesh made of 1000 elements up to time tend = 1 . The vertical velocity along y= 0 433 is plotted in Figure 20 against the exact solution. We observe a good agreement between both curves for all 434 three viscosities. Let us note that µ= 10−4 is the only simulation run using the ENO reconstruction so that 435 we completely avoid the small bump arising after the discontinuity if any limiting strategy is employed. In 436 the other cases, such reconstruction can be neglected due to the high physical viscosity considered. 437 Figure 20: Comparison of the vertical velocity u2 along the 1D cut y= 0 obtained for the rst Stokes problem using LADER scheme against the exact solution at tend = 1 . µ= 10−4 LADER-EN (left), µ= 10−3 LADER without limiters (center), µ= 10−2 LADER without limiters (right). 4.6. Viscous shock 438 Here we analyse a steady viscous shock with Ms= 2 the shock Mach number. Considering the particular 439 case Pr = 0.75 , with Pr the Prandtl number, it is possible to nd an exact solution of the compressible Navier440 Stokes equations, derived by Becker in 1923, see [108, 109, 110] for all the necessary details to setup this 441 test case. 442 The computational domain Ω=[−0.5,0.5] ×[0,0.1] is discretized with 12500 triangular elements of 443 characteristic mesh spacing h= 1/250 . The shock wave is centered at x= 0 . 444 The values of the uid in front of the shock wave are given by ρ0= 1 , u0=−2 , v0=w0= 0 , and 445 p0= 1/γ so that the corresponding sound speed is c0= 1 and the uid is moving into the shock from the 446 right to the left at shock Mach number Ms= 2 . The Reynolds number based on a unitary reference length 447 ( L= 1 ) and on the ow speed u0 is given by Res=ρ0c0MsL µ . The uid parameters are chosen as γ= 1.4 , 448 cv= 2.5 , µ= 2 ·10−2 and λ= 91 3·10−2 , hence the corresponding shock Reynolds number is Res= 100 . 449 The simulation with the new hybrid FV/FE scheme proposed in this paper is run until time tend = 0.025 , 450 setting cα= 3 . The comparison between the numerical solution and the exact solution is shown in Figure 21 451 for the density ρ , the velocity u , and the pressure p . For all quantities, one can note a very good agreement. 452 32 by Spanish MCIU under project MTM2017-86459-R and by FEDER and Xunta de Galicia funds under the 551 ED431C 2017/60 project. SB was also funded by INdAM via a GNCS grant for young researchers and by 552 a UniTN starting grant of the University of Trento. SB, LR and MD are members of the GNCS group of 553 INdAM. 554 References 555 [1] S. Busto, J. L. Ferrín, E. F. Toro, M. E. Vázquez-Cendón, A projection hybrid high order nite volume/nite element 556 method for incompressible turbulent ows, J. Comput. Phys. 353 (2018) 169192. 557 [2] S. Klainermann, A. Majda, Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible 558 limit of compressible uid, Comm. Pure Appl. Math. 34 (1981) 481524. 559 [3] S. Klainermann, A. Majda, Compressible and incompressible uids, Communications on Pure and Applied Mathematics 560 35 (1982) 629651. 561 [4] F. Harlow, J. Welch, Numerical calculation of time-dependent viscous incompressible ow of uid with a free surface, 562 Phys. Fluids 8 (1965) 21822189. 563 [5] A. Chorin, A numerical method for solving incompressible viscous ow problems, J. Comput. Phys. 2 (1967) 1226. 564 [6] A. Chorin, Numerical solution of the NavierStokes equations, Math. Comput. 23 (1968) 341354. 565 [7] S. V. Patankar, D. B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional 566 parabolic ows, Int J Heat Mass Transfer 15 (10) (1972) 17871806. 567 [8] V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. 568 [9] J. van Kan, A second-order accurate pressure correction method for viscous incompressible ow, SIAM Journal on 569 Scientic and Statistical Computing 7 (1986) 870891. 570 [10] J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, J. 571 Comput. Phys. 85 (2) (1989) 257283. 572 [11] C. W. Hirt, B. D. Nichols, Volume of uid (VOF) method for dynamics of free boundaries, J. Comput. Phys. 39 (1981) 573 201225. 574 [12] V. Casulli, A semiimplicit numerical method for the freesurface NavierStokes equations, International Journal for 575 Numerical Methods in Fluids 74 (2014) 605622. 576 [13] C. Taylor, P. Hood, A numerical solution of the Navier-Stokes equations using the nite element technique, Computers 577 and Fluids 1 (1973) 73100. 578 [14] A. Brooks, T. Hughes, Stream-line upwind/Petrov Galerkin formulstion for convection dominated ows with particular 579 emphasis on the incompressible Navier-Stokes equation, Computer Methods in Applied Mechanics and Engineering 32 580 (1982) 199259. 581 [15] T. Hughes, M. Mallet, M. Mizukami, A new nite element formulation for computational uid dynamics: II. Beyond 582 SUPG, Computer Methods in Applied Mechanics and Engineering 54 (1986) 341355. 583 [16] M. Fortin, Old and new nite elements for incompressible ows, International Journal for Numerical Methods in Fluids 584 1 (1981) 347364. 585 [17] R. Verfürth, Finite element approximation of incompressible Navier-Stokes equations with slip boundary condition II, 586 Numerische Mathematik 59 (1991) 615636. 587 [18] J. G. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier-Stokes Problem. I. Regularity 588 of solutions and second order error estimates for spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275311. 589 [19] J. G. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier-Stokes Problem. III. Smoothing 590 property and higher order error estimates for spatial discretization, SIAM J. Numer. Anal. 25 (1988) 489512. 591 [20] P. Lax, B. Wendro, Systems of conservation laws, Commun. Pur. Appl. Math. 13 (2) (1960) 217237. 592 [21] S. K. Godunov, A nite dierence method for the computation of discontinuous solutions of the equations of uid 593 dynamics, Mat. Sb. 47 (1959) 357393. 594 [22] P. Roe, Approximate Riemann solvers, parameter vectors, and dierence schemes, J. Comput. Phys. 43 (1981) 357372. 595 [23] S. Osher, F. Solomon, Upwind dierence schemes for hyperbolic conservation laws, Math. Comput. 38 (1982) 339374. 596 [24] A. Harten, P. Lax, B. van Leer, On upstream dierencing and Godunov-type schemes for hyperbolic conservation laws, 597 Vol. 25, 1983, pp. 3561. 598 [25] B. Einfeldt, C. Munz, P. Roe, B. Sjögreen, On Godunov-type methods near low densities, J. Comput. Phys. 92 (1991) 599 273295. 600 [26] C. D. Munz, On Godunovtype schemes for Lagrangian gas dynamics, SIAM Journal on Numerical Analysis 31 (1994) 601 1742. 602 [27] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the Harten-Lax-van Leer Riemann solver, 603 Journal of Shock Waves 4 (1994) 2534. 604 [28] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, 2002. 605 [29] E. F. Toro, Riemann solvers and numerical methods for uid dynamics: A practical introduction, Springer, 2009. 606 [30] V. Casulli, D. Greenspan, Pressure method for the numerical solution of transient, compressible uid ows, Int. J. Numer. 607 Methods Fluids 4 (1984) 10011012. 608 [31] A. Meister, Asymptotic single and multiple scale expansions in the low mach number limit, SIAM Journal on Applied 609 Mathematics 60 (1) (1999) 256271. 610 [32] C. Munz, R. Klein, S. Roller, K. Geratz, The extension of incompressible ow solvers to the weakly compressible regime, 611 Computers and Fluids 32 (2003) 173196. 612 33 [33] R. Klein, Semi-implicit extension of a godunov-type scheme based on low mach number asymptotics I: one-dimensional 613 ow, J. Comput. Phys. 121 (1995) 213237. 614 [34] R. Klein, N. Botta, T. Schneider, C. Munz, S.Roller, A. Meister, L. Homann, T. Sonar, Asymptotic adaptive methods 615 for multi-scale problems in uid mechanics, Journal of Engineering Mathematics 39 (2001) 261343. 616 [35] J. Park, C. Munz, Multiple pressure variables methods for uid ow at all Mach numbers, International journal for 617 numerical methods in uids 49 (8) (2005) 905931. 618 [36] E. Toro, M. Vázquez-Cendón, Flux splitting schemes for the Euler equations, Computers and Fluids 70 (2012) 112. 619 [37] F. Cordier, P. Degond, A. Kumbaro, An Asymptotic-Preserving all-speed scheme for the Euler and Navier-Stokes equa620 tions, J. Comput. Phys. 231 (2012) 56855704. 621 [38] P. Degond, M. Tang, All speed scheme for the low Mach number limit of the isentropic Euler equations, Comm. Comput. 622 Phys. 10 (1) (2011) 131. 623 [39] M. Dumbser, V. Casulli, A conservative, weakly nonlinear semi-implicit nite volume scheme for the compressible Navier624 Stokes equations with general equation of state, Applied Mathematics and Computation 272 (2016) 479497. 625 [40] S. Boscarino, G. Russo, L. Scandurra, All Mach number second order semi-implicit scheme for the Euler equations of 626 gasdynamics, J. Sci. Comput. 77 (2018) 850884. 627 [41] G. Dimarco, R. Loubère, V. Michel-Dansac, M. Vignal, Second-order implicit-explicit total variation diminishing schemes 628 for the euler system in the low mach regime, J. Comput. Phys. 372 (2018) 178  201. 629 [42] E. Abbate, A. Iollo, G. Puppo, An asymptotic-preserving all-speed scheme for uid dynamics and nonlinear elasticity, 630 SIAM Journal on Scientic Computing 41 (2019) A2850A2879. 631 [43] S. Avgerinos, F. Bernard, A. Iollo, G. Russo, Linearly implicit all mach number shock capturing schemes for the euler 632 equations, J. Comput. Phys. 393 (2019) 278  312. 633 [44] W. Boscheri, G. Dimarco, R. Loubère, M. Tavelli, M. Vignal, A second order all mach number imex nite volume solver 634 for the three dimensional euler equations, Journal of Computational Physics 415 (2020) 109486. 635 [45] W. Boscheri, G. Dimarco, M. Tavelli, An ecient second order all Mach nite volume solver for the compressible 636 NavierStokes equations, Computer Methods in Applied Mechanics and Engineering 374 (2021) 113602. 637 [46] W. Boscheri, L. Pareschi, High order pressure-based semi-implicit IMEX schemes for the 3D Navier-Stokes equations at 638 all Mach numbers, Journal of Computational PhysicsTo appear (2021). 639 [47] S. Shanmuganathan, D. L. Youngs, J. Griond, B. Thornber, R. Williams, Accuracy of high-order density-based com640 pressible methods in low mach vortical ows, International Journal for Numerical Methods in Fluids 74 (2014) 335358. 641 [48] N. Fleischmann, S. Adami, X. Hu, N. Adams, A low dissipation method to cure the grid-aligned shock instability, Journal 642 of Computational Physics 401 (2020) 109004. 643 [49] F. Bassi, S. Rebay, A high-order accurate discontinuous nite element method for the numerical solution of the com644 pressible Navier-Stokes equations, J. Comput. Phys. 131 (1997) 267279. 645 [50] C. Baumann, J. Oden, A discontinuous hp nite element method for convection-diusion problems, Comput. Methods 646 Appl. Mech. Eng. 175 (3-4) (1999) 311341. 647 [51] C. Baumann, J. Oden, A discontinuous hp nite element method for the euler and navier-stokes equations, Int. J. Numer. 648 Methods Fluids 31 (1) (1999) 7995. 649 [52] B. Cockburn, C. W. Shu, The local discontinuous Galerkin method for time-dependent convection diusion systems, 650 SIAM Journal on Numerical Analysis 35 (1998) 24402463. 651 [53] B. Cockburn, C. W. Shu, Runge-Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. 652 Comput. 16 (2001) 199224. 653 [54] F. Bassi, A. Crivellini, D. D. Pietro, S. Rebay, An implicit high-order discontinuous Galerkin method for steady and 654 unsteady incompressible ows, Computers and Fluids 36 (2007) 15291546. 655 [55] E. Ferrer, R. Willden, A high order discontinuous galerkin nite element solver for the incompressible navierstokes 656 equations, Computer and Fluids 46 (2011) 224230. 657 [56] N. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable discontinuous galerkin method for the incom658 pressible navier-stokes equations, J. Comput. Phys. 230 (2011) 11471170. 659 [57] S. Rhebergen, B. Cockburn, A space-time hybridizable discontinuous Galerkin method for incompressible ows on de660 forming domains, J. Comput. Phys. 231 (2012) 41854204. 661 [58] S. Rhebergen, B. Cockburn, J. J. van der Vegt, A space-time discontinuous Galerkin method for the incompressible 662 Navier-Stokes equations, J. Comput. Phys. 233 (2013) 339358. 663 [59] A. Crivellini, V. D'Alessandro, F. Bassi, High-order discontinuous Galerkin solutions of three-dimensional incompressible 664 RANS equations, Computers and Fluids 81 (2013) 122133. 665 [60] B. Klein, F. Kummer, M. Oberlack, A SIMPLE based discontinuous Galerkin solver for steady incompressible ows, J. 666 Comput. Phys. 237 (2013) 235250. 667 [61] V. Dolejsi, Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible ows, Comm. Comput. 668 Phys. 4 (2008) 231274. 669 [62] V. Dolejsi, M. Feistauer, A semi-implicit discontinuous Galerkin nite element method for the numerical solution of 670 inviscid compressible ow, J. Comput. Phys. 198 (2004) 727746. 671 [63] V. Dolejsi, M. Feistauer, J. Hozman, Analysis of semi-implicit DGFEM for nonlinear convection-diusion problems on 672 nonconforming meshes, Comput. Methods Appl. Mech. Eng. 196 (2007) 28132827. 673 [64] M. Tavelli, M. Dumbser, A staggered space-time discontinuous Galerkin method for the incompressible Navier-Stokes 674 equations on two-dimensional triangular meshes, Comput. Fluids 119 (2015) 235  249. 675 [65] M. Tavelli, M. Dumbser, A staggered space-time discontinuous Galerkin method for the three-dimensional incompressible 676 Navier-Stokes equations on unstructured tetrahedral meshes, J. Comput. Phys. 319 (2016) 294  323. 677 34 [66] M. Tavelli, M. Dumbser, A pressure-based semi-implicit space-time discontinuous Galerkin method on staggered unstruc678 tured meshes for the solution of the compressible Navier-Stokes equations at all Mach numbers, J. Comput. Phys. 341 679 (2017) 341  376. 680 [67] F. Fambri, M. Dumbser, Spectral semi-implicit and space-time discontinuous Galerkin methods for the incompressible 681 Navier-Stokes equations on staggered Cartesian grids, Applied Numerical Mathematics 110 (2016) 4174. 682 [68] F. Fambri, M. Dumbser, Semi-implicit discontinuous Galerkin methods for the incompressible Navier-Stokes equations 683 on adaptive staggered Cartesian grids, Computer Methods in Applied Mechanics and Engineering 324 (2017) 170203. 684 [69] S. Busto, M. Tavelli, W. Boscheri, M. Dumbser, Ecient high order accurate staggered semi-implicit discontinuous 685 galerkin methods for natural convection problems, Comput. Fluids 198 (2020) 104399. 686 [70] I. Peshkov, E. Romenski, A hyperbolic model for viscous Newtonian ows, Continuum Mechanics and Thermodynamics 687 28 (2016) 85104. 688 [71] M. Dumbser, I. Peshkov, E. Romenski, O. Zanotti, High order ADER schemes for a unied rst order hyperbolic 689 formulation of continuum mechanics: Viscous heat-conducting uids and elastic solids, J. Comput. Phys. 314 (2016) 824 690  862. 691 [72] S. Busto, S. Chiocchetti, M. Dumbser, E. Gaburro, I. Peshkov, High order ADER schemes for continuum mechanics, 692 Frontiers in Physiccs 8 (2020) 32. doi:10.3389/fphy.2020.00032. 693 [73] W. Boscheri, M. Dumbser, M. Ioriatti, I. Peshkov, E. Romenski, A structure-preserving staggered semi-implicit nite 694 volume scheme for continuum mechanics, Journal of Computational Physics 2021 (2010) 109866. 695 [74] I. Peshkov, M. Dumbser, W. Boscheri, E. Romenski, S. Chiocchetti, M. Ioriatti, Modeling solid-uid transformation in 696 non-newtonian viscoplastic ows with a unied ow theory, Computers and FluidsSubmitted. 697 [75] S. Godunov, An interesting class of quasilinear systems, Dokl. Akad. Nauk SSSR 139(3) (1961) 521523. 698 [76] S. Godunov, E. Romenski, Nonstationary equations of the nonlinear theory of elasticity in Euler coordinates., Journal of 699 Applied Mechanics and Technical Physics 13 (1972) 868885. 700 [77] S. Godunov, Symmetric form of the magnetohydrodynamic equation, Numerical Methods for Mechanics of Continuum 701 Medium 3 (1) (1972) 2634. 702 [78] E. Romenski, Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechanics, Mathe703 matical and computer modelling 28(10) (1998) 115130. 704 [79] S. Godunov, E. Romenski, Elements of Continuum Mechanics and Conservation Laws, Kluwer Academic/ Plenum Pub705 lishers, 2003. 706 [80] L. Li, J. Luo, H. Nishikawa, H. Luo, Reconstructed discontinuous Galerkin methods for compressible ows based on a 707 new hyperbolic navier-stokes system, Journal of Computational Physics (2021) 110058. 708 [81] A. Bermúdez, J. L. Ferrín, L. Saavedra, M. E. Vázquez-Cendón, A projection hybrid nite volume/element method for 709 low-Mach number ows, J. Comp. Phys. 271 (2014) 360378. 710 [82] S. Busto, G. Stabile, G. Rozza, M. Vázquez-Cendón, POD-Galerkin reduced order methods for combined Navier-Stokes 711 transport equations based on a hybrid FV-FE solver, Computers & Mathematics with Applications 79 (2) (2020) 256  712 273. 713 [83] A. Bermúdez, S. Busto, M. Dumbser, J. Ferrín, L. Saavedra, M. Vázquez-Cendón, A staggered semi-implicit hybrid fv/fe 714 projection method for weakly compressible ows, Journal of Computational Physics 421 (2020) 109743. 715 [84] S. Busto, E. F. Toro, M. E. Vázquez-Cendón, Design and analisis of ADERtype schemes for model advectiondiusion 716 reaction equations, J. Comp. Phys. 327 (2016) 553575. 717 [85] A. Bermúdez, S. Busto, J. L. Ferrín, L. Saavedra, E. F. Toro, M. E. Vázquez-Cendón, SEMA SIMAI Springer Series. 718 Computational Mathematics, Numerical Analysis and Applications, Springer, 2017, Ch. A projection hybrid nite volume719 ADER/nite element method for turbulent Navier-Stokes, pp. 201206. 720 [86] A. Bermúdez, S. Busto, M. Dumbser, F. Ferrín, V.-C. M. Saavedra, L., A staggered semi-implicit hybrid FV/FE projection 721 method for weakly compressible ows, J. Comput. Phys. 421 (2020) 109743. 722 [87] L. Brugnano, V. Casulli, Iterative solution of piecewise linear systems, SIAM Journal on Scientic Computing 30 (2007) 723 463472. 724 [88] L. Brugnano, V. Casulli, Iterative solution of piecewise linear systems and applications to ows in porous media, SIAM 725 Journal on Scientic Computing 31 (2009) 18581873. 726 [89] L. Brugnano, A. Sestini, Iterative solution of piecewise linear systems for the numerical solution of obstacle problems, 727 Journal of Numerical Analysis, Industrial and Applied Mathematics 6 (2012) 6782. 728 [90] V. Casulli, P. Zanolli, A nested Newtontype algorithm for nite volume methods solving Richards' equation in mixed 729 form, SIAM Journal on Scientic Computing 32 (2009) 22552273. 730 [91] V. Casulli, P. Zanolli, Iterative solutions of mildly nonlinear systems, Journal of Computational and Applied Mathematics 731 236 (2012) 39373947. 732 [92] V. V. Rusanov, The calculation of the interaction of non-stationary shock waves and obstacles, USSR Computational 733 Mathematics and Mathematical Physics 1 (1962) 304320. 734 [93] S. Busto, Contributions to the numerical solution of heterogeneous uid mechanics models, Ph.D. thesis, Universidade 735 de Santiago de Compostela (2018). 736 [94] E. F. Toro, R. C. Millington, L. A. M. Nejad, Godunov methods, Springer, 2001, Ch. Towards very high order Godunov 737 schemes. 738 [95] R. Millington, E. Toro, L. Nejad, Arbitrary high order methods for conservation laws i: The one dimensional scalar case, 739 Ph.D. thesis, Manchester Metropolitan University, Department of Computing and Mathematics (June 1999). 740 [96] M. Dumbser, C. Enaux, E. F. Toro, Finite volume schemes of very high order of accuracy for sti hyperbolic balance 741 laws, J. Comput. Phys. 227 (8) (2008) 3971  4001. 742 35 [97] W. Boscheri, M. Dumbser, A direct arbitrary-lagrangian-eulerian ADER-WENO nite volume scheme on unstructured 743 tetrahedral meshes for conservative and non-conservative hyperbolic systems in 3D, J. Comput. Phys. 275 (2014) 484523. 744 [98] T. Barth, D. Jespersen, The design and application of upwind schemes on unstructured meshes, Tech. rep. (1989). 745 [99] P. L. Roe, Modelling of Discontinuous Flows, Vol. 22, 1985. 746 [100] S. Clain, S. Diot, R. Loubère, A high-order nite volume method for systems of conservation lawsmulti-dimensional 747 optimal order detection (mood), J. Comput. Phys. 230 (2011) 40284050. 748 [101] V. Casulli, P. Zanolli, A nested newton-type algorithm for nite volume methods solving richards' equation in mixed 749 form, SIAM J. Sci. Comput. 32 (4) (2010) 22552273. 750 [102] M. Dumbser, D. S. Balsara, E. F. Toro, C.-D. Munz, A unied framework for the construction of one-step nite volume 751 and discontinuous Galerkin schemes on unstructured meshes, J. Comput. Phys. 227 (18) (2008) 82098253. 752 [103] L. Pareschi, G. Russo, Implicit-explicit Runge-Kutta schemes for sti systems of dierential equations, Advances in the 753 Theory of Computational Mathematics 3 (2000) 269288. 754 [104] G. A. Sod, A survey of several nite dierence methods for systems of nonlinear hyperbolic conservation laws, J. Comput. 755 Phys. 27 (1) (1978) 1  31. 756 [105] P. Woodward, P. Colella, The numerical simulation of two-dimensional uid ow with strong shocks, Journal of Compu757 tational Physics 54 (1984) 115173. 758 [106] V. A. Titarev, E. F. Toro, ADER schemes for three-dimensional non-linear hyperbolic systems, J. Comp. Phys. 204 (2) 759 (2005) 715736. 760 [107] H. Schlichting, K. Gersten, Boundary-layer theory, Springer, 2016. 761 [108] R. Becker, Stosswelle und Detonation, Physik 8 (1923) 321. 762 [109] A. Bonnet, J. Luneau, Aérodynamique. Théories de la dynamique des uides, Cepadues Editions, Toulouse, 1989, iSBN: 763 2.85428.218.3. 764 [110] M. Dumbser, I. Peshkov, E. Romenski, O. Zanotti, High order ADER schemes for a unied rst order hyperbolic 765 formulation of continuum mechanics: Viscous heat-conducting uids and elastic solids, Journal of Computational Physics 766 314 (2016) 824862. 767 [111] U. Ghia, K. Ghia, C. Shin, High-re solutions for incompressible ow using the Navier-Stokes equations and a multigrid 768 method, J. Comput. Phys. 48 (3) (1982) 387  411. 769 [112] M. van Dyke, An album of uid motion, The Parabolic Press, 2005. 770 [113] M. Dumbser, M. Käser, V. A. Titarev, E. F. Toro, Quadrature-free non-oscillatory nite volume schemes on unstructured 771 meshes for nonlinear hyperbolic systems, Journal of Computational Physics 226 (2007) 204243. 772 [114] F. Kemm, E. Gaburro, F. Thein, M. Dumbser, A simple diuse interface approach for compressible ows around moving 773 solids of arbitrary shape based on a reduced Baer-Nunziato model, Computers and Fluids 204 (2020) 104536. 774 [115] H. Schardin, in: Proc. VII Int. Cong. High Speed Photg., Darmstadt, O. Helwich Verlag, 1965, pp. 113119. 775 [116] V. Casulli, R. T. Cheng, Semi-implicit nite dierence methods for threedimensional shallow water ow, International 776 Journal for Numerical Methods in Fluids 15 (1992) 629648. 777 [117] V. Casulli, R. A. Walters, An unstructured grid, threedimensional model based on the shallow water equations, Inter778 national Journal for Numerical Methods in Fluids 32 (2000) 331348. 779 [118] S. C. Kramer, G. S. Stelling, A conservative unstructured scheme for rapidly varied ows, International Journal for 780 Numerical Methods in Fluids 58 (2008) 183212. 781 [119] M. Tavelli, M. Dumbser, A high order semi-implicit discontinuous Galerkin method for the two dimensional shallow water 782 equations on staggered unstructured meshes, Appl. Math. Comput. 234 (2014) 623644. 783 [120] P. K.G., Upwind and High-Resolution Schemes, Springer, 1997, Ch. An Approximate Riemann Solver for Magnetohy784 drodynamics. 785 [121] C. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, U. Voss, Divergence correction techniques for Maxwell solvers based 786 on a hyperbolic model, Journal of Computational Physics 161 (2000) 484511. 787 [122] D. Balsara, Second order accurate schemes for magnetohydrodynamics with divergence-free reconstruction, The Astro788 physical Journal Supplement Series 151 (08 2003). doi:10.1086/381377. 789 [123] D. Balsara, M. Dumbser, R. Abgrall, Multidimensional HLL and HLLC Riemann solvers for unstructured meshes-with 790 application to Euler and MHD ows, Journal of Computational Physics 261 (2014) 172208. 791 [124] M. Dumbser, D. Balsara, M. Tavelli, F. Fambri, A divergence-free semi-implicit nite volume scheme for ideal, viscous, 792 and resistive magnetohydrodynamics, International Journal for Numerical Methods in Fluids 89 (1-2) (2019) 1642. 793