Proceedings of ITC/ISHW2007

Derivation of jump conditions in multiphase incompressible flows
with singular forces
Caesar O. Harahap1) and Hideaki Miura2)
The Graduate University for Advanced Studies (SOKENDAI), Toki 509-5292, Japan2)
National Institute for Fusion Science, Toki 509-5292, Japan(Received 15 October 2007)
With the aim to numerically solve problems of plasma dynamics consisting of multiple phases and/or diﬀerent governing equations, an immersed interface technique is extended to solve a system with mass density andviscosity jumps. The jump conditions for velocity, pressure and their derivatives (sets of coupled equations) are
derived.
Keywords: dicontinuities, singular forces, jump conditions, immersed interface method
DOI: 10.1585/pfr.1.001

1Introduction

Clarifying key physics of complex behaviors of a hot
plasma in a magnetic confinement device such as the Large
Helical Device (LHD) is essential for understanding experimental results. For this purpose, various kinds of simulation studies such as magnetohydrodynamics (MHD), twofluid equations, Vlasov and gyrokinetic simulations havebeen conducted. (See, for example, Refs.[1, 2, 3]) TheMHD simulation class is the simplest among them in thesense that the system of equations consists of a relatively
few number of equations, and the number of dimensions of
the independent variable space is only three. Nevertheless,there still remain many diﬃculties mimicking the device
geometry in detail and taking various experimental and engineering parts into account. One of those diﬃculties is inconnecting the hot plasma region and the vacuum region,
which is complicated because the governing equations ofthese regions can be diﬀerent from each other.
In MHD simulations, the vacuum region is often described by the MHD equations with very low pressure
and/or mass density, or simply omitted from the simulation by imposing the boundary condition on the outermost
magnetic surface. For the purpose of studying the eﬀects
of plasma deformations around the hot plasma boundaryor to study peripheral regions and hot plasma core simultaneously, the former approach is preferable. In the former
simulations, very large resistivity and/or viscosity are imposed to the low pressure region. However, jumps of physical variables such as mass density, pressure, and temperature often causes numerical oscillations. Although suchoscillations may be avoided by adopting some numerical
techniques such as the Godunov, TVD or CIP scheme, thecomputation becomes complex and the numerical accuracy
can become ambiguous.
Recently, a class of numerical technique called the Imauthor’s e-mail: harahap.caesar@LHD.nifs.ac.jp

mersed Interface Method (IIM) is developed to simulate
the motion of a neutral fluid separated into two regions bya moving surface exerting forces. (See Xu and Wang [5]and references therein.) The IIM is a variant of the immersed boundary method, which was originally developedby Peskin[6] to simulate blood flows and has become oneof the most powerful tool to solve fluid flows with materialboundaries. The basic ideas of the IIM are (1)to utilize thegeneralized Taylor expansion[4] and (2)to derive the jumpconditions required in (1). A schematic figure is shown inFig.1. In a computational box, we have a material surface
S. The dashed lines are some typical grid lines. Utilizing
the generalized Taylor expansion, derivatives of a physical quantity u(x) can be approximated in a manner similar
to the standard finite diﬀerence method, provided that thejump conditions for u(x) and its derivatives on the material
surface1 are already derived. Then the motion of the fluids is marched by standard techniques such as the RungeKutta-Scheme and the motion of the discontinuous surface
can be tracked by a Lagrangean technique.

2Jump conditions

While IIM is an excelent technique, at present its application is limited to the incompressible Navier-Stokes equation (for neutral fluids). One important step to make thetechnique applicable to fusion plasma simulations is to derive jump conditions for various physical quantities (including magnetic field and current density) in the casewhere the viscosity and the mass density are dicontinuous accross a surface. In this article, we present thejump conditions for velocity, pressure, and their derivatives in such a case. For simplicity, we restrict ourselves tothe three-dimensional incompressible Navier-Stokes (NS)
1 Note that the words material surface mean a flow does not go acrossit. The surface is not necessariry a wall but can be, for example, a fluidsurface with surface tension.

Proceedings of ITC/ISHW2007
Computational Box

Fluid 2

α2 α1Fluid 1bτnu(x)material surface
Fig. 2 Coordinate lines α1 , α2 which cover the material surface
and the tangential, normal and binormal unit vectors.

xFig. 1 A schematic view of two fluids separated by a materialsurface.

equations:
 i∂u∂ui+ ujρ∂t∂x j==Fi=×∂X α , α , ti1∂t2−  i
∂u
∂
∂ui
∂p++µ∂xi ∂xi∂x j ∂x j+F i ,(1)0,
(2)



f i α1 , α2 , tS
δ x − X α1 , α2 , t dα1 dα2,(3)=  ui X α1 , α2 , t , t ,(4)where (α , α ) are Lagrangian parameters of points on theforce
surface at a reference time, f i is a Lagrangian density(due to surface tension, for example), δ x − X α1 , α2 , t
is a three dimensional delta function, and X α1 , α2 , t are
Cartesian coordinates of the surface points. (See Fig.2).
Taking the divergence of the momentum equation, a Poisson equation for pressure is obtained.

∂2 p
∂F i
∂2∂ui=+2µ∂xi ∂xi∂xi∂xi ∂x j ∂x j
  i
∂u
∂
∂ui− i ρ+ uj.
(5)
∂x
∂t
∂x j
1

n =J

∂u
∂xi

i

vector n normal to the surface (in Fig.1, pointing towardsfluid 2) is given by

2

Diﬀerentiating X we obtain two tangent vectors along
the coordinate lines of the surface.
∂Xτi =(i = 1, 2).
(6)∂αi
The three components of the tangential vectors τ i in theCartesian coordinates are denoted by (τ 1i , τ2i , τ3i ). The unit=1(τ1 × τ2 ) ,J|τ1 × τ2 | .

(7)
(8)

We also introduce three mutually orthogonal unit vectors
as τ = τ1 / |τ1 | , n, b = n × τ.
We assume that fluid 1 and fluid 2 have constant but
diﬀerent mass densities and viscosities. To show the jumpconditions, we define some symbols here. A jump of aphysical quantity accross the two phases is denoted by [.]where[.] X(α1 , α2 , t), t = (.)+ X(α1 , α2 , t), t− (.)− X(α1 , α2 , t), t .(9)The superscript + (− ) denotes the side towards which npoints, while − denotes the other side.
The jump conditions are derived as follows. First, weformulate some basic jump coditions following the workby Xu and Wang[4]. No slip condition on the dicontinuous
surface S leads to(10)
ui = 0.
and the acceleration jump condition∂ui∂ui= 0,+ uj∂t∂x j(11)Extensions from the works by Xu and Wang[4, 5] appear in
the jump conditions associated with the viscosity and themass density. The jump condition for the viscous stresstensor in eq.(1) is given by k
 if nk∂uk∂u jfi∂uj+−2µτk+=−nµ∂x j ∂xiJJ∂τ∂uk k ib n−2 µ(12)∂bProceedings of ITC/ISHW2007
The pressure jump condition is given by

One way to utilize eqs.(15) and (16) is by giving the source
terms in the right-hand-side by some appropriate interpolations or combinations of extrapolations from the twosides of each phases, forming a closed system of equations. Then the jump conditions can be obtained by inverting coeﬃcient matrices in eqs.(15) and (16), enabling
us to discretize the NS equations by the IIM method, andsolve problems of multiphase flows with singular forces.
Jump conditions for µ∂2 ui /∂x j ∂xk are given as follows.

f i ni∂uk∂uk k−2 µτk − 2 µb,(13)J∂τ∂bin which the jump of the viscosity is included. The jumpcondition for the pressure normal derivative 1∂p i1  ∂ ∂f2f
 1 + 2 n =∂xi
J ∂α∂α 2 i∂ u i ∂2 ui i ∂ui ∂ni ∂ui ∂ni+n + 2n ++2 µ∂τ ∂τ∂b ∂b∂τ2∂b

ii
∂u∂u−
ρ
(14)
ni .
+ uj ρ∂t∂x jp=

contains both the viscosity jump and the mass density
jump, making the formulation much more complicatedthan those in Refs.[4, 5].
After some manipulations, we have a set of coupledequations for the jump conditions as  ∂u1
 1  µ ∂x1
2300 0 00   τ1 τ1 τ1 0
  µ ∂u12 τ1 τ2 τ3 0
   ∂x100000
 222  µ ∂u3
 2n1 n2 n3 n2 0
30n00   ∂x
∂u2111
 00   µ ∂x1
0 0 τ1 τ2 τ3 0 0

2


 00  ·  µ ∂u
0 0 τ21 τ22 τ23 0 0
∂x2



 0 n1 0 n1 2n2 n3 0 n3 0   µ ∂u23   ∂x
0 0 000 τ 11 τ12 τ13   µ ∂u3
 0

 ∂x10 0 000 τ 21 τ22 τ23    ∂u3
 0  µ ∂x2  ∂u3
00 n1 00 n2 n1 n2 2n3µ ∂x3

1 ∂u  µ ∂α1 ∂u1 
 µ ∂α2  r1 u µ ∂u2 1∂α ∂u2  ,=  µ ∂α2  ru2 

3 
∂u 
 µ ∂α1 
∂u3 
 µ ∂α2 ru3
and 1 τ1 2
 τ1
 1n

τ12
τ22n2

   ∂p
τ13    ∂x1
  ∂pτ23   ∂x2
 n3  ∂p3∂x
 = ∂r p  1  ∂α ∂r p2  , ∂α r pm

 C11 C21

C31

C12C22C32

 B11−  B21

B31


















Since the expressions of the block matrices are lengthy,
here we show only the first line of the expression:

X1=C11=(15)∂uk∂uk kf i ni
−2 µ
τk − 2 µb,
(17)
J
∂τ
∂bfi
rui = − + R p ni ,(18)J 1
1  ∂ ∂f2f 1 + 2 r pn =J ∂α∂α
 2 i∂ u i ∂2 ui i ∂ui ∂ni ∂ui ∂ni+n + 2n ++2 µ∂τ ∂τ∂b ∂b∂τ2∂b
ii
∂u
∂u−
ρ
(19)
ni .
+ uj ρ∂t∂x j=  ∂2 u1  µ ∂x1 ∂x1   ∂2 u1  µ ∂x1 ∂x2   ∂2 u1  µ 1 3  ,(21)  ∂x∂2 ∂xu1
 µ ∂x2 ∂x2


2 1u
 µ ∂x∂2 ∂x  ∂2 u1 3 µ ∂x3 ∂x3  2
 2 τ12 τ11 τ21
2 τ11 τ3111 11 22 1
τ11 τ32 + τ31 τ12
 2 τ1 τ2 τ1 τ2 + τ1 τ2
  2
 2 τ12 τ12 τ222 τ12 τ32
2111221
 3 n τ2 n τ1 + n τ1 2 n1 τ31 + n3 τ111
11
 3 n τ2 2 n1 τ22 + n2 τ12 2 n1 τ32 + n3 τ12

100
 2
 2 2 τ21
2 τ21 τ31
2 τ31 
2 τ21 τ22 τ21 τ32 + τ31 τ22 2 τ31 τ32  2
 2 2 τ222 τ22 τ32
2 τ32  (22)
2 22 33 22 n τ1 n τ1 + n τ12 n3 τ31 
2 22 33 23 3 2 n τ2 n τ2 + n τ22 n τ2 

101(16)

whereRp A1 C13   X1  

C23   X2 = µ  A2 



C33
X3A3 E1 
B12 B13   D1  

B22 B23   .D2  +  E2  (20)

B32 B33
D3E3C12C13== 0 0 0 n2 τ1
1 n2 τ120
 0
 0 0 n3 τ1
1
 n3 τ12
0000n2 τ21n2 τ220000n2 τ31
n2 τ320000000000000000n3 τ21
n3 τ220000n3 τ31
n3 τ320000000000000


0 

0 

0 
 , (23)
0 

0 

0

0 

0 

0 
 , (24)
0 

0 
0

Proceedings of ITC/ISHW2007

A1=

D1=

B11=

B12=

B13=

Ei=

Ei6=

 ∂2 u1  ∂α1 ∂α1  ∂2 u1  ∂α1 ∂α2  ∂2 u1  ∂α2 ∂α2  , 0

 0

0
  ∂u1  µ ∂x1   ∂u1  µ ∂x2   ∂u1  µ 3  ∂x  ,
 0


 0



0
 ∂τ1

∂τ21
∂τ31
 11

000
11
∂α
∂α

 ∂α
∂τ11
∂τ21
∂τ31

 200022

∂α2
∂α3
 ∂α1
∂τ
∂τ
∂τ

22
 2000 , ∂α2 1 ∂α22 ∂α32 2 ∂n1 ∂n1 ∂n1 0 0 0 ∂α2∂α3 ∂α∂n1∂n∂n 2 ∂α0 0 0 
222
∂α
∂α

0000 0 0

0 0 0 0 0  0 00 0 0 0 0 

 00 0 0 0 0  ∂n2
 1 0 0 0 0 0  ,∂α ∂n2 2 0 0 0 0 0  ∂α
00 0 0 0 0

0 0 0 0 0 
 0
 00 0 0 0 0 
 00 0 0 0 0  ∂n3
 1 0 0 0 0 0  , ∂α3 ∂n 2 0 0 0 0 0 ∂α
00 0 0 0 0
 0  0 
 0  ∂ri  , u1  ∂αi  ∂ru  ∂α2 Ei6∂p∂ui∂ui+ ρ+ u1 ρ 1i
∂x
∂t
∂xii
∂u
∂u+u2 ρ 2 + u3 ρ 3 .
∂x∂x

conditions for many variables and coupled combinationsof the variables must be formulated.
(25)(26)(27)(28)(29)(30)(31)

The other block matrices
 are also derived easily. The jumpconditions for3

∂2 p
∂x j ∂xk

are given in similar expressions.

Concluding Remarks

We derived the jump conditions in multiphase NavierStokes flows with singular forces, the phases having different mass densities and viscosities. Numerical computations utilizing these formulations will be shown in ourfuture work. To make the IIM method applicable to MHD
simulations or other fluid models of fusion plasmas, jump

[1] H. Miura et al., J. Plasma Physics 72, 1095 (2006).
[2] H. Miura et al., Fusion Science and Technology 51, 8 (2007).
[3] T.-H. Watanabe et al., Nucl. Fusion 47, 1383 (2007).
[4] S. Xu and Z. Jane Wang, SIAM J. Sci. Comput. 27, 1948
(2006).
[5] S. Xu and Z. Jane Wang, Comput. Methods Appl. Mech.
(2007), doj:10.10.16/j.cma.2007.06.012.
[6] C.S. Peskin, J. Comput. Phys. 25, 220 (1977).