Computers and Structures 85 (2007) 879–890

www.elsevier.com/locate/compstruc

Numerical simulation of fluid–structure interaction by SPH

Carla Antoci *, Mario Gallati, Stefano Sibilla
Dipartimento di Ingegneria Idraulica e Ambientale, Università degli Studi di Pavia, via Ferrata, 1, 27100 Pavia, Italy

Received 31 May 2006; accepted 2 January 2007

Available online 15 February 2007

Abstract

A Lagrangian model for the numerical simulation of fluid–structure interaction problems is proposed in the present paper. In the
method both fluid and solid phases are described by smoothing particle hydrodynamics: fluid dynamics is studied in the inviscid approx-
imation, while solid dynamics is simulated through an incremental hypoelastic relation. The interface condition between fluid and solid is
enforced by a suitable term, obtained by an approximate SPH evaluation of a surface integral of fluid pressure.
The method is validated by comparing numerical results with laboratory experiments where an elastic plate is deformed under the
effect of a rapidly varying fluid flow.
 2007 Elsevier Ltd. All rights reserved.

Keywords: Fluid–structure interaction; Smoothed particle hydrodynamics; Interface conditions

1. Introduction ture are important). Structures are usually described by

Lagrangian formulations, whereas fluids are often
In many engineering applications, the forces exerted by described by Eulerian formulations. The coupling of the
a fluid flow on the confining solid boundaries do not mod- two media is usually obtained by an Arbitrary-Lagrang-
ify significantly the geometry of the boundaries. In this ian–Eulerian (ALE) formulation for the fluid. Significant
cases, the fluid flow can be studied as occurring within rigid contributions [1–3] have been proposed in the simulation
boundaries, and the forces applied on the solid boundaries of FSI problems in this context. Rugonyi and Bathe [1] per-
can be obtained after the characteristics of the fluid motion form a simplified stability analysis of the interface equa-
have been determined. tions and study the long-term dynamic stability of FSI
On the other hand, whenever the characteristic times of systems by use of Lyapunov characteristic exponents. They
the motion of the fluid flow and of the solid boundaries are also show the solution of some FSI problems, as the
comparable, it is necessary to couple the dynamics of the dynamics of spring-loaded valves in fuel pumps, that indi-
two media. These fluid–structure interaction (FSI) prob- cate the actual possibility to simulate complex coupled phe-
lems can be solved by employing either a simultaneous nomena. Recent developments in the simulation of viscous
(or direct) solution or a partitioned (or iterative) solution. incompressible and compressible fluid flows with structural
A description of the two procedures can be found in [1], interactions are discussed in [2]. Le Tallec and Mouro [3]
together with the explanation of their main advantages simulate the dynamics of an hydroelastic shock absorber
and drawbacks. The simultaneous technique is particularly adopting an ALE formulation for the fluid equations.
convenient when the interaction between the structure and An alternative approach to the numerical simulation of
the fluid is very strong (and the displacements of the struc- FSI problems consists in the description of both the fluid
and the structure motion by a Lagrangian formulation.
This can be especially effective when studying problems
*
Corresponding author. Tel.: +39 0382 985321; fax: +39 0382 985589. characterized by large displacements of the fluid–structure
E-mail address: [email protected] (C. Antoci). interface and by a rapidly moving fluid free-surface. An

0045-7949/$ - see front matter  2007 Elsevier Ltd. All rights reserved.
doi:10.1016/j.compstruc.2007.01.002
880 C. Antoci et al. / Computers and Structures 85 (2007) 879–890

example of these problems is the FSI inside safety valves tensor and the notation implies summation over repeated
for pressure reduction, where an elastic plate deforms indices.
owing to water pressure, allowing part of the fluid to flow The stress tensor can be decomposed into its isotropic
out at atmospheric conditions, thus causing a pressure and deviatoric parts:
relief in the connected pipe. In this kind of problems, the rij ¼ pdij þ S ij ; ð3Þ
use of Lagrangian techniques for both the solid and the
fluid part of the problem appears promising, as it permits where p ¼ rkk =3 is pressure, Sij is the deviatoric stress ten-
to easily follow in time the motion of the fluid–solid inter- sor and dij is the Kronecker tensor.
face and to simulate the free-surface of the fluid without Pressure can be formally defined in the same way for
any specific treatment. In particular, encouraging results both fluid and solid by the following linearized equation
have been recently obtained by the smoothed particle of state, which holds for small variations of density:
hydrodynamics (SPH) technique (see [4] for a recent review p ¼ c20 ðq  q0 Þ; ð4Þ
of the method), which allows to obtain numerical solutions qffiffiffiffi qffiffiffiffi
of the continuum equations by defining the variables at a where c0 ¼ qe0 for the fluid and c0 ¼ qK0 for the solid,
set of suitable moving points, reconstructing the continu- being e the compressibility modulus of the fluid and K
ous field by means of interpolation functions centred on the bulk modulus of the solid. Eq. (1) is strictly valid only
each moving point. for compressible flows, while for incompressible flows it re-
The SPH technique was first developed in astrophysics duces to the divergence-free condition for the velocity field.
by Lucy [5] and by Gingold and Monaghan [6]. It was then However, in order to avoid the complexity of the implicit
successfully applied to the study of various fluid dynamics computation of pressure in a meshless method, the incom-
problems, such as free-surface incompressible flows [7], and pressible fluid can be studied as weakly compressible, thus
viscous flows [8,9]. Since the early 1990s, SPH was applied retaining the validity both of (1) and of the equation of
also to the simulation of elasticity and fragmentation in state (4). However, since the stability of an explicit numer-
solids: in particular, Libersky et al. [10] modelled the elastic ical integration of Eqs. (1) and (2) depends on the Courant
response of solid structures by an incremental formulation condition, the maximum time step is inversely proportional
of Hooke’s law. to the sound speed c0. It is therefore often necessary to as-
SPH has been also used to simulate the interaction sign to the compressibility modulus a value which is lower
between different fluids [11,12], different solids [13] and than the real one, in order to limit the computational time.
between fluids and structures [14] in presence of explosions. This leads to errors that can be reduced if a proper value is
In some commercial codes, an SPH description of the fluid assigned to e. In particular Monaghan [7] suggests that, in
motion is coupled to a finite element formulation for the order to limit density fluctuations to 1%, the Mach num-
solid dynamics, in order to simulate FSI problems. ber, i.e. the ratio between the local flow velocity and c0,
The present paper discusses a FSI model where both the must be everywhere lower than 0.1. Many applications of
fluid and the solid parts are modelled by SPH. Aim of the weakly compressible SPH models (see, for instance
model is the analysis of FSI problems where large elastic [9,11,15,16]) confirm that incompressible flows can be sim-
displacements of the solid occur, while rapidly moving ulated with good precision in this way.
free-surfaces characterize the fluid motion. If the dynamics of the fluid flow is dominated by inertial
The reliability of the numerical results yielded by the forces, viscosity effects can be safely neglected, and Sij = 0
proposed SPH FSI model is checked against laboratory can be assumed for fluids. For solids, the linear elastic rela-
data obtained during a simple 2D interaction experiment. tion between stress and deformation tensors can be derived
in time in order to obtain an evolution equation for Sij. The
2. Numerical model use of the corotational, or Jaumann, time derivative guar-
antees that the formulation is independent from superposed
2.1. Equations of motion rigid rotations, resulting in the incremental formulation of
Hooke’s law corrected by the Jaumann rate:
The motion of a continuum subjected to the action of  
DS ij 1
gravity, in isothermal conditions, is described by the conti- ¼ 2l Dij  dij Dij þ S ik Xjk þ Xik S kj ; ð5Þ
nuity equation Dt 3

Dq ovi where
þq ¼ 0; ð1Þ  
Dt oxi 1 ovi ovj
Dij ¼ þ ð6Þ
and by the momentum equation 2 oxj oxi
Dvi orij is the rate of deformation tensor,
q ¼ qgi þ ; ð2Þ  
Dt oxj 1 ovi ovj
Xij ¼  ð7Þ
2 oxj oxi
where t is time, q is density, vi is the velocity vector, xi is the
position vector, gi is the gravity vector, rij is the stress is the spin tensor and l is the shear modulus.
 C. Antoci et al. / Computers and Structures 85 (2007) 879–890 881

8  2  3
Eq. (5) belongs to the class of hypoelastic constitutive > 3 j~xj 3 j~
xj j~
xj
>
relations, which relate the objective measure of stress rate >1  2 h
>
>
þ4
h
06
h
6 1;
to the rate of deformation. For this kind of relations a r<  3
W ð~
x; hÞ ¼ 2 j~
xj j~
xj ð9Þ
stored energy function can not be defined: therefore, in h >>
>
1
2  16 6 2;
the case of large deformations, energy is not conserved, >4
> h h
:
as happens instead when hyperelastic laws are considered 0 elsewhere;
[17]. On the other hand, the relation (5) is rate-indepen-
dent, incrementally linear and reversible: therefore, when where r is a normalization constant equal to 10/7p, is used
small increments about a state of finite deformation are in the present work.
considered, the increments in stress are proportional to Neglecting viscous stresses, the SPH approximation of
those in strain and are recovered after unloading [18]. This the momentum equation (2) for the generic fluid particle
means that, if (5) is integrated in time with a conveniently a reduces to
small time step, it can be adopted as a constitutive model Dvia X p pb

oW ab
a
when small finite deformations are considered. It can be ¼ mb 2 þ 2 dij þ gi ; ð10Þ
Dt b
q a q b oa x j
observed that the integration in time of Eq. (5) from an ini-
tial equilibrium configuration preserves the proportionality where the symmetric formulation of the pressure gradient
between stress and deformation, which characterizes the term guarantees that the action–reaction principle between
elastic response. particles a and b is satisfied.
The use of the constitutive equation (5) for the solid, The SPH approximation of the momentum equation (2)
which has been proposed in an SPH context by [10,19], for solids becomes
has the advantage to allow for a common description of   
both the fluid and the solid dynamics in terms of pressure Dvia X p p S ij S ij
¼ mb  a2 þ b2 dij þ 2a þ 2b
and velocity. Dt b
q a q b q a qb

oW ab
þ Pab dij þ Rijab f q þ gi ; ð11Þ
oa x j
2.2. SPH formulation
where the last two terms between brackets have been intro-
In the SPH formulation, the continuum is divided duced in order to solve numerical problems connected with
into pseudo-particles of constant mass m. This mass is the meshless nature of the SPH method. The first one is an
distributed around the centre of mass of each particle, artificial viscosity term which has been proposed by [21] in
according to a density distribution defined by a suitable order to smooth out the velocity oscillations which can
kernel function [4], which has non-zero values only on arise owing to non-uniform particle distribution in space
a circle of radius 2h centred on the particle itself. The when particles get too close to each other:
length scale h, defined as ‘‘smoothing length’’, characterizes 8
the dimension of the domain of definition of the kernel < acab lab ð~ va ~
vb Þ  ð~
xa ~
xb Þ < 0;
function, and hence the discretization of the numerical Pab ¼ qab ð12Þ
:
scheme. 0 ð~
va ~
vb Þ  ð~
xa ~
xb Þ > 0;
According to the Lagrangian approach, each particle is
followed along its trajectory. The density and the velocity where cab ¼ 12 ðca þ cb Þ, qab ¼ 12 ðqa þ qb Þ and lab ¼
hð~
va ~
vb Þð~
xa ~
xb Þ
of each particle are updated by explicit integration of the xb j2 þð0:1hÞ2
xa ~
j~
.
continuity (1) and the momentum equations (2). Particle This term introduces numerical dissipation in the model,
trajectories are then computed by integrating in time the which must be kept to a minimum to reduce unwanted
material derivative of velocity obtained from Eq. (2). oscillations in the velocities without affecting the solution
According to the continuity equation (1), the material in a sensible way. The sensitivity to parameter a of the
derivative for the density of the generic particle a, fluid numerical solution of a test elastic problem is discussed
or solid, can be obtained by SPH interpolation as in Section 2.4.
The second term plays the role of an artificial stress and
Dqa X was proposed by [19] to eliminate the effects of the so-called
¼ mb ð~
va ~
vb Þ  ra W ab ; ð8Þ
Dt b ‘‘tensile instability’’ [14,22,23]. This instability, which is
strictly related to the interpolation technique of the stan-
where the summation is extended to all the particles b dard SPH method [24], is especially noticeable when simu-
around a and the use of the velocity difference ~ va ~vb in lating tension states in solids: particles tend then to clump
the SPH approximation of the divergence guarantees that, together, causing non-physical fractures in the material.
for a constant velocity field, the material derivative of den- The ‘‘artificial stress’’ term in Eq. (11) becomes effective
sity is zero. In (8) ra W ab is the gradient of the kernel func- when particle i is in tension and acts as a repulsive force
tion W ð~ x; hÞ evaluated at ~ x ¼~xa ~xb . The cubic spline avoiding particle clumping. In the form proposed by [19]
kernel proposed by [20] it is evaluated as
882 C. Antoci et al. / Computers and Structures 85 (2007) 879–890
 q
  W ðj~
xa ~xb jÞ velocity oscillations, but introduces some energy dissipa-
Rijab f q ¼ Rija þ Rijb ; ð13Þ
W ðdÞ tion in the numerical solution. However, its use avoids
the inclusion of the artificial viscosity term (12) in (10).
where q is a parameter, d is the mean initial distance be- A sensitivity analysis of the numerical solution of an
tween particles and Rij is obtained as follows. For each par- unsteady hydraulic problem to the smoothing coefficient
ticle the stress tensor rij is diagonalised. Then an artificial h is discussed in Section 2.4.
stress term is evaluated for any of the diagonal components An Euler explicit scheme of integration in time is
i which are positive:
r adopted for both the fluid and the solid structure. The time
r
ia integration scheme is staggered, i.e. the solution of momen-
Ria ¼ e ; ð14Þ
q2a tum equation is shifted by half time step with respect to
other variables. The choice of the time step depends on
where e is a parameter. the Courant stability condition. Since there are two media
The artificial stress in the original coordinates system Rij in the domain, the smallest of the time steps required by
is then calculated by rotating the coordinates back. each medium is chosen.
Gray et al. derive optimal values from the dispersion The adoption of an implicit integration, although possi-
equations, suggesting to use e = 0.3 and q = 4. An analysis ble in principle, would require the inversion of a large
of the sensitivity of the numerical solution to the parameter sparse matrix, thus being computationally intensive.
e is also discussed in Section 2.4.
The deviatoric stress Sij in (11) is calculated by integrat-
2.4. Sensitivity to the corrective terms
ing in time, by an implicit scheme, the rate of deviatoric
stress from the incremental hypoelastic relation (5). The
The sensitivity to the artificial viscosity and artificial
SPH estimate of velocity gradient, needed to compute the
stress terms in Eq. (11) has been checked by simulating
spin and the rate of deformation, can be obtained by first
the free oscillations of an elastic plate [19,26] having one
introducing a Taylor expansion of the velocity of particle
end clamped and the other one free. The plate is initially
a truncated to the first order, which leads to the following
horizontal (Fig. 1) and the initial velocity distribution is
system in the unknown components of the velocity gradient
assigned according to the analytical expression of the free
tensor:
oscillations of a thin plate:
X mb oW ab X mb oW ab
vib ¼ via f ðxÞ
qb oa xj qb oa xj vy ðxÞ ¼ vL0 c0 ;
b b
  X f ðLÞ
ovi mb oW ab
þ ðxkb  xka Þ: ð15Þ where
oxk a b qb oa xj
f ðxÞ ¼ ðcosðkLÞ þ coshðkLÞÞðcoshðkxÞ  cosðkxÞÞ
The solution of system (15), although needing to be solved
at each velocity gradient evaluation, is preferred to a direct þ ðsinðkLÞ  sinhðkLÞÞðsinhðkxÞ  sinðkxÞÞ;
SPH evaluation of the velocity gradient because it accounts
while vL0 ¼ 0:01 determines the initial velocity of the free
for a non-uniform particle distribution, especially close to
end. Since we consider the fundamental mode, kL ¼ 1:875.
the boundaries. The plate oscillates around the initial position, with ampli-
tude and period depending on elastic and geometric
2.3. Velocity correction and time integration scheme properties.
The results discussed in the following are obtained
Particle velocities ~
v obtained by time integration of the for:L = 0.2 m, H = 0.02 m, q = 1000 kg/m3, K ¼ 3:25 
momentum equation are corrected in order to smooth 106 N=m2 , l = 715000 N/m2.
out unwanted numerical peaks. The correction is obtained The computed non-dimensional amplitude and period
by of the first oscillation (A/L = 0.124 and Tc0 =L ¼ 81:5) are
P mb
bq ð~
vb ~ va ÞW ab
v~a ¼ ~
ab
~ va þ ð1  hÞ P mb ; ð16Þ
bq
ab
W ab

where qab ¼ 12 ðqa þ qb Þ and h is a smoothing coefficient.

For solid particles, these corrected values ~ v~ are used,
according to the XSPH scheme [25], in (8) and (15) and
in the computation of the updated position of particles,
while the uncorrected velocities ~ v are used for time integra-
tion of Eq. (11) at the following step. For fluid particles, ~ v~
are used also in the time integration of the momentum Eq.
(10), according to the scheme proposed by [15]: this scheme
proves to be more effective than XSPH in smoothing out Fig. 1. Scheme of the plate.
 C. Antoci et al. / Computers and Structures 85 (2007) 879–890 883

very similar to the ones computed by Gray et al. [19]

(AG =L ffi 0:125 and T G c0 =L ffi 82), although they both differ
to some extent from the analytical solution (Aa =L ¼ 0:115
and T a c0 =L ¼ 72:39).
The main effect that has been verified while checking the
sensitivity to parameters e and a consists in a damping of
the free oscillations of the plate. Fig. 2, where the ratio
between the amplitude of the oscillations after five periods
(A5) and the amplitude of the first oscillation (A1) is plot-
ted, shows that this damping occurs when e < 0.15.
Also the introduction of the artificial viscosity term con-
tributes to elastic energy conservation (Fig. 3): actually, for
a = 0 the computation is unstable, while damping of the
Fig. 4. Dam-break: profiles for different values of h at t = 0.36 s (in
free oscillations reduces only for a > 1. This is probably simulations A, B and C the following values have been assigned to h
due to the fact that the local oscillations of velocity cause respectively: h = 0.96, h = 0.92 and h = 0.84).
disorder in the spatial distribution of solid particles, thus
compromising interpolated estimates. When this disorder
occurs, a decrease in time of the amplitude of the oscilla-
tions is observable. Artificial viscosity avoids this problem,
thus leading to a satisfactory conservation of elastic energy.
The amplitudes of the first oscillation, as well as the oscil-
lation periods, do not vary significantly with a. On the
other hand the higher is e, the larger the amplitude A1 of
the first oscillation (Fig. 2), since the artificial stress reduces
the effective stiffness of the plate, reducing the elastic
Fig. 5. Comparison between a snapshot by [15] at t = 0.36 s and an image
from a simulation with h = 0.92.

response of the material. Thus care must be taken in the

choice of e, in order not to affect the solution excessively.
The sensitivity of fluid flow simulations to the velocity
correction scheme (16) has been checked on the classical
dam-break problem: a mass of fluid at rest (here having
an initial height of 0.1 m) moves after the instantaneous
removal of the confining wall at t = 0. The dissipative effect
of the velocity correction scheme can be observed in Fig. 4,
where free-surface profiles for different values of h at
t = 0.36 s are compared with the experimental data
Fig. 2. Influence of parameter e on the amplitude of the free oscillations of obtained by [15]. In particular, high values of h lead to
the plate. wave fronts speeds and to free-surface profiles closer to
the experimental ones. The best value of h must therefore
be chosen in order to assure stability to the computation
without introducing excessive numerical dissipation.
In Fig. 5, a snapshot of the dam-break experiment at
t = 0.36 s is compared with the corresponding particle dis-
tribution from the SPH simulation with h = 0.92.

3. Fluid–structure interaction model

Two different sets of particles are used for solid and

fluid. All the particles located farther than 2h from the
interface interact, of course, only with particles of the same
species (particle a2 in Fig. 6). On the other hand, when SPH
interpolation is performed on particles closer to the inter-
Fig. 3. Influence of parameter a on the amplitude of the free oscillations face (like particle a1 in Fig. 6), particles of both media
of the plate. are involved. The simplest approach in this case consists
884 C. Antoci et al. / Computers and Structures 85 (2007) 879–890

Fig. 8. Definition of tangent and normal vectors for solid ‘‘boundary’’

particles.

vif ni ¼ vis ni ; ð19Þ

Fig. 6. Fluid (black) and solid (white) particles near the interface.
rijs nis njs ¼ pf : ð20Þ
The application of (19) and (20) requires the precise com-
in extending the summations in (10) and (11) to all the par-
putation of the position of the interface surface and of its
ticles b regardless of their nature. Although interpenetra-
normal direction.
tion of fluid and solid particles can occur in principle, it
is prevented by the XSPH correction of velocity, as
3.1. Definition of the interface and its normal
observed by [13] in the case of contact between solids.
This method is equivalent to the introduction of cou-
A way to define the interface and its normal can be
pling conditions on the interface between a solid and a vis-
found in SPH estimates of constant functions and of their
cous fluid, i.e.:
gradient, as proposed by [14] to define complex solid
vif ¼ vis ð17Þ boundaries. However, when no fragmentation takes place,
as in the problems here considered, solid particles maintain
and
the same regular pattern in time, making the identification
rijs njs ¼ rijf njf ; ð18Þ of the interface easier. The interface is therefore defined
here as a line distant d/2 from the row of solid particles
where the subscripts s and f identify the solid and the fluid, closest to the fluid (Fig. 8).
respectively, while the sign convention is the same as in For every solid particle closer than 2h to the interface
Fig. 7. (‘‘boundary’’ particles) the tangent unit vector is defined by
This simple approach, which is the correct one for vis-  
cous fluids as it imposes automatically a no-slip condition xaþ1  xa1 y aþ1  y a1
^ta ¼ ðtax ; tay Þ ¼ ; ; ð21Þ
on the contact surface, is not valid when the effect of vis- xaþ1 ~
j~ xaþ1 ~
xa1 j j~ xa1 j
cosity is neglected and an inviscid flow formulation is where the subscripts a + 1 and a  1 identify the particles
adopted. The following kinematic and dynamic interface which immediately precede and follow particle a along
conditions, which imply the continuity of the normal com- the row parallel to the boundary.
ponent of velocity and of normal stress, must be therefore The normal unit vector is computed accordingly
applied in this case:
^na ¼ ðtay ; tax Þ: ð22Þ
The position of the interface point closest to each particle a
is also defined
 
1
xinta ¼ ~
~ xa þ r þ d^na ; ð23Þ
2
where r is the number of rows between particle a and the
boundary.

3.2. Dynamic interface condition

The dynamic interface condition, which is equivalent to

Fig. 7. Normal convention. the action–reaction principle, is guaranteed by imposing
 C. Antoci et al. / Computers and Structures 85 (2007) 879–890 885

that the force exerted by the fluid on the solid has the same the two media separate (and C ceases to be an interface
modulus as the force exerted by the solid on the fluid, but surface between solid and fluid), introducing a smooth
opposite direction. To obtain the SPH evaluation of the decrease in the interface pressure value.
force applied by one phase on the other, one can first con- The term ~F f!s;a =qa is therefore added to (11), which, for
sider that the ‘‘smoothed’’ estimate of the gradient of pres- boundary particles, is modified as follows:
 
sure in a domain X bounded by the surface CX is Dvia X rija rijb q oW ab
Z ¼ mb þ 2 þ Pab dij þ Rijab f
Dt q2a qb oa x j
hrpð~xÞi ¼ pð~ 0
x ÞW ð~
x ~
x0; hÞjCX þ pð~ x0 Þr~x W ð~ x0 ; hÞdX0 ;
x ~ b2Xs

X F if!s;a
ð24Þ þ gi þ : ð28Þ
qa
which becomes, when discretized by SPH: Since, in general, the position of fluid and solid particles
Z X mb across the interface is not symmetric, the action–reaction
hrpð~
xÞi ¼ x0 ÞW ðj~
pð~ x0 j; hÞdC0 þ
x ~ p ra W ab : principle must be imposed through a linear interpolation
CX b
qb b procedure: an interpolated value ~ F f!s;a is evaluated at the
ð25Þ position a, symmetric to a* across the interface (Fig. 9)
and the reaction term ~ F s!f;a is calculated as
The surface term on the right-hand side of (24) vanishes if
~
F s!f;a ¼ ~
F f!s;a : ð29Þ
X coincides with the domain of definition of the kernel
function, since W is zero for j~ x ~x0 j ¼ 2h. On the other Thus the continuity of normal stresses through the inter-
hand, when computing the pressure gradient term for the face (20) is assured.
solid ‘‘boundary’’ particle a (Fig. 9), the summation in The momentum equation (10) for fluid ‘‘boundary’’ par-
(25) is extended only to the solid particles b in Xs and ticles is modified by adding the force per unit mass
CX ¼ CXs . In this case, the surface term does not vanish ~
F s!f;a =qa :
any more and it accounts for the force per unit volume
~ Dvia X p  p  oW a b Fi 
F f!s;a exerted by the fluid on particle a. If the pressure ¼ mb a2 þ b2 dij þ gi þ s!f;a : ð30Þ
on the interface is approximated by its SPH interpolation Dt b2X
qa qb oa xj qa
f
limited to the fluid particles in Xf:
X mb
pinta ¼ 2 p W ðj~
xinta ~
xb j; hÞ; ð26Þ 3.3. Kinematic interface condition
b2X
qb b
f

The velocity of the points on the interface can be

the surface term in (25) can be expressed as obtained as
Z P mb
~
F f!s;a ¼ pinta W ðj~ x0 j; hÞdC0
xinta ~ ð27Þ b2Xs qb vib W ðj~
xinta ~xb j; hÞ
C viinta ¼ P mb : ð31Þ
b2Xs qb W ðj~xinta ~ xb j; hÞ
where C is the portion of interface surface included in the
domain of definition of the kernel (Fig. 9). The interface surface represents a moving wall boundary
It should be noted that (26) contains in fact an approx- for the fluid flow. This moving boundary condition is en-
imate normalization for the interface pressure pinta : actu- forced by a modified version of the ‘‘ghost particle tech-
ally, the summation nique’’ [16], where the solid particle positions are used as
P on mthe right-hand side should have
a geometric support to assign a fluid velocity distribution
been divided by b
b2Xf qb W ðj~
xinta ~
xb j; hÞ, which is an
SPH approximation for 0.5 when the liquid and solid able to reproduce, by SPH interpolation, the required nor-
phases are in contact. However, it has been found that mal component of ~ vinta on the interface.
the approximate normalization (26) works better when
4. Deformation of an elastic plate subjected to
time-dependent water pressure

The proposed SPH interaction model has been validated

by comparison of the numerical results with data measured
during suitable laboratory experiments. In these experi-
ments an elastic gate, clamped at one end and free at the
other one, interacts with a mass of water initially confined
in a free-surface tank behind the gate.

4.1. Experimental set-up

A section of a plexiglas flume has been confined by two

Fig. 9. Dynamic condition for solid boundary particle a. vertical walls, thus creating a tank, as represented in
886 C. Antoci et al. / Computers and Structures 85 (2007) 879–890

of the plate and the water levels in the tank have been mea-
sured by digital image processing.

4.2. Details of the numerical simulation

Although the width B of the channel is of the same order

of magnitude of the plate height L and of the initial water
level H, the analysis of the experiments shows that the
water transient flow and the plate deformation can be stud-
ied as a two-dimensional phenomenon. In particular, it
must be noted that no bending of the rubber plate was
observed along the direction normal to the lateral walls.
A 2D simulation in the plane (x, y) shown in Fig. 11 has
therefore been realized.
Water is considered as a perfect fluid with density
qf = 1000 kg/m3. As explained in Section 2.1, the incom-
Fig. 10. Scheme of the tank and of the gate: frontal view, lateral and plan. pressible fluid is analysed as weakly compressible to retain
the explicit solution of the continuity equation according to
(1) and to the equation of state (4). However, the time step
Fig. 10. One of the walls consists in an upper rigid part and yielded by the Courant condition [7] when adopting the
in a lower deformable plate made of rubber. The rubber real value of the compressibility module (approximately
plate is free at its lower end, thus representing an elastic 2 · 109 N/m2) is too small to obtain numerical solutions in a
gate closing the tank. The geometric dimensions of the sys- reasonable computing time. Therefore, a lower value of the
tem and the physical characteristics of the elastic gate are compressibility modulus ðe ¼ 2  106 N=m2 Þ is adopted,
reported in Table 1. The Young modulus E has been while the appearance of unphysical compressibility effects
obtained by performing tension tests in the range of defor-
mations which occur in the phenomenon. Since the Young
modulus for rubber depends on deformation, the reported
value is intended as an average value.
Being the rubber plate clamped only along its upper
side, it is free to deform when subjected to the pressure
of the fluid behind it. In order to assure sealing and mini-
mize friction along the lateral plexiglas walls of the tank,
a thin layer of transparent fat has been inserted between
the lateral sides of the plate and the tank, paying care to
distribute it uniformly along the sides of the plate, in order
to minimize its effect on the motion.
While the rubber plate is held fixed by an external rigid Fig. 11. Initial configuration.
support, the tank is filled with water up to the desired level.
When the water in the tank is in hydrostatic conditions, the
rigid support is suddenly removed, thus allowing the plate
to deform while water flows under it. The experiments were
recorded by a digital video camera, at a frequency of
25 frames per second. The displacements of the free end

Table 1
Dimensions of the system and physical characteristics of the rubber plate
Dimensions
A (m) 0.1
H (m) 0.14
B (m) 0.1
B* (m) 0.098
L (m) 0.079
s (m) 0.005
Rubber
q (kg/m3) 1100
E (MPa) ffi10
Fig. 12. Simulation: initial configuration.
 C. Antoci et al. / Computers and Structures 85 (2007) 879–890 887

Fig. 13. Frames and images from simulation every 0.04 s from t = 0 s (a) until t = 0.4 s (k).
888 C. Antoci et al. / Computers and Structures 85 (2007) 879–890

is avoided by requiring that the local Mach number is placements, whereas, close to the free end, the plate seems
everywhere lower than 0.1. to move almost as a rigid body.
The rubber plate is discretized by solid particles having The displacements of the plate in the simulation are
density qs ¼ 1100 kg=m3 . Since some uncertainty occurs in slightly larger than those in the experiment. This is possibly
the estimate of an average Young modulus for rubber, sim- due to the fact that the stiffness of the calculated plate is
ulations with different values of E were run. The results of lower than the real one. However, the differences here high-
the simulation with E ¼ 1:2  107 N=m2 , which better lighted are consistent with those found in the examples dis-
reproduces the experimental phenomenon, are discussed cussed in Section 2.4.
in the following. In addition, it must be noted that during the experiment
The Poisson coefficient is set equal to 0.4. This value, some leakage of water occurs besides the gate (as it can be
although lower than the theoretical one for incompressible viewed in Fig. 13 at times larger than 0.12 s). Owing to the
rubber, allows us to use larger time steps while respecting leakage, the pressure behind the plate in the experiment
the Courant condition; it has been checked that the influ- might be slightly lower than the one predicted in the simu-
ence of the modification of the Poisson coefficient on the lation, leading to a hydrodynamic force on the plate lower
numerical results is negligible. than the simulated one.
The adopted values of Young modulus and Poisson coef- In Fig. 14, the horizontal and vertical displacements
ficient correspond to a bulk modulus K ¼ 2  107 N=m2 computed for the plate are compared with those measured
and to a shear modulus l ¼ 4:27  106 N=m2 . in the digitalized images acquired during the experiments.
According to the results of the sensitivity analysis dis-
cussed in Section 2.4, the artificial stress parameter e and
the velocity smoothing parameter h have been set respec-
tively equal to 0.3 and 0.92. Since the phenomenon simu-
lated is rapid, a value a = 1 has been preferred for the
artificial viscosity parameter to avoid excessive numerical
dissipation. Fig. 12 shows the initial particle distribution
in the simulation.
Slip boundary conditions are imposed to the fluid flow
on the rigid walls: these conditions are imposed through
the method of ghost particles [16] on the right tank wall
and on the bottom (continuous lines in Fig. 12) and by
layers of fixed ghost particles in the upper rigid part of
the left tank wall.
Fig. 14. Horizontal and vertical displacements of the free end of the plate.
The clamp condition is imposed to the plate through a
layer of fixed solid particles. This layer is orthogonal to
the end of the plate. Only densities and stresses are updated
in time for these particles, whereas their velocities are set
equal to zero and their positions are fixed in time.
Two different refinement levels (d 1 ¼ 8:33  104 m
close to the plate and d2 = 0.002 m far from the plate) have
been adopted in the discretization of the fluid mass,
whereas for the solid is d = d1. The resulting total number
of particles is np ¼ 6012. The time step in the simulation is
dt ¼ 8:34  106 s. At the initial time, the fluid is assumed
to be in hydrostatic conditions, while in the plate stresses
and deformations are equal to zero.

4.3. Discussion of results

The comparison between frames from the experiment

and images of SPH particle positions at corresponding
times are shown in Fig. 13. The dynamics of the phenom-
enon, in terms of deformation of the plate and of evolution
of the fluid free-surface, is well reproduced. It can be
observed that the calculated shape of the plate is similar
to the observed one: the deformation is maximum near
the clamp, with null second derivative of the horizontal dis- Fig. 15. Water level (m) just behind the gate (a) and 5 cm far from it (b).
 C. Antoci et al. / Computers and Structures 85 (2007) 879–890 889

Fig. 16. rxx distribution at t = 0.15 s: colour scale ranging from rxx = Fig. 18. rxy distribution at t = 0.15 s: colour scale ranging from rxy =
150 000 N/m2 (white) to rxx = 150 000 N/m2 (black). 200 000 N/m2 (white) to rxy = 200 000 N/m2 (black).

(Fig. 16–18). The stress distribution seems to be qualita-

tively correct, showing the internal part of the plate in ten-
sion (dark particles), whereas the external part is
compressed (light particles). A maximum stress value of
0.75 MPa is found for the ryy component near the clamp,
while the other components show lower maximum values.
All the stress components tend to zero when approaching
the free end of the plate.

5. Conclusions

A fluid–structure interaction model based on the SPH

method has been described. In the model, both the fluid
and the solid phases are discretized by SPH particles. Once
the fluid–solid interface has been defined, coupling condi-
tions are imposed to particles close to it. In particular,
the action of the fluid on the solid is computed through
Fig. 17. ryy distribution at t = 0.15 s: colour scale ranging from ryy = the evaluation of an approximated surface integral of fluid
750 000 N/m2 (white) to ryy = 750 000 N/m2 (black). pressure. On the contrary, the action on the fluid is com-
puted by linear spatial interpolation from the previous
one, in order to satisfy the action–reaction principle.
It can be seen that the time evolution of the phenomenon is The SPH model, thanks to its Lagrangian nature, has
well described by the SPH simulation, even if the differ- the advantage of allowing an easy definition of the fluid–
ences discussed above lead to a 10% overprediction of the solid interface, even in presence of large displacements of
maximum horizontal displacement of the plate. the structure, and does not need any specific treatment
The evolution of the free-surface is also well reproduced for the free-surface of the fluid. The contact between differ-
by the simulation. Fig. 15 shows the water level history ent media is also automatically established when particles
immediately behind the gate and in the middle of the tank. of different media are closer than twice as the smoothing
Consistently with the larger vertical displacement of the length of the kernel function.
gate, and hence with the larger gate opening, the computed In addition no time shift is introduced in the calculation
flow rate is slightly higher than the real one, leading to a of the dynamics of the two media, since the variables of
faster decrease of the water level in the first part of the tran- both media are updated simultaneously. Anyway, in prob-
sient. Later than t = 0.075 s, the computed and measured lems where the characteristic times of the dynamics of the
values of the water level evolve in the same way. two media differ too much, the model can be easily modi-
The stresses rxx, ryy, rxy in the elastic plate are plotted fied in order to integrate the equations for the two media
at t = 0.15 s, when the displacements are maximum with different time steps.
890 C. Antoci et al. / Computers and Structures 85 (2007) 879–890

The model has been tested by comparing numerical [8] Takeda H, Miyama SM, Sekiya M. Numerical simulation of viscous
results and data collected during a laboratory experiment flow by smoothed particle hydrodynamics. Progr Theor Phys 1994;
92:939–60.
concerning the outflow under an elastic sluice gate. This [9] Morris JP, Fox PJ, Zhu Y. Modeling low Reynolds number
experiment has been chosen because it concerns the interac- incompressible flows using SPH. J Comput Phys 1997;136:214–26.
tion between a long and thin elastic structure and a fluid [10] Libersky LD, Petschek AG, Carney TC, Hipp JR, Allahdadi FA.
mass, in presence of large displacements of the structure High strain Lagrangian hydrodynamics. J Comput Phys 1993;109:
and of a free surface flow. 67–75.
[11] Monaghan JJ, Cas RAF, Kos AM, Hallworth M. Gravity currents
The results have shown that realistic predictions, both of descending a ramp in a stratified tank. J Fluid Mech 1999;379:39–69.
the displacement of the elastic structure subjected to fluid [12] Colagrossi A, Landrini M. Numerical simulation of interfacial flows
pressure and of the resulting fluid flow, can be obtained by smoothed particle hydrodynamics. J Comput Phys 2003;191:
by the SPH model, although some improvement in the 448–75.
treatment of the elastic solid dynamics still needs to be [13] Vignjevic R, De Vuyst T, Campbell J. The use of an homogeneous
repulsive force for contact treatment in SPH, WCCM V, Fifth world
achieved. congress of computational mechanics, Vienna, Austria, July, 7–12,
2002.
Acknowledgements [14] Randles PW, Libersky LD. Smoothed particle hydrodynamics: some
recent improvements and applications. Comput Methods Appl Mech
We are also grateful to Roberto Allieri and Ivano Brivio Eng 1996;139:375–408.
[15] Gallati M, Braschi G. Simulazione Lagrangiana di flussi con
for their valuable help in the execution of the experiments. superficie libera in problemi di idraulica. L’acqua 2000;5:7–18.
We acknowledge Dresser Italia S.r.l. for the financial [16] Gallati M, Braschi G. Numerical description of rapidly varied flows
support to the present research. via SPH method. In: IASTED Int Conf, ASM, Creta 2002.
[17] Simo JC. Numerical analysis and simulation of plasticity, Handbook
of numerical analysis, vol. VI, Part. 3. Elsevier Science B.V.; 1998.
References [18] Belytschko T, Liu WK, Moran B. Nonlinear finite elements for
continua and structures. John Wiley & Sons; 2000.
[1] Rugonyi S, Bathe KJ. On finite element analysis of fluid flows fully [19] Gray JP, Monaghan JJ, Swift RP. SPH elastic dynamics. Comput
coupled with structural interactions. Comput Model Eng Sci 2001; Methods Appl Mech Eng 2001;190:6641–62.
2:195–212. [20] Monaghan JJ. Smoothed particle hydrodynamics. Annu Rev Astron
[2] Bathe KJ, Zhang H. Finite element developments for general fluid Astrophys 1992;30:543–74.
flows with structural interactions. Int J Numer Methods Eng 2004; [21] Monaghan JJ, Gingold RA. Shock simulation by the particle method
60:213–32. SPH. J Comput Phys 1983;52:374–89.
[3] Le Tallec P, Mouro J. Fluid structure interaction with large structural [22] Swegle JW, Hicks DL, Attaway SW. Smoothed particle hydro-
displacements. Comput Methods Appl Mech Eng 2000:1–29. dynamics stability analysis. J Comput Phys 1995;116:123–34.
[4] Monaghan JJ. Smoothed particle hydrodynamics. Rep Progr Phys [23] Morris JP. Analysis of smoothed particle hydrodynamics with
2005;68:1703–59. applications, Ph.D. thesis, 1996.
[5] Lucy LB. A numerical approach to the testing of the fission [24] Rabczuk T, Belytschko T, Xiao SP. Stable particle methods based on
hypothesis. Astron J 1977;82:1013–24. Lagrangian kernels. Comput Methods Appl Mech Eng 2004;193:
[6] Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: 1035–63.
theory and application to non-spherical stars. Mon Not R Astron Soc [25] Monaghan JJ. On the problem of penetration in particle methods.
1977;181:375–89. J Comput Phys 1989;82:1–15.
[7] Monaghan JJ. Simulating free surface flows with SPH. J Comput [26] Landau LD and Lifšits EM. Teoria dell’elasticità. Editori riuniti
Phys 1994;110:399–406. Edizioni Mir, 1979.