7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A SPHALE method to model multiphase flows with surface tension
J. Leduc*, F. Leboeuf*, Michel Lance*,
E. Parkinsont and J.C. Marongiut
Laboratoire de Mecanique des Fluides et d'Acoustique, Ecole Centrale de Lyon,
University of Lyon, Ecully, FRANCE
t ANDRITZ Hydro, Vevey, SWITZERLAND
julien.leduc @eclyon.fr and jeanchristophe.marongiu @andritzhydro.com
Keywords: SPHALE simulations, meshless methods, multiphase flows, Riemann solvers, preconditioning
Abstract
This work focuses on a new multiphase method for Smooth Particles Hydrodynamics (SPH) simulations which uses
many advantages of the arbitrary LagrangeEuler (ALE) formalism. The use of an acoustic Riemann solver in the
SPHALE method allows sharp interfaces between fluids and is fully compliant with a purely lagrangian description
of the interface. The mass fluxes between volumes of control associated with different fluids are controlled by the
interfacial velocity of the Riemann problem which allows to cancel these fluxes. One of the major advantages of
the proposed formulation is that it remains stable for very high density ratios (air/water density ratio=1000) without
adopting unphysical speeds of sound in both media. In the limit of low Mach number flows, preconditioned Riemann
solvers are used to limit the numerical diffusion linked with upwind schemes.
Surface tension effects are also included in this model. The proximity of the SPHALE formalism with the Finite
Volumes method allows an adaptation of the Continuum Surface Force (CSF) formalism. An other model of surface
tension is presented. This Local Laplace Pressure Correction model (LLPC) integrates the Laplace law in the Riemann
solver between volumes of control associated with different fluids. The acoustic solver and the momentum equations
are then modified to take into account the surface tension. The LLPC model recovers sharp variation of pressure at
the interface and decreases the intensity of parasitic currents.
This approach is successfully validated on well documented academic cases: gravity waves and static/oscillating
droplet.
Introduction
Roman symbols
g gravitational constant (m.s 2)
p pressure (N.m 2)
Greek symbols
p density (kg.m 3)
(a surface tension coefficient (N. m 1)
Subscripts
1 left particule
r right particle
High speed water flows as jets from Pelton hydraulic
turbines nozzles, can lead to important interface defor
mations without phase changes. These kinds of de
formations are really hard to catch with eulerian nu
merical methods since the diffusion of the interface
smoothes the sharpness of the physical interface pertur
bations. On an other hand meshless/lagrangian meth
ods as SPH offers the possibility to naturally simu
late flows with large deformations. Furthermore pre
vious work on SPHALE formalism brought improve
ment in term of precision, stability of single phase
flows (Marongiu 2008). Different approaches were de
veloped to simulate multiphase flows with SPH method
(Grenier 2008), (Colagrossi 2003) or (Hu 2006). Gre
nier et al. (Grenier 2008) developed a method based on
Nomenclature
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the SPHALE formalism using a volume fraction ap
proach. The goal of the present study is to develop a
model which can simulate multifluid flows with surface
tension effects using an SPHALE formalism without
volume fraction. The advantage of such an approach is
that control volumes always belong to one single phase,
which allows a sharp transition from one fluid to another
and a lagrangian tracking of interfaces.
The SPHmeshless method
The Smoothed Particle Hydrodynamics (SPH) is a
meshless method based on an interpolation technique.
A continuous field f is approximated from discrete data
points through the formula:
(f(x)) = w fj W(xz x, h) (1)
J D,.
W is a regular function (usually called kernel function)
different from 0 on its support D, h is the smoothing
length (linked with the size of Dx) and wj is the weight
of the particle j.
W
Figure 1: SPH approximation
One of the major advantage of this formalism is that
it allows to compute gradient of continuous field using
the gradient of the kernel function with a simple formula
(here written far from the domain boundary):
(V f(x)) = w fjVW(xj x, h) (2)
The traditional SPH formalism is a purely lagrangian
formalism and suffers usually from instability and lack
of precision. In this paper the numerical method is based
on an Arbitrary Lagrange Euler formalism in the SPH
framework (1) and (2) (Vila 1999)
The SPHALE formalism
This formalism considers conservative form of Euler
equations in a moving frame of reference at the veloc
ity vo.
d j dQ + J j (v vo).ndS=
dt 12 S
0 a s (3)
j Qs.ndS + QvdQd
where 2 is the volume of control, S its boundary,
vector of conservative variables and Qs and Qv the sur
face and volume source terms. If surface terms are re
stricted to pressure term, the equation (3) can be rewrit
ten:
LO () + div (FE() vo ) Q (4)
where LOo is the transport operator associated to vo and
FE is the flux vector of the Euler equations.
From this point Vila observed that the discretization of
the conservative formulation leads to the appearance of
one dimensional Riemann problems between each pair
of control volumes (Vila 1999). If we consider two con
trol volumes i and j, their interaction results in contribu
tions along the direction joining i and j, with a discon
tinuous state evolution at midpoint. It can be expressed
as:
at ) x (FE 4 D t), .)
j( (X ), 0) (Ii if x( ) < 0
^j if z("ij) > 0
(5)
where is the unit vector between i and j (from i to j),
xij is the midpoint between i and j, x"3 is the curvi
linear abscissa along the straight line between i and j,
whose origin is taken at xij and <( and
of conservative variables at i and j respectively. Vila
showed that the solution can be expressed as:
z(xci) + Xo(t), 4)4))
S t (6)
Xo (t) .,,, ,.
I Jo
where
steady Riemann problem for Euler equations expressed
in a fixed frame of reference, and vo(xij, t) is the veloc
ity of the interface between points i and j.
Vila showed that the final set of discrete equations in
the SPHALE formalism is:
d
t (Xi) = vo(X, t)
d
(Wi) = oi Yj (vo (xj) vo(xi))ViWij
d (wiPi) + wi p wj2pEij(VEij vo(Xi, t)).ViWi, =
jdtD
S(ipivi) + i 3 uj2[pE,ijVE,ij (E,i Vo(xj, t))+
dtjEDi
PE,ij]ViWij = Uwipig
(7)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where (PE,ij, VE,ij)t = ij (A) is the upwind solution
of the moving Riemann problem.
This final formalism has different interesting proper
ties: first it allows an upwind solution for the veloc
ity which stabilizes the model, mass fluxes exist be
tween control volumes(which are no more "particles")
and third, one equation is added in comparison with tra
ditional SPH which takes into account changes of vol
ume control distributions. Finally this formalism is re
ally close to the Finite Volume formalism. More de
tails on this formalism are presents in the work of Vila
(Vila 1999) and Marongiu et al. (Marongiu 2008).
The multifluid model
General solution to the Riemann problem
The resolution of the system (7) needs the solution of
the nonlinear Riemann problem of the form (5). In this
work the behaviour of the fluid is governed by the Tait
equation of state:
and the local speed of sound c is given by:
aPO
and the local speed of sound c is given by:
S (p + B)
Pressure and density derivatives are then linked by
speed of sound as:
Op
O p
at
adp
c20
C:c
SPI*
)v(1)
(2)
(2) *\
P,
vr(*
(2)
V,
Figure 2: Riemann problem scheme
if pi < pr, the Riemann problem will be structured by
a shock wave between the state I and I*, a rarefaction
wave between the state r and r*, and a contact wave
between 1* and r*. Considering low Mach number
flows and considering that the transport field is chosen
so that 0 < vo II
problem has to be found in the star region. It means
that the vector (,., 'E,j)t of the equation (7) is
equal to (p*,v() '' )t (p* is computed with the
Tait equation from p*, and v(2)* = v () if x/t < v1)*,
otherwise v(2)* = v2 ).
Linearized approximate Riemann solver for multi
phase simulations
Consider here that the contact discontinuity is the in
terface between two fluids, to determine the solution of
the Riemann problem, the method described by Murrone
is applied (Murrone 2008). It consist in a linearization
of the problem in each side of the contact discontinuity:
1 3.(q* q,) =0
1(q ) 0
It,.(q* o,
It results the following Riemann problem written fo
variable U (p, v(1), v(2))
BU BU
au + A(u) =S
at A(U)a O
The left eigenvectors of the matrix
(1/pc, 1, )t,12 (0,0, 1)t and 1 
(11) From (12), we gets two equations (13):
A are 11
S(1/pc1, .)t
where v(1) is the velocity normal to the interface
between i and j, v(2) is the tangential velocity and
V () pc2 0
A = 1/p v(1) 0 is the jacobian Matrix of
0 0 v(1)
the system. The eigenvalues of the jacobian matrix A
are: A1 v(1) c, A2 v() and A3 v(1) + c.
Since we are considering low Mach numbers in our
applications, A, represents a non linear wave moving
to the left and A3 a non linear wave moving to the
right. Fig. 2 shows the structure of the solution of the
Riemann problem using the Tait equation. It means that
1(v )) V (p+ )
prCr(V)* V(1)) ( Pr)
Since v(1) and the pressure are continuous at the in
terface between two fluids (if surface tension effects are
neglected), one gets:
{ (i)* V 1)* v*
p**r
P1 Pr P =
Then the state star (*) is computed as:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1mM
P1C1 + PrCr
CPyCl P(v(1)
pl cl + Pr cr
Plcl + prce
V(1))
v^ )
(15)
The tangential velocities v(2) are conserved through
the left and right wave and discontinuous through the
contact wave.
Solution of the multiphase Riemann problem
The ALE property is used to impose the velocity of
the interface vo equal to the resulting velocity from the
Riemann solver v* if the left and the right states are as
sociated with two different fluids. In this case the term
VE,ij v(xij, t) in the equation (7) is identically null
and the mass exchanges between the volumes of control
are blocked. From this point since the mass fluxes are
blocked at the interface by imposing the interface veloc
ity between two volumes of control, we then assume im
plicitly that the interface is described with a lagrangian
form of the ALE formalism. This solver has the advan
tage to be based on velocity and pressure which are con
tinuous at the interface without surface tension. The re
sulting density is computed from the pressure using the
Tait equation. The first terms of the equations (15) show
that the resulting velocity is balanced upwind on the side
of the heaviest fluid while the pressure is weighted on the
side of the lightest fluid, as expected across an interface
between a liquid and a gas for instance. In more prac
tical terms the heavy fluid (water) imposes the velocity
and the light fluid (air) imposes the pressure.
A special treatment is also applied to the equation for
the weight wc (7) when the two states correspond to dif
ferent fluids. This equation is in this case:
d (W'i)
dt
Ui wj2(vE,ij o(xi))ViW, (16)
j D,
The weight variations are no more supported sym
metrically by the two fluids. For an heavy fluid
.,, .) vo (xi) is nearly equal to 0 which implies that
the light fluid will support all the weight variations.
Gravity wave test case
The numerical developments are evaluated on the
gravity wave test case. The geometry is described on
figure 3 and the reference data are given in table 1. Wall
on the top and bottom of the domain and periodicity on
the side of the domain are used as boundary conditions.
Three different discretization are used to perform the
simulation: 60*60, 80*80, and 120*120.
Figure 3: Geometry for gravity waves test case
0.45
0.4
0.35
E 0.3
~0.25
0.2
a 0.15
0.1
0.05
0 0.5 1 1.5 2 2.5
time [s]
Figure 4: Evolution of the kinetic energy depending on
discretization
The error on frequencies are 2.66%, 0.791% and
0.50% for the different discretizations (Fig 5, Ax[m]
1/60, 1/80 and 1/120 respectively) These differ
ences are of the same order of magnitude in com
parison with the results obtained with eulerian codes
(Bonometti 2005). An important numerical diffusion
is responsible for the decay of the kinetic energy (it is
recalled that no physical viscosity is modelled in the
present study). This phenomenon is sensitive to the nu
merical discretization (see figure 4).
Preconditioning of Riemann Problems
The previous results show a satisfactory prediction of
the wave frequency but also an important numerical dif
fusion. Preconditioning techniques are widely used in
Finite Volumes methods to decrease this excessive nu
(1) (1)
pc* i + PrcrV
PiCl + PrCr
* PiClPr + Pr CrPtI
SPiCl + PTrCr
g,
.....
S 0.01
Figure 5: Frequency error depending on discretization
Table 1: Configuration gravity waves
domain length L [m] 1
initial deformation a/L [] 0.01
density fluid 1 [kg.m 3] 1000
density fluid 2 [kg.m 3] 1
sound speed 1 [m.s 1] 5
sound speed 2 [m.s ] 15
surface tension [N.m 1] 0
gravity [m.s 2] 9.81
medical diffusion obtained with Godunov type schemes
for low Mach number flows, and were introduced by
Murrone et al. (Murrone 2008) for multiphase flows.
The linear Riemann solver in (11) can thus be modified
by multiplying the Jacobian matrix A associated to the
Euler equations by the preconditionner of Turkel:
32 0 0
P 0 1 0 (17)
0 0 1
where 3 is of the order of magnitude of the Mach
number. For the following results, preconditioning
is applied for both single and multiphase interac
tions: interactions between two particles of the same
fluid are treated by a VFRoeTurkel Riemann solver
(Marongiu 2009) and those involving two particles of
different fluids by a preconditioned acoustic Riemann
solver.
The results on figure 6 show a positive effect of pre
conditioning on the numerical diffusion. This effect is
smaller for fine discretization. The preconditioning has
also a positive effect on frequency determination since
UV.4
d 0.3
0.2
0.1
0
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Non preconditioned 80*80 
Preconditioned 80*80 
Non preconditioned 120*120 ..........
Preconditioned 120*120
A
; l\ \*
0 0.5 1 1.5 2 2.5
time [s]
Figure 6: Comparison between nonpreconditioned and
preconditioned schemes
the errors are 0.61% for the 80*80 discretization (instead
of 0.79% without preconditioning) and 0.47% for the
120*120 discretization (instead of 0.50% without pre
conditioning).
Surface tension models
Surface tension effects can lead to different multiphase
phenomena such as coalescence or capillary waves.
This topic was considered in various studies with SPH
methods (Hu 2006). Two models are studied here,
one based on the Continuum Surface Force (CSF,
(Brackbill 1992)) method and one based on a novel
approach related to SPHALE formalism and the use of
Riemann solvers, the Local Laplace Pressure Correction
(LLPC).
Continuum surface force
This classical model arises from the work of
(Brackbill 1992) who rewrites this surface effect as a
volume force. A mathematical development allows to
rewrite this volume force as
(c Vc
F8 V.1
VcK1)
where c is the smoothed color function (c=0 or 1 far from
the interface and varies continuously near the interface).
The previous equation (18) considers that the force is
symmetrically applied to the two fluids. As it was done
by Brackbill (Brackbill 1992) the previous force will be
balanced between the two fluids when there is an impor
tant density ratio between the fluids. For airwater inter
face, it will model the fact that surface tension is mainly
arising from interactions between water particles.
To keep the conservation property of the previous
model, no fluxes of surface tension are exchanged be
tween volumes of control from different fluids and a
constant coefficient (2 P1 ) is added based on the
PD1 +PD2
reference density of each fluid. The final expression for
the surface tension effect on particle i in this formalism
and for a high density ratio is given by:
Fti = 2 PD1 : Wj (Fj Fi) VWi (19)
PD, +PD2 cD1
where D1 is the group of particles from the same fluid
as the particle i, pD, is the reference density of the par
ticle i and j (which are from the same fluid) and
Vci < Vci
F, v VcI (20)
Since the coefficient PD1 is constant, this
PD1 PD2
formalism is conservative. For low density ratios the
computation of surface tension effects is based on the
first formalism of the volumic force (18). The limit of
the density ratio to use the equation (18) or the equation
(19) is around ten but nor clearly defined since the use
of the formalism (18) allows to use the fluxes on each
side of the interface to compute the force on a particular
control volume. On the contrary the equation (19) uses
the fluxes between control volume from the same fluids.
Local Laplace pressure correction LLPC
The previous surface tension model applied to the
acoustic multiphase model has the major drawback that
a volume force representing a surface phenomenon is
applied to a sharp density interface. Indeed Brackbill
(Brackbill 1992) developed the original model on a
smoothed density transition through the interface. In
order to take benefit from the sharp interface obtained
with the lagrangian feature of SPH, a novel surface
tension model is developed which can reproduce the
sharp pressure jump predicted by Laplace law. In the
presence of surface tension, the difference of attraction
between particles is creating a resulting pressure jump
at the interface. Between two control volumes from
different fluids, the pressure continuity written in (14)
is then no more valid to describe this effect. As a con
sequence the pressure jump is directly included in the
acoustic Riemann solver between two control volumes
corresponding to 2 different fluids (Perigaud 2005):
S l)* (1)* *(21)
r (21)
where K, is the local surface curvature from the point of
view of I (Kr = Kl) and a the surface tension coeffi
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
cient. The second equation of (21) can be written as:
S Pi Pr
Pi  cTl = p + rKl
Pi + Pr Pi + Pr
We can introduce a pressure p*:
Pl *
P = pi p acl = Pr
P1 + Pr
Pr
P1 + Pr
c(Kr (23)
Introducing this last equation for the resolution of the
linearized Riemann problem, one gets:
(1) (1) 
P* 1 ;c1 + Crcr ) (Pr Prst) (Pi
PlCl + PrCr PlCl + Prcr
* PlCl(pr Prst) + PrCr(Pl Pst)
PiCl + PrCr
Plce PrCr
PlClprCr(V +) V 1))
plCl + PrCr
with:
Pist
Prst
P1
P1 + Pr
Pr
P + P Kr
P1 + Pr
Pist)
(24)
(25)
The equations (25) are balancing the effect of surface
tension from one side to another depending on the den
sity ratio. Indeed if we are considering a static liquid
droplet at rest in a stagnant gaz, the Laplace law predicts
a pressure jump /' ) through the interface. Figure
7 shows that the pressure terms in the Riemann solver
(24) are canceling: There are no flux momentum created
through the Riemann solver by the pressure jump for a
static droplet.
pressure
Kt < 0
P, = 1
Plst = 1g Kl
Pt Pist 
Pl
Pg
r > 0
i 9
aK
Figure 7: Example of pressure treatment for LLPC
model (see eq. 24)
In this model the fluxes linked with surface tension
are blocked at the interface (as for the CSF model for
large density ratios). As figure 7 shows, effects of sur
face tension are removed from pi and Pr before comput
ing the pressure through the Riemann solver (24). The
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 2: Configurations static dropletbubble simula
tions
droplet radius [mm] 10
surface tension coeff. [Nm 1] 0.072
density Fluid 1 [kg.m 3] 1000
density Fluid 2 [kg.m 3] 1
pressure components arising from surface tension terms
have to be added in the momentum equation when the
two volumes of control are from different fluids. The
momentum equation is then:
d
dt (wpivi)
wi wjuj2[PE,ijvE,ij ) (vE,ij o0) +PE,ij + Pijst] ViWj
3
= wipi9
(26)
where/ A(i) ' TKi and A(ij) 1 if i and j
are from different fluids, otherwise and A(ij) 0. Since
we are computing low Mach number test cases, we get
that the solution of the Riemann problem (i, i, VE,ij)
is equal to i / v*) computed for the interaction between
i and j. It has to be noted that this model is close from
the one developed by Francois et al. (Francois 2006)
Validation
Two test cases were chosen in order to validate the
different models: a static droplet surrounded by air and
a droplet with an initial deformation. The first test case
is useful to test the ability of the model to reproduce the
Laplace Law without developing high spurious currents.
The second test case allows to compare oscillation
frequency with the value predicted by the theory in the
case of small deformations.
Static droplet
As it could be expected, the pressure jump is more accu
rately reproduced with the LLPC method than with the
CSF method. On figure 9 the pressure jump occurs di
rectly from a particle of air to a particle of water. The
pressure jump is much smoother with the CSF method
since three layers of volumes of control are needed to
perform the interface transition (Figure 8).
Parasitic currents are also reduced by the use of
LLPC model in comparison with the CSF model. Figure
10 show the evolution of parasitic currents over a period
(0.5s) that can be considered as long with respect to a
typical time scale for flows in Pelton turbines (3.10 3s).
Pressure(x,y)[Pa]
Figure 8: Pressure field for CSF model
Pressure(x,y)[Pa]
4*"^fti'tm" *
Figure 9: Pressure field for LLPC model
Oscillating droplet
The test case of a free oscillating droplet is performed
with the LLPC model. At t=0s, an initially deformed
droplet begins to oscillate under the effect of surface ten
sion. The initial deformation is around 1% (27) (differ
ence rate between longer and smaller diameters).
r =o(l+ 1+ (3cos(20) +1)
4
The evolution of kinetic energy is shown on figure 11,
for preconditioned Riemann solvers. The theoretical pe
riod is 2.868101s. This model gives smaller periods
(higher frequencies) than expected by the linear theory.
The errors are 2.82% for the discretization (correspond
ing to 10 volumes of control in the diameter) and 0.674%
with 20 control volume along the diameter.
Conclusions
The previous study demonstrates the interest of using the
SPHALE formalism with an acoustic Riemann solver
to simulate multiphase flows. It gives a model with
a sharp density interface, stable for high density ratios
and which respects a physical ratio of sound speed be
tween the fluids. The Local Laplace Pressure Correction
le08
le09
Sle10
S lell
le12
le13
0 0.1 0.2 0.3 0.4 0.5
time [s]
Figure 10: Evolution of parasitic currents
3.5e08
3e08
2.5e08
2e08
1.5e08
le08
5e09 KV I '
0 0.1 0.2 0.3 0.4 0.5
time [s]
0.6 0.7
Figure 11: Oscillating Droplet: kinetic Energy evolu
tion
(LLPC) model for surface tension gives a sharp pres
sure variation at the interface when CSF formalism re
quires 3 volumes of control to establish the good level of
pressure inside a droplet. The parasitic currents are also
strongly reduced with the LLPC model. The use of pre
conditioned Riemann solvers enable a reduction of the
numerical diffusion. This methods will now be used on
more complicated academic test cases and also on three
dimensional industrial cases.
Acknowledgements
The authors wish to thank the ANRT (association na
tionale recherche technologies these 407/2007 for help
ing to fund this work.
L_ ." CSF 10
CSF 20
LLPC 10
LLPC 20
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
References
[Bonometti 2005] T. Bonometi, D6veloppement d'une
m6thode de simulation d'6coulements a bulles et a
gouttes, Phd Thesis, INP Toulouse, 2005.
[Brackbill 1992] J.U. Brackbill, D.B. Kothe, and
C.H. Zemach, A continuum method for modeling
surfacetension, J. Comp. Phys., 100(2), 335354,
1992.
[Colagrossi 2003] A. Colagrossi and M. Landrini,
Numerical simulation of interfacial flows by
smoothed particle hydrodynamics, J. Comp. Phys.,
191(2), 448475, 2003.
[Francois 2006] M. Francois, S. Cummins, E. Dendy,
D. Kothe, J. Sicilian and M. Williams, A balanced
force algorithm for continuous and sharp interfa
cial surface tension models within a volume track
ing framework, Journal of computational Physics,
213, 141173, 2006.
[Grenier 2008] N. Grenier, D. Le Touze, P. Ferrant and
J.P. Vila, Twophase simulations using a volume
fraction SPH scheme with a Riemann solver, Proc.
of 3rd Int. SPHERIC Workshop, Switzerland, 173
179, 2008.
[Hu 2006] X.Y. Hu, and N.A. Adams, A multiphase
SPH method for macroscopic and mesoscopic
flows, J. C. Physics, 213, 844, 2006.
[Marongiu 2008] J.C. Marongiu, F. Leboeuf,
E. Parkinson, Riemann solvers and efficient
boundary treatments: an hybrid SPHfinite volume
numerical method, Proc. of 3rd Int. SPHERIC
Workshop, Switzerland, 101108, 2008.
[Marongiu 2009] J.C. Marongiu, F. Leboeuf, J. Caro.
E. Parkinson, Low Mach number numerical
schemes for the SPHALE method. Application in
free surface flows in Pelton turbines, Proc. of 4th
Int. SPHERIC Workshop, Nantes, 324328, 2009.
[Murrone 2008] A. Murrone and H. Guillard, Behavior
of upwind scheme in the low Mach number limit:
III. Preconditioned dissipation for a five equation
two phase model, Computers and Fluid, 37, 1209
1224, 2008.
[Perigaud 2005] G. Perigaud and R. Saurel, A com
pressible flow model with capillary effects, Journal
of Computational Physics 209, 139178, 2005.
[Saurel 2003] R. Saurel, S. Gavrilyuk and F. Renaud,
A multiphase model with internal degrees of free
dom: application to shockbubble interaction,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Journal of Fluid Mechanics, vol. 495, 283321,
2003.
[Vila 1999] J.P. Vila, On particle weighted methods
and Smooth Particle Hydrodynamics, Mathemati
cal Models and Methods in Applied Sciences, 9(2),
161209, 1999.
