Numerical simulation of fluid forces on moving solid body by the vortex in cell method with volume penalization

Aerospace Science and Technology 94 (2019) 105360

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Numerical simulation of fluid forces on moving solid body by the vortex in cell method with volume penalization Tomohiro Degawa a , Qiang Gu b , Tomomi Uchiyama a,∗ , Kotaro Takamure a a b

Institute of Materials and Systems for Sustainability, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan Graduate School of Informatics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan

a r t i c l e

i n f o

Article history: Received 11 February 2019 Received in revised form 3 June 2019 Accepted 23 August 2019 Available online 29 August 2019 Keywords: CFD Fluid force Vortex in cell method Volume penalization method Moving body

a b s t r a c t This study proposes a method for calculating the fluid forces acting on moving solid bodies in flows simulated by the vortex in cell (VIC) method. The body is represented using the volume penalization (VP) method, which introduces the influence of the body on the flow as an external force in the form of a penalization term into the Navier-Stokes equation. The pressure field is computed by solving the Poisson equation, which is derived by taking the divergence of the Navier-Stokes equation. The pressure and viscous stress on the body surface are calculated using a quadratic function along the vector normal to the surface. This study simulates the flow around a cylinder oscillating in still water by the VICVP method, and calculates the fluid forces acting on the cylinder using the proposed method. The velocity profiles on four cross-sections at three different time points agree well with the results of a measurement. The in-line force acting on the cylinder also agrees with the existing results computed by a boundary-conforming formulation. These demonstrate the validity of the authors’ method for calculating the fluid forces on solid bodies in flows simulated by the VIC-VP method. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction For the analysis and prediction of the performance of turbomachinery, the numerical simulations of the flow through rotating rotors and blades have been conducted [1–4]. As turbomachinery have complex shape, simulation methods that can handle complex geometries, such as finite element methods and finite volume methods, are frequently used. These methods employ computational meshes constructed from computer aided design data, so the validity of the simulation results is strongly affected by the mesh quality. In addition, reconstruction of the mesh data sometimes requires reflecting the rotating motion of the rotor and blade in the mesh, which increases the computational load and computation time. To overcome the aforementioned limitations, immersed boundary methods, which can represent a solid body on Cartesian grids, have been proposed and are used extensively. An immersed boundary method was initially developed to simulate blood flows through a beating heart [5]. Subsequently, many researchers have improved and presented variants of the immersed boundary method [6]. Immersed boundary methods are used to

*

Corresponding author. E-mail address: [email protected] (T. Uchiyama).

https://doi.org/10.1016/j.ast.2019.105360 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

simulate many types of flows [7–12], including turbulent flows [13–16]. Immersed boundary methods are classified into two types, the direct forcing method and the feedback forcing method, based on the procedure used to compute fluid velocities on the surface of solid bodies. The direct forcing method computes the fluid velocity in a cell containing a body surface by using the velocity of the solid surface and the velocity of the fluid at nearby grid points. Although such treatment of the body surface is intuitively understandable, it is necessary to switch the treatment of velocities in the cells containing the solid body and the cells containing fluid only. The feedback forcing method considers the existence of bodies as external forces, and it adds an external force term to the governing equations. The external forces are determined such that they satisfy boundary conditions of velocity on the solid surface. The volume penalization (VP) method [17], a feedback forcing method, can simply represent solid bodies in fluids, and it has been used in conjunction with higher-order numerical methods [18,19]. The properties of the method have been investigated as well [20,21]. The aforementioned immersed boundary methods have been proposed for the Navier-Stokes equation in which velocity and pressure are primitive variables. Meanwhile, the vortex in cell (VIC)

2

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

method has been used successfully to simulate high-Reynoldsnumber flows. The VIC method is a vortex method [22], and it discretizes the vorticity of the flow field into vortex elements and traces the convection of each vortex element to simulate the time evolution of the flows. The VIC method solves the Poisson equation of the stream-function or vector potential on a spatial grid to determine the convection of vortex elements. Because the VIC method computes time evolutions of the vorticity of the fluid field, the time evolution of velocity in this method depends on the evolution of vorticity. In addition, pressure does not appear in the vorticity equation, and therefore, it is impossible to apply the direct forcing method to correct the velocity and pressure of the flow for representing solid bodies. A few immersed boundary methods [23–25] have been proposed for the vorticity equation, in which vorticity and stream-function are primitive variables. By contrast, the feedback forcing method can easily be combined with the vorticity equation. VIC methods combined with the boundary element method [26] or the volume penalization (VP) method [27–30] are proposed to simulate the flows around complex or moving bodies. In these methods, the fluid forces acting on the body are calculated according to a method based on momentum exchange. However, this method cannot estimate individually the pressure and the viscous components of the fluid forces acting on the bodies. Herein, the present authors propose a method for calculating the fluid forces by using the VIC method combined with the VP method for numerical simulation of the flow around moving bodies. In this method, the spatial distribution of pressure is computed by solving the pressure Poisson equation obtained by taking the divergence of the Navier-Stokes equation; then, the pressure and the viscous stress acting on the body surface are estimated from the pressure at nearby grid points and tangential velocity near the body surface. The proposed method is applied to simulate the fluid force acting on a cylinder oscillating in still water. The simulated flow around the oscillating cylinder and fluid forces on the cylinder agree with those of the existing studies [15,31]. These demonstrate the validity of the present method. 2. Governing equation and numerical method 2.1. Governing equation The mass and momentum conservation equations for twodimensional incompressible flow in the x– y plane are as follows:

∇ ·u = 0

(1)

∂u 1 + (u · ∇) u = − ∇ p + ν ∇ 2 u + λχ (u s − u ) ∂t ρ

(2)

χ=

⎪ ⎪ ⎩

1 2

l ≤ −1  ε



1 + εl + π1 sin π εl

εl < 1 l 1 ε ≥1



0



where l is the distance from the solid body surface, and ε is a parameter determining the slope of the function at the surface. ε is given as 1% of the characteristic length [28]. Taking the curl of Eq. (2) and substituting Eq. (1) into resulting equation, the following vorticity equation is derived:

∂ω + ∇ · (u ω) = ν ∇ 2 ω + λ∇ × χ (u s − u ) (4) ∂t where ω is the vorticity around the z-axis perpendicular to the x– y plane. According to the Helmholtz theorem, the velocity u is the sum of the curl of a vector potential ψ and the gradient of a scalar potential φ :

u = ∇ × ψ + ∇φ

(3)

(5)

The velocity calculated from Eq. (5) remains unaltered even if any gradient of a scalar function is added to ψ . To eliminate this arbitrariness, the following solenoidal condition is imposed on ψ ,

∇ ·ψ =0

(6)

Taking the curl of Eq. (5) and substituting Eq. (6) into the resulting equation, the following Poisson equation for ψ is derived:

∇ 2 ψ = −ω

(7)

Because the present study focuses only on two-dimensional flow, the z-component of the vector potential, ψz , is considered. ψz corresponds to the stream-function. Hereinafter, ψz is expressed as ψ. According to the vorticity equation, it is not necessary to compute pressure through the time evolution of the vorticity field. Hence, pressure is computed independently from the velocity field. Taking the divergence of Eq. (2) and substituting Eq. (1) into the resulting equation, the following pressure Poisson equation is obtained:

1

− ∇ 2 p = ∇ · [(u · ∇) u] − λ∇ · [χ (u s − u )]

ρ

where u is the velocity, p is the pressure, and ρ and ν are the density and the kinematic viscosity, respectively. The final term on the right-hand side of Eq. (2) is the penalization term for non-slip conditions on the surface of solid bodies. λ, χ , and u s are the penalization parameter, the mask function, and the velocity of the solid body, respectively. The mask function χ in Eq. (2) is a step function, which is 0 within the fluid and 1 inside the solid body. This discontinuity causes numerical oscillations of the simulated flow. To suppress the non-physical oscillations, the following sigmoid-type mask function [28] is used to guarantee continuity:

⎧ ⎪ ⎪ ⎨

Fig. 1. Definition of location of variable on staggered grid.

(8)

2.2. Computational grid A regular grid, which defines all physical variables on the same grid points, is usually used to divide the computational domain into the spatial grids in the finite difference method described by the stream-function-vorticity formulation, and it is used in the VIC method as well. One of the present authors [32] successfully clarified that dividing the computational domain by a staggered grid, which defines the physical variables at different points, can ensure consistency among the discretized equations and can prevent numerical oscillation of the flow field. In this study, a staggered grid is used for spatial discretization of the computational domain. As shown in Fig. 1, the vorticity and the stream-function are defined at each grid point. The velocity components are defined at the center of the cell sides. Scalar variables, including pressure and mask function, are defined at the cell center.

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

δu δ v + =0 δx δ y

2.3. Vortex in cell method It is postulated that the position vector and vorticity of the vortex element p are x p = x p , y p and ω p , respectively. The Lagrangian form of the vorticity equation, Eq. (4), is written as follows:

 ∂ xp = u xp ∂t     ∂ω = ν ∇ 2 ω x p + λ∇ × χ u s x p − u x p ∂t

(9) (10)

When the position and vorticity of a vortex element are known at time t, the values at t + t can be computed from the Lagrangian calculation  of Eqs. (9) and (10). If ω is defined at a position x g = x g , y g , the vorticity ω is assigned to x g , or a vortex element with vorticity ω is redistributed onto x g or the vortex element with vorticity ω is redistributed onto x g :

ω=

Nν



ωp W

p

xg − xp



x

 W

yg − yp



y

(11)

where N ν is the number of vortex elements, and x and y are the width of a cell. The following B-spline function [22] is used as the redistribution function W (σ ).

⎧ 3 2 |σ | ≤ 1 ⎪ ⎨ 1 − 2.5σ + 1.5 |σ | 2 W (σ ) = 0.5 (2 − |σ |) (1 − |σ |) 1 ≤ σ ≤ 2 ⎪ ⎩ |σ | ≥ 2 0

(12)

2.4. Discretization on staggered grids



δ δω δ δω + δx δx δy δy



 +λ

(16)

into the resulting equation for simplification, the following formulas of the convection term and the penalization term are obtained:

    x y x y δ δ x δu x δu y δv y δv u +v u +v + δx δx δy δy δx δy δu δu δx δx

xx

xy

yy

δ v δu δv δv + δx δ y δy δy   δ  y δ  x χ (u s − u ) + χ (v s − v ) δx δy   x δus δ v s δχ δχ + (u s − u ) + + (v s − v ) =χ δx δy δx δy =

+2

(17)

y

(18)

For these reformulations, an identical equation introduced by Kajishima [34] is used. The time derivative and viscous diffusion terms become 0. The discretized form of the pressure Poisson equation is derived as follows:

 xx xy yy δ δp δu δu δ v δu δv δv δ δp = − + +2 + ρ δx δx δ y δ y δx δx δx δ y δy δy     x y δus δ v s δχ δχ + (u s − u ) −λ χ + + (v s − v ) δx δy δx δy 1



(19)

where ·xx and ·y y represent the interpolation operations twice in the x and y directions, respectively. Similarly, ·xy represents an interpolation operation in the x and y directions. 2.5. Method for calculating fluid forces acting on solid bodies

Equation (10) is discretized on the staggered grid. When the second-order central-finite-difference method is used for spatial difference calculations, Eq. (10) is discretized as follows:

∂ω =ν ∂t

3

δ χ y ( v s − v ) δ χ x (u s − u ) − δx δy



(13) where δ · /δ x, δ · /δ y and ·x , ·y are the differencing and interpolation operators [33], respectively. Here, the subscript representing the stencil width is omitted because it is always 1 when second-order differencing and interpolation are performed. The time derivative term in Eq. (13) is discretized using a time integration method, such as the forward Euler method, and then Eq. (13) is solved after redistributing the vorticity of vortex elements. There is no appropriate discretized form of the pressure Poisson equation (8) on the staggered grid. In this study, a discretized form of the pressure Poisson equation was obtained by taking the divergence of Eqs. (14) and (15), which are discretized appropriately on the staggered grid.

  x y δ δu ∂u δu δu 1 δp δ δu + ux + vx =− +ν + ∂t δx δy ρ δx δx δx δ y δ y + λχ x ( u s − u )   x y δ δv ∂v 1 δp δ δv y δv y δv +u +v =− +ν + ∂t δx δy ρ δy δx δx δ y δ y + λχ y ( v s − v )

(14)

(15)

The differencing operators δ ·/δ x and δ ·/δ y are applied to Eqs. (14) and (15), respectively, and then adding those equations. Substituting the mass conservation equation discretized on the staggered grid

After the pressure field is obtained, the pressure and velocity on the solid body are estimated to calculate the fluid forces acting on the body. The fluid forces are calculated using the following equation:



F=

   − p S − p  n + τ · n dS

(20)

ΓS

where p S is the pressure on the body surface, p  is a correction term when the integral surface Γ S is not a closed surface, and n is the outward unit vector normal to the surface. τ is the viscous tensor given as:



τ = μ ∇ u + (∇ u )T



(21)

where the superscript T represents the transpose of matrix. Because the surface of the bodies does not coincide with the defined point of pressure in almost all cases, the value of pressure at the surface is estimated from nearby values by means of interpolation. To compute the pressure on the surface, p S , an interpolation polynomial is constructed [12].

p (l) − p  = αl2 + β l + γ

(22)

where l is the distance from the body surface in the normal direction. As shown in Fig. 2, two points I1 , I2 which are at lengths l1 , l2 away from the surface are defined in the normal direction to calculate the parameters α , β , and γ . The pressure on the surface is obtained at l = 0 by using Eq. (22). The parameter β is defined as follows:

β=

∂ p

∂ n l→0

(23)

4

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

and Yang and Balaras [15] calculated β from the inner product of the normal vector n and the Navier-Stokes equation without the viscous term. According to Yang and Balaras [15], in the present study, β can be calculated as follows:

  ∂u ∂ p

·n = − ρ + u · ∇) u − λ χ u − u ( ) ( s ∂ n l→0 ∂t

(24)

where the diffusion term disappears because the velocities around the defined point of pressure are calculated by a linear interpolation method. After obtaining β , the other two parameters, α and γ , are calculated by solving the following simultaneous linear equations



p (l1 ) − p  = αl21 + β l1 + γ p (l2 ) − p  = αl22 + β l2 + γ

Fig. 2. Interpolation point of pressure and shear force.

where p (l1 ) and p (l2 ) are the pressures at I1 and I2 , respectively, and they are calculated from nearby values by means of linear interpolation. For the calculation of shear force, a Lagrangian polynomial constructed from the tangential velocity at points I1 and I2 was proposed [35]. This polynomial assumes that the velocity at the body surface (corresponding to point I in Fig. 2) is 0. In the present research, a new Lagrangian polynomial is constructed to consider the velocity of bodies. The modified polynomial is

u T (l1 ) l22 − u T (l2 ) l21 − (l1 + l2 ) (l2 − l1 ) u T (0) ∂uT = ∂n l1l2 (l2 − l1 )

Fig. 3. Computational domain.

(25)

(26)

where u T is the tangential velocity. If I1 and I2 are placed at equal intervals (l2 = 2l1 ), Eq. (26) coincides with the second-order onesided difference formulation. When the velocity at a certain time step n is known, the fluid force acting on bodies is calculated according to the following procedure.

Fig. 4. Streamlines and velocity distribution in the fourth period.

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

5

1. Determine a point I and calculate the pressures p (l1 ) and p (l2 ) at I1 and I2 , respectively, by means of linear interpolation. 2. Calculate the surface pressure from Eq. (22) using parameters α , β , and γ computed from Eqs. (23), (24), and (25). 3. Calculate the tangential velocity u T (l1 ) and u T (l2 ) by means of linear interpolation. 4. Calculate the shear stress on the surface from Eq. (26). 5. Calculate fluid forces acting on bodies by numerical integration of Eq. (20) over the entire integral surface Γ S . 3. Application to a cylinder oscillating in still water The proposed method is applied to simulate the flow around a cylinder oscillating in still water. The flow and fluid forces were investigated experimentally and numerically by Dütsch et al. [31]. As shown in Fig. 3, a cylinder with the diameter D is placed at the center of a tank (15D × 8.75D) and is forced to oscillate. The position of the cylinder center (xc , y c ) is given by



xc (t ) = − A sin(2π f t )

(27)

y c (t ) = 0

where A denotes the amplitude of the cylinder oscillation. f is the oscillation frequency. A = 0.8D and f = 0.2 1/s are selected. The Reynolds number Re and the Keulegan-Carpenter number KC are defined as Re = (U max D ) /ν and KC = U max /( f D ), respectively, where U max is the maximum velocity of the cylinder and ν is the kinematic viscosity of water. The values of Re and KC are 100 and 5, respectively. These conditions are the same as the measurement and simulation conditions [31]. The computational domain is divided into 1500 × 875 square cells so that the cylinder diameter contains ten cells. Simulation with the time step t ∗ (= tU max / D ) = 8.68 × 10−4 is conducted from t ∗ (= tU max / D ) = 0 to 5, which corresponds to ten periods of oscillating motion. The time step is set so as to satisfy the stability conditions for the convection, diffusion and penalization terms [19]. The fluid forces acting on the cylinder are calculated at intervals of t ∗p = 64 t ∗ . The pressure is computed 90 times in one period, which is adequate for capturing the time variation of the fluid forces. The pressure gradient along the normal vector on the cylinder surface is calculated as follows:

∂ p

∂n

l→0

= −ρ ac · n 

(28)



where ac = ∂ 2 xc /∂ t 2 , ∂ 2 y c /∂ t 2 is the acceleration of the cylinder. Non-slip boundary condition is imposed on all boundaries; thus, the following conditions are imposed on the stream-function and pressure, respectively:

ψ =0

(29)

∂p =0 ∂n

(30)

The scalar potential φ is set to 0 over the entire domain because the inlet flow does not exist. In this simulation, the integral path along the cylinder surface is a closed curve; thus, the correction term p  vanishes. The distances from the surface I to points I1 and I2 are l1 = 3.75 x and l2 = 7.5 x, respectively.

Fig. 5. Pressure and vorticity distributions in four different phases at four different time points in fourth period of oscillation.

4. Numerical results The oscillating motion of the cylinder induces a flow around the cylinder. Fig. 4 shows the streamlines and the velocity distribution around the cylinder in the fourth period of the oscillation. The streamlines correspond to contour lines of the stream-function. The positions of the cylinder are xc / D = 0 and 0.4 in case of the motion toward the +x direction and xc / D = 0.4 in case of the motion toward the −x direction. The direction of movement of the cylinder is indicated by an arrow in the figure. As the cylinder moves, the vortex pair is shed symmetrically in the y direction from the sides of the cylinder. The rotating direction of the vortices varies as the direction of motion of the cylinder changes. Fig. 5 shows the distributions of pressure and vorticity around the cylinder at four different time points in the fourth period of the oscillation. When the cylinder moves toward the −x direction, the pressure acting on the forward direction of the cylinder increases, and a pair of vortices is formed in the wake of the cylinder. When the cylinder reaches the left end of the oscillation (x/ D = −0.8), the pressure on the cylinder surface increases due to the collision of the cylinder with the vortices formed in the wake of the cylinder. Thereafter, the cylinder changes its direction of motion from the −x to +x directions and reaches the right end of the oscillation. Similar phenomena were observed over a series

6

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

of cylinder oscillations. The pressure and vorticity distributions of the simulation results qualitatively agree well with the existing results obtained using the boundary-conforming formulation [31]. To quantitatively validate the simulation results, the simulated velocity profiles were compared with the values measured using a laser Doppler velocimeter [31]. The measurement locations of the four cross-sections are shown in Fig. 6. As the period number of the measured results is not mentioned in the literature [31], the simulated results in the fourth period are used for comparison. The comparisons are shown in Fig. 7. The velocity profiles show good agreement with the experimental results, demonstrating the validity of the simulation. The fluid forces F acting on the cylinder was calculated from the pressure and velocity fields. Fig. 8 shows the time histories of the total fluid force, pressure component and shear component of the force during 10 periods of the oscillation. Here, only the

Fig. 6. Points for measurement of velocity distribution.

component acting along the x-direction is shown, so F represents an in-line force. At the beginning of the simulation (0 ≤ t f ≤ 1) the in-line force is higher than the force simulated using the boundary-conforming formulation [31]. This difference can be ascribed to the pressure part of the results obtained using the proposed method. Because the starting time of the simulation is different from that in [31], the evolution of transient flow differs from that in [31]. The in-line force F in a period is shown in Fig. 9. This distribution is in the tenth period (9 ≤ t f ≤ 10) and it started from the time when F was almost 1. The results obtained using the boundary-conforming formulation [31] and the immersed boundary method [15] are also shown in Fig. 9. The present results agree well with the results of the existing simulations. However, F is slightly lower than the existing results when the cylinder changes its moving direction. When the direction of motion of the cylinder changes, the surface pressure increases owing to the collision of the cylinder with the vortices formed behind the cylinder, as shown in Fig. 5. In the present work, the two interpolation points, I1 and I2 , for pressure and velocity are considered far from the surface. It was confirmed that these differences decrease as the interpolation points are gradually moved closer to the surface. However, the in-line force F may be affected by the changes in the mask function χ expressed in Eq. (3) near the cylinder surface. To improve the results of F , it is necessary to narrow the range of change in χ , |l/ε |, and to move the interpolation points closer to the surface.

Fig. 7. Comparison of velocity distribution at four cross-sections over three different time points.

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

7

Fig. 8. Time histories of fluid dynamic force acting on cylinder and its components.

Declaration of competing interest All authors declare no competing interest regarding the research, authorship, and publication of this article. References

Fig. 9. History of fluid force over the tenth period.

5. Conclusions This study proposed a method for calculating the fluid forces acting on solid bodies in flows simulated by the VIC method with the volume penalization. The pressure field over the entire domain is computed from the pressure Poisson equation which is obtained by taking divergence of the Navier-Stokes equations. The pressure and the velocity on the body surface are estimated from a Lagrangian polynomial. This study also applied the proposed method to the simulation of the flow around an oscillating cylinder. The cylinder is forced to oscillate in still water in a rectangular tank. The Reynolds number and the KeuleganCarpenter number are 100 and 5, respectively. The simulation results were validated through successful comparison between the simulated velocity profiles and those obtained by means of laser Doppler velocimeter measurements. The fluid forces acting on the cylinder agreed well with the value reported in the existing literatures. These demonstrate the validity of the author’s method.

[1] L.H.L. Whitacker, J.T. Tomita, C. Bringhenti, An evaluation of the tip clearance effects on turbine efficiency for space propulsion applications considering liquid rocket engine using turbopumps, Aerosp. Sci. Technol. 70 (2017) 55–65. [2] K.N. Kiran, S. Anish, An investigation on the effect of pitchwise endwall design in a turbine cascade at different incidence angles, Aerosp. Sci. Technol. 71 (2017) 382–391. [3] H. Zhu, W. Hao, C. Li, Q. Ding, Simulation on flow control strategy of synthetic jet in an vertical axis wind turbine, Aerosp. Sci. Technol. 77 (2018) 439–448. [4] M. Righi, V. Pachidis, L. Konozsy, L. Pawsey, Three-dimensional through-flow modelling of axial flow compressor rotating stall and surge, Aerosp. Sci. Technol. 78 (2018) 271–279. [5] C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252–271. [6] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. [7] J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, J. Comput. Phys. 171 (2001) 132–150. [8] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries, J. Comput. Phys. 156 (1999) 209–240. [9] D. Kim, H. Choi, Immersed boundary method for flow around an arbitrarily moving body, J. Comput. Phys. 212 (2006) 662–680. [10] A. Gilmanov, F. Sotiropoulos, A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies, J. Comput. Phys. 207 (2005) 457–492. [11] Y. Yuki, S. Takeuchi, T. Kajishima, Efficient immersed boundary method for strong interaction problem of arbitrary shape object with the self-induced flow, J. Fluid Sci. Technol. 2 (2007) 1–11. [12] J. Choi, R.C. Oberoi, J.R. Edwards, J.A. Rosati, An immersed boundary method for complex incompressible flows, J. Comput. Phys. 224 (2007) 757–784. [13] E. Fadlun, A. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, J. Comput. Phys. 161 (2000) 35–60. [14] Y.-H. Tseng, J.H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, J. Comput. Phys. 192 (2003) 593–623. [15] J. Yang, E. Balaras, An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries, J. Comput. Phys. 215 (2006) 12–40.

8

T. Degawa et al. / Aerospace Science and Technology 94 (2019) 105360

[16] T. Ikeno, T. Kajishima, Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations, J. Comput. Phys. 226 (2007) 1485–1508. [17] P. Angot, C.-H. Bruneau, P. Fabrie, A penalization method to account obstacles in incompressible viscous flows, Numer. Math. 81 (1999) 497–520. [18] K. Schneider, M. Farge, Numerical simulation of the transient flow behaviour in tube bundles using a volume penalization method, J. Fluids Struct. 20 (2005) 555–556. [19] D. Kolomenskiy, K. Schneider, A Fourier spectral method for the Navier–Stokes equations with volume penalization for moving solid obstacles, J. Comput. Phys. 228 (2009) 5687–5709. [20] A. Sarthou, S. Vincent, J.P. Caltagirone, P. Angot, Eulerian–Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects, Int. J. Numer. Methods Fluids 56 (2008) 1093–1099. [21] I. Ramière, P. Angot, M. Belliard, A general fictitious domain method with immersed jumps and multilevel nested structured meshes, J. Comput. Phys. 225 (2007) 1347–1387. [22] G.-H. Cottet, P.D. Koumoutsakos, Vortex Methods: Theory and Practice, Cambridge University Press, New York, 2000. [23] D. Russell, Z.J. Wang, A Cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow, J. Comput. Phys. 191 (2003) 177–205. [24] M. Linnick, H.F. Fasel, A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains, J. Comput. Phys. 204 (2005) 157–192. [25] Z. Wang, J. Fan, K. Cen, Immersed boundary method for the simulation of 2D viscous flow based on vorticity-velocity formulations, J. Comput. Phys. 228 (2009) 1504–1520.

[26] G. Morgenthal, J.H. Walther, An immersed interface method for the vortex-incell algorithm, Comput. Struct. 85 (2007) 712–726. [27] J.T. Rasmussen, G.-H. Cottet, J.H. Walther, A multiresolution remeshed vortexin-cell algorithm using patches, J. Comput. Phys. 230 (2011) 6742–6755. [28] M. Gazzola, P. Chatelainm, W.M. van Rees, P.D. Koumoutsakos, Simulations of single and multiple swimmers with non-divergence free deforming geometries, J. Comput. Phys. 230 (2011) 7093–7114. [29] M. Gazzola, C. Mimeau, A.A. Tchieu, P.D. Koumoutsakos, Flow mediated interactions between two cylinders at finite Re numbers, Phys. Fluids 24 (2012) 043103. [30] C. Mimeau, I. Mortazavi, G-H. Cottet, A vortex penalization method for low and moderate Reynolds number flows, in: Proc. 7th Int. Conf. Vortex Flows and Vortex Models, Rostock, Germany, 2016. [31] H. Dütsch, F. Durst, S. Becker, H. Lienhart, Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers, J. Fluid Mech. 360 (1998) 249–271. [32] T. Uchiyama, H. Yoshii, H. Hamada, Direct numerical simulation of a turbulent channel flow by an improved vortex in cell method, Int. J. Numer. Methods Heat Fluid Flow 24 (2015) 103–123. [33] Y. Morinishi, T.S. Lund, O.V. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, J. Comput. Phys. 143 (1998) 90–124. [34] T. Kajishima, Conservation properties of finite difference method for convection, Trans. Japan. Soc. Mech. Eng., Ser. B 60 (1994) 2058–2063. [35] T. Nonomura, J. Onishi, A comparative study on evaluation methods of fluid forces on Cartesian grids, Math. Probl. Eng. (2017) 8314615.