Eulerian weakly compressible smoothed particle hydrodynamics (SPH) with the immersed boundary method for thin slender bodies

Journal of Fluids and Structures 84 (2019) 263–282

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Eulerian weakly compressible smoothed particle hydrodynamics (SPH) with the immersed boundary method for thin slender bodies ∗

A.M.A. Nasar , B.D. Rogers, A. Revell, P.K. Stansby, S.J. Lind School of Mechanical Aerospace and Civil Engineering, The University of Manchester, Manchester, M13 9PL, UK

highlights • • • • •

Lagrangian and Eulerian weakly compressible SPH are tested for FSI. Immersed boundary method used to model FSI with rigid thin body. Well-known treatments used to improve accuracy/stability of Lagrangian SPH. Results show Lagrangian SPH fails to accurately capture flow features. Eulerian SPH accurately captures flow highlighting advantages for FSI problems.

article

info

Article history: Received 24 November 2017 Received in revised form 7 June 2018 Accepted 8 November 2018 Available online xxxx

a b s t r a c t Fluid–structure interaction (FSI) problems involving thin slender bodies are a difficult modelling problem particularly for highly flexible bodies. An accurate, numerically stable approach is presented for fixed bodies as a first stage where weakly compressible smoothed particle hydrodynamics in Eulerian form (EWCSPH) is coupled with the immersed boundary method (IBM) and applied to the problems of an accelerating square box, as a nonslender case, and an impulsively started flat plate as a slender body case. The results for the box case show that vortical flow structures are virtually identical to those predicted by industry standard finite-volume solvers for lower Reynolds numbers (up to 150). For a higher Reynolds number of 450 there is some evidence that EWCSPH provides more accurate vortical structures for equivalent mesh/particle resolutions. Computational resource requirements are similar. The results for the plate case demonstrate that to simulate strong dynamic shear layers and vortex shedding the Eulerian form of SPH offers significantly greater stability and accuracy in comparison to the conventional Lagrangian form. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Weakly compressible smoothed particle hydrodynamics (WCSPH) is well known for its ability to model non-linear fluid flows such as multi-phase flows (Mokos et al., 2015), sediment transport (Fourtakas et al., 2013) and a wide range of freesurface flow applications (Gómez-Gesteira and Dalrymple, 2004, Longshaw and Rogers, 2015 and Antoci et al., 2007). WCSPH has also been applied to a number of FSI applications ranging from the violent interaction of waves with rigid floating bodies (Canelas et al., 2013) to the analysis of complex fluid behaviour within flexible plant cells (Karunasena et al., 2014 and Van Liedekerke et al., 2010). For many industrial and biomedical FSI problems, vortex shedding plays a crucial role as ∗ Corresponding author. E-mail address: [email protected] (A.M.A. Nasar). https://doi.org/10.1016/j.jfluidstructs.2018.11.005 0889-9746/© 2018 Elsevier Ltd. All rights reserved.

264

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

it induces structural vibration which feeds back and alters the flow field. This is especially evident in problems involving highly flexible structures such as flow through flexible valves either in industrial applications or biological applications such as the human cardiovascular system. These are particularly challenging modelling problems yet to be resolved. Quinlan et al. (2006), Fatehi et al. (2008) and Basa et al. (2009), have shown that the accuracy and stability of the SPH approximation is strongly degraded by the irregularity of particle distributions with non-uniform inter-particle spacing. This has been controlled by the particle regularisation techniques of shifting (see Xu, 2010, Lind et al., 2012 and Khayyer et al., 2017) and several diffusion methods in WCSPH: Shepard filter (Panizzo, 2004), XSPH (Monaghan, 1989) and δ -SPH (Marrone et al., 2011). δ -SPH in combination with shifting, as used by Luis Cercos-Pita (0000), is investigated here. These methods generally enable convergence between first and second order. Hieber and Koumoutsakos (2008) investigated flow past bluff bodies in 2-D and 3-D and also flow past slender swimmers. They used a variant of WCSPH (r-SPH) which employs a re-meshing technique to maintain uniform particle distributions. This improves the accuracy and stability of the method in comparison to fully Lagrangian WCSPH (LWCSPH). They also employ the immersed boundary method (IBM), initially developed by Peskin (1972), to model the boundary effects and their combined r-SPH IBM method captures the flow structures in the wake of slender swimmers. However, the flow structures appear diffused in comparison to the reference results used for comparison, possibly due to the additional smoothing required for the re-meshing procedure. This is confirmed by the work of Cherfils (2011) where the r-SPH technique is compared with well-established flow solvers and a moving least square immersed boundary method for WCSPH (IBMLS-SPH) showing that the vorticity field is indeed diffused over a larger area when using r-SPH in comparison to reference results (Figure 6.22 in Cherfils (2011)). Eulerian Incompressible SPH (ISPH) is a recent proposition made by Noutcheuwa and Owens (2012) and by Lind and Stansby (2016). Without a free surface it has been shown that high-order convergence is possible, through high-order kernels, and that special treatment does not seem necessary for numerical stability. This has been demonstrated for simple geometries to date but it is expected that extension to complex geometries will follow with improvements in boundary conditions to high-order (Nasar et al., 2018). It is also anticipated that particle position generation with adaptivity for complex geometries will be readily automated (Fourtakas et al., 2018) and this is a particular blockage for mesh/cell based methods. The downside is the number of particle connections within a kernel, typically 20 to 30 in 2D and 200 to 300 in 3D. While there are developments here also, these numbers are much greater than cell connections in mesh-based methods. However, particle methods are ideally suited to parallel processing particularly with modern, inexpensive and energy efficient processing available through GPUs (Crespo et al., 2011; Chow et al., 2018). The main aim of this work is to investigate and establish the capabilities of WCSPH coupled to the immersed boundary method for slender and bluff bodies immersed in a flow. Both Eulerian and Lagrangian WCSPH are evaluated. The motion of a square box through a quiescent fluid within a sealed tank is investigated as a standard test case (SPHERIC Benchmark 6) for comparison of EWCSPH with industry-standard CFD. An impulsively started rigid plate problem is also investigated as it involves a most severe test with shear layer mixing from a thin surface, this test is used to validate EWCSPH and for comparison with Lagrangian forms of WCSPH. In practice this type of flow occurs in applications involving very thin structures such as cilia or thin flexible valves which are subjected to fluid loading in the human body. This paper is structured as follows: in the following Section 2 the formulation for LWCSPH and EWCSPH is presented incorporating techniques used to limit the effects of irregular particle distributions on the accuracy and stability of LWCSPH. In Section 3, the formulation for the IBM is presented in the form suitable for WCSPH. Section 3 also covers the formulation for the MLS kernel correction which was used to improve the accuracy of the IBM kernel interpolations for cases involving disordered particle distributions. Results are presented in Section 4 and in Section 5 while some conclusions are made in Section 6. 2. Lagrangian and Eulerian weakly compressible SPH 2.1. Formulation The Navier–Stokes equations can be written in Lagrangian form as du dt dρ

1

= − ∇ P + ∇ · (ν∇ u) , ρ

(1)

= −ρ∇ · u, (2) dt where u is the fluid velocity vector, ρ is fluid density, t is time, P is fluid pressure and ν is kinematic fluid viscosity. The discrete Lagrangian WCSPH Navier–Stokes equations can be written as Monaghan (2005) dui dt

=− SPH

dρi dt

j=1

= SPH

N ∑

N ∑ j=1

( mj

Pi

ρi2

+

Pj

ρj2

) ∇i Wij + Πij ,

mj (ui − uj ) · ∇i Wij + δi ,

(3)

(4)

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

265

where m is the mass of an SPH particle, N is the total number of fluid particles and subscripts i and j reference the discrete properties evaluated at the positions of the fluid particle of interest and its neighbours, respectively. Πij denotes the viscous terms in the discrete SPH formulation and ∇i Wij is the gradient of the smoothing kernel where Wij = W (xj −xi , h) for the interaction between fluid particles i and j where h is the smoothing length and Wij has compact support. δi is an artificial density diffusion term used to limit spurious density fluctuations as discussed in Section 2.2. A stiff equation of state (EOS) is used for closure of the Navier–Stokes equations. The EOS relates fluid pressure to changes in density via P =

co2 ρo

((

γ

ρ ρo

)γ

) −1 ,

(5)

where γ = 7 is the fluid polytropic index used herein such that small variations in density lead to large variations in pressure, greatly limiting fluid compressibility, and co is the reference numerical speed of sound. co is typically set to co = 10Umax where Umax is the maximum expected flow velocity such that the overall flow Mach number Ma ≤ 0.1. Given Ma2 ∼ (ρ−ρo )/ρo , for Ma ≤ 0.1 the expected density variation is less than 1%. Fluid compressibility was found to detrimentally affect vortical flow structures due to pressure wave propagation as discussed in Sections 4 and 5. Hence, the values used in this paper range from co = 50Umax to co = 100Umax so that the estimated fluid compressibility would be minimal (less than 0.04%) without rendering the computational cost unfeasible. The governing equations for WCSPH can be simply converted into Eulerian form. This requires re-writing the momentum equation such that N

∑ ∂ ui =− mj ∂ t SPH

(

Pi

ρi2

j=1

+

)

Pj

∇i Wij + Πij − ((u · ∇ u)SPH )i ,

ρj2

(6)

where (u·∇ u)SPH is the SPH discrete form of the advective flux term arising from re-writing the total derivative du/dt in Eq. (1) as the sum of partial derivatives (∂ u/∂ t + u·∂ u/∂ x). The convective term in Eq. (6) can be discretised using the smoothing kernel’s derivative as

((u · ∇ u)SPH )i = ui ·

N ∑ [

(uj − ui ) ⊗ ∇i Wij Vj ,

]

(7)

j=1

where ⊗ denotes the tensor product. The laminar viscosity model, introduced for SPH by Morris et al. (1997) is used for the investigations herein and is written as

( Πij = ∇ ·

µo ∇ ui ρ

) = SPH

N ∑

( mj

2µo rij · ∇i Wij

ρi ρj |r|2ij

j=1

) ,

uij

(8)

where µo is the dynamic fluid viscosity and subscript o denotes the reference value of a fluid property. The Eulerian form of the continuity equation can be written as

∂ρ = −ρ∇ · u − u · ∇ρ. ∂t

(9)

Since weakly compressible SPH is used herein for flows without rapid variations in density ∇ρ ≈ 0. The term u·∇ρ can thus be omitted to give N

∑ ∂ρi ≈ mj (ui − uj ) · ∇i Wij . ∂ t SPH

(10)

j=1

The simulations outlined in this paper were performed with and without incorporating the term u·∇ρ into the continuity equation. Pressure drag coefficient measurements for the test in Section 4 were practically identical, within 1×10−3 % of each other, confirming that the term may be neglected. Time integration of the system of equations is performed via a second-order predictor corrector scheme (Batchelor, 2000) which allows for a variable time step size. The size of the maximum allowable time step is dependent on the maximum fluid velocity, the maximum force acting on any fluid particle in the flow field and the fluid viscosity, in the work herein this is calculated using

[

(

√

)]

∆tmax = CFL min (∆x/co ), min ( h/fi ), (0.5h /ν ) i∈fluid

2

,

(11)

where fi is the force acting on fluid particle i at the current time step and CFL denotes the Courant–Friedrichs–Lewy coefficient, typically set to 0.1 for LWCSPH.

266

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

2.2. δ -SPH The conventional WCSPH approach can suffer from spuriously fluctuating pressure and velocity fields. In the δ -SPH method (Molteni and Colagrossi, 2009) an artificial diffusion parameter is used to diffuse the spurious fluctuations and minimise their effect on the flow field. As suggested by Molteni and Colagrossi (2009), a density diffusion term δi is added to the right hand side of the continuity equation (Eq. (4)) where N ∑ mj

ψij · ∇ Wij , ρj ( ) ρj xi − xj ψij = 2 −1 ⏐ , ⏐ ρi ⏐xi − xj ⏐2 + 0.1h2 δi = ad hc0

(12)

j=1

(13)

and ad is the diffusion parameter typically set to 0.1. 2.3. Diffusion-based particle shifting Particle shifting is an increasingly standard technique for Lagrangian SPH which is used to prevent particle clustering and/or disordered particle arrangements which can cause numerical instability. As the flow field evolves, SPH particles are displaced slightly into a more uniformly distributed configuration. A number of techniques such as those of Xu et al. (2009) and Shadloo et al. (2012) were tested but the method which provided the most stable results for the impulsively started plate test was found to be the diffusion based method of Lind et al. (2012). With this method a concentration gradient is used to evaluate a directional measure of particle disorder through

∇ Ci =

N ∑

mj ∇ Wij /ρj .

(14)

j=1

The concentration gradient is then used to diffuse the particles into a more uniform distribution following Fick’s law

∆xs = −Di ∆t ∇ Ci ,

(15)

where Di is a numerical diffusion coefficient and ∆xs is a shifting displacement vector. From a stability analysis highlighted in Lind et al. (2012); Skillen et al. (2013) the maximum theoretical limit for the diffusion coefficient required to avoid instability is Di = 0.5h2 /∆t .

(16)

The numerical diffusion coefficient can be re-written as Di = λh2 /∆t ,

(17)

where λ is a dimensionless parameter which defines the maximum shifting magnitude for a given concentration gradient per time step. From comparing Eq. (16) and Eq. (17) λ < 0.5 is required to avoid instability. 2.4. Kernel gradient correction To improve the accuracy of the SPH approximation of derivatives for cases with disordered particle distributions a correction for the kernel gradient was proposed by Bonet and Lok (1999). This is evaluated as,

∇ c Wij = Li ∇ Wij ,

(18)

c

where ∇ Wij is the corrected kernel gradient and

⎛ Li = ⎝

N ∑

⎞−1 ∇ Wij ⊗ rij Vj ⎠

.

(19)

j=1

As can be inferred from the formulae above this procedure adds a significant computational cost to the solver algorithm as it involves an additional summation which must be performed over the entire fluid field prior to the summations in Eq. (3) and Eq. (4) and the inversion of N matrices where N here is the total number of fluid particles. Also, the corrected kernel gradient is non-symmetric hence SPH schemes using kernel gradient corrections are non-conservative.

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

267

Fig. 1. (a) 2-D definition sketch highlighting a planar boundary surface Γ discretised into a set of markers (star shapes) immersed in a discrete fluid domain (blue circles denote fluid particles). The sketch illustrates the radius of influence (2h) of an immersed boundary marker b (highlighted in red). (b) A membrane-like surface is fully immersed in SPH fluid particles. The arrows represent the velocity vectors of the fluid particles and the nodes used to discretise the structure, which are also used as boundary markers.

3. The Immersed Boundary Method The Immersed Boundary Method (IBM) has proven to be successful for modelling FSI problems involving complex geometries and/or high levels of boundary deformation. Initially developed for mesh-based methods, it has been used as an alternative to re-meshing which is a complex and computationally expensive process and can lead to high levels of mesh skewness potentially corrupting the computations. Since its inception many modifications have been proposed to the original IBM of Peskin (1972). Mittal and Iaccarino (2005) have listed some of the most prevalent variants highlighting their advantages and limitations. Uhlmann’s diffuse interface direct forcing IBM formulation (Uhlmann, 2004) was developed to model rigid particle laden flows but it is also useful for moving flexible boundaries (Favier et al., 2014). The formulation is robust since it does not involve the use of empirical stiffness parameters which are used in some continuous forcing diffuse interface IBMs which can lead to instability or the requirement for very low CFL numbers when modelling rigid boundaries (Mittal and Iaccarino, 2005). Hence Uhlmann’s approach is used as the formulation for the WCSPH-IBM proposed herein. 3.1. Formulation The principle of the IBM is that two superimposed stencils are used to discretise the fluid and the structure separately as shown in Fig. 1 while forcing and spreading procedures are used to ensure continuity of stresses at the interface. Forces are evaluated in the region where the two stencils overlap (the immersed boundary) in a manner which approximates the no-slip condition such that the resulting forces are applied to both the fluid and the structural solvers as external forces in a conservative formulation. Note: The no-slip and no-penetration conditions are only approximately enforced since numerical experiments have shown that the volume of fluid enclosed within an elastic immersed boundary decreases with simulation time (Griffith, 2012). Smaller time steps and/or the use of semi-implicit time integration can reduce fluid penetration across the immersed boundary (Griffith, 2012). Lagrangian boundary markers Fig. 2 shows a conceptual flow diagram of the main steps that will now be explained. To account for the immersed boundary in the fluid solver the position of the structure/boundary is tracked via a set of Lagrangian boundary markers shown as star shapes in Fig. 1(a) which also shows an example radius of support centred on a boundary marker. The markers are used to interpolate the force required (commonly referred to as ‘force density’ in IBM formulations) to approximately enforce the no-slip condition at the position of the markers. This also enables the hydrodynamic load exerted on the structure to be computed as being equal, and opposite in direction, to the load required to approximately enforce the no-slip condition at the position of the Lagrangian markers.

268

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

For fully coupled FSI problems, where the structure is free to suffer deflections and/or full body displacements due to fluid loading, the nodes used to discretise the structure and solve the structural mechanics governing equations are typically used as the Lagrangian markers. Fig. 1(b) shows an illustrative snapshot from a 2-D simulation performed using the SPH-V Model depicting the deformation of a flexible membrane immersed within a fluid (Nasar et al., 2016). IBM for fluid and solid The full description of the IBM method for SPH now follows. Throughout this section vector variables written in lower case font such as position x, velocity u and the right-hand side of the momentum equation rhs refer to properties belonging to the fluid. Upper case symbols X and U are used to signify the position and velocity of the Lagrangian markers located on the structure. Finally, upper case RHS denotes acceleration acting on the Lagrangian markers which are coincident with and rigidly attached to the structural nodes used to discretise the structural governing equations. Recalling the Lagrangian Navier–Stokes momentum equation (Eq. (1)) in the absence of external forces the terms on the right-hand side of Eq. (1) can be grouped into 1 rhsL = − ∇ P + ∇ · (ν∇ u)

ρ

(20)

In similar fashion, rearranging the Eulerian formulation gives 1 rhsE = − ∇ P + ∇ · (ν∇ u) − u · ∇ u

ρ

(21)

The following three steps are then followed for the implementation of the IBM formulation. The first two steps are the velocity interpolation and forcing step and then the spreading step. Since these form a distinct part of the IBM formulation they are discussed separately. Step 3 highlights how the IBM forces are incorporated within the time integration process. Hence the description of the time integration algorithm is presented within the description of Step 3. Step 1: Velocity interpolation and forcing (a) Structure The equations of motion governing the behaviour of the structure are integrated one step forward in time to evaluate the velocity the structure would possess if it were not immersed in the fluid. For a Lagrangian marker, denoted by subscript b, located on a structural discretisation point, the desired velocity, denoted by Ud , at the end of the time step would be evaluated as n+1/2

Udb = Unb + RHSb

∆t ,

(22)

where RHSb is the total acceleration acting on each individual boundary marker b. This can be due to external forces, such as gravity for example, or due to internal structural forces for cases involving flexible structures or due to both internal and external forces. For the tests in this paper the immersed body motion is prescribed as a function of time with RHSb = (0, 0) for the impulsively started plate test and RHSb = (ax , 0) for the accelerated box test where ax is evaluated via Eq. (29) . (b) Fluid Next, using the rhs values from Eq. (20) for Lagrangian SPH or from Eq. (21) for Eulerian SPH, the fluid governing equations are integrated to time n + 1/2 in the predictor step of the time integration scheme which gives the predictor values p

uj = unj + rhsnj ∆t /2,

(23)

p xj

/2,

(24)

∆t /2,

(25)

=

xnj

+

ρjp = ρjn +

p uj ∆t n

dρ dt

j

where superscript p denotes predictor and ρ p j is the predicted density of particle j at time step n + 1/2. Note: the E and L subscripts are dropped for this part of the formulation as the procedure is the same for both Lagrangian and Eulerian WCSPH. Using the predictor values for the fluid variables the forces arising from the Navier–Stokes equations are evaluated without the inclusion of any external force due to the immersed boundary. This gives the acceleration rhsn+1/2 acting on the fluid which is used to evaluate a preliminary velocity field at time n+1 for all points j in the fluid n+1/2

u∗j = unj + rhsj

∆t .

(26) ∗

The intermediate velocity u is the velocity that the fluid would possess at time n+1 in the case that the immersed boundary was not present in the flow field. The time integration scheme is discussed later in Step 3. With knowledge of u∗ for the fluid particles a smoothing kernel interpolation can be used to evaluate the preliminary velocity at the position of Lagrangian marker b as U∗b =

Nf ∑ j

u∗j Wbj Vj ,

(27)

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

269

where Nf is the total number of fluid particles, subscript b denotes a boundary marker and subscript j denotes neighbouring fluid particles. Fig. 1(a) illustrates the stencil used for the velocity and force interpolations. Based on the desired velocity Ub d and the preliminary velocity Ub ∗ , the force density required to maintain Ub d at the position of a marker b (approximately enforcing no-slip condition) is evaluated as n+1/2

Fb

= ρo

Udb − U∗b

∆t

,

(28)

where ρo is the reference fluid density. kernel interpolations can be used to evaluate the fluid density at the ∑ nSmoothing n+1/2 +1/2 position of the marker as ρb = ρj Wbj Vj but since the fluid is only weakly compressible the deviation from the reference density is small, less than 1% for cases with the recommended co = 10Umax where Umax is the maximum fluid velocity. For most cases in this paper co > 10Umax so the deviation from ρo is even smaller than 1%. Since Fb represents the force per unit volume required to enforce the no-slip condition for the fluid at the position of the Lagrangian marker, the force per unit volume acting on a Lagrangian marker due to its interaction with the surrounding fluid is −Fb and the acceleration vector arising from this force is evaluated as n+1/2

ab

n+1/2

= −Fb

/ρs ,

(29)

where ρs is the structural material density. This acceleration is added as a collection of point loads acting on each individual point b in the discrete structure for fully coupled FSI cases as presented in Nasar et al. (2016). Step 2: Force Spreading The force density in Eq. (28) must be spread back onto the fluid in a manner which enforces the no-slip condition, to a good level of accuracy. This is achieved via the spreading procedure where the force evaluated at the position of a marker b is spread back onto neighbouring fluid points n+1/2

fj

=

Nb ∑

n+1/2

Fb

Wjb Vb ,

(30)

b

where Nb is the total number of boundary markers. Note that since the summation in the right hand side of Eq. (30) is an approximation of an integral (Uhlmann, 2005) conservation of momentum can only be guaranteed if Nb ∑

Wjb Vb = 1.

(31)

b

The MLS-IBM variant of the IBM, discussed in Section 3.2 and used for tests in Section 5, adopts moving least square corrections (Vanella and Balaras, 2009) which ensures Eq. (31) is satisfied by correcting the shape of the kernel function to account for imperfect kernel support and ensure first-order consistency. The velocity interpolation and spreading steps (steps 1 and 2) form the basis of the interpolation procedures employed within the IBM. The remainder of the IBM formulation and the implementation algorithm are heavily intertwined with the time integration procedure because they involve the use of the IBM forces to update the fluid properties at various points in time. Hence, details of the remaining steps are incorporated into the following description of the time integration procedure which is labelled step 3. Step 3: Time integration of the governing equations to time step n + 1 Following the approach highlighted in steps 1 and 2, the forces required to approximate the no-slip/no-penetration conditions at the position of the boundary can be evaluated such that the solution may be marched through time in a manner where the no-slip and no-penetration conditions are satisfied to a good level of accuracy. As the WCSPH formulation is used the WCSPH-IBM system is integrated through time via an explicit scheme. The predictor–corrector time integration scheme (Batchelor, 2000) is suited for WCSPH as it provides a good compromise between stability, accuracy and computational expense and importantly allows for the use of a variable time step size. The description of the predictor– corrector scheme as used for the studies herein is as follows: Using the predicted state of the fluid at time step n + 1/2 (from Eq. (23) to Eq. (25)) the continuity equation (Eq. (4)) and the momentum equation (Eq. (3) for LWCSPH and Eq. (6) for EWCSPH) are evaluated to give rhsn+1/2 . With knowledge of rhsn+1/2 , u∗ is evaluated for all the fluid particles via Eq. (26) and interpolated onto the Lagrangian markers via Eq. (27) to give U∗ which is used to evaluate the Lagrangian force density Fn+1/2 as in Eq. (28). The forces are spread back onto the fluid using Eq. (30) to evaluate fn+1/2 which is then added to the momentum equation for LWCSPH

ρ

du n+1/2

n+1/2

+ fn+1/2 ,

(32)

∂ u n+1/2 n+1/2 = ρ rhsE + fn+1/2 . ∂ t SPHj

(33)

dt SPHj

= ρ rhsL

and for EWCSPH

ρ

270

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

Fig. 2. (a) Flow chart illustrating the time integration scheme for SPH-IBM solver. (b) Flow chart for IBM forcing and spreading steps.

The fluid corrector step is then performed to give ucj = unj + ucj

=

unj

du n+1/2 dt SPHj

n+1/2

∆t /2 and xcj = xnj + uj

∂ u n+1/2 + ∆t /2 for Eulerian SPH . ∂ t SPHj

∆t /2 for Lagrangian SPH , (34)

Finally, unj +1 = 2ucj − unj ;

xnj +1 = 2xcj − xnj .

(35)

n

Note that for EWCSPH xn + 1 = x since the fluid particle positions are fixed in space. The three steps discussed above are summarised in Fig. 2 which shows a flow chart description of the IBM methodology incorporated within the SPH time integration scheme. For the case with prescribed structural motion at constant velocity investigated herein the velocity of the boundary is not influenced by the fluid and hence Ub d = Ub n = Ub n+1 . A more detailed description of the IBM method discussed above can be found in Nasar (2016) where the methodology description is extended to cover incorporation of the time integration scheme for flexible structures into the IBM framework. 3.2. MLS-IBM For some of the Lagrangian tests a moving least squares correction is used to correct the shape of the smoothing function to improve the accuracy of the IBM interpolations near the boundary in the case of disordered particle arrangements and recover first-order consistency. Following the approach of Cherfils et al. (2010) the IBM smoothing kernel function is corrected in 2-D via WbjMLS = βb · 1, xbj , ybj

(

(

))

Wbj ,

(36)

where W MLS replaces W for the forcing and spreading interpolations (Eq. (27) and Eq. (30)), xbj = (xb − xj ) and ybj = (yb − yj ) are the distances between fluid particle j and boundary marker b in the x and y directions, βb is a corrective vector evaluated

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

271

as

( βb = Ab

−1

1 0 0

) ,

(37)

and Ab is evaluated based on the distribution of particles near the point of interpolation b such that

⎡ 1 N ∑ ⎢ Ab = ⎣ xbj j=1

ybj

xbj

ybj

⎤

xbj ybj ⎦ Wbj Vj .

x2bj

⎥

xbj ybj

(38)

y2bj

−1

Ab in Eq. (37) is constructed using the determinant and cofactors of Ab . 4. Accelerating square box in a confined domain In this section the ability of EWCSPH to model incompressible flow in a sealed domain is validated and assessed via the SPHERIC Benchmark (BM) 6 (Colagrossi, 0000) in Section 4.1. A comparison of EWCSPH-IBM computational requirements with an industry standard CFD solver is presented in Section 4.2. 4.1. Test setup and results The SPHERIC BM 6 test involves accelerated motion of a 2-D square box, initially at rest, through a 2-D rectangular sealed tank filled with a quiescent fluid. This benchmark test was developed to validate the capability of different WCSPH formulations to accurately portray incompressible pressure fields. Given the complexity of the box geometry which involves right angle corners it is also useful for validating the accuracy of wall boundary conditions. This case has been previously investigated using remeshed-SPH (Cherfils et al., 2009) and a number of LWCSPH variants (the works of Cherfils et al. (2010) and Vacondio et al. (2013) are a few examples). Herein, the case is investigated using the EWCSPH-IBM and the ‘‘immersed solid’’ module (CFX-IBM) available for the ANSYS CFX incompressible fluids solver which uses a control volume finite element model. Results acquired using a level set-Navier Stokes Finite Difference solver (LS-NS) (Colagrossi, 0000) are also used for comparison. For the EWCSPH-IBM simulation the particle spacing was set to ∆x = L/100 where here L is the side length of the square box but for these tests the ratio h/∆x was lowered to 1.4 to reduce the computational requirements while still providing a good level of accuracy. The same resolution was used to discretise the fluid for CFX-IBM while the results for the LSNS solver (Colagrossi, 0000) were acquired using ∆x = L/60. The box was therefore discretised into 100×100 Lagrangian markers for both EWCSPH-IBM and CFX-IBM. Four tests were performed using EWCSPH-IBM and CFX-IBM where the fluid viscosity was varied (µ = 1/5, 1/10, 1/15, and 1/45 Pa s) while the reference fluid density was kept constant at ρo = 10 kg/m3 for all 4 tests. The external acceleration imposed on the box for the LS-NS tests (Colagrossi, 0000) varied with time according to ax (t) =

e−25(t −0.5275)

√

π /25

2

.

(39)

Integration of the above with respect to time and using suitable integration boundary conditions gives Ubox (t) = 0.5 (1 + erf (5t − 211/80)) ,

(40)

where Ubox (t) is the box velocity. Note, at t = 1 s Eq. (39) gives a(t) ≈ 0 m/s2 and hence the final box velocity of Ubox max = 1 m/s and is reached at t = 1 s, the box continues to move through the fluid with a velocity of 1 m/s after that point in time. Eq. (40) is used to prescribe the box velocity for all the EWCSPH-IBM and CFX tests therefore the Reynolds number is matched to the LS-NS tests with Re = 50, Re = 100 and Re = 150. Colagrossi (0000) did not perform the test for Re = 450 which corresponds to the lowest value used for viscosity for the tests herein. For EWCSPH-IBM the CFL number was set to CFL = 0.5 and the numerical speed of sound was set to co = 50 m/s, the time step size was fixed for these computations based on the particle spacing and the numerical speed of sound ∆t = 0.5∆x/co = 1×10−4 s. WCSPH simulations of this problem (as presented herein and also others in literature (Cherfils et al., 2010; Luis Cercos-Pita, 0000) predict the creation of pressure waves due to the box acceleration. These pressure waves travel faster than the box and reflect off the tank walls multiple times before the box has traversed the length of the tank. If left unchecked, pressure wave interference leads to large deviations from the pressure field expected for incompressible flow which eventually corrupt the numerical solution when the pressure fluctuations become unresolvable by the discrete domain (when the fluctuation wavelength is smaller than ∆x). This has been alluded to previously by Cherfils et al. (2009, 2010) who use a modified continuity equation to help dissipate the effects of the interfering pressure waves. Herein, δ -SPH is used for this purpose with δ = 0.25 for the Re = 50 test and δ = 0.5 for the higher Re tests. For the CFX tests the time step was set to ∆t = 1×10−3 s and the normalised residual target for the momentum and continuity equations was set to 1×10−6 . This lead to convergence of the implicit time integration scheme, to the residual

272

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

Fig. 3. Vorticity contours for SPHERICIC BM 6 simulations with various Re. The top half of each image shows EWCSPH-IBM results while the bottom half shows CFX-IBM results. For (a) and (b) Re = 50, (c) and (d) show results for Re = 100, for (e) and (f) Re = 150 and (g) and (h) show results for Re = 450.

target, within 2 to 9 iterations per time step, the number of iterations required for convergence varied for each time step and was also dependent on Re. A Cartesian mesh with constant control volume size was used in order to minimise any mesh-related errors arising from variations in element size and mesh skewness and also to allow for a direct comparison with the EWCSPH-IBM results. The same mesh was used for the tests with Re = 50, Re = 150 and Re = 450. Vorticity contours are plotted in Fig. 3 highlighting the development of the flow field at different points in time for both the EWCSPH-IBM and CFX tests. The contours show that the flow patterns from Re = 50 and Re = 150 are virtually identical for the two solver results especially for lower Re. The lowest value vorticity contours in the wake of the box at t = 5 s are generally slightly longer for the CFX-IBM results than for the EWCSPH-IBM results. However, as shown in Fig. 4, the solution appears to be underresolved for the CFX-IBM results with Re = 450 where a secondary vortex sheet, counter-rotating with respect to the shear layer rotation is predicted to be attached to the frontal shear layer which is formed on the box’s leading face. A convergence study presented in Section 5.1 shows that this is an unphysical numerical artefact arising from underresolved flow. The unphysical vortex sheets predicted by CFX-IBM may also be due to the rotating flow near the box leading

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

273

Fig. 4. Magnified vorticity contours for SPHERIC BM 6 with Re = 450. The top half shows EWCSPH-IBM results while CFX-IBM results are presented in the lower half.

Fig. 5. Pressure drag coefficient results for SPHERIC BM 6 with Re = 50 (a), Re = 100 (b), Re = 150 (c) and Re = 450 (d).

face being misaligned with the computational mesh used for the CFX-IBM tests which is Cartesian, that is aligned with the x and y axes. Smoothing kernel interpolations are utilised for EWCSPH-IBM and these are much less sensitive to grid/flow misalignment than FE and FV interpolations due to the spherical nature of the stencil used for interpolation. Also, the number of neighbours used for each EWCSPH-IBM fluid to fluid interpolation is 20 in comparison to only 4 neighbours used for the CFX-IBM interpolations. Similar numerical artefacts can be noted for the Re = 150 CFX-IBM results near the edges of the main frontal shear layer but to a much smaller extent, see Fig. 3 (e) and (f). Pressure fields from a weakly compressible solver inevitably show fluctuations due to compressible pressure waves not present in incompressible solvers. However, pressure drag variation is shown to be in reasonable agreement in Fig. 5 with some fluctuation in EWCSPH due to pressure waves. This effect will be less with open boundaries or a free surface as the pressure waves disperse from the domain. Fig. 5 shows the time variation of pressure drag coefficient for the EWCSPH-IBM, ANSYS CFX and LS-NS (Colagrossi, 0000), results for the latter at Re = 450 were not available and are thus not shown. The drag coefficient was computed as p

p

Cd = −

2Fd

max ρ Ubox L

,

(41)

274

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

where, p

∫

Fd =

Γ

Pnx dy ≈

∑

Pi nx ∆y,

(42)

i

where Γ is the box surface, nx is the wall normal vector component in the x-direction, ∆y is the particle spacing in the y-direction and Pi =

∑

/

j∈Ω + Pj Wij Vj

∑

j∈Ω + Wij Vj

.

(43)

Note that Ω + denotes fluid particles located outside the box at any point in time. Fluid particles lying inside the box are excluded from the pressure interpolation in Eq. (43) to match the sampling procedure used by ANSYS CFD-Post to calculate the forces acting on the box leading and trailing faces which only uses element nodes lying outside the immersed solid to interpolate fluid variables. The results show that EWCSPH-IBM slightly over predicts the peak value of Cd p in comparison to both the LS-NS and CFXIBM solvers. Also, the mean quasi-steady state Cd p value is roughly constant for both CFX-IBM and LS-NS-FD with a relatively small-amplitude high-frequency oscillation about the mean, here quasi-steady state refers to the period of time after the box has reached its final velocity. The EWCSPH-IBM solver predicts this high-frequency low-amplitude oscillation but also an additional lower-frequency higher-amplitude oscillation about the mean. This discrepancy and the slightly over-predicted mean pressure values are directly a result of the pressure wave reflections discussed previously in this section which appear to have less of an effect with reduced viscosity. In comparison to other results available in the literature, gathered using LWCSPH variants, it may be noted that previous LWCSPH models predict additional noise for Cd p (for example see Fig. 9 in Cherfils et al. (2010)). This additional noise is due to disordered particle distributions which are difficult to eliminate when using LWCSPH, as highlighted further by the results in Section 5.2 for the impulsively started plate test which shows that these spurious oscillations heavily corrupt the vorticity field. The EWCPSH-IBM captures smooth vortical flow structures for all tests presented in this paper which appear unaffected by the pressure waves arising from compressibility. Still, the results of the SPHERIC BM 6 tests suggest some limitations for the use of any weakly-coupled WCSPH model for incompressible FSI in sealed domains. A stronger coupling for the governing equations can be achieved by employing Eulerian ISPH so it is expected that an Eulerian ISPH-IBM solution of SPHERIC BM 6 would produce pressure fields in better agreement with the CFX-IBM and LS-NS results. Since investigation of Eulerian ISPH techniques is beyond the scope of this paper it will be the focus of future work.

4.2. Computational requirement analysis The average computational time required to complete 1 time step was evaluated for both the EWCSPH-IBM and the CFX-IBM as the time required to complete 20 time integration steps divided by 20. For the CFX-IBM stability is not affected by the size of the time step. However, for a given problem, larger time steps require more iterations to achieve convergence to the desired normalised tolerance value for the momentum and continuity equations, which for these tests was set to 1 × 10−6 to improve the overall accuracy for CFX-IBM. As such, the following figure shows the computational time required to complete ∆tstep = 1 × 10−3 s of simulation time for both solvers (which is 1 time step for CFX-IBM and 10 time steps for EWCSPH-IBM) plotted against the number of processors used for the computations. R All the results for the computational cost comparison are for the Re = 150 tests which were performed using an Intel⃝ R Xeon⃝ E5-2620 v4 CPU with a clock speed of 2.1 GHz with 64 GB of RAM available. Hyperthreading is known to adversely affect the speed of parallel ANSYS CFX computations (ANSYS user manual) and hence it was disabled. As shown in Fig. 6 the time required to complete 1 × 10−3 s of time integration for EWCSPH-IBM is not too far off from the CFX-IBM compute requirements, in serial the time required is 24.28 s and 21.49 s for EWCSPH-IBM and CFX-IBM respectively. Considering that the EWCSPH-IBM code is a self-written C++ code, which is far from optimised, and that R R ANSYS⃝ CFX⃝ is an industry leading highly optimised solver with virtually linear scaling (ANSYS user manual) it is fair to say that the EWCSPH-IBM computational requirements for this test could be lowered even further than those for CFX-IBM at the same spatial resolution. Also, on average the CFX-IBM only requires about 4 iterations to reach convergence meaning that only 4 sub-iterations are performed per ∆tstep while the EWCSPH-IBM performs 10 time steps per ∆tstep . This result is assuring as it means that employing an optimally parallelised ISPH solver, which would be more computationally expensive in comparison to the Eulerian WCSPH based solver but also eliminate the oscillatory pressure, could be a viable alternative to industry leading solvers offering greater accuracy at the same spatial resolution. The additional expense would be mostly due to solution of the Pressure Poisson Equation which comes at a significant computational overhead. The work of Guo et al. (2013) and Chow et al. (2018) show that ISPH can be parallelised nearly as effectively as WCSPH (Crespo et al., 2011) via solvers utilising MPI CPU parallelisation and GPU parallelisation, respectively. This however must be validated before concrete claims can be made.

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

275

Fig. 6. Average computational time required to complete 1 time step vs number of processors.

Fig. 7. Schematic illustrating the computational setup for the impulsively started plate test.

5. Investigation of EWCSPH and LWCSPH for FSI involving a thin rigid impulsively started plate The impulsively started plate test case has been investigated by a number of authors (see Favier et al., 2014, Koumoutsakos and Shiels, 1996, Yoshida and Nomura, 1985). It involves high levels of viscous shear formed due to the presence of large pressure gradients and the no-slip condition affected on the surface of a very thin plate. Previous LWCSPH studies involving FSI with immersed bluff bodies have been limited to circular and square bodies accelerated through a fluid. Hence, it is interesting to investigate this case to evaluate the ability of SPH to model the mixing of shear layers from a point and impulsively started flow. 5.1. Convergence study and validation of EWCSPH-IBM The test involves a rigid flat plate of length L = 0.3 m moving at a constant velocity through a fluid initially at rest. For this case the fluid domain is 6.66L long and 3L wide enclosed within periodic boundaries, Fig. 7 illustrates the setup for this case. The case is modelled at Reynolds number Re = Up L/ν = 126 and the plate is impulsively started at Up = 0.1 m/s so that it moves through the fluid at a constant velocity leaving a pair of counter rotating vortices in its wake which grow in size as the simulation progresses. For the tests in this section the Wendland smoothing kernel (Wendland, 1995) is used with h = 2∆x. The vortices are fed by strong shear layers which are formed on the two surfaces of the plate and act in opposite directions creating a source of vorticity at each tip of the plate as seen in Fig. 8 which shows snapshots of the vorticity field at different points in time plotted using vorticity contour levels of (−2z ,2z ) where z = (−8, −7, . . ., 3, 4). An issue encountered with this test is that the shock of the impulsive start can create a chequer board pressure field which develops due to the destructive interference of fast travelling pressure waves crossing the periodic boundaries. Hence, δ -SPH with ad = 0.1 is used for the full simulations presented in this section to diffuse the noise created by the chequer board effect

276

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

Fig. 8. Eulerian impulsively started flat plate test wake vorticity contours for Re = 126 with ∆x/L = 120. (a), (b), (c) and (d) show the contours at tU p /L = 0.5, 1.0, 2.0 and 4.0 respectively. The thick dark lines just behind the plate are a result of the vorticity contours merging into each other at the interface between the main vortex bubble and the shear layer in the wake of the plate.

which if left untreated could corrupt the flow field. Note; the aforementioned chequerboard effect is not the result of SPH numerical instability. The source of the chequer board pattern is discussed in detail in Nasar (2016) and also briefly in Section 4. For this case co = 50Umax where Umax here is the maximum observed fluid velocity (≈0.3 m/s) and the CFL coefficient was set to 0.5 leading to an average time step size of roughly 3.5 × 10−5 s for the test with ∆x/L = 120. The convergence behaviour of the EWCSPH-IBM formulation is assessed using multiple resolutions beginning with ∆x/L = 60 and ending with ∆x/L = 160. In the absence of a reference solution, results obtained using ∆x/L = 160 are used as the reference solution for the convergence study. For the convergence study the time step size is kept constant at ∆t = 2.5×10−5 s for all resolutions to keep temporal error small compared to spatial error. δ -SPH is removed only for the convergence study to prevent diffusive effects influencing the study as suggested by Roache (Roache, 1997). For all the resolutions the x-direction velocity u is measured at Ny = 32 intervals along the line x = 2.22 m which lies in the wake of the plate and the error is evaluated as

  Ny (  ∑ )2 1 p p e=√ unp =160 − unp . Ny

(44)

p=1

The measurements for u are taken at the same instant in time (t = 0.1 s) for all resolutions. The pressure wave interference issue only begins to affect the flow structures near the plate once the wave has re-entered the computational domain, which occurs after this point in time. Hence, δ -SPH was not required for the convergence study. As shown in Fig. 9 the error (e) follows second-order convergence. This agrees with the theoretical accuracy for both EWCSPH and the IBM interpolations using second-order accurate smoothing kernels and a uniformly distributed set of interpolation points. Higher-order approximations can be made using fourth-order or sixth-order Gaussian kernels with large radii of support (such as 6h) (Lind and Stansby, 2015). Fig. 10 shows snapshots of the vorticity field at t = 0.1 s for selected refinement levels (∆x/L = 80, 100, 120 and 140 markers). It shows that a region of counter rotating vorticity connected to the back of the shear layer in the wake (along the line x ≈ 2.23 m) disappears with smaller particle spacing suggesting that the lower resolutions were underresolved. A similar phenomenon was observed with under-resolved tests performed using ANSYS CFX-IBM for the flow past an accelerating square box presented in Section 4. It can also be noted from Fig. 10 that the shape of the vortex bubble does not change significantly for resolutions finer than ∆x/L = 120 hence this resolution was chosen to measure the development of the wake through a long duration simulation for comparison with reference results. Fig. 11 shows the development of the length of the vortex bubble s plotted against normalised time for the case with ∆x/L = 120 for Lagrangian markers. The length of the vortex bubble was measured as the point in the wake of the flat plate at which x-direction flow reversal occurs along the line defined by y = 0 m. As can be seen in Fig. 11 the results are in agreement with reference numerical results obtained using a vortex particle method by Koumoutsakos and Shiels (Koumoutsakos and

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

277

Fig. 9. Error e for EWCSPH with np = (60 to 110) plotted against particle spacing ∆x using log–log axis.

Fig. 10. Effect of resolution refinement on the vorticity measurements for the Eulerian impulsively started flat plate test. (a), (b), (c) and (d) correspond to np = 80, 100, 120 and 140 respectively. The vorticity contour levels are (−2z , 2z ) for z = (−8, −7, . . ., 3, 4).

Fig. 11. Evolution of the normalised length of the vortex bubble s/L in the wake of the plate plotted against normalised time for the impulsively started plate test with ∆x/L = 120.

Shiels, 1996), a finite element method by Yoshida & Nomura (Yoshida and Nomura, 1985) and experimental measurements obtained by Taneda and Honji (Taneda and Honji, 1971). These results further highlight the ability of EWCSPH-IBM to capture high vorticity flow fields to a good degree of accuracy and show that when uniformity of the particle arrangement is maintained SPH can be used to model flows involving high

278

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

Fig. 12. Impulsively started flat plate test. Magnified pressure and vorticity fields about the top half of plate. Images (a), (b), (c) and (d) show the fields at tU p /L = 0.04 for the case using only the MLS correction, only shifting, the MLS correction and particle shifting and EWCSPH with no corrections, respectively. Images (e), (f), (g) and (h) show the results for the case using only the MLS correction, only shifting, the MLS correction and particle shifting, and EWCSPH with no corrections, respectively at tU p /L = 0.24. In all above images the vorticity contour levels are (−2z , 2z ) where z = (−1, 0, 1, 2, 3, 4). Note that different scales are used for each of the pressure field snapshots to ensure that the scale reflected the pressure range predicted by each method at that point in time.

levels of viscous shear and recirculation with confidence in the numerical convergence of the method even for such a demanding test.

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

279

Fig. 13. Snapshots of the vorticity field for the impulsively started plate test after the plate has moved a significant distance from its initial position, tU p /L = 1.17. (a) Vorticity contours for the EWCSPH predictions. (b) Vorticity contours for the results obtained using the MLS-IBM Lagrangian WCSPH formulation with δ -SPH and shifting. The vorticity contour levels are (−2z , 2z ) where z = (−1, 0, 1, 2, 3, 4).

5.2. Comparison of the EWCSPH-IBM and the LWCSPH-IBM for the impulsively started plate problem For the Lagrangian WCSPH tests highly irregular particle distributions were expected due to the circulatory nature of the wake. Using different combinations of the corrective techniques discussed in Section 2 for WCSPH and Section 3.2 for the IBM an investigation was performed into which are most effective in suppressing any errors or instabilities caused by the disordered arrangement. In addition to those techniques the corrected kernel gradient, Eq. (18), was used for the SPH fluid–fluid interpolations to assess its effect on the results. However, the correction did not have a significant influence on the flow field. The size of the domain for the Lagrangian tests was reduced to save computational effort (width = 1.0 m and length = 1.0 m). Since the domain is smaller the periodic boundaries can affect the flow field near the plate more significantly than in the previous test. Hence the velocity of the fluid particles was initialised to U∞ = (0.01, 0.0) m/s and the velocity of particles re-entering the domain from the left boundary was constrained to U∞ = (0.01, 0.0) m/s. The prescribed re-entry velocity acts to damp pressure waves as they leave the domain along the x-coordinate and also maintains a uniform inlet velocity which is required due to the smaller domain size. The velocity of the plate was initialised and maintained at Up = (−0.09, 0.0) m/s to ensure that the relative normal velocity between the fluid inflow and the plate was (0.1, 0.0) m/s. Other variations of the ratio |Up |/|U∞ | were trialled and the flow fields were initially similar for all cases. However, decreasing |Up |/|U∞ | leads to less plate motion and more fluid motion which means that as the simulation progresses the wake ends up being dragged by the motion of fluid into the domain boundaries which stretches the vortex bubble. In cases with lower |Up |/|U∞ | values the initial position of the plate must be moved to be closer to the inlet periodic domain. The numerical speed of sound for these tests was set to co = 100Umax and the CFL number was set to 0.04 for the initial 6000 time integration steps in an attempt to prevent large voids forming in the computational domain due to the impulsive start of the test case, the CFL number was then increased to 0.1 for the remainder of the time integration steps. The particle spacing was set to ∆x/L = 60 and the test was re-performed using EWCSPH and the same computational setup to provide a benchmark for the Lagrangian tests. Both EWCSPH and LWCSPH are expected to perform better in general with spatial resolution refinement but the aim here is not to achieve converged results. The aim of this test is purely to highlight the capabilities of EWCSPH in comparison to LWCSPH hence this resolution is sufficient. A background pressure was also trialled as a means of regularising the particle distribution and to prevent the creation of voids as suggested by Violeau and Leroy (2014) and Adami et al. (2012) where the values used for the background pressure ranged from 10 Pa to 1000 Pa. However, the background pressure had little effect on improving the accuracy and stability of the scheme when using low values within the range tested while higher values led to increased spurious fluctuations in the pressure field. The first test was performed using δ -SPH with ad = 0.3 to suppress pressure fluctuations which occur due to the violent nature of the impulsive start and the disordered particle arrangement which is strongly influenced by vorticity; similar values were used by Vacondio et al. (2013) in their simulation of SPHERIC Benchmark 6. The MLS-IBM proposed by Cherfils et al. (2010) for Lagrangian SPH is also used to ensure first-order consistency for the IBM interpolations. For LWCSPH with δ -SPH and MLS-IBM corrections, despite the high value for the speed of sound and the low CFL number, particles in the wake separate from the plate creating a void as shown in Fig. 12(a). The fluid particles eventually catch up with the plate as shown in Fig. 12(e) however spurious pressure fluctuations are visible in the vortex near the tip and the fluctuations prevent the formation of a stable wake. For the second test the MLS-IBM was not used. Instead, in addition to δ -SPH, the particle distribution was regularised using the diffusion based particle shifting technique in order to minimise errors caused by particle disorder. A wide range of values was trialled for the diffusion controlling parameter λ in Eq. (17) (from 0.01 to 0.35) and the most stable wake was captured using λ = 0.25. Values lower than λ = 0.25 were not sufficient to minimise particle disorder and prevent

280

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

the formation of voids leading to large pressure fluctuations while values higher than λ = 0.25 led to excessive shifting corrupting the simulation. Considering the impulsive start and the high level of vorticity encountered during the test it is understandable why a high value for λ was required to maintain regularity. Also, it was found that correcting the particle properties after shifting using the non-conservative Taylor series extrapolation discussed in Lind et al. (2012) led to errors quickly eroding the quality of the results hence it was not used for this case. As seen in Fig. 12(b and f) the simulation using a combination of shifting and δ -SPH initially captures a relatively welldefined wake in good agreement with the EWCSPH results in Fig. 12(d and h). However, as the simulation progresses and the vortex bubble grows larger the shifting technique is no longer able to maintain a stable bubble and the wake disintegrates. To investigate if the wake disintegration was due to errors in the IBM interpolations the test was re-performed using, δ -SPH, the MLS-IBM correction and the particle shifting technique. Fig. 12(c) shows a snapshot of the results obtained using this setup where, initially, the pressure field is identical to the pressure field obtained using only shifting shown in Fig. 12(b). However, the snapshots at tU p /L = 0.24 Fig. 12(f and g) show that there is a small difference between the results where the results obtained using δ -SPH, MLS-IBM and shifting show a pressure field in slightly better agreement with the EWCSPH results in Fig. 12(h). This leads to a vorticity field in better agreement with the EWCSPH results. Spurious pressure fluctuations can be observed in all the LWCSPH results in Fig. 12 while the EWCSPH results are virtually noise free near the plate. Analysing these results, it can be concluded that the MLS-IBM corrections do slightly improve the pressure field. However, the major improvement is due to the particle shifting technique which regularises the particle distribution reducing the overall kernel derivative interpolation error. Hence the major source of error with SPH-IBMs is most likely related to errors in the SPH kernel derivative approximations for the interactions between fluid particles, Eqs. (3) and (4), and is not due to the accumulation of small errors in the IBM kernel interpolations, Eqs. (27) and (30), between the fluid and the boundary, as suggested by Cherfils et al. (2010) (that is they are less sensitive to particle disorder). Fig. 13 shows snapshots of the test performed using δ -SPH, MLS-IBM and shifting at a later point in time during the simulation (tU p /L = 1.17) illustrating the irregularity of the wake using LWCSPH compared to the stable wake predicted by the EWCSPH formulation. While the combined use of δ -SPH, the MLS-IBM correction and shifting reduces the particle distribution anisotropy and prevents the formation of voids which can completely corrupt the simulation from the outset (as shown in Fig. 12(a and e)), the errors introduced by the large magnitude of shifting required to keep the particle distribution uniform accumulate and eventually corrupt the simulation. Persistent pressure fluctuations in the shear layers attached to the plate also contribute to the corruption of the wake. The combination of δ -SPH, MLS-IBM and shifting was found to be the best combination of techniques for this type of test but at the same resolution for particle spacing the EWCSPH approach produces significantly more accurate and stable results. 6. Conclusions This paper has presented an accurate, numerically stable method for coupling WCSPH to the IBM. The impulsively started flat plate case was investigated using Eulerian WCSPH-IBM and the results are in good agreement with reference experimental and specialised numerical tests showing that the dynamic mixing of shear layers from the plate tip can be accurately captured when uniformity of the particle distribution is maintained. A grid convergence study was presented showing that the technique is second-order accurate in space which is in agreement with the theoretical order of accuracy. The results of studies performed using several Lagrangian WCSPH-IBM formulations show that irregularity in the particle distributions severely affects the methods’ ability to capture such flows due to spurious fluctuations in the flow properties. Using a combination of well-known treatments such as particle shifting, density diffusion terms and MLS smoothing kernel corrections improves the quality of the flow structures captured but eventually errors accumulate and corrupt the simulation leading to vortical structures with ill-defined contours in comparison to the results predicted using Eulerian WCSPH and the same resolution for the particle spacing. Further validation of the EWCSPH-IBM for the accelerated square box (SPHERIC Benchmark 6) showed that vorticity was accurately captured although pressure waves were evident due to compressibility, a feature of all weakly compressible solvers. For confined domains these waves reflect off the wall boundaries and interfere with the pressure field leading to small oscillatory forces. EWCSPH-IBM has similar computational requirements as an industry leading solver. There is some evidence that vorticity structures at higher Reynolds number investigated are better resolved. Importantly the method is inherently numerically stable and hence robust. Acknowledgements The lead author would like to thank the University of Tripoli, formerly the University of Al-Fateh in Tripoli, Libya and the Libyan Ministry of Higher Education for funding the PhD research project in which this work was developed. The authors would also like to thank the reviewers for their suggestions which significantly helped in improving this paper.

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

281

References Adami, S., Hu, X., Adams, N., 2012. Simulating three-dimensional turbulence with SPH. In: Proceedings of the Summer Program. Stanford, USA. Antoci, C., Gallati, M., Sibilla, S., 2007. Numerical simulation of fluid–structure interaction by SPH. Comput. Struct. 85 (11), 879–890. Basa, M., Quinlan, N.J., Lastiwka, M., 2009. Robustness and accuracy of SPH formulations for viscous flow. Internat. J. Numer. Methods Fluids 60 (10), 1127–1148. Batchelor, G.K., 2000. An Introduction to Fluid Dynamics. Cambridge University Press. Bonet, J., Lok, T.-S., 1999. Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Comput. Methods Appl. Mech. Engrg. 180 (1), 97–115. Canelas, R., Ferreira, R.M., Crespo, A., Domínguez, J., 2013. A generalized SPH-DEM discretization for the modelling of complex multiphasic free surface flows. In: The 8th International SPHERIC Workshop. Cherfils, J.-M., 2011. Développements et Applications de la Méthode SPH Aux écoulements Visqueux à Surface Libre (Ph.D. Thesis), Université du Havre, France. Cherfils, J.-M., Blonce, L., Pinon, G., Rivoalen, E., 2009. Towards the simulation of wave-body interactions with SPH. In: 4th International SPHERIC Workshop. Nantes, France. Cherfils, J.-M., Blonce, L., Pinon, G., Rivoalen, E., 2010. IB-SPH simulations of wave-body interactions. In: 5th International SPHERIC Workshop. Manchester, UK. Chow, A.D., Rogers, B.D., Lind, S.J., Stansby, P.K., 2018. Incompressible SPH (ISPH) with fast Poisson solver on a GPU. Comput. Phys. Comm.. Colagrossi, A., 0000. Bencmark Test 6: 2D Incompressible flow around a moving square inside a rectangular box. SPHERIC. Crespo, A.C., Dominguez, J.M., Barreiro, A., Gómez Gesteira, M., Rogers, B.D., 2011. GPUs, a new tool of acceleration in CFD: efficiency and reliability on smoothed particle hydrodynamics methods. PLoS One 6 (6), e20685. Fatehi, R., Fayazbakhsh, M., Manzari, M., 2008. On discretization of second-order derivatives in smoothed particle hydrodynamics. In: Proceedings of World Academy of Science, Engineering and Technology. Favier, J., Revell, A., Pinelli, A., 2014. A Lattice Boltzmann–Immersed Boundary method to simulate the fluid interaction with moving and slender flexible objects. J. Comput. Phys. 261, 145–161. Fourtakas, G., Nasar, A.M.A., Vacondio, R., Stansby, P.K., Lind, S.J., Guo, X., Rogers, B.D., 2018. Towards high-order 3-D Eulerian incompressible SPH for arbitrary geometries with generalised particle distributions. In: 13th SPHERIC Workshop. Galway, Ireland. Fourtakas, G., Rogers, B.D., Laurence, D.R., 2013. Modelling sediment resuspension in industrial tanks using SPH. Houille Blanche (2), 39–45. Gómez-Gesteira, M., Dalrymple, R.A., 2004. Using a three-dimensional smoothed particle hydrodynamics method for wave impact on a tall structure. J. Waterw. Port Coastal Ocean Eng. 130 (2), 63–69. Griffith, B.E., 2012. On the volume conservation of the immersed boundary method. Commun. Comput. Phys. 12 (2), 401–432. Guo, X., Lind, S., Rogers, B., Stansby, P.K., Ashworth, M., 2013. Efficient Massive Parallelisation for Incompressible Smoothed Particle Hydrodynamics (ISPH) with 10^8 particles. Hieber, S.E., Koumoutsakos, P., 2008. An immersed boundary method for smoothed particle hydrodynamics of self-propelled swimmers. J. Comput. Phys. 227 (19), 8636–8654. Karunasena, H., Senadeera, W., Gu, Y., Brown, R.J., 2014. A coupled SPH-DEM model for micro-scale structural deformations of plant cells during drying. Appl. Math. Model. 38 (15), 3781–3801. Khayyer, A., Gotoh, H., Shimizu, Y., 2017. Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context. J. Comput. Phys. 332, 236–256. Koumoutsakos, P., Shiels, D., 1996. Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate. J. Fluid Mech. 328, 177–227. Lind, S.J., Stansby, P.K., 2015. Investigations into high-order incompressible SPH. In: SPHERIC. Parma, Italy. Lind, S., Stansby, P., 2016. High-order Eulerian incompressible smoothed particle hydrodynamics with transition to Lagrangian free-surface motion. J. Comput. Phys. 326, 290–311. Lind, S.J., Xu, R., Stansby, P.K., Rogers, B.D., 2012. Incompressible smoothed particle hydrodynamics for free-surface flows: a generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. J. Comput. Phys. 231 (4), 1499–1523. Longshaw, S.M., Rogers, B.D., 2015. Automotive fuel cell sloshing under temporally and spatially varying high acceleration using GPU-based Smoothed Particle Hydrodynamics (SPH). Adv. Eng. Softw. 83, 31–44. Luis Cercos-Pita, J., 0000. SPHERIC test case number 6 (Flow around a moving square). AQUAgpusph. Marrone, S., Antuono, M., Colagrossi, A., Colicchio, G., Le Touzé, D., Graziani, G., 2011. δ -SPH model for simulating violent impact flows. Comput. Methods Appl. Mech. Engrg. 200 (13), 1526–1542. Mittal, R., Iaccarino, G., 2005. Immersed boundary methods. Annu. Rev. Fluid Mech. 37, 239–261. Mokos, A., Rogers, B.D., Stansby, P.K., Domínguez, J.M., 2015. Multi-phase SPH modelling of violent hydrodynamics on GPUs. Comput. Phys. Comm. 196, 304–316. Molteni, D., Colagrossi, A., 2009. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. Comput. Phys. Comm. 180 (6), 861–872. Monaghan, J.J., 1989. On the problem of penetration in particle methods. J. Comput. Phys. 82 (1), 1–15. Monaghan, J.J., 2005. Smoothed particle hydrodynamics. Rep. Progr. Phys. 68 (8), 1703. Morris, J.P., Fox, P.J., Zhu, Y., 1997. Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys. 136 (1), 214–226. Nasar, A.M.A., 2016. Eulerian and Lagrangian Smoothed Particle Hydrodynamics as Models for the Interaction of Fluids and Flexible Structures in Biomedical Flows (Ph.D. Thesis), The University of Manchester, UK. Nasar, A.M.A., Fourtakas, G., Lind, S.J., Rogers, B.D., Stansby, P.K., 2018. Towards higher-order boundary conditions for eulerian SPH. In: 13th SPHERIC Workshop. Galway, Ireland. Nasar, A.M.A., Rogers, B.D., Revell, A., Stansby, P.K., 2016. An eulerian weakly compressible smoothed particle hydrodynamics immersed boundary method for fluid structure interaction. In: 11th International Smoothed Particle Hydrodynamics European Research Interest Community (SPHERIC) Workshop. Munich, Germany. Noutcheuwa, R.K., Owens, R.G., 2012. A new incompressible smoothed particle hydrodynamics—immersed boundary method. Int. J. Numer. Anal. Model. 3 (2), 126–167. Panizzo, A., 2004. Physical and Numerical Modelling of Subaerial Landslide Generated Waves (Ph.D. thesis), Universita degli Studi di L’Aquila. Peskin, C.S., 1972. Flow patterns around heart valves: a numerical method. J. Comput. Phys. 10 (2), 252–271. Quinlan, N.J., Basa, M., Lastiwka, M., 2006. Truncation error in mesh-free particle methods. Internat. J. Numer. Methods Engrg. 66 (13), 2064–2085. Roache, P.J., 1997. Quantification of uncertainty in computational fluid dynamics. Annu. Rev. Fluid Mech. 29 (1), 123–160. Shadloo, M.S., Zainali, A., Yildiz, M., Suleman, A., 2012. A robust weakly compressible SPH method and its comparison with an incompressible SPH. Internat. J. Numer. Methods Engrg. 89 (8), 939–956.

282

A.M.A. Nasar, B.D. Rogers, A. Revell et al. / Journal of Fluids and Structures 84 (2019) 263–282

Skillen, A., Lind, S., Stansby, P.K., Rogers, B.D., 2013. Incompressible smoothed particle hydrodynamics (SPH) with reduced temporal noise and generalised Fickian smoothing applied to body–water slam and efficient wave–body interaction. Comput. Methods Appl. Mech. Engrg. 265, 163–173. Taneda, S., Honji, H., 1971. Unsteady flow past a flat plate normal to the direction of motion. J. Phys. Soc. Japan 30 (1), 262–272. Uhlmann, M., 2004. New results on the simulation of particulate flows. In: Centro de Investigaciones Energeticas Medioambientales y Tecnologicas (CIEMAT). Madrid (Spain. Uhlmann, M., 2005. An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209 (2), 448–476. Vacondio, R., Rogers, B.D., Stansby, P.K., Mignosa, P., Feldman, J., 2013. Variable resolution for SPH: a dynamic particle coalescing and splitting scheme. Comput. Methods Appl. Mech. Engrg. 256, 132–148. Van Liedekerke, P., Tijskens, E., Ramon, H., Ghysels, P., Samaey, G., Roose, D., 2010. Particle-based model to simulate the micromechanics of biological cells. Phys. Rev. E 81 (6), 061906. Vanella, M., Balaras, E., 2009. A moving-least-squares reconstruction for embedded-boundary formulations. J. Comput. Phys. 228 (18), 6617–6628. Violeau, D., Leroy, A., 2014. On the maximum time step in weakly compressible SPH. J. Comput. Phys. 256, 388–415. Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math. 4 (1), 389– 396. Xu, R., 2010. An Improved Incompressible Smoothed Particle Hydrodynamics Method and its Application in Free-surface Simulations (Ph.D. Thesis), The University of Manchester, UK. Xu, R., Stansby, P.K., Laurence, D.R., 2009. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. J. Comput. Phys. 228 (18), 6703–6725. Yoshida, Y., Nomura, T., 1985. A transient solution method for the finite element incompressible Navier–Stokes equations. Internat. J. Numer. Methods Fluids 5 (10), 873–890.