N. Mitsume Advanced Modeling and Simulation
in Engineering Sciences(2023)10:11
https://doi.org/10.1186/s40323-023-00248-w

Advanced Modeling and Simulation
in Engineering Sciences

RESEARCH ARTICLE

Open AccessCompatible interface wave–structureinteraction model for combining
mesh-free particle and finite element
methods
Naoto Mitsume1* Correspondence:

mitsume@kz.tsukuba.ac.jp
1
Faculty of Engineering,
Information and Systems,
University of Tsukuba, Tsukuba,
Japan

Abstract
This study presents a novel wave–structure interaction model, which is a compatible
interface wave–structure interaction model that is based on mesh-free particle
methods for free-surface ﬂow analysis; the FEM for structural analysis. We adopt the
explicitly represented polygon (ERP) wall boundary model, which is a polygon wallboundary model for mesh-free particle methods, to express the ﬂuid–structureinterfaces. The ﬂuid–structure interfaces in the proposed model are geometricallycompatible because the ERP model has advantages in dealing with complex-shaped
and moving boundaries and it enables the direct use of surface meshes and shapefunctions given by ﬁnite element models of structures. This allows the automatic
generation of polygons for free-surface ﬂow analysis from ﬁnite elements, andtherefore, it greatly increases the ﬂexibility of the analysis. Based on the compatible
interface between the polygons and ﬁnite elements, we propose a strong couplingalgorithm based on an iterative partitioned scheme that deﬁnes the interacting models
between the ﬂuid and the structure. We model the force exerted by ﬂuid particles on
structures such that the kinetic boundary condition on the ﬂuid–structure interface issatisﬁed. We perform the veriﬁcation and validation tests of the proposed model bysolving two benchmark problems.
Keywords: Finite element method, Mesh-free particle method, Fluid–structureinteraction, Free-surface ﬂow, Partitioned iterative coupling method

Introduction
A large tsunami can inﬂict devastating damage, as evidenced by the Great East Japan
Earthquake and Tsunami that occurred on March 11, 2011. Severe accidents triggeredby tsunamis such as the Fukushima Daiichi nuclear disaster have the potential to cause
catastrophic harm to people; this has led to an increased awareness regarding the vital
importance of designing safe, tsunami-resistant structures (including the buildings andmachinery of electric power, energy, and chemical plants) and evacuation facilities incoastal regions.
© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit
to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The
images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise
in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyrightholder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

0123456789().,–: volV

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Thus far, several experimental studies have evaluated pressure forces generated bytsunamis and proposed practical experimental formulas [1–3]. These Experimental formulas for estimating collision forces attributed to ﬂoating debris [4,5] have also been presented. Although these formulas are very useful in the design process of defense structures
such as levees, numerical analyses can be considered viable alternatives. Two-dimensional
wave models, as represented by the shallow-water equations and Boussinesq-type wave
models, have been used to understand large-scale wave behaviors [6–8] because thesemodels can save computational costs compared to that incurred by running direct simulations of the three-dimensional Navier–Stokes equations. The inundation area for tsunamis
can be predicted using an appropriate representation of bottom friction, Coriolis force,and any other parameter [9]. In contrast with the two-dimensional numerical simulations, direct simulations based on Navier–Stokes equations [10–12] have been conducted
to assess wave forces exerted by tsunamis. Although these simulations can consider thethree-dimensional nature of waves accurately, they often incur considerably high computational costs for computing large-scale wave behaviors. A multiscale approach to wave
modeling is employed to reduce the computational burden; this approach exploits boththe two- and three-dimensional computations. For example, a combination of the twodimensional shallow-water equations and three-dimensional Navier–Stokes equations[13–16] has recently attracted considerable research attention.
Although existing studies have played an important role in tsunami disaster prevention, it is sometimes economically impossible to mitigate all the deleterious eﬀects ofsuch disasters. Therefore, structural design that is resilient to the catastrophic damage
incurred by tsunamis, i.e., structural design that can maintain functions even for beyondthe-design-basis conditions [17] has recently gained increasing attention. Damage createdby a tsunami is attributed to many factors: the direct force of the wave, elastic and plastic
deformation, collisions with ﬂoating debris, and buoyancy [18]. All these factors must be
considered when developing tsunami-resilient designs; however, it is diﬃcult to predict
them using conventional processes that are based only on empirical knowledge and wave
simulations. An integrative strategy that includes the wave–structure interaction analysis
is a powerful and important part of the tsunami-resilient design processes.
The wave–structure interaction problem can be considered a type of ﬂuid–structureinteraction (FSI) problem with free surfaces. Considerable research has been conductedto solve FSI problems using the ﬁnite element method (FEM). Given the high accuracy
and widespread application of the FEM, it is the de facto standard method for structural
analysis. However, the analysis of free-surface ﬂow using FEM is not straightforward.
Interface capturing methods such as the volume of ﬂuid method [19] and the level set
method [20,21] use an artiﬁcial scalar ﬁeld to describe the interfaces implicitly. These
approaches can perform robust calculations even if the interfaces are subjected to severe
motion caused by topological changes. However, the interface capturing methods suffer from mass conservation and interface smearing problems. In addition to the above
approaches, which are Eulerian or ALE formulations in ﬂuids, the particle FEM (PFEM)[22–26] uses Lagrangian formulations in the ﬂuid and eliminates the problems of freesurfaces and moving boundaries. However, PFEM requires frequent remeshing.
In contrast with mesh-based methods, mesh-free particle methods such as smoothed
particle hydrodynamics (SPH) [27,28] and moving particle semi-implicit/simulation(MPS) [29] are inherently Lagrangian methods wherein a continuum is discretized asPage 2 of 31

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

a group of moving particles. These methods continue to clearly deﬁne moving boundaries
that include free surfaces even during long analyses because of their Lagrangian formulations. Coupled analyses that include a ﬂuid-rigid body and ﬂuid–structure interactions
have been conducted using mesh-free particle methods [30–33]. However, these methods
are less accurate on solid surfaces and have limited applicability to structural analysis,
such as with the ﬁnite diﬀerence method, because the Neumann boundary conditions on
solid surfaces must be explicitly imposed and the derivatives needs to be approximatedby particles within the circular support domain unlike the FEM, which solves weak form
equations.
Hybrid coupled models that use mesh-free particle methods and FEM have been proposed to exploit the advantages of both methods. Early studies on coupled analysis related
to these models [34–39] were conducted to deal with large deformation, penetration,
crack propagation, and fragmentation, wherein structure–structure interaction problems
were solved. Existing studies on FSI analyses using the FEM and the particle methods [40–43], in which commercial programs such as the PAM-CRASH and LS-DYNA are often
used, have also been conducted. In recent years, existing studies on modeling the hybrid
FSI model [42,44–49] have focused on approaches to formulate interactions between thefree-surface ﬂow computation based on particle methods and the structural computation
based on the FEM. Lee et al. [45] developed a partitioned coupled FSI model between the
MITC4 shell element based on structural computation and the MPS based free-surface
ﬂow computation, wherein wall particles and nodes of the shell elements overlapped to
exchange physical values. Fourey et al. [42] proposed a model between the SPH methodwith the ghost particle model as a wall boundary model and the FEM. They obtained force
from ﬂuid to structure via the surface integration of pressures on the ﬂuid–structureinterface interpolated by the weighted average of pressures on the neighboring particles.
Yang et al. [49,50] developed a model using the SPH method and the FEM. As with theLee’s study, the wall particles and nodes of the ﬁnite elements were overlapped at the same
positions, and the force exerted on the surface of the structure was obtained considering
the reaction forces in the ﬂuid computation. They conducted computations for quantitative validation compared with experimental results. Thiyahuddin et al. [48] developed
a model using the SPH method and the FEM to solve for a car colliding with water-ﬁlled
barriers. They modeled the ﬂuid–structure interface with a ﬁnite thickness for preventing
the penetrations of particles into structures. Recently, Li et al. [46,47] proposed a partitioned coupled FSI model between the FEM and SPH-ALE method [51]. They discussed
the nonphysical creation and dissipation of energy on the ﬂuid–structure interface anddeveloped a formulation to prevent it. However, these hybrid coupled models have mainly
employed weak coupling schemes based on staggered algorithms such as the conventional
serial staggered (CSS) scheme, and do not implicitly treat ﬂuid–structure interfaces.
Similar to this trend, we developed MPS-FE methods [52–54] that adopts the MPS
method for ﬂuid computation involving free surfaces and the FEM for structure computation. These two methods are coupled with a partitioned approach, i.e., the CSS scheme[55], which can set diﬀerent time step sizes for the ﬂuid and structure computations. The
method combines advantages of both methods for achieving eﬃciency and robustness.
However, the conventional MPS-FE method in which the MPS wall boundary particles andﬁnite elements are overlapped to exchange information on ﬂuid–structure interfaces cannot deal with complex-shaped ﬂuid–structure boundaries because the wall particles need

Page 3 of 31

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 4 of 31

to be set in an orthogonal and uniform grid manner for the accurate execution of the MPS
computation. This requires cumbersome interpolation for the exchange of physical values
based on the node-particle correspondence relation. Therefore, the MPS-FE method lacks
versatility and reduces software modularity. Further, forces on ﬂuid–structure interfaces
are not balanced when the pressure on the walls is calculated in this manner.
This study presents a novel wave–structure interaction model, which is a compatible
interface wave–structure interaction model that is based on mesh-free particle methods
for free-surface ﬂow analysis; the FEM for structural analysis. We adopt the explicitly
represented polygon (ERP) wall boundary model [56], which is a polygon wall boundary model for mesh-free particle methods, to express the ﬂuid–structure interfaces. The
ﬂuid–structure interfaces in the proposed model are geometrically compatible because
the ERP model has advantages in dealing with complex-shaped and moving boundaries
and it enables the direct use of surface meshes and shape functions given by ﬁnite element
models of structures. This allows the automatic generation of polygons for free-surface
ﬂow analysis from ﬁnite elements, and therefore, it greatly increases the ﬂexibility of theanalysis. Detailed formulations of the ERP model introduced to the explicit MPS method,which is a mesh-free particle method with explicit pressure calculation, is described in“Free-surface ﬂow analysis by explicit moving particle simulation (MPS) method with
explicitly represented polygon (ERP) wall boundary model”. In order to map physicalquantities between particle methods for strong-form equations and ﬁnite element methods for weak-form equations, we propose a novel formulation based on the Dirac delta
function, where the force exerted on the structure by each ﬂuid particle is considered as apoint load. It is formulated such that the kinetic boundary condition on the ﬂuid–structureinterface is satisﬁed. In addition, we introduce a strong coupling algorithm based on an
iterative partitioned scheme to discuss its eﬀect on numerical results and the convergenceof iterations, which has not been discussed much in existing studies on the hybrid coupledmodels. The details on these procedures are introduced in “Veriﬁcation and validation”.
We perform the veriﬁcation and validation tests of the proposed model by solving two
benchmark problems. First, we solve a dam break problem with an elastic obstacle; theresult of the proposed model are compared those of one of the several existing studies.
Then, we simulate the experiment of a sloshing tank with an elastic beam; these test cases
are then introduced in “Conclusions”.

Governing equations of fluid–structure interaction (FSI) problem with freesurfaces
Deﬁnition of spaces

Figure 1 shows the spaces and boundaries in an FSI problem involving free surfaces. The
whole area of the ﬂuid–structure interaction analysis, , which is partitioned into areasof free-surface ﬂow analysis and structural analysis, F and S , respectively, is given as = F ∪ S .

(1)

The boundaries of the free-surface ﬂow and structural analyses are denoted as F andS , respectively.
The domain of the free-surface ﬂow analysis is partitioned into the ﬂuid domain, F ,F ,and the other domain, F ,V . A free-surface boundary, F ,FS , is deﬁned between theseN. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 5 of 31

fluid domain

structural domain

Fig. 1 Schematic of the ﬂuid–structure interaction problem with free surface

two domains. In this study, we consider the open boundary, F ,O , and wall boundary,F ,W , on the boundary of the free-surface ﬂow analysis. We do not consider inﬂow
boundaries in this study. The domain of structural analysis is partitioned into boundaries
with displacement and traction boundary conditions, S ,D and S ,L , respectively.
Let I represent the ﬂuid–structure interface. Then, the FSI boundaries in the freesurface ﬂow and structural analyses denoted as, F ,I and S ,I , can be considered diﬀerent.
Governing equations and boundary conditions

The governing equations in both the ﬂuid and the structure are Cauchy’s 1st law of motion,D(ρv) = ∇ · σ + ρg,
Dt

(2)

D
where Dt
, ρ, v, σ, and g represent the material derivative, density, velocity vector, Cauchy
stress tensor, and gravitational acceleration vector, respectively.
In our FSI problem, the constitutive laws in the ﬂuid and structure in terms of stress are
diﬀerent. Let σ F and σ S represent the stress tensors in the ﬂuid and structural analyses,
respectively.

σ F = −pI + μ ∇v + (∇v)TσS =1F SF T
det F

in F ,F ,(3)

in S ,

(4)

where p, μ, F , and S represent the pressure, viscosity, deformation gradient tensor, and2nd Piora–Kirchhoﬀ stress tensor, respectively.
On the ﬂuid–structure interface, the following geometric and kinetic boundary conditions need to be imposed.

F F

vF − vS = 0

on I ,(5)S S

on I ,

(6)

σ n +σ n =0

where v F and v S represent the velocities in the ﬂuid and structure domain, and nF and nS
represent the outward unit normal vectors in the ﬂuid and structure domains.

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 6 of 31

Table 1 Deﬁnition of indices for the formulation of compatible interface wave–structure interaction
model
Fluid particle:i(i ∈ particle = {1, · · · , Nparticle })Neighboring particle within eﬀective radius:
Polygon wall:
Polygon wall on the ﬂuid–structure (FS) interface:jss̃(s ∈ pol = {1, · · · , Npol })˜ pol ⊂ pol )(s̃ ∈ Global node of the polygon wall:
Global node of the polygon wall on the FS interface:
Finite element (FE):LL̃e(L ∈ pnode = {1, · · · , Npnode })˜ pnode ⊂ pnode )(L̃ ∈ (e ∈ elem = {1, · · · , Nelem })FE surface:
FE surface on the FS interface:
Global node of FE:bb̃
K

Global node of the FE on the FS interface:
Timestep:K̃n(b ∈ surface = {1, · · · , Nsurface })˜ surface ⊂ surface )(b̃ ∈ (K ∈ fnode = {1, · · · , Nfnode })˜ fnode ⊂ fnode )(K̃ ∈ (1 ≤ n ≤ Ntimestep )Rule of indices

The FSI model has many ingredients such as ﬂuid particles, polygon walls, ﬁnite elements,
and nodes. For the ease of formulation, we adopt a rule in which each character used as
an index with a physical value speciﬁes an ingredient the value is deﬁned on, as indicatedin Table 1. In this table,  represents index set.

Free-surface flow analysis by explicit moving particle simulation (MPS) methodwith explicitly represented polygon (ERP) Wall Boundary Model
Governing equations

Substituting Eq. (3) into Eq. (2) produces the Navier–Stokes equations in a Lagrangianreference frame asDv
1= − ∇p + ν∇ 2 v + g,
Dt
ρ

(7)

where p represents the pressure and ν represents the kinetic viscosity deﬁned as ν = μρ .
Diﬀerential operator models in MPS method

In the MPS discretization, the diﬀerential operators acting on particle i are evaluated usingneighboring particles, j, located within an eﬀective radius, re .
The neighboring particles are weighted by a function with respect to the relative distance, r. In this study, we use the quadratic function⎧2⎪⎪⎨ r −1(0 ≤ r < re )w(r) =
re⎪⎪⎩0(re ≤ r).

(8)

The normalization factor, called the particle number density, is deﬁned as
Nparticleni =w(|xj − xi |),
j=1
j =i

(9)

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 7 of 31

where Nparticle represents the number of particles. The particle number density is proportional to the density of the ﬂuid.
In the EMPS method, diﬀerential operators in the governing equations are computed asd∇p i =

Nparticle

n0

j=1
j =i

2d
λ0 n0

∇ 2 v i =Nparticle
pj + p i x j − x iw(|xj − xi |) ,|xj − xi |2
v j − v i w(|xj − xi |) .

(10)

(11)

j=1
j =i

Here, d, n0 , and  represent the number of dimensions, initial value of the particle number
density given by Eq. (9) and calculated for the initial particle geometry, and discretizationby the MPS diﬀerential model, respectively. Further, λ0 represents a correction parameter,
which ensures that the increase in the variance is equal to that of the analytical solution;
like n0 , it is calculated for the initial geometry as

λi =

1
ni

Nparticle|xj − xi |2 w(|xj − xi |).

(12)

j=1
j =i

Time-stepping algorithm in the EMPS method

In the original MPS method [29], the fractional step method is applied for time discretization; this method splits each time step into two pseudo-time steps of prediction
and correction. Therefore, the velocity in the new time step n + 1 is given by the sum of
the predicted velocity, v ∗ , and the velocity correction value, v , as v n+1 = v ∗ + v . The
predicted velocity v ∗ is explicitly calculated using viscosity terms and source terms asv ∗i − v ni
= ν∇ 2 v ni + g.
t

(13)

Next, the velocity at n + 1 step is given byv n+1− v ∗i1i= − ∇p n+1
.
i
t
ρ

(14)

Although the original MPS method solves the pressure Poisson equation implicitly,
the EMPS method [57,58] explicitly calculates the pressure by introducing a pseudocompressibility. The pressure is calculated as
pin+1 = c2

ρ ∗(n − n0 ),
n0 i

(15)

where c represents the speed of sound, which can be set arbitrarily to satisfy the conditionof stability and incompressibility.

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Wall particle

Fluid particle

Page 8 of 31

Fluid particle

Polygon walliiWall Particle Model

Polygon Wall Model

Ωwall
i

Fig. 2 Schematic of the conventional wall model that represents wall boundaries as particles (wall particle
model) and the polygon wall boundary model

Free-surface criterion

For particles on a free surface, F ,FS , we apply the Dirichlet condition for the pressure,p = 0. This is a simple approximation of the free-surface boundary condition.
σ F nF = 0

on F ,FS .

(16)

The condition for free-surface particle recognition is given asn∗i < βn0 ,(17)where β represents a threshold coeﬃcient. In the EMPS model, we set β = 1.0 to prevent
instability caused by tensile stress.
ERP wall boundary model
Polygon detection

The ERP model [56] expresses the wall boundaries as polygons, which are explicitly represented without using the signed distance function; we assumes that each ﬂuid particle isaﬀected by only the nearest polygon. Here, we denote the nearest points on each polygon
wall
near
wall s from the particle i as xnear
i,s . The closest point, x i , for all x i,s is deﬁned as
xwall= argmin |xi − xnear
ii,s | (s = 1, 2, · · · , Npol ),
xnear
i,s

(18)

where Npol represents the number of polygons. The outward unit normal vector at xnear
i,s
is deﬁned as
=
nwall
i

xi − xwalli|xi − xwall
i |

.

(19)

Corrections of particle number density

Figure 2 shows that the conventional wall model for the MPS method [29] represents wallboundaries as particles (wall particles). Although the ERP model does not require thesewall particles, the particle number density inside the wall region shown in Fig. 2, walli ,needs to be computed. The wall weight function, z(r), is computed to interpolate thisN. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 9 of 31

value in the preprocessing stage (For details on the wall weight function, see [56,59]).
The particle number density in the ERP model, nERP
, which is used instead of Eq. (9), isi
computed as= ni + z(|xwall− xi |).
nERP
i
i

(20)

Corrections of diﬀerential operators

The ERP model can impose the pressure Neumann and slip/no-slip boundary conditions
in a similar manner to a ghost particle approach [60–64] that does not use ghost particles.
We deﬁne the reﬂectional transformation matrix for particle i across its nearest polygon
inv
whose unit normal vector is nwall
as Rref
i , and for the inverse transformation matrix, R ,
i
as
wall
Rref
⊗ nwalli = I − 2ni
i ,

(21)

Rinv = −I,(22)

and the mirroring position corresponding to particle i, xi , as
− xi ),
xi = xi + 2(xwall
i(23)

where I represents the second-order identity tensor, and xi represents the mirroring
position of xi .
For the no-slip boundary condition, the following diﬀerential operator model for the
Laplacian of the velocity is used instead of Eq. (11).
∇ 2 v ERP= ∇ 2 v i + Rinv ∇ 2 v i ,
i





wall
wall
wall
v i = Rinv v i − 2 v wall
−
n
·
v
n
,
i
i
i
i

(24)

where v wall
represents the velocity of the wall at xwallii .
The following diﬀerential operator model for the pressure gradient is used instead ofEq. (10) to impose the pressure Neumann boundary condition on the wall; this condition
corresponds to the slip boundary condition.
rep∇p ERP= ∇p i + Rrefii ∇p i + f i ,(25)

pi = p i ,

(26)
rep

where the second and third terms are additional terms in the ERP model, and f i represents the repulsive force term for preventing penetrations at curved edges. An earlier
study [56] used the repulsive force
⎧

1 0
⎪
⎪
⎨−α rep 2 l − 1 nwall (0 ≤ r < 1 l 0 )
rep
i
2
r
,
fi =
⎪
⎪
1 0
⎩0
( 2 l ≤ r)which has a parameter α rep .

(27)

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 10 of 31

Structural analysis using the finite element method (FEM)
Weak form of governing equation

In the structural analysis, we determine the displacement, u, and velocity, v, as
u = x − X,

(28)

Du
,
v=
Dt

(29)

where X represents a point in the reference conﬁguration and x represents one in thecurrent conﬁguration. The deformation gradient tensor, F , is deﬁned as the linear operator
between inﬁnitesimal ﬁbers in the reference and the current conﬁgurations.
F=dx= ∇X x.
dX(30)

Based on the weighted residual method with the weight function, W , the weak form of
the governing equation (Eq. ( 2)) in the current conﬁguration, can be written asD(ρv) · W d +
S Dt




S

σ F : (W )d =(W ) =S

ρg · W d +

I


1
∇W + (∇W )T .
2

σ S nS · W d, (31)(32)

Here, we assume that the traction boundary excluding the ﬂuid–structure interface is
imposed on the traction-free conditionσ S nS = 0

on S ,L \ I .

(33)

This study uses the following Saint Venant–Kirchhoﬀ model as a constitutive law ofstructures.
∂W (E),∂EλS
W (E) =(trE)2 + μS tr(E 2 ),
2

1
E=
∇X u + (∇X u)T + (∇X u)(∇X u)T ,
2
S=

(34)
(35)
(36)

where E, W (E), and λS and μS represent the Green–Lagrange strain tensor, strain energy
density function, and Lamé’s constants, respectively.
Spatial discretization based on the Galerkin method

Equation (31) is discretized using the Galerkin method. The weight function and displacement are approximated as
NfnodeW (x)W K NK (x),(37)

uK NK (x),

(38)

K =1
Nfnodeu(x)K =1

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 11 of 31

where W K and uK represent values of the weight function and the displacement at nodeK , respectively; Nfnode represents the number of nodes, and NK (x) represents the shapefunction of node K . We assume that the shape function used here satisﬁes the partitionof the unity condition at any x ∈ S ,NfnodeNK (x) = 1.

(39)

K =1

This study adopts an isoparametric 2nd-order triangle element whose shape function
satisﬁes this condition.
Time integration scheme

Newmark’s β method is used as the time integration scheme in this study [65].

1
un+1 = un + tv n + ( − β)t 2 an + βNB t 2 an+1 ,
2
n+1
n
v
= v + (1 − γNB )tan + γNB tan+1 ,

(40)
(41)

where βNB and γNB are the coeﬃcients. Although the coeﬃcients are decided within
the ranges of the Von-Neumann stability criteria, other parameter sets, which generate
numerical damping, are sometimes used for stability in a positive manner [66–68]. This
study adopts βNB = 0.3025 and γNB = 0.6, which generates numerical damping for stablecomputations.
Nonlinear and linear solvers

The discretized equation given by Eq. (31) with Eqs. (37), (38), (40), and (41) is still
nonlinear. This study uses the Newton–Raphson method for the linearization of thisequation, and the conjugate gradient method as the linear solver.

Compatible interface wave–structure interaction model
We present an FSI model—the compatible interface wave–structure interaction model—
based on the EMPS method for free-surface ﬂow computation and the FEM for structure
computation. The ERP model is adopted to represent the ﬂuid–structure interface as
illustrated in Fig. 3. This model is formulated to exploit the advantages of both methods
when dealing with complicated phenomena including moving boundaries and structureswith various types of constitutive laws. In the proposed model, ﬂuid–structure interfaces
are geometrically compatible because the ERP wall boundary model allows us to use
surfaces of ﬁnite elements directly as wall boundaries in computation of the free-surface
ﬂow. In addition, the force exerted by ﬂuid particles on structures is modeled such thatthe kinetic boundary condition on the ﬂuid–structure interface is satisﬁed.

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 12 of 31

Polygon walls (ERP model)

Fluid particles
Finite elements

Fig. 3 Schematic of the compatible interface wave–structure interaction model between the ﬁnite element

Approximations of the ﬂuid and structure domains
Fluid domain represented by ﬂuid particles in particle methods

The ﬂuid domain at t = 0, 0F ,F , is approximated by ﬂuid particles, each of which has its
own volume and space, as
Nparticle



0

F ,F

0i .

(42)

i=1

However, unlike the mesh-based methods, the particle methods do not have connectivity
information, and therefore, the area of the particle i, i , is not clearly speciﬁed.
Although diﬀerent volumes of particles can be used [69,70], computations using particle
methods assume that all particles have the same volume. This study also adopts thesame assumption. Assuming that ﬂuid particles are set in a uniform grid in the initial
conﬁguration and the initial particle spacing is l 0 , the volume and mass of a particle i are
given asVi = l 0

d

(43)

,

mi = ρVi .

(44)

Wall boundary represented by the polygon walls of the ERP model

The ERP wall boundary model approximates the wall boundary at t = 0, 0F ,F by polygon
walls as
Npol
0
F
,W



s0 .

(45)

i=1

The shape function in the FEM is deﬁned on each polygon and node.
Structural domain represented by ﬁnite elements

When computing a structure based on the FEM, the structural domain at t = 0, 0S , isapproximated by ﬁnite elements as
0S

N
elem

0e ,

i=1

where Nelem represents the number of ﬁnite elements.

(46)

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 13 of 31

Correspondence relation between the polygon walls and ﬁnite elements

In the structural computation, ﬁnite elements represent the ﬂuid–structure interface,S ,I , as their surfaces,
0
F
,I



s̃0 .

(47)

˜ pol
s̃∈

The polygon walls for the ERP model can be obtained by the surface traction of the ﬁnite
elements, which approximates the ﬂuid–structure interface, F ,I , asS0 ,I˜ surfaceb̃∈b̃0 .

(48)

If the same shape functions are used, the polygon walls generated by the surface extraction
are compatible to the ﬁnite elements
I

S ,I

F ,I .

(49)

The distance evaluation between ﬂuid particles and curved polygons requires a high computational cost, and therefore, the polygon walls are approximated as ﬂat triangles in thisstudy even if the original elements are high-order isoparametric elements.
Here, we deﬁne a map, NCR, that represents the node correspondence relation between
˜ pnode andnodes of the polygons and the elements on the ﬂuid–structure interface, L̃ ∈ ˜K̃ ∈ fnode , asK̃ = NCR(L̃).

(50)

Since the number of nodes in both the ﬂuid and structure are the same, NCR is a bijective
mapping. Therefore, the inverse mapping, NCR−1 , exists asL̃ = NCR−1 (K̃ ).

(51)

These node correspondence relation data are necessary for the partitioned coupling analyses conducted in this study.
FSI model from ﬂuid to structure

The EPR model deﬁnes the force exerted on the particle i by walls, f wall
i , as= Vi −∇p ERP+ μ∇ 2 v ERP
f wall
.
i
i
i

(52)

As shown in Fig. 4, we assume that this force is located at the closest point, xwall
i , deﬁned
in Eq. (18), as a point load. The traction force of the ﬂuid can be approximated as
Nparticlewall
δ(x − xwalli )f i ,σ F nF
i=1

where δ represents the Dirac delta function.

(53)

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 14 of 31

equivalent nodal force
point loadiequivalent nodal force

Polygon Walls

Finite Elements

Fig. 4 Point load on polygon wall exerted by particle i and equivalent nodal forces corresponding to the load

The equivalent nodal force exerted from the ﬂuid to the structure at the global node of
the polygon wall, L̃, which corresponds to the second term of right-hand side in Eq. (31),is obtained asσ S nS NNCR(L̃) (x)d.
(54)
f FL̃ →S =
I

Substituting the force boundary condition on the ﬂuid–structure interface (Eq. (6)) into
Eq. (54), the equivalent nodal force can be ﬁnally computed as
f FL̃ →S =−σ F nF NNCR(L̃) (x)dI⎛⎞
Nparticle
wall ⎠−⎝
δ(x − xwall
NNCR(L̃) (x)d
=i )f i
I

(55)

(56)

i=1

Nparticle 

=−
i=1

I

wall
δ(x − xwalli )f i NNCR(L̃) (x)d,(57)

Nparticlewall
f walli NNCR(L̃) (x i ).

=−

(58)

i=1

These forces from the ﬂuid to the structure on each polygon node are projected to thecorresponding global node, K̃ = NCR(L̃), in the FEM as→S= f L̃F →S .
f FNCR(L̃)(59)

FSI model from the structure to the ﬂuid

The contribution from the structure to the ﬂuid can be simply represented by renewing
the position, velocity, and acceleration values at the global node of polygon walls L̃ using
the displacement, velocity, and acceleration values at the corresponding global node ofﬁnite elements, K̃ = NCR(L̃), obtained by the FEM computation as
xNCR−1 (K̃ ) = x0K̃ + uK̃ ,(60)

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 15 of 31

v NCR−1 (K̃ ) = v K̃ ,(61)

aNCR−1 (K̃ ) = aK̃ .

(62)This procedure satisﬁes the boundary condition of the velocity on the ﬂuid–structureinterface given by Eq. (5).
The predictors of the displacement, velocity, and acceleration of polygon walls can be
determined based on the known results of structural computations, and they are usedas the wall boundary condition in the ﬂuid computation. We adopt the acceleration at anode K̃ in the nth step as the predicted acceleration, ân+1
, as
K̃
ân+1
K̃

anK̃ .

(63)

Substituting this into the equations of the time integration based on the Newmark’s beta
and v̂ n+1, can
method, Eqs. (40) and (41), the predicted displacement and velocity, ûn+1
K̃K̃
be expressed as1(t) = unK̃ + tv nK̃ + t 2 anK̃ ,
ûn+1
K̃
2
n
n
v̂ n+1
(t)
=
v
+
ta
.
K̃
K̃
K̃

(64)
(65)

Their equations are second-order interpolation for the displacement, and ﬁrst-order forthe velocity from t n to t n+1 .
Partitioned coupling scheme

Partitioned coupling schemes [55,71–73] solve each physical ﬁeld separately while
exchanging information. Although a smaller time step width is often required for sufﬁcient coupling to be achieved in the case of strong interaction problems, partitioned
approaches have the advantage of software modularity, which allows the use of existingsolvers.
For the ﬂuid and structural computations as abstract functions, F and S , respectively,
these computations in the partitioned coupled analyses can be written as(v  , p ) = F (x ),(x ) = S (v  , p ),

(66)
(67)

where x ,v  , and p represents the position, velocity, and pressure on the ﬂuid–structureinterface, I , respectively. A residual vector whose components are the diﬀerence between
positions before and after a coupling computation is given as
r = x − S (F (x )).

(68)

Partitioned coupling schemes can be categorized into iterative schemes that achievea strong coupling and staggered schemes (so called weak coupling). In analogy with theiterative methods for linear systems, the iterative schemes adopt an iteration loop withina time step that can be selected to treat the entire problem including the coupling eﬀects

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 16 of 31

0.292
0.012
0.08
0.292
0.584
Fig. 5 Dam break problem with an elastic obstacle: initial conﬁguration

Table 2 Dam break with an elastic obstacle: analysis conditions for the ﬂuid computation based on
the EMPS method
Time step width
Number of particles
Particle spacing

2.5 × 10−6 [s]20,000 and 80,000
1.46 × 10−3 and 7.30 × 10−4 [m]Eﬀective radius
Density
Kinetic viscosity

3.1l 0 [m]1.0 × 103 [kg/m3 ]1.0 × 10−6 [m2 /s]Gravitational acceleration
Sound speed coeﬃcient
Repulsive coeﬃcient in the ERP model

10.0[m/s2 ]17.1[m/s]
1.0 × 107 [N/m3 ]implicitly. They skip the iteration loop if the L2 norm of the residual vector in the ith
iteration is less than a small constant .
r i 
<(
r 0 

0).

(69)

Thus far, various iterative partitioned schemes have been proposed; for example, theblock Gauss–Seidel method [74], block Newton method [72,73], and Broyden method[75]. Staggered schemes such as the conventional parallel staggered scheme, CSS scheme,
and improved serial staggered scheme [55] skip the iteration in a coupling step that ignores
the norm of the residual.

Verification and validation
Dam break with an elastic obstacle

The dam break problem with an elastic obstacle is solved to verify the proposed FSI model.
This problem was initially proposed in a previous study [76] using the space-time FEM,
and it has been investigated in other studies [26,33,47,77] for veriﬁcation.
Figure 5 illustrates the initial conﬁguration of the dam break problem with an elastic
obstacle. The elastic obstacle is a linear elastic body that has a width of 0.012[m] anda height of 0.08[m] with a bottom constraint. The water column is 0.146[m] wide and0.292[m] high. The details of the analysis conditions for the ﬂuid and structure are presented respectively in Tables 2 and 3.

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 17 of 31

Table 3 Dam break with an elastic obstacle: analysis conditions for structural computation basedon the FEM
Time step width

2.5 × 10−6 [s]Number of elements
Element type
Young’s modulus

48
Isoparametric 2nd-order
1.0 × 106 [kg/s2 m]

Poisson’s ratio
Density
Gravitational acceleration

0.0
2.5 × 103 [kg/m3 ]10.0[m/s2 ]Table 4 Dam break with an elastic obstacle: details of existing studies
Fluid

Structure

Coupling type

Ryzhakov et al. [26]Walhorn et al. [76]Idelsohn et al. [77]Particle FEM
Stabilized space-time FEM
Particle FEM

FEM
Space-time FEM
FEM

Monolithic scheme
Monolithic scheme
Monolithic scheme

Raﬁee et al. [33]Li et al. [47]

SPH method
SPH-ALE method

SPH method
FEM

Staggered partitioned scheme
Staggered partitioned scheme

Figure 6 shows the displacement of the upper-left corner of the elastic obstacle obtainedin existing studies [26,33,47,76,77]; the numerical methods for ﬂuid and structure, andthe coupling type used in their studies are presented in Table 4. As shown in the ﬁgure,the qualitative behaviors are similar. The elastic obstacle is deformed to the left a little by
the ﬂuid hitting the lower part of the obstacle at approximately t = 0.15[s]. The obstacle
deﬂects to the right because of the force from the ﬂuid, and it reaches the maximum
deﬂection at approximately t = 0.25[s]. In terms of the maximum deﬂection, the existingstudies have a range from 0.04[m] to 0.05[m]. Then, the deﬂected obstacle pushes theﬂuid back while the ﬂuid climbs over the obstacle, and then, the results of the existingstudies become diﬀerent. The existing studies based on mesh-based methods for both
ﬂuid and structure; Ryzhakov et al. [26], Walhorn et al. [76], and Idelsohn et al. [77];exhibit relatively similar trends, although studies in which mesh-free particle methods areused for ﬂuid; Raﬁee et al. [33] and Li et al. [47] show relatively larger deﬂections.
The results obtained by the proposed wave–structure model are shown in Figs. 7, 8, and9; they are similar to the ones in the Ryzhakov’s study [26]. Therefore, we compare ourresults with that in the following discussion.
Figure 7 shows the deﬂection of the upper-left corner of the obstacle, as in Fig. 6, obtainedby the proposed model based on the CSS without predictor. Cases with two spatial resolutions for the ﬂuid computations, 20,000 and 80,000 ﬂuid particles, are compared. Figure 8
shows the deﬂection obtained by the proposed FSI model based on the CSS with the
2nd-order predictor, which is given by Eqs. (64) and (65). Both the results with and without predictor are like the Ryzhakov’s result though the peaks are slightly higher thanthat. There are a few diﬀerences when the spatial resolution is increased in terms of the
maximum deﬂections; the results of the proposed model exhibit small ﬂuctuations fromt = 0.5[s] to 0.6[s]. This is caused by a collision of a splash of ﬂuid to the obstacle after
the ﬂuid impacts the right wall; this may be very sensitive to the free-surface and wallboundary treatment.
Figure 9 shows the deﬂection obtained by the proposed FSI model based on the block
Gauss-Seidel method in cases with 20,000 and 80,000 ﬂuid particles. The snapshots of theN. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 18 of 31

0.05

Ryzhakov et al.
Walhorn et al.
Idelsohn et al.
Rafiee et al.
Li et al.

0.04

Deflection [m]

0.03

0.02

0.01

0

-0.01

-0.02
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time [s]Fig. 6 Dam break with an elastic obstacle: deﬂection of the upper-left corner of the obstacle obtained in existingstudies ([26], [76], [77], [33], and [47])case with 80,000 particles are shown in Fig. 11. The convergence criterion for the iterationin Eq. (69) is set as  = 1.0×10−4 . A Strong coupling in this case is like the nonlinear solver
for the boundary nonlinearity of the surfaces of the structure. Although higher accuracy
and stability is expected by adopting the iteration, there are slight diﬀerences between
the cases with the block Gauss–Seidel (BGS) method and CSS. The result indicate that
this problem may not have a large coupling eﬀect, and the time step width is suﬃciently
small for weak coupling because that used for this computation is limited by the ﬂuid
computation with the fully explicit algorithm.
Figure 10 shows the number of iterations averaged every 200 steps in the simulation
with the BGS method. We limit the number of iterations to 200; there are 200,000 steps inthis simulation. Although the number of iterations is less than 10 in most steps, 210 outof all calculation steps do not reach our convergence criteria, and a rapid increase in thenumber of iterations is observed between t = 0.45[s] to 0.6[s]. This result coincides withthe small ﬂuctuation shown in Figs. 7 and 8, and it is believed to be associated with theisolated droplets of the free-surface particles colliding with the obstacle. In such scenarios,
the repulsive force can be dominant, and the initial residual norm can be very small, andthis can result in poor convergence properties.
Sloshing tank with an elastic beam

Iglesias et al. [78,79] performed a series of experiments to simulate the free-surface ﬂow
interacting with deformable structures as a SPHERIC benchmark test case [80]. Further,
they showed the results of the simulations based on PFEM [79]. Yang et al. [49] computed
one of the experiments using SPH for free-surface ﬂow and FEM for structure.
Figure 12 illustrates the geometric speciﬁcations of the experiment. The length, height,
and width of the tank are 0.609[m], 0.3445[m], and 0.039[m], respectively. An elastic beam

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

0.05

Page 19 of 31

Ryzhakov et al.
Proposed (CSS without predictor; 100 x 200)
Proposed (CSS without predictor; 200 x 400)

0.04

0.03

Deflection [m]

0.02

0.01

0

-0.01

-0.02

-0.03
0

0.1

0.2

0.3

0.4
0.5
Time [s]

0.6

0.7

0.8

0.9

Fig. 7 Dam break with an elastic obstacle: deﬂection of the upper-left corner of the obstacle obtained by the
proposed model based on the conventional serial staggered (CSS) scheme without a predictor
0.05

Ryzhakov et al.
Proposed (CSS with predictor; 100 x 200)
Proposed (CSS with predictor; 200 x 400)

0.04

0.03

Deflection [m]

0.02

0.01

0

-0.01

-0.02

-0.03
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time [s]Fig. 8 Dam break with an elastic obstacle: deﬂection of the upper-left corner of the obstacle obtained by the
proposed model based on the conventional serial staggered (CSS) scheme with a second-order predictor

is initially located at the center of the tank, and it is immersed in sunﬂower oil. There aregaps of 0.0029[m] between the beam and the tank. The tank moves in an oscillatory
manner around the center point at the bottom, as shown in Fig. 13. A local coordinate
system moving along the tank is deﬁned, and deﬂection at the top of the beam is measured
in the local coordinate system. Figure 14 shows the time history of the roll angle provided
in the SPHERIC benchmark website [80].

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 20 of 31

0.05

Ryzhakov et al.
Proposed (BGS; 100 x 200)
Proposed (BGS; 200 x 400)

0.04

0.03

Deflection [m]

0.02

0.01

0

-0.01

-0.02

-0.03
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time [s]Fig. 9 Dam break with an elastic obstacle: deﬂection of the upper-left corner of the obstacle obtained by the
proposed model based on the block Gauss–Seidel (BGS) method

The average number of iterations

20

BGS; 200 x 400

15

10

5

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time [s]Fig. 10 Dam break with an elastic obstacle: number of iterations in the simulation of the proposed FSI model
averaged every 200 steps in the case with the iterative coupling scheme

We conduct two-dimensional simulations of the proposed FSI model to reproduce theexperiment with the initial conﬁguration shown in Fig. 12 and we use the roll angle shownin Fig. 14. Tables 5 and 6 present the analysis conditions for the ﬂuid and structural
computations. Further, we use the same physical values of the oil and the beam (density,
kinetic viscosity, and Young’s ratio) as that provided in the report [80].

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

t = 0.0[s]t = 0.12[s]t = 0.24[s]t = 0.36[s]t = 0.48[s]t = 0.60[s]t = 0.72[s]t = 0.84[s]Fig. 11 Dam break with an elastic obstacle: snapshots of the simulation conducted by the proposed model
based on the block Gauss–Seidel (BGS) method from t = 0.0[s] to t = 0.84[s]; pressure contours are indicated bycolor (min: 0.0[N/m2 ], max: 2500.0[N/m2 ])In the following considerations, we use the result in the case with a time spacing of
1.0 × 10−4 [s], initial particle spacing of 2.5 × 10−3 [m], predictor, and CSS scheme.
Figure 15 shows the displacement in the x direction at the top of the beam observed in thelocal coordinate system as computed by the proposed FSI model. The result is comparedwith that of the experiment, and with the two-dimensional numerical results of the PFEM
(Souto-Igresias et al. [79]) and SPH-FEM (Yang et al. [49]). The PFEM computation isconducted from 0.0[s] to 5.0[s]. Although the result of the proposed model overestimates

Page 21 of 31

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 22 of 31

Table 5 Sloshing tank with an elastic beam: analysis conditions for the ﬂuid computation
Time step width
Number of particles

1.0 × 10−4 and 5.0 × 10−5 [s]44,528 and 278,300

Particle spacing
Eﬀective radius
Density

2.5 × 10−3 and 5.0 × 10−4 [m]3.1l 0 [m]9.0 × 102 [kg/m3 ]Kinetic viscosity
Gravitational acceleration
Sound speed coeﬃcient

5.0 × 10−5 [m2 /s]
9.81[m/s2 ]10.59[m/s]Repulsive coeﬃcient in the ERP model

1.0 × 107 [N/m3 ]Table 6 Sloshing tank with an elastic beam: analysis conditions for the structural computation
Time step width

1.0 × 10−4 and 5.0 × 10−5 [s]Number of elements
Element type
Young’s modulus

200
Isoparametric 2nd-order
3.3 × 106 [kg/s2 m]

Poisson’s ratio
Density
Gravitational acceleration

0.49
1.1 × 103 [kg/m3 ]9.81[m/s2 ]the peak displacements, it agrees well with the result of the PFEM. We will discuss theoverestimation later.
We conduct sensitivity tests in terms of the eﬀect of the predictor, time spacing, andspatial resolution. Figure 16 shows the results with and without the predictor, and Fig. 17
shows the ones with t = 1.0 × 10−4 [s] and 5.0 × 10−5 [s]. There are almost no diﬀerence
between the results, and therefore, in this test case, these eﬀects can be considered trivial.
Figure 18 shows the diﬀerence between the cases with l 0 = 2.5 × 10−3 [m] (44,528 ﬂuid
particles) and l 0 = 5.0 × 10−4 [m] (278,300 ﬂuid particles). The peak displacements in thehigher spatial resolution are slightly small compared with those of the coarser ones.
Figure 19 shows the results obtained using the proposed FSI model with the CSS coupling and iterative partitioned coupling scheme with the BGS method. The number of
iterations averaged every 100 steps in the case with the BGS method is presented inFig. 20. We set the small constant used as a convergence criterion in equation (69) to
1.0 × 10−4 . The snapshots of that simulation are presented in Fig. 21. The result with the
BGS method exhibits smaller peak displacements compared to that with the CSS scheme.
A comparison of the second positive peak displacement of the numerical results of the
proposed model with the experiment results indicates diﬀerences of about 13.0% in thecase with the BGS method and 25.3% with the CSS scheme. Though the result provided
by Yang et al. are in good agreement with the experiment results, as indicated in Fig. 15,
the overestimations seen in the results obtained by the proposed model and the PFEM
seem reasonable. This is because the experimental setting has gaps of 0.0029[m] between
the tank and the beam, as illustrated in Fig. 12. The area without the gaps is approximately
14.9% larger than the one with the gaps. Focusing on the pressure diﬀerence between
the left and right regions, displacement of the structure is proportional to surface area inthe case of a linear elastic body. Although it is diﬃcult to discuss rigorously because of
the nonlinear deformation, overestimations can occur because two-dimensional simulations implicitly assume no gap. In addition, the two-dimensional simulations do not have

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

0.039

Page 23 of 31

0.609

z

z

y

x

0.34450.34450.11480.00290.03320.30250.004

Fig. 12 Sloshing tank with an elastic beam: initial conﬁgurationyx

Fig. 13 Sloshing tank with an elastic beam: coordinate system for the tankﬂuid ﬂowing in and out between the left and right regions through the gap compared tothe experiment. Therefore, the water level diﬀerence is expected to be more pronounced
than in the experiment, and the pressure diﬀerence due to the water level diﬀerence isalso expected to peak later. The displacement of the elastic body is expected to follow
this peak. It also contributes to overestimations of displacement in two-dimensional simulations and is thought to be the cause of the phase diﬀerence from the experimental
values as seen in Fig. 19. However, there are also three-dimensional eﬀects or dynamic
pressures caused by ﬂow through the gaps; these three-dimensional computations should
be conducted as part of a future work for further discussion of this test case.

Conclusions
This study presented a compatible interface wave–structure interaction model based ona mesh-free particle method for free-surface ﬂow analysis, and the FEM for structural
analysis.
The EMPS method is employed for free-surface ﬂow analysis. We adopted the ERP
wall boundary model [56] in the EMPS method. The proposed model was geometricallyN. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 24 of 31

Experiment4
32

Angle [deg]

1
0
-1
-2
-3
-4
0

1

2

3

4

5

6

7

8

9

Time [s]Fig. 14 Sloshing tank with an elastic beam: roll angle in the experiment
Experiment
Iglesias et al.
Yang et al.
Proposed (CSS with predictor)

0.1

Deflection [m]

0.05

0

-0.05

-0.1
0

1

2

3

4

5

6

7

8

9

Time [s]Fig. 15 Sloshing tank with an elastic beam: displacements of the x direction at the top of the beam obtained by
the proposed FSI model compared to experiments and numerical results in the existing studies ([79] and [49])compatible on the interfaces because the ERP model could directly use FE surfaces andtheir shape functions on ﬂuid–structure interfaces.
The force exerted by ﬂuid particles on structures was modeled as point loads using
the Dirac delta function based on the compatible interface between the polygons and theﬁnite elements. This procedure was relatively simple and satisﬁed the kinetic boundary
condition on the ﬂuid–structure interface. We formulated a predictor for the displacement and velocity to verify its eﬀectiveness when used with the staggered (weak coupling)N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 25 of 31

Experiment
Proposed (CSS with predictor)
Proposed (CSS without predictor)

0.1

Deflection [m]

0.05

0

-0.05

-0.1
0

1

2

3

4

5

Time [s]Fig. 16 Sloshing tank with an elastic beam: displacements of the x direction at the top of the beam obtained by
the proposed FSI model with and without the second-order predictor
Experiment
Proposed (CSS with predictor; Δ t=0.0001)
Proposed (CSS with predictor; Δ t=0.00005)

0.1

Deflection [m]

0.05

0

-0.05

-0.1
0

1

2

3

4

5

Time [s]Fig. 17 Sloshing tank with an elastic beam: displacements of the x direction at the top of the beam obtained by
the proposed FSI model with diﬀerent time spacings

schemes. Further, we developed a strong coupling algorithm based on an iterative partitioned scheme with the BGS method.
We performed veriﬁcation and validation tests and sensitivity tests for the eﬀect of thepredictor, time spacing, and spatial resolution by solving two benchmark problems. The
results obtained for the dam break problem with an elastic obstacle using the proposed
wave–structure model exhibited good agreement with the ones reported in Ryzhakov’s

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 26 of 31

Experiment
Proposed (CSS with predictor; l0=0.0025)
Proposed (CSS with predictor; l0=0.0005)

0.1

Deflection [m]

0.05

0

-0.05

-0.1
0

1

2

3

4

5

Time [s]Fig. 18 Sloshing tank with an elastic beam: displacements of x direction at the top of the beam obtained by the
proposed FSI model with diﬀerent spatial resolutions

Experiment
Proposed (CSS with predictor)
Proposed (BGS)

0.1

Deflection [m]

0.05

0

-0.05

-0.1
0

1

2

3

4

5

Time [s]Fig. 19 Sloshing tank with an elastic beam: displacements of the x direction at the top of the beam obtained by
the proposed FSI model with the CSS scheme and iterative coupling scheme based on the block Gauss–Seidel(BGS) method

study [26], which adopts a strong coupling scheme and the particle FEM. In this test
case, there are little diﬀerence in the cases with the BGS method (strong coupling) andCSS method (weak coupling). The result indicate that this problem may not have a large
coupling eﬀect, and the time step width is suﬃciently small for weak coupling because
the time step width used for this computation is limited by the ﬂuid computation with

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Page 27 of 31

10

Proposed (BGS)

The average number of iterations

8

6

4

2

0
0

1

2

3

4

5

Time [s]Fig. 20 Sloshing tank with an elastic beam: number of iterations in the simulation of the proposed FSI model
averaged every 100 steps in the case with the iterative coupling schemet = 1.84 [s]t = 2.12 [s]t = 2.32 [s]t = 2.56 [s]0.0

1200.0Pressure [N/m2]

Fig. 21 Sloshing tank with an elastic beam: snapshots of the simulation conducted by the proposed model
based on the block Gauss–Seidel method at t = 1.84, 2.12, 2.32, and 2.56[s]; pressure contours are indicated bycolor (min: 0.0[N/m2 ], max: 1200.0[N/m2 ])N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

the fully explicit algorithm. However, this result is valuable because we cannot discuss thecoupling eﬀect without attempting the strong coupling scheme. Next, we simulated anexperiment of a sloshing tank with an elastic beam. Although the results of the proposedmodel overestimated the peak displacements compared with those of the experiment, it
agreed well with the result of the existing study that uses particle FEM (Souto-Igresias etal. [79]). The experimental setting has gaps between the tank and the beam. This study
conducted two-dimensional simulations, which implicitly assumed no gap, and therefore,
overestimations were considered to occur.
Incorporating the three-dimensional and parallel simulator for the EMPS method and
ERP wall boundary model [59] to achieve three-dimensional wave–structure interaction
analysis will be considered as a future task for practical applications. The 3D simulator
would also allow for a more quantitative discussion of the issues in “Sloshing tank with anelastic beam”.
Acknowledgements
This work was supported by a JSPS Grant-in-Aid for Scientiﬁc Research (B) (Grant Number 22H03601) and JST FOREST
Program (Grant Number JPMJFR215S).
Author contributions
NM developed the proposed model in this study. NM conducted the veriﬁcation and validation test cases andinterpreted the results. NM wrote this manuscript.
Funding
Japan Society for the Promotion of Science, Grant-in-Aid for Scientiﬁc Research (B), 22H03601. Japan Science and theTechnology Agency, Fusion Oriented REsearch for disruptive Science and Technology (FOREST), JPMJFR215S.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable
request.

Declarations
Competing interests
They have no competing interests.
Received: 1 December 2022 Accepted: 28 May 2023

References
1. Asakura R, Iwase K, Ikeya T, Takao M, Kaneto T, Fujii N. et al. Experimental study on wave force acting on onshore structures due to overﬂowing tsunamis. In: Proceedings of Coastal Engineering. vol. 47. Japan Society of Civil
Engineers; 2000;. p. 911–915. (in Japanese).
2. Asakura R, Iwase K, Ikeya T, Takao M, Kaneto T, Fujii N. et al. The tsunami wave force acting on land structures. In:
Proceedings of 28th International Conference on Coastal Engineering; 2002; p. 1191–1202.
3. Ikeno M, Matsuyama M, Tanaka H. Shoaling and soliton ﬁssion of tsunami on a shelf and wave pressure for tsunamiresistant design of breakwater by large wave ﬂume-experiments. In: Proceedings of Coastal Engineering. vol. 45; 1998;
p. 366–370. (in Japanese).
4. Ikeno M, Mori N, Tanaka H. Experimental study on drift behavior and its impulse force by tsunami. Annual J Coast Eng.
2001;48:846–50 (in Japanese).
5. Matsutomi H. A practical formula for estimating impulsive force due to driftwoods and variation features of theimpulsive forces. J Hydraulic Coast Environ Eng. 1999;621:111–27 (in Japanese).
6. Grilli ST, Harris JC, Bakhsh TST, Masterlark TL, Kyriakopoulos C, Kirby JT, et al. Numerical simulation of the 2011 Tohoku
tsunami based on a new transient FEM co-seismic source: comparison to far-and near-ﬁeld observations. Pure Appl
Geophys. 2013;170(6–8):1333–59.
7. Oishi Y, Imamura F, Sugawara D. Near-ﬁeld tsunami inundation forecast using the parallel TUNAMI-N2 model: application to the 2011 Tohoku-Oki earthquake combined with source inversions. Geophys Res Lett. 2015;42(4):1083–91.
8. Westerink JJ, Luettich RA, Feyen JC, Atkinson JH, Dawson C, Roberts HJ, et al. A basin-to channel-scale unstructured
grid hurricane storm surge model applied to southern Louisiana. Mon Weather Rev. 2008;136(3):833–64.
9. Dao M, Tkalich P. Tsunami propagation modelling? A sensitivity study. Natural Hazards Earth Syst Sci. 2007;7(6):741–54.
10. Hamzah MA, Mase H, Takayama T. Simulation and experiment of hydrodynamic pressure on a tsunami barrier. In:
Proceedings of Coastal Engineering. vol. 47; 2000; p. 371–375 (in Japanese).

Page 28 of 31

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

11. Murotani K, Koshizuka S, Tamai T, Shibata K, Mitsume N, Shinobu Y, et al. Development of hierarchical domain
decomposition explicit MPS method and application to large-scale tsunami analysis with ﬂoating objects. J Adv
Simul Sci Eng. 2014;1(1):16–35.
12. Yoneyama N, Nagashima H. Development of a three dimensional numerical analysis method for the drift behavior intsunami. J Jpn Soc Civil Eng Ser B2 (Coast Eng). 2009;65(1):266–70.
13. Pringle WJ, Yoneyama N, Mori N. Two-way coupled long wave-RANS model: solitary wave transformation and breaking
on a plane beach. Coast Eng. 2016;114:99–118.
14. Sitanggang K, Lynett P. Multi-scale simulation with a hybrid Boussinesq-RANS hydrodynamic model. Int J Numer
Meth Fluids. 2010;62(9):1013–46.
15. Takase S, Moriguchi S, Terada K, Kato J, Kyoya T, Kashiyama K, et al. 2D–3D hybrid stabilized ﬁnite element method fortsunami runup simulations. Comput Mech. 2016;58(3):411–22.
16. Mitsume N, Donahue AS, Westerink JJ, Yoshimura S. Coupling methods between ﬁnite element-based Boussinesqtype wave and particle-based free-surface ﬂow models. Int J Numer Meth Fluids. 2018;88(3):141–68.
17. Sato S, Okayasu A, Yeh H, Fritz HM, Tajima Y, Shimozono T. Delayed survey of the 2011 Tohoku Tsunami in the former
exclusion zone in Minami-Soma, Fukushima Prefecture. Pure Appl Geophys. 2014;171(12):3229–40.
18. Yoshimura S, Nakamura H, Kanayama H, Aoki T, Himeno T, Sakai Y. et al. Understanding the mechanism of tsunamiinduced damage to machines and structures based on a discipline of mechanics. In: Report of the JSME ResearchCommittee on the Great East Japan Earthquake Disaster. The Japan Society of Mechanical Engineers; 2014; p. 24–57.
19. Hirt CW, Nichols BD. Volume of ﬂuid (VOF) method for the dynamics of free boundaries. J Comput Phys. 1981;39(1):201–25.
20. Fedkiw SOR. Level set methods and dynamic implicit surfaces. 2003.
21. Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi
formulations. J Comput Phys. 1988;79(1):12–49.
22. Idelsohn S, Oñate E, Pin FD, Calvo N. Fluid-structure interaction using the particle ﬁnite element method. Comput
Methods Appl Mech Eng. 2006;195(17):2100–23.
23. Idelsohn SR, Oñate E, Pin FD. The particle ﬁnite element method: a powerful tool to solve incompressible ﬂows with
free-surfaces and breaking waves. Int J Numer Meth Eng. 2004;61(7):964–89.
24. Oñate E, Celigueta MA, Idelsohn SR, Salazar F, Suárez B. Possibilities of the particle ﬁnite element method for ﬂuidsoil-structure interaction problems. Comput Mech. 2011;48(3):307–18.
25. Oñate E, Idelsohn SR, Del Pin F, Aubry R. The particle ﬁnite element method—an overview. Int J Comput Methods.
2004;1(02):267–307.
26. Ryzhakov P, Rossi R, Idelsohn S, Oñate E. A monolithic Lagrangian approach for ﬂuid–structure interaction problems.
Comput Mech. 2010;46(6):883–99.
27. Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not
R Astron Soc. 1977;181(3):375–89.
28. Lucy LB. A numerical approach to the testing of the ﬁssion hypothesis. Astron J. 1977;82:1013–24.
29. Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible ﬂuid. Nucl Sci Eng.
1996;123(3):421–34.
30. Antoci C, Gallati M, Sibilla S. Numerical simulation of ﬂuid–structure interaction by SPH. Comput Struct.
2007;85(11):879–90.
31. Chikazawa Y, Koshizuka S, Oka Y. A particle method for elastic and visco-plastic structures and ﬂuid–structure interactions. Comput Mech. 2001;27(2):97–106.
32. Hwang SC, Khayyer A, Gotoh H, Park JC. Development of a fully Lagrangian MPS-based coupled method for simulationof ﬂuid–structure interaction problems. J Fluids Struct. 2014;50:497–511.
33. Raﬁee A, Thiagarajan KP. An SPH projection method for simulating ﬂuid-hypoelastic structure interaction. Comput
Methods Appl Mech Eng. 2009;198(33):2785–95.
34. Aktay L, Johnson AF. FEM/SPH coupling technique for high velocity impact simulations. In: Advances in Meshfree
Techniques. Springer; 2007; p. 147–167.
35. Attaway SW, Heinstein MW, Swegle JW. Coupling of smooth particle hydrodynamics with the ﬁnite element method.
Nucl Eng Des. 1994;150(2):199–205.
36. De Vuyst T, Vignjevic R, Campbell J. Coupling between meshless and ﬁnite element methods. Int J Impact Eng.
2005;31(8):1054–64.
37. Gu Y, Zhang L. Coupling of the meshfree and ﬁnite element methods for determination of the crack tip ﬁelds. Eng
Fract Mech. 2008;75(5):986–1004.
38. Johnson GR. Linking of Lagrangian particle methods to standard ﬁnite element methods for high velocity impact
computations. Nucl Eng Des. 1994;150(2):265–74.
39. Lu Y, Wang Z, Chong K. A comparative study of buried structure in soil subjected to blast load using 2D and 3D
numerical simulations. Soil Dyn Earthq Eng. 2005;25(4):275–88.
40. Anghileri M, Castelletti LM, Tirelli M. Fluid–structure interaction of water ﬁlled tanks during the impact with theground. Int J Impact Eng. 2005;31(3):235–54.
41. Comas-Cardona S, Groenenboom P, Binetruy C, Krawczak P. A generic mixed FE-SPH method to address hydromechanical coupling in liquid composite moulding processes. Compos A Appl Sci Manuf. 2005;36(7):1004–10.
42. Fourey G, Oger G, Le Touzé D, Alessandrini B. violent ﬂuid–structure interaction simulations using a coupled SPH/FEM
method. In: IOP Conference Series: Materials Science and Engineering. vol. 10. IOP Publishing; 2010; p. 012041.
43. Groenenboom PH, Cartwright BK. Hydrodynamics and ﬂuid–structure interaction by coupled SPH-FE method. J
Hydraul Res. 2010;48(S1):61–73.
44. Hu D, Long T, Xiao Y, Han X, Gu Y. Fluid–structure interaction analysis by coupled FE-SPH model based on a novel
searching algorithm. Comput Methods Appl Mech Eng. 2014;276:266–86.
45. Lee CJK, Noguchi H, Koshizuka S. Fluid–shell structure interaction analysis by coupled particle and ﬁnite element
method. Comput Struct. 2007;85(11):688–97.

Page 29 of 31

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

46. Li Z, Leduc J, Combescure A, Leboeuf F. Coupling of SPH-ALE method and ﬁnite element method for transientﬂuid–structure interaction. Comput Fluids. 2014;103:6–17.
47. Li Z, Leduc J, Nunez-Ramirez J, Combescure A, Marongiu JC. A non-intrusive partitioned approach to couple smoothed
particle hydrodynamics and ﬁnite element methods for transient ﬂuid-structure interaction problems with large
interface motion. Comput Mech. 2015; 1–22.
48. Thiyahuddin M, Gu Y, Gover R, Thambiratnam D. Fluid–structure interaction analysis of full scale vehicle-barrier impact
using coupled SPH-FEA. Eng Anal Boundary Elem. 2014;42:26–36.
49. Yang Q, Jones V, McCue L. Free-surface ﬂow interactions with deformable structures using an SPH-FEM model. Ocean
Eng. 2012;2012(55):136–47.
50. Yang Q. SPH Simulation of ﬂuid–structure interaction problems with application to hovercraft. Virginia Polytechnic
Institute and State University; 2012;
51. Vila J. On particle weighted methods and smooth particle hydrodynamics. Math Models Methods Appl Sci.
1999;9(02):161–209.
52. Mitsume N, Yoshimura S, Murotani K, Yamada T. Improved MPS-FE ﬂuid–structure interaction coupled method with
MPS polygon wall boundary model. Comput Model Eng Sci. 2014;4:229–47.
53. Mitsume N, Yoshimura S, Murotani K, Yamada T. MPS-FEM partitioned coupling approach for ﬂuid–structure interaction with free surface ﬂow. Int J Couput Methods. 2014;11(4):1350101.
54. Mitsume N, Yoshimura S, Murotani K, Yamada T. Inundation simulation coupling free surface ﬂow and structures.
In: High-performance computing for structural mechanics and earthquake/tsunami engineering. Springer; 2016; p.
179–199.
55. Farhat C, Lesoinne M. Two eﬃcient staggered algorithms for the serial and parallel solution of three-dimensional
nonlinear transient aeroelastic problems. Comput Methods Appl Mech Eng. 2000;182(3):499–515.
56. Mitsume N, Yoshimura S, Murotani K, Yamada T. Explicitly represented polygon wall boundary model for explicit MPS
method. Comput Particle Mech. 2015;2(1):73–89.
57. Shakibaeinia A, Jin YC. A weakly compressible MPS method for modeling of open-boundary free-surface ﬂow. Int J
Numer Meth Fluids. 2010;63(10):1208–32.
58. Oochi M, Koshizuka S, Sakai M. Explicit MPS algorithm for free surface ﬂow analysis. transactions of the Japan Society
for Computational Engineering and Science. 2010;(20120013) (in Japanese).
59. Mitsume N, Yoshimura S, Yamada T. Parallel analysis system for free-surface ﬂow using MPS with explicitly represented
polygon wall boundary model. Comput Particle Mech. 2020;7(2):279–90.
60. Colagrossi A, Landrini M. Numerical simulation of interfacial ﬂows by smoothed particle hydrodynamics. J Comput
Phys. 2003;191(2):448–75.
61. Delorme L, Colagrossi A, Souto-Iglesias A, Zamora-Rodriguez R, Botia-Vera E. A set of canonical problems in sloshing,Part I: pressure ﬁeld in forced roll-comparison between experimental results and SPH. Ocean Eng. 2009;36(2):168–78.
62. Libersky LD, Petschek AG, Carney TC, Hipp JR, Allahdadi FA. High strain Lagrangian hydrodynamics: a threedimensional SPH code for dynamic material response. J Comput Phys. 1993;109(1):67–75.
63. Morris JP, Fox PJ, Zhu Y. Modeling low Reynolds number incompressible ﬂows using SPH. J Comput Phys.
1997;136(1):214–26.
64. Oger G, Doring M, Alessandrini B, Ferrant P. Two-dimensional SPH simulations of wedge water entries. J Comput Phys.
2006;213(2):803–22.
65. Newmark NM. A method of computation for structural dynamics. J Eng Mech Div. 1959;85(EM3):67–94.
66. Ishihara D, Kanei S, Yoshimura S, Horie T. Eﬃcient parallel analysis of shell-ﬂuid interaction problem by using monolithic
method based on consistent pressure Poisson equation. J Comput Sci Technol. 2008;2(1):185–96.
67. Minami S, Kawai H, Yoshimura S. Parallel BDD-based monolithic approach for acoustic ﬂuid–structure interaction.
Comput Mech. 2012;50(6):707–18.
68. Yamada T, Yoshimura S. Line search partitioned approach for ﬂuid–structure interaction analysis of ﬂapping wing.
Comput Model Eng Sci. 2008;24(1):51.
69. Vacondio R, Rogers B, Stansby P. Accurate particle splitting for smoothed particle hydrodynamics in shallow waterwith shock capturing. Int J Numer Meth Fluids. 2012;69(8):1377–410.
70. Vacondio R, Rogers B, Stansby P, Mignosa P, Feldman J. Variable resolution for SPH: a dynamic particle coalescing andsplitting scheme. Comput Methods Appl Mech Eng. 2013;256:132–48.
71. Felippa CA, Park K, Farhat C. Partitioned analysis of coupled mechanical systems. Comput Methods Appl Mech Eng.
2001;190(24):3247–70.
72. Matthies HG, Niekamp R, Steindorf J. Algorithms for strong coupling procedures. Comput Methods Appl Mech Eng.
2006;195(17):2028–49.
73. Matthies HG, Steindorf J. Partitioned strong coupling algorithms for ﬂuid–structure interaction. Comput Struct.
2003;81(8):805–12.
74. Le Tallec P, Mouro J. Fluid structure interaction with large structural displacements. Comput Methods Appl Mech Eng.
2001;190(24):3039–67.
75. Minami S, Yoshimura S. Performance evaluation of nonlinear algorithms with line-search for partitioned couplingtechniques for ﬂuid–structure interactions. Int J Numer Meth Fluids. 2010;64(10–12):1129–47.
76. Walhorn E, Kölke A, Hübner B, Dinkler D. Fluid–structure coupling within a monolithic model involving free surface
ﬂows. Comput Struct. 2005;83(25):2100–11.
77. Idelsohn SR, Marti J, Limache A, Oñate E. Uniﬁed Lagrangian formulation for elastic solids and incompressible ﬂuids: application to ﬂuid–structure interaction problems via the PFEM. Comput Methods Appl Mech Eng.
2008;197(19):1762–76.
78. Botia Vera E, Souto Iglesias A, Bulian G, Lobovskỳ L. Three SPH Novel Benchmark Test Cases for free surface ﬂows. In:
Proceedings of the 5th ERCOFTAC SPHERIC workshop on SPH applications; 2010.
79. Iglesias AS, Idelsohn S, Martí J, Rodriguez RZ, Ibáñez EO. Modeling of free surface ﬂows with elastic bodies interactions.
In: Proceedings of the 27th Symposium on Naval Hydrodynamics; 2008; p. 1–17.
80. Spheric benchmarks;. http://canal.etsin.upm.es/ftp/SPHERIC_BENCHMARKS/, Accessed 23 Nov 2022.

Page 30 of 31

N. Mitsume Advanced Modeling and Simulation in Engineering Sciences(2023)10:11

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional aﬃliations.

Page 31 of 31