An Eulerian Lagrangian Scheme For Solving Large Deformation Fluid Structure Interaction Problems
An Eulerian Lagrangian Scheme For Solving Large Deformation Fluid Structure Interaction Problems
An Eulerian Lagrangian Scheme For Solving Large Deformation Fluid Structure Interaction Problems
Abstract
1 Introduction
The work presented here describes a theoretical and numerical approach for
“full physics” dynamic simulations of fluid structure interactions involving
large deformations and material transformations (like phase change). By “full
physics” we refer to problems involving strong interactions between the fluid
field and solid field temperatures and velocities, with a full Navier Stokes
There are two classes of theoretical models in common use for full physics
FSI problems: one considers separated materials, and the other considers
their behavior on average. The former is separated in the sense that each
material is permitted to occupy some part of space, to the total exclusion
of any other material. Interactions among materials can take place only at
boundaries between materials, across which there is no flow. In the averaged
model, any material can occupy a point in space with some finite probabilty;
hence interactions among materials can take place anywhere. For brevity, we
shall call the first model the “separated” model, and the second will be called
the “averaged” model.
The relative merits of the two models is frequently judged by the robustness
of their very different numerical solutions, rather than their accuracy or
their ability to address a wide range of physical scenarios. In the separated
model, accurate solutions for the structure can be obtained using commercially
available simulation tools based on finite element methods. A simple code
coupling routine performs the boundary mating with Arbitrary Lagrangian
Eulerian (ALE) codes for the fluid part of the domain. Examples of such
approaches are widespread, (see [?,?,?,?,?]). Among the handicaps of this
type of approach is that a large distortion of the solid domain can lead
to catastrophic mesh tangling. Physical transformations, e.g., phase change,
2
between separated materials is generally beyond the capabilities of this type
of approach.
The approach described in this paper uses the averaged model, and addresses
the issue of stress transport by integrating the state of the solid field in a frame
of reference that follows the center of mass motion of the solid field. In this
“material” frame, the transport of the stress is error free. Furthermore, the
Material Point Method permits the solid field to undergo arbitrary distortion.
Because the fluid state is integrated on the Eulerian frame it too can undergo
arbitrary distortion. If the rate of momentum exchange is infinitely fast, the
single velocity limit is obtained, and the interface region between material
fields can be limited to at most a few grid cells. Hence in the differential
limit, the separated model is recovered. This means that with sufficient grid
resolution, the accuracy of the separated model and the robustness of the
averaged model can be obtained simultaneously.
3
2 The Multifield Model
The model equations used in this work are given here. Section 2.1 develops
the model equations for the rate of change of mass, linear momentum, and
total energy of the interacting materials considered. Section 2.2 shows how
the differing frames of reference are reconciled in the FSI approach used here.
Section 2.3 discusses one peculiarity of the approach that arises from the
practice of using different frames of reference for the fluid and solid materials.
The genesis of the Multifield model is given briefly here, a more thorough
version can be found in [?]. The Multifield model is just an ensemble average
of the exact conservation equations for a collection of materials, any one of
which can reside at a point in space time (xo , t). To be general, suppose there
are N materials. The state at (xo , t) is given by the density, velocity, total
energy, stress, and material type indicator [ρo , uo , Eo , σo , αr ], r = 1, 2, ..., N .
The subscript “o” refers to the value of the quantity at a point in space. For
each of the N materials, αr is Saffman’s H function[?], defined everywhere in
space: it is equal to one if r-material is found at (xo , t), and is zero otherwise.
The H function is just a ledger containing the observer’s record of what
material happened to be at (xo , t) in any one of the infinite members of the
ensemble. The exact evolution equations in conservative differential form with
gravitational acceleration g are
∂ρo
= −∇ · ρo uo (1)
∂t
∂ρo uo
= −∇ · ρo uo uo + ∇ · σo + ρo g (2)
∂t
∂ρo Eo
= −∇ · ρo Eouo + ∇ · σo · uo − ∇ · qo + ρo uo · g (3)
∂t
These equations are closed by appropriate constitutive models for the stress
(or stress rate), equations of state, and the heat flux vector qo due to thermal
conduction. An exact equation for the evolution of αr is not available, except
when the materials never change identity (that is, no transformations in
material type occur, in which case αr is invariant). For the purposes here
we just need to be clear that
∂αr Do α r
+ uo · ∇αr ≡ , (4)
∂t Dt
4
where we use the subscript o on the material derivative to signify that the
material trajectory is uo . The equivalence sign is used to emphasize that this
is just a definition, not a closure relation. Equation 4 is used as a placeholder
for what represents the change in material type due to physical processes
(mass exchange between fields).
Let angle brackets signify the ensemble average, which is like an arithmetic
average over the ensemble; this is formally an integral over a weighted
probability function. The real quantities of interest are the ensemble averaged
ones, defined with the help of Saffman’s H function. The data at (xo , t) are
discontinuous (called weak solutions); application of the ensemble averaging
operator creates continuous ones at general locations in space time (x, t). The
continuous (averaged) variables are refered to as “field” variables. They are
Evolution equations for the field variables are formally obtained by taking the
partial time derivative of the definition, and rearranging with the help of the
exact equations. For the mass conservation equation, the result is
∂ρr Do α r
+ ∇ · ρ r ur = ρ o
∂t Dt
!
∂ρr Do α r
Z Z
+ ∇ · ρr ur dV = ρo dV . (5)
∂t Dt
V V
The left side is the total rate at which r-mass changes in a volume moving
with the r-field mean velocity; it is equal to the r-field material derivative of
the r-mass Dr Mr /Dt (where the subscript r on the material derivative means
that the transport velocity is ur ). Accordingly, the right side is the total rate
at which r-mass is created or destroyed in V due to conversion to/from other
material types, while moving along the instantaneous trajectory uo . Because
any field can contribute/receive mass to/from the r-field, we can express this
by a sum over all materials. Let Γrs be the rate at which r-mass is converted
to/from s-mass, averaged in V . Then we can write
5
Dr M r Do α r
Z
PN
= ρo dV = V s=1 Γrs = V Γr (6)
Dt Dt
V
which signifies the net rate at which r-material is generated (or depleted) in
an arbitrary volume, as a result of conversion to/from all other fields. This is
the total mass exchange rate, per unit of volume. These steps define a volume
average of the ensemble average of the differential equation for the mass.
Using the foregoing notation, we have the (volume & ensemble) averaged
equations
1 Dr M r
= Γr (7)
V Dt
1 Dr (Mr ur )
= − ∇ · hρo αr u0r u0r i + θr ∇ · σ + ∇ · θr (σr − σ)
V Dt PN PN
+ ρr g + s=1 frs + s=1 urs Γrs (8)
1 Dr (Mr er ) Dr v r
= − ρr p − ∇ · hρo αr e0r u0r i + θr τ r : ∇ur − ∇ · θr qr
V Dt DtP
N PN
+ ρr εr + s=1 hsr + s=1 (e + pv)rs Γrs (9)
where the right side of each equation is the averaged rate occuring in V . Also,
the averaged r-field internal energy comes from Er = er + 21 u2r ; in which er
typically includes all internal modes (translational + vibrational + rotational
+ chemical).
Equations 7 - 9 are the averaged model equations for mass, momentum, and
internal energy of r-material. The nonsubscripted σ is the mean mixture stress,
taken here to be isotropic, so that σ = −pI in terms of the hydrodynamic
pressure p. In this, the deviatoric part of the material stress is given by
τ r = σr − 13 IT r(σr ).
In Eq. 8 the quantity − hρo αr u0r u0r i is the r-material fluctuational momentum
flux, otherwise known as the multifield Reynolds Stress. Similarly, in Eq. 9
− hρo αr e0r u0r i is the r-material fluctuational flux of internal energy, while ρr εr
is the r-material fluctuational dissipation rate. Each of these multiphase
turbulence terms, which must be modeled, was neglected in the current work,
they are included above for completeness.
In Eq. 8 the term N s=1 frs signifies a model for the momentum exchange among
P
materials. This term results from the deviation of the r-field stress from the
mean stress, averaged. This is typically modeled as a function of the relative
velocity between materials at a point. (For a two material problem this term
might look like f12 = K12 θ1 θ2 (u1 −u2 ) where the coefficient K12 determines the
6
rate at which momentum is transferred between materials). Likewise, in Eq. 9,
PN
s=1 hsr represents an exchange of heat energy among materials and for a two
material problem typically takes the form h12 = H12 θ1 θ2 (T2 − T1 ) where Tr is
the r-material temperature and the coefficient Hrs is analogous to a convective
heat transfer rate coefficient. The heat flux is given by qr = −ρr br ∇Tr where
the kinematic thermal diffusion coefficient br is an effective one that contains
both molecular and turbulent effects (when the turbulence is included).
In Eqs. 8 and 9, the velocity urs and the enthalpy (e + pv)rs are the velocity
and enthalpy carried with the material that is converted into that of the r-
field, from material that was formerly accounted for in the s-field. For practical
purposes these are simply the mean values associated with the field donating
its mass to r-field. The rate at which the conversion occurs is given by a model
for Γrs which can by any of the well known forms that satisfy the law of mass
action. For a specific example, consider a solid propellant (field-1) converted by
burning into a reaction product gas (field-2). Typically Γ12 = −C(T1 , p)M1 /V
where the phenomenological coefficient depends on the propellant temperature
and the system pressure. Generally, once burning commences it proceeds until
the mass is depleted.
Given appropriate constitutive models for the mass exchange rate, the force
density, heat flux, and heat exchange, Eqs. 7-9 are still unclosed (despite
ignoring the turbulence effects.) The temperature Tr , specific volume vr ,
volume fraction θr , and hydrodynamic pressure p are all yet to be specified;
these quantities are related to the mass density ρr and internal energy er by
way of equations of state. The four relations for the four quantites (Tr , vr , θr , p)
are
er = er (vr , Tr ) (10)
vr = vr (p, Tr ) (11)
θr = ρ r vr (12)
PN
0 = 1 − s=1 ρs vs (13)
Equations 10 and 11 are, respectively, the caloric equation of state and the
thermal equation of state. Equation 12 defines the ensemble averaged H
function as the volume fraction, θ and with that definition, Eq. 13 follows. We
shall refer to Eq. 13 as the multifield equation of state. The hydrodynamic
pressure p is the one value permitting arbitrary masses of the multiple
materials to occupy the entirety of the volume V (for this reason it has been
called the “equilibration” pressure).
For the present purposes, the final closure relation needed is for the material
stress σr . In symbolic form this is
7
σr = Φ[F, σr , ∇ur , ∗r ] (14)
The operator Φ is the general material response law that may depend on the
deformation gradient, the stress itself, the averaged r-material rate of strain
and a list of variables indicated by the ∗r which may include history variables
such as a damage parameter, as well as the rotation of the polar axes relative
to the axes in a laboratory frame of reference.
Equations 7-14 form a set of eight equations for the state vector [Mr , ur , er , Tr , vr , θr , σr , p],
in any arbitrary volume of space V . The Eulerian Lagrangian approach for
large deformation FSI problems requires obtaining approximate solutions to
this equation system. To clarify some subtle aspects of the numerical solution
method used here, a discussion of the different frames of reference employed
is given next.
The averaged field equations displayed in the previous section specify the
Lagrangian rate of change for the state of the materials momentarily residing
in an arbitrary volume of space. The coordinates of that volume may change
with time. The key feature of the FSI approach in this paper is the choice of
reference frame on a material by material (field by field) basis. To affect this
flexibility a generalized Reynolds Transport Theorem is used [?]. Because of
its importance to this work, the theorem is developed here.
Z " #
D ∂q
Z Z
[q]dV = dV + [n̂ · qu]dS , (15)
Dt ∂t
V (t) V (t) S(t)
where q is a quantity which varies within the volume. The integral of q over
the volume is the total of that quantity Q. (For q = 1 this is just a kinematic
expression for the rate of volume change, along the trajectory u.) In the
multifield case there are N different trajectories associated with the averaged
velocity of the N fields in the problem. Let the subscript r designate one of
those fields. Leibnitz’s Rule for any quantity associated with that field is
" #
Dr Dr Q r ∂qr
Z Z Z
[qr ]dV = = dV + [n̂ · qr ur ]dS . (16)
Dt Dt ∂t
Vr (t) Vr (t) Sr (t)
8
Now consider a second volume (that of a computational grid, say), moving
with a completely arbitrary velocity ug which may be time dependent, time
steady, or zero. The corresponding rule, operating again on qr is
Z " #
Dg Z Dg Q r ∂qr Z
[qr ]dV = = dV + [n̂ · qr ug ]dS . (17)
Dt Dt ∂t
Vg (t) Vg (t) Sg (t)
Let us now consider an instant in time at which the surfaces bounding the two
volumes is exactly overlapping. At that instant the integrals cover precisely
the same space, so we can subtract one from the other thereby eliminating the
term involving the partial derivative. That is,
Dg Q r Dr Q r
Z
+ [n̂ · qr (ur − ug )]dS = . (18)
Dt Dt
Sg (t)
Dg Q r Dr Q r
+ V ∇ · qr (ur − ug ) = , (19)
Dt Dt
in which we implicitly assume that the gradient is constant in the volume V .
Equation 19 is the generalized Reynolds Transport Theorem. It relates the
change in a total quantity (Q) in a volume moving with velocity ug to the
rate of change along the material motion ur . The difference between the two
rates is the second term on the left, called the advection term. Equation 19 is
valid, so long as the r-Lagrangian volume and the grid volume are precisely
overlapping at the instant for which the evaluation is made.
In the FSI approach used here, ug = ur for the field representing a solid
structure and for the fluid field ug = 0. Hence we say that the solid field state
is integrated in the Lagrangian (material) frame of reference, and the fluid field
state is integrated in the Eulerian (laboratory) frame. Thus, the approach uses
mixed frames of reference. At the instant the states are to be updated in time,
the solid field state is interpolated to the center of each grid volume in the
computational domain. In this way, the gradients and state differences that
measure the change in state are evaluated in the same volume of space, as is
required by Eq. 19.
When mixed reference frames are used it becomes necessary to integrate the
state relations Eqs. 10-13 in the differential form. The reason for this was
9
shown in detail by Kashiwa et al.[?] by considering the model equations in the
incompressible limit. This section illustrates the issue, and its resolution.
Consider the multiphase equation of state for the two field case, and multiply
by the grid volume V in which the equation is to be satisfied. Because
ρr = Mr /V this is
0 = V − Ms vs (p, Ts ) − Mf vf (p, Tf )
where the subscripts designate solid and fluid respectively. Now let vs =
constant, which corresponds to an isothermal incompressible material. Let
the fluid field be a perfect gas with gas constant Rf , so that vf = Rf Tf /p. The
solution for p is
p = Mf Rf Tf /(V − Ms vs )
which can always be found, providing that the solid field does not occupy all
of the space in V , that is, Ms vs < V . If the field equations for both solid and
fluid are integrated in the Eulerian frame, this condition can be guaranteed
by a physical stress on the solid field that is compressive (sometimes called a
“configuration” stress because the specific volume of the solid is constant). If
instead the solid field is integrated in the Lagrangian frame, the mass Ms in
V is approximated by interpolation from mass points near V , in which case
the condition Ms vs < V can not be guaranteed in general.
The apparent pathology is removed by using Eqs. 11-13 and the generalized
Reynolds Transport Theorem for the r-material specific volume to derive an
expression for the specific volume, or volume per unit of mass. The specific
volume can be considered to be a dynamic (averaged) variable of the state, to
be integrated in time from an initial condition. When this is done, the total
volume associated with the materials is
PN
Vt = r=1 Mr v r
The development of the rate equation of the r-material specific volume is given
in the appendix. The final result is:
N
" #
Dr (Mr vr ) Dr T r Ds T s
= V vr Γr + frθ ∇ · u + θr βr − frθ
X
θ s βs . (20)
Dt Dt s=1 Dt
10
is discussed in the next section, where the numerical solution algorithm
is described. In the limit κ → 0 all materials are incompressible, and
Eq. 51 becomes the equation for p whose numerical approximation requires
no correction. Hence in the incompressible limit the discrete approximation
is once again exact; for compressible materials, we exhibit the goodness of
the approximation by way of examples comparing exact and approximate
solutions.
3 Solution Algorithm
Our use of the cell centered ICE method employs time splitting: first, a
Lagrangian step updates the state due to the physics of the conservation laws
(i.e., right hand side of Eqs. 7-9); this is followed by an Eulerian step, in which
the change due to advection is evaluated. For solution in the Eulerian frame,
the method is well developed and described in [?].
11
5
v gas phase
condensed phase
condensed phase (altered)
4
0
0 100000 200000
P
Fig. 1. Specific volume vs pressure for a gas phase material and a condensed phase
material. Light dashed line reflects an altered condensed phase equation of state to
keep all materials in positive equilibration pressure space.
The difficulty, and the modification that resolves it, can be understood
by considering a solid material in tension coexisting with a gas. For solid
materials, the equation of state is just the bulk part of the constitutive
response (that is, the isotropic part of the stress versus specific volume and
temperature). If one attempts to equate the isotropic part of the stress with
the fluid pressure, there exist regions of pressure-volume space for which Eq. 13
has no physical solutions (because the gas pressure is only positive). This can
be seen schematically in Fig. 1, which sketches equations of state for a gas
phase (solid line) and a condensed phase (heavy dashed line), at an arbitrary
temperature. Recall that the isothermal compressiblity is the negative slope
of the specific volume versus pressure. Imbedded structures considered here
are condensed phase and possess a very much smaller compressibility than
the gasses in which they are submerged, at low pressure. Nevertheless the
variation of condensed phase specific volume can be important at very high
12
pressures, where the compressibilities of the gas and condensed phase materials
can become comparable (as in a detonation wave, for example). Because the
speed of shock waves in materials is determined by their equations of state,
obtaining accurate high pressure behavior is an important goal of our FSI
studies.
A new feature of the approach presented here is the ability to neatly couple a
Lagrangian description for selected materials within the overarching multifield
method. Specifically, the underlying theory is the multifield equations, solved
on an Eulerian mesh. However, materials with solid like behavior and
history dependent response laws are more conveniently treated by Lagrangian
methods. The unification of the Lagrangian method within the multifield
solution is described in Sec. 3.3.
The Lagrangian frame representation adopted here describes the state vari-
ables of the material on particles, or “material points.” The specific numerical
technique used is known as the Material Point Method (MPM). Originally
described by Sulsky, et al., [?,?], the MPM is a particle method for struc-
tural mechanics simulations. MPM is an extension to solid mechanics of FLIP
[?], which is a particle in cell (PIC) method for fluid flow simulation. The
method typically uses a regular structured grid as a computational scratchpad
for computing spatial gradients. This same grid also functions as an updated
Lagrangian grid that moves with the particles during advection and thus elim-
inates the diffusion problems associated with advection on an Eulerian grid.
At the end of a timestep, the grid is reset to the original regularly ordered
position.
13
In explicit MPM, the equations of motion are cast in the form [?]:
The solution procedure begins by interpolating the particle state to the grid,
to form Mg , Fextg , and to determine a velocity on the grid vg . In practice,
a lumped mass matrix is usually used. These quantities are calculated by the
following equations, where the represents a summation over all particles:
P
p
Sip mp vp
P
X p X
Mi = Sip mp , vi = , Fexti = Sip Fextp (22)
p Mi p
At this point, a velocity gradient, ∇vp is computed at the particle using the
velocities interpolated to the grid:
X
∇vp = Gip vi (23)
p
where Gip is the gradient of the ith nodes shape function, evaluated at xp .
X
Finti = Gip σp Vp , (24)
p
14
Equation 21 can then be solved for ag . An explicit backward Euler method is
used for the time integration:
vLi = vi + ai dt (25)
and the particle position and velocity are explicitly updated by:
X
vp (t + dt) = vp (t) + Sip ai dt (26)
i
Sip vLi dt
X
xp (t + dt) = xp (t) + (27)
i
3.3 Integration of the Material Point Method within the Eulerian multifield
ICE Formulation
The key feature of this work is the development, demonstration, and validation
of a tightly coupled, multifield “averaged” solution approach that allows
simulations of an arbitrary number of materials undergoing a variety of
physical and chemical processes. Fundamental to the proposed approach is the
representation of a material field in either an Eulerian or Lagrangian reference
frame, integrated within a general multifield description. This allows treating
specific phases in their traditionally preferred frame of reference, Lagrangian
for solid, Eulerian for fluid. The MPM algorithm, described in the previous
section, is used to time advance those materials that are best described in a
Lagrangian reference frame. However, by choosing the background mesh used
to update the MPM materials, to be the same mesh used in the multi-material
Eulerian description, all interactions among materials can be enforced to occur
in the common framework. This results in a robust and tightly coupled solution
for interacting materials with very different response behaviors.
All materials, solid or fluid, will have a presence in the Eulerian multifield
description. This common reference frame is used for all physics that involve
15
mass, momentum, or energy exchange among materials. This allows for a
tight coupling between the fluid and solid phases. The coupling occurs directly
within the terms in the multifield equations. Since a common Eulerian frame
of reference is used for interactions among materials, typical problems with
convergence and stability of solutions for separate domains communicating
only through boundary conditions are alleviated. Thus, all solid phase
materials thus have a dual representation - in both the Lagrangian (MPM)
and Eulerian framework.
16
(2) Compute the equilibrium pressure: While Eq. 13 and the
surrounding discussion describes the basic process, a few specifics warrant
further explanation. In particular, the manner in which each material’s
volume fraction is computed is crucial. Typically, in multi-material CFD
calculations, this would simply be the material volume divided by the
cell volume. However, because the solid and fluid materials are advected
in different manners (see below) the total volume of material in a cell
is not necessarily equal to the volume of a computational cell. Because
of this, it is important to have an accurate accounting of the volume of
material in each cell. This is done by solving the evolution equation for
each material’s specific volume given by Eq. 53 and described below in
step 11.
With the materials’ masses and specific volumes in hand, material
volume can be computed and summed to find the total material volume.
The volume fraction θr is then computed as the volume of r-material per
total material volume. With this, the solution of Eq. 13 can be carried out
at each cell using a Newton-Raphson technique[?], which results in new
values for the equilibrium pressure, peq , volume fraction, θr and specific
volume, vr .
(3) Compute face centered velocities, u∗r , for the Eulerian advec-
tion: At this point, fluxing velocities are computed at each cell face.
The expression for this is based on a time advanced estimate for the cell
centered velocity. A full development can be found in [?] and [?] but here,
we only state the result:
!
ρ r ur + ρ r R ur R 2vrL vrR ∆t peqR − peqL
u∗r = L L − (28)
ρ rL + ρ rR v rL + v rR ∆x
The first term above is a mass weighted average of the logically left and
right cell centered velocities and the second term is a pressure gradient
acceleration term. Not shown explicitly here is the necessary momentum
exchange at the face centers. This is done here in the same manner as
it is described subsequently in step 10 for the cell centered momentum
exchange.
(4) Multiphase chemistry: Compute sources of mass, momentum, energy
and specific volume as a result of phase changing chemical reactions
for each r-material, Γr , ur Γr , er Γr and vr Γr . Specifics of this step are
model dependent, but in general the most important consideration
is to accurately account for the changes in state for all materials
involved, particularly the reactant material. Failure to do so can result
in unreasonable values for temperature and velocity as the mass of the
reactant material is consumed. When the reactant material is described
by particles, this is somewhat simplified as decrementing the particle
17
mass automatically decreases the momentum and internal energy of that
particle by the appropriate amount. Otherwise, care must be taken to
reduce the momentum and internal energy of the reactant by amounts
proportional to the mass consumed each timestep. This mass, momentum
and internal energy is transferred to the product material’s state.
(5) Compute an estimate of the time advanced pressure, p: Based
on the volume of material being added to (or subtracted from) a cell in
a given timestep, an increment to the cell centered pressure is computed
using:
N N
v r Γr − ∇ · (θu∗ )r
P P
r=1 r=1
∆p = ∆t N
(29)
(θκ)r
P
r=1
p = peq + ∆p (30)
where the subscripts L and R refer to the logically left and right cell
centered values respectively, and ρ is the sum of all material’s densities
in that cell. This will be used subsequently for the computation of the
pressure gradient, ∇p∗ .
(7) Material Stresses: For the solid, we calculate the velocity gradient at
each particle based on the grid velocity (Eq. 23) for use in a constitutive
18
model to compute particle stress. Fluid stresses are computed on cell
faces based on cell centered velocities.
(8) Accumulate sources of mass, momentum and energy at cell
centers: These terms are of the form:
N
X
∆(m)r = ∆tV Γs (32)
s=1,s6=r
N
∆(mu)r = −∆tV θr ∇p∗ + ∇·θr (σr − σ) +
X
us Γ s (33)
s=1,s6=r
N N
∆(me)r = −∆tV frθ p ∇ · (θu∗ )s +
X X
e s Γs (34)
s=1 s=1,s6=r
(10) Momentum and heat exchange: The last step in the Lagrangian
phase is the exchange of momentum and heat between materials.
N
(mu)Lr = (mu)L− θr θs Krs (uLs − uLr )
X
r + ∆tmr (38)
s=1
N
(me)Lr = (me)L− θr θs Rrs TsL − TrL
X
r + ∆tmr cvr (39)
s=1
19
results in driving contacting materials to the same velocity. Intermaterial
heat exchange is usually modeled at a lower rate. Again, note that the
same operation must be done following Step 3 above when computing the
face centered velocities.
(11) Specific volume evolution: As discussed above in step 2, in order
to correctly compute the equilibrium pressure and volume fraction, it is
necessary to keep an accurate accounting of the specific volume for each
material. Here, we compute the evolution in specific volume due to the
changes in temperature and pressure, as well as phase change, during the
foregoing Lagrangian portion of the calculation, according to:
N N
" #
frθ ∇ θs∗ u∗s frθ
X X
∆(mv)r = ∆tV vr Γr + · + θr βr Ṫr − θs βs Ṫs (40)
s=1 s=1
(mv)Lr = (mv)nr + ∆(mv)r (41)
−T L t
where β is the constant pressure thermal expansivity and Ṫ = T ∆t is
the rate of change of each material’s temperature during the Lagrangian
phase of the computation.
(12) Advect Fluids: For the fluid phase, use a suitable advection scheme,
such as that described in [?], to move fluid material between cells. This
includes the transport of mass, momentum, internal energy and specific
volume. As this last item is an intensive quantity, it is converted to
material volume for advection, and then reconstituted as specific volume
for use in the subsequent timestep’s equilibrium pressure calculation.
(13) Advect Solids: For the solid phase, interpolate the time advanced grid
velocity and the corresponding velocity increment back to the particles,
and use these to advance the particle’s position and velocity, respectively.
This constitutes advection of the solid phase material.
This completes one timestep. In the preceding, the user has a number of
options in the implementation. The approach taken here was to develop a
working MPM code and a separate working multifield ICE code. We note,
however, that the fluid structure interaction methodology should not be looked
at in the context of a “marriage” between an Eulerian CFD code and the
MPM. The underlying theory is a multifield description that has the flexibility
to incorporate different numerical descriptions for solid and fluid fields within
the overarching solution process. To have flexibility in treating a widest range
of problems, it was our desire that in the coupling of the two algorithms, each
of the original codes be able to function independently, or in concert.
Thus, for instance, in step 8 above, for the solid material we’ve chosen not to
aggregate all of the momentum sources at the cell centers. Rather, the pressure
gradient and mean pressure are computed at the cell centers, and interpolated
20
back to the nodes, where the integration of the MPM equations normally take
place. The time advanced solid material momentum is then aggregated at the
cell centers for momentum exchange, and the resulting change in momentum
is propagated back to the nodes and eventually to the particles. Also, as
described here, this method is fully explicit in time. To make this implicit
with respect to the propagation of pressure waves, a Poisson equation is solved
in the calculation of ∆p, which is in turn used to iteratively update the face
centered velocities. This is described in [?].
4 Validation Tests
21
Air P=101.3 kPa
Solid
P r2 r2 P r2 r2
! !
σradial = 2 i 2 1 − o2 σθ = 2 i 2 1 + o2 (42)
(ro − ri ) r (ro − ri ) r
1 X ||σr + σθ | − 23 P |
Error = 2 (43)
N N 3
P
Results from that are shown in Fig. 3b. When plotted in this log-linear manner,
the slope of the straight line that best approximates these points indicates the
order of accuracy. Here, a slope of 0.89 is found.
22
-2
σr
7
σθ
+0
-2.2
0E
σr
3.
σθ -2.4
7
+0
-2.6
0E
2.
-2.8
log(Error)
07
E+
-3
0
1.
Stress (Pa)
-3.2
0
+0
0E
-3.4
0.
-3.6
07
E+
-3.8
.0
-1
-4
0.5 0.6 0.7 0.8 0.9 -5 -4.5 -4 -3.5 -3
07
log(dx)
E+
Radius (m)
.0
-2
a b
Fig. 3. Results from pressurization of an annulus, (a) radial and circumferential
stresses, compared to the theoretical solutions, of the row of particles nearest the
x-axis, and (b) natural log of solution error, averaged over all particles, vs. natural
log of grid spacing. The slope of the line fit to the four data points indicates an
order of accuracy of 0.89.
γ
V1
P2 = P 1 (44)
V2
A large momentum exchange coefficient, (Krs = 1015 ) was chosen. This drives
the velocity of the gas and the piston to the same value in cells where
both gas and piston are present (i.e., at the interface between the gas and
the solid). This enforces a no-slip, no-interpenetration condition between the
piston and the gas. For this simulation, the compression ratio of the gas (initial
volume/final volume) reached a maximum of 5.9.
23
14
Simulation Results
Exact solution
12
10
8
Pressure
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Time
The final validation test presented here is a symmetric 20o wedge moving
at Mach 3 through initially stationary air. The main purpose of this test is
to demonstrate the ability of the compressible multifield ICE formulation to
capture the oblique shocks and expansion fan that result from the interaction
of the solid wedge with the surrounding air. Traditionally, tests such as this are
done by keeping the geometry fixed and describing it as a boundary condition
to a CFD code which moves the fluid over the stationary item.
24
β=38ο
M=3
θ=20ο
as it interacts with the flow field. We believe that these are capabilities unique
to this approach.
25
a b
c d
26
reference, thus alleviating some of the problems typically encountered when
attempting to treat both phases in a single reference frame. Several simple
validation tests were performed, each of which was intended to bring a degree
of confidence to one or more elements of the physics involved in multifield
problems. Finally, an example simulation was presented to demonstrate the
complex types of problems that this method is capable of treating.
Dr (Mr vr ) Dr M r Dr v r
= vr + Mr
Dt Dt Dt
Dr v r
= v r V Γr + M r (45)
Dt
The Lagrangian rate for the change of vr is obtained by differentiation of the
state equation, Eq. 11. That is
Dr v r Dr p Dr T r
= vr −κr + βr (46)
Dt Dt Dt
where the partial derivatives are
" #
1 ∂vr
κr = − isothermal compressibility ; (47)
vr ∂p T
27
" #
1 ∂vr
βr = constant pressure volumetric expansivity . (48)
vr ∂T p
In the problems shown here, ideal materials with constant specific heats are
assumed.
The rate equation for the hydrodynamic pressure begins with the differenti-
ation of the multifield equation of state, Eq. 13, and proceeds with the use
of conservation of mass, the definition of the material derivative of vs , the
product rule and Eqs. 46 and 12.
N N
" # " #
∂ X X ∂vs ∂ρs
1− (ρs vs ) = 0 = ρs + vs
∂t s=1 s=1 ∂t ∂t
N
" #
X ∂vs
= ρs − vs ∇ · (ρs us ) + vs Γs
s=1 ∂t
N
Ds v s
X
= ρs − ρs us · ∇vs − vs ∇ · (ρs us ) + vs Γs
s=1 Dt
N
Ds v s
X
= ρs − ∇ · (ρs us vs ) + vs Γs
s=1 Dt
N N N
Ds p Ds T s
X X X
= θs −κs + βs −∇· θ s us + v s Γs
s=1 Dt Dt s=1 s=1
PN
κ= s=1 θ s κs , (49)
PN
u= s=1 θ s us . (50)
N
Dp p Ds T s
X
κ = −∇ · u + θ s βs (51)
Dt s=1 Dt
PN
up = ( s=1 θs κs us ) /κ .
28
N
" # " #
Dr p Dp p ∂p ∂p X
κ −κ = κ + κur · ∇p − κ + θs κs us · ∇p
Dt Dt ∂t ∂t s=1
N
X N
X
= ur θs κs · ∇p − θs κs us · ∇p
s=1 s=1
= κ (ur − up ) · ∇p
or:
Dr p Dp p
= + (ur − up ) · ∇p (52)
Dt Dt
Substitution of Eq. 52 into Eq. 46 gives a rate equation for the specific volume
of the r-material:
N
!
M r Dr v r Dr T r Ds T s
= frθ ∇ · u + θr βr − frθ
X
θ s βs
V Dt Dt s=1 Dt
N
− frθ
X
vs Γs − κr θr (ur − up ) · ∇p (53)
s=1
where frθ = θr κr /κ. Finally, Eq. 53 is used with Eq. 45 to give a rate equation
for the volume of r-material
N
!
1 Dr (Mr vr ) Dr T r Ds T s
= frθ ∇ · u + θr βr − frθ
X
θ s βs
V Dt Dt s=1 Dt
N
!
frθ
X
+ v r Γr − v s Γs
s=1
− θr κr (ur − up ) · ∇p (54)
Acknowledgements
29
References
[4] S. Rugonyi, K. Bathe, On finite element analysis of fluid flows fully coupled
with structural interactions, CMES 2 (2001) 195–212.
[5] K. Bathe, H. Zhang, S. Ji, Finite element analysis of fluid flows coupled with
structural interactions, Computers & Strctures 72 (1999) 1–16.
30
[17] D. Sulsky, S. Zhou, H. Schreyer, Application of a particle-in-cell method to solid
mechanics, Computer Physics Communications 87 (1995) 236–252.
[25] J. Guilkey, J. Weiss, Implicit time integration for the material point method:
Quantitative and algorithmic comparisons with the finite element method,
International Journal for Numerical Methods in Engineering 57 (2003) 1323–
1338.
31