An immersed boundary method in OpenFOAM : Verification and validation

Accepted Manuscript

An Immersed Boundary Method in OpenFOAM : verification and validation E. Constant, J. Favier, M. Meldi, P. Meliga, E. Serre PII: DOI: Reference:

S0045-7930(17)30275-X 10.1016/j.compfluid.2017.08.001 CAF 3556

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

8 March 2017 10 June 2017 4 August 2017

Please cite this article as: E. Constant, J. Favier, M. Meldi, P. Meliga, E. Serre, An Immersed Boundary Method in OpenFOAM : verification and validation, Computers and Fluids (2017), doi: 10.1016/j.compfluid.2017.08.001

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.

ACCEPTED MANUSCRIPT

Highlights • A modified PISO algorithm integrating an efficient Immersed Boundary method is proposed, using an improved direct forcing approach.

CR IP T

• A rigorous characterization of the Immersed Boundary method is performed using an original verification technique. • The numerical errors at various stages of the algorithm are precisely estimated.

AC

CE

PT

ED

M

AN US

• A thorough validation of the solver on relevant literature test-cases is provided

1

ACCEPTED MANUSCRIPT

An Immersed Boundary Method in OpenFOAM : verification and validation E. Constanta , J. Faviera,∗, M. Meldib , P. Meligaa , E. Serrea a Aix-Marseille Univ., CNRS, Centrale Marseille, M2P2 Marseille, France Institut P 0 , CNRS – Universit´e de Poitiers – ENSMA, 11 Bd Marie et Pierre Curie, 86962 Futuroscope Chasseneuil Cedex, France

CR IP T

b

Abstract

AC

CE

PT

ED

M

AN US

The present work proposes a modified Pressure-Implicit Split-Operator (PISO) solver integrating the recent Immersed Boundary Method (IBM) proposed by Pinelli et al. [1] in order to perform reliable simulations of incompressible flows around bluff bodies using the open source toolbox OpenFOAM version 2.2 (ESI-OpenCFD [2]). The (IBM) allows for a precise representation of fixed and moving solid obstacles embedded in the physical domain, using uniform or stretched Cartesian meshes. Owing to this feature, the maximum level of accuracy and scalability of the numerical solvers can be systematically achieved. An iterative scheme based on sub-iterations between (IBM) and pressure correction has been implemented in the native (PISO) solver of OpenFOAM. This allows one to use fast optimized Poisson solvers while satisfying simultaneously the divergence-free flow state and the no-slip condition at the body surface. To compute the divergence of the momentum equation (in the PISO loop) and the interpolation of the fluxes, we propose an hybrid calculation with an analytical resolution (using the kernel function equation) of the quantities involving the force term (singular quantities). A careful and original verification study has been carried out which allows to estimate three different errors related to the discretization and to the (IBM). Various 2D and 3D well-documented test cases of academic flows around fixed or moving solid bodies (cylinder and sphere) have been simulated and carefully validated against existing data from the literature in a large range of Reynolds numbers, Re = 30 − 3900 and in the frame of DNS and DDES OpenFOAM native models. Keywords: Immersed Boundary Method (IBM), OpenFOAM, Bluff body, ∗

Corresponding author Email address: [email protected] (J. Favier)

Preprint submitted to Computers & Fluids

August 4, 2017

ACCEPTED MANUSCRIPT

Incompressible flows 1. Introduction

AC

CE

PT

ED

M

AN US

CR IP T

Flows around bluff body are of practical interest for many engineering applications. Their simulations require to develop and implement accurate and efficient numerical methods in versatile and powerful numerical tools that can be further used to study real industrial problems. The immersed boundary methods have shown in the past few years to be efficient to encompass the presence of fixed and moving solid obstacles in a computational mesh without conforming to their boundaries. A wide spectrum of methods included in this family have proven efficient to simulate complex and moving geometries, such as Lagrangian multipliers [3], level-set methods [4], fictitious domain approaches and surface [5] and volume penalization approaches [6, 7]. A common feature of all (IBM) techniques is that the Navier-Stokes equations are discretized over a simple structured Cartesian grid, which significantly improves the computational efficiency and the stability. The complex geometry is then immersed into a larger computational domain, and the boundary conditions are represented by the addition of an ad-hoc body force in the momentum equations. Such a body force is meant to mimick the effect of classical no-slip boundary conditions at the physical surface of the obstacle. The present work deals with the (IBM) primarily proposed in the seminal work of Peskin [8], who introduced this method to simulate fluid-structure interactions into a cardio-vascular system (see the late, seminal paper by Peskin [9] for the mathematical foundation). The body force is evaluated on Lagrangian markers that can move over the Eulerian mesh to capture the motion or the deformation of the body. Approximations of the Delta distribution by smoother functions allow to interpolate between the two grids [8]. Since then, other formulations as well as several improvements and extensions of Peskin’s method have been proposed in the literature (see a review in Mittal and Iaccarino [10]). Among them, the direct forcing approach initially proposed by Uhlmann [11] and improved by Pinelli et al. [1] has been identified as the best candidate to be implemented in the OpenFOAM solver. Due to a new efficient quadrature for the spreading step and its capacity to deal with non-uniform and curvilinear meshes, this improved (IBM) offers modularity, stability, computational efficiency and accuracy, in particular in the analysis of moving/deformable configurations which are targeted in this work. OpenFOAM is now extensively used both in academic research (see among others the papers by Tabor and Baba-Ahmadi [12], Meldi et al. [13], Lysenko et al. [14], Komena and Shamsa [15]) and for industrial flows analysis (Selma et al. [16], Flores 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

et al. [17], Gao et al. [18]). It provides an efficient coding and a suitable environment for the implementation and the rapid dissemination of new algorithms to the users community, with the target to be used further to investigate real industrial problems. Several recent papers in the literature deal with the implementation of new numerical techniques or models in OpenFOAM (see among others the papers by Flores et al. [17], Towara et al. [19], Vuorinen et al. [20] ). In OpenFOAM, immersed bodies are primarily accounted by the use of wall-boundary conditions, although an (IBM) was recently implemented by Jasak et al. [21]. However, the (IBM) proposed by Pinelli et al. [1] proposes some improvements for the study of unsteady/deforming structures, as it relies only on the accuracy of the interpolation and spreading steps, which are independent of the complexity of the geometry. In addition, the implementation of this new (IBM) in the code gives us the opportunity to propose an original and careful verification and validation study, which is essential to provide reliable simulations and very useful for the OpenFOAM users community, in the absence of the quality certification and the lack of high-quality documentation and references compared to other commercial codes. Although it is not systematically mentioned explicitly in the literature, the application of direct forcing approaches in the context of incompressible flow solvers with predictor-corrector schemes is not straightforward. In fact, it is a two-constraints problem: on the one hand, the force term needed to impose the no-slip condition at the solid boundary must be calculated, and on the other hand, a divergence-free velocity at the boundaries must be satisfied. This means that enforcing divergence free conditions on the velocity affects the accuracy of the immersed boundary force at the wall. Although this issue has been claimed to be negligible by Fadlun et al. [22], it may actually lead to significant differences depending on the configuration considered. It has been shown to systematically introduce a first-order error in time on the actual boundary values [23]. A solution has been proposed by Ikeno and Kajishima [24] which changes the matrix structure of the Poisson problem solved to compute the value of the projector term (i.e., pressure or pressure correction), by directly imposing Neumann type conditions on the immersed boundary on the corresponding matrix terms. In order to avoid changing the matrix structure, Taira and Colonius [25] have suggested to use Lagrangian multipliers associated to boundary values to impose the expected velocity condition on the immersed boundary. Those Lagrangian multipliers are obtained solving a system derived from an algebraic splitting of the full spatial operator of the Navier-Stokes equations. In the present work, we choose an iterative scheme based on sub-iterations between (IBM) and pressure correction. This allows to use fast optimized Poisson solvers while keeping control of the error made on both the velocity at the immersed boundary and the divergence 4

ACCEPTED MANUSCRIPT

AN US

CR IP T

of the velocity field. The immersed boundary method has been widely used to solve problems with moving interface. In literature, the typical test-case used to check convergence is often the Taylor-Green vortex case. The drawback of this configuration is that it is difficult to evaluate precisely the convergence on variables of practical interest, such as the force coefficient or the no-slip condition on one boundary point. In this work, we propose an original method which circumvents this issue to evaluate convergence, by defining various errors at each step of the PISO solver. The objectives of the present work are threefold: (i) propose a modified PISO algorithm integrating the improved direct forcing approach of Pinelli et al. [1], (ii) propose an original verification technique allowing to estimate the errors related to the discretization and the (IBM), (iii) validate the solver on well-documented testcases of the literature. The paper is structured as follows: the numerical model is presented in section 2 including all amendments carried out. An original verification work based on the analysis of the convergence of various errors is performed in section 3. In section 4 the solver is validated comparing flow simulations around fixed and moving circular cylinders to available data of the literature. Finally, concluding remarks and perspectives are provided in section 5. 2. Numerical model

AC

CE

PT

ED

M

2.1. The geometrical and mathematical model Flows over bluff bodies are solved in a computational domain made of a parallelipedic volume of height H, width W and length Li + L0 , as shown in Figure 1. At the inlet, a steady uniform velocity is imposed along the streamwise direction x together with a zero pressure gradient. A mass conservation condition is imposed at the outlet. We assume periodic conditions in the spanwise direction z. Free-slip boundary conditions on the velocity are applied at the top and bottom of the domain. The incompressible flow of a viscous, Newtonian fluid is described by the dimensionless Navier Stokes equations written in a Cartesian frame of reference (x, y, z): 1 2 ∂u + ∇ · (uu) = −∇p + ∇u ∂t Re ∇ · u = 0,

(1) (2)

where u is the velocity vector, p is a normalized pressure by the fluid density ρ and Re = u∞ D/ν is the Reynolds number with ν the kinematic velocity, and D a characteristic length of the flow. The incompressibilty of the fluid is guaranteed by

5

ACCEPTED MANUSCRIPT

the continuity equation (2). Taking the divergence of equation (1) and using the continuity equation (2) yields classically a Poisson equation for the pressure : ∇2 p = −∇ · (u∇u)

(3)

AN US

CR IP T

Equations (1) and (3) couple pressure and velocity in an elliptic manner that requires specific numerical algorithms.

Figure 1 Computational domain. Example of a flow configuration around a cylinder of diameter D.

AC

CE

PT

ED

M

2.2. The numerical discretization The governing equations are discretized into OpenFOAM on fixed and blockstructured meshes composed of hexaedral elements [26]. Two neighboring blocks are connected by unstructured elements to ensure a better continuity between the blocks. The time discretization is based on the first-order implicit Euler scheme which has been chosen for its simplicity but could be staithforwardly extended to more accurate second-order OpenFOAM schemes, as the IBM implementation does not depend on the time discretization scheme. The velocity-pressure coupling is solved by the built-in solver pisoFoam in which the stabilizing extra term on the mass flux in the pressure correction loop is set to zero (see details in the recent paper by Vuorinen et al. [20]). The solver is thus Pressure Implicit with Splitting of Operators (PISO) algorithm described in the paper by Ferziger and Peric [27]. Three and one iterations were set for a PISO loop and for non orthogonal corrections, respectively [28]. Linear algebraic systems are solved using the Diagonal Incomplete LU Preconditioned Biconjugate Gradient DILUPBG (for the momentum equation (1)) and the Diagonal Incomplete Cholesky Preconditioned Conjugate Gradient DICPCG (for the Poisson equation (3)). For the present simulations involving low to moderate Reynolds numbers and regular block-structured meshes, no preconditionning has been used. For 6

ACCEPTED MANUSCRIPT

all independent variables, the required accuracy is 10−7 at each time step. For the numerical tests and the meshes presented in the paper, this solver has provided a good compromise between stability, accuracy and numerical cost.

CR IP T

2.3. The immersed boundary approach The predictor-corrector solver pisoFOAM has been modified to integrate the direct forcing approach proposed by Pinelli et al. [1]. The Navier-Stokes equations are discretized on a fixed mesh (Eulerian) while the solid boundary is discretized by a set of Lagrangian markers free to move over the Eulerian mesh, depending on the motion of the solid. In order to satisfy the two-constraints problem formed by no-slip and divergence-free conditions, we propose the new following procedure at each time step n :

AN US

1. Predictor step: b is obtained by solving the momentum Navier– (a) An estimate velocity u Stokes equations without any force term, and using the pressure p computed at the previous time step n − 1 :

ED

M

1 2 ∂b u b ) = −∇p + b + ∇ · (b uu ∇u (4) ∂t Re b (Eq. (b) The (IBM) force Fs is calculated on the Lagrangian markers using u (A.1)), and its values spread on the Eulerian mesh to calculate f (Eq. (A.4)), as detailed in Appendix A. (c) A new velocity u?,1 is calculated from the Navier-Stokes equations accounting now the immersed boundary force term f : ∂u?,1 1 2 ?,1 + ∇ · (u?,1 u?,1 ) = −∇p + ∇ u + f (b u) ∂t Re

(5)

CE

PT

u?,1 is the guess value of the velocity in the iterative PISO loop 2. PISO loop for m = 1 to M − 1, and up to convergence: (a) At each iteration, a pressure field p?,m is calculated from the following Poisson equation : ∇2 p?,m = −∇ · (u?,m ∇u?,m ) + ∇ · f (b u)

(6)

AC

(b) The velocity field is thus corrected using: u?,m+1 = g (u?,m , ∇p?,m , f (b u))

(7)

g as well as all discretized operators used in the algorithm are defined in Appendix B. 7

ACCEPTED MANUSCRIPT

3. Final step: the velocity and the pressure are updated : un+1 = u?,M −1 pn+1 = p?,M −1

CR IP T

(8) (9)

M

AN US

The discretization of the structure boundary leads to non-negligible errors due to the sharpness of the body force term (ideally discretized over 3 markers). This is illustrated in Figures 2 and 3, which show the interpolation of the IBM force term and its influence on the derivative of the force on one given Lagrangian marker.

PT

ED

Figure 2 Example of Eulerian discretization of the IBM force term f : comparison between 2nd order discretization using a centered scheme and the kernel analytical solution (left). Zoom on the mesh discretization (right). Flow simulation past a fixed cylinder at Re = 30.

AC

CE

Thus, the correct calculation of both the velocity divergence around the structure and the flux requires an appropriate integration of the body force term into the momentum equation. We propose here to compute the divergence of the momentum equation (in the PISO loop) and the interpolation of the fluxes by an hybrid calculation involving an analytical resolution (with the kernel function Eq. (A.3)) of the quantities involving the force term (singular quantities). Another solution would be to enlarge the stencil introducing additional points to derivate and interpolate the force term but this would lead to a more diffuse, and thus less accurate definition of the boundary. Numerical results in Figure 4 show the efficiency of this correction on the derivative, which can provide an error on the mass conservation per cell which is even smaller than the one obtained with a body fitted mesh (by a factor of 15%). 8

CR IP T

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

Figure 3 Derivative of the force term using three different kinds of interpolation. The results of a 2nd order discretization using a centered scheme are shown in empty grey symbols and refer to the left scale 10 times higher than the other results. The full grey symbols represent the new derivative calculated with the kernel function. The black line represents the derivative theoretical value for the flow simulation past a fixed cylinder at Re = 30.

AC

CE

Figure 4 Error on the mass conservation defined by ρij div(uij )Sij , with Sij the area of the 2D cell (i, j), for the flow around a 2D cylinder at Re = 30 using 312 Lagrangian markers. Classical no-slip boundary condition (left), IBM calculations without (center) and with (right) correction of the derivative. The divergence is calculated at the end of the PISO loop. The error magnitude varies from −10−5 (black) to 10−5 (white).

3. Solver verification A careful verification study is provided allowing to distinguish three different kind of errors coming from the discretization and the (IBM). This verification has been 9

ACCEPTED MANUSCRIPT

performed in 2D using a polynomial manufactured solution for the velocity and the pressure. Polynomial functions f (x, y), g(x, y) and h(x, y) have been chosen for the two components of the velocity ua (divergence free), and the pressure pa , respectively such that:

CR IP T

ua =

 2 2 2 2    f (x, y) = (1 − 0.01x ) (1 − 0.03y )(1 − 0.01y )

−0.02(1 − 0.01x2 )2 (y − 0.01)(y − 0.01y 3 )    g(x, y) = 0.5 + 0.04x(1 − 0.01x2 )(y − 0.01y 3 )(1 − 0.01y 2 ) pa = h(x, y) = f (x, y)g(x, y)

(10)

(11)

Different steps in the solver are verified according the definition of three errors:

AN US

• eFIBM is the error on the estimate of the (IBM) force term defined by equation (A.1) (during Step 2 of the IBM/PISO solver) and integrated on the body, hence computed as: X eFIBM =| (Fk − Fa )k | (12) k∈Dj

where :

Udk − Ua (13) ∆t and Ua is the value of the analytical solution ua , defined above by equation (10), evaluated on the Lagrangian markers.

ED

M

Fa =

AC

CE

PT

• enoslip is the error on the estimate of the no-slip condition at the boundary of the obstacle. This error is evaluated during the calculation of the IBM force term on the Eulerian mesh (end of Step 2 of the IBM/PISO solver). It is defined as the L∞ norm of the difference between the velocity on one Lagrangian marker (equation (A.2)), and the Eulerian velocity that has been spread and re-interpolated, i.e. : enoslip =k Us − I[S[Uk ]]s k∞

(14)

• eutot is the error on the velocity at the end of the PISO loop (Step 6 of the IBM/PISO solver). It is calculated in terms of both the L2 and L∞ norms: eutot /L2 =k u − ua k2 eutot /L∞ =k u − ua k∞ 10

(15) (16)

ACCEPTED MANUSCRIPT

The verification is made in five steps summarized below: 1. Computation of u and p according to: (17) (18)

CR IP T

∂u 1 2 + ∇ · (uu) = ∇p + ∇ u + Sa ∂t Re ∂ua 1 2 + ∇ · (ua ua ) − ∇pa − ∇ ua Sa = ∂t Re

AN US

2. Computation on the Lagrangian markers of the analytical values of the IBM force term Fa using ua and of the IBM force term Fs using the interpolated velocity u from the former step. 3. Calculation of eFIBM and enoslip using equations (12) and (14). 4. Spreading of the residual force Fs − Fa on the Eulerian grid. 5. Execution of steps 3 to 6 of the PISO algorithm and calculation of eutot /L2 and eutot /L∞ .

PT

ED

M

The errors eFIBM , enoslip and eutot have been calculated for two 2D flows past a circular and a square cylinder (to quantify the impact of the geometry on the method) of diameter and side L/5, respectively, L being the size of the computational domain. For a (10 × 10)-domain, four uniform grids have been tested, corresponding to ∆x = ∆y = 5 × 10−2 , 3.3 × 10−2 , 2.5 × 10−2 and 1.25 × 10−2 . To properly evaluate the interpolation error, the square cylinder is not aligned with the cells centers. Thus, Lagrangian markers have been chosen to be located between two Eulerian nodes. Results are shown on Figure 5. All errors decrease when the mesh is refined. The error enoslip exhibits a second-order rate of convergence whereas eutot /L2 and eutot /L∞ only exhibit a rate of convergence between 1 and 2 for both geometries. Finally, eFIBM exhibits a rate of convergence that depends on the geometry (as could have been expected), namely 1 for the square cylinder and nearly 2 for the circular cylinder.

CE

4. Solver validation

AC

The solver validation is performed using two- and three-dimensional (2D/3D) simulations of flows past a circular cylinder and a sphere of diameter D and at various Reynolds numbers (Re = U∞ D/ν) ranging between 30 and 3900. Such flows are of great interest since it constitutes a generic configuration with many applications in fields such as external aerodynamics, offshore engineering or environmental sciences. Vortex shedding is physically challenging because the separation point on the surface of the cylinder is not fixed by the geometry. Since the pressure field changes 11

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

CR IP T

rapidly near the separation and reatachment points, pressure prediction is therefore decisive for a correct estimation of the drag and lift coefficients. Numerically, in contrast to its square counterpart which can be tackled using Cartesian grids, this configuration requires curvilinear body-fitted (or unstructured) grids. The sphere is also a challenging configuration for (IBM) because the gemetry is intrinsically three-dimensionnal. All these features reveal these flows to be excellent test cases to evaluate our novel (IBM) solver fin OpenFOAM.

AC

Figure 5 Log-log plots of eFIBM , enoslip and eutot as a function of the mesh refinement, for flows past a circular (full thick line) and square (dashed thick line) cylinder. The dashed black and grey lines show the slopes of order 1 and 2, respectively.

12

ACCEPTED MANUSCRIPT

The Strouhal number, drag and lift coefficients are defined by: CD =

2Fd 2Fl Dfv , CL = 2 , St = , 2 ρu∞ D ρu∞ D u∞

(19)

CR IP T

where fv is the shedding frequency and Fd and Fl are the drag force and lift force per unit length, respectively, computed by integrating the immersed boundary force term in the Lagrangian space.

AC

CE

PT

ED

M

AN US

4.1. 2D and 3D flows around a cylinder at low Re 4.1.1. Computational details The center of the cylinder is the origin of the domain at (0, 0). The dimensions of the computational domain are those proposed by Pinelli et al. [1] and Vanella and Balaras [29], namely [−16D, 48D]×[−16D, 16D]×[−5.12D, 5.12D] in the streamwise (x), vertical (y) and spanwise (z) directions (Figure 6).

:

Figure 6 Computational domain decomposition and grid spacings. (x, y)-plane (top) and spanwise direction (bottom).

13

ACCEPTED MANUSCRIPT

CR IP T

The grid is uniform in the neighborhood of the cylinder, i.e. in the region −D ≤ x ≤ D and −D ≤ y ≤ D. For 3D computations, the 2D mesh has been extruded in the spanwise direction. Details pertaining to the resolution, as well as the number of Lagrangian markers and their relative spacing with respect to the Eulerian mesh are given in Table 1. Outside this region, the mesh size is stretched with a factor of 2.0 on five grid levels in the (x, y)-plane (as shown in Figure 6). Table 1 Mesh resolutions in the neighborhood of the cylinder: 2D cases 1 and 2 [−D, D] × [−D, D], 3D case 3 [−D, D] × [−D, D] × [−5.12D, 5.12D]. The α parameter defines the ratio of the distance between Lagrangian markers over the local Eulerian grid size [1].

Resolution

Lagrangian markers

α

1

∆x = ∆y = 0.02D

147

1.061

2

∆x = ∆y = 0.01D

312

1.004

3

∆x = ∆y = 0.02D, ∆z = 0.16D

AN US

Case

9792

1.004

M

All 2D and 3D simulations have been performed on 12 and 96 cpu, respectively. The CFL has been fixed to 0.5 and the number of PISO loop to 3. Simulations time varies from 24 hours (2D simulations) to 168 hours (3D simulations) depending on the mesh size and the flow regime on the middle size cluster of the Aix-Marseille computing facilities.

L/D

a/D

b/D

θo

CD

PT

ED

Table 2 Geometrical parameters of the wake and drag coefficients for the configuration of a fixed cylinder at Re = 30. L is the length of the recirculation, a is the distance between the cylinder and the recirculations centers, b is the vertical distance between the two recirculation centers, and θ is the separation angle measured from the rear stagnation point. Numerical and experimental data from literature are provided for comparison.

∆x = ∆y = 0.02D

1.66

0.556

0.53

47.80

1.78

∆x = ∆y = 0.01D

1.64

0.55

0.53

48.40

1.77

1.70

0.56

0.52

48.05

1.80

Blackburn and Henderson [30] (Num.)

-

-

-

-

1.74

Coutanceau and Bouard [31](Expe.)

1.55

0.54

0.54

50.00

-

Tritton [32](Expe.)

-

-

-

-

1.74

Present (Re = 30)

AC

CE

Pinelli et al. [1] (Num.)

4.1.2. 2D steady flow past a fixed cylinder At Re = 30, the flow is characterized by a steady recirculating region located just behind the cylinder. All characteristic geometrical parameters compare well with the 14

ACCEPTED MANUSCRIPT

data of the literature reported in Table 2, with differences less than 6% for the most refined grid.

M

AN US

CR IP T

4.1.3. 2D unsteady flows past a fixed cylinder Simulations in 2D unsteady regimes with vortex shedding have been performed at Re = 100 and 185, i.e. above Rec = 40 for the transition to unsteadiness according to Williamson [33] and Norberg [34]. The vorticity contours shown in Figure 7 for Re = 185 exhibit the well-known Karman vortex street featuring the periodic shedding of vortices, convected and diffused away from the cylinder. The topology of the solution compares well with that reported in several reference studies, see for instance the papers by Guilmineau and Queutey [35] and Pinelli et al. [1]. The corresponding time evolutions of CD and CL are plotted in Figure 8 and show that the amplitude of the lift and drag fluctuations increase with Reynolds number, in good agreement with the paper by Guilmineau and Queutey [35]. The Strouhal number, the mean drag (computed over 10 time periods) and the rms lift coefficients compare well with the literature data summarized in Table 3.

AC

CE

PT

ED

Figure 7 Vorticity countours evidencing the shedding of large-scale vortices in 2D flow past a fixed circular cylinder at Re = 185 (b). Vorticity magnitude varies from -1 (black) to 1 (white).

Figure 8 Temporal evolutions of CD (full line) and CL (dashed line) for the 2D flow past a fixed cylinder at Re = 100 (left); Re = 185 (right).

15

ACCEPTED MANUSCRIPT

Table 3 Mean drag, rms lift coefficients and Strouhal number for 2D flow past a fixed cylinder at Re = 100 and Re = 185. Numerical and experimental data from the literature are provided for comparison.

CLrms

St

o θmean

∆x = ∆y = 0.02D

1.38

-

0.165

118.9

∆x = ∆y = 0.01D

1.37

-

0.165

118.9

Blackburn and Henderson [30](Num.)

1.35

-

-

-

Williamson [36](Expe.)

-

-

0.164

-

Henderson [37](Num.)

1.35

-

-

-

∆x = ∆y = 0.02D

1.387

0.436

0.198

110.8

∆x = ∆y = 0.01D

1.379

0.427

0.198

110.8

∆x = ∆y = 0.005D

1.430

0.423

0.196

-

∆x = ∆y = 0.01D

1.509

0.428

0.199

-

1.377

0.461

-

-

Guilmineau and Queutey [35](Num.)

1.287

0.443

0.195

-

Williamson [33](Expe.)

-

-

0.193

-

Present (Re = 185) Pinelli et al. [1]

Vanella and Balaras [29](Num.)

AN US

Present (Re = 100)

CR IP T

CD

AC

CE

PT

ED

M

4.1.4. 2D unsteady flow past an oscillating cylinder In order to assess the capacity of the solver to encompass moving obstacles, a 2D simulation of flow past a sinusoidally moving cylinder is performed at Re = 500. Following Blackburn and Henderson [30], the cylinder is forced to oscillate in the vertical direction at a fixed amplitude ratio of A = 0.25 and with a frequency ratio of F (= fo /fv ) = 0.975 (with fo the frequency of the forced oscillation). The computational domain is the same as for the fixed simulation, with ∆x = ∆y = 0.02D around the cylinder. The shedding frequency fv is obtained from a preliminary flow simulation past a fixed cylinder at Re = 500. We start moving the cylinder once the flow has settled down to the established 2D shedding regime. A detailed description of the flow is provided in Figure 10. The five leftmost figures show vorticity contours and streamlines plotted at five instants spreading over half of the vortex shedding cycle, during which the lift force acts in the upwards direction, starting and ending at times of zero lift. Comparison with the results of Blackburn and Henderson [30] shown in the rightmost figures, provides good evidence that the spatial dynamics of the separation bubbles, but also the temporal evolution of the reattachment and separation points is very well predicted. Figure 9 shows the evolution of the lift coefficient as a function of the body displacement over the 10 last periods of oscillations. Again, the results are seen to match well the 16

ACCEPTED MANUSCRIPT

AN US

CR IP T

reference data of Blackburn and Henderson [30].

ED

M

Figure 9 Lift coefficient CL as a function of the cylinder displacement for the 2D flow past an oscillating cylinder at Re=500 : present results (grey line) vs. results obtained by Blackburn and Henderson [30] (black line).

AC

CE

PT

4.1.5. 3D unsteady flows past a fixed cylinder Additional simulations have been performed at Re = 200 and 300, i.e., above the critical value Rec = 190 for the transition to 3D flow, and within the range of Reynolds numbers where the 3D pattern transitions from mode A to mode B, according to the reference study of Williamson [36]. The present simulations predict well the occurrence of 3D vortex shedding (Figure 11). When increasing Reynolds number from Re = 200 to Re = 300, the solution shows a strong decrease of the spanwise wavelength λz , from λz /D ' 4.5 to λz /D ' 1.25 as previously observed by Williamson [36] at the transition between mode A and mode B. The temporal evolution of CD and CL in Figures 12 show a modulated behaviour characteristic of these 3D flows, all values being in agreement with the literature data, as seen from Table 4.

17

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 10 Instantaneous contours of vorticity for the 2D flow past an oscillating cylinder at Re=500: present results (left columns) vs. results obtained by Blackburn and Henderson [30] (right column) at five instants spreading over half of the shedding cycle. The vorticity magnitude is ranging from -1 (grey) to 1 (black).

18

CR IP T

ACCEPTED MANUSCRIPT

ED

M

AN US

Figure 11 Iso-surfaces of the instantaneous Q-criterion (−0.8 < Q < 0.8) at Re = 200 (left) and Re = 300 (right).

Figure 12 Temporal evolution of CD (full line) and CL (dashed line) for the 3D flow at Re = 200 (left); Re = 300 (right).

AC

CE

PT

4.2. Turbulent flow at moderate Re Envisioning application of the proposed IBM in OpenFOAM for the analysis of turbulent regimes, we consider here the flow around a circular cylinder in the subcritical regime at Reynolds number Re = 3900 (see in Breuer [38]). The von Karman vortex street at this Reynolds number already exhibits most of the characteristic features of industrial applications. Even though this test case is defined by a simple geometry, it is fully three-dimensional and unsteady, including transition regions to turbulence as well as flow separations along the sidewall. Therefore, it is identified as a relevant test case for the assessment of the IBM solver to perform simulations of complex turbulent flows, while being well-documented in the literature for validation.

19

ACCEPTED MANUSCRIPT

Table 4 Mean drag, rms lift coefficients and Strouhal number for 3D flow past a fixed cylinder at Re = 200 and Re = 300. Numerical and experimental data from the literature are provided for comparison.

CLrms

St

1.384

0.346

0.1802

Rajani et al. [39](Num.)

1.338

0.4216

0.1936

Qu et al. [40](Num.)

1.24

0.339

0.1801

-

-

0.1800

1.371

0.163

0.1915

1.43

0.453

0.198

1.28

0.499

0.195

1.26

0.38

0.203

-

-

0.203

-

0.435

0.203

1.22

-

-

Present (Re = 200)

∆x = ∆y = 0.02D &

∆Z = 0.16D

Williamson [36] (Expe.) Pinelli (private Communication) Present (Re = 300)

∆x = ∆y = 0.02D & ∆Z = 0.16D

Rajani et al. [39](Num.) Williamson [36](Expe.) Norberg [42](Expe.) Wieselsberger [43](Expe.)

AN US

Mittal and Balachandar [41](Num.)

CR IP T

CD

PT

ED

M

At this moderate Reynolds number, the flow is subcritical i.e., the boundary layers at the cylinder exhibit laminar separation and the transition takes place in the free shear layers. Therefore any DES model, which work in RANS mode in the near wall region, is expected to provide reliable representation. Here, the Delayed Detached Eddy Simulation with the OpenFOAM native implementation is used. It is formulated for the Spalart Almaras model [44]. Its implementation in the IBM framework simply requires to evaluate the distance to the wall (distance between the cells center and the Lagrangian points) for each cell of the domain.

AC

CE

4.2.1. Computational details The center of the cylinder is the origin of the domain at (0, 0). The dimensions of the computational domain are [5D, 15D] × [−10D, 10D] × [−1.57D, 1.57D] in the streamwise (x), vertical (y) and spanwise (z) directions, Figure 13. The flow periodicity was assumed to be in the spanwise direction. Inflow, outflow, upper and lower sides of the domain are defined as shown Figure 1.

20

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

Figure 13 Computational domain decomposition and grid spacings (with k the cell-to-cell stretching ratio). (x, y)-plane (top) and spanwise direction (bottom).

AC

CE

PT

A rather coarse mesh of 4 millions grid cells has been considered (most of the CFD references of the literature consider between 5 to 10 millions of grid cells) in order to assess the capabilities of the code in capturing the main coherent structures developing in the wake. The domain has been discretized in 15 elements in the spanwise direction, leading to very elongated computational cells in the z-direction. No wall function is used, so the mesh is refined near the cylinder wall (y + . 4 in wall units) and in the wake, as shown on Figure 13. The time-step equal to 10−3 D/u∞ is chosen so as to satisfy a CFL condition of 0.3 here to improve the accuracy. Computations have been run over a rather long time interval of 460D/u∞ that corresponds to about 90 shedding periods to achieve meaningful statistics. The simulation takes14 days on 96 cpu of the regular custer of the Aix-Marseille University computing facilities. 21

ACCEPTED MANUSCRIPT

Table 5 Rms lift, mean drag coefficients, Strouhal number and mean separation angle for 3D flow past a fixed cylinder at Re = 3900. Numerical and experimental data from the literature are provided for comparison.

CDmean

St

o θmean

Present (Re = 3900)

0.2593

1.0389

0.2108

89.89

Parnaudeau et al. [45] PIV

-

-

0.21

-

Lourenco and Shi [46] PIV

-

0.99

0.21

-

D’Alessandro et al. [47] [OpenFOAM SA-IDDES]

0.1458

1.0235

0.222

87

D’Alessandro et al. [47] [OpenFOAM NLDES]

0.3832

1.1751

0.217

88.99

D’Alessandro et al. [47] [OpenFOAM SA-DES]

0.4248

1.2025

0.215

89.28

D’Alessandro et al. [47] [OpenFOAM v2-f DES]

0.1088

0.9857

0.214

86.40

Kravchenko and Moin [48] LES

-

1.04

0.21

88

AN US

CR IP T

CLrms

AC

CE

PT

ED

M

4.2.2. DDES results The analysis of the main flow features indicates that a reasonably high level of precision in the flow prediction has been achieved. Although being rapidly dissipated in the wake due to the stretching of the mesh , the computations exhibit the first vortices of the well-known von Karman vortex street with periodic vortex shedding as shown on Figure 14. The picture clearly shows the transition to turbulence takes place in the free shear layers in the very near wake.

Figure 14 Snapshot of instantaneous contour of the vorticity magnitude ranging from -10 to 10 at Re = 3900. IBM DDES (left), Body fitted DDES (right).

Typical time histories of the computed lift coefficient CL and drag coefficient CD 22

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 15 Lift CL and drag CD coefficients time histories at Re = 3900.

AC

CE

PT

ED

M

are plotted on Figure 15 and show cyclic oscillations due to the vortex shedding phenomenon and high-frequency turbulent fluctuations. The mean drag and rms lift coefficients are in agreement with numerical and experimental data available in the literature, as shown on Table 5. First and second-order flow statistics are also shown on Figures 16, together with experimental data and other OpenFOAM results of the literature. The mean streamwise velocity profiles (Figure 16) are rather well predicted, both in the wake centerline and at three locations in the wake. The agreement is particuarly satisfactory with the PIV measurements of Lourenco and Shi [46], with a good estimate of the mean streamwise length of the recirculation zone (corresponding to negative velocity values) and of the U-shape of the profiles. It is noteworthy to notice that the present results match well the body-fitted results obtained with the same turbulence model meaning that the IBM correctly mimicks the solid boundary. The same overall agreement is obtained on the streamwise Reynolds stress in the near wake, Figure 16. The overestimate of the streamwise velocity fluctuations at x = 1.06D and around y = 0 , and so of the turbulence intensity, should be related to an early transition to turbulence in the thin shear layers behind the cylinder leading to their shortening. This seems supported by the snapshot of vorticity on Figure 14.

23

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

.

Figure 16 3D flow past a fixed cylinder at Re = 3900. Mean streamwise velocity profiles in the wake centerline (top) and at three locations in the near wake (middle). Resolved streamwise Reynolds stresses at threet locations in the wake24 (bottom). (–) Present simulation IBM-DDES, (–) Present simulation DDES Body fitted, () PIV [45], (4)PIV [46], (- -)OpenFOAM DES [47], (- -) OpenFOAM IDDES [47]

ACCEPTED MANUSCRIPT

CD =

CR IP T

4.3. Flow around a sphere The performance of the IBM solver to accurately represent three dimensional configurations is assessed by the analysis of the laminar flow around a sphere of diameter D. Different Reynolds numbers in the range [100, 300] are considered. Because of the increased dimensionnality with respect to the cylinder configuration, the flow past a sphere exhibits the emergence of more complex dynamics and vortex interactions. The drag, lift and side coefficients are defined as: 4Fd 4Fl 4Fs , CL = 2 , CS = 2 , 2 2 2 ρu∞ πD ρu∞ πD ρu∞ πD2

(20)

AN US

where Fd , Fl and Fs are the drag force, lift and side force per unit length, respectively, computed by integrating the immersed boundary force term in the Lagrangian space.

CE

PT

ED

M

4.3.1. Computational details The center of the sphere is the origin of the domain at (0, 0, 0). The dimensions of the computational domain are those used above for the cylinder, namely [−16D, 48D] × [−16D, 16D] × [−H/2, H/2] in the streamwise (x), vertical (y) and spanwise (z) directions as shown in Figure 17. To evaluate the effect of the domain size in the spanwise direction, two values of H have been considered, H = 10 nd H = 32, denoted hereafter (dom1) and (dom2), respectively. The grid is uniform in the neighborhood of the sphere, i.e. in the region −1.2D ≤ x ≤ 2D, −1.2D ≤ y ≤ 1.2D and −1.2D ≤ z ≤ 1.2D. The body is represented via 7652 lagrangian markers (α = 1.012). Outside this region, the mesh size is coarsened with a factor of 2.0 on four grid levels in the (x, y)-plane and (x, z)-plane (as shown Figure 17) with ∆x = ∆y = ∆z = ∆. As for the cylinder, all simulations have been performed on 96 cpu of the AMU computing facilities. The CFL has been fixed to 0.2 and the number of PISO loops to 3.

AC

4.3.2. 3D steady axisymmetric flow For Reynolds numbers between 20 and approximately 210, the flow is separated, steady, axisymmetric, and topologically similar [49]. The flow is computed here at Re = 100. It is characterized by a steady axisymmetric recirculating region located just behind the sphere. All characteristic geometrical parameters compare well with the data of the literature, Table 6. Differences lower than 6% are observed.

25

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 17 Computational domain decomposition and grid spacings. (x, y)-plane and (x, z)-plane

M

Table 6 Geometrical parameters of the wake and drag coefficient for the flow past a sphere at Re = 100. L is the length of the recirculation and θ is the separation angle measured from the rear stagnation point. Numerical and experimental data from literature are provided for comparison.

L/D

θo

CD

0.92

53.03

1.14

Taneda [50] (Expe.)

0.89

-

-

Nakamura [51] (Expe.)

-

53

-

Johnson and Patel [49] (Num.) ∆x = ∆y = 0.005D

0.88

53

1.08

Giacobello et al. [52] (Num.)

0.88

53

-

PT

ED

Present (Re = 100) ∆x = ∆y = ∆z = 0.02D

AC

CE

4.3.3. 3D steady non-axisymmetric flow For Reynolds numbers in the range [211, 270], the flow no longer exhibits axial symmetry [49] but however, remain steady. Despite being non-axisymmetric, the flow does contain a plane of symmetry. Computations are performed here at Re = 250. The location of the symmetry plane was allowed to arise naturally, emerging by the numerical perturbation of the solver only. For a clear presentation of the results, the flow field has been rotated such that the symmetry plane coincides with the (x, y)-plane. Results obtained in the large domain in the spanwise direction (dom2) 26

ACCEPTED MANUSCRIPT

(Table 7) are in better agreement with reference data, with differences less than 3%. Table 7 Drag and lift mean coefficients for the flow past a fixed sphere at Re = 250. Data (numerical) from the literature are provided for comparison.

∆x = ∆y = ∆z = 0.02D − (dom1) ∆x = ∆y = ∆z = 0.02D − (dom2)

Johnson and Patel [49] (Num.) ∆x = ∆y = 0.005D Giacobello et al. [52] (Num.)

CLmean

0.76

-0.057

CR IP T

Present (Re = 250)

CDmean 0.72

-0.062

0.70

-0.061

0.702

-0.061

AC

CE

PT

ED

M

AN US

The presence of a plane of symmetry in the flow can be clearly observed on Figures 18 and 19 showing the three-dimensional particle path out of the (x,y)-plane upstream of the sphere and the snapshot of 3D streamwise, respectively. Similarly to the simulation of Johnson and Patel [49], it can be seen that the upper spiral in the (x,y)-plan is actually fed by fluid from upstream, while the lower spiral releases fluid into the wake after sending it up and around the upper spiral.

Figure 18 Snapshot of 3D particle path for (x,y)-view (a), (x,z)-view (b) and (y,z)-view (c). Present results (left) and results of Johnson and Patel [49] (right). Flow past a sphere at Re = 250.

Present computations predict also well the pressure coefficient distribution as 27

ACCEPTED MANUSCRIPT

AN US

CR IP T

shown by the comparions with numerical results of Johnson and Patel [49] on Figure 20. The pressure field in the (x,z)-plane is completely symmetric but the pressure contours in the (x, y)-plane are not symmetric. The pressure minimum in the region of the lower vortex is lower than that in the region of the upper vortex which correspond to the phenomenon observed by Johnson and Patel [49].

AC

CE

PT

ED

M

Figure 19 Snapshot of 3D streamwise for (x,y)-view (top), (x,z)-view (bottom) compare to the photographs of the dye-filled of Johnson and Patel [49]. Present results (a), numerical (b) and experimental (c) results of Johnson and Patel [49]. Flow past a sphere at Re = 250

Figure 20 Snapshot of instantaneous flow field pressure coefficient contours for a non-axisymmetric flow. Flow past a sphere at Re = 250. Present results (left) and results of Johnson and Patel [49] (right).

28

ACCEPTED MANUSCRIPT

CR IP T

4.3.4. 3D unsteady non-axisymmetric flow For Reynolds numbers larger than Re = 270 ,the flow past a sphere is expected to become unsteady [49]. The flow shows thus a highly organized periodic structure dominated by vortex shedding. Computations are performed here at Re = 300. Present results are summarized in Table 8 together with literature data. As for the 3D steady case, results obtained in the larger (spanwise) domain (dom 2) are in closer agreement with the data of the literature. The agreement between numerical results is very good, with differences less than 4% for all quantities. Regarding experimental results, only few measurements are avalaible, and differences with present results are of about 6% and 8% for the St and the CDmean , respectively.

AN US

Table 8 Drag and lift mean coefficient for the flow past a sphere at Re = 300. Numerical and experimental data from literature are provided for comparison.

CDmean

Present (Re = 300)

∆x = ∆y = ∆z = 0.02D − (dom1) ∆x = ∆y = ∆z = 0.02D − (dom2)

Johnson and Patel [49] (Num.) Giacobello et al. [52] (Num.) Tomboulides et al. [53] (Num.)

10−3

0.656 ± 3.5 × 10−3

St 10−2

0.13

10−2

0.139

0.0659 ± 1.89 × 0.066 ± 2.03 ×

0.136

0.067

0.134

-

0.136

0.629

-

-

-

-

0.148-0.165

0.658

10−3

AC

CE

PT

Johnson and Patel [49] (Expe.)

CLmean

0.069 ± 1.6 × 10−2

0.671 ± 2.8 ×

ED

Roos and Willmarth [54] (Expe.)

0.679 ± 3.9 ×

M

∆x = ∆y = 0.005D

0.705 ± 3.3 ×

10−3

Figure 21 Average streamwise velocity (left) and r.m.s (right). Present result (–), Johnson and Patel [49] (- -). Flow past a sphere at Re = 300.

29

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 22 Streamwise velocity at every quarter period along the axis of motion from the rear of the sphere. φ = 0 (–), φ = π/2 (- -), φ = π (-.-), φ = 3π/2 (...). Present results (up), Johnson and Patel [49] (down). Flow past a sphere at Re = 300.

CE

PT

ED

M

The near wake dynamics (x ≤ 5) is well predicted as shown by Figure 21 both for the mean streamwise velocity and r.m.s. quantities. For (x > 5) the mean streamwise velocity remains well predicted while the r.m.s. becomes overestimated 15% probably due to the coarsening of the mesh in the far wake. Figure 22 shows the streamwise velocity at every quarter period from the rear of the sphere. We can observe the same convecting structures in the wake as for the work of Johnson and Patel [49]. A small delay is observed, which is probably due to the mesh refinement in the wake. The traveling wave with a peak at x = 5 instead of xref = 4.5 for φ = 0 move to x = 9 instead of xref = 8.5 for φ = 3π/2. The location of zero velocity is well captured for every quarter of phase going from x = 2 to x = 1.5 for φ = π/2 to φ = 3π/2 respectively. 5. Concluding remarks

AC

This study has proposed a new immersed boundary method in OpenFOAM to simulate incompressible bluff body fluid flows. This IBM, originally proposed by Pinelli et al. [1], is accurate and versatile for the study of unsteady/deforming structures, as it relies only on the accuracy of the interpolation and spreading steps, which are independent of the complexity of the geometry. The paper details the integration 30

ACCEPTED MANUSCRIPT

AN US

CR IP T

of the method in the finite-volume formalism and the changes in the native PISO algorithm to satisfy the two-constraints problem involved by the imposition of both the no-slip and the divergence free conditions on the velocity at the solid boundary. A careful and original verification study has been provided using a manufactured solution, which may be applied in a more general context of algorithms using IBM. The efficiency and the accuracy of the algorithm has been studied on various 2D and 3D simulations of flows around fixed and moving cylinders, as well as around a fixed sphere, and for Reynolds numbers ranging from 30 to 3900, i.e. from 2D steady to turbulent regimes. Careful comparisons with available numerical and experimental results of the literature have shown very good agreements. Analysis of the computational cost, numerical behavior and accuracy of the numerical method show that the global properties of the OpenFOAM solver are not alterated. This work, based on a careful verification and validation study, offers new capabilities to OpenFOAM to perform reliable simulations.

AC

CE

PT

ED

M

Acknowledgements This work was granted access to the HPC resources of AixMarseille University nanced by the project Equip@Meso (ANR-10-EQPX-29-01) and of IDRIS, under the allocations i20170242 made by GENCI. This work has been carried out thanks to the support of the FUI 15 URABAILA granted in the frame of the Energy Climate Program of the French government. The financial support of the European Commission through the PELskin FP7 European project (AAT.2012.6.31-Breakthrough and emerging technologies) is greatly acknowledged. EC thanks Aix-Marseille University and Direction G´en´erale de l’Armement (DGA) for his PhD grant.

31

ACCEPTED MANUSCRIPT

References [1] A. Pinelli, I. Naqavi, U. Piomelli, J. Favier, Immersed-boundary methods for general finite-difference and finite-volume Navier–Stokes solvers, J. of Comp. Phys. 229 (24) (2010) 9073 – 9091.

CR IP T

[2] ESI-OpenCFD, http://wwww.openfoam.com, OpenFoam . [3] R. Glowinski, T.-W. Pan, T. Hesla, A distributed Lagrange multiplier/fictitious domain method for particulate fows, Int. J. of Multiphase Flow 25 (1999) 755– 794.

AN US

[4] Y. Cheny, O. Botella, Set Method for the Computation of Incompressible Viscous Flows in Complex Moving Geometries with Good Conservation Properties, J. of Comp. Phys. 229 (2010) 1043–1076. [5] C. S. Peskin, Flow patterns around heart valves: A numerical method, J. of Comp. Phys. 10 (2) (1972) 252 – 271. [6] M. Minguez, R. Pasquetti, E. Serre, High-order large-eddy simulation of flow over the Ahmed body car model, Phys. Fluids 20.

M

[7] L. Isoardi, G. Chiavassa, G. Ciraolo, Penalization modeling of a limiter in the Tokamak edge plasma, J. of Comp. Phys. 229 (2010) 2220–2235.

ED

[8] C. S. Peskin, Numerical analysis of blood flow in the heart, J. of Comp. Phys. 25 (3) (1977) 220 – 252. [9] C. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 1 – 39.

PT

[10] R. Mittal, G. Iaccarino, Immersed boundary methods, Ann. Rev. of Fluid Mech. 37 (2005) 239 – 261.

CE

[11] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. of Comp. Phys. 209 (2) (2005) 448 – 476.

AC

[12] G. Tabor, M. Baba-Ahmadi, Inlet conditions for large eddy simulation: A review, Comp. & Fluids 39 (4) (2010) 553 – 567. [13] M. Meldi, M. Salvetti, P. Sagaut, Quantification of errors in large-eddy simulations of a spatially evolving mixing layer using polynomial chaos, Phys. Fluids 24 (3) (2012) 035101. 32

ACCEPTED MANUSCRIPT

[14] D. Lysenko, I. Ertesvg, K. Rian, Modeling of turbulent separated flows using OpenFOAM, Comp. & Fluids 80 (2013) 408–422. [15] E. Komena, A. Shamsa, Quasi-DNS capabilities of OpenFOAM for different mesh types, Comp. & Fluids 96 (2014) 87–104.

CR IP T

[16] B. Selma, M. Dsilets, P. Proulx, Optimization of an industrial heat exchanger using an open-source CFD code, Applied Thermal Eng. 69 (12) (2014) 241 – 250. [17] F. Flores, R. Garreaud, M. R.C., CFD simulations of turbulent buoyant atmospheric flows over complex geometry: Solver development in OpenFOAM, Comp. & Fluids 82 (2013) 1–13.

AN US

[18] L. Gao, J. Xu, G. Gao, Numerical Simulation of Turbulent Flow past Airfoils on OpenFOAM, Comp. & Fluids 31 (2012) 756–761. [19] M. Towara, M. Schanen, U. Naumann, MPI-Parallel Discrete Adjoint OpenFOAM, Comp. & Fluids 51 (2015) 19–28.

M

[20] V. Vuorinen, J. Keskinen, D. C., On the implementation of low-dissipative Runge-Kutta projection methods for time dependent flows using OpenFOAM, Comp. & Fluids 93 (2014) 153–163.

ED

[21] H. Jasak, D. Rigler, Z. Tukovic, Design and implementation of Immersed Boundary Method with discrete forcing approach for boundary conditions, In proceedings of 6th European Congress on Computational Fluid Dynamics - ECFD VI Barcelona, Spain ISBN: 978-849428447-2.

PT

[22] E. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined ImmersedBoundary Finite-Difference Methods for Three-Dimensional Complex Flow Simulations, J. of Comp. Phys. 161 (1) (2000) 35 – 60.

AC

CE

[23] F. Domenichini, On the consistency of the direct forcing method in the fractional step solution of the Navier-Stokes equations, J. Comp. Phys. 227 (12) (2008) 6372–6384. [24] T. Ikeno, T. Kajishima, Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations, J. of Comp. Phys. 226 (2) (2007) 1485–1508.

33

ACCEPTED MANUSCRIPT

[25] K. Taira, T. Colonius, The immersed boundary method: A projection approach, J. Comp. Phys. 225 (10) (2007) 2118–2137. [26] H. Jasak, Error analysis and estimation for the Finite Volume method with applications to fluid flows, PhD. Thesis, Imperial College, University of London .

CR IP T

[27] J. Ferziger, M. Peric, Computational Methods in Fluid Dynamics, New-York : Springer-Verla . [28] E. Devilliers, The potential of large eddy simulation for the modeling of wall bounded flows, PhD thesis, Imperial college of science, technology and medicine .

AN US

[29] M. Vanella, E. Balaras, A moving-least-squares reconstruction for embeddedboundary formulations, J. Comp. Phys. 228 (18) (2009) 6617 – 6628. [30] H. Blackburn, R. Henderson, A study of two-dimensional flow past an oscillating cylinder, J. Fluid Mech. 385 (1999) 255–286.

M

[31] M. Coutanceau, R. Bouard, Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow, J. Fluid Mech 79 (1977) 231–256.

ED

[32] D. Tritton, Experiments on the flow past a circular cylinder at low Reynolds numbers, J. Fluid Mech 6 (1959) 547 – 567.

PT

[33] C. Williamson, Defining a universal and continuous StrouhalReynolds number relationship for the laminar vortex shedding of a circular cylinder, Phys. Fluids 31 (10) (1988) 2742–2744.

CE

[34] C. Norberg, An experimental investigation of the flow around a circular cylinder: influence of aspect ratio, J. Fluid Mech 258 (1994) 287 – 316. [35] E. Guilmineau, P. Queutey, A numerical simulation of vortex shedding from an oscillating circular cylinder, J. Fluids Struct. 16 (6) (2002) 773 – 794.

AC

[36] C. Williamson, Vortex dynamics in the cylinder wake, Annu. Rev. Fluid. Mech. 28 (1996) 477 – 539. [37] R. Henderson, Details of the drag curve near the onset of vortex shedding, Phys.Fluids 7 (9) (1995) 2102 – 2104. 34

ACCEPTED MANUSCRIPT

[38] M. Breuer, Large eddy simulation of the subcritical flow past a circular cylinder : numerical and modeling aspects, Phys. Fluids 12 (2) (2000) 403–417. [39] B. Rajani, A. Kandasamy, M. Sekhar, Numerical simulation of laminar flow past a circular cylinder, Appl. Math. Modelling 33 (2009) 1228–1247.

CR IP T

[40] L. Qu, C. Norberg, L. Davidson, S. Peng, F. Wang, Quantitative numerical analysis of flow past a circular cylinder at Reynolds number between 50 and 200, J. Fluids Struc. 39 (2013) 347–370. [41] R. Mittal, S. Balachandar, Generation of streamwise vortical structures in bluff body wakes, Phys. Rev. Letter 75 (1995) 1300–1303.

AN US

[42] C. Norberg, Pressure forces on a circular cylinder in cross flow, Bluff-Body Wakes, Dynamics and Instabilities (1993) 275–278. [43] C. Wieselsberger, New data on the law of hydro and aerodynamic resistance, NACA TN 84. [44] P. Spalart, S. Deck, M. Shur, K. Squires, M. Strelets, A. Travin, A new version of detached-eddy simulation, resistant to ambiguous grid densities, Theor. Comp. Fluid Dyn. 20 (2006) 181 – 195.

ED

M

[45] P. Parnaudeau, J. Carlier, D. Heitz, E. Lamballais, Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900, Phys. Fluids 20 (6) (2008) 441–453. [46] L. Lourenco, C. Shi, Characteristics of the plane turbulent near wake of a circular cylinder, a particle image velocimetry., 1993.

CE

PT

[47] V. D’Alessandro, S. Montelpareb, R. Riccia, Detachededdy simulations of the flow over a cylinder at Re = 3900 using OpenFOAM, Comp. & Fluids 136 (2016) 152–169. [48] A. Kravchenko, P. Moin, Numerical studies of flows over a circular cylinder at ReD = 3900, Phys. Fluids 12 (2) (2000) 403–417.

AC

[49] T. Johnson, V. Patel, Flow past a sphere up to a Reynolds number of 300, J.Fluid Mech. 378 (1998) 19–70. [50] S. Taneda, Experimental investigation of the wake behind a sphere at low Reynolds numbers, J. Phys. Soc. Japan 11 (1956) 1104. 35

ACCEPTED MANUSCRIPT

[51] I. Nakamura, Steady wake behind a sphere, Phys Fluids 19 (1976) 5. [52] M. Giacobello, A. Ooi, S. Balachandar, Wake structure of a transversely rotating sphere at moderate Reynolds numbers, J. Fluid Mech. 621 (2009) 103–130.

CR IP T

[53] A. Tomboulides, S. Orszag, G. Karniadakis, Direct and large-eddy simulation of axisymmetruc wakes, AIAA Paper 93 (1993) 0546. [54] F. W. . Roos, W. W. Willmarth, Some experimental results on sphere and disk drag, AIAA J. (4) (1971) 285291. [55] Z. Li, J. Favier, U. D’Ortona, S. Poncet, A numerical approach to combine immersed boundary method and lattice Boltzmann model for single- and multicomponent fluid flows, J. Comp. Phys. in press.

AC

CE

PT

ED

M

AN US

[56] A. Roma, C. Peskin, M. Berger, An adaptive version of the immersed boundary method, J. Comp. Phys. 153 (1999) 509 – 534.

36

ACCEPTED MANUSCRIPT

Appendix A. Calculation of the body force terms in the (IBM)

CR IP T

The Navier–Stokes equations are discretized on a fixed mesh (Eulerian) while the solid boundary is discretized by a set of Lagrangian markers free to move over the Eulerian mesh, depending on the motion of the solid. As in traditional direct forcing methods, the target velocity U d is directly imposed at the boundary nodes. This velocity is equal to either the local fluid velocity or zero depending on wether the solid moves or is at rest. Appendix A.0.1. Calculation of the body force term on the Lagrangian markers: the interpolation step

AN US

The body force is computed in the Lagrangian space, i.e. at all Lagrangian markers. On the sth Lagrangian marker and at time step (n + 1), the force term Fn+1 , is given s by: u]s Uds − I[ˆ (A.1) ∆t where Uds is the target velocity on the sth Lagrangian marker. I[ˆ u]s stands for the interpolation on the sth Lagrangian marker of the fluid velocity known on the Eulerian mesh at time step n, and computed without any force term. As presented in Li et al. [55], the discrete expression of the interpolation operator is given by : X I[un ]s = unj δh (xj − Xs )∆v (A.2)

M

Fn+1 = s

ED

j∈Ds

AC

CE

PT

where the j-index refers to the discrete value of the fluid velocity on the Eulerian mesh, Xs refers to the coordinates of the sth Lagrangian marker and ∆v refers formally to an Eulerian quadrature, i.e. ∆v = ∆x∆y∆z for the case of a Cartesian uniform mesh. The interpolation kernel δh is the discretized delta function used in Roma et al. [56] :    √ 1  2+1  −3r 1 + 0 ≤ r ≤ 0.5   3 h i (A.3) δh (r) 1 5 − 3r − p−3(1 − r)2 + 1 0.5 ≤ r ≤ 1.5   6   0 otherwise

It is centered on each Lagrangian marker s and takes non-zero values inside a finite domain Ds , called the support of the sth Lagrangian marker. This δh can be calibrated as proposed in Pinelli et al. [1] to take into account the non uniform meshes. 37

ACCEPTED MANUSCRIPT

Appendix A.0.2. Calculation of the body force term on the Eulerian mesh: the spreading step

CR IP T

Once the force term is computed from equation (A.1), one needs to transfer its value to the Eulerian mesh. This is done by the spreading step, which is the inverse operation of the interpolation. The value of the force term evaluated on the Eulerian mesh, f n+1 (xj ), is given by: X S[Fn+1 ] = f n+1 (xj ) = Fn+1 δh (xj − Xk )k (A.4) k k k∈Dj

AN US

The k-index refers to a loop over the Lagrangian markers whose support contains the Eulerian node j. k is the Lagrangian quadrature, which is calculated solving a linear system : A = 1 (A.5) where the vectors  = (1 , . . . , Ns )T and 1 = (1, . . . , 1)T have a dimension of Ns , Ns being the number of Lagrangian markers, and A is the matrix defined by the product between the k th and the lth interpolation kernels such that: X Akl = δh (xj − Xk )δh (xj − Xl ) (A.6)

M

j∈Dl

ED

For further details, see the papers by Pinelli et al. [1], Li et al. [55]. Appendix B. Discretized operators in the modified PISO Loop At each PISO loop iteration m, let’s define [S] such that:

PT

1 2 ?,m ∂u?,m + ∇ · (u?,m u?,m ) − ∇u ∂t Re with A its diagonal part and H its extra-diagonal part, i.e.:

CE

[S] =

[S] = [A] − [H]

(B.1)

(B.2)

AC

with :



 [A] =  

ai−1;j−1 u?,m i−1;j−1

0

0

0

ai;j u?,m i;j

0

0

ai+1;j+1 u?,m i+1;j+1

0

38



  = {aii } × [U ?,m ] ii 

(B.3)

ACCEPTED MANUSCRIPT

and 

 [H] =  

?,m ai;j−1 u?,m i;j−1 ai+1;j−1 ui+1;j−1

0 ai−1;j u?,m i−1;j

0

ai+1;j u?,m i+1;j

?,m ai−1;j+1 u?,m i−1;j+1 ai;j+1 ui;j+1

0

   

(B.4)

[A] − [H] = −∇p?,m + f (b u)

CR IP T

Taking into account the (IBM) force term spread on the Eulerian mesh f (b u) and the unknown pressure p?,m , we get :

AN US

[A] = {aii } × [Uii?,m ] = −∇p?,m + f (b u) + [H] ?,m [Uii?,m ] = {a−1 + f (b u) + [H]) ii }(−∇p

(B.5) (B.6)

(B.7)

Taking the divergence of equation (B.7) and using the discrete continuity condition ∇ · [Uii?,m ] = 0 we get the Poisson equation for the pressure: (B.8)

M

?,m ∇ · ({a−1 ) = ∇ · ({a−1 u)) + ∇ · ({a−1 ii }∇p ii }f (b ii }[H])

ED

The term ∇ · {a−1 u) corresponds to the divergence of the force term. The ii }f (b velocity is then corrected, using the new pressure p?,m such that:

with

PT

[Uii?,m+1 ] = g(u?,m , p?,m , f (b u))

g ≡ {aii−1 }(−∇p?,m + f (b u) + [H])

CE

In order to update the flux, the force should be interpolated on the surface, introducing the issue discussed in section 2.2. Thus, the equation is written as: F = S · [{a−1 u) + [H])]f aces ii }(−∇p + f (b

(B.9)

AC

with f (b u)f aces the immersed boundary force is calculated analytically on the surface.

39