A sharp-interface immersed boundary method for moving objects in compressible viscous flows

A sharp-interface immersed boundary method for moving objects in compressible viscous flows

Journal Pre-proof

A sharp-interface immersed boundary method for moving objects in compressible viscous flows Francesco De Vanna, Francesco Picano , Ernesto Benini PII: DOI: Reference:

S0045-7930(19)30373-1 https://doi.org/10.1016/j.compfluid.2019.104415 CAF 104415

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

31 July 2019 3 December 2019 24 December 2019

Please cite this article as: Francesco De Vanna, Francesco Picano , Ernesto Benini , A sharp-interface immersed boundary method for moving objects in compressible viscous flows, Computers and Fluids (2020), doi: https://doi.org/10.1016/j.compfluid.2019.104415

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Elsevier Ltd. All rights reserved.

Highlights • Numerically methodology for moving objects in compressible flows. • Boundary oscillations can be reduced by energy-preserving schemes. • Numerical benchmarks for moving objects in compressible flows.

1

A sharp-interface immersed boundary method for moving objects in compressible viscous flows Francesco De Vanna1,, Francesco Picano1 , Ernesto Benini1 5

Department of Industrial Engineering, University of Padova, Via Venezia 1, 35131 Padova, Italy

Abstract

10

15

20

25

Sharp-interface Immersed Boundary Methods for moving solids in compressible viscous flows often exhibit spurious noise propagating from the moving boundary. In compressible flows, filtering or upwinding discretization techniques are often used to overcome the problem. In the present work, we show that using a conservative energy-preserving finite difference method for convection in combination with a Ghost-Point-Forcing-Method (GPFM) is able to keep controlled the pressure-velocity spurious noise. In order to deal with compressible flow in a wide range of Mach numbers the shock-dynamics appearing in high-speed flows conditions has been addressed hybridising the latter scheme with a fifth-order weighted-essentially-non-oscillatory (WENO) scheme. The latter, in the path of keeping the numerical dissipation minimal, has been confined around the shock locations using a proper detector. In this respect, the method appears a suitable alternative for direct and large-eddy simulation of moving objects in compressible flows. The entire methodology was found to be robust ranging from weakly compressible to highly supersonic flows including shocks. In order to prove the cleanliness and the robustness of the method, several well-documented benchmarks and test cases, in a wide range of both Reynolds and Mach numbers, are presented. Keywords: Immersed boundary method, Moving boundary, Split convective operators, Spurious force oscillations, Compressible flows

Email address: [email protected] (Francesco De Vanna) Preprint submitted to Computers and fluids

January 8, 2020

1. Introduction

30

35

40

45

50

55

60

Numerical simulation of compressible viscous flows around moving objects is a crucial target in many engineering applications. These problems range from the control of aeroelastic structures up to the study of particulate flows and surely today the success of landing on dense atmosphere planets [1, 2] or the efficiency improvement of high-bypass turbofan engines [3] are even more affected by the non-linear dynamics between the structural components and the supersonic flow field. In order to simulate the interactions between a moving object inside a fluid two main approaches can be identified: The Body-Fitted Grid Method (BFGM) and the Immersed Boundary Method (IBM). Even if the former is very accurate for the treatment of the boundary layer in complex geometries, dealing with moving bodies the Eulerian mesh needs to be computed and regenerated at each time step, resulting in a very demanding effort. On the other hand, the Immersed Boundary Method allows the body surface to cut the computational cells making it possible to use a fixed Cartesian mesh even with complex moving objects. The IBM has been originally developed for biological incompressible flows by Peskin [4] to deal with elastic moving membranes. After this pioneering work several formulations have been proposed to simulate different flow problems, see e.g. Fadlun et al. [5], Iaccarino and Verzicco [6], Mittal and Iaccarino [7], Eshghinejadfard et al. [8] for application to turbulent incompressible flows, Uhlmann [9], Breugem [10], Kempe and Fr¨ohlich [11], Picano et al. [12], Schwarz et al. [13], Eshghinejadfard and Th´evenin [14], Abdol Azis et al. [15] for interface resolved particulate flow. More recently, the IBM has been applied also to thermal flow problems by Ren et al. [16], Wang et al. [17] and to compressible flows by Merlin et al. [18], Luo et al. [19], Piquet et al. [20], Bernardini et al. [21], Boukharfane et al. [22]. Different families of the IBMs have been proposed in order to comply with the different peculiarities of the applications. Restricting to the class of the discrete forcing methods [6] - where the technique is developed at discrete level and integrated with the flow solver algorithm - the most common approaches are the Direct-Forcing Method (DFM) and the Ghost-Point-Forcing Method (GPFM). The DFM consists in adding a forcing term to the Navier-Stokes system of equations that indirectly imposes the boundary conditions induced by the presence of a solid body. The forcing process can be carried out using 3

65

70

75

80

85

90

95

100

various numerical techniques, both explicit and implicit. However, while the DFM is straightforward when dealing with Dirichlet boundary condition, it often involves complicated implementation strategies to accurately impose Neumann boundary conditions. Hence, the DFM is usually implemented in incompressible flow solvers, in which the velocity components are usually assigned at the boundary location, see e.g.[7, 10]. On the other hand, in the GPFM, to impose the boundary condition on the object surface, the mesh nodes inside the body are employed as ghost points. Hence, the body surface becomes a sharp interface between the outer fluid and the inner ghost domain. Using proper techniques to determine the ghost point values, it is possible to deal both with Dirichlet and Neumann boundary conditions in a similar manner. Hence, the method is more suitable for problems where the natural boundary conditions are specified as the gradient of some fluid variables, as in thermal and compressible flows problems. Despite its greater flexibility, the GPFM needs a more sophisticated data management to deal with general geometries, since the method needs to discern whether a computational node is a fluid, a ghost or a solid one. In order to pursue this target, some automatic solid-detection techniques must be employed [6]. Even if the GPFM is more suitable for compressible boundary problem, employing a sharp interface introduces spurious oscillations, a frequent phenomenon especially in the case of moving IBM. Previous studies [23, 9, 24, 10] proved that the problem is affected by pressure instabilities at the wall boundaries, resulting in a noisy behaviour in terms of lift, drag and pressure coefficients. The problem is even more persistent in case of moving objects in compressible flow conditions. The co-located structure of the numerical methods for Euler equations, in combination with the difficulties associated with the thermodynamical nature of pressure, makes compressible simulations very delicate if combined with moving boundaries. In this respect, many authors adopted an ad-hoc numerical technique to reduce the boundary noise. However most of the works in the field of moving boundaries have been focused on incompressible or slightly-compressible flow conditions, while few methodologies have been formulated for strong compressible simulations [25]. In their pioneering work, Guilmineau and Queutey [23] made use of the Consistent Physical Interpolation (CPI) method previously developed by Deng et al. [26]. In [23], the main problem of moving boundary has arisen and the authors understood the need for increasing the pressure-velocity coupling near the solid interface in order to control the number of numerical oscillations. 4

105

110

115

120

125

130

135

140

On the other hand, Uhlmann [9] decided to relax the fluid-solid interface making use of a smooth Dirac delta interpolation able to slow the transition between the Cartesian mesh and the boundary points. As Breugem [10] noticed, the latter method required the solid to be reduced in terms of its size in order to accurately reproduce its force exchange with the fluid. Yang and Balaras [27] proposed a sharp interface method in combination with an incompressible large-eddy model. Here, to reduce spurious pressure and velocity noise a field-extension procedure was employed. Always in the field of incompressible flows, Chi et al. [28] recently developed a directional ghostcell immersed boundary method in which fictitious discrete boundary forcing terms have been involved instead of the ghost values in the governing equations. The method was proved to deliver a lower truncation error compared to standard IBM incompressible procedures. Schneiders et al. [24] developed a sophisticated technique combining a sharp IBM with a dissipation function regulating the fluid-solid transition. The method was successfully applied to compressible flow simulation in low-Mach conditions. Ehsan Khalili et al. [29] presented a GPFM for compressible flows combined with a high-order low-pass filter. Thus, filtering all the conservative variables, they were able to cut down the numerical oscillations and improve the stability of the whole solver. While these techniques do not impact on the quality of the results in highly viscous and in inviscid conditions, they should be avoided when dealing with turbulent direct numerical and large-eddy simulations. In the present computational method, no artificial stabilisation techniques in terms of upwinding or filtering are implemented. The main aim of the present paper is to show that employing the fully-split formulation developed by KennedyGr¨ uber and Pirozzoli [30] in the discretisation of the convective derivatives is a suitable strategy to drastically reduce the oscillations at the boundary location. The peculiarity of the scheme, in fact, is the stronger coupling between the pressure and the velocity fields; the problem representing one of the most critical issues related to the boundary noisy behaviour. Previous studies ranging on different purposes [31, 32, 33, 34, 35] have already proved the method to be able to stably evolve strong gradients in terms of thermodynamic and mechanical flow properties without introducing any additional numerical dissipation. Here we show that using the split-convective formulation for the smooth part of the flow hybridised with a standard WENO scheme for the correct shock-capturing combined with GPFM is an effective recipe to accurately and stably evolve the complex dynamic of moving 5

145

150

objects in strong compressible flows. Avoiding any filtering or upwinding techniques and introducing a minimal numerical viscosity just around the shock locations, the proposed methodology appears a valid strategy for accurate direct and large-eddy simulations of moving objects in compressible turbulent flows. The present paper is organised as follow: In the first section, we will briefly describe the numerical algorithms and the methodology we employed in the discretisation of the Navier-Stokes equations; in the second part, the GPFM will be presented in terms of its data management and numerics; finally, a collection of the main results, performed over the most documented aerodynamic test cases available in literature, will be shown and discussed. 2. Numerical Methodolody

155

2.1. Governing equations The solver used in the present study, URANOS, an acronym standing for Unsteady Robust All-around Navier-StOkes Solver, has been recently developed at the Department of Industrial Engineering, University of Padova. The code implements the three-dimensional Navier-Stokes equations written for an ideal gas in conservative formulation: ∂Fj (U) ∂Fvj (U) ∂U =− + , ∂t ∂xj ∂xj

j = 1, . . . , 3

(1)

where 

ρ







ρui



0



             , Fj = ρui uj + pδij  , Fvj =   σ ρu U= ij i             ∂T σij uj + λ ∂xi ρE (ρE + p)ui

160

165

(2)

denote the vector of conservative variables, the vector of convective and viscous fluxes in the jth direction, respectively. In particular, ρ is the density, ui is the velocity component in the ith direction, p is the thermodynamic pressure, T is the temperature, E = e + ρui ui /2 is the total energy per unit mass, e is the internal energy per unit mass, R is the gas constant, γ = cp /cv is the specific heat ratio, λ is the thermal conductivity of the fluid. The viscous stress tensor σij is expressed as: 6

 2 ∂us ∂ui ∂uj + − . (3) σij = µ ∂xj ∂xi 3 ∂xs The above system of conservation equations is completed by the ideal gas equation of state 

p = ρRT

(4)

and a constitutive equation for the total energy pR ρ + (ui ui ) (5) (γ − 1) 2 Finally viscosity µ is assumed to obey to the Sutherland’s law   T0 + S 3/2 µ(T ) = T (6) T +S where T0 = 273.15K and S = 110.4K are empirical parameters. The flow solver relies on a co-located finite-difference (FD) conservative approach as described in the following. ρE =

170

175

180

185

2.2. Discretisation schemes The main scheme adopted for the discretization of the convective derivatives is the sixth-order fully-split-convective energy-preserving scheme by Kennedy, Gr¨ uber and Pirozzoli (KGP6) [30]. Dealing with supersonic conditions with shocks such scheme has been hybridised with a fifth-order weightedessentially non-oscillatory scheme (WENO5) [36]. For the sake of clarity we explain the two schemes and the shock detector used to hybridise the schemes separately. The main feature of the KGP6 method, compared to other central-finite difference approximations, consists in granting the semidiscrete preservation of the total kinetic energy in the limit of inviscid, incompressible flows. The property was found to be effective also in case of non-uniform Cartesian grids [31], thus allowing for a robust spatial discretisation method for convective terms without the addition of numerical dissipation in the form of upwinding or filtering. An extensive description of the method can be found in [30] and [35] and here a brief overview will be presented. Being  1 ˆ ∂ρuj φ ˆi−1/2 ' f − f (7) i+1/2 ∂xj i ∆xj 7

190

a convective derivative associated to a generic transported scalar quantity (i.e. unity for the continuity equations, {ui }3i=1 for the momentum equation and H = γ/(γ − 1)p/ρ + ui ui /2 for the total energy equation) and fˆi+1/2 the numerical flux at the cell face i + 1/2, a fully split approximation of equation (7) is expressed by   ∂ρuj φ ∂ρuj φ ∂ρφ ∂uj φ ∂ρuj = k1 + k2 uj +ρ +φ + ∂xj ∂xj ∂xj ∂xj ∂xj   ∂φ ∂uj ∂ρ (1 − k1 − k2 ) ρuj + ρφ + uj φ (8) ∂xj ∂xj ∂xj

195

As shown by Pirozzoli [30] a conservative approximation for (8) is granted only if k1 = k2 = 1/4 and the numerical flux associated with the split formulation can be written as fˆi+1/2 = 2

L X l=1

where

200

205

210

al

l−1 X

^ (ρ, u, φ)i−m,l

(9)

m=0

1 ^ (ρ, u, φ)i,l = (ρi + ρi+l )(ui + ui+l )(φi + φi+l ) (10) 8 is the two-point, three-variables discrete averaging operator associated with the transported variable φ. The coefficient al maximise the formal order of accuracy of a central approximation of a 2L-size stencil. Compared to other split formulations available in literature (see e.g. Morinishi [37]) the one proposed by Pirozzoli [30] - and in particular the equation (10) - takes into account separately the thermo-mechanical effect of density variation from the other variables, and in this sense it can be considered a densityweighted conservative discretisation. In this respect, the method was found to be extremely robust and accurate in a wide range of applications ranging on very different purposes see e.g. [31, 32, 33, 34, 35]. Even if the KGP6 method is considerably more robust than a standard central-finite difference approximation, the method can not be employed to deal with discontinuities as shocks. For these regions, the present work adopts a WENO5 reconstruction. WENO scheme, originally developed by Jiang and Shu [36] improves the fundamental suggestion of the essentially nonoscillatory (ENO) schemes, which were introduced by Harten et al. [38]. 8

215

220

The key idea behind ENO is to employ the smoothest stencil among several candidates in order to discretise the fluxes at cell interface. The process grants both a high-order of accuracy and avoids spurious oscillations near shocks. After its first implementation, a large variety of WENO versions have been developed to improve the scheme resolution (see e.g. Balsara and Shu [39], Balsara et al. [40]) and in the present work, the so-called WENOZ method by Castro et al. [41] was employed. Considering a conservative approximation of a convective flux in the form of  1 ˆ ∂Fj ' fi+1/2 − ˆ fi−1/2 ∂xj ∆xj j

225

230

235

240

(11)

the aims of WENO is to provide a high order reconstruction of the numerical flux ˆ fi+1/2 at the cell-bound. This goal is pursued using the following procedure: The first step computes the Roe’s [42] averaged state at the cell face i + 1/2 and determines the left and right eigenmatrices Li+1/2 , Ri+1/2 associated to the Jacobian of flux component [43]. Thus, in order + to split the fluxes associated with each s-characteristic direction fs,i+1/2 − and fs,i+1/2 [44] the local-Lax-Friedrichs (LLF) flux splitting, according to ± fi+1/2 = 1/2Li+1/2 (Fi ± |λmax |Ui ) is employed, and here λmax denotes the maximum eigenvalue of the Jacobian matrix associated to the flux component Fj . For every s-direction the WENO interpolation procedure is performed. + In the following, we will explain how fs,i+1/2 is computed and, for the sake of clarity, we drop the + sign in the superscript. The formulas for the negative part of the split flux are symmetric in respect of the i+1/2. Thus, interpolat+ ing fi+1/2 with three different polynomia on three different stencils according P (j) to fi+1/2 = 2j=0 wj qi+1/2 , the key point of WENO is the determination of the non-linear weights {wj }2j=0 . The latter are expressed by ! |β0 − β2 |2 αj , where αj = γj 1 + (12) wj = P (βj + ε)2 j αj and ε is a small coefficient, in our computations set to 10−10 , whose role is to avoid a vanishing denominator. The {βj }2j=0 functions denote the smoothness indicators and their formulation, as far as the equations for the (j) three interpolation polynomia qi+1/2 , can be found in Jiang and Shu [36]. Finally, the coefficients {γj }2j=0 are the linear weights of the scheme resulting ± in γj = 1/10 {3, 6, 1}T . Once the two flux components, fi+/2 , have been 9

245

250

255

260

265

270

computed the procedure ends with their recombination in the physical space + − according to ˆ fi+1/2 = Ri+1/2 (fi+1/2 + fi+1/2 ) equation that represents a fifth order approximation of the numerical flux at the left cell bound. The WENO procedure is well-known to be robust and stable but, especially in low-Mach conditions and/or in turbulent flow regions adds some artificial dissipation. Moreover, the WENO reconstruction represents a computational demanding effort and because of these reasons, a hybridisation strategy is desirable. In the present work, the WENO5 scheme was confined only around shock locations, while in the rest of the flow domain the KGP6 method was used. The latter strategy improves the accuracy of the results and minimises the non-linear spatial operations per time-step, optimising the efficiency of the entire algorithm. Here the shock detection technique proposed by Ducros [45] was employed and, at each time-step, the shocks were located through a modified Ducros sensor ! div(u) ,0 (13) θ = max − p (div(u))2 + rot(u)2 + ε

where ε = (u∞ /c)2 and c is the local speed of sound. Here a threshold ¯ of θ = 0.05 was set for equation (13), so that, if θ > θ¯ the flow is considered shocked and the WENO5 reconstruction procedure is activated in the actual cell and in all the three proximal cells, otherwise the KGP6 scheme is employed. To avoid odd-even decoupling, the viscous terms were also discretised using a sixth-order central finite difference approximation after being expanded in Laplacian form. Also, in this case, such an arrangement improves the spectral properties of the scheme, increasing the wavenumber resolution and the scheme stability [46]. The time is advanced using a third-order Total-Variation-Diminishing (TVD) low-storage Runge-Kutta scheme following the one proposed by Gottlieb and Shu [47]. Being N (U) a high-order representation of the the nonlinear spatial differential operator applied to the conserved variables U the

10

scheme reads: U(1) =Un + ∆tN (Un ) 3 1 1 U(2) = Un + U(1) + ∆tN (U(1) ) 4 4 4 1 n 2 (2) 2 (n+1) U = U + U + ∆tN (U(2) ) 3 3 3 275

280

285

290

295

300

(14a) (14b) (14c)

The method was proved to be the optimal third-order TVD Runge-Kutta scheme, i.e. the third-order TVD method with the higher value of the CourantFriedrich-Lewy parameter. In particular it was proved by Gottlieb and Shu [47] that the stability of the present Runge-Kutta scheme is recovered up to CF Lmax = 1. In our numerical tests we employed a CF L number equals to 0.5 everywhere (apart where differently stated). The value represents the best compromise in terms of speed and stability of the scheme and was largely applied by other authors in literature [24, 22]. 2.3. The Immersed Boundary method In this section, we will explain the details of the adopted IBM through three main parts: The first one explains the data structure employed in the determination of a grid point status (i.e. fluid, ghost or solid); the second one shows the numerics behind the GPFM and, finally, the third part explains how the boundary conditions are enforced. 2.3.1. Immersed boundary data structure The implementation of the data structure is a critical aspect of dealing with GPFM. In order to discretise a two-dimensional object the .msh format has been employed so that the immersed solid was represented by a series of n−1 on a Cartesian mesh. In the present closed points {xl }nl=1 and edges {el }l=1 approach, we decided to link the number of points n employed in the solid discretisation to the minimum grid step using a criterion based on mutual distances. The procedure makes it possible to minimise the geometrical discretisation memory storage, saving just the points whose mutual distance is comparable to the minimum grid step. Once the strictly necessary points have been stored, the geometry is represented using a many-edges-polygon, a suitable model for the most common algorithms in Computational Geometry. Following O’Rourke [48] every Cartesian grid point is flagged in respect of the geometry position, determining whether it is inside or outside the polygon, 11

305

or if it is a ghost point. Thus, being D the whole computational domain, this is divided into three portions: Ωsolid , Ωf luid and Ωghost and the information is stored inside a marker variable ζ defined as   0 if x = (x, y) ∈ Ωsolid ζ(x) = 1 if x = (x, y) ∈ Ωghost (15)   2 if x = (x, y) ∈ Ωf luid

thus the set of Cartesian points belonging to the solid and fluid regions respectively are defined by Ωsolid = {x ∈ D | ζ(x) = 0}

Ωf luid = {x ∈ D | ζ(x) = 2}

(16) (17)

It is important to highlight that Ωsolid ∪ Ωf luid = D so that the ghost region is a function of both the solid and fluid portions defined as Ωghost = {(xi , yj ) ∈ Ωsolid if ∃(xm , yl ) ∈ Ωf luid for m = i−3 . . . , i+3, l = j−3, . . . , j+3} (18) In particular, the ghost region includes three layers of points which are needed for the high-order numerical scheme stencils. A distinction must be made between steady and moving boundary simulations. In the case of steady objects, the geometrical computation of the mapping function ζ = ζ(x) and the grid point status determination is performed once, at the beginning of the simulation. On the other hand, dealing with moving objects the polygon tracing and the points status must be recomputed every time step. In addition, in case of moving solids, a fourth value of the marker function ζ must be taken into account considering the so called fluid emerging points. The fluid emerging points at time level t = tn represent le grid locations belonging to the following set  Ωf e (tn ) = x ∈ R3 | x(tn−1 ) ∈ Ωghost ∧ x(tn ) ∈ Ωf luid (19) Because at the time level tn these points do not have a Runge-Kutta history, the time integration process could fail in terms of accuracy, resulting in a drastical drop of the local numerical error. To reduce this source of instabilities and spurious oscillations, the present method treats the fluid emerging points with an ad-hoc interpolation procedure that is explained in 12

the following. In terms of data management a fourth value of the marker variable ζ = ζ(x) must be taken into account. Employing the definition of Ωf e the marker variable associated with the grid flagging procedure can be compared time-step per time-step. The comparison highlights the differences between the status (i.e. solid, ghost or fluid) of the grid points and univocally determines the fluid emerging set as Ωf e (tn ) = {x ∈ R3 | ζ(x, tn−1 ) = 1 ∧ ζ(x, tn ) = 2}

(20)

The values of the primitives variables (i.e. u, v, w, p, T ) are interpolated employing a discrete Dirac delta function, thus being φ a scalar field the latter reads Z ∞ φ(x)δ(x − x0 )dx. (21) φ(x0 ) = −∞

where

310

  3 Y 1 ¯ xi − xi,0 φ δ(x − x0 ) = ∆x ∆xi i i=1

(22)

¯ Here φ(r) denotes the kernel of the Dirac delta distribution for which the smoothest formulation of Yang et al. [49] was employed. Figure 1 reports a sketch of the interpolation strategy.

13

R

Figure 1: Sketch of the interpolation procedure around a fluid emerging point. Here the dashed black line represents the position of the boundary at time t(n−1) while the black solid line is the fluid-solid interface at time t(n) . The orange points are the fluid emerging points and the blue circle represents the support employed for the interpolation.

14

2.3.2. The Ghost Point Forcing Method

Figure 2: Sketch of the boundary interpolation procedure. Here the blue points represent the ghost points, the light blue points represent the fluid nodes while x1 , x2 , x3 , x4 represent the points retained for the interpolation.

315

320

325

Figure 2a shows a graphical representation of the different interpolation procedures adopted. Once a mapping function ζ = ζ(x) is known, for each point flagged as ghosts a normal probe is built in order to locate the correspondent image point xip in the fluid domain. In the case of analytical boundary curves ψ (or surfaces) the image point associated to a ghost point is unique and it is defined as xip = {x1 ∈ Ωf luid such that ∃! x2 ∈ Ωghost | kx1 − xb k = kxb − x2 k with (x1 − xb ) and (xb − x2 )⊥ψ}. In a discrete path instead an approximation of the boundary curve is needed. In particular since the solid is discretised by a many-edges-polygon, so being xg a ghost point, we can easily find the nearest edge ne in respect of xg . Once ne is found we compute the line passing though xg and orthogonal to ne . This line cross the edge in a point xb that, in our computation, is considered a good approximation of the boundary point associated with xg . Knowing the distance d = kxb − xg k between the bound and the ghost point, the image point xip is located as xip = xg + 2nd

(23)

where n represents the local normal direction through the bound associated with xg , defined as n= 330

xb − xg kxb − xg k

(24)

Once the image point xip is found, the flow field variables are bilinearly interpolated from surrounding values. A generic flow variable φ is assumed to obey to 15

φ(x, y) = c1 xy + c2 x + c3 y + c4

335

340

345

350

(25)

with the coefficients {ci }4i=1 to be determined. Being {φi }4i=1 the values of φ at the surrounding coordinates {(xi , yi )}4i=1 , the interpolation coefficients are determined by solving the linear system,      x 1 y 1 x 1 y 1 1 c1 φ1 x2 y2 x2 y2 1 c2  φ2       (26) x3 y3 x3 y3 1 c3  = φ3  . x 4 y 4 x 4 y 4 1 c4 φ4

Some exceptions must be pointed out concerning the main procedure. For the case where one (or two) points retained for the interpolation are flagged as ghosts - so that when the interpolation cell is cut by the boundary curve - the number of fluid points for the variable reconstruction reduce to three or two. In this case, the numerical accuracy drastically drops and the interpolation procedure becomes a source of numerical instabilities and spurious oscillations. Following the sketch reported in Figure 2b, Figure 2c, in this situation we directly impose the wall value at the correspondent matrix line using the coordinate of the boundary interception point. A distinction between Dirichlet and Neumann boundary condition is needed. For instance, if (x2 , y2 ) ∈ Ωghosts and we are dealing with Dirichlet boundary the correspondent system becomes:      x 1 y 1 x 1 y 1 1 c1 φ1  xb yb xb yb 1 c2  φb       (27) x3 y3 x3 y3 1 c3  = φ3  x 4 y 4 x 4 y 4 1 c4 φ4 otherwise, when dealing with Neumann boundary condition the system becomes:      φ1
x1 y 1 x 1 y 1 1 c1 yb nx + xb ny nx ny 0 c2   ∂φ      =  ∂n b  (28)  x3 y 3 x3 y3 1 c3   φ3  x4 y 4 x 4 y 4 1 c4 φ4 where the line corresponding to (x2 , y2 ) is replaced by ∇φb · nb .

16

355

2.3.3. Boundary conditions Once the flow variable φip is interpolated at the image point position the wall condition is enforced reflecting its value in the correspondent ghost point. In particular for Dirichlet boundary condition, being φb is the desired value of φ at the bound location, the value of φg in the ghost position is expressed by φg = 2φb − φip .

(29)

For inhomogeneous Neumann boundary conditions instead the value of φg is given by   ∂φ , (30) φg = φip − ∆l ∂n b 360

∂φ is the desired value of the φ gradient at the Here ∆l = |xip − xg | while ∂n solid-fluid interface. In our computations, we enforced the boundary conditions on primitive variables {ui }3i=1 , T and p following the method proposed by Piquet et al. [20] so that the equations (29) and (30) become

ui,g = 2ui,b − ui,ip , Tg = Tip − ∆l pg = pip − ∆l

365

i = 1, . . . , 3 



∂T ∂n ∂p ∂n





(31) (32)

b

(33)

b

In particular, for adiabatic wall (∂T /∂n)b = 0 for both steady and moving solids, while the pressure gradient (∂p/∂n)b is expressed by (   0 if ub = 0 ∂p
= (34) Dub ∂n b − ρ Dt · nb otherwise 2.4. Parallelisation strategy The entire numerical methodology can take the advantages of recent massive-parallel numerical architectures. In particular, the core of the numerical solver makes use of a fully three-dimensional domain decomposition 17

370

375

380

385

390

395

400

405

sharing the halo-nodes thought non-blocking MPI communications. The latter process has been efficiently addressed thought point-to-point directives and has been partially desynchronised concerning the computational kernels. As far as the coupling with the Immersed Boundary Method, the initialisation process reads the geometry (i.e. from an external .msh file at the beginning of the simulation) and stores it inside a global boundary array. This set of data is known by all the computational units. Consequently, the parallelisation strategy employs the following steps: 1. The global boundary array is distributed thought the processors looking at the geometry’s point location in respect to the processors’ bounds; thus for every computational unit a local boundary array is saved in a memory buffer. 2. To discern the solid and the fluid domains, the ray-tracing algorithm is employed. To reduce the complexity of a fully parallel ray-tracing procedure, the present method looks at the global boundary array but flags just the Cartesian points belonging to a local computational unit. This implementation grants parallel scalability concerning the number of grid points but does not provide any speed-up associated with the geometry splitting. Anyway, because the number of points which describe the geometry is always orders of magnitude less than the number of Cartesian locations, the strategy was found very efficient and avoids any compliances related to MPI communications inside the ray-tracing computational kernel. 3. Once the solid location has been determined for every computational unit, the ghost nodes can be easily found employing the definitions of the ghost-region. The integer marker variable ζ = ζ(x) is then locally stored to distinguish the different computational areas and the halopoints information is shared through standard MPI procedures. 4. Finally, once the grid points have been flagged (i.e. as solid, ghost or fluid) the ghost points are employed to detail the boundary interface. The normal direction n and the boundary point location xb associated with each ghost node xg are then computed and stored inside local memory buffers. The latter procedure is performed once, at the beginning of the simulation, if the boundary is stationary. Thus, in the case of steady objects, the geometry computation does not impact the solver performances. On the other hand, if the solid is moving, this algorithm is repeated every Runge-Kutta sub-step. 18

Moreover, in the case of moving objects, the tracing of the fluid emerging points introduces an additional overload to the computation. 3. Results and discussions 410

415

420

425

430

435

440

The present section aims at validating and discussing the results obtained using our methodology. To highlight the cleanliness offered by the fullysplit convective scheme in shockless flows regions, several tests have been performed in low-Mach regimes. In these conditions, the convective terms are discretised using KGP6 method coupled with the GPFM, do to avoid any additional uncertainty associated with the hybridisation scheme with WENO (not needed for shockless flows). The results description will continue addressing the supersonic flow regime and these configurations require the WENO reconstruction procedure to be activated around the shocks’ locations. 3.1. Subsonic flow past a confined cylinder at various Re numbers and D/Ly ratios Standard Navier-Stokes Characteristic Boundary Conditions (NSCBC) (Poinsot [50], Lodato et al. [51], Pirozzoli and Colonius [52]) were imposed at the inlet and outlet boundary while no-slip wall conditions were enforced at the top and at the bottom boundary respectively. The computation have been performed at various Reynolds numbers Re = u∞ D/ν = {1, 10, 20, 40, 100}T and at D/Ly = {0.5, 0.7}T . A parabolic inflow velocity profile was imposed at the left boundary and here u∞ denotes the pick value of inlet speed, D is the cylinder diameter and ν the kinematic viscosity of the fluid. The initial values for the non-dimensional pressure and temperature were set to p/p0 = T /T0 = 1 respectively, while initial velocity components are set to zero. In order to avoid any kind of compressibility effect, the u∞ value was set in such a way that the free stream Mach number M∞ equals to 0.01 at the inlet position. Here a note is demanding as far as the case at Re = 1. This low Reynolds number could be a source of uncertainties concerning the validity of the Navier-Stokes system of equations. Remaining inside the Continuum Mechanics Hypothesis in fact, the latter formalism requires the Knudsen number (Kn) - i.e. the ratio between the molecular mean free path (λ) and a characteristic length of the system (L0 ) - to be much smaller of the unit. Practical applications

19

445

450

455

set to 0.01 the upper limit of Kn. On the other hand, theoretical considerations link Kn to pthe Mach and the Reynolds numbers through the expression Kn = M a/Re γπ/2. In this respect, the case at Re = 1 shows Kn = 0.0148 which can be considered in the limit of the Navier-Stokes model applicability. The computational grid consisted of a structured uniform grid featuring Nx × Ny = 1600 × 400 = 640000 computational nodes. The time-step was evaluated from the Courant-Friedrichs-Levy stability condition with a CFL number equals to 0.5. Interpolating the computed pressure and velocity fields bilinearly around the boundary positions the drag coefficient was calculated. In particular, the following formulation was implemented: H H pb nb · ˆidSb + µσb nb · ˆidSb (35) cD = 1 2 γp∞ M∞ 2 here pb , nb , σb represents respectively the pressure, the normal outer direction and the viscous stress tensor at the boundary position while ˆi is the unit vector associated to the u∞ flow direction. In Figure 3, the computed values of the drag coefficients are compared with the results obtained by Sahin and Owens [53] showing a very good agreement for all the Reynolds numbers investigated.

20

Figure 3: Comparison of the computed drag coeffient for a subsonic flow past a confined at M = 0.01 with various Reynolds number and D/Ly ratios. In figure 3a the drag coeffient obtained with the present method (black points) is reported as a function of the Reynolds number for the cases at D/Ly = {0.5, 0.7}T . The results are compared with the simulations of Sahin and Owens [53](solid lines). Figure 3b shows the stationary condition for velocity field at Re = 100 and D/Ly = 0.5. The simulation was performed on a uniform Cartesian grid consisted in Nx × Ny = 1600 × 400 grid points for a domain of [8D; 2D] and two diameters were reserved upstream the cylinder centre location.

460

465

470

3.2. Subsonic flow past a free cylinder at various Re numbers A similar domain configuration as described the previous paragraph was employed to simulate the subsonic flow past a cylinder in open domain. Again NSCBC were imposed at the inlet and outlet portion of the domain while symmetry conditions were enforced at the top and the bottom boundary. It should be noticed that, considering that the transverse dimension of the domain Ly /D = 20 the choice of the lateral boundary does not produce any significant confinement phenomenon. The streamwise dimension was set equal to Lx /D = 50. Here we simulated the flow at Re = u∞ D/ν = {10, 20, 40, 70, 100, 300}T . A uniform velocity profile is imposed at the inlet boundary and here u∞ denotes the value of the inflow speed, D is the circle diameter and ν the kinematic viscosity of the fluid. Following the previous tests, ambient conditions were initialised at t/t0 = 0 while the u∞ value was set in order to ensure M = 0.05 at the inlet flow position. The computational grid consisted of a structured non-uniform grid, stretched around the position 21

475

of the cylinder and featuring Nx ×Ny = 800×320 = 256000 grid points. Grid stretching was applied along both x and y directions but, near the cylinder, a uniform grid spacing, approximately with D/∆x = 50, was employed. The CFL number was set equal to 0.5. In Figure 4a a global view of the computational grid is shown.

Figure 4: Cylinder subsonic laminar aerodynamics in open domain configuration. Figure 4a shows a global view of the computational grid employed for the simulation. The latter consisted in a non-uniform Cartesian mesh featuring Nx × Ny = 800 × 320 computational points for a domain of [50D; 20D]; 13D were reserved upstream the cylinder centre location and a hyperbolic tangent distribution was employed to cluster the node around the body surface. In panel 4b the Mach field and the wake recirculation for the case at Re = 40 and M = 0.05 is reported.

480

485

A large number of results, both numerical and experimental, have been produced regarding the free cylinder aerodynamics in subsonic flow regime (see e.g. [54, 55, 56, 57, 58, 22]) and here we present our results compared with the most accurate data available in the literature. We divided the test campaign into two main groups. Because of the intrinsic unsteadiness of the cases at Re = 100 and Re = 300 we will treat them separately. Focusing on the configurations at Re = 10, 20, 40, 70, in table Table 1 the values of the drag coefficient are compared with the one in literature; in Figure 4b the velocity field and the wake recirculation are reported while 22

in Figure 5 the circular distribution of the pressure coefficient, i.e. cp (θ) = 2(p(θ) − p∞ )/(γp∞ M 2 ) is shown compared to the one obtained by Dennis and Chang [54].

Dennis and Chang [54] Tritton [55] De Palma et al. [56] Linnick and Fasel [57] Schneiders et al. [58] Boukharfane et al. [22] Present

cD (Re = 10) 2.85 3.07 2.83 2.99 2.93

cD (Re = 20) 2.05 2.09 2.05 2.06 2.15 2.07 2.09

cD (Re = 40) 1.52 1.59 1.55 1.54 1.55 1.52 1.55

cD (Re = 70) 1.21 1.30 1.25

Table 1: Drag coefficient literature comparison for a laminar flow past a free cylinder at various Reynolds numbers.

Figure 5: Comparison of the pressure coeffiecient circular distribution around a steady cylinder at M = 0.05 and various Reynolds numbers. The figures show the circular distribution of the pressure coefficient for a cylinder in open domain at Re = {10, 20, 40, 70}T and M = 0.05. Present results (black points) are compared with the one of Dennis and Chang [54] (blue dashed lines).

490

As far as the tests at Re = 100 and Re = 300 are concerned, the intrinsic unsteadiness of the flow needs to be treated statistically. In particular, we focused our attention on three parameters: the time-averaged drag coefficient, the root mean square of the lift coefficient and the Strouhal number of the wake defined as St = f D/u∞ . Here f denotes the wake frequency, D 23

495

500

the cylinder diameter and u∞ the free stream velocity. In Figure 6 the nondimensional vorticity field and the time history of the lift and drag coefficient are reported. As can be appreciated from the snapshots the vorticity fields do not show oscillations or noise. The time histories of CL and CD are also very clean. Values obtained by present computation always lie in the range reported by other studies. In Table 2 the results of the present computation are compared with the one available in the literature.

Figure 6: Von Karman wake behind a steady cylinder in open domain configuration at Re = {100, 300}T and M = 0.05. Panels 6a and 6b reports the results in term of instantaneous non-dimensional vorticity contours (ωz ) while figures 6c and 6c show the lift and drag history in a statistically-steady range of time. The simulation has been performed over the same grid in §3.2.

Re Ghias et al. [59] De Palma et al. [56] Rajani et al. [60] Boukharfane et al. [22] Present

cLrms 100 0.32 0.23 0.17 0.25 0.22

300 0.67 0.60 0.62 0.61

c¯D 100 1.36 1.32 1.33 1.36 1.32

300 1.40 1.28 1.26 1.34

St 100 0.16 0.16 0.15 0.16 0.16

300 0.21 0.21 0.21 0.21

Table 2: Literature comparison of force coefficient statistics for a steady cylinder at Re = 100 and Re = 300 and M = 0.05.

24

505

510

515

520

3.3. Subsonic flow past a confined square at various Re numbers A similar domain configuration as in §3.1 was used for the flow prediction of a confined square at various Reynolds numbers. NSCBC boundary conditions were imposed at the inlet and outlet portions of the domain while no-slip wall conditions were set at the top and bottom sides. The computation has been run at Re = u∞ D/ν = {10, 20, 30, 40, 50}T . Here u∞ denotes the maximum of the parabolic inflow profile, D is the square diameter and ν the kinematic viscosity of the fluid. The flow was initialised with p/p0 = T /T0 = 1 respectively, while the initial velocity components were set to zero. To avoid any uncertainty due to the flow compressibility, the inflow Mach number M∞ has been fixed at 0.01. A grid sensitivity campaign, consisting in two non-uniform Cartesian meshes was performed. In particular, the two computational grids featured Nx × Ny = [200 × 50]; [1000 × 200] grid points for a domain of Lx × Ly = [40D; 8D]. With this settings, the grid stretching granted D/h = {25; 100}T number of points per diameter for a blockage ratio of B = D/Ly = 1/8 as in Breuer et al. [61]. The time step was evaluated setting the CF L number equal to 0.5. In figures 7a and 7b the drag coefficient and the non-dimensional recirculation length are reported as a function of both the Reynolds number and the grid resolution. A very good agreement is recovered compared to the finite-volume solution of Breuer et al. [61] even for the coarser grid.

25

Figure 7: Comparison of the computed drag coefficient (7a) and the recirculation length (7b) for a subsonic flow past a confined square at M = 0.01 with various Reynolds numbers and with a blockage ratio B = 1/8. In panel 7c an enlargement of the recirculation bubble for the case at Re = 50 is reported. The results are compared with the simulations of Breuer et al. [61] (solid blue lines). The simulation was performed on two non-uniform Cartesian meshes stretched around the square with increasingly resolution Nx × Ny = [200 × 50]; [1000 × 200] for a domain of Lx × Ly = [40D; 8D].

525

530

535

3.4. Subsonic unsteady flow past a square To boost the reader in the validity of the present methodology even in case of non-smooth geometries here we consider the unsteady flow past a square in open domain. Due to the effect of the sharp-corners in respect to the wake dynamics, the test represents a very difficult benchmark for any nonconformal Immersed Boundary Method. Again NSCBC boundary conditions were imposed at the inlet and outlet locations while symmetry conditions were enforced at the top and bottom sides. The simulation was performed at Re = u∞ D/ν = 100. Here u∞ denotes the value of the uniform free-stream speed, while D and ν follow the usual conventions. As in previous tests the flow was initialised with p/p0 = T /T0 = 1 at t/t0 = 0 while the velocity components were set equal to the free stream value u∞ in all the computational domain. The latter was expressed as a function of the free-stream Mach number ensuring M∞ = 0.05. To avoid any confinement phenomena, the domain fit Lx ×Ly = [50D ×20D] and 15 diameters were reserved upstream the 26

540

545

550

square centre. A grid sensitivity campaign was performed consisting of two structured non-uniform Cartesian meshes stretched around the square location. The meshes featured respectively Nx × Ny = [200 × 80]; [400 × 160] grid points and granted D/h = {25; 50}T computational nodes over the square diameter. The time was advanced setting a CF L number equal to 0.5. As in all the subsonic tests cases, to avoid any help related to WENO5 scheme’s inherent diffusivity, just the KGP6 scheme has been employed for the discretisation of the convective terms. In figure 8a the non-dimensional vorticity contours (ωz ) are reported, while panels 8b and 8c show the time history of the drag and lift coefficient in a statistically-steady range of time. As we can notice, both the contours and the signals of the force coefficients are very clean and noise-free, proving the effectiveness of the fully split formulation in combination with the IBM procedure. Finally, table 3 resumes the force coefficients statistics and the Strouhal number of the wake as a function of the grid resolution and compared to well-documented literature benchmarks, concluding a very high level of accuracy of the present methodology even for such difficult benchmark.

Figure 8: Von Karman wake behind a steady cylinder in open domain configuration at Re = 100 and M∞ = 0.05. Panel 8a reports the results in term of the instantaneous nondimensional vorticity contours (ωz ) while figures 8b and 8c shows respectively the time history of the drag and the lift coefficients in a statically-steady range of time. The figures refer to the finer grid level (D/h = 50).

27

Re = 100 Sohankar et al. [62] Darekar and Sherwin [63] Singh et al. [64] Sahu et al. [65] Sen et al. [66] Present (D/h = 25) Present (D/h = 50)

c¯D 1.4778 1.4855 1.5040 1.4895 1.5287 1.5059 1.4861

cLrms 0.1530 0.1864 0.1620 0.1880 0.1928 0.2029 0.1571

St 0.1460 0.1463 0.1480 0.1485 0.1452 0.1427 0.1482

Table 3: Literature comparison of force coefficient statistics for a steady square cylinder at Re = 100 and M∞ = 0.05. 555

560

565

570

575

3.5. Subsonic flow past a moving cylinder with prescribed motion To validate the IBM for moving objects we simulated a transversely oscillating cylinder in a free-stream configuration and low-Mach conditions. Such case is quite well documented both numerically and experimentally (see e.g.[23, 49, 24, 29]). Here the same domain configuration described in the previous paragraph was employed. NSCBC were used at the inlet and outlet portion of the domain while symmetry conditions were employed at the top and the bottom boundaries. To avoid any confinement effect due to the motion of the cylinder the transverse dimension of the domain was doubled compared to the steady case so that Ly /D = 40; the streamwise dimension remains the same. We simulated the flow at Re = u∞ D/ν = 185 where u∞ , D and ν follow the same convections explained in §3.2. The value of u∞ was selected in such a way that a Mach number M∞ = 0.25 is ensured in the free stream conditions following the computation performed by Ehsan Khalili et al. [29]. A grid sensitivity campaign, consisting in three non-uniform Cartesian meshes stretched around the cylinder, was performed. The three meshes feature Nx × Ny = [400 × 320]; [800 × 640]; [1600 × 1280] grid points are denoted as coarse, medium and fine, respectively. Grid stretching was applied in both streamwise and spanwise directions but, near the cylinder, a uniform grid spacing was employed. In particular, the three meshes were characterised by a minimum grid step equals to ∆x = ∆y = [D/25]; [D/50]; [D/100], respectively. The initial values for the non-dimensional pressure and temperature were set to p/p0 = T /T0 = 1 respectively, while the initial value of the velocity components was set equal to the free stream speed in all the domain.

28

580

The CF L number was set equal to 0.5. The body motion was prescribed using a single harmonic oscillation in the form of yc (t) = A sin(2πfe t)

585

590

595

600

(36)

and here t is the dimensionless simulation time, yc (t) is the location of cylinder centre, A is the amplitude associated to the harmonic oscillation and fe is its frequency. Following the available literature the amplitude was set equals to A/D = 0.2 while the frequency fe was expressed as a function of the cylinder natural shedding frequency f0 = St0 D/u∞ . Here St0 denotes the Strouhal number of a steady cylinder under the same Mach and Reynolds regime and experimental and numerical results suggest that St0 is in the range of [1.85 ÷ 1.95]. In our computation we set St0 = 0.195. Following the experiments of Gu et al. [67] and the simulations of Guilmineau and Queutey [23] and Schneiders et al. [24] a spectrum of six frequencies was considered, i.e. fe /f0 = {0.80, 0.90, 1.00, 1.10, 1.12, 1.20}T . In Table 4 the force coefficient statistics for a moving cylinder with a reduced frequency of fe /f0 = 0.8 are reported in comparison with the available literature and as function of the grid spacing, while in Figure 9 the results over the complete frequency spectrum are shown. From the obtained results we can see that the average drag coefficient (¯ cD ) is closer to the results of Yang et al. [49], the root mean square of the lift coefficient (cLrms ) is more aligned with Guilmineau and Queutey [23] and the root mean square of the drag coefficient (cDrms ) features well both the two references. Hence we believe that, within the error bar, our method is in a good agreement with previous studies on moving objects.

Guilmineau and Queutey [23] Yang et al. [49] Schneiders et al. [24] Ehsan Khalili et al. [29] Present (coarse) Present (medium) Present (fine)

M 0 0 0.10 0.25 0.25 0.25 0.25

c¯D cLrms 1.195 0.080 1.281 0.076 1.279 0.082 1.287 0.079 1.248 0.063 1.263 0.058 1.272 0.065

cDrms 0.036 0.042 0.042 0.045 0.044 0.043 0.043

Table 4: Literature comparison of force coefficients statistics for a moving cylinder in low Mach regime, Re = 185 and with a reduced frequency fe /f0 = 0.8.

29

Figure 9: Variation of the force coefficients statistics for a transversely oscillating cylinder at Re = 185 and M = 0.25 respect of fe /f0 . In the figure our computation (black circles) is compared with Guilmineau and Queutey [23] (blue dashed triangles) and Yang and Balaras [27] (red dashed squares).

605

To highlight the low level of noise introduced by the present method, in Figure 10 the drag coefficient for the case at fe /f0 = 0.8 has been plotted as a function of the location of the cylinder. Magnifying the drag history over a single period we can see that even without introducing any additional dissipation the signal is very clean. In Figure 11 the time history of the lift and drag coefficient is reported for all the reduced frequency spectrum, while in Figure 12 the non-dimensional z-vorticity contours are shown.

30

Figure 10: Resolution comparison at different grid spacing. Here the drag coefficient for the case at fe /f0 = 0.8 is plot as a function of the cylinder centre location over a single period.

610

615

620

625

According to the results of our simulations and in particular looking at Figure 9 some comments are demanding to discuss the complex dynamics related to a moving cylinder in a transitional Reynolds regime. Looking at the case at fe /f0 < 1 we can easily conclude that force coefficients statistics depend linearly to this ratio. On the other hand, increasing the wavenumber the trend is not confirmed and something happens in the flow domains. Looking closer to the Figure 12 we can observe that, up to fe /f0 > 1 the interaction between the body motion and the vortex shedding becomes strongly non-linear, resulting in a vortex an opposite circulation detached from the lower side of the cylinder. This fact causes the force exchange to suddenly invert its behaviour, showing a typical knot effect in the system dynamics. The phenomenon is clearly highlighted also in Figure 11. Here the force coefficient time history, up to fe /f0 = 1, exhibits a secondary frequency, higher than the fundamental one, resulting in a secondary vortex shedding. Above the good agreements with the results available in the literature, from this analysis, we can conclude that our methodology is properly predictive in a wide range of body motions and also the strong non-linear effects related to fluid-structure interactions are properly caught and resolved.

31

Figure 11: Time history of the lift and drag coefficient past a moving cylinder in a cross flow configuration with Re = 185, M = 0.25 and a 0.8 ≤ fe /f0 ≤ 1.2. Here the results are associated to the mesh with ∆x = ∆y = D/50.

32

Figure 12: Non-dimensional vorticity contours (ωz ) of a cross flow oscillating cylinder at M = 0.25 and Re = 185. Here the results are presented in function of the oscillation frequency 0.8 ≤ fe /f0 ≤ 1.2.

630

635

640

645

3.6. Shock wave diffraction against a steady cylinder Concerning the validation of solid objects in strong compressible flow conditions, here we present the diffraction of a planar shock wave over a cylinder and over a wedge. The problem related to the shock wave diffraction over steady obstacles represents a widely studied topic in the literature and during the last decades it has been addressed by many authors both experimentally and numerically (see e.g. Bryson and Gross [68], Chang and Chang [69], Chaudhuri et al. [70], Piquet et al. [20] and Boukharfane et al. [22]). As far as the first benchmark, here we will present two versions of the problem: the classical one, where the cylinder is steady and a shock is travelling trough it and a moving version that sees a steady shock and the cylinder impinging to it. Regarding the shock wave diffraction over a wedge just the steady version will be shown and the results are reported in §3.8. The flow field was initialised in such a way that the Rankine-Hugoniot conditions ensure a travelling shock wave moving from left to right with a Mach number equals to 2.81. In particular, the shock, initially located two diameters upstream the cylinder, divided the flow field into two regions and the left part consisted in a gas initially at rest with p/p0 = T /T0 = 1. The Reynolds number of the flow Re = us D/νa was set equal to 2000. Here us denotes the shock speed, D the cylinder diameter and νa the ambient kinematic viscosity of the 33

650

655

660

665

fluid. Extrapolation boundary conditions were ensured in all the four domains edges. The domain size was set equal to Lx × Ly = [10D; 8D] and the cylinder was located at zero coordinates so that four diameters were reserved upstream the cylinder. The computational grid consisted in a structured uniform mesh featuring Nx × Ny = 2000 × 1600 = 3.2 · 106 nodes. A number of points equal to D/∆x = 200 was employed around the cylinder. The CFL number was set equal to 0.5. The complexity of the wave interactions is depicted in Figure 13a. Here the numerical Schlieren of the density field is reported at a reduced time equals to t/t0 = 2. The flow field shows a strong symmetry concerning the centre line of the domain and it evidences some characteristics features such as Mach stems, triple points, shocks, slip lines and vortices. In particular, the trajectory of the two main triple points was extensively reported by previous studies and here, in Figure 13b the results of our computation are compared with the one reported in [68], [70], [20] and [22]. As we notice, in the limit of the uncertainty related to the Bryson and Gross [68]experiments and looking at the numerical results provided by Chaudhuri et al. [70], Piquet et al. [20] and Boukharfane et al. [22] our simulations show a very good agreement compared to all the available data, concluding that our methodology provides very good predictions also in strong compressible flow conditions.

34

Figure 13: Shock-wave-diffraction over a steady cylinder at Ms = 2.81 and Res = 2000. In figure 13a the shock positions has been highlighted thought their density variation, thus the Schlieren density field is reported. Panel 13b shows the WENO-activated-cells while in figure 13c a literature comparison of the triple points trajectory for the problem is reported. Present results are compared with the experimental data of Bryson and Gross [68] and the numerical simulations of Chaudhuri et al. [70], Piquet et al. [20] and Boukharfane et al. [22]. The simulation was performed over a uniform Cartesian grid featuring Nx × Ny = 2000 × 1600 nodes for a domain of [10D; 8D]. Four diameters were reserved upstream the cylinder centre location

3.7. Order of accuracy analysis of the overall methodology To quantify the scaling of the numerical error of the overall methodology the diffraction of a planar shock wave against a steady cylinder was found a significative and challenging benchmark. The problem resumes all the ingredients of the present numerical method (i.e. scheme hybridisation, shock-wave automatic detection, viscous effects and time-dependency). The simulation has been repeated over a set of five uniform Cartesian meshes {Mi }4i=0 , featuring Nx × Ny ' [180 × 180] ∗ {1.5i }4i=0 grid points for a domain of Lx × Ly = [8D × 8D]. The shock Mach and Reynolds numbers were imposed equal to Ms = 2.81 and Res = 1000 while the CF L number - to avoid 35

any uncertainty due to the error contribution of the time-integration - was set equal to 0.1. To assess the error scaling the surface pressure coefficient cp (θ) has been monitored as a function of the grid resolution and the one computed with the i-th grid-level was compared with the M4 results. The complexity of the flow - which features shocks and acoustics wave behind the cylinder - makes the error evaluation quite difficult. In particular, is well known that the WENO5 scheme reduces to zero-order of accuracy in the region where the flow is not smooth. For this reason, the ±80◦ shock-free frontal circular sector was selected as the control zone in which to monitor the error behaviour [70, 22]. The comparative process has been carried out over n = 64 control points (i.e. every 2.5◦ ) inside the control region and the analysis has been performed looking at the results at t/t0 = 1.5. The planar shock wave has been initialised two diameters upstream of the cylinder centre. The Lp and the L∞ norms defined as i LM p =

n p 1 X Mi Cp (θj ) − CpM4 (θj ) n j=1

!1/p

n Mi M4 i LM ∞ = max Cp (θj ) − Cp (θj ) j=1

(37) (38)

have been retained as metrics of the error and in particular, the p = 1 and p = 2 cases have been used in present analysis. Whatever the type of norm q i retained, the numerical error should asymptotically scales as LM (·) = Cε /Ni , where Cε is a constant parameter and q denotes the error order of accuracy. Mi−1 i we obtain Thus, comparing the ratio of the error decay LM (·) /L(·) Mi L(·) M

L(·)i−1

=



Ni−1 Ni

q

 q 1 = r

The latter expression implicitly defines the order of accuracy as   Mi−1 i log LM /L (·) (·) q= log(1/r)

670

(39)

(40)

with r = 1.5 representing the grid ratio of the present case. In figure 14a the pressure coefficient in the control zone is reported as a function of the grid resolution while panel 14b shows the evolution of the metrics associated 36

675

with the numerical error. Finally, table 14 reports the computed errors and the related order od accuracy. For the obtained results we conclude that the present immersed boundary method - in combination with the fluid dynamics solver - keeps the accuracy in between the second and the third-order resulting in an excellent goal. In particular, the latter result is fully consistent with the nominal second-order of accuracy of the boundary interpolation method.

Figure 14: Surface pressure coefficient and its error scaling as a function of the grid resolution. Figure (14a) shows the pressure coefficient for a M = 2.81, Re = 1000 planar shock wave against a steady cylinder. Panel 14b reports the error scaling of the pressure coefficient. The data were collected at t/t0 = 1.5 for four uniform Cartesian meshes whose resolution was set ensuring a grid ratio r = 1.5. The case with Nx ×Ny = 180×180∗{1.54 } was considered as the reference solution. From the obtained results the overall accuracy of the entire algorithm is in between of the second and the third-order, consistently with the nominal second order of accuracy of the bilinear interpolation.

37

N 180 270 405 607

L1 0.827335E-01 0.471348E-01 0.144535E-01 0.683969E-02

L2 0.102169E+00 0.617132E-01 0.188992E-01 0.809366E-02

L∞ 0.190200E+00 0.138030E+00 0.441100E-01 0.146200E-01

L1 order 1.388 2.915 1.845

L2 order 1.243 2.919 2.092

Table 5: Error scaling of the present methodology for the surface pressure coefficient.

3.8. Shock wave diffraction against a wedge

Figure 15: Shock-wave-diffraction over a wedge at Ms = 1.3 and Res = 106 . In figure the numerical Schlieren density field is reported at two different time level and compared with the experimental results of Chang and Chang [69]. The simulation was performed over a uniform Cartesian grid featuring Nx × Ny = 2000 × 2000 nodes for a domain of [7L0 ; 7L0 ] where L0 is the triangle edge height. Three reference lengths were reserved upstream the wedge vertical edge.

680

685

To represent the dynamics of the shock wave diffraction against a steady wedge a similar configuration to the one explained in §3.6 was employed, thus the flow field was initialised in such a way that the Rankine-Hugoniot conditions ensure a travelling shock wave moving from left to right with a Mach number equals to 1.3. The numerical setting follows the experiments of Chang and Chang [69]. The left part of the domain consisted in a gas initially at rest with p/p0 = T /T0 = 1 while the Reynolds number of the flow Re = us L0 /νa was set equal to 106 . The choice of such high Reynolds number, in combination with no particular treatments for wall-turbulence model, makes the simulation practically inviscid in respect to its macroscopic behaviour. Here us denotes the shock speed, L0 the wedge edge height and νa the ambient 38

L∞ order 0.791 2.814 2.724

690

695

700

kinematic viscosity of the fluid. Extrapolation boundary conditions were ensured in all the four domains edges. The domain size was set equal to Lx × Ly = [7L0 ; 7L0 ] and the vertical edge of the wedge was located in such a way that three reference length would be reserved upstream. The mesh consisted in a structured uniform mesh featuring Nx × Ny = 2000 × 2000 nodes and the CFL number was set equal to 0.5. The computation was performed up to t/t0 = 5. In figure 15 the Schlieren density at two different time levels is reported. Here, as far as in the cylinder diffraction case, lots of discontinuities appear behind the wedge, detailing the complexity of the flow. In particular, three main features can be highlighted, consisting of two main triple points and a vortex-pair (a detailed description of the flow field can be found in [69]). Figure 16 reports the trajectories of those features compared to the one available in literature resulting in an excellent agreement also for this validation benchmark.

Figure 16: Shock-wave-diffraction over a wedge at Ms = 1.3 and Res = 106 . Literature comparison of the triple points trajectory and vortex centre location for the problem is reported. L0 is the wedge height while b = L0 tan(π/3) is the edge length. Present results are compared with the experimental data of Chang and Chang [69] and the numerical simulations of Chaudhuri et al. [70] and Boukharfane et al. [22].

39

705

710

715

720

725

3.9. Moving cylinder against a steady shock wave In order to show the performance of the present method with moving object, the shock wave diffraction over a cylinder (described in §3.6) is repeated in a different reference frame. Here we focus on the frame fixed with the shock so that in this configuration the cylinder is moving from left to right. Due to the Galileian invariance, the physics in the two frame is identical once the translational speed is accounted for. However, it should be remarked that the intrinsic non-linearity of WENO scheme does not guarantee a perfect Galileian invariance transformation. A similar set up as the previous paragraph was employed. A Cartesian grid featuring Nx × Ny = 700 × 400 grid points was employed while the domain sized Lx × Ly = [14D; 8D]. The Reynolds and the Mach numbers were set equal to 2000 and 2.81 and followed the same convections in §3.6. In Figure 17, the drag coefficient history is plotted in the two references. As we can notice, a very good agreement is found with the small differences barely noticeable. In Figure 18 two instantaneous fields in terms of density contours are shown comparing the steady with the moving cylinder concluding that our method is properly predictive also in case of moving objects in strong compressible conditions with shock waves. Here we want to highlight that the same simulation - setting k1 = 1 and k2 = 0 in the equation (8) - made the simulation to blow up, confirming the robustness improvement of a fully expanded formulation for convective derivatives.

Figure 17: Comparison between the time history of the drag coefficient in two different reference frames. The dashed line reports the results for a moving cylinder against a steady shock while the solid line shows the results for a moving shock against a steady cylinder. The cylinder (the shock) is travelling at M ach = 2.81, while the Reynolds number is set equals to 2000.

40

Figure 18: Density field in two homologous instants of a moving cylinder against a steady shock (left panel) and of a moving shock against a steady cylinder (right panel).

730

735

740

745

750

3.10. Supersonic flow past a confined moving square Finally, to conclude with a robustness benchmark for the entire methodology we select the supersonic flow past a moving confined square as a good and delicate benchmark of a compressible non-smooth moving object. The presented results provide reasonable qualitative behaviour since we did not find previous experimental data or numerical simulations to compare with. The simulation consisted in a M∞ = 3 and Re = 3000 flow over a single harmonic oscillating square along the crossflow direction. Slip-wall conditions were enforced at the top and at the bottom side of the domain while a supersonic inflow condition was enforced at the left side of the domain. The outlet bound was treated characteristically. The movement was analytical superimposed as y(t) = A sin(2πfe t) and here A = 0.2D is the amplitude associated to the harmonic oscillation, while fe is its effective oscillation frequency. The latter was set equal to 0.5f0 where f0 = St0 u∞ = 0.147u∞ represents a physical consistent value whose order of magnitude is the same of the natural shading frequency of a square in similar reference viscosity conditions (i.e. µ∞ ∼ M∞ /Re∞ ). The simulation was carried out over a uniform Cartesian mesh featuring Nx × Ny = [800 × 200] grid points for a domain of Lx × Ly = [32D; 8D]. Four reference lengths were reserved upstream the squared centre and the problem was run till a statistically steady convergence of the flow. The CFL number was set equal to 0.5. For this setting, a mean drag coefficient equal to c¯D = 1.62 has been recovered. The 41

instantaneous temperature fields in five different time levels are reported in figure 19. As we notice, the frontal shock is dynamically stretched following the square centre square.

Figure 19: Instantaneous temperature fields for a confined moving square at high-Mach number conditions.

4. Conclusions 755

760

The main result of the present study is that, using the fully-split-convective formulation of Kennedy-Gr¨ uber and Pirozzoli [30] in combination with IBM and a WENO procedure for shock-capturing, represents a suitable numerical strategy able to stably evolve moving objects in compressible flow. In particular, in order to reduce the typical noise oscillations introduced by the IBM filtering or upwinding techniques are usually used. We have shown that using the described numerical strategy it is possible to avoid any artificial viscosity injection even in the case of moving boundaries. The capability of 42

765

770

the method to accurately reproduce a large set of flow conditions, from subsonic laminar flows (e.g. flow over a confined and free objects at low Mach numbers), to transitional flows (e.g. vortex shedding downstream a cylinder and a square) up to high-speed shock-wave interactions has been quantitatively proved both in steady and moving cases reproducing well-documented datasets and showing very good agreements for the entire validation campaign. In particular, the low-Mach computations of moving cylinder and the results obtained with moving circle in strong compressible conditions have highlighted the vanishing level of noise introduced by the present approach without the need of any additional artificial dissipation. Hence the present methodology appears a promising strategy for accurate direct and large-eddy simulations of compressible turbulent flows with moving boundaries.

775

Appendix A. Sharp corner treatment A standard Immersed Boundary procedure can fail if the geometry is characterised by sharp corners because the bound normal direction associated with a ghost point is not unique. To deal with this singularities in the present work an ad-hoc numerical methodology has been employed. Being xg a ghost point; if xg has more then four fluid surrenders fluid the point is flagged as a corner-points. Consequently, a corner-region (Ωcorner ) is identified looking at the following domain n o  j+2 Ωcorner = (xi , yj ) ∈ Ωghost if (xm , yl )i+2 is a corner-point l=i−2 m=j−2 (A.1) and the information is stored inside a marker variable ζcorner defined as ( 1 if x ∈ Ωcorner ζcorner (x) = (A.2) 0 elsewhere The boundary interpolation method is slightly modified for the ghost points in the corner region. In particular we suppose that a two-dimensional sharp corner is characterised by two outer-normal directions, one for each Cartesian direction. In this way two boundary points (xbx , xby ) and two image points (xip x , xip y ) can be associated to the same ghost-corner-point (xg c ). Thus, the flow is bilinearly interpolated around both xip x and xip y and the value of 43

a general flow variable φ associated to a ghost-corner-point xg c is assumed equal to φip x ηx + φip y ηy (A.3) φ(xg c ) = ηx + ηy Here ηx and ηy are respectively the inverse of the square distances between xg c and the two image points xip x , xip y , thus ηx = 1/d2x and ηy = 1/d2y . Figure A.20 shows a sketch of the corner interpolation procedure.

Figure A.20: Sketch of the corner correction of the present methodology. The greenshaded region denotes the corner region Ωcorner while the two red points xip x , xip y are the image points associated to a ghost-corner-point xg c . 780

Appendix B. Quality improvement of a fully expanded formulation for convective derivative

785

790

In this appendix we discuss the quality improvement due to the local energy preservation character of the employed numerical discretisation. The higher robustness and the noise reduction due to a fully split formulation for convective derivative have been largely proved by many authors in recent publications (see e.g. [37, 30, 34, 35]). To asses this peculiarity an inviscid flow over a steady cylinder at M = 3.5 has been simulated with two different hybridisation strategies. The first wants the convective derivatives to be discretised with the WENO5+KGP6 method; the second seas the smooth regions of the flow to be solved by a standard sixth-order central finite difference method, thus setting k1 = 1 and k2 = 0 in the equation (8). The comparison 44

Figure B.21: Quality comparison of the non-dimensional vorticity contours (ωz ) for an inviscid steady cylinder at Mach 3.5. Panel B.21a shows the results obtained with the hybrid WENO5+KGP6 scheme, while panel in panel B.21b the results have been obtained setting k1 = 1 and k2 = 0, thus employing a standard sixth-order finite difference approximation for the convective fluxes in the smooth regions of the flow. Figures B.21c and B.21d report an enlargement of the wake region.

795

800

805

between the two methodologies highlights the characteristics of the energy preserving method in term of both cleanliness and time-stability. Looking at figure B.21 the non-dimensional vorticity contours of a steady M∞ = 3.5 cylinder are reported and plotted in a very narrow range of values. Panel B.21a shows the results computed with the hybrid WENO5+KGP6 scheme while panel B.21b reports the same simulation in which the local energy preservation character of the scheme has been de-activated. As we notice, the macroscopic behaviour of the flow is the same but the energy preserving method denotes a much cleaner results in term of acoustical waves and wake prediction. In particular, the flow behind the cylinder in the second panel is characterised by spurious oscillations in the region were the central-scheme is activated. Moreover, the core of the wake is very noisy compared to the simulation in B.21a. References [1] Robert C. Blanchard and Prasun N. Desai. Mars Phoenix Entry, Descent, and Landing Trajestory and Atmosphere Reconstruction. Journal of Spacecraft and Rockets, 2012. ISSN 0022-4650. doi: 10.2514/1.46274. ¨ ur Karatekin, Stephen R. Lewis, Fran¸cois Forget, [2] Francesca Ferri, Ozg¨ 45

Alessio Aboudan, Giacomo Colombatti, Carlo Bettanini, Stefano Debei, Bart Van Hove, Veronique Dehant, Ari Matti Harri, Mark Leese, Teemu M¨akinen, Ehouarn Millour, Ingo Muller-Wodarg, Gian Gabriele Ori, Andrea Pacifici, Sebastien Paris, Manish Patel, Mark Schoenenberger, Jeffrey Herath, Tero Siili, Aymeric Spiga, Tetsuya Tokano, Martin Towner, Paul Withers, Sami Asmar, and Dirk Plettemeier. ExoMars Atmospheric Mars Entry and Landing Investigations and Analysis (AMELIA). Space Science Reviews, 215(1):1–21, 2019. ISSN 15729672. doi: 10.1007/s11214-019-0578-x. URL http://dx.doi.org/10.1007/ s11214-019-0578-x. [3] Nicol´as Garc´ıa Rosa, Guillaume Dufour, Roger Bar`enes, and G´erard Lavergne. Experimental Analysis of the Global Performance and the Flow Through a High-Bypass Turbofan in Windmilling Conditions. Journal of Turbomachinery, 137(5): 051001, 2015. ISSN 0889-504X. doi: 10.1115/1.4028647. URL http://turbomachinery.asmedigitalcollection.asme.org/ article.aspx?doi=10.1115/1.4028647. [4] Charles S. Peskin. Flow patterns around heart valves: A numerical method. Journal of Computational Physics, 10(2):252–271, 1972. ISSN 10902716. doi: 10.1016/0021-9991(72)90065-4. [5] 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(1): 35–60, 2000. ISSN 00219991. doi: 10.1006/jcph.2000.6484. [6] Gianluca Iaccarino and Roberto Verzicco. Immersed boundary technique for turbulent flow simulations. Applied Mechanics Reviews, 2003. ISSN 00036900. doi: 10.1115/1.1563627. [7] Rajat Mittal and Gianluca Iaccarino. Immersed Boundary Methods. Annual Review of Fluid Mechanics, 37(1):239–261, 2005. ISSN 0066-4189. doi: 10.1146/annurev.fluid.37.061903.175743. URL http://www.annualreviews.org/doi/10.1146/annurev.fluid.37. 061903.175743. [8] Amir Eshghinejadfard, Abouelmagd Abdelsamie, Seyed Ali Hosseini, and Dominique Th´evenin. Immersed boundary lattice Boltzmann sim46

ulation of turbulent channel flows in the presence of spherical particles. International Journal of Multiphase Flow, 2017. ISSN 03019322. doi: 10.1016/j.ijmultiphaseflow.2017.07.011. [9] Markus Uhlmann. An immersed boundary method with direct forcing for the simulation of particulate flows. Journal of Computational Physics, 209(2):448–476, 2005. ISSN 00219991. doi: 10.1016/j.jcp.2005. 03.017. [10] Wim Paul Breugem. A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. Journal of Computational Physics, 231(13):4469–4498, 2012. ISSN 00219991. doi: 10.1016/j.jcp.2012.02.026. URL http://dx.doi.org/10.1016/j. jcp.2012.02.026. [11] Tobias Kempe and Jochen Fr¨ohlich. An improved immersed boundary method with direct forcing for the simulation of particle laden flows. Journal of Computational Physics, 2012. ISSN 10902716. doi: 10.1016/ j.jcp.2012.01.021. [12] Francesco Picano, Wim Paul Breugem, and Luca Brandt. Turbulent channel flow of dense suspensions of neutrally buoyant spheres. Journal of Fluid Mechanics, 764:463–487, 2015. ISSN 14697645. doi: 10.1017/ jfm.2014.704. [13] S. Schwarz, T. Kempe, and J. Fr¨ohlich. A temporal discretization scheme to compute the motion of light particles in viscous flows by an immersed boundary method. Journal of Computational Physics, 2015. ISSN 10902716. doi: 10.1016/j.jcp.2014.10.039. [14] A. Eshghinejadfard and D. Th´evenin. Numerical simulation of heat transfer in particulate flows using a thermal immersed boundary lattice Boltzmann method. International Journal of Heat and Fluid Flow, 2016. ISSN 0142727X. doi: 10.1016/j.ijheatfluidflow.2016.04.002. [15] Mohd Hazmil Abdol Azis, Fabien Evrard, and Berend van Wachem. An immersed boundary method for flows with dense particle suspensions. Acta Mechanica, 2019. ISSN 00015970. doi: 10.1007/s00707-018-2296-y.

47

[16] Weiwei Ren, Chang Shu, and Wenming Yang. An efficient immersed boundary method for thermal flow problems with heat flux boundary conditions. International Journal of Heat and Mass Transfer, 64: 694–705, 2013. ISSN 00179310. doi: 10.1016/j.ijheatmasstransfer.2013. 05.020. URL http://dx.doi.org/10.1016/j.ijheatmasstransfer. 2013.05.020. [17] Y. Wang, C. Shu, and L. M. Yang. Boundary condition-enforced immersed boundary-lattice Boltzmann flux solver for thermal flows with Neumann boundary conditions. Journal of Computational Physics, 306: 237–252, 2016. ISSN 10902716. doi: 10.1016/j.jcp.2015.11.046. URL http://dx.doi.org/10.1016/j.jcp.2015.11.046. [18] Cindy Merlin, Pascale Domingo, and Luc Vervisch. Immersed boundaries in large eddy simulation of compressible flows. Flow, Turbulence and Combustion, 90(1):29–68, 2013. ISSN 13866184. doi: 10.1007/s10494-012-9421-0. [19] Kun Luo, Zhenya Zhuang, Jianren Fan, and Nils Erland L. Haugen. A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions. International Journal of Heat and Mass Transfer, 92:708–717, 2015. ISSN 00179310. doi: 10.1016/j.ijheatmasstransfer.2015.09.024. URL http://dx.doi. org/10.1016/j.ijheatmasstransfer.2015.09.024. [20] A. Piquet, O. Roussel, and A. Hadjadj. A comparative study of Brinkman penalization and direct-forcing immersed boundary methods for compressible viscous flows. Computers and Fluids, 136:272– 284, 2016. ISSN 00457930. doi: 10.1016/j.compfluid.2016.06.001. URL http://dx.doi.org/10.1016/j.compfluid.2016.06.001. [21] Matteo Bernardini, Davide Modesti, and Sergio Pirozzoli. On the suitability of the immersed boundary method for the simulation of high-Reynolds-number separated turbulent flows. Computers & Fluids, 130:84–93, may 2016. ISSN 00457930. doi: 10.1016/j.compfluid.2016.02.018. URL http://dx.doi.org/10.1016/j.compfluid.2016.02.018https: //linkinghub.elsevier.com/retrieve/pii/S0045793016300421.

48

[22] Radouan Boukharfane, F´abio Henrique Eugˆenio Ribeiro, Zakaria Bouali, and Arnaud Mura. A combined ghost-point-forcing / direct-forcing immersed boundary method (IBM) for compressible flow simulations. Computers and Fluids, 162:91–112, 2018. ISSN 00457930. doi: 10.1016/j.compfluid.2017.11.018. URL https://doi.org/10.1016/j. compfluid.2017.11.018. [23] E. Guilmineau and P. Queutey. A numerical simulation of vortex shedding from an oscillating circular cylinder. Journal of Fluids and Structures, 16(6):773–794, 2002. ISSN 08899746. doi: 10.1006/jfls.2002.0449. [24] Lennart Schneiders, Daniel Hartmann, Matthias Meinke, and Wolfgang Schr¨oder. An accurate moving boundary formulation in cut-cell methods. Journal of Computational Physics, 235:786–809, 2013. ISSN 00219991. doi: 10.1016/j.jcp.2012.09.038. URL http://dx.doi.org/ 10.1016/j.jcp.2012.09.038. [25] Yegao Qu, Ruchao Shi, and Romesh C. Batra. An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids. Journal of Computational Physics, 354:672–691, 2018. ISSN 10902716. doi: 10.1016/j.jcp.2017.10.045. URL https://doi. org/10.1016/j.jcp.2017.10.045. [26] G. B. Deng, J. Piquet, P. Queutey, and M. Visonneau. Incompressible flow calculations with a consistent physical interpolation finite volume approach. Computers and Fluids, 23(8):1029–1047, 1994. ISSN 00457930. doi: 10.1016/0045-7930(94)90003-5. [27] Jianming Yang and Elias Balaras. An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries. Journal of Computational Physics, 215(1):12–40, 2006. ISSN 00219991. doi: 10.1016/j.jcp.2005.10.035. [28] Cheng Chi, Abouelmagd Abdelsamie, and Dominique Th´evenin. A directional ghost-cell immersed boundary method for incompressible flows. Journal of Computational Physics, page 109122, 2019. ISSN 0021-9991. doi: 10.1016/j.jcp.2019.109122. URL https://doi.org/10.1016/j. jcp.2019.109122.

49

[29] M. Ehsan Khalili, Martin Larsson, and Bernhard M¨ uller. Immersed boundary method for viscous compressible flows around moving bodies. Computers and Fluids, 170:77–92, 2018. ISSN 00457930. doi: 10.1016/j.compfluid.2018.04.033. URL https://doi.org/10.1016/j. compfluid.2018.04.033. [30] Sergio Pirozzoli. Generalized conservative approximations of split convective derivative operators. Journal of Computational Physics, 229(19): 7180–7190, 2010. ISSN 00219991. doi: 10.1016/j.jcp.2010.06.006. URL http://dx.doi.org/10.1016/j.jcp.2010.06.006. [31] Sergio Pirozzoli. Stabilized non-dissipative approximations of Euler equations in generalized curvilinear coordinates. Journal of Computational Physics, 2011. ISSN 10902716. doi: 10.1016/j.jcp.2011.01.001. [32] Francesco Salvadore, Matteo Bernardini, and Michela Botti. GPU accelerated flow solver for direct numerical simulation of turbulent flows. Journal of Computational Physics, 235:129–142, 2013. ISSN 00219991. doi: 10.1016/j.jcp.2012.10.012. URL http://dx.doi.org/10.1016/j. jcp.2012.10.012. [33] Davide Modesti and Sergio Pirozzoli. Reynolds and Mach number effects in compressible turbulent channel flow. International Journal of Heat and Fluid Flow, 59:33–49, 2016. ISSN 0142727X. doi: 10.1016/j.ijheatfluidflow.2016.01.007. URL http://dx.doi.org/10. 1016/j.ijheatfluidflow.2016.01.007. [34] Davide Modesti and Sergio Pirozzoli. Direct numerical simulation of supersonic pipe flow at moderate Reynolds number. International Journal of Heat and Fluid Flow, 2019. ISSN 0142727X. doi: 10.1016/j. ijheatfluidflow.2019.02.001. [35] G. Coppola, F. Capuano, S. Pirozzoli, and L. de Luca. Numerically stable formulations of convective terms for turbulent compressible flows. Journal of Computational Physics, 382:86–104, 2019. ISSN 10902716. doi: 10.1016/j.jcp.2019.01.007. URL https://doi.org/10.1016/j. jcp.2019.01.007. [36] G S Jiang and C W Shu. Efficient implementation of weighted WENO schemes. Journal of Computational Physics, 126(126):202–228, 1996. 50

[37] Yohei Morinishi. Skew-symmetric form of convective terms and fully conservative finite difference schemes for variable density low-Mach number flows. Journal of Computational Physics, 2010. ISSN 10902716. doi: 10.1016/j.jcp.2009.09.021. [38] Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R. Chakravarthy. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of Computational Physics, 1987. ISSN 10902716. doi: 10.1016/0021-9991(87)90031-3. [39] Dinshaw S. Balsara and Chi Wang Shu. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics, 160(2):405–452, 2000. ISSN 00219991. doi: 10.1006/jcph.2000.6443. [40] Dinshaw S. Balsara, Sudip Garain, and Chi Wang Shu. An efficient class of WENO schemes with adaptive order. Journal of Computational Physics, 326:780–804, 2016. ISSN 10902716. doi: 10.1016/j.jcp.2016.09. 009. URL http://dx.doi.org/10.1016/j.jcp.2016.09.009. [41] Marcos Castro, Bruno Costa, and Wai Sun. High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws. Journal of Computational Physics, 230(5):1766–1792, 2011. ISSN 0021-9991. doi: 10.1016/j.jcp.2010.11.028. URL http://dx.doi.org/ 10.1016/j.jcp.2010.11.028. [42] Culbert B. Laney. Computational Gasdynamics. Cambridge University Press, 1998. ISBN 9780521570695. [43] A. Rohde. Eigenvalues and eigenvectors of the Euler equations in general geometries. AIAA Paper 2001-2609, (1):1–6, 2001. ISSN 1538-3598. doi: doi:10.2514/6.2001-2609. [44] Eleuterio F. Toro. Riemann solvers and numerical methods for fluid dynamics: A practical introduction. Springer, 2009. ISBN 9783540252023. doi: 10.1007/b79761. [45] F. Ducros. Large-eddy simulation of the shock/turbulence interaction. Journal of Computational Physics, 152:517—-549, 1999. ISSN 00011452. doi: 10.2514/1.6002. 51

[46] Sergio Pirozzoli. Numerical Methods for High-Speed Flows. Annual Review of Fluid Mechanics, 43(1):163–194, 2011. ISSN 00664189. doi: 10.1146/annurev-fluid-122109-160718. URL http://www. annualreviews.org/doi/10.1146/annurev-fluid-122109-160718. [47] Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing RungeKutta schemes. Mathematics of Computation of the American Mathematical Society, 67(221):73–85, 1998. ISSN 0025-5718. doi: 10. 1090/S0025-5718-98-00913-2. URL http://www.ams.org/jourcgi/ jour-getitem?pii=S0025-5718-98-00913-2. [48] Joseph O’Rourke. Computational Geometry in C., volume 64. Cambridge University Press, 1998. doi: 10.2307/2153463. [49] Xiaolei Yang, Xing Zhang, Zhilin Li, and Guo Wei He. A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations. Journal of Computational Physics, 228(20):7821–7836, 2009. ISSN 00219991. doi: 10.1016/j.jcp. 2009.07.023. URL http://dx.doi.org/10.1016/j.jcp.2009.07.023. [50] T.J. Poinsot. Boundary conditions for direct simulations of compressible viscous flows. Journal of Computational Physics, 99(2):352, 2004. ISSN 00219991. doi: 10.1016/0021-9991(92)90227-p. [51] Guido Lodato, Pascale Domingo, and Luc Vervisch. Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows. Journal of Computational Physics, 2008. ISSN 10902716. doi: 10.1016/j.jcp.2008.01.038. [52] Sergio Pirozzoli and Tim Colonius. Generalized characteristic relaxation boundary conditions for unsteady compressible flow simulations. Journal of Computational Physics, 248:109–126, 2013. ISSN 00219991. doi: 10.1016/j.jcp.2013.04.021. URL http://dx.doi.org/10.1016/j.jcp. 2013.04.021. [53] Mehmet Sahin and Robert G. Owens. A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder. Physics of Fluids, 16(5):1305–1320, 2004. ISSN 10706631. doi: 10.1063/1.1668285.

52

[54] S. C. R. Dennis and GAU-ZU Chang. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. Journal of Fluid Mechanics, 42(3):471–489, jul 1970. ISSN 0022-1120. doi: 10.1017/S0022112070001428. URL http://link.springer.com/ 10.1007/BF02942594https://www.cambridge.org/core/product/ identifier/S0022112070001428/type/journal{_}article. [55] D. J. Tritton. Experiments on the flow past a circular cylinder at low Reynolds numbers. Journal of Fluid Mechanics, 6(04):547, 2006. ISSN 0022-1120. doi: 10.1017/s0022112059000829. [56] P. De Palma, M. D. de Tullio, G. Pascazio, and M. Napolitano. An immersed-boundary method for compressible viscous flows. Computers and Fluids, 35(7):693–702, 2006. ISSN 00457930. doi: 10.1016/j. compfluid.2006.01.004. [57] Mark Linnick and Hermann Fasel. A High-Order Immersed Boundary Method for Unsteady Incompressible Flow Calculations. Number January, pages 1–17, 2012. doi: 10.2514/6.2003-1124. [58] Lennart Schneiders, Daniel Hartmann, Matthias Meinke, and Wolfgang Schr¨oder. An accurate moving boundary formulation in cut-cell methods. Journal of Computational Physics, 235:786–809, 2013. ISSN 00219991. doi: 10.1016/j.jcp.2012.09.038. URL http://dx.doi.org/ 10.1016/j.jcp.2012.09.038. [59] Reza Ghias, Rajat Mittal, and Thomas Lund. A Non-Body Conformal Grid Method for Simulation of Compressible Flows with Complex Immersed Boundaries. In 42nd AIAA Aerospace Sciences Meeting and Exhibit, number January, pages 1–10. American Institute of Aeronautics and Astronautics, 2013. doi: 10. 2514/6.2004-80. URL https://pdfs.semanticscholar.org/5baf/ 145473caf7cc495442c48e189ef41b18721b.pdf. [60] B. N. Rajani, A. Kandasamy, and Sekhar Majumdar. Numerical simulation of laminar flow past a circular cylinder. Applied Mathematical Modelling, 33(3):1228–1247, 2009. ISSN 0307904X. doi: 10.1016/j.apm. 2008.01.017. URL http://dx.doi.org/10.1016/j.apm.2008.01.017.

53

[61] M. Breuer, J. Bernsdorf, T. Zeiser, and F. Durst. Accurate computations of the laminar flow past a square cylinder based on two different methods: Lattice-Boltzmann and finite-volume. International Journal of Heat and Fluid Flow, 21(2):186–196, 2000. ISSN 0142727X. doi: 10.1016/S0142-727X(99)00081-8. [62] A. Sohankar, C. Norberg, and L. Davidson. Low-Reynolds-number flow around a square cylinder at incidence: Study of blockage, onset of vortex shedding and outlet boundary condition. Technical report, 1998. [63] Rupad M. Darekar and Spencer J. Sherwin. Flow past a square-section cylinder with a wavy stagnation face. Journal of Fluid Mechanics, 2001. ISSN 00221120. doi: 10.1017/S0022112000002299. [64] A. P. Singh, A. K. De, V. K. Carpenter, V. Eswaran, and K. Muralidhar. Flow past a transversely oscillating square cylinder in free stream at low reynolds numbers. International Journal for Numerical Methods in Fluids, 2009. ISSN 02712091. doi: 10.1002/fld.1979. [65] Akhilesh K. Sahu, R. P. Chhabra, and V. Eswaran. Two-dimensional laminar flow of a power-law fluid across a confined square cylinder. Journal of Non-Newtonian Fluid Mechanics, 2010. ISSN 03770257. doi: 10.1016/j.jnnfm.2010.03.011. [66] Subhankar Sen, Sanjay Mittal, and Gautam Biswas. Flow past a square cylinder at low Reynolds numbers. International Journal for Numerical Methods in Fluids, 2011. ISSN 02712091. doi: 10.1002/fld.2416. [67] W. Gu, C. Chyu, and D. Rockwell. Timing of vortex formation from an oscillating cylinder. Physics of Fluids, 6(11):3677–3682, 1994. ISSN 10706631. doi: 10.1063/1.868424. [68] A. E. Bryson and R. W. F. Gross. Diffraction of strong shocks by cones, cylinders, and spheres. Journal of Fluid Mechanics, 10(1): 1–16, feb 1961. ISSN 0022-1120. doi: 10.1017/S0022112061000019. URL https://www.cambridge.org/core/product/identifier/ S0022112061000019/type/journal{_}article. [69] Se Myong Chang and Keun Shik Chang. On the shock-vortex interaction in Schardin’s problem. Shock Waves, 10(5):333–343, 2000. ISSN 09381287. doi: 10.1007/s001930000061. 54

[70] Arnab Chaudhuri, Abdellah Hadjadj, and Ashwin Chinnayya. On the use of immersed boundary methods for shock/obstacle interactions. Journal of Computational Physics, 230(5):1731–1748, 2011. ISSN 00219991. doi: 10.1016/j.jcp.2010.11.016. URL http://dx.doi.org/ 10.1016/j.jcp.2010.11.016.

55

Declaration of Competing Interest The authors declare that they have not any conflicts of interest.

56

CRediT author statement Francesco De Vanna: Software, Methodology, Investigation, Formal analysis, Writing-Reviewing & Editing. Francesco Picano: Conceptualization, Investigation, Supervision, Project Administration, Writing-Reviewing & Editing. Ernesto Benini: Supervision, Project administration, Writing-Reviewing & Editing.

57