A Level Set Approach to Eulerian-Lagrangian
Coupling
Marco Arienti, Patrick Hung, Eric Morano and Joseph E. Shepherd
Graduate Aeronautical Laboratories, MS 205-45, California Institute of Technology,
Pasadena, CA 91125 USA.
Received June 2001
DRAFT
September 9, 2002, 11:39am
DRAFT
2
ARIENTI, HUNG, MORANO AND SHEPHERD
We present a numerical method for coupling an Eulerian compressible
∞ow solver with a Lagrangian solver for fast transient problems involv-
ing ∞uid-solid interactions. Such coupling needs arise when either spe-
ciØc solution methods or accuracy considerations necessitate that diÆerent
and disjoint subdomains be treated with diÆerent (Eulerian or Lagrangian)
schemes.
The algorithm we propose employs standard integration of the Eulerian
solution over a Cartesian mesh. To treat the irregular boundary cells that
are generated by an arbitrary boundary on a structured grid, the Eulerian
computational domain is augmented by a thin layer of Cartesian ghost cells.
Boundary conditions at these cells are established by enforcing conservation
of mass and continuity of the stress tensor in the direction normal to the
boundary. The description and the kinematic constraints of the Eulerian
boundary rely on the unstructured Lagrangian mesh. The Lagrangian mesh
evolves concurrently, driven by the traction boundary conditions imposed
by the Eulerian counterpart.
Several numerical tests designed to measure the rate of convergence and
accuracy of the coupling algorithm are presented as well. General problems
in one and two dimensions are considered, including a test consisting of an
isotropic elastic solid and a compressible ∞uid in a fully coupled setting
where the exact solution is available.
Keywords:
65N99 Partial diÆerential equations, boundary value prob-
lems; 74F10 Fluid-solid interaction; 76L05 Shock
waves and
blast
waves.
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
3
1. INTRODUCTION
There is a variety of numerical methods today that can be used to tackle multi-
material mechanics problems involving both ∞uid and solid motions as well as their
interactions. Traditionally, ∞uid dynamicists have favored Eulerian methods while
solid mechanicians prefer Lagrangian methods.
It is clear that a large class of applications is neither ideally suited for a pure
Lagrangian nor a pure Eulerian approach. In a Lagrangian calculation, the mesh
points correspond to elements of mass in the material and their trajectories follow
the particle paths of the material elements. Thus, the position of a boundary is au-
tomatically calculated. Since the initial accuracy of the approximation is generally
maintained throughout the computation, Lagrangian schemes have proven to be
very accurate for a constant number of mesh points as long as the approximating
mesh remains regular. However, as the computation evolves, the stretching of the
computational grid may drastically reduce the stable time step (which is propor-
tional to the minimum size of an element of a triangular grid, say). Furthermore, if
the mesh becomes highly distorted, the calculation becomes increasingly inaccurate.
Remeshing is then required, at the price of increased complexity and computational
eÆort. Conversely, a pure Eulerian calculation allows the development of a com-
plex ∞ow, at the price of a loss in accuracy when treating an arbitrarily varying,
time-dependent boundary (which is intrinsically of a Lagrangian type). Examples
of the limitations of a pure Eulerian or Lagrangian approach can be found in [1].
Hybrid methods such as the ALE (Arbitrary Lagrangian Eulerian) method ad-
dress multi-material applications involving ∞uid-solid interactions. References [2, 3,
4, 5] report on recent advances in ALE coupling. In ALE methods, the conservation
equations are expressed in a control volume formulation, where the control volume
is bounded by a surface
S
a
(
t
) moving with arbitrary local velocity
u
a
. The control
volume velocity
u
a
has two important limit values. In the Eulerian limit,
u
a
=
0
DRAFT
September 9, 2002, 11:39am
DRAFT
4
ARIENTI, HUNG, MORANO AND SHEPHERD
and the control volumes are Øxed. In the Lagrangian limit
u
a
=
u
, the local ∞ow
Øeld velocity and the control volumes coincide with material volumes.
In ALE methods, a Lagrangian phase is followed by a coordinate transforma-
tion due to the mesh motion (remapping or advection phase). The advection step
performs an incremental rezone in which nodes are moved only a small fraction of
a typical length of the surrounding elements. Monotonic advection algorithms are
used to prevent the advection step from creating new minimum or maximum values
for the solution variables.
In ∞uid-structure interactions with ALE methods, the equations of the structural
elements are usually expressed using a purely Lagrangian scheme, so that the nodes
follow the motion of material particles. The interaction with a ∞uid can be modeled
through intermediate regions in which the mesh moves with a spatially varying
velocity. A grid-rezoning technique is used within the bulk of the ∞uid domain
to respect the movement of Lagrangian interfaces by simultaneously minimizing
the grid distortion. The speciØcation of
u
a
is key to the success of ALE methods.
Unfortunately, this process often requires
a priori
knowledge of the solution when
modeling the problem. For systems where the Lagrangian domain suÆers large
deformation or where the Eulerian ∞ow has high rotations, ALE methods will often
fail to give a solution. Finally, the ALE method does not appear to be very suitable
for loosely coupling separate Eulerian and Lagrangian software packages since this
introduces a third solution algorithm and increases the complexity of the coupling
process rather than simplifying it.
The scheme we propose provides an alternative to ALE methods. Unlike ALE
schemes, there is no mixed Eulerian-Lagrangian region and coupling occurs only
through interaction at the boundary of the Eulerian and Lagrangian regions. The
two subdomains are integrated separately by two independent (synchronized) Eu-
lerian and Lagrangian solvers. We will refer to this scheme using the acronym GEL
(Ghost-∞uid Eulerian-Lagrangian).
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
5
Since the Eulerian region mesh is usually Cartesian and stationary, the problem
of how to treat moving and irregular ∞ow domain boundaries arises. When an ar-
bitrary boundary is embedded in a structured grid, any cell whose interior contains
this boundary will be a \cut cell". By Ølling in the cut cells and a small layer
of neighboring regular cells with an appropriate \ghost ∞uid", cell updates can be
performed in the same standard way as in the bulk of the computational domain.
Additionally, the time step is not constrained by the geometry of the cut cells.
The treatment of the cut cells is crucial to any Cartesian grid method, since a
straightforward approach would reduce the stable time step to an arbitrarily small
value due to the reduced size of the irregular cells. Also, accuracy and conservation
at the boundary must be addressed, but in such a way that extension of the coupling
scheme to two and three dimensions is relatively easy. Early work by Noh [1]
made use of redistribution and cell-merging, probably the most popular approach
in Cartesian grid methods. An extensive review of these methods can be found in
[6]. Numerical results typically present only Ørst-order accuracy at the boundary,
independent of the accuracy of the ∞ow solver in the bulk of the computational
domain.
Pember
et al.
[6] propose an adaptive Cartesian grid method where the boundary
is reconstructed with a shock tracking approach and treated as a stationary re∞ect-
ing wall. The coupling scheme is an explicit two-step method, enforcing conserva-
tion in the cut cells and in the neighboring regular cells through a redistribution
algorithm. Numerical results for a Prandtl-Meyer expansion show a degradation of
the accuracy of the scheme to Ørst-order at the boundary. Recently, Falcovitz
et al.
[7] proposed a coupling scheme that maintains conservation across the boundary
cells when both the ∞uid and the boundary undergo uniform motion. No indication
is given on how the algorithm performs in actual dynamical problems.
The method we propose has its origin (and name) in the Ghost Fluid Method
(GFM) originated by Fedkiw
et al.
[8, 9, 10]. GFM originated as an algorithm
DRAFT
September 9, 2002, 11:39am
DRAFT
6
ARIENTI, HUNG, MORANO AND SHEPHERD
for handling Eulerian multi-phase multi-∞uid problems where interfaces separate
regions of diÆerent ∞uids, e.g., air and water. The original GFM is designed to
capture discontinuous interfaces with an Eulerian solver on each side. Within a
prescribed distance of an interface, an Eulerian grid point is a real node to one
solver and a ghost node to the other. The prescription for populating a state of
a ghost node in the GFM is to replace pressure and normal velocity from the real
node, while extrapolating in the normal direction a second thermodynamic variable
(entropy) and the tangential velocity.
GEL diÆers from the original formulation in the way the solution variables in
the ghost region are populated. In this respect, GEL is more akin to the local
mirroring extrapolation technique presented by Forrer and Berger [11] and Forrer
and Jeltsch [12], since it treats the Eulerian-Lagrangian interface as an impermeable
wall at any new iteration of the Eulerian solver. The main diÆerence with respect
to [11, 12] lies in the way the Eulerian-Lagrangian interface is tracked; in our case, a
level set-based approach [13, 14] is followed. Also, GEL does not use boundary cell
averaging and the order of the mirroring extrapolation is lower than in [11, 12], at
least for a smooth ∞ow. While less sophisticated, GEL is expected to be more robust
when dealing with arbitrarily complex boundaries and ∞uid-solid shock interactions,
particularly in three dimensions. We note that results presented by Forrer and
Berger for moving boundaries are limited to ∞ow interactions with rigid bodies,
whereas we address fully coupled problems in the sense described by Noh [1].
In Section 1, a thorough treatment of the coupling algorithm is provided. Section
2 is an overview, stating the class of problems we intend to address and outlining
the coupling algorithm. As the Eulerian and Lagrangian solvers that are used in
this paper are well established, their descriptions are brief by intent, with further
references provided for readers desiring more information. Section 3 focuses on the
decomposition of the solution domain into Eulerian and Lagrangian subdomains,
which leads to the deØnition of the Eulerian-Lagrangian interface. The dynamical
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
7
and kinematic constraints that need to be satisØed by this interface are described
there. Section 4 addresses the numerical implementation of the coupling scheme.
Other pertinent issues, including the temporal evolution of the subdomains, possible
enhancements, and areas of future research are also discussed.
The second part of the paper presents several numerical tests designed to measure
the rate of convergence and accuracy of the algorithm. Section 5 describes several
one-dimensional veriØcation tests. The behavior of the interface under shock
wave
transmission is examined in the \transparency tests". Section 6 discusses several
two-dimensional tests. VeriØcation tests proposed in [11] are adopted in order to
investigate issues of mass conservation. The Ønal example we consider is a veriØ-
cation test consisting of a shock load over an isotropic elastic solid (superseismic
loading problem). This is a fully coupled ∞uid-solid interaction problem where the
exact solution is known.
2. OVERVIEW
The work on Eulerian-Lagrangian coupling described in this paper is a part of
the research eÆort at Caltech to develop a Virtual Test Facility (VTF) that pro-
vides a problem-solving environment for three-dimensional parallel simulations of
the dynamic response of materials in multi-physics problems. GEL coupling is used
extensively and is the workhorse behind the VTF [15]. A typical example applica-
tion of the VTF is modeling the cylinder test used for high-explosive performance
studies [16]. In this test, sketched in the left panel of Fig. 1, a cylinder of high
explosive is detonated inside of a metal tube. The subsequent motion of the metal
is used as a measure of the performance of the explosive. The motion of the metal
and the progress of the detonation
wave can be
strongly coupled if the chemical
reaction processes in the detonation
wave are
su±ciently slow.
The approach we follow divides the problem into a ∞uid mechanics and a solid
mechanics portion (see Fig. 1). This division is both theoretical and computational.
DRAFT
September 9, 2002, 11:39am
DRAFT
8
ARIENTI, HUNG, MORANO AND SHEPHERD
Metallic confinment
Reaction products
HE
reaction zone
stress fronts
FIG. 1.
Cylinder test and domain decomposition. The cylinder test problem (left) is de-
composed into a Cartesian domain, here with Adaptive Mesh ReØnement, for the Fluid Mechanics
solver (center); and a Lagrangian domain for the Solid Mechanics solver (right).
On the one hand, it allows research and development on the two diÆerent subject
areas to progress separately and concurrently; on the other hand, this separation
lets us combine the strengths of Eulerian and Lagrangian solvers.
Following standard practice in gas dynamics, shock physics, and high explosives
modeling, we approximate the motion in the Eulerian region as hydrodynamic
in nature and neglect viscous eÆects. We discuss only nonreactive ∞ow in this
paper but that is not an essential limitation and the method has been applied to
detonation problems. Nonreactive, inviscid ∞uid dynamics is described by the Euler
equations [17]
@Ω
@t
+
r¢
(
Ω
u
)= 0
;
@
@t
(
Ω
u
)+
r¢
(
Ω
uu
+
I
P
)= 0
;
@
@t
∑
Ω
μ
e
+
u
2
2
∂∏
+
r¢
∑
Ω
u
μ
e
+
P
Ω
+
u
2
2
∂∏
=0
:
(1)
Here
Ω
is the density,
u
the velocity vector,
P
the pressure, and
e
the speciØc
internal energy. These equations may be rewritten in conservative form
U
t
+
r¢
F
(
U
)=0
;
(2)
by using the vector of conserved variables
U
=(
Ω; Ω
u
;ΩE
)
T
, the ∞ux vector
F
(
U
)=(
Ω
u
;Ω
uu
+
I
P; Ω
u
(
E
+
P=Ω
))
T
, and the speciØc total energy
E
=
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
9
e
+1
=
2
k
u
k
2
. The equation set (1) must be supplemented by an Equation of State
(EoS) that describes the equilibrium thermodynamic state of the material. For our
purposes, it is su±cient to have a relationship between pressure, internal energy,
and mass density,
P
=
P
(
e; Ω
). The examples presented in the paper all use the
approximation of a perfect gas with a constant ratio of speciØc heats
∞
,
P
=(
∞
°
1)
Ωe :
(3)
Examples of a multi-species reactive ∞ow with the Mie-Grƒ
uneisen EoS for the solid
explosive and the Jones-Wilkins-Lee EoS for the reaction products are presented in
[15].
A number of simple Lagrangian solvers were developed to verify and validate
Eulerian-Lagrangian coupling. These solvers are a one- and two-dimensional in-
tegrator of motion for rigid bodies with assigned trajectory or traction boundary
conditions, an Euler-Bernoulli beam module, a one-dimensional Lagrangian gas
dynamics solver, and a two-dimensional explicit Finite Element (FE) solver. Al-
though these Lagrangian solvers are rather simple, the methods described in this
paper have also been applied [15] with a rather sophisticated [18] FE solver.
Each one of these discretizations leads to a set of ordinary diÆerential equations
in the generalized displacements
X
of the form
M
d
2
X
dt
2
+C
d
X
dt
+K
X
=
f
:
(4)
We recognize in Eq. (4) the mass matrix M, the stiÆness matrix K, the dissipation
matrix C, and the vector
f
of generalized external forces.
In the rigid body module,
X
is the vector of coordinates of the center of mass. We
have C = 0, K = 0 and
f
(
t
) is computed by integrating the Eulerian pressure Øeld
over the boundary. The solution is advanced in time either by a third-order Runge-
Kutta or an explicit Newmark scheme. Alternatively, what is called an essential
boundary condition,
i.e.
a speciØcation of the trajectory of the body, can be used.
DRAFT
September 9, 2002, 11:39am
DRAFT
10
ARIENTI, HUNG, MORANO AND SHEPHERD
The beam module is based on the classical Euler-Bernoulli (sometimes known
as Bernoulli-Euler, or Coulomb) beam theory, attributing the resistance to ∞exure
entirely to extension and contraction of longitudinal Ølaments [19]. See, for exam-
ple, [20] for a list of the theory's main assumptions and limitations. The derivation
of the element matrices is standard and will not be reproduced here. A consistent
mass matrix is used for the integration of the transient beam equation. See [21, 22]
and [23] for further details.
A one-dimensional Lagrangian gas dynamics module has been developed to con-
duct a class of veriØcation problems referred to as transparency tests. The user can
specify any combination of linear (von-Neumann) artiØcial viscosity [24], quadratic
artiØcial viscosity [25, 26], and artiØcial heat conduction [27] for solving the Euler
equations with a perfect gas EoS. Nodal variables are nodal displacements and their
derivatives. Cell variables are pressure, internal energy, speciØc volume, artiØcial
viscosity, and heat conduction. Properties of the Ønite elements (or cell variables)
are staggered spatially with respect to the nodal variables (in this one-dimensional
setting). Cell variables and nodal variables are also staggered temporally as shown
in Fig. 2. Initial conditions corresponding to
n
= 0 are displacement and velocity
at
t
=
°
0
:
5 and pressure and speciØc volume at
t
= 0. Pressure at
t
= 0 is used
to update velocity to
t
=0
:
5 by the momentum equation. Velocity at
t
=0
:
5is
used to update speciØc volume to
t
=1
:
0 by the continuity equation. The energy
equation and EoS are used to compute the pressure at
t
=1
:
0. This completes one
cycle, and
n
is now at 1. The Ønite diÆerence numerical method follows essentially
the scheme of von-Neumann and Richtmyer [24]. For a more recent treatment of
artiØcial viscosity and a detailed description of its implementation in two and three
dimensions, see [28].
The two-dimensional FE module is a displacement-based solver used for com-
puting plane stress/plane strain problems. Linear (3 nodes) or quadratic (6 nodes)
triangular plane elements, with linearized kinematics and explicit time integration
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
11
FIG. 2.
Finite diÆerence space-time discretization for one-dimensional Lagrangian solver.
of a Hookean elastic material, are implemented. The solution is advanced by evalu-
ating the force balance at time
t
n
in order to compute the acceleration of the system
between
t
n
and
t
n
+1
. This is equivalent to computing the velocity and acceleration
via central diÆerences, and the method is, therefore, second-order accurate
O
(¢
t
2
)
when equal time steps are used. The advantage of this scheme is low computational
cost and low storage when a diagonal mass matrix is used [23].
In this paper, we restrict our attention to the class of problems involving the
coupling of an inviscid ∞uid with a solid whose interface is assumed to be imperme-
able, nonreactive, adiabatic, and unable to support surface tension. The VTF, for
which the method is designed, targets problems (e.g., the cylinder test) which fall
under this category. With these assumptions, the Eulerian-Lagrangian interface
(EL-interface, to be deØned in the following Section) can be treated as a contact
discontinuity with the following properties: 1) No mass ∞ux; 2) No jump in normal
velocity; 3) Free-slip boundary condition for the tangential velocity; 4) No jump in
the normal stress. Jumps in entropy (or density) across the interface are admitted.
These properties are enforced by the coupling scheme through the application of
boundary conditions at discrete times, as will be seen in the next Section.
DRAFT
September 9, 2002, 11:39am
DRAFT
12
ARIENTI, HUNG, MORANO AND SHEPHERD
3. TIME AND SPACE DISCRETIZATION
In our procedure, the entire solution domain is decomposed into subdomains of
two diÆerent types, Eulerian and Lagrangian. We assume that coupling exists only
between subdomains that share a boundary, and only through boundary conditions.
The GEL coupling scheme provides these boundary conditions to the solvers in a
manner to be described next. The cases of Eulerian-Eulerian coupling (e.g., multi-
∞uid simulations) and Lagrangian-Lagrangian coupling (e.g., contact mechanics)
will not be discussed but their importance is noted.
3.1. Spatial Discretization and Interface Representation
Given an IBVP (Initial Boundary Value Problem) on ≠
£
T
, the domain ≠ (and
the time coordinate on the interval
T
, treated in the next Section) needs to be
discretized for numerical solution. Loosely speaking, part of the domain will be
covered by a Lagrangian mesh and the rest by an Eulerian mesh.
The discretization ≠
L
of ≠ that is associated with the Lagrangian solver is the
Lagrangian domain
. Generally, ≠
L
is an unstructured grid. The discretization ≠
E
that is associated with the Eulerian solver is the
Eulerian domain
. In this paper, we
further specialize ≠
E
to be a collection of structured grids, namely Cartesian grids.
Note that, in general, ≠
L
Ω
≠ and ≠
E
Ω
≠. Figure 3 is a sketch of a Lagrangian
domain made of triangular elements which is partially superimposed on a Cartesian
grid.
The boundary representation of ≠
L
is an oriented surface denoted by
@
≠
L
.We
are particularly interested in the subset of
@
≠
L
whose points lie in ≠
E
. We call
this subset the
EL-interface
@
≠
EL
. All the coupling that takes place between the
Eulerian domain and the Lagrangian domain is assumed to occur at this interface.
The level set function (denoted by the scalar Øeld
'
) contains an implicit repre-
sentation of
@
≠
EL
on ≠
E
. It is deØned as the signed distance from
@
≠
EL
, evaluated
at the center of a cell (in a Ønite volume scheme) or at a grid vertex (in a Ønite
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
13
diÆerence scheme). In this paper, we will refer to such a point as a
node
when
the distinction between vertex and cell center is unimportant. The level set is
discretized as the set of points
'
i; j
=
'
(
x
i
;y
j
), where the indices
i
,
j
span the
Cartesian grid. The level set was originally applied to computations on Cartesian
grids by Osher and Sethian in [29] and has been successfully employed to resolve
sharp interfaces between materials with diÆerent properties or diÆerent equations
of state [8].
The equation
'
= 0, or the zero level of the discretized Øeld
'
i; j
, identiØes
the EL-interface
@
≠
EL
. At each time step, the level set is reconstructed on the
Cartesian grid from the current Lagrangian description of the boundary. To this
end, ray-intersection, a popular approach for determining whether a point lies inside
or outside of a surface (a process known as point-classiØcation), is implemented [30].
The complexity of the overall algorithm is of order
O
(
M
¢
N
) where
M
is the number
of points in the Lagrangian boundary representation and
N
is the number of points
in the Eulerian grid. Recently, an algorithm to reconstruct (in three dimensions)
the closest point transform
'
with optimal complexity has been developed [31].
Throughout this paper, we will use the sign convention introduced in [9],
≠
R
=
f
≠
E
:
'
i; j
∑
0
g
:
(5)
The set ≠
R
designates the real (or ∞ow) part of the Eulerian domain, where the
∞ow Øeld is computed, as opposed to the ghost region, where boundary conditions
are set. Equation (5) means that the level set in the real part of the Eulerian
domain is negative. Since
'
deØnes contour levels of the signed distance function,
the gradient
r
'
at
@
≠
EL
must be perpendicular to
@
≠
EL
itself. A corollary of the
sign convention is that the normal vector
r
'=
kr
'
k
is oriented from the Eulerian
to the Lagrangian domain. In practice,
r
'=
kr
'
k
is numerically approximated by
central diÆerencing of
'
i; j
.
DRAFT
September 9, 2002, 11:39am
DRAFT
14
ARIENTI, HUNG, MORANO AND SHEPHERD
FIG. 3.
Eulerian (regular grid) and Lagrangian discretizations (triangular mesh) showing
overlap of domains and a layer of ghost cells (Ølled points). The EL-interface is indicated by the
thicker line on the Lagrangian boundary.
To use regular Cartesian grid cells, the collection of cells ≠
R
needs to be aug-
mented so that each computational node has a complete stencil. This is done by
adding to ≠
R
a thin layer of cells ≠
G
Ω
≠
E
at the EL-interface. We call this subset
the
ghost region
. The Eulerian solver operates exclusively on ≠
R
[
≠
G
μ
≠
E
.In
Fig. 3, these cells are marked by a Ølled circle. The ghost region ≠
G
is required to
be \big enough" and the required size depends on the details of the Eulerian solver.
A more precise statement will be made in the next Section.
3.2. Ghost cells
The concept of ghost cells as used here is an extension for arbitrary boundaries of
the commonplace guard (or ghost) cells that surround a computational patch,
i.e.
,
a rectangular Cartesian grid. This practice allows for the application of boundary
conditions (BC). Additional rows of guard cells are used to \complete" the stencil
of the external cells of a patch so that the solver does not need to be aware of the
boundary of its computational domain. The actual number of rows depends only
on the stencil of the Eulerian scheme.
The introduction of a level set function allows for the extension of this idea to
an arbitrary boundary. The Eulerian node (
i; j
)isa
ghost node
if the level set at
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
15
that point satisØes
≠
G
=
f
≠
E
:0
<'
i; j
∑
'
S
g
:
(6)
The corresponding cell can be partially or fully covered by the Lagrangian dis-
cretization ≠
L
.
The parameter
'
S
depends exclusively on the stencil of the numerical scheme that
is used to compute the ∞uxes. For example, a second-order accurate Essentially-
Non-Oscillatory (ENO) scheme requires four nodes, two on each side of the node
that is being updated, to compute the numerical ∞ux. The extent of the ghost
region should be at least 2
¢
¢
x
, since we cannot expect a cell boundary to lie
exactly at
'
= 0. We must account for the motion of the EL-interface because a
ghost cell can be \exposed" and become a real cell at the next time step. Thus,
in this example, the buÆer area must be increased at least to
'
S
=3
¢
¢
x
. This
one-dimensional argument can be extended to higher dimensions as long as the
solver implements a dimension-by-dimension integration.
According to Eq. (6), the ghost region has to be initialized (or populated) up to
a distance
'
S
; the details will be discussed in the next Section. Note, however, that
the coupling procedure is completely independent of the patch integrator. We can
think of it as setting the proper boundary conditions before advancing the solution
by one time step. Indeed, given a generic patch integrator, we require only two
additions to the code: 1) A test to compute the numerical ∞ux only if
'
∑
'
S
;2)
A test to update the solution only if
'<
1
¢
¢
x
. The second point above is important
because a ghost cell within one ¢
x
from the EL-interface can become a real cell
after a time step. No more than one ghost cell layer needs to be updated since
the Courant-Friedrichs-Levy (CFL) condition is applied over all ≠
R
[
≠
G
when
estimating a stable time step for the ∞uid solver. The CFL condition prevents the
contact discontinuity at the EL-interface from sweeping more than a fraction of a
Cartesian cell when the solution is advanced.
DRAFT
September 9, 2002, 11:39am
DRAFT
16
ARIENTI, HUNG, MORANO AND SHEPHERD
3.3. Time Discretization and Temporal Coupling
The time coordinate
T
=[
t
i
;t
f
] is discretized (partitioned) by
μ
=
f
t
0
;t
1
;t
2
;:::;t
n
g
as follows
T
=
N
[
1
ø
n
(7)
where
ø
n
=(
t
n
°
1
;t
n
].
It is clear that
t
0
=
t
i
(the initial time) and
t
n
=
t
f
(the Ønal time). The set
μ
represents instants in time when the
coupling is performed
. Note that this can be
diÆerent from the temporal discretizations used for the Eulerian or the Lagrangian
solvers. In fact, Eulerian solvers and Lagrangian solvers usually employ diÆerent
time integrators, possibly with multiple steps.
The current implementation of GEL deØnes the common time increment ¢
t
which is allowable from stability considerations for the Eulerian and Lagrangian
grids, say ¢
t
E
and ¢
t
L
. This increment can be computed at run time from the
previous cycle of computation,
i.e.
,¢
t
n
+1
= min (¢
t
n
E
;
¢
t
n
L
). However, Noh [1]
observes that several applications typically have ¢
t
L
ø
¢
t
E
;
and allows ¢
t
E
to
be a (stable) multiple of ¢
t
L
.
To simplify the discussion on temporal coupling, we will refer to the set of equa-
tions integrated by the Eulerian solver as
E
. Similarly,
L
refers to the set of equa-
tions integrated by the Lagrangian solver. The solution of
E
depends on boundary
conditions provided by
L
, which depends on the state of
L
, and vice versa. The
solution is assumed to be available at discrete times up to
t
n
.
One can advance the solution by integrating
E
using
L
(
t
n
) and by integrating
L
using
E
(
t
n
); this simple approach is called
concurrent time coupling
and it is used
in the VTF as it is more suitable for solution by parallel computers than staggered
methods. An example of a staggered method is as follows:
E
is integrated using
L
(
t
n
) and
L
is subsequently integrated using
E
(
t
n
+1
). With staggered methods,
only one set of equations can be solved at any given time.
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
17
A viable alternative is called
PC-Heun time coupling
[32]. After the Ørst integra-
tion, one can re-integrate
E
using a combination (average) of
L
(
t
n
) and
L
(
t
n
+1
)to
get a new
E
(
t
n
+1
) and, concurrently, re-integrate
L
using a combination (average)
of
E
(
t
n
) and
E
(
t
n
+1
) to get a new
L
(
t
n
+1
). Using a predictor-corrector scheme
such as this is often not practical because of its high overhead in both CPU and
memory requirements. Numerical experiments using this scheme are presented in
Section 5.2.
4. COUPLING SCHEME
GEL is a boundary condition coupling scheme for Eulerian and Lagrangian solvers
which are sharing portions of their boundaries (their EL-interface). In this Section,
we will describe GEL, and variations of it, in detail.
In formulating the boundary condition exchange, we make the following assump-
tions: 1) The EL-interface is deØned by the boundary as geometrically determined
by the Lagrangian solver. It is identiØed by
'
= 0; 2) the Lagrangian solver uses a
natural (pressure) boundary condition at the EL-interface; 3) The Eulerian solver
requires a no-∞ux boundary condition (with free-slip) at the EL-interface.
A consequence of the Ørst assumption is that, as stated earlier, the EL-interface
location is recomputed as the Lagrangian boundary moves.
As for the second assumption, either displacements or force boundary conditions
could be applied to the boundary of a Lagrangian solver. The second assumption
indicates that only a force boundary condition is used. In our implementation, pres-
sure is linearly interpolated in the Eulerian domain at the location of the Lagrangian
pressure control points and used to enforce the traction boundary condition. It is
interesting to note that applying a velocity boundary condition by using the Eule-
rian velocity at the boundary does not work. This can be seen with the following
one-dimensional experiment. Imagine a shock
wave in ≠
E
traveling towards an
initially stationary solid ≠
L
. Since the solid is initially stationary, it will act as a
DRAFT
September 9, 2002, 11:39am
DRAFT
18
ARIENTI, HUNG, MORANO AND SHEPHERD
re∞ective boundary to the Eulerian solver at the Ørst integration, but this implies
that the Eulerian ∞ow velocity at the boundary remains zero also after marching
by one time step. This information is fed back to the solid, which, therefore, re-
mains still. Thus, there is no shock transmission in the Lagrangian domain when
we would expect one.
The last assumption implies that the Eulerian solver sees the domain computed
by the Lagrangian solver only as a moving boundary (completely ignoring all states
in the interior of the Lagrangian domain). The only information needed from ≠
L
is the velocity vector evaluated at the boundary. This conclusion is similar to the
one obtained in [33] for coupling compressible to incompressible ∞ow.
4.1. Populating the Ghost Cells: Fluid Solver Boundary Condition
Population of a ghost cell states requires an extrapolation algorithm operating
on the real ∞ow. In the following, the state that is extrapolated from ≠
R
is marked
with the subscript
E
, whereas the state that is evaluated from
@
≠
EL
is denoted by
W
. The populated state in ≠
G
is denoted by
G
.
The level set
'
introduces a vector Øeld of normals
n
=
r
'=
kr
'
k
. We deØne
t
to be the unit vector normal to
n
. The projection of a ghost node over ≠
EL
is
x
W
=
x
G
(
i; j
)+
'
(
i; j
)
n
(8)
and the corresponding boundary velocity
V
W
can be found by interpolation of
Lagrangian boundary values at
x
W
.
For the extrapolation algorithm, the approach of capturing the EL-interface as
a contact discontinuity (proposed in [10] and [9]) suggests that pressure and nor-
mal velocity should be continuous across the interface. However, the choice of the
extrapolation scheme is, in general, not unique. In this paper, we experimented
with diÆerent algorithms for extrapolating density, pressure, and the ∞ow Øeld ve-
locity. These techniques are described in the following Sections, and they consist of
one-sided constant extrapolation (injection), a variation on constant extrapolation
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
19
(re∞ection), and linear extrapolation (mirroring). Results for each of these schemes
can be found in the second part of this paper for one- and two-dimensional prob-
lems. The current Section is closed by a few considerations on the implementation
of the boundary conditions as a Riemann problem.
4.1.1. Constant Extrapolation or Injection
Constant extrapolation of a scalar quantity
I
can be achieved through advection
by integrating the Eikonal equation
I
t
+
n
¢r
I
=0
;
(9)
subject to the boundary condition
I
=
I
W
on
@
≠
EL
:
The numerical discretization of Eq. (9) can be implemented as a Ørst-order up-
wind space discretization with Ørst-order accurate time integration [34]. The equa-
tion has to be solved for a number of pseudo time steps until the ghost region has
been fully populated [9, 10]. Our experience is that this scheme is robust even
for irregular interfaces and shock interactions, since sharp variations in the ad-
vected quantities are smoothed by the Ørst-order advection algorithm. To reduce
computational cost, advection needs to be performed only on a tiny strip of the
computational domain, enclosing the ghost region and the closest strip of real ∞ow
region (say
°
¢
x
∑
'
∑
'
s
).
The prescription for populating a ghost node is
0
B
B
B
B
B
@
Ω
G
V
G
P
G
1
C
C
C
C
C
A
=
0
B
B
B
B
B
@
Ω
E
(
V
W
¢
n
)
n
+
V
E
°
(
V
E
¢
n
)
n
P
E
1
C
C
C
C
C
A
;
(10)
which assigns the normal velocity of the Lagrangian boundary to the normal velocity
of the ghost node. We call this process the normal velocity treatment by
injection
.
DRAFT
September 9, 2002, 11:39am
DRAFT
20
ARIENTI, HUNG, MORANO AND SHEPHERD
4.1.2. Extrapolation by Re∞ection
The extrapolation scheme by
re∞ection
is a variation of the extrapolation by
injection. In the standard treatment of an impermeable wall, the normal component
of the velocity in a ghost cell is set to be antisymmetric with respect to the normal
velocity in the reference frame moving at (
V
W
¢
n
)
n
. Based on this argument, we
propose an alternate formula for
V
G
0
B
B
B
B
B
@
Ω
G
V
G
P
G
1
C
C
C
C
C
A
=
0
B
B
B
B
B
@
Ω
E
(2
V
W
¢
n
°
V
E
¢
n
)
n
+(
V
E
¢
t
)
t
P
E
1
C
C
C
C
C
A
:
(11)
The normal velocity treatment by re∞ection relies on advection as the previous
treatment. It is not a linear extrapolation because it does not account for the
distance from the wall to the ghost node. An advantage over linear extrapolation
is that this prescription avoids the overshoot of an extrapolated velocity when the
interface is very close to a real node. An advantage over injection is that the
antisymmetric treatment of normal velocity is a more accurate implementation of
an impermeable wall.
4.1.3. Linear Extrapolation or Mirroring
Linear extrapolation in the ghost region can be performed by computing the
mirror image
^
x
in ≠
R
of the (
i; j
) ghost node,
x
G
[11, 35]. This operation requires
the normal distance of
x
G
from
@
≠
EL
, which is (with the correct positive sign) the
value
'
(
i; j
). Thus, the relation between
^
x
and
x
G
is provided by
^
x
=
x
G
(
i; j
)+2
'
(
i; j
)
n
:
(12)
The prescription for populating the ghost state is formally identical to the one in
Eq. (11), only now
V
E
,
Ω
E
,
P
E
are interpolated values (we use bilinear interpola-
tion) at
^
x
in ≠
R
.
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
21
It should be noted that extrapolation by mirroring is based on the assumption
that the EL-interface is moving at constant speed within a given time increment; in
general, this is not true. A typical counter-example is an explicit FE solver, where
all nodes assume a constant
acceleration
in any given time step. One can include
acceleration eÆects by imposing a suitable pressure gradient in the ghost region. In
one dimension, if ƒ
x
W
(
t
) is the acceleration of the wall, we would have
P
x
(
x
W
;t
)=
°
Ω
(
x
W
;t
)ƒ
x
W
:
(13)
However, complications arise when the ∞ow is not smooth (a shock, even when
smeared, is associated with high acceleration) and implementation is di±cult in
two and three dimensions. Furthermore, one may argue that since the state on
the EL-interface is only interpolated to Ørst-order at best, it is of little advantage
trying to predict the ∞ow behavior at the interface to better than Ørst-order.
4.1.4. Coupling as a Riemann Problem
In this Section, we suggest an improvement of GEL via formulation of a Riemann
problem at the EL-interface. The implementation requires that the equations of
state of the materials be communicated to a software module which solves the
Riemann problem at the interface. In the simplest case of a linearized solver,
knowledge of the stiÆness (or the impedance) of the material on both sides of the
EL-interface is su±cient. In any case, this approach involves a more complicated
coupling algorithm.
The Riemann problem can be understood as follows. Consider an initial jump
discontinuity at
x
= 0 separating two constant states,
U
(
x; t
=0)=
8
>
<
>
:
U
l
for
x<
0
U
r
for
x>
0
:
(14)
Equation (14) is an example of a shock-tube problem, where a hypothetical \mem-
brane" situated at
x
= 0 separates the two states at
t
= 0. Without loss of
DRAFT
September 9, 2002, 11:39am
DRAFT
22
ARIENTI, HUNG, MORANO AND SHEPHERD
generality, we take ≠
R
to be on the right of the membrane and ≠
L
to be on the
left. If we consider only the subdomain ≠
L
, this becomes a piston problem where
the entire domain is initially a constant state with a boundary condition at
x
=0.
If the original shock-tube problem is to be solved as a coupled Eulerian-Lagrangian
problem, then the Lagrangian scheme solves the piston problem for an interval
dt
subject to a natural (pressure) boundary condition
P
p
.
There is a unique driving pressure
P
p
for the piston problem that corresponds to
the original shock-tube problem. This is the boundary condition that needs to be
applied to the Lagrangian solver if it were to simulate the original system Eq. (14)
correctly. It is natural to ask whether the pressure
P
l
corresponding to the left state
U
l
and the pressure boundary condition
P
p
are equal. Rephrased in the context
of GEL coupling, is the pressure boundary condition to the Lagrangian solver the
pressure of the ∞uid given by the Eulerian solver? The answer is no, in general. In
the case of a shock-tube Ølled with perfect gases having the ratio of speciØc heats
∞
, sound speed
a
and pressure
P
, with subscripts
l
and
r
denoting the left and right
states,
P
l
and
P
p
are related by the shock-tube [17] equation
P
l
P
r
=
P
p
P
r
"
1
°
(
∞
l
°
1)(
a
r
=a
l
)(
P
p
=P
r
°
1)
p
2
∞
r
[2
∞
r
+(
∞
r
+ 1)(
P
p
=P
r
°
1)]
#
°
2
∞
l
=
(
∞
l
°
1)
:
(15)
Thus, when there is a discontinuity in pressure across the EL-interface, the pres-
sure boundary condition (
P
p
) assigned to the Lagrangian solver is not the pressure
given by the Eulerian solver (
P
l
). In practice, large jumps of pressure are often
smeared and
P
l
=P
r
º
1. In this limit,
P
l
is well approximated by
P
p
and solving
the Riemann problem at the interface to calculate
P
p
is not necessary. Furthermore,
these pressure discontinuities do not remain at the interface as they get convected
away
with the ∞ow, and the resulting error is usually small.
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
23
5. ONE-DIMENSIONAL TESTS
In this Section, we will restrict ourselves to the analysis of one-dimensional test
problems. Two-dimensional problems are discussed in Section 6.
The following results are obtained by using an Essentially-Non-Oscillatory Local
Lax-Friedrichs (ENO-LLF) ∞ow solver [36, 37] on a Cartesian grid. A Local Lax-
Friedrichs form [36] is used to avoid entropy violating shocks near sonic points,
making the scheme very robust at the price of perhaps too much dissipation to
accurately capture all the ∞ow features. Time marching is achieved through the
Total Variation Diminishing (TVD) third-order Runge-Kutta method devised by
Shu and Osher [36]. The time step ¢
t
is selected to account for the CFL stability
condition.
We Ørst consider a prescribed boundary motion to examine issues of conservation
of mass and entropy in isentropic ∞ows. The next step is to couple the boundary
motion to the ∞ow in a simple fashion and examine the problem of free expansion
of a piston. The simplicity of this test case allows us to explore diÆerent temporal
couplings and coupling strategies. Then we study the oscillations of a spring-mass
system in a compressible gas. The problem is reducible to a nonlinear oscillator,
and we compare trajectories of the piston in the underdamped and overdamped
cases with the numerical solutions of the corresponding ordinary diÆerential equa-
tion. Finally, we examine
wave in
teractions in a setting where the Lagrangian and
the Eulerian materials have identical acoustic impedance. This test, called the
transparency test, assesses how transparently
waves can be
transmitted between
domains that diÆer only in their numerical modeling. The Lagrangian modules
were described previously in Section 2 of this paper.
5.1. Prescribed Rigid Boundary Motion
Conservation is not guaranteed [7] by a scheme which does not take into account
cut cells and the motion of cell interfaces during a time step ¢
t
. In this Section, the
DRAFT
September 9, 2002, 11:39am
DRAFT
24
ARIENTI, HUNG, MORANO AND SHEPHERD
convergence properties of the treatments by re∞ection and injection are compared.
Results for mirror ∞ow extrapolation are also computed for completeness. In all
cases, we show that the losses in mass and entropy decrease at least linearly with
the mesh size,
i.e.
, with Ørst-order convergence. For this smooth one-dimensional
problem, quadratic reduction of these errors could be achieved by an appropriate
treatment of how boundary cells are cut by the interface [11, 12]. As noted earlier,
these schemes are potentially cumbersome in higher dimensions and they tend to
provide only linear convergence in multi-dimensional simulations involving shocks
[11].
We will consider a perfect gas conØned between two rigid walls at
x
r
=1
:
0 and at
x
l
=0
:
5+
v
l
t
+
a
l
=
2
t
2
. This problem has been previously considered by Forrer and
Berger [11], who demonstrated second-order convergence. Initial conditions are
Ω
(
x;
0) = 1 + 0
:
2 cos (2
º
(
x
°
0
:
5))
;
v
(
x;
0) = 2 (1
°
x
)
v
l
;
P
(
x;
0) =
Ω
(
x;
0)
∞
:
(16)
If the left wall is moving leftward (
v
l
<
0), an expansion takes place and the ∞ow
Øeld is isentropic for all times,
s
(
x; t
)=
s
(
x;
0) =
p
(
x;
0)
Ω
(
x;
0)
∞
=1
:
0. The
numerical value of the entropy
s
can, therefore, be monitored for error analysis.
We denote an initial value with the superscript
i
and a Ønal value with the
superscript
f
. Discrepancies in mass ¢
m
and entropy ¢
s
at the Ønal time are
given by
¢
m
=
μ
P
n
i
C
j
=1
Ω
i
j
j
C
j
j°
P
n
f
C
j
=1
Ω
f
j
j
C
j
j
∂
=
P
n
i
C
j
=1
Ω
i
j
j
C
j
j
;
¢
s
=
P
n
f
C
j
=1
Ø
Ø
Ø
s
f
j
°
1
Ø
Ø
Ø
j
C
j
j
=
P
n
i
C
j
=1
j
C
j
j
;
(17)
where
j
C
j
j
is the length of that part of a computational cell
C
j
that lies in ≠
R
. The
Ønal time
t
f
is chosen so that the walls are located exactly at the interface between
two grid nodes. Since, at that time, there are no cut cells, the two equations above
are exact estimates of the errors of a piecewise constant solution. Thus, there
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
25
are exactly
n
i
C
cells belonging to ≠
R
at time
t
i
and
n
f
C
cells at time
t
f
. We also
introduce a measure of the error of the entropy at the left moving wall
¢
s
W
=
Ø
Ø
Ø
s
f
W
°
1
Ø
Ø
Ø
;
(18)
where
W
is the index of the left boundary cell.
Equipped with these error estimators, we now consider two cases and study the
convergence of the results as the grid spacing
h
decreases. In all the computations,
a Øxed ratio
dt=h
=0
:
32 is used, corresponding to a CFL number approximately
equal to 0
:
6.
Case A. We set the left wall velocity
v
l
=
°
0
:
5 and the acceleration
a
l
=0.
Results for diÆerent grid reØnements at
t
f
=0
:
5 are shown in log-log plots in
Fig. 4 (left column). Reference lines for linear and quadratic convergence are also
displayed.
For all the extrapolation schemes, the convergence rate is linear for ¢
m
and
quadratic for ¢
s
, whereas ¢
s
W
displays an intermediate behavior. Mirroring gives
the best performance in all the three error indicators, particularly the value ¢
m
which is an order of magnitude smaller than the value obtained by using injec-
tion or re∞ection. This result is expected, as the linear extrapolation described in
paragraph 4.1.3 is designed to implement an impermeable boundary for the case of
constant velocity of the interface. However, we notice that mirroring is slower in
achieving linear convergence of ¢
m
and that the rate of decrease of ¢
s
W
is only
linear, and not superlinear, as for the other two extrapolation schemes.
Case B. The left wall velocity is initially zero, while the acceleration is constant
a
l
=
°
2
:
0. Results at
t
f
=0
:
5 are shown in log-log plots in Fig. 4 (right column).
The convergence is again Ørst-order for ¢
m
and second-order for ¢
s
. The error
in entropy at the wall, ¢
s
W
, decreases linearly for re∞ection and mirroring, and
superlinearly for injection. In this situation of constant wall acceleration, the mass
loss is slightly smaller with the re∞ection treatment. This result suggests that
DRAFT
September 9, 2002, 11:39am
DRAFT
26
ARIENTI, HUNG, MORANO AND SHEPHERD
1/h
∆
m
10
-7
10
-6
10
-5
10
-4
10
-3
10
-2
injection
reflection
mirroring
reference linear convergence
200
100
400
800
1600
3200
1/h
∆
m
10
-6
10
-5
10
-4
10
-3
10
-2
10
-1
injection
reflection
mirroring
reference linear convergence
200
100
400
800
1600
3200
1/h
∆
s
10
-9
10
-8
10
-7
10
-6
10
-5
10
-4
injection
reflection
mirroring
reference quadratic convergence
200
100
400
800
1600
3200
1/h
∆
s
10
-8
10
-7
10
-6
10
-5
10
-4
10
-3
injection
reflection
mirroring
reference quadratic convergence
200
100
400
800
1600
3200
1/h
∆
s
W
10
-8
10
-7
10
-6
10
-5
10
-4
10
-3
injection
reflection
mirroring
reference linear convergence
reference quadratic convergence
200
100
400
800
1600
3200
1/h
∆
s
W
10
-7
10
-6
10
-5
10
-4
10
-3
10
-2
injection
reflection
mirroring
reference linear convergence
reference quadratic convergence
200
100
400
800
1600
3200
FIG. 4.
Convergence study for prescribed rigid motion. Left column: wall moving with
constant speed (case A). Right column: wall moving with constant acceleration (case B). Solid
reference line indicates linear convergence; dashed reference line indicates quadratic convergence.
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
27
t
x
−
C
−
C
Undisturbed
u = 0
c = c
0
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
8888888888
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
x(t)
.
x(t)
C
+
FIG. 5.
Wave diagram for continuous piston withdrawal.
none of the extrapolation schemes considered here can be expected to minimize all
error indicators for all possible tests. The role of interface acceleration is clearly
important for simulations of shock interactions at an EL-interface, and it is further
investigated in the following Sections.
5.2. Free Expansion
The free-expansion experiment is the simplest nontrivial test of GEL coupling. It
is simple because an exact solution is available; it is nontrivial because the physics
of the problem has the ∞uid and the solid tightly coupled together.
The setup of the free-expansion problem consists of a frictionless piston in a tube
with a vacuum to the right of the piston and an initially constant state to the left
at pressure, density, and sound speed given by
P
o
,
Ω
o
,
c
=
p
∞P
0
=Ω
0
,
∞
being the
ratio of speciØc heats (see Fig. 5).
The solution to the problem can be simpliØed with the deØnition of the following
time constant
ø
=
2
mc
o
P
o
(1 +
∞
)
(19)
DRAFT
September 9, 2002, 11:39am
DRAFT
28
ARIENTI, HUNG, MORANO AND SHEPHERD
where
m
is the mass per unit area of the frictionless piston. The speed of the piston
and the pressure at the interface are determined by the method of characteristics
[17]. Writing
t
§
=
t=ø
, the solution is given by
u
c
o
=
2
∞
°
1
"
1
°
∑
1
1+
t
§
∏
∞
°
1
∞
+1
#
(20)
and
P
P
o
=
∑
1
1+
t
§
∏
2
∞
∞
+1
:
(21)
The exact velocity and pressure are plotted in Fig. 6.
The simplicity of this test case allows us to explore diÆerent temporal couplings
and coupling strategies. Four cases are presented in this Section.
1. Concurrent integration, normal velocity treatment by injection.
2. Concurrent integration, normal velocity treatment by re∞ection.
3. PC-Heun integration, normal velocity treatment by injection.
4. PC-Heun integration, normal velocity treatment by re∞ection.
In all cases, a third-order ENO solver is used for the solution of the ∞uid problem.
For cases using concurrent integration (cases 1 and 2), a third-order TVD time
integration is used for the ∞uid and explicit integration is used for the piston motion.
For cases using PC-Heun integration (cases 3 and 4), a second-order predictor-
corrector method is used for the ∞uid. The motion of the solid is also written as a
Ørst-order system and integrated with a second-order predictor-corrector method.
The grid size is compared against a characteristic length
L
,
L
=
c
o
ø
=
2
c
2
o
m
P
o
(1 +
∞
)
:
(22)
Values ¢
x=L
of 0
:
2, 0
:
1, 0
:
05, and 0
:
025 are used for the convergence study, cor-
responding to 20, 40, 80, and 160 grid cells. The total time of the simulation is
t
§
= 2. A CFL number of 0.1 is used for all simulations.
DRAFT
September 9, 2002, 11:39am
DRAFT
A LEVEL SET APPROACH TO EULERIAN-LAGRANGIAN COUPLING
29
0
0.5
1
1.5
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Nondimensional time t
*
= t/
τ
Nondimensional value
Nondimensional Velocity
Nondimensional Pressure
FIG. 6.
Nondimensional velocity of piston of the piston-air system and the nondimensional
pressure at the piston-air interface.
DRAFT
September 9, 2002, 11:39am
DRAFT