A weakly-compressible Cartesian grid approach for hydrodynamic flows

Accepted Manuscript A weakly-compressible Cartesian grid approach for hydrodynamic flows P. Bigay, G. Oger, P.-M. Guilcher, D. Le Touzé

PII: DOI: Reference:

S0010-4655(17)30192-3 http://dx.doi.org/10.1016/j.cpc.2017.06.010 COMPHY 6244

To appear in:

Computer Physics Communications

Received date : 16 March 2017 Revised date : 9 June 2017 Accepted date : 12 June 2017 Please cite this article as: P. Bigay, G. Oger, P.-. Guilcher, D. Le Touzé, A weakly-compressible Cartesian grid approach for hydrodynamic flows, Computer Physics Communications (2017), http://dx.doi.org/10.1016/j.cpc.2017.06.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Computer Physics Communications Computer Physics Communications 00 (2017) 1–22

A weakly-compressible Cartesian grid approach for hydrodynamic flows P. Bigaya,b , G. Ogera,∗, P.-M. Guilcherb , D. Le Touz´ea a Ecole

Centrale Nantes, LHEEA research dept. (ECN and CNRS), Nantes, France b Nextflow Software, Nantes, France

Abstract The present article aims at proposing an original strategy to solve hydrodynamic flows. In introduction, the motivations for this strategy are developed. It aims at modeling viscous and turbulent flows including complex moving geometries, while avoiding meshing constraints. The proposed approach relies on a weakly-compressible formulation of the Navier-Stokes equations. Unlike most hydrodynamic CFD (Computational Fluid Dynamics) solvers usually based on implicit incompressible formulations, a fully-explicit temporal scheme is used. A purely Cartesian grid is adopted for numerical accuracy and algorithmic simplicity purposes. This characteristic allows an easy use of Adaptive Mesh Refinement (AMR) methods embedded within a massively parallel framework. Geometries are automatically immersed within the Cartesian grid with an AMR compatible treatment. The method proposed uses an Immersed Boundary Method (IBM) adapted to the weakly-compressible formalism and imposed smoothly through a regularization function, which stands as another originality of this work. All these features have been implemented within an in-house solver based on this WCCH (Weakly-Compressible Cartesian Hydrodynamic) method which meets the above requirements whilst allowing the use of high-order (> 3) spatial schemes rarely used in existing hydrodynamic solvers. The details of this WCCH method are presented and validated in this article. Keywords: Weakly-compressible formulation, Cartesian grid, Fully-explicit scheme, WENO5, Local Mesh Refinement, Numerical diffusion, High-order schemes, Immersed-Boundary Method (IBM), Mobile geometries

1. Motivations When solving free-surface hydrodynamic flows of marine engineering interest, i.e. flows around large bodies in the sea, people usually use an incompressible flow assumption (div ~u = 0), implying a null M ach number assumption (pressure waves have an infinite propagation speed). This particularity entails that most of state-of-the art hydrodynamic solvers use an implicit solution of the Navier-Stokes equations. In the context of High Performance Computing, good parallel performances can be obtained with implicit solvers but with possibly consequent implementation efforts. Furthermore, in a context where more and more unsteady problems need to be solved, the time-steps used for marching the solution in time need to be smaller and smaller, increasing the pertinence a fully explicit resolution. For free-surface problems where Volume-Of-Fluid or Level-Set techniques are used, the situation is accentuated since explicit advection is used for the volume fraction/distance function. As a result, the time-steps used in these implicit solvers get ∗ Corresponding

author Email address: [email protected] (G. Oger)

1

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

2

closer and closer to the ones that would be imposed by the CFL (Courant-Friedrichs-Lewy) condition of an explicit resolution for the same problem. Additionally, hydrodynamic solvers require an accurate treatment of fixed and moving boundaries. Most of the usual CFD methods described in the literature use complex meshing strategies to set-up: unstructured meshes, structured meshes with body-fitted meshes (with all the complexity that arises when dealing with boundary layers). To treat mobile geometries, various methods may be used: sliding grids, grid morphing, overlapping grids (with the use of interpolations), re-meshing strategies. These methods display some advantages such as an enhanced accuracy in near body regions, but they also yield some drawbacks such as interpolation issues, high computational costs related to mesh modifications, or human resource costs related to meshing operations. Furthermore, mainly due to this mesh complexity, these state-of-the-art hydrodynamic solvers rely on second-order accurate schemes (at least in space). This rather low order prevents solving efficiently problems such as, e.g., turbulent wake interactions. More generally, in the trend of modeling high-Reynolds number flows by solving a larger and larger part of the turbulent spectrum, the numerical diffusion of the scheme should be as low as possible, in order to correctly capture turbulent effects without bringing additional numerical viscosity. For instance when using Large-Eddy Simulation (LES) with a low-order scheme, the intrinsic numerical viscosity may become higher than the sub-grid viscosity added to the scheme to model the viscous dissipation at small scales. This analysis led us to propose in this paper a Weakly-Compressible Cartesian Hydrodynamic (WCCH) solver, relying on ingredients which differ from the standard (implicit) Navier-Stokes solvers usually adopted in naval hydrodynamics. We thus propose to combine: • a weakly-compressible form of the Navier-Stokes equations, • an explicit solution, • a purely Cartesian mesh with local mesh refinement capabilities, • high-order (> 3) space-time schemes, • the use of an Immersed Boundary Method (IBM) to describe the presence of solid bodies. These different technical choices are driven by the will to find a good compromise between accuracy, performance, adaptability and algorithmic simplicity to answer the above points, in the context of more and more intensive calculations involving larger and larger meshes. First of all, the WCCH solver is based on the Finite Volume approach, for its strict conservation properties (mass, momentum and energy). A Cartesian grid is used exclusively, allowing for simplified algorithmic and a fast evaluation of the various spatial derivative operators. In addition, Cartesian grids are naturally well adapted to AMR methods [5, 18]. The choice of local mesh refinement has been made in order to reach high accuracy and to capture very local phenomena, while minimizing the total number of cells and concentrating fine spatial resolutions in the most important flow regions, while making it possible to preserve a high-order scheme all over the domain. With this meshing and explicit solution strategy, it is possible to use a high-order scheme in a straightforward way. A 5th order WENO (Weighted Essentially Non-Oscillatory) scheme is used in the present work. The proposed WCCH method strategy is also based on a fully-explicit solution, thanks to a weaklycompressible formulation of the Navier-Stokes equations. Contrary to an implicit solution, only variables from the current instant tn are necessary for their time updating at tn+1 . In this explicit solving context, computational times are efficiently optimized by using parallel computing on distributed memory architectures (MPI) involving an arbitrary high (thousand) number of computational cores. In return, the main limitation of the explicit schemes lies on the CFL condition for slow or viscosity-driven flows, imposing much smaller time steps than those of implicit schemes. This restriction is mitigated by the use of a weaklycompressible approach allowing to close the system through the use of an equation of state linking pressure to density, but not using the true sound speed of the considered medium, and still allowing a fully-explicit solution of the problem. Not having to enforce the incompressible feature of the flow prevents from solving a Poisson equation, and the subsequent linear system. A similar approach is used for instance in the 2

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

3

weakly-compressible SPH method [25, 27]. In addition to conferring its explicit feature to the solver, actual compressible flows can also be treated if required. Concerning the treatment of complex geometries, the strategy adopted here resides in the use of a fixed Cartesian grid without any modification (apart from local static refinement), even in the case of mobile boundaries. This method is based on the Immersed Boundary Method (IBM), for which the boundary Features the terms method: conditionsMain are enforced by of source acting on each Cartesian cell inside the considered boundary. The treatment of moving boundaries inside a Cartesian grid is especially addressed in this article. Figure 1 summarizes the different characteristics of the proposed WCCH solver. Finite Volume method Good conservation properties

Immersed Boundary Method

Cartesian grid Precision Performance HPC

High-Order methods

Automatic geometry handling Adaptive Mesh Refinement (AMR)

Weakly-compressible solution (Ma<0.1) Explicit resolution, Intrinsic multiphysic capabilities

Figure 1: Scheme of the main characteristics of the WCCH solver. In the present article, the WCCH solver is firstly introduced by discussing its weakly-compressible approach based with a preliminary 2nd -order method using the MUSCL (Monotonic Upwind Scheme for Conservation Laws) scheme. The intrinsic numerical diffusion is addressed and evaluated on different test cases, involving steady and convected vortices. To reduce the inherent numerical diffusion, a higher-order method based on a 5th -order WENO scheme is presented and discussed. The second part of this study focuses on validations onto laminar viscous flows. Lastly, the Direct Forcing Immersed Boundary Method applied for imposing boundary conditions on complex geometries is presented and discussed. This method is finally validated on low-Reynolds test cases involving fixed and moving geometries. 2. Hyperbolic solver The problem to solve resides on the Navier-Stokes equations ∂ρ ~ + ∇.(ρ~u) = ∂t ∂ρ~u ~ ~ + ∇.(ρ~u ⊗ ~u) + ∇P ∂t

=

0,

(1)

~ τ¯ + ρf~ , ∇.

(2)

where f~, ρ, P , ~u and τ¯ are respectively the body forces, density, pressure, velocity and viscous stress tensor. The strategy adopted here consists in decomposing the problem into an hyperbolic part (Euler equations) to which a viscous part resolution is then added if needed. 2.1. Euler equations in weakly compressible approach For a compressible inviscid flow, the following Euler equations are addressed: ∂ρ ~ + ∇.(ρ~u) ∂t

=

0,

(3)

∂ρ~u ~ ~ + ∇.(ρ~u ⊗ ~u) + ∇P ∂t 3

=

0.

(4)

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

4

In the present work this system is thermodynamically closed with the following barotropic equation of state:  γ  ρ ρ0 c20 −1 , (5) P = γ ρ0 with γ, ρ0 and c0 the polytropic constant, nominal density and nominal speed of sound respectively. Values commonly used for γ are 1.4 for air and 7 for water. This equation of state does not involve the fluid energy, which has not to be solved. In the field of naval hydrodynamics, flows are usually assumed to be incompressible (i.e. div(~v ) = 0), which corresponds to a null M ach number so that pressure waves propagate at infinite speed. Instead, the solver presented hereafter is based on the compressible Navier-Stokes equations, leading to a weaklycompressible approach. Such an approach is used for instance in the SPH method [25, 27]. It allows a fully explicit resolution in time with the following main advantages: • to avoid the resolution of a linear system (needed in an incompressible approach), • to ease and simplify significantly the implementation of the scheme in a parallel framework together with high parallel performances (see [28] for more details). The sound speed c0 can be chosen arbitrarily provided that the density variations remains small enough, that is for instance: |ρ − ρ0 | ≤ 0.01, (6) ρ0 which corresponds to the following constraint on pressure variations using the equation of state (5): γP ≤ 1.01γ − 1. (7) ρ0 c20 As outlined by Marrone et al. [19], the choice of the sound speed is ruled by three predominant pressure scales: the static pressure scale due to the gravity ρgHmax (where Hmax is the highest liquid height value 2 encountered in the fluid domain at any instant), the pressure scale related to the kinetic energy ρUmax (where Umax is the maximum velocity value encountered in the fluid domain at any instant), and the acoustic pressure scale ρcUmax due to compressibility effects occurring namely during impact events. The respect of the condition (7) for each of these three pressure scales leads to the following constraints:

c0 c0 c0

r

γgHmax , 1.01γ − 1 r γ ≥ Umax , 1.01γ − 1 γUmax ≥ , 1.01γ − 1

≥

(8) (9) (10)

for which setting for instance γ = 1 leads exactly to the conditions presented in [19]. The last condition (10) is therefore the more restrictive one. Note that in the presence of impacts at large velocity the following choice might be finally preferred:   γUmax ϕ c0 = min , c , (11) 1.01γ − 1 0 where cϕ 0 is the physical sound speed of the considered fluid. Indeed, intense impacts can be responsible for physical compressibility effects requiring the use of the actual sound speed. In such cases the compressible approach adopted in this work is also particularity advantageous. The present methodology therefore requires an a priori estimation of the maximum velocity and height values Umax and Hmax in the flow. 4

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

5

2.2. Finite-Volume discretization and flux computation Euler equations can be expressed in the following conservative way:

with the conservative variables vector 

 ~ c (W ) =  and ∇.Ψ  |

  ρu  ρu2 + P   +   ρuv ρuw ,x {z } | F (W ),x

∂W ~ c (W ) = 0, + ∇.Ψ ∂t   ρ  ρu   W =  ρv , ρw   ρv ρw  ρwu ρvu   +  ρwv ρv 2 + P  ρvw ρw2 + P ,y {z } | {z G(W ),y

H(W ),z

(12)

   

,z

}

where F (W ), G(W ) and H(W ) are the flux vectors in direction x, y and z respectively. Equation (12) can be expressed in the integral form on a control volume V (t): ZZZ ZZZ ∂W ~ c (W ) dV = 0. ∇.Ψ dV + V (t) V (t) ∂t

(13)

The solver is purely Eulerian, therefore with a fixed control volume V . The divergence theorem leads to the following expression: ZZZ ZZ d W dV + Ψc (W ).~n dS = 0, (14) dt V S with ~n the unit normal vector to surface S. The above equation is classically discretized as follows: d Wijk 1 1 1 1 1 − H = − Fi+ 12 ,j,k ) + − Gi,j+ 12 ,k ) + (F 1 (G (H i,j,k+ 21 ), dt ∆x i− 2 ,j,k ∆y i,j− 2 ,k ∆z i,j,k− 2

(15)

with Wijk =

1 1 1 ∆x ∆y ∆z

ZZZ

W (x, y, z, t) dx dy dz,

(16)

Cijk

the averaged value in the conservative variables vector in cell Cijk = [xi− 21 , xi+ 12 ] × [yj− 21 , yj+ 12 ] × [zk− 12 , zk+ 12 ]. | {z } | {z } | {z } ∆x

∆y

∆z

The flux computation can then be performed by using any exact or approached Riemann solver (for instance HLL, HLLC...), and this flux is then computed from the left and right states at each cell edge [32]. These states are reconstructed from the surrounding cell values, with a given reconstruction spatial order which rules the final order of the global scheme. Two different types of reconstructions are used in WCCH: the first one is Total Variation Diminishing (MUSCL scheme [37]), while the second one is Essentially Non Oscillatory (5th order WENO scheme [17], implemented following Titarev & Toro [31]). dW The temporal derivative dtijk is finally integrated using a 4th order Runge-Kutta scheme. The explicit nature of the scheme implies the following CFL condition: ∆tacous = CF Lacous with the Courant number set to CF Lacous = 0.7.

5

mini (∆xi ) , c0

(17)

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

6

2.3. Numerical diffusion evaluation: the 2D Taylor-Green vortex We propose here to estimate the numerical diffusion of the scheme on a multi-dimensional academic case: the Taylor-Green vortex. The analytical solution for the velocity and pressure fields in a bi-dimensional periodic square domain of size L is: 

−8π 2 νt y x  2   u(x, y, t) = U0 e L 2 sin 2π L cos 2π L ,

−8π νt x (18) v(x, y, t) = −U0 e L2 cos 2π L sin 2π Ly ,  2 


 ρ0 U02 −16π2 νt y x cos 4π L + cos 4π L , p(x, y, t) = 4 e L

with ν the kinematic viscosity, U0 the reference flow velocity and ρ0 the nominal flow density. From a theoretical point of view, steady vortices are expected in the inviscid case (ν = 0). From a numerical point of view, a decay of the vortex intensity is expected due to the intrinsic diffusion of the scheme. The simulations presented hereafter have been performed with a nominal sound speed c0 = 10m/s, and with ρ0 = 1000kg/m3 , U0 = 1m/s and L = 1m. Both MUSCL and 5th order WENO schemes were tested. For the MUSCL scheme, two different limiters have been tested: Minmod [30] and Van Leer [36] limiters. Figure 2 compares the pressure fields obtained at time t = 3s with a 100x100 grid resolution in the MUSCL cases, and a 64x64 grid resolution in the WENO 5 case.

Figure 2: 2D Taylor-Green vortex test case. Pressure fields at time t= 3s with MUSCL scheme using Minmod (left) and Van Leer (middle) limiters with a 100x100 grid resolution, and for the WENO 5 scheme with a 64x64 grid resolution (right). As expected the Minmod limiter appears very diffusive compared to the Van Leer limiter, for which the pressure fields are better preserved. The lesser diffusive solution is obtained with the WENO 5 scheme, for which an undissipated solution is already obtained with a coarser spatial resolution. For a more quantitative evaluation, an exponential regression from the local velocity decay time evolution allowed to determine an equivalent numerical viscosity, by assuming ν = νnum in Equation (18). Figure 3 shows the equivalent kinematic viscosity obtained for each scheme. Note that performing this regression on a more global variable such as the kinetic energy would have also been possible. After verifications both methods revealed to provide sensibly the same results. The equivalent viscosity is nearly two orders of magnitude smaller for the WENO 5 scheme than for the MUSCL scheme. For a 100x100 grid resolution, the extrapolated equivalent kinematic viscosity is 10−5 m2 /s for the WENO 5 scheme compared to 2 × 10−4 m2 /s and 2 × 10−3 m2 /s for the MUSCL scheme with Van Leer and Minmod limiters respectively. On this test case, the WENO 5 scheme proved to be quantitatively less dissipative, thus meeting our needs for a low-diffusive solver. 2.4. Single vortex convection test case The ability to properly convect vortices is critical for modeling turbulent flows and wakes. The second diffusion test case presented in this paper refers to a benchmark proposed by the CERFACS institute [1]. L It consists in a single vortex of radius Rc = 20 convected with a velocity U0 along ~x within a periodic 6

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

7

Figure 3: 2D Taylor-Green vortex test case. Equivalent kinematic viscosity computed with MUSCL (Minmod and Van Leer limiters) and WENO 5 schemes. For information purpose, an equivalent order of numerical diffusion rate is plotted. square domain of size L. In the particular case of an inviscid fluid, the expected solution is the initial vortex convected without deformation. The equations relative to this test case are the following:   u = U0 + ∂Φ  ∂y  v = − ∂Φ ∂x (19) (x−xc )2 +(y−yc )2   2 Rc  P = − ρ0 Γ22 e− (

2Rc

p if (x − xc )2 + (y − yc )2 < 4Rc Γe 0 otherwise The parameters are defined as: U0 = 35m/s, L = 0.3112m, ρ0 = 1.7170kg/m3 , Γ = 0.0359157m2 .s−1 . The simulation is performed under the weakly-compressible assumption with a nominal sound speed c0 = 350m/s and γ = 1.4. Figure 4 presents a comparison for a 128x128 grid resolution between the expected pressure field solution and the results obtained with MUSCL (with Van Leer limiter) and WENO 5 schemes at t = UL0 (vortex displacement of one domain length). These snapshots clearly emphasize the limitation of the MUSCL scheme in vortex convection situations, as a strong vortex diffusion is observed together with the presence of dispersive effects. Figure 5 shows the evolution of the velocity and pressure profiles as a function of the grid resolution at t = 40 UL0 (vortex displacement of 40 domain lengths). For a 128x128 grid resolution, errors in the velocity amplitudes are lower than 10% in the WENO 5 case instead of 66% for the MUSCL scheme. WENO 5 scheme results with 32x32 are almost equivalent to those obtained with the MUSCL scheme with a grid 128x128. This interesting test case also enables to compute the heuristic order of convergence of the solver. This convergence study has been performed on the velocity profiles U − U0 on the x = 0 line. The reference solution used for computing the order is chosen as the initial velocity profile. Results are computed for different error norms and are presented in Table 1. For MUSCL, convergence orders of about 1.6 are obtained, slightly lower than the theoretical order of 2. Same results were obtained when computing the convergence order on the pressure profiles. It seems 7 with Φ(x, y) =

(x−xc )2 +(y−yc )2 − 2 Rc

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

8

Figure 4: Single vortex convection test case. Pressure field at t = UL0 for a 128x128 grid resolution: expected pressure field (left), MUSCL scheme solution with Van Leer limiter (middle) and WENO 5 solution (right).

Figure 5: Single vortex convection test case. Grid convergence on pressure (WENO 5 only) and U − U0 (WENO 5, and MUSCL for 128x128) fields along x = 0 line at t = 40 UL0 . Norm L1 L2 Inf

MUSCL spatial convergence 1.66 1.57 1.30

WENO 5 spatial convergence 3.36 3.24 3.14

Table 1: Single vortex convection test case. Spatial convergence orders obtained with MUSCL and WENO 5 schemes on velocity profiles U (y) along x = 0 line. therefore that performing turbulence simulations with reasonable computational costs cannot be achieved with this scheme because of its strongly diffusive behavior, namely for wake turbulence studies. Conversely, the WENO 5 scheme exhibits an overall scheme convergence order higher than 3 for any error norm. Finally, Figure 6 displays a cost vs accuracy study underlining the benefits of using a higher order scheme compared to the classical MUSCL scheme on this convection case. This cost/benefit study shows that for a given computational cost, the WENO 5 scheme is much more accurate than the MUSCL scheme. As a consequence, the WENO 5 scheme has to be recommended for 8

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

10

100

100

Relative error (%) on maximum velocity value

32x32

Equivalent sequential CPU time (s) 1000

10000

9

100000

128x128 128x128

64x64

64x64

256x256

10

512x512 128x128

1

256x256

0,1

MUSCL Minmod

MUSCL VanLeer

WENO5

Figure 6: Single vortex convection test case. Precision as a function of the computational cost for different solvers and grid resolutions. simulations in which the numerical diffusion displays significant consequences on the flow prediction. 3. Treatment of viscous terms 3.1. Viscous flux solver The whole Navier-Stokes equations are addressed here by adding the viscous terms to the Euler equations considered previously. By using the weakly-compressible approach (compressibility effects are neglected), the expression of the viscous terms is simplified by setting div ~u = 0. In an integral form, and by using an analogy to Equation (14), the Navier-Stokes equations can be rewritten as: ZZZ ZZ ZZ d W dV + Ψc (W ). ~ndS = Ψv (W ). ~ndS, (20) dt V S S with Ψv (W ) the viscous flux tensor defined as: 



  0 0 0     2µu,x µ(v,x + u,y ) µ(w,x + u,z )   Ψv (W ) =   2µv,y µ(w,y + v,z )   µ(u,y + v,x )   2µw,z  µ(u,z + w,x ) µ(v,z + w,y )  | {z }| {z }| {z } Fv (W )

Gv (W )

Hv (W )

with Fv (W ), Gv (W ) et Hv (W ) the viscous fluxes in directions x,y and z. 9

(21)

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

10

The above velocity gradients are discretised with the scheme proposed by Berger et al. [4]. This discretisations is briefly presented here, by showing the example of the computation of the derivatives of u around a face of a cell of normal ~x in two dimensions. Two kinds of derivatives can be distinguished for which the computational method differs: the normal derivative to the face u,x and the derivatives in transverse directions (here u,y ).

Figure 7: Schematic of the computation of u velocity derivatives on a face of normal ~x in two dimensions. (

u

−u

i,j (ui+ 12 ,j ),x = i+1,j ∆x u −u +ui,j+1 −ui,j−1 ) (ui+ 12 ,j ),y = 12 ((ui,j ),y + (ui+1,j ),y ) = 12 ( i+1,j+1 i+1,j−1 2∆y

(22)

This solver shows the advantages of directly computing the first derivatives. These viscous terms are then treated similarly to the hyperbolic fluxes in Equation (15). Finally, as for the hyperbolic part the explicit nature of the temporal scheme implies the respect of the following viscous CFL condition: ∆tvisc = CF Lvisc

ρ(mini (∆xi ))2 , µ

(23)

in which the diffusion number CF Lvisc is set to 0.125, as usually adopted in weakly-compressible approaches (see for instance [26, 2]). Finally, the time step considered for time integration is chosen as: ∆t = min(∆tacous , ∆tvisc ).

(24)

This additional stability condition can become preponderant for very small cell sizes and/or low Reynolds numbers. Nevertheless, in practice all test cases performed up to now were constrained by the acoustic CFL condition only, due to the low viscosity of the fluids involved. 3.2. Wall boundary conditions A ghost cell method is used to impose wall boundary conditions onto the domain limits. The ghost cells variables (velocity and pressure) are deduced from the near boundary fluid cells as:   Pg = Pf ~vg = 2~vb − ~vf no-slip condition (25)  ~vg = ~vf − 2~n.(~vf − ~vb )~n free-slip condition

where subscripts g, f and b refer respectively to ghost cell, fluid cell and boundary. Note that this method exclusively addresses plane boundaries in the present study (it differs from the treatment of boundary conditions on complex geometries which is described further in this paper).

10

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

11

3.3. Validation on the lid-driven cavity at Re = 1000 The lid-driven cavity at Re = 1000 [12] is chosen as a validation test case for the viscous solver. A square shaped fluid cavity of length L = 1m is driven by the upper plane moving at a constant imposed velocity Ue = 1m.s−1 . No-slip conditions are applied on each boundary, and the fluid features are the following: ρ0 = 1000kg.m−3 and µ = 1kg/m.s. The sound speed used here is c0 = 10m.s−1 (i.e. 10 times the maximal expected flow velocity). The fluid is initially at rest, and 40 seconds are necessary to reach the steady solution. Figure 8 shows a snapshot of the flow obtained with the WENO 5 scheme for a 128x128 grid resolution.

Figure 8: Lid-driven cavity test case at Re = 1000. Streamlines and velocity magnitude fields obtained with the present WENO 5 scheme for a 128x128 grid resolution. For this case the presence of primary, secondary and tertiary vortices is expected, for which the positions are studied. Figure 9 presents the results obtained with MUSCL (Van Leer limiter) and WENO 5 schemes, compared to a highly accurate reference solution from Botella et al. [7] resulting from a Chebyshev collocation method with a polynomial of degree N = 160. In this figure, the relative errors Err(%) on the vortex center locations are presented, with |X − Xref | , (26) L where X and Xref represent respectively the computed and the reference vortex center locations. As expected, some differences appear between MUSCL and WENO 5 scheme results, the latter presenting lower errors due to its lower numerical diffusion. In addition to the above results, a convergence study performed on the velocity profile V (x) along the y = 0.5m axis is presented in Figure 10. MUSCL and WENO 5 scheme results are compared to both Botella et al. and Ghia et al. [12] solutions. The solution provided by the WENO 5 scheme with a 128x128 resolution is by far closer to the reference solutions than the MUSCL scheme for this same grid as well as for a 512x512 grid, emphasizing the interest of a higher order spatial scheme also for this case. Results do converge towards the reference solution but the use of the WENO 5 scheme together with a very fine spatial resolution (512x512) is necessary to approach the full convergence. It is further observed that the solution from Ghia et al. is not fully converged, confirming the need for the highly accurate reference solution from Botella et al. [7]. Note that the WENO 5 solution on a 512x512 grid is very close to this reference profile. Finally, convergence orders of the two respective spatial schemes tested are presented in Table 2. Err(%) = 100

11

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

Vortices Primary

Re=1000

Botella et al.

MUSCL

Coordinates

N=160

128x128

x (m)

0,5308

0,5358

0,4987

0,5312

0,0400

0,5309

0,0100

0,5657 0,0829 0,0772 0,8640 0,1120 0,9914 0,0082

0,0500

0,5653

0,0100

0,0429

0,0830

0,0300

0,0870

0,0778

0,0300

0,0000

0,8641

0,0100

0,0200

0,1119

0,0100

0,0920

0,9923

0,0020

0,0550

0,0074

0,0250

Err (%)

WENO5

12

Err (%)

128x128

y (m)

0,5652

0,5632

0,1953

Secondary left

x (m)

0,0833

0,0809

0,2400

y (m)

0,0781

0,0721

0,5962

Secondary right

x (m)

0,8640

0,8747

1,0660

y (m)

0,1118

0,1121

0,0300

Tertiary (right)

x (m)

0,99232

0,9870

0,5320

y (m)

0,00765

0,0130

0,5300

WENO5

Err (%)

512x512

Figure 9: Lid-driven cavity test case at Re = 1000. Comparison of the vortex position obtained by Botella et al. (Chebyshev collocation method with a polynomial of degree N = 160) taken as reference and with present MUSCL and WENO5 schemes. Spatial scheme MUSCL WENO5

Spatial convergence order 1.59 2.44

Table 2: Lid-driven cavity test case at Re = 1000. Spatial convergence order for MUSCL and WENO 5 schemes. As expected, a convergence order lower than 2 is obtained for the MUSCL solver. Concerning the WENO 5 scheme, the convergence order is lower than the expected theoretical value. This could be explained by the order of the viscous solver used here (2nd order), but this assumption should be further confirmed by performing a convergence study with a higher order viscous solver. As expected, the choice of the spatial scheme adopted within the hyperbolic part has a leading role for the simulation of viscous flows. The differences between MUSCL and WENO solvers are here explained µ by the previously performed study by looking at the very low physical ratio. The WENO 5 solver should µnum therefore be preferred for cases involving higher Reynolds numbers. 3.4. Parallelization and Mesh Refinement strategies To reduce the computational times while preserving a sufficient accuracy, an octree-based hierarchical structure has been developed in order to refine the spatial resolution in desired areas of interest while derefining everywhere else. Note that static mesh refinements are used exclusively in the present work, as the implementation of efficient AMR refinement criteria will be the topic of future investigations. The solver has also been integrated in a parallel environment by the mean of MPI libraries using a domain decomposition strategy. Ordering, refinement level management and load-balancing of octree elements throughout the process network are performed through the use of a Z-order space filling curve, similarly to the method adopted in [18]. The explicit feature of WCCH is leveraged to overlap the MPI communication latencies through a “Local-Inner-Outer” framework, similarly to the strategy proposed for the SPH method in [28]. A scalability study of WCCH has been performed on the CRIHAN computational cluster, equipped with Westmere EP cores together with Infiniband switch interconnects. Figure 11 shows the speedup obtained in the 3D case on up to 1024 cores, with a static mesh refinement. Both strong and weak scalability speedups are plotted in this figure, which illustrates the high parallel performances obtained in terms of MPI communication latency overlapping. In order to complete this study, a performance index has been measured in the above conditions, leading to an average CPU time TCP U = 6.8 × 10−5 s.core/cell/∆t. A super-linear speedup is obtained in all cases tested in strong scalability analysis. This can be explained by beneficial cache effects occurring as the buffer size decreases in proportion to the number of cores involved. Good weak-scalability performances are also observed, with a parallel 12

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

13

Figure 10: Lid-driven cavity test case at Re = 1000. Grid convergence on the V(x) velocity profile at y = 0.5m for the MUSCL and WENO solvers. Top left: global view. Top right: zoom around x ∈ [0.05; 0.25]. Bottom: zoom around x ∈ [0.75; 1]. efficiency greater than 95% in any cases involving from 32, 000 cells per core to higher process loads. Note that the parallel strategy adopted in WCCH is not affected by mesh refinements, at least in a static case. Further scalability studies should be performed to analyze the parallel performances obtained on mesh refinements adapted dynamically. 4. Boundary conditions on complex geometries This part is dedicated to complex geometry handling with a Cartesian grid, which is one of the fundamentals of the solver presented in this paper. The desired properties for our solver are the following: • robustness of the method • algorithmic simplicity of the method, compatible with the use of a parallel mesh refinement framework including the 5th order WENO scheme • flow simulation results independent from the position of the geometries relatively to the Cartesian grid 13

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

14

Figure 11: Strong (left) and weak (right) scalability studies performed on a 3D case (WENO 5 scheme) on up to 1024 cores and involving a static mesh refinement. • easy extension to 3-D computations • low overall computational costs • possibility to perform viscous and turbulent simulations 4.1. Quick review of the related literature Various approaches are possible for handling boundary conditions on complex geometries. The most common method is the use of body-fitted grid, offering the advantages of a good compromise between accuracy and number of cells involved. While it is commonly used in free-surface hydrodynamics, this method also presents some drawbacks such as mesh creation complexities, impossibility to treat relative motions of mobile geometries, difficulties to manage adaptive mesh refinement and/or high order spatial schemes. An alternative to body-fitted methods is the cut-cell method, which consists in explicitly inserting the geometry in the Cartesian grid, leading to cells cut into polygonal shapes. The main advantage of this method is the accurate preservation of the boundary position and its conservativeness. In return, some difficult issues appear due to its algorithmic complexity (namely in the 3-D case) and the need for cell merging procedures to avoid CFL restriction on small cells [15, 16, 21, 39]. Some other methods have been designed to avoid any grid conformation, such as the ghost cell method. This method imposes the boundary conditions through a combination of near-boundary fluid cell variables and boundary kinematic properties, which is applied in the cells located inside the geometry (ghost cells) [6, 20, 22]. Another possible approach is the Immersed Boundary Method (IBM), for which the boundary is not explicitly embedded in the grid but is represented by adding a fictitious force term modifying the momentum equations in the near-boundary region to enforce the desired boundary condition [29, 23]. The main advantages of this method reside in its relative simplicity together with the absence of limitation on the time-steps. Force terms can be added in a continuous way as in its original version from Peskin [29], with possible improvements as in the “feedback forcing” method proposed by Goldstein et al. [13], or alternatively in a discrete way as in the so-called “direct forcing” method [10, 24, 34, 35, 38]. Note that while an assimilation is often made by authors between these two methods, we make here a distinction between the ghost cell method and the IBM. Both methods use flux completion thanks to the cells located inside the geometry, but the main difference resides on the imposition of variables of ghost 14

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

15

cells, which strongly differs from the IBM in which a source term is added in the momentum equation. This difference is particularly important for (weakly-)compressible schemes (as in the present work) where mass fluxes act in the mass conservation equation. 4.2. Direct forcing IBM proposed The method proposed in this paper is an Immersed Boundary Method with explicit direct forcing that is close to the methods proposed in [10] and [34]. The boundary condition is enforced by adding a forcing source term within the momentum equation in all cells located inside the geometry and in its near outside vicinity. An originality of the method proposed is the addition of a regularization function to balance smoothly the momentum fluxes and the forcing source term in the near-boundary region. Another peculiarity of the method proposed resides on the application of the IBM method onto a (weakly-)compressible flow, contrasting with the IBM methods usually encountered in the literature. To this end, no special treatment is applied to the mass conservation equation, leaving the density (i.e. the pressure) unconstrained throughout the fluid domain, including inside the geometry. The pressure inside the geometry therefore adapts and tends to naturally impose the desired boundary condition. This method is therefore both simple and computationally inexpensive. We first introduce a Heaviside regularization function H(δ), with δ the non-dimensional normal distance to the boundary, normalized with respect to the characteristic size ∆x of the grid in the near boundary vicinity. H(δ) tends towards 0 far from the body (δ > 0) and tends towards 1 inside the geometry (δ < 0). This Heaviside function is smoothed using a hyperbolic tangent function as follows: 1 (1 + tanh (−Kδ)) , (27) 2 with K a dimensionless coefficient allowing to adjust the regularization width ζ, as presented in Figure 12. H(δ) =

Figure 12: Regularization function used for the direct forcing method. At instant tn , the acceleration of a given cell (ijk) located anywhere in the fluid domain (including inside or onto a solid geometry) can be decomposed as n  n  n  d ~vijk d ~vijk d ~vijk = (1 − H(δijk )) + H(δijk ) , (28) dt dt dt f luid body  n d~ vijk by weighting the fluid and solid contributions on this cell. The fluid contribution is obtained dt f luid

through

15

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22



d ~vijk dt

n

=

f luid

1 ρnijk



d (ρ~v )ijk dt

n

n − ~vijk



d ρijk dt

n 

−−−→ = RHS ijk ,

16

(29)

−−−→ wherethe right n hand side term RHS ijk corresponds to the hyperbolic and elliptic contributions deduced d Wijk from detailed in sections 2 and 3. dt n  d~ vijk corresponds to the forcing source term that imposes the body accelThe solid contribution dt body

eration and which is approximated as



d ~vijk dt

n

body

≈

n+1 n ~vijk − ~vijk

∆t

,

(30)

with ∆t = t n+1 − tn the current time step (i.e. actually the sub-timestep of each Runge-Kutta stage). We assume n+1 ~vijk ≈ v~b n+1 ,

(31)

with v~b n+1 the body velocity at t n+1 , so that a cell (ijk) fully embedded within the geometry sees 

d ~vijk dt

n

body

≈

n v~b n+1 − ~vijk ~ijk . =S ∆t

(32)

The direct forcing immersed boundary method is therefore applied here at time tn by weighting a forcing −−→ ~ijk and the right hand side terms − source term S RHS ijk within the computational grid through: n  −−−→ d ~vijk ~ijk . = (1 − H(δijk ))RHS ijk + H(δijk )S (33) dt ~ijk is Note that Eq.(33) is finally integrated in time to update the velocity. The forcing source term S therefore applied in a continuous and progressive manner through the geometry boundary. This boundary condition is a no-slip boundary condition both from the hyperbolic and viscous part points of view. The bigger ζ, the larger the transition zone between the fluid and the geometry. The influence of the boundary is thus spread on several cells. As a result, the boundary can be considered as distributed in the Cartesian ~ijk over all cells represents the total force exerted by the body grid. Note that the sum of the terms H(δijk )S on the fluid. Note finally that the above formulation is first order accurate only, since the approximation (31) does not take into account the relative position between the solid boundary and the cell centroids. The use of the Heaviside regularization function helps in maintaining an admissible order and in preventing the scheme from spurious pressure oscillations, but our future investigations will aim at increasing the order of the method by getting rid off this limitation, as proposed for instance in [3]. Meanwhile, we present the results obtained with this formulation on the validation test cases proposed hereafter. 4.3. Validation on the flow around a circular cylinder at Re = 40 4.3.1. Validation of the present results with reference values This direct forcing Immersed Boundary Method is applied here to the case of a flow around a fixed cylinder at Re = 40. The sound speed used here is c0 = 20U0 , where U0 denotes the incident x-velocity. The flow is initially set as uniform with U0 as x-velocity and a null pressure field. The steady solution is fully reached at tU D = 50, D denoting the diameter of the cylinder. The results obtained are compared with the literature in terms of forces and geometric characteristics of the flow as shown on Figure 13. The left plot of Figure 14 displays the convergence of the drag coefficient Cd for three grid resolutions. The drag coefficient converges towards the reference value Cd = 1.48. The right plot of Figure 14 compares the pressure profile obtained on the finest grid with reference pressure profiles extracted from the literature ([9], [34], [14] and [6]), showing a very good agreement. Note that some oscillations can be observed 16

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

17

b/2D

θ a/D

L/D

Figure 13: Characteristics of the flow past a cylinder at Re = 40.

Figure 14: Flow past a fixed cylinder at Re = 40. Left: convergence on the drag coefficient Cd. Right: comparison of the pressure profiles around the cylinder with regard to reference results. between the angles 60◦ and 120◦ . These effects are due to the relatively small regularization width used here (ζ=2), as complementary numerical tests showed that increasing ζ remarkably improves the pressure profile around the cylinder by removing such oscillations. Nevertheless, this also leads to significantly smear the fluid-solid interface and tends to decrease the convergence rate of the scheme. Table from Figure 15 shows the different quantitative results obtained (Cd, geometrical features of the flow) for the different numerical configurations tested, and shows the influence of the different parameters used for this method. The values obtained from the simulations are compared to experimental values ([33] and [8]), and the range of results usually obtained in the literature ([11], [9], [34], [14] and [6]) are also presented. For the cases used for convergence (first three cases) all quantities converge towards the experimental values, in terms of forces, vortex positions and separation angles. Only the vortex length is different from the experimental values but is close to the values from numerical computations found in the literature. A case with a resolution ∆Dx = 30 was performed with the MUSCL scheme to evaluate the effects of the spatial order on the solution, leading to a 1.55% error which is not very different from the 0.45% error with WENO 5, as expected for these low Reynolds number. An other computation was performed without using the regularization function H(δ), leading to a 6% error on the drag forces with respect to the experimental values and to major differences in terms of flow features, emphasizing the interest of a regularization function. Figure 16 compares the pressure fields obtained for ∆Dx = 30 and ∆Dx = 60 in the cylinder region. 17

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

18

Tritton (Exp.)

1.48

-

-

-

-

Error on (%) -

Bouard & Coutanceau (Exp.)

-

53.8

2.13

0.76

0.59

-

Numerical Experiments (Range) WENO D/Δx= 15 WENO D/Δx = 30 WENO D/Δx = 60 MUSCL D/Δx = 30 WENO D/Δx = 30

1,491,62

52,755,6

2,142,35

0,710,76

0,590,60

-

1.522

51.5

2.54

0.84

0.63

2.82

1.487

52

2.33

0.75

0.61

0.45

1.486

53.7

2.27

0.72

0.6

0.39

1.503

51.85

2.34

0.77

0.6

1.55

1.567

50.3

2.22

0.7

0.59

5.90

Case

References

Mesh Convergence WENO

MUSCL No Regularization

Figure 15: Flow past a fixed cylinder at Re = 40. Comparison of forces and geometrical flow features with respect to the reference solutions.

Figure 16: Flow past a fixed cylinder at Re = 40. Effect of the spatial resolution in the cylinder region. Pressure fields are obtained with ∆Dx = 30 (left) and ∆Dx = 60 (right). The pressure field inside the cylinder (forced cells) is displayed. The effects of convergence are visible, namely by improving the pressure field regularity along the cylinder boundary. This figure also displays the pressure field inside the cylinder, illustrating the natural pressure adaptation that imposes the desired boundary condition, since no forcing term is applied to the mass conservation equation for inner cells. 4.3.2. Comparison between fixed and moving case solutions The comparison of the cases of the flow past a fixed cylinder with the expected equivalent case of a moving cylinder in the fluid at rest is now performed . Inlet/outlet boundary conditions are applied respectively on 18

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

19

the left and right domain boundaries for the fixed cylinder case, while the moving cylinder should perform a long displacement through the grid before reaching the Re = 40 steady state. A very long domain is thus required for this latter case, implying the need for the refinement features of WCCH, as shown in Figure 17.

Figure 17: Grids used for the comparison between fixed (left grid) and moving (right grid) cylinder. Comparisons are performed while the moving cylinder has reached the origin coordinates. For this comparison the grid resolution in the cylinder area is only ∆Dx = 30. Indeed, the moving cylinder case requires a long domain that would result in large computational costs for a higher spatial resolution (not especially needed for this fixed/moving geometry comparison purpose). Figure 18 shows a comparison at time tU D = 50 between the pressure fields of the mobile and fixed cases, with a regularization width ζ = 2.

Figure 18: Comparison of the pressure fields of the fixed (left) and mobile (right) cylinder cases with a direct forcing method with ζ = 2. Both cases use a spatial resolution ∆Dx = 30 in the cylinder region. Very small differences are observed between both pressure fields, showing a very good agreement between fixed and mobile cases. This very good agreement is confirmed by the left part of Figure 19, comparing the pressure profiles for both simulations. 19

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

1

20

Fixed cylinder Moving cylinder

Cp

0.5

0

-0.5

0

50

100

150

Angle (°)

Figure 19: Comparison of the pressure profiles between the fixed and mobile cylinder simulations (left). Comparison of drag forces for various regularization widths ζ (right). A spatial resolution ∆Dx = 30 is used in the cylinder region. The use of a regularization function H(δ) allows to smooth the oscillations in the pressure fields and therefore on the resulting forces. The right plot of Figure 19 shows a zoom of the drag force time history obtained for various values of regularization width ζ. For a regularization width large enough (ζ = 2), the oscillations around the mean drag force value are strongly damped. For smaller values of ζ, oscillations arising from the cells transitioning from solid to fluid cells during the cylinder displacement are responsible for force variations from 7% of the mean temporal value up to 20% according to the value chosen for ζ. 5. Conclusion The WCCH method dedicated to solving complex hydrodynamic flows has been introduced in the present paper. The long-term targeted applications are viscous and turbulent flows in presence of complex moving geometries while releasing the meshing constraints. This original strategy is based on a fully-explicit scheme through a weakly-compressible approach, which strongly differs from the implicit incompressible formulations usually adopted in most hydrodynamic CFD solvers. The current implementation relies on a purely Cartesian grid facilitating the implementation of a massively parallel framework together with local mesh refinement capabilities. The discretization of the hyperbolic part of the Navier-Stokes equations has been first studied, with a particular attention paid on the numerical diffusion of the spatial scheme. The use of a low-order upwind method induces a large numerical diffusion that might become much larger than the physical flow viscosity. A 2nd -order (MUSCL) scheme has been first studied on steady and convection vortex test cases, for which MUSCL proved to be too much diffusive. A 5th -order (WENO) spatial scheme has been then used, showing a strong minimization of the numerical diffusion. The study then focused on viscosity modeling and its validation on a Re = 1000 lid-driven cavity. Results obtained were very close to the reference solutions and showed that the final accuracy is mainly correlated to the ratio between physical and numerical viscosity. The paper finally addressed the imposition of boundary conditions of complex solid geometries immersed within the Cartesian grid. The method proposed consists in a direct forcing immersed boundary technique allowing for a simple algorithmic that is fully compatible with a high-order spatial scheme (WENO5) and with the required mesh refinement and parallel features. Validations based on flows past a fixed and a mobile cylinder at low Reynolds number Re = 40 were performed and showed some good convergence properties, together with a large equivalence between mobile and fixed case, as expected. This method will need to be further studied, namely in cases involving unsteady mobile geometries. The future researches will be 20

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

21

mainly oriented towards the design of an efficient AMR refinement criteria in order to fully exploit the potential of this method, and towards modeling turbulent flows for high Reynolds numbers around complex mobile geometries. Moreover, a formulation based on the Volume Of Fluid (VOF) approach is also under development for the simulations of complex two-phase flows involving high density ratios. References [1] http://elearning.cerfacs.fr/numerical/benchmarks/vortex2d/. [2] M. Antuono, A. Colagrossi, and S. Marrone. Numerical diffusive terms in weakly-compressible SPH schemes. Computer Physics Communications, 183:2570–2580, 2012. [3] E. Balaras. Modeling complex boundaries using an external force field on fixed Cartesian grids in large-eddy simulations. Computers and Fluids, 33:375–404, 2004. [4] M. Berger, M.J. Aftosmis, and S. Allmaras. Progress towards a cartesian cut-cell method for viscous compressible flow. In Proceedings of the AIAA Conference Paper2012-1301, 2012. [5] M. J. Berger and P. Colella. Local Adaptive Mesh Refinement for shock hydrodynamics. Journal of Computational Physics, 82:64–84, 1989. [6] P.A. Berthelsen and O.M. Faltinsen. A local directional ghost-cell approach for incompressible viscous flow problems with irregular boundaries. Journal of Computational Physics, 227:4354–4397, 2008. [7] O. Botella and R. Peyret. Benchmark spectral results on the lid-driven cavity flow. Computers and Fluids, 27:421–433, 1998. [8] R. Bouard and M. Coutanceau. The early stage of development of the wake behind an impulsively started cylinder for 40 < Re < 104 . J. Fluid Mech., 101:583–607, 1980. [9] S. Dennis and G. Chang. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J. Fluid Mech., 42:471–489, 1970. [10] E.A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined immersed-boundary finite-difference methods for three dimensional complex flow simulations. Journal of Computational Physics, 161:35–60, 2000. [11] R. Gautier, D. Biau, and E. Lamballais. A reference solution for the flow over a circular cylinder at Re=40. Computers and Fluids, 75:103–111, 2013. [12] U. Ghia, K.N. Ghia, and C.T. Shin. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 48:387–411, 1982. [13] D. Goldstein, R. Handler, and L. Sirovich. Modeling a no-slip flow boundary with an external force field. Journal of Computational Physics, 105:354–366, 1993. [14] A.S Grove, F.H. Shair, E.E. Petersen, and A. Acrivos. An experimental investigation of the steady separated flow past a circular cylinder. J. Fluid Mech., 19:60–80, 1964. [15] X. Hu, B. Khoo, N. Adams, and F. Huang. A conservative interface method for compressible flows. Journal of Computational Physics, 219:553–578, 2006. [16] M. Kirkpatrick and S. Armfield. A representation of curved boundaries for the solution of the Navier-Stokes equations on a staggered three-dimensional cartesian grid. Journal of Computational Physics, 184:1–36, 2003. [17] X.-D. Liu, S. Osher, and T. Chan. Weighted essentially non-oscillatory schemes. Journal of Computational Physics, 115:200–212, 1994. [18] P. MacNeice, K. M. Olson, C. Mobarry, R. de Fainchtein, and C. Packer. PARAMESH: A parallel adaptive mesh refinement community toolkit. Computer Physics Communications, 126:330–354, 2000. [19] S. Marrone, A. Colagrossi, A. Di Mascio, and D. Le Touz´ e. Prediction of energy losses in water impacts using incompressible and weakly compressible models. Journal of Fluids and Structures, 54:802–822, 2015. [20] D.D. Marshall and S.M. Ruffin. An embedded boundary Cartesian grid scheme for viscous flows using a new viscous wall boundary condition treatment. In AIAA-2004-0581, 2004. [21] H. Udaykumar R. Mittal, P. Rampunggoon, and A. Khanna. A sharp interface cartesian grid method for simulating flows with complex moving boundaries. Journal of Computational Physics, 174:345–380, 2001. [22] R. Mittal, H. Dong, M. Bozkurttas, F.M. Najjar, A. Vargas, and A. von Loebbecke. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. Journal of Computational Physics, 227:4825–4852, 2008. [23] R. Mittal and G. Iaccarino. Immersed boundary methods. Annual Review of Fluid Mechanics, 37:239–261, 2005. [24] J. Mohd-Yosuf. Combined immersed boundary/B-spline methods for simulation of flow in complex geometries. Annu. Res Briefs, Cent. Turbul. Res. pp 317-28, 1997. [25] J. J. Monaghan. Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys., 30:543–574, 1992. [26] J.P. Morris, P.J. Fox, and Y. Zhu. Modeling low Reynolds number incompressible flows using SPH. Journal of Computational Physics, 136:214–226, 1997. [27] G. Oger, S. Marrone, D. Le Touz´ e, and M. de Leffe. SPH accuracy improvement through the combination of a quasiLagrangian shifting transport velocity and consistent ALE formalisms. Journal of Computational Physics, 313:76–98, 2016. [28] G. Oger, D. Le Touz´ e, D. Guibert, M. de Leffe, J. Biddiscombe, J. Soumagne, and J.-G. Piccinali. On distributed memory MPI-based parallelization of SPH codes in massive HPC context. Computer Physics Communications, 200:1–14, 2016. [29] C.S. Peskin. Flow patterns around heart valves: a numerical method. Journal of Computational Physics, 10:252–271, 1972.

21

P. Bigay et al. / Computer Physics Communications 00 (2017) 1–22

22

[30] P.L. Roe. Characteristic-based schemes for the Euler equations. Annual Review of Fluid Mechanics, 18:337–365, 1986. [31] V.A Titarev and E.F Toro. Finite-Volume WENO schemes for three dimensional conservation laws. Journal of Computational Physics, 201:238–260, 2004. [32] E.F Toro. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media, 1999. [33] D.J. Tritton. Experiments on the flow past a circular cylinder at low Reynolds numbers. J. Fluid Mech., 6:547–567, 1959. [34] Y.H. Tseng and J.H. Ferziger. A ghost-cell immersed boundary method for flow in complex geometry. Journal of Computational Physics, 192:593–623, 2003. [35] M. Uhlmann. An immersed boundary method with direct forcing for the simulation of particulate flows. Journal of Computational Physics, 209:448–476, 2005. [36] B. van Leer. Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme. Journal of Computational Physics, 14:361–370, 1974. [37] B. van Leer. Towards the ultimate conservative difference scheme V. A second order sequel to Godunov’s method. Journal of Computational Physics, 32:101–136, 1979. [38] R. Verzicco, J. Mohd-Yusof, P. Orlandi, and D. Haworth. Large Eddy Simulation in complex geometric configurations using boundary body forces. AIAA journal, 38:427–433, 2000. [39] T. Ye, H.Udaykumar, and W. Shyy. An acurate cartesian grid method for viscous incompressible flows with complex immersed boundaries. Journal of Computational Physics, 156:209–240, 2.

22