A.M.Salehizadeh and A.R.Shafiei

Received:07 November 2021/Accepted:09 February 2022

?Harbin Engineering University and Springer-Verlag GmbH Germany,part of Springer Nature 2022

Abstract A fully Lagrangian algorithm for numerical simulation of fluid-elastic structure interaction (FSI) problems is developed based on the Smoothed Particle Hydrodynamics (SPH) method. The developed method corresponds to incompressible fluid flows and elastic structures. Divergence-free (projection based) incompressible SPH (ISPH) is used for the fluid phase, while the equations of motion for structural dynamics are solved using Total Lagrangian SPH (TLSPH) method.The temporal pressure noise can occur at the free surface and fluid-solid interfaces due to errors associated with the truncated kernels. A FSI particle shifting scheme is implemented to produce sufficiently homogeneous particle distributions to enable stable, accurate, converged solutions without noise in the pressure field. The coupled algorithm,with the addition of proposed particle shifting scheme, is able to provide the possibility of simultaneous integration of governing equations for all particles, regardless of their material type. This remedy without need for tuning a new parameter, resolves the unphysical discontinuity beneath the interface of fluid-solid media. The coupled ISPH-TLSPH scheme is used to simulate several benchmark test cases of hydro-elastic problems.The method is validated by comparison of the presented results with experiments and numerical simulations from other researchers.

Keywords Smoothed particle hydrodynamics; Incompressible SPH; Total Lagrangian SPH; Fluidelastic structure interaction;FSI particle shifting scheme

1 Introduction

The analysis of multi-physics problems is important in various areas of engineering. Fluidstructure interaction(FSI) is such a multi-physics problem where the effect of the fluid on the structure could result in a significant struc‐tural deformation. This deformation affects the pressure and velocity of the fluid. Fluid-structure interaction prob‐lems are analyzed in numerical methods using two distinct approaches, a “direct” approach and a “partitioned” ap‐proach. In the first approach, the fluid and structure rela‐tions are expressed and solved simultaneously, done by Langer and Yang (2018). In the second approach, first, the fluid equations are solved for the determination of stress(traction or pressure)on solid boundaries,and then the cal‐culated stress is used as a boundary condition for the gov‐erning equations of the structure to determine the deforma‐tion of the solid media. In fact, in the “partitioned” ap‐proach, the information is transmitted from one media to another, so that the governing equations for both media converge simultaneously (Rao and Wan 2018; Salehizadeh and Shafiei 2021).

The main drawback in solving such problems with the standard Lagrangian Finite Element Method (FEM) is that severe mesh distortions, which accompany large deforma‐tions, lead to inaccurate results. In the case of the fluidstructure interaction, the coupling of fluid and structure media is usually obtained by Arbitrary Lagrangian-Eulerian (ALE) formulations for the fluid region, while the structure is modeled with the fully Lagrangian for‐mula, such as Aldlemy et al. (2020).Although re-meshing techniques and ALE formulations can overcome the mesh distortion issue, the remapping of state variables creates difficulties with history-dependent materials and the accu‐racy of results is questionable.Also,precision definition of the interface and free surfaces requires a particular remedy(Aldlemy et al.2020).

In the last decades, the mesh-free methods have been proposed to simulate fluid-structure interaction problems,which use a fully Lagrangian formulation to describe fluid flow and structure displacement(Idelsohn et al. 2008a;2008b; Nasar et al. 2019). These methods allow the mo‐tion of the boundary of fluid-structure interaction to be de‐termined without any special algorithm, such as the remeshing stage (He et al. 2018). Several studies have been carried out on the development of fluid-structure interac‐tion solvers using Lagrangian mesh-free methods.The par‐ticle modeling capability in the simulation of the breaking wave impact on a moving float is observed using Moving Particle Semi-implicit (MPS) by Koshizuka et al. (1998).Also, Lee et al. (2010) Described the MPS method for simulating the interaction of violent free-surface flow with structures. Khayyer et al. (2018b) coupled a fluid model based on MPS method with an MPS-based structure.They simulated the dam breaking with an elastic gate problem,experimented by Antoci et al. (2007) and hydro-elastic slamming of a marine panel.

On the other hand, the Smoothed Particle Hydrodynam‐ics (SPH) method has been developed as a mesh-free La‐grangian method by Lucy (1977) and Gingold and Monaghan (1977) to solve astronomical physics problems.This method has been used in simulating various applica‐tion problems such as turbulent flows(Kazemi et al.2017;2020), free-surface flows (Zheng et al. 2017; Salehizadeh and Shafiei 2019) and multi-material interactions (Sale‐hizadeh and Shafiei 2021).

Due to potential advantages of SPH method,a consider‐able number of studies have been conducted for develop‐ment of hydroelastic FSI solvers by coupling a SPH-based fluid model either with a mesh-based structure model such as SPH-FEM (Hermange et al. 2019), SPH-SFEM (Zhang et al. 2021b), SPH-ESFEM (Long et al. 2020) and SPHNMM (Xu et al. 2019). On the other hand, different FSI solvers based on pure particle methods have been devel‐oped. Zhang et al. (2021a) presented an open source SPH code named SPHinXsys to treat various FSI problems in‐cluding the dam break through an elastic gate. Sun et al.(2021) developed a FSI-SPH model that is effective to solve some challenging FSI problems in three-dimensional space.

The traditional SPH solver was formulated based on a weakly compressible fluid, which relates the pressure to the variation of density through an artificial equation of state, referred as “weakly compressible” approach(WCSPH), described by Liu and Liu (2003). The standard WCSPH suffers from the computational error of density variation causes non-physical oscillations of pressure that can produce numerical instability(Lee et al.2008).On the resolving of high-frequency pressure oscillations, exten‐sive prevention measures have been proposed. Huang et al. (2018) employed a artificial diffusive term (δ-SPH model), proposed by (Antuono et al. 2010), to treat wave interactions with structure in field of ocean engineering.This approach requires a small time step with regard to the“sound speed”, which is at least ten times higher than the maximum velocity of fluid, thus increasing computation time. To prevent such issues, the concept of “truly incom‐pressible fluid”was achieved for SPH method.The Incom‐pressible smoothed particle hydrodynamics method(ISPH) is a two-step method that first calculates the pre‐dicted velocity,and then the predicted values are corrected using the pressure derived from the solution of the Poisson equation to apply a divergence-free velocity field, done by Cummins and Rudman(1999).

The ISPH method can be accurate and effective in re‐ducing the pressure fluctuations in the simulation of the fluid-solid interaction problems. Rafiee and Thiagarajan(2009) coupled the ISPH method for a fluid with an SPHbased structure. In their model, the fluid and solid pres‐sures were determined by solving the pressure Poisson equation (PPE) with a simple scheme. They assumed that the pressure time variations are negligible between two computational time steps, which is not reasonable.Khayyer et al. (2018a) presented an enhanced ISPH-SPH coupled method for simulating fluid-structure interaction.They used a dynamic stabilizer scheme,introduced by Tsu‐ruta et al.(2013),for fluid that involves additional calcula‐tions to determine parallel and normal vectors of predicted relative distances.Recently,their model have been demon‐strated to be capable of reproducing hydroelastic FSI cor‐responding to composite structures(Khayyer et al.2021).

The SPH method is applied to simulate the structural dy‐namics, such as in the Libersky et al. (1993) and Gray et al. (2001) studies. Antoci et al. (2007) used the WCSPH method to solve hydro-elastic problems. The original for‐mulation of WCSPH method suffers from a series of insta‐bilities and low accuracy issues when applied to solid me‐chanics (Swegle et al. 1995; Rabczuk and Belytschko 2004).The first issue SPH suffered from was the so-called“tensile instability”which leads to particle sticking and nu‐merical fracture. To reduce the effects of tensile instabili‐ties in solving such problems, several treatments have been proposed such as normalizing kernel functions(John‐son and Beissel 1996) and applying artificial stress (Gray et al. 2001). Nevertheless, some techniques add numerical parameters that need to be carefully selected. To eliminate instabilities related to rank deficiency in particle-based FSI modeling, Rabczuk et al. (2010) applied stress point integration (Rabczuk and Belytschko 2004) into FSI for fracturing shell structures. The stress-points implemented in mesh-free methods require a background grid for nodal integration;hence the method is not truly mesh-free.In the stability analysis of numerical mesh-free methods, Be‐lytschko et al.(2000)found tensile instability to be caused by the use of Eulerian kernel functions. They proposed a fully Lagrangian particle method in which the initial con‐figuration is considered as a reference, the smoothing ker‐nel function and its derivatives are calculated based on the initial distribution of the particles,making Lagrangian ker‐nel functions.In the present study,a totally Lagrangian for‐mulation proposed by Lin et al.(2015)is used to the simu‐lation of two-dimensional elastic structures. It can directly impose the interaction values such as force and velocity to the structure particles in the interaction with the fluid.The performance of the structure model is validated by simulat‐ing a cantilever beam under the concentrated force and the dynamic response of a free oscillating plate. The numeri‐cal results are compared to their analytical solution. Sun et al. (2019) introduced a coupled FSI-SPH solver by combining a multi-phaseδ-SPH scheme and a Total-Lagrangian SPH (TLSPH) solver for elastic structures. In their study, the numerical results of TLSPH demonstrated sufficient accuracy and therefore, this structural solver could be further implemented in the FSI simulations. Re‐cently, Zhan et al. (2019) and O’Connor and Rogers(2021) developed a GPU-accelerated coupled total La‐grangian and weakly compressible SPH(TL-WC SPH)ap‐proach for 3D fluid-structure interaction modeling.

The temporal pressure noise due to irregular particle dis‐tribution through clustering or stretching, leads to numeri‐cal instability. It also causes inaccurate behavior of fluid and structure in interaction with each other. A study con‐ducted by Xu et al. (2009), enhanced particle uniformity by shifting particles slightly away from the streamlines.The hydrodynamic properties are modified after they are shifted to account for their new spatial positions by using Taylor expansion.Their method provides a stable and con‐verged solution for fluid flow problems especially at rela‐tively high Reynolds numbers. However, this algorithm works inaccurately in the simulation of free-surface flows due to errors associated with the truncated kernels.Lind et al. (2012) proposed a relation based on the reference con‐centration gradient to control normal diffusion for freesurface particles. This remedy requires the setting of the constant parameters that encounter numerical challenges in long-term simulations. To achieve an accurate and con‐sistent implementation of particle shifting for free-surface flows, Khayyer et al. (2017) proposed a so-called Opti‐mized Particle Shifting (OPS) scheme through a careful implementation of only tangential shifting for free-surface and free surface vicinity particles.

In this research,for the first time,a novel FSI algorithm is introduced such a way that the proposed FSI particle shifting method added to FSI coupling scheme based on ISPH approach in order that the interaction term is ob‐tained and then imposed on the particles. This algorithm improves the interface behavior between fluid and struc‐ture.To the best of our knowledge,this is the first applica‐tion of the particle shifting treatment scheme, together with the coupled algorithm into the context of simulations of FSI problems by combining an ISPH scheme and a Total-Lagrangian Particle (TLP) solver for deformable elastic structures. The developed SPH coupling method was used to simulate the interaction of fluid-elastic struc‐tures for various problems and was compared and vali‐dated with experimental results. The considered problems of fluid-structure interaction in this paper include the de‐flection of an elastic plate due to the hydrostatic weight of water column, done by Fourey et al. (2017), the dambreaking flow through an elastic gate, done by Antoci et al.(2007),breaking dam flow on a hypo-elastic baffle,ex‐perimented by Liao et al.(2015),the deflection of an elas‐tic baffle due to fluid sloshing in a rolling tank, experi‐mented by Idelsohn et al.(2008a).

2 Numerical procedure

In the SPH formulation, the particle approximation of variableA(r)is determined by a summation of particles in‐side the support domain of the particle located atra:

whereVbis the volume of particleb,his the smoothing length which represents the discretization scale of SPH ap‐proximations and|rab|= |ra?rb| is the distance between particlesaandb.In this paper,the cubic B-spline kernel is used for all cases, reviewed by Liu and Liu (2003). A smoothing length ofh=1.2Δxis typically used, where Δxis the initial particle spacing. HereafterW(|rab|,h)will be simply written asWab.

The SPH discretization of gradient, divergence and La‐placian operator for an arbitrary scalar functionfor tensor functionFare,respectively,calculated as:

and?ais a renormalization tensor for the second deriva‐tive given by Fatehi and Manzari(2012):

where

and?defines the dyadic product of two vectors.

2.1 The projection-based incompressible SPH algorithm

In the SPH method, the Navier-Stokes equations are solved in the projection-based framework, as expressed in the Lagrangian form:

whereuis the velocity field,ρis the fluid density,pis pressure,vis kinetic viscosity,gis the gravity acceleration andvtis the eddy viscosity due to spatial filtering. In Eq.(9),the accelerationas→fis due to the interaction force on the fluid by the structure. We used the eddy viscosity ap‐proach based on the sub-particle scale turbulence model presented by Gotoh et al.(2000):

whereε?ijis the strain rate tensor;csis the Smagorinsky constant which is considered to be 0.1 for all simulations and Δ is the filter width, which is defined to be the SPH kernel smoothing lengthh.The ISPH method has been in‐troduced by Cummins and Rudman (1999) and has been further developed by other researchers(Shao and Lo 2003;Hu and Adams 2007; Lind et al. 2012). In this approach,both the mass conversation equation and the equation of state are replaced by a Poisson equation, which is deter‐mined using the projection method. In general, to deter‐mine an intermediate velocity field, the momentum equa‐tion is solved without the effect of a pressure gradient.Then a Poisson equation for pressure is obtained by setting the divergence of intermediate velocity equal to the diver‐gence of pressure gradient.The final velocity at the end of the time step is achieved in such a way that it results in in‐compressible conditions.

The ISPH solver in its original formulation suffers from the density error accumulation and errors in velocity diver‐gence field. Several studies have been conducted to pre‐vent such issues. Shao and Lo (2003) proposed an algo‐rithm by imposing invariant density instead of zero veloc‐ity divergence that, despite its stability, is not accurate. Hu and Adams (2007) developed an iterative algorithm by si‐multaneously enforcing invariant density and zero velocity divergence that, despite its stability and accuracy, is very time-consuming. Xu et al. (2009) presented a new method to prevent the instability resulting from the intense cluster‐ing of particles. In their method, at the end of each time step, the particles are shifted away from the streamline leading to the uniform distribution of particles. Lind et al.(2012)proposed a generalized version of this scheme to al‐low extended applications to free-surface flows. Skillen et al.(2013)proposed a new method in the ISPH simulations for reducing the temporary disturbances caused by free surface effects. Their method provided a smoothing Pois‐son equation using a proposed relationship for adjacent free-surface particles.Also,the standard GPU implementa‐tion have been presented for the ISPH simulations, where the pressure Poisson equation is also solved on the GPU(Chow et al.2018).

Here a first-order time marching scheme is applied,where both the density and mass of particles are constant;The particle positions are predicted using current velocity field:

The intermediate velocity field is calculated based on the momentum equation,regardless of pressure gradient term at the positionr:

where

The SPH discretized form of Eq.(14)is calculated as:

Skillen et al. (2013) proposed a suitable method for ob‐taining a smoothing pressure field at the free-surface. By identifying the particles on the free-surface and adjacent to the free-surface, a variable correction coefficient is used which leads to effect of free-surface gradually. The freesurface particles are detected by calculating the divergence of the position vector of particles. By detecting the freesurface particles, these particles are given a zero pressure over a time step (Dirichlet condition in the Poisson equa‐tion of pressure). After the discretization of the Eq.(15),the Poisson equation of pressure is rewritten in the follow‐ing form:

whereAabare the Poisson’s coefficients of pressure and the summation is applied to all particles located in the support domain of the particle a. It is noted that the calculation of new pressures(pn+1)requires a limited number of iterations to reach the smoothing converged pressures. The coeffi‐cient appearing in Eq.(16) is proposed by Skillen et al.(2013) for the implementation of the effect of free surface conditions,which follows:

Then,by calculating the pressure of particles,the veloc‐ity at the end of the time step is determined from the inter‐mediate velocity and the negative version of pressure gra‐dient:

In the end,the position of the particles is updated as fol‐lows:

In order to stabilize the simulation and form a uniform distribution of particles, after the particle positions are ad‐vanced in time by Eq.(19), following Xu et al. (2009), the particles are shifted slightly, then the hydrodynamic vari‐ables are corrected by the Taylor series approximation:

whereAis a general variable;aandα′are the particle’s old position and new position respectively;δraa′is the displace‐ment vector between the particle’s new position and its old position.By modifying the particle shifting magnitude,ζ,in relation to the particle convection distance and the parti‐cle size,the position shift reads:

whereCis a constant in the range of 0.01?0.1;ζis the shifting magnitude which is equal to the maximum particle convection distance |u|maxΔt, with |u|maxthe maximum particle velocity, and Δtthe time step;Rais solved by fol‐lowing equation:

whereNais the number of neighboring particles around particlea,is the distance between particle a and particleb,andnabis the unit displacement vector between particle a andb;rˉais the average particle spacing in the neighbor‐hood ofa,and reads:

Lind et al. (2012) observed that due to the motion of free-surface particles,the distance of particle displacement would be much greater than the smoothing length which results in an irrational error; therefore, an upper limit for the particle shifting distance was defined as 0.2 h. The in‐complete kernel support and increased kernel interpolation error in free-surface regions, lead to error in determination of the displacement vector of particles.

According to Khayyer et al. (2017), in the area near the free surface,the shifting distance of the fluid particles was forced to be equal to zero in the normal direction of the free surface.To achieve a reliable physical particle shifting for free-surface and nearby particles, the particle shifting displacement vector for these particles is modified by Khayyer et al.(2017)as follows:

where Δrais obtained by Eq.(21).For the free-surface par‐ticles,the normal vector is determined by:

where

Figure 1 The particle shifting process for a free-surface particle(Khayyer et al.2017)

Figure 1 shows the particle shifting process for a free surface particle. The particles as depicted in Figure 2, are determined in the simulation as internal particles, particles adjacent to the free surface, free-surface particles and splash particles based on the divergence of the position vec‐tor of particles.At first, the free-surface particles are deter‐mined based on the divergence of position vector ??r＜xand the normal vector is obtained for them.Then the adja‐cent particles to the free surface are determined based on the establishment of two criteria of the divergence of posi‐tion vectorx＜ ??rand the spatial distance|ra′b|＜h(a′nearest neighbor particle in the free surface to the particle b)and the normal vector of the nearest neighbor particle located on the free surface is assigned to these particles. Splash par‐ticles have the divergence of position vector ??r＜xwhich do not have any neighboring particles in their support do‐main. By detecting free-surface and nearby particles and also splash particles, the internal particles are distin‐guished. It should be noted that for splash particles, the correction coefficient is zero. In the present study, the value of x was chosen to be 1.4.

Figure 2 Particle status identification to calculate the vector of particle displacement

Within discontinuous regions on the surface of the struc‐ture, the large displacements lead to instability, which re‐sult in forming jamming clumps due to disrupt of the fluid particles adjacent to the surface of the structure. To solve this problem, a novel method is proposed to shift the fluid particles adjacent to the surface of the structure.An inter‐mediate boundary layer is characterized in the fluid do‐main on surface of the structure as shown in Figure 3.

Figure 3 Definition of intermediate boundary layer in fluid domain

When the fluid particles are located in this intermediate boundary layer, these particles are displaced based on a predefined nonlinear function. This method results in a preserved solution with smoothing the velocity field.Relo‐cation of fluid particles located in the intermediate bound‐ary layer is executed using an exponential function.

whereδraa′?nais prescribed displacement vector by Eq.(21) in local normal directionnof the boundary surface anddis a flow dependent vector in local normal directionn.Considering a large depth for the intermediate boundary layer,|d|,enhance the stability of the solution with increas‐ing the number of particles within the layer that need to be relocated. Applying this relocation is determined with de‐fining limits for |d|. If the value ofδraa′|nis assumed to be equal to|u|maxΔtandδraa′|newis limited to the shifting dis‐tance proposed in Eq.(21) in local normal directionn, Eq.(27)becomes:

By applying second order Taylor series to scalar form of Eq.(28),it becomes:

whereCis a constant parameter in the range of 0.01?0.1,limits for|d|can be determined:

Finally, the hydrodynamic variables of modified at the new position as follows:

2.2 Total lagrangian formulation for elastic structures

The Total Lagrangian SPH(TLSPH)method is applied by Lin et al.(2015)to model the elastic structures;in this formu‐lation,the density is considered constant and to solve the mo‐mentum equation in the reference framework,the first Piola-Kirchhoff stress tensorPsreplaces the Cauchy stress tensorσ.Therefore,Eq.(9)is given in the following form:

In the above equation,αf→sis the acceleration due to the interaction force on the structure from the fluid side.The relation betweenPsandσis as follows:

It is necessary to obtain first-order consistency to avoid numerical errors caused by the truncated kernel function on the structure.For this purpose,the shape tensor is deter‐mined by Rabczuk et al. (2004) to reproduce the SPH ker‐nel gradients correctly,as follows:

where ?0Wabis the gradient of kernel function corresponds to particles a and b in their initial configuration as the sub‐script 0 has indicated, andXis a coordinate in the refer‐ence configuration.Nbrepresents the number of structure particles in the support domain of particlea.The deforma‐tion gradient tensorFis approximated as follows:

wherexis coordinate vectors in the current configuration.The Cauchy stressσij(iandjrepresent coordinate vectors)based on Eulerian strainεijis determined in two dimen‐sions by Hooke’s law for an isotropic elastic material:

For the plane-strain problem and,

For the plane-stress problem. Accordingly,vsis Pois‐son’s ratio and,Esis Young’s moduli.The Eulerian strains are defined using the deformation gradient tensor:

Regarding non-linear geometric behavior, the Green-Lagrange strain tensor is used,as follows:

In the above relation, the displacement gradient tensor is based on the displacement vector:

The TLSPH scheme suffers from rank-deficiency of zero energy modes due to its particle arrangement nature.Therefore, the discrete form of Eq. (33) is written in the initial configuration using the first Piola-Kirchhoff stress tensor and the interaction acceleration term.

In the present study, the 4th-order Runge-Kutta scheme as a multi-stage method is used to integrate the momentum equation (Eq. (42)) in time. This time integration scheme is fully explicit and is expressed as:

wherekumandkrmare solved simultaneously for each stage:

where Δtis time step anda=du/dt.

2.3 The fluid-structure interaction treatment for ISPH-TLSPH model

Although the SPH method is used to solve the govern‐ing equations for both fluid and structure media, two nu‐merical separated approaches for these media are consid‐ered. Therefore, they must be coupled together in a physi‐cally manner and applicable mathematics. In the follow‐ing, the adopted algorithm for coupling is expressed and can be seen in Figure 4.

Throughout the interface between fluid-structure, the normal components of the velocity vector must be continu‐ous. For this purpose, the use of the velocity field of the structure provides a suitable boundary condition for the fluid flow. In other words, the structure particles are used as moving boundaries for fluid governing equations to complete the support domain of fluid particles.

Essentially,the normal component of Cauchy stresses in fluid and structure are continuous on the boundary of fluidstructure interaction.In assuring the normal stress continu‐ity at the fluid-tructure interface, The pressure (an isotro‐pic section of Cauchy stress)is considered continuous dur‐ing the interaction between fluid and structure. Using the projection-based particle method, the ISPH algorithm is applied to all particles adjoining each other, regardless of their nature (belonging to the structure or fluid media). In other words,the position and velocity of structure particles as a boundary condition are used to calculate the fluid pressure field using the pressure Poisson equation (Figure 4b). The present algorithm accompanied by the proposed FSI particle shifting scheme ensures that the fluid velocity field remains free-divergence. When the fluid pressure field is calculated, the interaction term in the equation of fluid and structure momentum is determined based on the pressure gradient between the fluid and structure particles.Hence, the components of the interaction force, perpen‐dicular to the interface of the fluid and structure, are equal in the opposite direction.According to this:

Figure 4 Schematic View of the coupling model between fluid and structure

Regarding Neumann boundary condition for the pres‐sure field across the interface of fluid-structure interaction and considering a no-slip boundary condition at the fluidstructure interface, the Eq.(46) shows the relation of nor‐mal stress continuity at the interaction boundary. By using optimized particle shifting method, the fluid particles are then shifted slightly, also, the pressure field of all particles regardless of their nature are corrected by the Taylor series approximation,

The modified pressure field of particles is used to calcu‐late pressure gradient at a typical structure particle and then imposed on the structure particle as interaction termαf→s. This proposed coupling scheme works in a robust and efficient manner so that,no internal iteration for a con‐verged particle position is needed.The interaction term ob‐tained on the basis of the pressure gradient for a structure particle is as follows:

wherekdenotes all the particles of fluid and structure in the support domain of structure particle(Figure 4(c)).

The exchange of kinematic information throughout the interaction thus leads to a coupled process for a problem with material discontinuity. The pseudo-code of the pres‐ent solution algorithm is shown in Algorithm1 in the Ap‐pendix.

2.4 Time-stepping condition

In SPH method,the time step is generally determined in the Courant-Friedrich-Levy conditions:

Where Cr is Courant number andUis the characteristic ve‐locity of the model. The SPH-based structure model and the fluid model utilize different characteristic velocities in their formulation.The sound speed is used in the structure model, while for fluid, the maximum velocity of the fluid is used in each time-step. Consequently, the time steps al‐lowed for a structure should be quite smaller than the fluid.Therefore,significant improvements are achieved us‐ing separate and different time steps for both media. In this study,the time step of the structure Δtswas set to 0.01 Δtf. Therefore, the fluid governing equations are updated for every 100 steps of the structure. Here, the Courant number is taken as Cr=0.1 for all cases.

The proposed scheme is programmed in Fortran and an Intel?Core?i7-4770 CPU is used to run the computations.

3 Validations and comparisons

3.1 Benchmark tests for the structural analysis

To validate the structural analysis method,static and dy‐namic simulations were performed on an elastic structure.In a static simulation, a concentrated load is applied to clamped beam at its end.The geometric dimensions of the clamped beam are shown in Figure 5, whichL=100 mm is the length of the beam,b=1 mm the width andT=10 mm thickness of the beam. The load is gradually increased ranging from 0 to 1750 N withint=1.5 ms; then the load remains constant until the simulation ends att=3 ms.Young’s moduli and Poisson’s ratio for the beam were con‐sideredE=210 GPa andv=0.3 respectively.

Figure 5 The sketch of static problem of a clamped beam under concentrated load

The maximum displacement of the beam under the con‐centrated load at its end is obtained by:

A sensitivity analysis of the particle discretization is car‐ried out.The Different particle resolutions have been used,starting fromNy=5 particles through the thickness untilNy=17 particles.For every case,the solution ΔTLSPH/Δ is calcu‐lated and evolved vs number of particles as depicted in Figure 6. As one can observe, the predicted TLSPH enddeflection converges to the analytical solution as the num‐ber of particles is over ten, the error between the end de‐flection obtained using the TLSPH model and the analyti‐cal value is less than 1.04 percent,which is acceptable.

Figure 6 The ratio of predicted TLSPH end-deflection to the analytical solution vs the number of particles through the thickness

Table 1 presents the Root Mean Square Error(RMSE)of the results obtained from the simulations by TLSPH model presented in Figure 6.

Table 1 RMSE corresponding to deflection of a clamped beam subjected to a? concentrated load

Figure 7 shows the vertical displacement of the beam.The presented result has an acceptable correspondence with the analytical solution.

Figure 7 The deformed configuration of the clamped beam and the vertical displacement distribution of the beam

The present structure model was also used to simulate the free oscillation of an elastic plate of lengthL=0.2 m and thicknessT=0.02 m with a clamped edge and free edge(Fig‐ure 8).The initial vertical velocity distribution for the elastic plate is considered as follows(Gray et al.2001):

where

Figure 8 Schematic sketch of the dynamic problem of a free oscillating plate

wherekL=1.875.V0L=0.01 is the initial velocity at the free end andc0=K/ρis the sound speed,whichKis the bulk modulus of elastic plate. For comparison with the results of Antoci et al. (2007), the properties of the elastic plate were considered to beρ=1 000 kg/m3,K=3.25 MPa andv=0.397 5. For resolution test, the particle spacings of Δx=0.002 m, 0.001 m and 0.000 5 m are chosen which corre‐spond to 10, 20 and 40 particles in vertical direction. The simulation results are compared with analytical solution in Figure 9. The plate tip’s vertical displacement in Figure 9 shows that the present simulation is well-resolved forT/Δx=20.

Table 2 provides RMSE and Normalized Root Mean Square Error(NRMSE)corresponding to the results in Fig‐ure 9, quantitatively verifying the convergence property of the TLSPH model.

In Table 3,the simulation results are compared with ana‐lytical solution and numerical results by other researchers who had performed this test using other approaches.

Figure 10 shows the deformed configuration of the elas‐tic plate at the maximum vertical displacement of the plate at the instantt=0.65 s(Figure 10(a)).

Figure 9 Time histories of oscillations of the free end of the elastic plate reproduced by TLSPH with a set of different particle spacing with corresponding analytical solution

Table 2 RMSE and NRMSE corresponding to dynamic response of a free oscillating plate for 0.7 s of computation—dynamic response of a free oscillating plate.

Table 3 Comparison of dynamic test results with analytic solution and other researchers’results

In order to investigate the energy conservation property of 2D TLSPH model, time histories of kinetic, elastic strain, total energies untilt=0.7 s are presented in Figure 11. The presented time histories correspond to par‐ticle spacing of Δx=5.0e-4 m. From the presented figure,kinetic and elastic strain energies consistently vary in time with a phase difference of half of the period of plate’s os‐cillatory motion

3.2 Benchmark tests for the fluid-structure interaction

The proposed algorithm for analyzing fluidstructure in‐teraction problems is evaluated through several examples of applied problems.

Figure 10 Dynamic response of elastic plate under free oscillation

Figure 11 Time histories of kinetic,elastic strain,total energies by 2D TLSPH—dynamic response of a free oscillating elastic plate

3.2.1 Hydrostatic column of water on aluminum elastic plate

Figure 12 shows the initial configuration of the prob‐lem. An elastic aluminum plate with a widthL=1m and a thicknessT=0.05 m clamped on both sides and the column of water to the depth ofH=2 m with a similar width at the top of the plate. Water pressure variation is due to the ef‐fect of gravity. For the fluid, density and dynamic viscos‐ity wereρw=1000 kg/m3andμω=1.0×10-3Pas respectively,and the physical parameters of the aluminum plate were consideredρs=2700 kg/m3,E=67.5 GPa andv=0.34.

First,the aluminum plate bends due to the weight of the water column and the weight of the plate. Then, with the energy transfer between the potential energy and the ki‐netic energy, the elastic plate vibrates around its equilib‐rium state. The average mean descent of the elastic plate with the assumption of no damping in the structure and no viscosity in fluid,is calculated based on static analysis:

Figure 12 Schematic sketch of hydrostatic water column on an aluminum elastic plate

Figure 13 presents a set of time histories of plate’s midspan deflectionUycorresponding to three different spatial resolutions (initial particle spacing) by refining the resolu‐tion fromT/Δx=8 toT/Δx=12. The accuracy of deflection time history increases and eventually, forT/Δx=12 is quite well consistent with theoretical static displacement. In the beginning, the aluminum plate falls under sudden load caused by the weight of the water column. In the follow‐ing, the plate vibrates in response to this extensive load.Progressively, the plate damped and remains close to its static point.The acceptable stability was obtained by using the fluid-structure interaction algorithm based on the SPH method and a noise-free pressure field was provided.

Figure 13 The time histories of deflection at the midspan of an elastic plate under hydrostatic water column simulated by ISPHTLSPH FSI solver with respect to the analytical solution

Table 4 presents the RMSE corresponding to numerical results shown in Figure 13.

Figure 14 shows the vertical displacement of the center of the elastic plate Uy during the time forT/Δx=12 as com‐pared to the Fourey et al.(2017)results.

Figure 15 shows the fluid pressure field and stress distri‐bution of the structure reproduced by ISPH-TLSPHcoupled method att=0.3s.Thecomputational conditions for the correspondingsimulationareT/Δx=12 and Δt=8.0e-6 s.

Table 4 RMSE (Root Mean Square Error) corresponds to deflection of the mid-span of plate under the water hydrostatic column for 1.0 s of computation

Figure 14 The time histories of deflection at the mid-span of an elastic plate under hydrostatic water column

Figure 15 The fluid pressure field and stress distribution of the structure at t=0.3 s

In addition, Table 5 shows the required CPU time for ISPH fluid model and SPH structure model for one second of simulation of the hydrostatic water column on elastic plate for the resolution ofT/Δx=12.

Table 5 Required CPU time for simulation of hydrostatic water column on an elastic plate by using ISPH-TLSPH FSI solver until t=1 s for the resolution of T/Δx=12 (s)

3.2.2 Dam breaking through an elastic gate

The problem of water dam-breaking through an elastic gate was considered in comparison with the experimental results presented in Antoci et al. (2007) and the numerical results (Antoci et al. 2007; Rafiee and Thiagarajan 2009;Khayyer et al. 2018a; Fourey et al. 2017). The initial con‐figuration of the problem is shown in Figure 16.The water is located inside a reservoir that can flow by pushing a de‐formable gate that is clamped from above to a rigid wall.The elastic rubber gate is of thicknessT=5 mm with den‐sityρgate=1 100 kg/m3, Young modulusE=12 MPa, and Poisson ratiov=0.4.

Figure 16 Schematic sketch of the test case: breaking-dam flow through an elastic gate

The comparison between experimental results and nu‐merical simulation is shown in Figure 17 for every 0.04 seconds.The hydrostatic pressure of water causes the gate to gradually open, resulting in water flowing. By increas‐ing the gate opening,the amount of water leakage is inten‐sified. With the passing of time at the instantt≈0.14 s, the pressure applied on the back of the gate decreases due to the lowering of the water surface. Therefore, the gate slowly returns to its initial position.

The distribution of stressesσxxandσyyfor the elastic structure at its maximum vertical displacement at the in‐stantt≈0.14 s is shown in Figure 18.

Figure 17 Qualitative comparison between simulation snapshots and experimental photographs, obtained by Antoci et al. (2007) for every 0.04 second in breaking-dam flow through an elastic gate

Figure 18 The distribution of stresses σxx and σyy for the elastic structure reproduced by ISPH-TLSPH coupled method at its maximum displacement

Figures 19(a)and 19(b)show the horizontal and vertical displacement of the free end of the gate in time to compare with the results of Antoci et al. (2007) (experimental,WCSPH), Rafiee and Thiagarajan (2009) (ISPHWCSPH), Khayyer et al. (2018a) (ISPH-SPH) and Fourey et al.(2017)(SPH-FEM).As can be seen,at the momentt≈0.14 s, the elastic plate reaches its maximum vertical dis‐placement. The results show that the proposed coupled al‐gorithm for structure-fluid interaction analysis provides ac‐ceptable accuracy for determining the displacement of the elastic gate. It should be noted that the difference between the results of the present study and the experiments of An‐toci et al. (2007) is due to the nonlinear behavior of the rubber gate in the experiment. Comparison of the results of the present study with other studies shows that the dif‐ferent treatment of incompressible constraint of fluid in weakly compressible approach and truly incompressible approach (based on the SPH method) and different struc‐tural models lead to a difference in the determination of the amount of displacement of the gate.

Figure 19 Time histories of the displacement of the gate’s free end; present results by ISPH-TLSPH model vs. experimental results(Antoci et al.2007)and numerical results(Antoci et al.2007;Rafiee and Thiagarajan 2009;Khayyer et al.2018a;Fourey et al.2017)

The energy conservation property of ISPH-TLSPH FSI solver is investigated in the benchmark test of dam break‐ing through an elastic gate by considering the energy com‐ponents corresponding to fluid and structure.The effect of elastic energy of fluid due to compressibility is neglected by considering an incompressible fluid within the context of ISPH. The time evolution of the mechanical energy of the coupled system in each sub-domain is presented in Fig‐ures 20(a)?20(e) as well as the energy accumulated at the fluid-solid interface in Figure 20(f).

Figures 20(a), 20(b) and 20(c) show the time variations of gravitational potential energy of structure, the kinetic energy of structure and elastic strain energy of the struc‐ture, respectively. Figure 20(d) illustrates the time evolu‐tion of the mechanical energy of the coupled system in structure sub-domain. Figure 20(e) shows the time evolu‐tion of the mechanical energy of the coupled system in fluid sub-domain. From Figure 20(e), this drop of energy is mainly linked with the drop in gravitational potential en‐ergy of fluid.Figure 20(f)presents the time histories of to‐tal energy of the FSI system.A slight decreasing of the to‐tal system energy is observed that can be interpreted as a slight numerical dissipation. Indeed, within the context of ISPH, the fluid is considered to be incompressible and thus,the elastic energy of fluid due to compressibility is ig‐nored. For the considered benchmark test, the proposed ISPH-SPH FSI solver is shown to possess acceptable en‐ergy conservation properties.

3.2.3 Breaking dam on a hypo-elastic baffle

In order to investigate the interaction of fluid-structure with free-surface and the proposed algorithm for the par‐ticle shifting scheme concerning the free surface, numeri‐cal simulations of the collision of a water column with an elastic baffle are discussed. The computational condition of the benchmark test corresponds to the experiment by Liao et al. (2015). The schematic sketch of the problem is illustrated in Figure 21. The elastic plate is ofT=0.004 m thick with Young’s modulus and density ofE=3.5 MPa andρ=1 161.54 kg/m3,respectively.

Figure 22 presents a set of snapshots corresponding to the simulation of water dam-breaking on an elastic baffle versus the experimental results at 0.25 s＜t＜0.75 s.By comparing the present numerical results with experi‐mental data, done by Liao et al. (2015), It is observed that the developed coupled algorithm based on SPH method determines the main characteristics of the phe‐nomenon of breaking dam on an elastic baffle. As shown in Figure 22, in the early stages of the impact of the fluid on the elastic structure, an acceptable accuracy is obtained with respect to the baffle movement. With time passing, track of the blade tip shows oscillation in the experiment; however, numerical results do not show clearly these vibrations.

Figure 20 Time histories of Mechanical energy in structure and fluid sub-domain and energy balance for the whole coupled system

Hence, the horizontal displacement of the blade tip con‐tinues to the right and then moves to its stable point by fluid accumulation in two sides of the obstacle.The appli‐cation of the proposed particle regularization scheme is effective for elimination of pressure oscillations, in par‐ticular for interaction regions. The configurations of freesurface fluid flow and deformed baffle are qualitatively well reproduced compared with those of experimental re‐sults by Liao et al. (2015). Despite the local fluctuations in the pressure field, the reproduced pressure/stress fields obtained from the simulation are consistent with the nu‐merical results of Khayyer et al. (2019). Figure 23 shows

Figure 21 Schematic sketch of the test case: breaking dam on a hypo-elastic baffle based on Liao et al.(2015)experiment

Figure 22 Qualitative comparison between simulation snapshots and experimental photographs, obtained by Liao et al. (2015) in breaking-dam flow on an elastic baffle:fluid pressure field and stress component σyy in structure

Figure 23 Time histories of the horizontal displacement of the free end of baffle; present results by ISPH-TLSPH model vs.experimental results (Liao et al. 2015) and numerical results (Liao et al.2015;Khayyer et al.2019;Sun et al.2019)-Breaking dam on an elastic baffle the time histories of the displacements of the baffle’s free endUxreproduced by proposed SPH-based FSI solver along with corresponding experimental data (Liao et al. 2015) and the numerical results by FDM-FEM (Liao et al. 2015),multi-resolution MPS (Khayyer et al. 2019), single-phaseδ-WCSPH-TLSPH scheme (Sun et al. 2019) and multi-phaseδ-WCSPH-TLSPH scheme(Sun et al.2019).

It can be seen that the proposed SPH-based coupled method has provided an almost acceptable level of accu‐racy with the experimental result. To study this test case more precisely, the proposed algorithm is used for simula‐tion this problem with different geometry followed by other numerical studies. This configuration has been stud‐ied by other researchers using finite element method(FEM), simulated by Walhorn et al. (2005), particle finite element method (PFEM), done by Idelsohn et al. (2008b)and ISPH-WCSPH coupled method, simulated by Rafiee and Thiagarajan (2009).Although there are no experimen‐tal results for this configuration, however, the free surface profile and the elastic baffle deformation correspond to the physics of the problem. The initial configuration of the problem is presented in Figure 24. A column of water placed in hydrostatic equilibrium in a reservoir, and an elastic baffle with densityρb=2500 kg/m3,Young modulusE=1 MPa,and Poisson ratiov=0 in the middle of the reser‐voir is clamped from below. The geometry of the elastic baffle is widthB=1.2 cm and heightH=20/3B. The sound speed of the fluidcw=50 m/s and the initial distance be‐tween the particles Δx=2.4 mm were considered.

Figure 24 Schematic sketch of breaking dam on an hypo-elastic baffle based on numerical studies(Idel-sohn et al. 2008b;Rafiee and Thiagarajan 2009;Walhorn et al.2005)

Figure 25 Numerical simulation results of the test case:breaking-dam flow on an elastic baffle-fluid pressure field and stress component σyy in structure

Figure 25 shows the fluid pressure field and the stress componentσyyof the elastic structure. When the fluid en‐counters a baffle, because of the impact of the fluid on the lower part of the obstacle,firstly,its free end is diverted to the left.And by increasing the amount of water behind the obstacle, the baffle moves to the right. The maximum de‐viation of the obstacle to the right occurs when water jets pass through it,and the left side of the baffle is completely under the pressure of the fluid. In other words, the baffle starts to oscillate due to the pressure of the fluid, subse‐quently by fluid accumulation on both sides of the baffle,the baffle motion is gradually damped.

Figure 26 shows the horizontal displacement of the free left end of the elastic baffleUxin terms of time in order to compare with the results of Walhorn et al. (2005) (FEM method), Idelsohn et al. (2008b) (PFEM method), and Rafiee and Thiagarajan (2009) (ISPH-WCSPH method). It is observed that the curves have a relative correspondence to the first peak. Walhorn’s results predict the maximum displacement by about 18% less than the results of the present study.With passing the time, the results have been different.However,the second peak is relatively similar in time to other studies.

Figure 26 Comparison between present results by ISPH-TLSPH model vs. numerical results (Idelsohn et al. 2008b; Walhorn et al.2005; Rafiee and Thi-agarajan 2009) for time history of the horizontal displacement of the upper left corner of the baffle

The force that the fluid exerts on the front wall, after passing through the elastic baffle, depends on the Young modulus of the baffle. By increasing the Young modu‐lus, the amount of force applied to the wall decreases.The reason is that an elastic baffle with less flexibility makes more resistant to bending, therefore, more fluid momentum changes direction from horizontal to verti‐cal, meaning that there is less horizontal momentum available in the flow when it impacts the side wall. The effect of the baffle on the amount of applied force on the front wall was studied by selecting different Young’s modulus. Figure 27 shows the amount of applied forceFxon the wall in terms of time for three different Young’s moduli, which is qualitatively concordant with the physics of the problem.

Figure 27 Comparison of impact load for different Young modulus

3.2.4 Flow in sloshing tank interacting with a clamped elastic baffle

The proposed ISPH-TLSPH FSI solver is applied in the simulation of sloshing in a rolling tank with a bottom clamped elastic baffle as compared with Idelsohn et al.(2008a). The initial configuration is illustrated in Figure 28.The fluid used is sunflower oil with density and kinematic viscosity ofρoil=917 kg/m3andvoil=5×10?5m2/s, respectively.The elastic baffle with a density ofρb=1100 kg/m3and Young’s modulus ofE=6 MPa, is clamped in the middle of the tank.A sinusoidal oscillation with a maximum amplitude of 4°and a period oft=1.211 s is prescribed to the tank.

Figure 28 Schematic sketch of sloshing with an elastic baffle clamped at the bottom of a rolling tank partially filled by oil(Idelsohn et al.2008a)

Figure 29 Qualitative comparison between simulation snapshots and experimental photographs, obtained by Idelsohn et al. (2008a)corresponding to sloshing with elastic baffle clamped at bottom wall of a tank partially filled by oil.

The experimental results are presented in Idelsohn et al.(2008a), inclusive of the local displacement of the tip of the elastic baffle and four photographs taken during the experiment to illustrate the deformation of the elastic body and the free-surface flow of fluid. Figure 29 shows qualitative comparisons between simulation snapshots and experimental photographs at the same instants. In general,the reproduced free-surface profiles and elastic baffle dis‐placements well agree with the experiment. The positions of the wave trough and crest show good agreement between the simulation and the experiment. For timet=1.84 s, the free-surfaces matches with experimental result while struc‐tural deformation shows phase lag behind the expected.The minor observed discrepancies are expected to be due to an imprecise reproducibility of the working fluid’s vis‐cosity as well as considered assumptions in the mathemati‐cal model for structure. According to the presented fig‐ures, there is no un-physical gap in between fluid and structure and in specific, in the vicinity of fluid-structure interface, both the pressure and velocity divergence fields are well reproduced, thanks to the implemented fluidstructure coupling and the consistent particle shifting scheme.

Figure 30 shows the time histories of displacements of elastic baffle’s free end for comparison of simulation re‐sults with other published results (Idelsohn et al. 2008a;Khayyer et al.2018a).The displacement curve shows peri‐odic motion with respect to the sloshing tank as in the ex‐periment.The good agreement in the amplitude of the dis‐placement indicates that the coupling model can accu‐rately perform in a long term simulation of viscous fluid‐elastic structure interaction.

Figure 30 Time histories of the displacement of the elastic baffle free end; present results by ISPH-TLSPH model vs. experimental results (Idelsohn et al. 2008a), PFEM by Idelsohn et al. (Idelsohn et al. 2008a) and ISPH-SPH (Khayyer et al. 2018a)-Sloshing with an elastic baffle clamped at the bottom wall of a tank

4 Conclusion

The focus of this study is on the development, valida‐tion and application of a mesh-free computational method for simulating the flow of incompressible fluid in interact‐ing with an elastic deformable structure. Hence, a numeri‐cal coupling model based on the SPH method was devel‐oped to investigate the fluid-structure interaction. The ISPH method was used to model the fluid in the determina‐tion of pressure to create a free-divergence velocity field,while the elastic structure was analyzed by TLSPH. In the ISPH method, the regular particle distribution is effective in the reproducing of a noise-free pressure field, to give a stable and accurate simulation.Xu et al.(2009)proposed a suitable algorithm for the particle displacement and correc‐tion of hydrodynamic variables to prevent numerical insta‐bilities. However, due to incompleteness of the kernel do‐main of free-surface particles and adjacent to free-surface particles, numerical errors in their simulation lead to the lack of proper formation of the boundary of the fluid with the free surface. The proposed coupled method, without tuning a new parameter, results in a smooth and uniform pressure field in the fluid, as well as in the interface with the structure.The structure model was validated using two static and dynamic problems for the elastic structure and compared with its analytical results. Then, the numerical FSI coupled model combining with particle shifting scheme was evaluated by simulating several benchmark test cases in relation to the fluidstructure interaction. The simulation results in this study were compared with the analytical, experimental, and numerical results of other re‐searchers. The good agreement of the results presented in this study with data from the other researchers showed the ability of the proposed model to simulate the phenomenon of fluid-structure interaction. The stability and capability of the coupled FSI-SPH solver by combining an ISPH scheme for fluid flows and a Total Lagrangian SPH (TLSPH) solver for elastic structures were demonstrated by comparing the flow profile of the fluid at the interface with the structure and the free surface with experimental results.

Appendix: Pseudo-code of the proposed algorithm

Algorithm 1 Summary of the ISPH-TLSPH algorithm for each time-step n do find the neighboring particles;for each fluid particle do convect to the intermediate position r?a,by Eq.(10)calculate the intermediate velocity u?a,from Eq.(11)end for for each all particles do calculate the pressure using Eq.(13)end for for each fluid particle do update the velocity using Eq.(16)update the position by Eq.(17)evaluate δraa'based on status of particles if particle a away from surface of the structure then by Eq.(22)else by Eq.(25)end if shift the position by δraa'correct the velocity and pressure using Eqs.(29)-(30)end for for each structure particle do calculate acceleration term αf →s from Eq.(46)find the structure neighboring particles;calculate shape matrix Ka from Eq.(33)for each stage of 4th-order Runge-Kutta scheme do compute the deformation gradient from Eq.(34)calculate the displacement gradient by Eq.(39)evaluate the Lagrangian strains by Eq.(38)compute the Eulerian strains εa by Eq.(37)calculate the Cauchy stress σa using Eq.(36)compute first Piola-Kirchhoff stress by Eq.(32)evaluate the acceleration of particle from Eq.(40)update the velocity using Eq.(42)advance the position by Eq.(41)end for end for end for