A Fully Lagrangian Mesh-Free Framework for Fluid–Structure Interaction Based on WC-MPS and Hybrid TL–UL Formulations

1  Introduction

Fluid–Structure Interaction (FSI) encompasses complex multiphysics between fluid dynamics and structural (solid) dynamics and is crucial in various fields [1]. The design and assessment of fluvial and offshore structures [2,3], underwater vehicles and equipment [4,5], aircraft bodies and wings [6], and artificial heart valves [7] are a few examples of the various applications that comprise FSI. In civil engineering, FSI has a crucial impact on the behavior of structures, especially in extreme events such as strong waves, flooding, storm surges, and high-speed winds. Neglecting FSI in design or safety evaluation can lead to catastrophic failures [8].

Analytical and experimental approaches, while valuable, have often been limited in their ability to resolve highly nonlinear behavior, large deformations, and evolving interfaces [9]. Recent advances in numerical techniques and computational power have made accurate simulation of FSI more feasible [10–13].

Various numerical methods have been developed for FSI analyses. Traditional mesh-based Eulerian methods, such as the Finite Element Method (FEM) and the Finite Volume Method (FVM), have long been the basis of FSI numerical analysis [14–17]. To address free surface flows, these methods often rely on interface-capturing methods, such as the level set method [18] or the Volume of Fluid (VOF) method [19].

The Arbitrary Lagrangian–Eulerian (ALE) method [20] combines Eulerian and Lagrangian descriptions to better handle moving interfaces and deformable boundaries. It has been used for simulating rigid-body interactions and flow-induced vibrations [21,22]. However, mesh distortion under large deformations requires frequent re-meshing, thus increasing computational cost and complexity [23]. To address these challenges, the Particle Finite Element Method (PFEM) [24] was developed, offering a fully Lagrangian framework that represents free surfaces explicitly as domain boundaries. While PFEM eliminates the need for volumetric re-meshing, its frequent boundary re-meshing and contact detection steps can still be computationally demanding. To overcome the re-meshing challenge, the Immersed Boundary Method (IBM) has been developed [25,26]. IBM allows for a combination of Eulerian (fluid) and Lagrangian (structure) descriptions [27].

Numerical modeling of FSI in civil engineering faces significant challenges due to the presence of violent free-surface flows, complex boundaries, and topological changes, combined with nonlinear structural responses [28,29]. Despite their strengths, mesh-based methods struggle with large deformations and fragmentation, especially at interfaces and boundaries. These issues are commonly encountered in many FSI problems. In contrast, mesh-free Lagrangian (particle) methods, such as Smoothed Particle Hydrodynamics (SPH) [30–32] and the Moving Particle Semi-Implicit (MPS) method [33], are naturally suited to simulating large interfacial deformations and fragmentation due to their flexibility and Lagrangian formulation [34,35]. In particle methods, the computational nodes (particles) can move and adapt to any deforming boundary or interface, making them particularly efficient in simulating violent flows [36,37], significant deformations [38] and cracks [39], and mass fragmentation [40]. These distinctive attributes make these methods well-suited for FSI problems.

The MPS and SPH methods discretize the continuum using a collection of free-to-move computational nodes, called particles, without connectivity. While SPH and MPS share many algorithms, they differ in their derivative-approximation techniques. SPH relies on the derivatives of the kernel function [31], whereas the standard MPS method employs a weighted averaging of gradients between particle pairs using the kernel function itself [33,41]. However, both approaches yield closely related governing equations and stabilization mechanisms, differing mainly in the choice of density or particle number density as the primary variable. This distinction can offer specific advantages in terms of preserving local interaction symmetries in certain high-deformation regimes. Recent advances in particle methods have led to the development of reliable, entirely Lagrangian solvers for complex hydro-elastic interactions [42–44].

MPS has exceptional capabilities in modeling fluid flows characterized by high viscosity and surface tension [34,45]. In MPS, two common techniques are used for pressure calculation. The original method, known as Fully Incompressible MPS (FI-MPS), uses a semi-implicit scheme, whereas the Weakly Compressible MPS (WC-MPS) [37] uses an explicit Equation of State (EOS). MPS has also been applied to various types of solid mechanics and fracture mechanics problems, such as bending, buckling, free vibration analysis, sensitivity analysis, and topology optimization [46–48].

Errors arising from kernel approximation and discretization in particle methods can compromise stability, accuracy, and convergence [49]. These errors may lead to issues such as non-conservation of momentum, spurious pressure fluctuations, particle clustering, and tensile instability [50,51]. To mitigate these challenges, a variety of stabilization techniques have been developed. Regularization methods, including pairwise Particle Collision (PC) and Particle Shifting (PS), improve the uniformity of particle distributions, thereby enhancing accuracy [52]. Additionally, artificial diffusion, artificial viscosity, and modifications to pressure and density formulations are employed to suppress numerical oscillations and promote convergence [10]. The sensitivity of particle methods to parameter selection has motivated the use of higher-order spatial discretization schemes and adaptive particle resolution [53,54]. Moreover, multi-resolution techniques allow efficient resolution of fine-scale features by enabling local or adaptive refinement in regions with strong gradients or other critical features [55,56]. For solid mechanics applications, the MPS method also suffers from limitations such as inconsistency, tensile instability, and hourglass modes, which can affect both accuracy and numerical stability [57].

Furthermore, while this study focuses primarily on MPS and related SPH approaches, it is acknowledged that other Lagrangian mesh-free families—particularly collocation-based methods such as the Generalized Finite Difference Method (GFDM) and Radial Basis Function-Finite Difference (RBF-FD) [58,59]—offer efficient alternatives for FSI by successfully circumventing the need for complex ghost particle treatments entirely.

FSI problems exhibit complex temporal behavior due to the strong coupling between fluid and solid phases. Iterative coupling strategies are commonly used to synchronize the solvers, either through monolithic schemes, where fluid and solid equations are solved simultaneously [60], or segregated (partitioned) schemes, where the fluid and solid sub-problems are solved sequentially [60]. The segregated approach applies fluid forces to the structure and updates the fluid flow based on the deformed structure in successive time steps.

FSI coupling strategies can be broadly classified into hybrid mesh-free/mesh-based and fully Lagrangian mesh-free approaches. Hybrid methods, such as MPS–FEM [61] or SPH–FEM [62], combine mesh-free fluids with mesh-based structures to leverage each method’s strengths, but they require careful interface handling to ensure stability and accuracy [61]. Fully Lagrangian methods discretize both the fluid and solid domains into particles [63–65]. Using either MPS or SPH, these methods enable consistent FSI treatment and straightforward interface handling [66,67]. These approaches have proven highly effective across a diverse range of complex scenarios. For instance, they have been successfully applied to simulate violent free-surface impacts, such as dam-break flows interacting with deformable structures [68,69]. Furthermore, fully Lagrangian frameworks are well-suited for analyzing confined internal flows, as demonstrated by numerical investigations of sloshing dynamics within tanks equipped with elastic baffles [70–72]. Beyond these specific applications, particle-based solvers continue to be widely adopted for general hydro-elastic problems, demonstrating their inherent versatility in handling severe structural deformations coupled with fluid dynamics [73].

An important aspect of fully Lagrangian coupling is the choice of the solid formulation. Within mesh-free solid modeling, Updated Lagrangian (UL) and Total Lagrangian (TL) formulations are widely used. Antoci et al. [74] proposed a unified SPH model to simulate FSI, using a Weakly Compressible SPH (WC-SPH) solver for the fluid and an Updated Lagrangian SPH (UL-SPH) solver for the solid. While UL can capture plastic deformations effectively (e.g., Bui et al. [75]), it has some limitations, such as a lack of consistency, tensile instabilities, and zero-energy modes [61]. TL formulations, which reference all quantities to the initial configuration, are more suitable for stable, large-deformation solid modeling in FSI problems [1]. Recent advances also include multi-resolution schemes, GPU acceleration, and thin-shell models for fully Lagrangian FSI [76–79].

Building on previous developments in particle methods, this work presents a fully Lagrangian mesh-free framework for fluid–structure interaction based on the MPS methodology. The objective is not to introduce a new particle formulation, but to improve robustness and consistency by systematically integrating stabilization and correction mechanisms that are often treated separately. For the fluid phase, the weakly compressible MPS formulation is enhanced with artificial density diffusion, a modified gradient operator, and particle regularization schemes, aimed at reducing pressure oscillations and improving numerical stability.

For the solid phase, corrected kernel gradient operators and rotationally consistent particle interactions are introduced to address rank deficiency and suppress spurious zero-energy deformation modes. A hybrid Total–Updated Lagrangian formulation is adopted for hyperelastic solids, where the strain is evaluated with respect to the initial configuration while momentum balance is enforced in the current configuration using the Cauchy stress divergence. This approach preserves a stable reference for strain computation while maintaining natural compatibility with the spatial description required for fluid–structure coupling.

To ensure stable interaction between phases, a normal-flux surface integration scheme and a velocity-consistent particle grouping strategy are incorporated. In addition, a surface particle identification criterion based on the maximum angular gap between same-phase neighbors is proposed to improve boundary detection for deforming interfaces. The overall contribution lies in assembling these elements into a coherent and computationally efficient MPS-based FSI framework that enhances stability and consistency for large-deformation problems without departing from the fundamental structure of Lagrangian particle methods.

While individual components of this framework, such as the WC-MPS and hybrid TL-UL formulations, are established in the literature, the primary novelty of this work lies in their consistent and targeted integration to address highly dynamic FSI challenges. Specifically, this framework introduces: (1) a fully consistent particle-based coupling strategy using identical interaction operators across all phases. This intrinsically mitigates spurious pressure fluctuations at the fluid-solid interface, bypassing the need for ad-hoc normal flux corrections, penalty terms, or artificial damping; (2) a hybrid TL-UL formulation tailored for FSI that maintains direct compatibility with fluid forces under large deformations; (3) an improved interface detection approach utilizing particle number density and angular spacing; and (4) the iterative corrected Particle Shifting (icPS) scheme to preserve stability and avoid order-dependent bias.

2  Methodology

2.1 Governing Equations

The mathematical model that describes the internal motion of a continuum in isothermal conditions under the effect of external forces, such as gravity, is governed by the conservation of mass (continuity) and momentum [80].

DρDt=−ρ∇⋅u,(1)

DuDt=1ρ∇⋅σ+f,(2)

where t represents time, ρ is the density of the material, and u denotes the velocity vector field, f represents the external accelerations applied to the material (e.g., gravity g), and σ denotes the stress tensor, describing internal response stresses within the material.

2.1.1 Fluid Dynamic Equations

The stress tensor for a Newtonian incompressible fluid is described as:

σF=−pI+τ=−pI+2μFε˙,(3)

where p is pressure, τ is the viscous stress tensor. μF is the fluid’s dynamic viscosity, ε˙ is the rate of strain tensor, I is the identity tensor, subscript “F” denotes fluid phase. Therefore, the governing equations describing incompressible fluid flow in a Lagrangian framework can be written as:

DρFDt=−ρF(∇⋅u),(4)

DuDt=−1ρF∇p+μFρF∇2u+g+fS→F,(5)

where fS→F is the solid force on fluid. To calculate pressure, the weakly-compressible (WC) MPS method is used herein [34,37], where an equation of state (EOS) gives the pressure as an explicit function of density. In the Lagrangian framework of fluid analysis, an equation of motion updates for translational movement, defining how position evolves in response to velocity:

DrDt=u,(6)

where r denotes the position vector.

2.1.2 Solid Dynamic Equations

Here, we assume the solid material to be homogeneous, isotropic, and purely elastic. The solid stress tensor, σS is given by:

σS=Ktr(ε)I+2G(ε−13tr(ε)I)=λStr(ε)I+2μSε,(7)

where subscript “S” denotes the solid phase, ε is the solid strain tensor, K and G denote the bulk and shear modulus of solid material, respectively. λS and μS are the Lame’s constants, describing the mechanical properties of the solid material, calculated as a function of Young’s modulus (Elasticity), E, and the Poisson’s ratio, νS, as [81]:

λS=EνS(1+νS)(1−2νS)andμS=E2(1+νS).(8)

The solid model is based on Hooke’s generalized law for linear elasticity, as described by Hwang et al. [63]. So, the governing equations are reshaped into the following formula employing Eqs. (2) and (7):

DuDt=1ρS∇⋅(λStr(ε)I+2μSε)+g+fF→S,(9)

where ρS represents the density of the solid material. fF→S denotes the fluid forces acting on the solid. For clarity in the present hybrid TL–UL solid implementation used later in Section 2.2.3, strain measures are evaluated with respect to the initial configuration, while momentum balance is enforced in the current configuration using the divergence of the Cauchy stress.

In solid mechanics, the conservation of angular momentum is integrated into the solution, as described by:

DDt(Imω)=DDt(r^×(mu))=r^×DDt(mu)=r^×F^,(10)

where ω is the angular velocity and denotes the moment of inertia of the rotating mass, m. The term mu represents the linear momentum, with u being the linear velocity vector. The symbol F^ represents the total force vector responsible for driving the object’s rotation, and r^ corresponds to the moment arm of the momentum and represents the position vector from the center of rotation to the application point of F^. The equations of motion describe how the position and orientation of a solid body change over time in response to its respective velocities. The translational and rotational motion of the solid body is given by:

DrDt=uandDθDt=ω,(11)

respectively, where r is the body’s position vector and θ is the rotational degree of freedom or orientation angle.

2.2 Numerical Method

2.2.1 MPS Fundamentals

Here, we use MPS mesh-free Lagrangian numerical methods for approximation and solution of the governing equation. In the MPS method, the continuum is represented by a set of particles (computational nodes) with no connectivity, which can move in response to the material’s motion. The particles are distributed uniformly, with an initial distance (or size) of l∘. MPS approximation of quantities and spatial derivatives is based on the interaction of each particle, i, and its neighbor particles, j≠i, within an effective radius re, and is evaluated with a kernel (weight) function. The same particle set and kernel concept are used for both phases; the phase-specific discretizations are detailed in Sections 2.2.2 and 2.2.3, and the coupling terms are defined in Section 2.2.4.

The kernel assigns weights based on the Euclidean distance between i and j particles, |rij|=|rj−ri| where r is the position vector. The most commonly used kernel function in MPS-based simulations is the rational function given by [33]:

W(r,re)={rer−1,0≤r<re0,re≤r(12)

This study uses the recommended [33] effective radius of 3.0l∘ for the Laplacian operator and other calculations. In MPS context, the particle number density at i, symbolized as ⟨n⟩i, is defined as a dimensionless measure of the material density ρi and is calculated as:

⟨n⟩i=∑j≠iNiWij ,Wij=W(|rij|,re).(13)

Using the notation “⟨ ⟩i” denotes the MPS kernel approximation. The basic MPS approximation of the gradient of a scalar function φ and the divergence and Laplacian of a vector field ϕ can be defined as:

⟨∇φ⟩i=Dn0∑j≠i(φj−φi)eij|rij|Wij(14)

⟨∇⋅ϕ⟩i=Dn0∑j≠i(ϕj−ϕi)⋅eij|rij|Wij(15)

⟨∇2ϕ⟩i=2Dn0λ∘∑j≠i(ϕj−ϕi)Wij(16)

where eij=rij/|rij| defines the normalized or unit vector associated with rij, D is the number of spatial dimensions, and n0 is the initial particle number density, which is used as a normalization factor. The normalization parameter λ∘ ensures that the increase in variance is equal to that of the analytical solution and is given by:

λ∘=∑j≠i|rij|2Wij∑j≠iWij(17)

The standard operators defined above are used for the fluid phase unless otherwise specified. To enhance stability and accuracy, modified formulations are employed in certain parts of the model, as specified in the corresponding sections.

2.2.2 MPS Discretization for Fluid Phase

The fluid solver is based on an enhanced weakly-compressible MPS method [50]. The governing equations of a fluid particle i can be written as:

{DniDt=ni(−⟨∇⋅u⟩i+Dmi)DuiDt=−1ρF⟨∇p⟩i+μFρF⟨∇2u⟩i+g+⟨fS→F⟩iDriDt=uipi=f(ni)(18)

where Dmi is a density diffusion term proposed by Jandaghian and Shakibaeinia [50], in the so-called DMPS method, to compensate for numerical diffusion and reduce pressure oscillation in the MPS method. The diffusion term is defined as:

Dmi=δMPSΔtF(cF)2n0⟨∇2n⟩i(19)

where δMPS is the diffusion coefficient (δMPS=0.50 is used for all the simulations in this study), ΔtF is the time step of fluid solution, and cF is a numerical sound speed (artificial parameter).

The Laplacian of particle number density is approximated as:

⟨∇2n⟩i=2Dn0∑j≠i(nj−ni|rij|−⟨∇n⟩i+⟨∇n⟩j2⋅eij)1|rij|Wij(20)

where

⟨∇n⟩i=Dn0∑j≠inj−ni|rij|(Cieij)Wij(21)

Ci=(Dn0∑j≠ieij⊗eijWij)−1(22)

Here, the divergence operator for the velocity field can be defined as [50]:

⟨∇⋅u⟩i=Dn0∑j≠i(ninj)(uj−ui)⋅eij|rij|Wij(23)

The gradient of the pressure term is approximated by [50]:

⟨∇p⟩i=Dn0∑j≠i((ninj)pj+(ninj)pi)eij|rij|Wij.(24)

Also, the Laplacian of the velocity vector is given by:

⟨∇2u⟩i=2Dn0λ∘∑j≠i(uj−ui)Wij.(25)

The coupling force, fS→F, will be defined later in the Section 2.2.4. The weakly-compressible MPS method explicitly calculates pressure as a function of particle number density (pi=f(ni)) using Tait’s equation of state [37]:

pi=ρF(cF)2γ[(⟨n⟩in0)γ−1],γ=7.0(26)

The numerical sound speed cF in the WC–MPS method controls the compressibility of the fluid. It is chosen to be sufficiently large to limit density variations to an acceptable level (typically below 1%), while remaining small enough to avoid excessive reduction of the time step required for numerical stability.

2.2.3 MPS Discretization of Solid Phase

The solid phase is modeled assuming co-rotational, linearly elastic behavior under any input (often small) strain conditions, for simplicity. The strain tensor is evaluated using particle interactions defined with respect to the initial configuration. The methodology for solid phase discretization employed in this section builds upon the particle-based elastic dynamics frameworks established by Song et al. [81] and Hwang et al. [63].

In the present implementation, a hybrid Total–Updated Lagrangian (TL–UL) strategy is adopted. Strain measures are evaluated in the reference configuration, while momentum balance is enforced in the current configuration using the Cauchy stress tensor. This formulation allows deformation to be measured relative to a fixed reference configuration while maintaining compatibility with the spatial description used in the fluid solver. As a result, pressure forces from the fluid can be directly transferred to the structure without requiring transformations between reference and spatial stress measures. The approach also facilitates the use of the standard MPS interaction operators, which are formulated in the current configuration. The choice of the hybrid TL–UL formulation is made to facilitate implementation within the particle-based framework and does not restrict the generality of the proposed method, as the formulation could equivalently be expressed using alternative stress measures.

The strain tensor is computed as:

εi=12(∇di+(∇di)T)(27)

The gradient operator used in strain evaluation employs the initial inter-particle configuration to improve numerical stability and reduce accumulated distortion under large deformation.

Fig. 1 illustrates how the current positions of particles j and i relative to each other are evaluated. The vectors rij0 and rij represent the position vectors of particle i to particle j in the initial state (rij0=rj0−ri0) and the current state (rij=rj−ri), respectively. The relative rotation, achieved by fixing i and a pure angular motion of j, is calculated as θij=0.5[θi+θj]. Consequently, the vector Rijrij0 signifies the virtual position of particle j if particle j only rotates by θij. The rotated position vector of the initial position vector is calculated as Riri0, where, as an example, for a 2D rotation, Ri is given by:

Ri=R(θi)=[cos⁡θi−sin⁡θisin⁡θicos⁡θi](28)

Figure 1: The model for relative motion between two structure particles i and j.

The vector of pure translational displacement due to strain, dij, is obtained by subtracting Rijrij0 from rij:

dij=rij−Rijrij0(29)

The displacement vector, dij, is divided into normal (subscript n) and tangential (subscript t) components as:

dij,n=(dij⋅eij)eijanddij,t=dij−dij,n.(30)

The normal strain, εij,n, and the shear strain, εij,t, between particles i and j are

εij,n=dij,n|rij0|andεij,t=dij,t|rij0|.(31)

The stress tensor is computed from Hooke’s law as:

σi=λStr(εi)I+2μSεi(32)

The solid momentum equation is discretized as:

DuiDt=1ρS∇⋅σi+g+fF→S(33)

Note that the divergence operator applied in Eq. (33) acts on the Cauchy stress tensor in the current configuration. Therefore, while strain is evaluated with respect to the reference configuration, equilibrium is enforced in the spatial configuration. This hybrid TL–UL treatment differs from a classical pure Total Lagrangian formulation, which would require divergence of the first Piola–Kirchhoff stress tensor with all operators expressed in the reference configuration.

The isotropic pressure corresponds to the trace of strain and is discretized using the MPS divergence operator:

⟨tr(ε)⟩i=D⟨n0⟩i∑j≠iεij⋅rij|rij0|2Wij0 ,Wij0=W(|rij0,re|),(34)

where using the local ⟨n0⟩i instead of a global n0 allows more stable approximation [63]. The latter simplifies into:

⟨tr(ε)⟩i=D⟨n0⟩i∑j≠i|dij,n||rij||rij0|3Wij0.(35)

Translational acceleration of particles is obtained from shear stress and isotropic pressure. Isotropic pressure is obtained at the particle locations, but stresses are obtained between the particles. Translation of particles is described by:

⟨1ρS∇⋅(λStr(ε)I+2μSε)⟩i=λSρS⟨∇⋅tr(ε)I⟩i+2μSρS⟨∇⋅ε⟩i(36)

where, tr(ε)I=0.5[⟨tr(ε)⟩i+⟨tr(ε)⟩j]I. The first and second terms on the right-hand side of Eq. (36) are given by:

λSρS⟨∇⋅tr(ε)I⟩i=12λSρSD⟨n0⟩i∑j≠i([⟨tr(ε)⟩i+⟨tr(ε)⟩j]I)⋅rij|rij0|2Wij0(37)

2μSρS⟨∇⋅ε⟩i=2μSρSD⟨n0⟩i∑j≠i(rij−Rijrij0)⋅rij|rij0|2Wij0(38)

MPS discretization Eq. (10) can be written as:

IiDωiDt=⟨r^×F^⟩i=−μS1ρSD⟨n0⟩i∑j≠irij×(Rijrij0)1|rij0|2Wij0(39)

where Ii=mil∘2/6 represents the moment of inertia for a particle i with mass mi in the case of a square (or cube) shaped particle. F^ is the force applied by a neighbor particle on particle i. Note that the interaction force driving the rotation is proportional to the tangential displacement (rij−Rijrij0). Because the cross product of the vector rij with itself is zero, the moment expression simplifies strictly to −rij×(Rijrij0). The cross product only considers the contribution of shear stress (as described in Eq. (7)) to the rotational momentum, while any contributions from forces acting in the i-to-j direction (i.e., eij) are disregarded. Although the present framework is implemented in 2D—where rotational motion reduces to a scalar about the Z-axis—consistent vector notation (ω and θ) is strictly maintained throughout the formulation to facilitate future extensions to 3D fully Lagrangian applications. In addition, the momentum rotating the pair is divided equally between the two interaction particles.

Suppressing spurious zero-energy deformation modes and hourglass-type instabilities is primarily achieved through the use of the Laplacian model for the estimation of the shear term, as highlighted by Wu et al. [82] in the TL-SPH formulation. Furthermore, the corrective matrix (Eq. (22)) is specifically utilized in the approximation of the gradient of the particle number density (Eq. (21)) to ensure uniform consistency.

Overall, the governing equations for a structure particle i can be summarized using MPS into:

{DuiDt=λSρS⟨∇⋅tr(ε)I⟩i+2μSρS⟨∇⋅ε⟩i+g+⟨fF→S⟩iDωiDt=1Ii⟨r^×F^⟩iDriDt=uiDθiDt=ωi(40)

2.2.4 Fluid–Structure Coupling

The interactions in FSI are modeled through a fully Lagrangian particle-based coupling strategy, where both fluid and solid phases are discretized into particles and interact through interface forces.

Interface forces are computed using the same MPS interaction operators defined previously, ensuring consistency between intra-phase and inter-phase interactions. No additional mesh-based interpolation or projection procedures are needed.

At the interface, the continuity of velocity and traction is enforced. The interaction force between fluid and solid particles is computed solely from the pressure (shear contributions being neglected). The force exerted by the fluid on the solid, fF→S, represents the pressure force originating from the fluid domain, which is consistently transferred to the solid phase. The force exerted by the solid on the fluid, fS→F, is induced by boundary motion and appears in the fluid momentum equation. In physical terms, the interaction force fS→F generates a localized pressure gradient in the neighbor fluid, resulting in positive or negative pressure variations as the boundary moves toward or away from the fluid phase, respectively (See Fig. 2). These interaction forces are strictly equal and opposite (fS→F=−fF→S), ensuring action–reaction consistency between phases.

Figure 2: Presentation of the fluid-structure interactions (indicated by dashed lines) occurring at the interface (IF). The tinted particles serve as the query point, and the big solid circles define the effective radius: (a) computation method for the coupling force of the structure surface to the fluid; (b) calculation process of the force applied by fluid pressure on the structure interface.

This symmetric force exchange guarantees local conservation of linear momentum at the interface. Angular momentum conservation is preserved through consistent pairwise interaction forces aligned along particle connection vectors.

During each fluid time step (ΔtF), only the positions of fluid particles are updated, while non-fluid particles are kept fixed, act as solid boundary particles, and exert forces on the fluid. These forces are incorporated into the fluid equations as a coupling force term, denoted by fS→F, in Eq. (18). This term represents the influence of the structure on the fluid. Specifically, the coupling force fS→F is defined based on the net acceleration of the structure interface (See Fig. 2a).

fS→F=ddtuS, Interface(41)

Considering only solid particles on the interface (IF) gives us

⟨fS→F⟩i=∑j∈IF(DujDt)Wij∑j∈IFWij(42)

In the solid solver, during the solid time-step (ΔtS), only solid (structure) particles will be updated, and otherwise, boundary, support, and fluid particles are kept fixed. The boundary and support particles define the boundary condition. The fluid particles primarily exert pressure on the solid surface, which initiates the motion of the structure interface. Therefore, the fluid effect is considered as a force term from the fluid to the structure, fF→S in the governing equations (See Eq. (40)). The formula to calculate the external force on the structure surface is given by integration of fluid pressure on the solid surface as [63]:

FF→S=∫∫InterfacepF dA(43)

where A denotes the solid interface area (i.e., the interface, IF), and pF represents the fluid pressure. Using MPS method, the acceleration force, fS→F, at a surface solid particle, i, is sourced from pressure of nearby fluid particles, j, and can be given by (See Fig. 2b) [83]:

⟨fF→S⟩i=1ρiVi∑j∈FpjW(|rj−ri′|)∑j∈FW(|rj−ri′|)Ai(44)

Here, Ai=l0D−1 and Vi=l0D are the surface area and volume of the target particle i, respectively. The point ri′ is expected as the pressure’s application point on the interface for particle i. This is to avoid invalid kernel values and extra cut-off of the pressure force by the effective radius, which leads to inconsistent separation of structure and fluid at the interface. Here, ri′ is calculated as:

ri′=ri+0.5l∘si(45)

where si is the surface normal vector of the interface towards the fluid at particle i.

2.2.5 Temporal Stability

Accurate time step calculation is crucial for stable and efficient simulation of both fluid and elastic materials. In this study, we calculate the fluid time step, ΔtF, based on stability conditions for advection and diffusion processes [73,84]. This approach ensures compliance with the Courant–Friedrichs–Lewy (CFL) condition:

ΔtF=CCFLmin{l∘vmax,F+cF,l∘2(μF/ρF)} .(46)

Here, vmax,F=max(|u|) represents the maximum fluid velocity, ensuring stability by controlling the influence of convective and diffusive terms [85]. We selected CCFL=0.1 for the fluid simulations, allowing optimal stability without significantly compromising computational efficiency.

The same CFL condition applies to the time step of the solid phase, ΔtS, which is essential for the stability of simulations involving elastic boundaries [81]. Specifically, we have

ΔtS=CCFLmin{l∘gE/ρS+cS,l∘2(G/ρS)} .(47)

In this case, the term gE/ρS represents the speed of an elastic wave in the solid, dependent on Young’s modulus, E, and the density of the material, ρS [81]. For the elastic material, a more conservative CCFL value of 0.05 is recommended to effectively capture the higher stiffness and slower wave propagation in solids. Similarly to the concept of diffusion in fluids, we consider the analogy of the elastic body as a viscous fluid with an effective viscosity equal to the shear modulus G [84]. This allows the diffusion-based time step term to be incorporated, akin to a fluid’s treatment in the second term of Eq. (46), enhancing the stability of our fully Lagrangian framework.

In addition, to achieve a synchronized simulation and a more stable coupling between the two separate solvers with varying time steps, the ratio of ΔtF to ΔtS is chosen as a natural number.

ΔtFΔtS∈N.(48)

2.2.6 Boundary Treatments

To accurately detect interface particles within highly dynamic FSI scenarios involving deformed boundaries and free-surface flows, an additional particle registration limit is applied. Interface particles are identified based on their neighborhood phase configuration using two distinct criteria. These criteria are based on 1) the particle number density condition (similar to Antoci et al. [74]) and 2) the angular spacing between radial lines of neighbor particles, as illustrated in Fig. 3. Specifically, a particle is classified as an interface particle when its effective neighborhood contains particles belonging to both phases, based not only on the interaction radius but also on the angular portion of the circular neighborhood they occupied. This criterion allows robust identification of evolving fluid-air (gap), solid-air, and fluid–solid (elastic/fixed-boundary) boundaries without explicit geometric reconstruction.

Figure 3: Illustration of the angular opening between radial lines from a fluid particle to neighbors. δamax is the maximum radial opening among others.

A fluid particle is classified as a free surface particle when its number density (n) falls below β1n0 and the maximum opening angle, represented by δamax, exceeds β2(360) degrees within the effective radius. In our simulations, we have set β1 to 0.95 and β2 to 0.25.

Accurate kernel interpolation requires full kernel support, which is compromised near interfaces. For the fluid solver, the simulation setup comprises one layer of fixed particles serving as a container wall, thus imposing rigid wall boundary conditions (See Fig. 4). The wall particles have their pressure values determined by solving Eq. (26) concurrently with the fluid particles. To augment the kernel domain for wall particles, three additional layers of particles are introduced beyond the wall particles. These supplementary particles, known as ghost particles, help in calculating numerical entities of wall and near-wall fluid particles (incomplete kernel domain, truncated-kernel error) within the fluid equations. Values are not computed for ghost particles; instead, interpolation techniques are employed for this purpose. To ensure consistency in the fluid solver, the outermost layer of structure particles is treated as wall particles, while those residing within the structure are designated as ghost particles, facilitating their interaction with the fluid solver. The pressure field of fluid is considered for the structure particles by imposing a Neumann boundary condition for the interface.

Figure 4: Definition of particle types for the fluid solver.

In this context, pseudo-wall and pseudo-ghost particles (See Fig. 4) refer to standard solid particles located at the interface that temporarily carry a dual set of physical values. This dual representation allows them to actively participate in fluid kernel approximations without requiring the generation or regeneration of independent, auxiliary fluid particle sets at every time step.

In case of a large deformation of the boundary particles (large deformation of the structure) or complex geometries (narrow edges or free fluid surface), we adapted a ghost particle boundary method proposed by Long et al. [86] and Zhang et al. [87], where the interceptive area of the kernel support domain is divided into subareas based on the boundary segment, and these sections are later used to produce the corresponding ghost particles. Using these ghost particles, the effect of the truncated kernel domain due to geometric boundaries can be calculated. This ghost particle method can treat complex boundaries and deforming interfaces of fluid and the structure.

For the solid solver component of our model, all fluid particles are excluded from calculations, and their influence, in the form of forces, is imparted to the system. In scenarios involving boundary collisions, due to extreme displacements or interactions with another solid object, a straightforward linear elastic collision technique, proposed by [34], is applied. This technique involves the exchange of linear momentum between sets of colliding particles. The structure particles in the exterior layer with more than two fluid particles within their effective radius are considered in contact with the fluid or interface particles (IF). It should be noted that, unlike the fluid phase, the solid phase discretization does not require the deployment of ghost particles. Boundary conditions, such as clamped supports, are strictly enforced by constraining the kinematic properties (i.e., fixing displacement and velocity to zero) of the designated boundary solid particles.

To improve stability and prevent numerical penetration at the interface, a normal-flux surface integration approach is incorporated. The normal-flux evaluation ensures that interface forces are primarily transferred along the local normal direction of the boundary, reducing spurious tangential oscillations and improving traction continuity [88].

A velocity-consistent particle grouping strategy is also employed to improve interaction consistency across phases. This grouping strategy ensures that interacting fluid–solid particle pairs share compatible velocity updates within each time step, reducing artificial slip and enhancing interface stability under large deformation [87].

No additional penalty stiffness or artificial interface damping is introduced; stability is achieved through consistent particle interaction, corrected gradient operators in the solid phase, and density diffusion in the fluid phase.

2.2.7 Particle Regularization for Fluid Particles

To maintain a uniform particle distribution—essential for the accuracy and stability of MPS kernel approximations—this study employs a combined regularization framework. We utilize the Dynamic pairwise Particle Collision (DPC) technique [50] to alleviate particle clustering, specifically in violent free-surface and splash regions. Building on the work of Jandaghian et al. [52], who integrated DPC with corrected Particle Shifting (cPS) to stabilize Lagrangian simulations of highly dynamic flows (e.g., dam breaks and sloshing), we propose a newly developed iterative corrected Particle Shifting (icPS) method. The combined framework is hereafter referred to as DPC+icPS.

Standard Particle Shifting Techniques (PST) [52,89,90] typically calculate and apply a single displacement vector once per fluid time step ΔtF. While computationally inexpensive, a single iteration often fails to fully resolve clusters in high-shear or high-density zones, potentially leading to boundary penetration or unphysical particle overlapping. To address this, the proposed icPS acts as a local spatial relaxation procedure. The shifting displacement is calculated and applied k=5 times sequentially within a single fluid time step. For each iteration m∈{1,…,k}, a partial shifting vector δrim=(1/k)δri,cPS is recalculated based on the updated particle positions from the previous sub-step:

rit,m=rit,(m−1)+δrim(49)

A critical feature of the icPS algorithm is that the new positions for all particles are updated simultaneously (altogether) at the end of each sub-step, rather than sequentially one after another. This overall update strategy ensures that the algorithm remains strictly independent of particle indexing order, preserving the symmetry of particle interactions and preventing the accumulation of directional bias.

As illustrated in Fig. 5, this multi-pass relaxation ensures the particle distribution progressively reaches a stable, uniform equilibrium. Unlike standard one-off shifting methods (Fig. 5b), which can distort interfaces under extreme deformation, the proposed DPC+icPS effectively preserves the integrity of the body shape without promoting numerical artifacts, such as unphysical penetration into boundaries or overlapping with adjacent particles (Fig. 5a). Ultimately, this approach provides a computationally efficient alternative to full Arbitrary Lagrangian-Eulerian (ALE) formulations. The icPS also significantly improves the smoothness of the pressure field and overall stability under large deformation compared to standard one-off shifting methods.

Figure 5: Schematic representation of (a) the iterative corrected Particle Shifting (icPS) vs. (b) prior standard shifting techniques. The multi-step relaxation (k=5 in this study) maintains the body shape and prevents particle clustering or boundary penetration. The orange particle represents the initial position, the green particle indicates the intermediary update stages, and the blue particle shows the final relaxed position.

2.2.8 Solution Algorithm

The coupling algorithm proceeds sequentially within each time step. First, fluid particles are updated using the WC-MPS formulation. Second, interaction forces are evaluated at the interface. Third, the solid phase is advanced using the hybrid TL–UL formulation.

This segregated coupling strategy ensures numerical stability while avoiding the complexity of a fully monolithic solution scheme. The chosen time integration remains explicit for both phases, with stability governed by CFL-type constraints associated with the artificial speed of sound and elastic wave speed.

Fig. 6 shows the coupling process and the fluid and solid solvers within the present model throughout the simulation time. Fluid and solid solvers operate sequentially and independently with distinct time intervals, namely, ΔtF and ΔtS for fluid and solid phases, respectively. Typically, ΔtF>ΔtS as the wave speed in the fluid is smaller than that of the solid (cF<cS). For the fluid, a symplectic time integration scheme is employed as described in Jandaghian and Shakibaeinia [50].

Figure 6: Calculation order coupling scheme of fluid and solid solvers.

The solution algorithm is implemented as a GPU-accelerated parallel code developed in C++ using NVIDIA’s Compute Unified Device Architecture (CUDA) framework. To reduce computational time, an improved neighbor-search method is used.

3  Numerical Validations

The developed numerical model is validated through four distinct test cases to systematically evaluate the framework’s capabilities. First, the individual solvers are assessed independently: the solid solver is validated in Case 1 (dynamic response of a cantilever beam) and the fluid solver in Case 2 (dam-break flow). Next, the fully coupled fluid-solid interaction mechanism is quantitatively validated in Case 3 using the standard Turek-Hron FSI3 benchmark. Finally, the framework is applied to a complex engineering scenario in Case 4, simulating a dam-break flow interacting with an elastic gate.

3.1 Case 1: Initially Disturbed Cantilever Elastic Beam in the Absence of External Forces

To evaluate the accuracy of the proposed algorithm in capturing large elastic deformations, a dynamic test of a cantilever beam is considered. As shown in Fig. 7, a thin plate is clamped at one end and free at the other. An initial vertical velocity is applied, causing harmonic vibration. This widely used benchmark has been used in past studies for model validations (e.g., [63,74,84,91]). The initial vertical velocity distribution along a beam of length L is analytically defined as:

uy0(x)=0.01cSf(x)f(L)(50)

where 0.01cS corresponds to the prescribed initial velocity at the free end of the beam (x=L), and the velocity decreases along the beam according to the modal shape function f(x). Also,

f(x)=(cos⁡(kL)+cosh⁡(kL))(cosh⁡(kLxL)−cos⁡(kLxL))+(sin⁡(kL)−sinh⁡(kL))(sinh⁡(kLxL)−sin⁡(kLxL))(51)

Figure 7: Initial configuration of case 1; the red dot (D1) shows the reference point to evaluate the beam’s deflection.

In beam theory, k (the wavenumber or modal parameter) and L (the length) are separate variables. Here, the kL product is equal to 1.875 for the fundamental mode. The parameter cS=K/ρS denotes the speed of sound in the structure. The actual physical longitudinal wave speed in a thin 1D beam is E/ρ. For this test case, the bulk modulus K=3.25 MPa and the solid density ρS=1000 kg/m3. Also, shear modulus G is set to 0.715 MPa, which results in νS=0.395 and E=2.01 MPa.

Fig. 8 shows the time history of the displacement at the free end of the plate (marked by a red dot in Fig. 7). The plate undergoes periodic motion; therefore, the simulated dimensionless period and amplitude are compared against both the analytical solution and previously published numerical results. These comparisons are summarized in Table 1, which shows that the present results closely match the analytical solution. The current study’s dimensionless period is 73.30, which is slightly higher than the Analytic solution (72.39). This 1.2% deviation is likely due to the discretization of the clamped boundary condition. Although minor deviations are observed in the numerical results, the agreement of our results remains strong overall, as shown in the table.

Figure 8: Time history of the deflection at D1 in case 1 using the analytical solution and the current numerical (Num.) model.

For error quantification and convergence study, different metrics are evaluated based on the specific test case. In Case 1, the displacement at the free end of the beam is compared against the analytical solution of linear Euler–Bernoulli beam theory. The relative error is computed as:

Error=∑k=1N(δkNum.−δkRef.)2∑k=1N(δkRef.)2(52)

where δkNum. and δkRef. denote the numerical and reference displacements, respectively, and N is the total number of recorded time steps.

Simulations were performed using systematically halved particle resolutions (l∘, l∘/2, and l∘/4) to rigorously assess convergence behavior. As shown in Table 2 and the log-log error plot in Fig. 9, the error reduces monotonically. The slope of the error curves approaches the 𝒪(h2) reference line, confirming that the spatial discretization converges at a nearly second-order rate.

Figure 9: Log-log plot demonstrating the spatial convergence rate of the proposed framework. The dashed line represents second-order convergence.

3.2 Case 2: Dam Break Flow

This case is a dam-break problem created by the sudden release of a water column confined by a dam (gate) in a tank. It is a widely used benchmark for evaluating particle-based simulation of highly dynamic free-surface flow problems [37]. In this test case, after gate removal, the water column collapses, and the water flows across the horizontal bed, subsequently impacting a downstream wall, leading to the formation of a returning wave (jet). The primary challenge with simulating this case is producing a stable solution with a noise-free pressure field and accurately representing the evolution of the highly deformed free surface. Therefore, the experiment conducted by Lobovský et al. [92] is selected for this validation. The simulated results are compared with experimental data and previous numerical findings from [50,93,94].

The water column has a height of H2=0.60 m, length of L2=2H2, and the tank has a length of LTank=5.366H2 as shown in Fig. 10. The fluid is water with a density ρF=1000 kg/m3 and a dynamic viscosity μF=0.001 kg/(m ⋅ s). The initial spacing between particles (particle size), l∘, is set to 0.005 m, corresponding to a total of approximately 30,000 fluid particles. The validation process includes assessing the propagation of the main wavefront and monitoring pressure at the impact wall (gauge S1 in Fig. 10).

Figure 10: Initial configuration of case 2; the red dot (S1) shows the reference point to compute the pressure, which is the pressure sensor in the original experiment done by [92].

Fig. 11 illustrates the time evolution of the simulated dam-break wave. A key highlight is the smoothness of the pressure field, which demonstrates the capability of the proposed method in addressing one of the known challenges in particle methods, i.e., pressure noise and instability in dynamic free-surface flows.

Figure 11: Snapshots of the dam break flow of case 2 at dimensionless times (T=tgH2) and showing dimensionless pressure contours (p/ρgH2).

A snapshot of the wave impact near the wall is presented in Fig. 12, where the simulated free-surface profile is compared against the results from a Mixed-Eulerian-Lagrangian Boundary Element Method (BEM-MEL) [93]. The comparison shows excellent agreement in both shape and position of the wave front, validating the physical accuracy of the model in reproducing highly nonlinear free-surface dynamics.

Figure 12: The magnification of the returning water jet at T=6.0 in case 2. The solid purple line represents the digitized free surface computed by the BEM-MEL solver [93].

Further validation is presented in Fig. 13, which shows the pressure time history at gauge S1. These results are compared with both experimental measurements by Lobovský et al. [92] and the numerical results of Jandaghian et al. [85], showing an overall good agreement. Finally, Fig. 14 shows the progression of the water front. The simulated results are benchmarked against analytical solutions (shallow flow theory), experimental observations, and existing numerical data. The propagation trends and impact characteristics align well with all references, reinforcing the credibility of the proposed approach. Minor deviations, particularly when compared with experimental results, may be attributed to unmodeled physical factors—such as the roughness of the flume bed, which can influence wave speed at the leading edge, effects that are difficult to resolve without detailed modeling.

Figure 13: Pressure time history at S1 gauge showing a near-noise-free fluid model (Exp. denotes experimental data from [92], and Num. denotes the numerical method. Numerical reference data from [85]).

Figure 14: Water wavefront progression vs. time (Exp. denotes experimental data from [92], and Num. denotes the numerical method. Numerical reference data from [84]).

Overall, the results confirm that the developed model is capable of accurately capturing the key fluid dynamics, including pressure evolution, wavefront progression, and nonlinear free-surface deformation. For the convergence study, the time evolution of the water wavefront position is compared with the reference measurements. The relative error in the front location is computed by normalizing the absolute deviation against the initial water column length (L2 = 2H2). Results (Table 2) confirm the model convergence (with an order of around 1.8).

3.3 Case 3: Flow around a Flexible Flag

To quantitatively validate the coupling mechanism of the model, the standard Turek-Hron FSI3 benchmark [95] is simulated. This benchmark is widely adopted in the literature to evaluate the accuracy of FSI solvers in capturing self-sustained, large-deformation oscillations. The problem consists of a laminar channel flow interacting with an elastic flag attached to a fixed rigid cylinder.

As depicted in Fig. 15, the computational domain is a rectangular channel with dimensions 2.5m×0.41m. A rigid cylinder of radius r=0.05m is centered at (0.2m,0.2m). An elastic flag of length L=0.35m and thickness h=0.02m is clamped to the downstream side of the cylinder. The fluid is characterized by a density ρF=1000kg/m3 and kinematic viscosity νF=0.001m2/s. A parabolic velocity profile is prescribed at the inlet with a mean velocity of u¯x=2.0m/s, resulting in a Reynolds number of Re=200. The solid flag has a density ρS=1000kg/m3, Poisson’s ratio νS=0.4, and shear modulus G=2.0MPa.

Figure 15: Initial geometric configuration for the Turek-Hron FSI benchmark [95].

Under these conditions, the flag undergoes periodic, self-sustained flapping induced by the vortex shedding from the upstream cylinder. To evaluate the solver’s performance, the displacement of the flag’s trailing edge (tip) is tracked over time. Table 3 compares the amplitude and frequency of the tip displacement in both the longitudinal (X) and transverse (Y) directions against the reference data provided by Turek & Hron [95].

For the spatial convergence analysis of the coupled solvers (l∘=0.0025 m), the amplitude of the transverse tip displacement (Y-displacement) serves as the primary convergence metric, as it is highly sensitive to the coupled hydro-elastic dynamics and interface resolution. The wave speed in the solid is roughly cS≈100 m/s. Using a conservative CFL condition, the solid time step should be around ΔtS=2.0×10−6 s. Because the max fluid velocity is 3.0 m/s and the numerical sound speed cF is usually 30 m/s, the fluid time step can be safely set to ΔtF=1.0×10−5 s. This perfectly accommodates the requirement that ΔtF/ΔtS∈N (in this case, exactly 5).

Fig. 16 represents the time history of the elastic flag’s tip displacement in the fully developed periodic state for the FSI3 benchmark. Fig. 16b shows the transverse (Y) displacement oscillating at ~5.4 Hz, and Fig. 16a demonstrates the longitudinal (X) retraction oscillating at the frequency ~10.7 Hz. The simulated tip displacement exhibits highly stable, self-sustained periodic oscillations. The longitudinal displacement oscillates at exactly twice the frequency of the transverse displacement, accurately reflecting the physical bending mechanics, in which the flag tip retracts towards the cylinder twice per full flapping cycle. The absence of numerical damping or amplitude decay over time highlights the robustness of the fully Lagrangian coupling strategy.

Figure 16: Tip displacement in fully developed state for (a) longitudinal displacement (X-axis), and (b) transverse displacement (Y-axis). The reference values are from Turek and Hron [95].

The computed displacement amplitudes and oscillation frequencies show excellent agreement with the reference finite element solutions. Stable evaluation of the hydrodynamic forces driving the flag’s motion was achieved without artificial normal flux corrections, thereby confirming the robustness of the proposed coupling algorithm for FSI problems.

3.4 Case 4: Dam with an Elastic Gate

Fluid flow and its interaction with a structure are simulated for the case of a dam with an elastic gate. The experiment conducted by Antoci et al. [74] featured a tank filled with water, confined by a rubber gate fixed at one end and free at the other. After the gate is suddenly freed, the water flows through a narrow opening created by the gate’s deformation. The hydrostatic pressure from the reservoir initiates the gate’s displacement, and the subsequent gate motion has a significant influence on the water flow patterns.

Fig. 17 shows the initial configuration of Case 4, with a water head in the reservoir measuring H3=0.14 m and a reservoir width of 0.10 m. The gate, with a thickness of 0.005 m and a length of L3=0.079 m, is clamped at the bottom of the reservoir wall, holding back the water. The gate is made of a flexible rubber material with a density of 600 kg/m3, Poisson’s ratio of 0.4, and Young’s modulus of 12 MPa. The fluid is water. The initial particle distance (l∘) is set at 0.001 m, resulting in approximately 18,000 particles.

Figure 17: Initial configuration of case 4.

Fig. 18 presents the simulation results for Case 4 at different time instances, alongside corresponding experimental snapshots for comparison. The color fields represent the fluid pressure and the first principal component of the strain tensor within the rubber gate. The first principal strain component offers a scalar measure of the maximum tensile strain direction, providing valuable insight into the deformation patterns and stress distribution in the elastic material.

Figure 18: Snapshots of case 4 from the developed numerical model and experiments [74]. Fluid colors represent dimensionless pressure (p/ρgH3); solid colors show strain.

Upon release, the rubber gate is forced open by the water pressure, which is proportional to the initial head in the full tank. This scenario highlights the dynamic fluid-structure interaction between the hydrostatic pressure exerted by the fluid and the elastic restoring forces of the gate. As the water discharges and the fluid level drops, the pressure acting on the gate gradually decreases. Once the pressure falls below a critical threshold, the gate begins to recover its original shape and move toward closure. In our simulation, this recovery process was observed to begin at approximately 0.12 s.

The results of the simulation show good agreement with experimental observations, which can be attributed to two key factors: (1) special attention was given to the treatment of geometric boundaries to find the surface normal vectors near the edges of the gate especially when it deforms, minimizing numerical artifacts and improving the simulation’s fidelity; (2) synchronization of the numerical sound speed used in the fluid and the solid material enhancing the stability of the coupling.

Figs. 19 and 20 present the time history of the displacement of the rubber gate’s free end. Overall, the simulation results align well with the experimental data, although some discrepancies exist. These differences are likely due to the assumption of a linear stress–strain relationship and may not fully capture the complex, highly nonlinear hyperelastic behavior of rubber under extreme strains, as well as the fluid pressure fluctuations in WC-MPS, which change the loading distribution on the gate. Similar discrepancies were observed in the results by Antoci et al. [74], and Hwang et al. [63]. However, our results match the experimental data more closely. Antoci et al. [74] overestimated the experimental data due to a two-level refinement of l∘=0.8 mm near the fluid-structure boundary and l∘=2.0 mm elsewhere. Hwang et al. [63] slightly underestimated and oscillated around the experimental data using a similar refinement to ours of l∘=1 mm, but by solving the PPE to update fluid pressure and using a Higher-order Laplacian (HL) to discretize the Laplacian of PPE.

Figure 19: Time history of displacement of the free end of the rubber gate: vertical displacement over time (Exp. denotes experimental data from [74] and Num. denotes the numerical method. Reference numerical data from [63,74]).

Figure 20: Time history of displacement of the free end of the rubber gate: horizontal displacement over time (Exp. denotes experimental data from [74] and Num. denotes the numerical method. Reference numerical data from [63,74]).

Fig. 21 shows the time variation of the water level, measured at the midpoint of the tank width. The observed trends are consistent with the experimental data and the results by Antoci et al. [74] and Hwang et al. [63].

Figure 21: Time variation of water level in the middle of the tank (Exp. denotes experimental data from [74] and Num. denotes the numerical method. Reference numerical data from [63,74]).

For the convergence study, the maximum tip displacement of the elastic gate (δmaxNum.) is compared with available numerical and experimental references (δmaxRef.). The peak displacement error is quantified as:

max(Error)=|δmaxNum.−δmaxRef.|L3+δmaxRef.(53)

The results (Table 2) show a good model convergence (~1.4-order convergence), confirming the stability of the hybrid TL–UL solid formulation and the fully Lagrangian coupling strategy under large deformation.

4  Conclusion

This study presented a fully Lagrangian, mesh-free framework for fluid–structure interaction based on the Moving Particle Semi-implicit (MPS) method. The proposed approach couples an enhanced weakly compressible MPS formulation for fluids with a hybrid Total–Updated Lagrangian formulation for hyperelastic solids, in which strains are evaluated in the reference configuration while momentum balance is enforced in the current configuration.

Several stabilization and consistency corrections were incorporated, including artificial density diffusion, particle regularization schemes, corrected kernel gradients, and suppression of spurious zero-energy modes. An optimized neighbor-search strategy was also introduced to improve computational efficiency. Stable and conservative coupling between phases was achieved through normal-flux surface integration and a velocity-consistent particle grouping strategy, supported by an interface particle identification criterion that improves boundary detection in deforming configurations.

The benchmark simulations demonstrate that the unified framework provides stable and consistent predictions for large-deformation FSI problems while maintaining the simplicity of Lagrangian particle methods.

Quantitative validation against four benchmark problems—cantilever beam vibration, dam-break flow, the Turek-Hron flexible flag benchmark, and a dam-break flow interacting with an elastic gate—demonstrated monotonic error convergence and excellent agreement with analytical, experimental, and numerical references. This confirms the framework’s ability to accurately capture large deformations, free-surface evolution, and nonlinear coupling. The results highlight the method’s efficiency and robustness for complex FSI applications.

It is important to note that the current framework is restricted to two-dimensional (2D) simulations. Specifically, the angular opening criterion utilized for interface detection is not directly extendible to 3D without modification. Extending this detection criterion using solid-angle-based measures and subsequently expanding the overall framework to 3D applications represent key directions for future research. Other future work will extend the framework to using multi-resolution schemes and explore nonlinear models for structural material to further enhance its predictive capability for real-world hydro-elastic and multiphase problems.

Acknowledgement: During the preparation of this work, the authors used Grammarly and ChatGPT to check grammatical correctness and to improve the readability and clarity of the text. After using these tools, the authors reviewed and edited the content as necessary and took full responsibility for the final content of the published article.

Funding Statement: This work was funded by Fonds de recherche du Québec—Nature et technologies (FRQ-NT) through the strategic research network CEISCE (Centre d’études interuniversitaire des structures sous charges extrêmes) and the Canada Research Chairs (CRC) Program.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Saeed Tavakoli and Ahmad Shakibaeinia; methodology, Saeed Tavakoli, Ahmad Shakibaeinia and Najib Bouaanani; software, Saeed Tavakoli; validation, Saeed Tavakoli and Ahmad Shakibaeinia; formal analysis, Saeed Tavakoli; investigation, Saeed Tavakoli; resources, Saeed Tavakoli and Ahmad Shakibaeinia; data curation, Saeed Tavakoli and Ahmad Shakibaeinia; writing—original draft preparation, Saeed Tavakoli; writing—review and editing, Saeed Tavakoli, Ahmad Shakibaeinia and Najib Bouaanani; visualization, Saeed Tavakoli and Ahmad Shakibaeinia; supervision, Ahmad Shakibaeinia and Najib Bouaanani; project administration, Ahmad Shakibaeinia and Najib Bouaanani; funding acquisition, Ahmad Shakibaeinia and Najib Bouaanani. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the Corresponding Author, Saeed Tavakoli, upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

1. Khayyer A, Gotoh H, Shimizu Y, Nishijima Y, Nakano A. 3D MPS-MPS coupled FSI solver for simulation of hydroelastic fluid-structure interactions in coastal engineering. J Japan Soc Civil Eng. 2020;76(2):I37–42. doi:10.2208/kaigan.76.2_I_37. [Google Scholar] [CrossRef]

2. Hassoon O, Tarfaoui M, Alaoui A, El Moumen A. Experimental and numerical investigation on the dynamic response of sandwich composite panels under hydrodynamic slamming loads. Compos Struct. 2017;178:297–307. doi:10.1016/j.compstruct.2017.07.014. [Google Scholar] [CrossRef]

3. Barrette P, Richard M, Poirier L, Babaei H, Frederking R. Assessing ice action on bridges in the context of climate change: prospective approach. Ottawa, ON, Canada: National Research Council of Canada; 2017. doi:10.4224/40000398. [Google Scholar] [CrossRef]

4. Rajendran R, Narasimhan K. Deformation and fracture behaviour of plate specimens subjected to underwater explosion—a review. Int J Impact Eng. 2006;32(12):1945–63. doi:10.1016/j.ijimpeng.2005.05.013. [Google Scholar] [CrossRef]

5. Ming F, Zhang A, Xue Y, Wang S. Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions. Ocean Eng. 2016;117(9–10):359–82. doi:10.1016/j.oceaneng.2016.03.040. [Google Scholar] [CrossRef]

6. Williamson C, Govardhan R. Vortex-induced vibrations. Annu Rev Fluid Mech. 2004;36(1):413–55. doi:10.1146/annurev.fluid.36.050802.122128. [Google Scholar] [CrossRef]

7. Mao W, Li K, Sun W. Fluid-structure interaction study of transcatheter aortic valve dynamics using smoothed particle hydrodynamics. Cardiovasc Eng Tech. 2016;7(4):374–88. doi:10.1007/s13239-016-0285-7. [Google Scholar] [PubMed] [CrossRef]

8. Szabó G, Györgyi J, Kristóf G. Three-dimensional FSI simulation by using a novel hybrid scaling–application to the Tacoma Narrows Bridge. Period Polytech Civ. 2020;64(4):975–88. doi:10.3311/PPci.15586. [Google Scholar] [CrossRef]

9. Lu W, Lubbad R, Løset S, Kashafutdinov M. Fracture of an ice floe: local out-of-plane flexural failures versus global in-plane splitting failure. Cold Reg Sci Technol. 2016;123(1):1–13. doi:10.1016/j.coldregions.2015.11.010. [Google Scholar] [CrossRef]

10. Khayyer A, Gotoh H, Shimizu Y. On systematic development of FSI solvers in the context of particle methods. J Hydrodyn. 2022;34(3):395–407. doi:10.1007/s42241-022-0042-3. [Google Scholar] [CrossRef]

11. Hou G, Wang J, Layton A. Review article: numerical methods for fluid-structure interaction—a review. Commun Comput Phys. 2012;12(2):337–77. doi:10.4208/cicp.291210.290411s. [Google Scholar] [CrossRef]

12. Zhang A, Sun P, Ming F, Colagrossi A. Smoothed particle hydrodynamics and its applications in fluid-structure interactions. J Hydrodyn. 2017;29(2):187–216. doi:10.1016/S1001-6058(16)60730-8. [Google Scholar] [CrossRef]

13. Liu M, Zhang Z. Smoothed particle hydrodynamics (SPH) for modeling fluid-structure interactions. Sci China Phys Mech. 2019;62(8):1–38. doi:10.1007/s11433-018-9357-0. [Google Scholar] [CrossRef]

14. Takizawa K, Tezduyar T. Multiscale space-time fluid-structure interaction techniques. Comput Mech. 2011;48(3):247–67. doi:10.1007/s00466-011-0571-z. [Google Scholar] [CrossRef]

15. Martínez-Ferrer P, Qian L, Ma Z, Causon D, Mingham C. An efficient finite-volume method to study the interaction of two-phase fluid flows with elastic structures. J Fluids Struct. 2018;83(8):54–71. doi:10.1016/j.jfluidstructs.2018.08.019. [Google Scholar] [CrossRef]

16. Mitsume N, Yoshimura S, Murotani K, Yamada T. Improved MPS-FE fluid-structure interaction coupled method with MPS polygon wall boundary model. Comput Model Eng Sci. 2014;101(4):229–47. doi:10.3970/cmes.2014.101.229. [Google Scholar] [CrossRef]

17. Long T, Huang C, Hu D, Liu M. Coupling edge-based smoothed finite element method with smoothed particle hydrodynamics for fluid structure interaction problems. Ocean Eng. 2021;225(11–14):108772. doi:10.1016/j.oceaneng.2021.108772. [Google Scholar] [CrossRef]

18. Zhang Y, Wan D, Hino T. Comparative study of MPS method and level-set method for sloshing flows. J Hydrodyn. 2014;26(4):577–85. doi:10.1016/S1001-6058(14)60065-2. [Google Scholar] [CrossRef]

19. Elahi R, Passandideh-Fard M, Javanshir A. Simulation of liquid sloshing in 2D containers using the volume of fluid method. Ocean Eng. 2015;96(2):226–44. doi:10.1016/j.oceaneng.2014.12.022. [Google Scholar] [CrossRef]

20. Sawada T, Hisada T. Fluid–structure interaction analysis of the two-dimensional flag-in-wind problem by an interface-tracking ALE finite element method. Comput Fluids. 2007;36(1):136–46. doi:10.1016/j.compfluid.2005.06.007. [Google Scholar] [CrossRef]

21. Hu H, Patankar N, Zhu M. Direct numerical simulations of fluid–solid systems using the arbitrary Lagrangian–Eulerian technique. J Comput Phys. 2001;169(2):427–62. doi:10.1006/jcph.2000.6592. [Google Scholar] [CrossRef]

22. Longatte E, Bendjeddou Z, Souli M. Methods for numerical study of tube bundle vibrations in cross-flows. J Fluids Struct. 2003;18(5):513–28. doi:10.1016/j.jfluidstructs.2003.08.010. [Google Scholar] [CrossRef]

23. Tezduyar T, Sathe S, Pausewang J, Schwaab M, Christopher J, Crabtree J. Interface projection techniques for fluid–structure interaction modeling with moving-mesh methods. Comput Mech. 2008;43(1):39–49. doi:10.1007/s00466-008-0261-7. [Google Scholar] [CrossRef]

24. Oñate E, Celigueta M, Idelsohn S, Salazar F, Suárez B. Possibilities of the particle finite element method for fluid–soil–structure interaction problems. Comput Mech. 2011;48(3):307–18. doi:10.1007/s00466-011-0617-2. [Google Scholar] [CrossRef]

25. Peskin C. The immersed boundary method. Acta Numer. 2002;11:479–517. doi:10.1017/S0962492902000077. [Google Scholar] [CrossRef]

26. Yang J, Stern F. A non-iterative direct forcing immersed boundary method for strongly-coupled fluid–solid interactions. J Comput Phys. 2015;295(3):779–804. doi:10.1016/j.jcp.2015.04.040. [Google Scholar] [CrossRef]

27. Zhao P, Xu J, Ge W, Wang J. A CFD-DEM-IBM method for Cartesian grid simulation of gas-solid flow in complex geometries. Chem Eng J. 2020;389(3):124343. doi:10.1016/j.cej.2020.124343. [Google Scholar] [CrossRef]

28. Meng Z, Zhang A, Yan J, Wang P, Khayyer A. A hydroelastic fluid–structure interaction solver based on the Riemann-SPH method. Comput Method Appl M. 2022;390(3):114522. doi:10.1016/j.cma.2021.114522. [Google Scholar] [CrossRef]

29. Gotoh H, Khayyer A, Shimizu Y. Entirely Lagrangian meshfree computational methods for hydroelastic fluid-structure interactions in ocean engineering. Appl Ocean Res. 2021;115(1):102822. doi:10.1016/j.apor.2021.102822. [Google Scholar] [CrossRef]

30. Gingold R, Monaghan J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc. 1977;181(3):375–89. doi:10.1093/mnras/181.3.375. [Google Scholar] [CrossRef]

31. Monaghan J. Smoothed particle hydrodynamics. Annu Rev Astron Astr. 1992;30(8):543–74. doi:10.1088/0034-4885/68/8/R01. [Google Scholar] [CrossRef]

32. Monaghan J. Simulating free surface flows with SPH. J Comput Phys. 1994;110(2):399–406. doi:10.1006/jcph.1994.1034. [Google Scholar] [CrossRef]

33. Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng. 1996;123(3):421–34. doi:10.13182/NSE96-A24205. [Google Scholar] [CrossRef]

34. Shakibaeinia A, Jin Y. MPS mesh-free particle method for multiphase flows. Comput Method Appl M. 2012;229–232(9):13–26. doi:10.1016/j.cma.2012.03.013. [Google Scholar] [CrossRef]

35. Sun P, Le Touze D, Oger G, Zhang A. An accurate FSI-SPH modeling of challenging fluid-structure interaction problems in two and three dimensions. Ocean Eng. 2021;221:108552. doi:10.1016/j.oceaneng.2020.108552. [Google Scholar] [CrossRef]

36. Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J Comput Phys. 2003;191(2):448–75. doi:10.1016/S0021-9991(03)00324-3. [Google Scholar] [CrossRef]

37. Shakibaeinia A, Jin Y. A weakly compressible MPS method for modeling of open-boundary free-surface flow. Int J Numer Meth Fl. 2010;63(10):1208–32. doi:10.1002/fld.2132. [Google Scholar] [CrossRef]

38. Robb D, Gaskin S, Marongiu J. SPH-DEM model for free-surface flows containing solids applied to river ice jams. J Hydraul Res. 2016;54(1):27–40. doi:10.1080/00221686.2015.1131203. [Google Scholar] [CrossRef]

39. Zhang N, Zheng X, Ma Q, Hu Z. A numerical study on ice failure process and ice-ship interactions by smoothed particle hydrodynamics. Int J Nav Archit Oc. 2019;11(2):796–808. doi:10.1016/j.ijnaoe.2019.02.008. [Google Scholar] [CrossRef]

40. Long X, Liu L, Liu S, Ji S. Discrete element analysis of high-pressure zones of sea ice on vertical structures. J Mar Sci Eng. 2021;9(3):348. doi:10.3390/jmse9030348. [Google Scholar] [CrossRef]

41. Bonet J, Lok T. Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Comput Method Appl M. 1999;180(1–2):97–115. doi:10.1016/S0045-7825(99)00051-1. [Google Scholar] [CrossRef]

42. Shimizu Y, Khayyer A, Gotoh H. An SPH-based fully-Lagrangian meshfree implicit FSI solver with high-order discretization terms. Eng Anal Bound Elem. 2022;137(1):160–81. doi:10.1016/j.enganabound.2021.10.023. [Google Scholar] [CrossRef]

43. Shimizu Y, Khayyer A, Gotoh H. An implicit SPH-based structure model for accurate fluid–structure interaction simulations with hourglass control scheme. Eur J Mech B-Fluid. 2022;96(1):122–45. doi:10.1016/j.euromechflu.2022.07.007. [Google Scholar] [CrossRef]

44. Khayyer A, Shimizu Y, Gotoh H, Hattori S. A 3D SPH-based entirely Lagrangian meshfree hydroelastic FSI solver for anisotropic composite structures. Appl Math Model. 2022;112(1):560–613. doi:10.1016/j.apm.2022.07.031. [Google Scholar] [CrossRef]

45. Khayyer A, Gotoh H, Falahaty H, Shimizu Y, Nishijima Y. Towards development of a reliable fully-Lagrangian MPS-based FSI solver for simulation of 2D hydroelastic slamming. Ocean Syst Eng. 2017;7(3):299–318. doi:10.12989/ose.2017.7.3.299. [Google Scholar] [CrossRef]

46. Daxini S, Prajapati J. A review on recent contribution of meshfree methods to structure and fracture mechanics applications. Sci World J. 2014;2014(1):247172. doi:10.1155/2014/247172. [Google Scholar] [PubMed] [CrossRef]

47. Chen Y, Lee J, Eskandarian A. Meshless methods in solid mechanics. Vol. 9. New York, NY, USA: Springer; 2006. doi:10.1007/0-387-33368-1. [Google Scholar] [CrossRef]

48. Benson D, Bazilevs Y, Hsu M, Hughes T. Isogeometric shell analysis: the Reissner–Mindlin shell. Comput Method Appl M. 2010;199(5–8):276–89. doi:10.1016/j.cma.2009.05.011. [Google Scholar] [CrossRef]

49. Francomano E, Paliaga M. Highlighting numerical insights of an efficient SPH method. Appl Math Comput. 2018;339(3):899–915. doi:10.1016/j.amc.2018.07.060. [Google Scholar] [CrossRef]

50. Jandaghian M, Shakibaeinia A. An enhanced weakly-compressible MPS method for free-surface flows. Comput Method Appl M Eng. 2020;360(1):112771. doi:10.1016/j.cma.2019.112771. [Google Scholar] [CrossRef]

51. Tanaka M, Masunaga T. Stabilization and smoothing of pressure in MPS method by quasi-compressibility. J Comput Phys. 2010;229(11):4279–90. doi:10.1016/j.jcp.2010.02.011. [Google Scholar] [CrossRef]

52. Jandaghian M, Krimi A, Zarrati A, Shakibaeinia A. Enhanced weakly-compressible MPS method for violent free-surface flows: role of particle regularization techniques. J Comput Phys. 2021;434(8):10902716. doi:10.1016/j.jcp.2021.110202. [Google Scholar] [CrossRef]

53. Zhang C, Wei Y, Dias F, Hu X. An efficient fully Lagrangian solver for modeling wave interaction with oscillating wave surge converter. Ocean Eng. 2021;236(11–14):109540. doi:10.1016/j.oceaneng.2021.109540. [Google Scholar] [CrossRef]

54. Zhang C, Rezavand M, Hu X. A multi-resolution SPH method for fluid-structure interactions. J Comput Phys. 2021;429(3):110028. doi:10.1016/j.jcp.2020.110028. [Google Scholar] [CrossRef]

55. Shibata K, Koshizuka S, Matsunaga T, Masaie I. The overlapping particle technique for multi-resolution simulation of particle methods. Comput Method Appl M. 2017;325(3):434–62. doi:10.1016/j.cma.2017.06.030. [Google Scholar] [CrossRef]

56. Khayyer A, Shimizu Y, Gotoh H, Hattori S. Multi-resolution ISPH-SPH for accurate and efficient simulation of hydroelastic fluid-structure interactions in ocean engineering. Ocean Eng. 2021;226(11):108652. doi:10.1016/j.oceaneng.2021.108652. [Google Scholar] [CrossRef]

57. Liu G. An overview on meshfree methods: for computational solid mechanics. Int J Comput Meth. 2016;13(5):1630001. doi:10.1142/S0219876216300014. [Google Scholar] [CrossRef]

58. Gavete L, Gavete M, Benito J. Improvements of generalized finite difference method and comparison with other meshless method. Appl Math Model. 2003;27(10):831–47. doi:10.1016/S0307-904X(03)00091-X. [Google Scholar] [CrossRef]

59. Shankar V, Wright GB, Kirby RM, Fogelson AL. A radial basis function (RBF)-finite difference (FD) method for diffusion and reaction–diffusion equations on surfaces. J Sci Comput. 2015;63(3):745–68. doi:10.1007/s10915-014-9914-1. [Google Scholar] [PubMed] [CrossRef]

60. Shimizu Y, Khayyer A, Gotoh H. Towards development of Lagrangian meshfree hydroelastic FSI solvers by incorporating implicit structure solvers. Proc Soc Civil Eng B2. 2019;75(2):799–804. doi:10.2208/kaigan.75.i_799. [Google Scholar] [CrossRef]

61. Zhang G, Chen X, Wan D. MPS-FEM coupled method for study of wave-structure interaction. J Mar Sci Appl. 2019;18(4):387–99. doi:10.1007/s11804-019-00105-6. [Google Scholar] [CrossRef]

62. Peng Y, Zhang A, Wang S. Coupling of WCSPH and RKPM for the simulation of incompressible fluid–structure interactions. J Fluids Struct. 2021;102(11):103254. doi:10.1016/j.jfluidstructs.2021.103254. [Google Scholar] [CrossRef]

63. Hwang S, Khayyer A, Gotoh H, Park J. Development of a fully Lagrangian MPS-based coupled method for simulation of fluid-structure interaction problems. J Fluids Struct. 2014;50(11):497–511. doi:10.1016/j.jfluidstructs.2014.07.007. [Google Scholar] [CrossRef]

64. Khayyer A, Gotoh H, Falahaty H, Shimizu Y. Towards development of enhanced fully-Lagrangian mesh-free computational methods for fluid-structure interaction. J Hydrodyn. 2018;30(1):49–61. doi:10.1007/s42241-018-0005-x. [Google Scholar] [CrossRef]

65. Khayyer A, Gotoh H, Falahaty H, Shimizu Y. An enhanced ISPH–SPH coupled method for simulation of incompressible fluid–elastic structure interactions. Comput Phys Commun. 2018;232(1):139–64. doi:10.1016/j.cpc.2018.05.012. [Google Scholar] [CrossRef]

66. Amaro Junior R, Gay Neto A, Cheng L. Three-dimensional weakly compressible moving particle simulation coupled with geometrically nonlinear shell for hydro-elastic free-surface flows. Int J Numer Meth Fl. 2022;94(8):1048–81. doi:10.1002/fld.5083. [Google Scholar] [CrossRef]

67. Han L, Hu X. SPH modeling of fluid-structure interaction. J Hydrodyn. 2018;30(1):62–9. doi:10.1007/s42241-018-0006-9. [Google Scholar] [CrossRef]

68. Sun P, Le Touzé D, Zhang A. Study of a complex fluid-structure dam-breaking benchmark problem using a multi-phase SPH method with APR. Eng Anal Bound Elem. 2019;104(21):240–58. doi:10.1016/j.enganabound.2019.03.033. [Google Scholar] [CrossRef]

69. Sun Z, Djidjeli K, Xing J, Cheng F. Modified MPS method for the 2D fluid structure interaction problem with free surface. Comput Fluids. 2015;122:47–65. doi:10.1016/j.compfluid.2015.08.017. [Google Scholar] [CrossRef]

70. Hu T, Wang S, Zhang G, Sun Z, Zhou B. Numerical simulations of sloshing flows with an elastic baffle using a SPH-SPIM coupled method. Appl Ocean Res. 2019;93:101950. doi:10.1016/j.apor.2019.101950. [Google Scholar] [CrossRef]

71. Hwang S, Park J, Gotoh H, Khayyer A, Kang K. Numerical simulations of sloshing flows with elastic baffles by using a particle-based fluid–structure interaction analysis method. Ocean Eng. 2016;118(11):227–41. doi:10.1016/j.oceaneng.2016.04.006. [Google Scholar] [CrossRef]

72. Yang C, Zhang H, Su H, Shen Z. Numerical simulation of sloshing using the MPS-FSI method with large eddy simulation. China Ocean Eng. 2018;32(3):278–87. doi:10.1007/s13344-018-0029-6. [Google Scholar] [CrossRef]

73. Liu M, Shao J, Li H. Numerical simulation of hydro-elastic problems with smoothed particle hydrodynamics method. J Hydrodyn. 2013;25(5):673–82. doi:10.1016/S1001-6058(13)60412-6. [Google Scholar] [CrossRef]

74. Antoci C, Gallati M, Sibilla S. Numerical simulation of fluid–structure interaction by SPH. Comput Struct. 2007;85(11–14):879–90. doi:10.1016/j.compstruc.2007.01.002. [Google Scholar] [CrossRef]

75. Bui H, Fukagawa R, Sako K, Ohno S. Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic–plastic soil constitutive model. Int J Numer Anal Met Geomech. 2008;32(12):1537–70. doi:10.1002/nag.688. [Google Scholar] [CrossRef]

76. Khayyer A, Tsuruta N, Shimizu Y, Gotoh H. Multi-resolution MPS for incompressible fluid-elastic structure interactions in ocean engineering. Appl Ocean Res. 2019;82(1):397–414. doi:10.1016/j.apor.2018.10.020. [Google Scholar] [CrossRef]

77. Bao T, Hu J, Wang S, Huang C, Yu Y, Shakibaeinia A. An entirely SPH-based FSI solver and numerical investigations on hydrodynamic characteristics of the flexible structure with an ultra-thin characteristic. Comput Method Appl M. 2024;431:117255. doi:10.1016/j.cma.2024.117255. [Google Scholar] [CrossRef]

78. O’Connor J, Rogers B. A fluid-structure interaction model for free-surface flows and flexible structures using smoothed particle hydrodynamics on a GPU. J Fluids Struct. 2021;104(21):103312. doi:10.1016/j.jfluidstructs.2021.103312. [Google Scholar] [CrossRef]

79. Zhan L, Peng C, Zhang B, Wu W. A stabilized TL–WC SPH approach with GPU acceleration for three-dimensional fluid-structure interaction. J Fluids Struct. 2019;86(21):329–53. doi:10.1016/j.jfluidstructs.2019.02.002. [Google Scholar] [CrossRef]

80. Batchelor G. An introduction to fluid dynamics. Cambridge, UK: Cambridge University Press; 2000. [Google Scholar]

81. Song M, Koshizuka S, Oka Y. Dynamic analysis of elastic solids by MPS method. Nippon Kikai Gakkai Ronbunshu A. 2005;17(1):16–22. doi:10.1299/kikaia.71.16. [Google Scholar] [CrossRef]

82. Wu D, Zhang C, Tang X, Hu X. An essentially non-hourglass formulation for total Lagrangian smoothed particle hydrodynamics. Comput Method Appl M. 2023;407(3):115915. doi:10.1016/j.cma.2023.115915. [Google Scholar] [CrossRef]

83. Ren D, Park J, Hwang S, Jeong S, Kim H. Failure simulation of ice beam using a fully Lagrangian particle method. Int J Nav Archit Oc. 2019;11(2):639–47. doi:10.1016/j.ijnaoe.2019.01.001. [Google Scholar] [CrossRef]

84. Rafiee A, Thiagarajan K. An SPH projection method for simulating fluid-hypoelastic structure interaction. Comput Method Appl M. 2009;198(33–36):2785–95. doi:10.1016/j.cma.2009.04.001. [Google Scholar] [CrossRef]

85. Jandaghian M, Siaben H, Shakibaeinia A. Stability and accuracy of the weakly compressible SPH with particle regularization techniques. Eur J Mech B-Fluid. 2022;94(3):314–33. doi:10.1016/j.euromechflu.2022.03.007. [Google Scholar] [CrossRef]

86. Long T, Hu D, Wan D, Zhuang C, Yang G. An arbitrary boundary with ghost particles incorporated in coupled FEM–SPH model for FSI problems. J Comput Phys. 2017;350:166–83. doi:10.1016/j.jcp.2017.08.044. [Google Scholar] [CrossRef]

87. Zhang C, Zhu Y, Lyu X, Hu X. An efficient and generalized solid boundary condition for SPH: applications to multi-phase flow and fluid-structure interaction. Eur J Mech B-Fluid. 2022;94(3):276–92. doi:10.1016/j.euromechflu.2022.03.011. [Google Scholar] [CrossRef]

88. Adami S, Hu X, Adams N. A generalized wall boundary condition for smoothed particle hydrodynamics. J Comput Phys. 2012;231(21):7057–75. doi:10.1016/j.jcp.2012.05.005. [Google Scholar] [CrossRef]

89. Lind S, Xu R, Stansby P, Rogers B. 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. 2012;231(4):1499–523. doi:10.1016/j.jcp.2011.10.027. [Google Scholar] [CrossRef]

90. Sun P, Colagrossi A, Marrone S, Antuono M, Zhang A. A consistent approach to particle shifting in the δ-Plus-SPH model. Comput Method Appl M. 2019;348(13):912–34. doi:10.1016/j.cma.2019.01.045. [Google Scholar] [CrossRef]

91. Gray J, Monaghan J, Swift R. SPH elastic dynamics. Comput Method Appl M. 2001;190(49–50):6641–62. doi:10.1016/S0045-7825(01)00254-7. [Google Scholar] [CrossRef]

92. Lobovský L, Botia-Vera E, Castellana F, Mas-Soler J, Souto-Iglesias A. Experimental investigation of dynamic pressure loads during dam break. J Fluids Struct. 2014;48(10):407–34. doi:10.1016/j.jfluidstructs.2014.03.009. [Google Scholar] [CrossRef]

93. Antuono M, Colagrossi A, Marrone S. Numerical diffusive terms in weakly-compressible SPH schemes. Comput Phys Commun. 2012;183(12):2570–80. doi:10.1016/j.cpc.2012.07.006. [Google Scholar] [CrossRef]

94. Rezavand M, Zhang C, Hu X. A weakly compressible SPH method for violent multi-phase flows with high density ratio. J Comput Phys. 2020;402(1):109092. doi:10.1016/j.jcp.2019.109092. [Google Scholar] [CrossRef]

95. Turek S, Hron J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. In: Fluid-structure interaction. Berlin/Heidelberg, Germany: Springer; 2006. p. 371–85. doi:10.1007/3-540-34596-5_15. [Google Scholar] [CrossRef]