LINEBURG


<< . .

 18
( 23)



. . >>

Benson [14, 15].
While Eulerian numerical methods are superior for high-speed gas ¬‚ows,
they can perform poorly for solid mechanics calculations. The capturing
17 nature of Eulerian methods generally leads to algorithms that are not ac-
curate or robust enough to track time-history variables or material response
to loading and damage. On the other hand, Lagrangian numerical methods
Solid-Fluid Coupling are extremely accurate and well tested in this area.
Many researchers agree that it is preferable to use Eulerian numerical
methods for ¬‚uids and Lagrangian numerical methods for solids. Then there
are two popular approaches for treating the solid-¬‚uid interface. First, one
can smear out the nature of the numerical approximations using a La-
grangian method in the solid and an Eulerian method in the ¬‚uid with a
“mushy” region in between where the grid moves with an intermediate ve-
locity. That is, the grid velocity is smoothly varying between the Lagrangian
mesh velocity and the zero-velocity Eulerian mesh. This is the fundamen-
tal idea behind arbitrary Lagrangian-Eulerian (ALE) numerical algorithms;
see, e.g., [14]. The problem with this approach is that the variable velocity
mesh has not been well studied, and the numerical algorithms employed
on it tend to be low-order accurate and suspect. The second approach for
treating the solid-¬‚uid interface is to keep the mesh representation sharp so
that the Eulerian and Lagrangian meshes are in direct contact. The problem
17.1 Introduction with this approach is that the Lagrangian mesh moves, causing Eulerian
mesh points to appear and disappear. In addition, the Eulerian cells tend
Solid-¬‚uid interaction problems are still rather di¬cult for modern compu-
to have irregular shapes, referred to as cut cells. These cut cells can lead to
tational methods. In general, there are three classical approaches to such
numerical errors and sti¬ time-step restrictions; see, e.g., [14] and the ref-
problems: One can treat both the solid and the ¬‚uid with Eulerian numer-
erences therein, speci¬cally Hancock [78], Noh [121], and McMaster [112],
ical methods, the ¬‚uid with an Eulerian numerical method and the solid
which discuss the PISCES, CEL, and PELE programs, respectively.
with a Lagrangian numerical method, or both the solid and the ¬‚uid with
In [62], the author took the second approach for the treatment of the
Lagrangian numerical methods.
solid-¬‚uid interface. However, problems with cut cells were avoided by us-
Many ¬‚uid ¬‚ows, e.g., high-speed gas ¬‚ows with strong shocks and large
ing ghost cells for the Eulerian mesh. These ghost cells are covered (or
deformations, are di¬cult to simulate with Lagrangian numerical methods
partially covered) by the Lagrangian mesh, but are used in the Eulerian
that use arti¬cial viscosity to smear out shock pro¬les over a number of
¬nite di¬erence scheme in order to circumvent small time-step restrictions.
zones in order to reduce postshock oscillations or ringing. The arti¬cial
The ghost cells are de¬ned in a way consistent with a contact discontinuity
viscosity can be both problem- and material-dependent. Flows with sig-
so that the interface boundary conditions or jump conditions are properly
ni¬cant deformations can cause large mesh perturbations and subsequent
captured. This method also avoids the blending problems associated with
numerical errors that can be removed only with complicated remeshing
covering and uncovering of grid points, since covered real grid nodes are
and/or mesh generation procedures that tend to be low-order accurate.
subsequently treated as ghost nodes, and uncovered ghost nodes are sub-
In particular, vorticity can cause the mesh to tangle and sometimes in-
sequently treated as real grid nodes. Moreover, the numerical treatment of
vert, in which case the calculation needs to be stopped. Eulerian numerical
the solid-¬‚uid interface does not compromise the solution techniques for
methods intrinsically avoid these mesh-associated di¬culties by using a
the solid or the ¬‚uid. That is, once the ghost cells™ values are speci¬ed, a
stationary mesh. Furthermore, these schemes capture shocks in a straight-
standard Eulerian program can be used to advance the ¬‚uid (and its ghost
forward fashion using conservation and robust limiters, eliminating the need
nodes) in time. A standard Lagrangian program can be used to advance the
for problem-dependent arti¬cial viscosity formulations. This allows shocks
solid in time as well as to acquire boundary conditions from the Eulerian
to be modeled with as few as one grid cell without oscillations, whereas La-
mesh using both the real grid nodes and the ghost nodes in a standard
grangian numerical methods usually su¬er from some amount of postshock
17.2. Lagrange Equations 203 204 17. Solid-Fluid Coupling

one-fourth the zonal mass. The nodal mass M p is de¬ned as the sum of the
interpolation procedure. Aivazis et al. [3] successfully used this method
to couple an Eulerian hydroprogram to a highly sophisticated Lagrangian (at most four) neighboring subzonal masses. Once again, the nodal, zonal,
materials program designed by Ortiz. and subzonal masses all remain ¬xed throughout the calculation. The in-
dependent variables are updated with equations (17.1), (17.2) and (17.3)
with x, u, and F n replaced with X, V , and F , respectively. Either force or
17.2 Lagrange Equations velocity boundary conditions are applied to the nodes on the boundary.
For three spatial dimensions we refer the reader to Caramana et al. [25].
While there are a number of sophisticated Lagrangian programs, we par-
ticularly like the approach of Caramana et al. [24] which allows for
implementation of arbitrary forces in a straightforward fashion.
The Lagrange equations are written in nonconservative form with po-
17.3 Treating the Interface
sition, velocity, and internal energy as the independent variables. In one
spatial dimension, both x and u are de¬ned at the grid nodes, while e is
Boundary conditions need to be imposed on both the Eulerian and La-
de¬ned at the cell centers located midway between the nodes. To initialize
grangian grids where the Lagrangian grid partially overlaps the Eulerian
the calculation, the mass of each zone, M z , is determined, and then the
grid. First the interface itself needs to be de¬ned, and since the Lagrangian
subzonal masses mz are de¬ned as half the zonal mass. The nodal mass M p
grid nodes move at the local material velocity, these nodes can be used
is de¬ned as the sum of the neighboring subzonal masses. The nodal, zonal,
to determine the position of the interface. This interface divides the Eule-
and subzonal masses all remain ¬xed throughout the calculation. At each
rian mesh into separate regions, a region populated by real grid nodes and
time step, the location of each grid node is updated according to
a region populated by ghost nodes. Interface boundary conditions for the
xn+1 ’ xn Eulerian mesh are imposed by de¬ning mass, momentum, and energy in
= un , (17.1)
t the ghost nodes. Interface boundary conditions for the Lagrangian mesh are
imposed by either specifying the velocity of the grid nodes on the boundary
where t is the size of the time step. The velocity at each node is updated
or by specifying the force applied to that boundary.
using
Since the interface moves with the local material velocity, it can be
un+1 ’ un Fn treated as a contact discontinuity. The pressure and the normal velocity are
= p, (17.2)
t M continuous across the interface, while the entropy and tangential velocities
where F n is the net force on the grid node. The internal energy in each are uncoupled across the interface. The interface values of the uncoupled
variables are captured by extrapolating these variables across the interface
zone is updated with
into the ghost cells. The continuous or coupled variables are determined
en+1 ’ en Hn
using the values from both the Eulerian and the Lagrangian meshes.
= z, (17.3)
t M The interface normal velocity can be determined by applying any num-
where H n is the heating rate of the zone. One can apply either force or ber of interpolation techniques to the Eulerian and Lagrangian mesh values.
velocity boundary conditions to the grid nodes on the boundary. Velocity However, one should be careful to de¬ne the interface normal velocity in a
boundary conditions are enforced by setting the velocity of a boundary way that is consistent with the material in the Lagrangian mesh. Pertur-
node to the desired boundary velocity instead of solving equation (17.2). bations to the velocity of the Lagrangian grid nodes can provide enormous
Force boundary conditions are applied by adding the boundary force to the stress due to resistive forces such as material strength. For this reason, in
net nodal force F n in equation (17.2). order to determine an accurate (and Lagrangian mesh-consistent) value of
In two spatial dimensions, both X = x, y and V = u, v are de¬ned the normal velocity at the interface, only the Lagrangian mesh is used to
at the grid nodes, which are connected in the same fashion as an Eulerian determine the interface velocity, as was done in [78], [112], and [121]. Both
grid, producing quadrilateral zones. Each quadrilateral zone is split into calculations use this interface normal velocity, so that [VN ] = 0 is enforced.
four subzones by connecting the midpoints of opposite edges of the zone. That is, the Lagrangian mesh uses the computed velocities of its boundary
The internal energy is de¬ned at the zone center located at the intersec- nodes, while the Eulerian calculation captures this interface normal veloc-
tion of the four subzones. To initialize the calculation, the mass of each ity by assigning to each ghost node the interface normal velocity of the
zone, M z , is determined, and then the subzonal masses mz are de¬ned as nearest Lagrangian mesh point on the interface.
17.3. Treating the Interface 205 206 17. Solid-Fluid Coupling

Since the interface normal velocity is de¬ned as the velocity of the nodes Eulerian mesh. This linear interpolation requires valid pressure values in
on the Lagrangian mesh boundary with no contribution from the Eulerian both the real and the ghost nodes. Therefore, the pressure extrapolation
mesh, velocity boundary conditions cannot be enforced on the Lagrangian step needs to be carried out before this linear interpolation step. This Eu-
lerian interface pressure makes a contribution of ±p to the net force on the
mesh at the interface. Instead, force boundary conditions are applied by
interpolating the Eulerian grid pressure to this Lagrangian interface. In this Lagrangian boundary node, depending on whether the Lagrangian mesh
fashion, the interface pressure is determined using only the Eulerian grid lies to the right or to the left of the interface, respectively.
values, ignoring contributions from the Lagrangian mesh, as in [78], [112], With boundary conditions speci¬ed on both the Eulerian and Lagrangian
and [121]. Both calculations use this interface pressure, so that [p] = 0 is meshes, both can be advanced one Euler step in time. Both the Eulerian
enforced. The interface pressure is captured by the Eulerian calculation by real grid nodes and a band of Eulerian ghost nodes are advanced in time.
extrapolating the pressure across the interface into the ghost cells which These ghost nodes are advanced in time so that they have valid values in
is similar to the treatment of entropy and tangential velocity. Then the case they are uncovered by the Lagrangian mesh.
interface pressure is interpolated from the Eulerian grid in order to apply The two-dimensional interface is de¬ned by the line segments of the La-
force boundary conditions to the Lagrangian mesh. grangian mesh boundary that are adjacent to grid nodes of the Eulerian
Noh [121] suggested that it might be better to use some average of the mesh. This interface is used to construct a signed distance function de¬ned
Lagrangian and Eulerian grid values for determining the pressure at the at every Eulerian grid node. Once again, φ is examined on the computa-
interface. For Lagrangian calculations with arti¬cial viscosity and mate- tional boundaries of the Eulerian mesh to ensure that enough ghost nodes
rial strength, the jump condition implies that the net stress in the normal are present, and the size of the Eulerian mesh is increased when necessary.
direction is continuous, not just the pressure. Therefore, this advocated The Eulerian ghost nodes are de¬ned by ¬rst extrapolating S, p, and V
averaging procedure would need to take place between the pressure in the across the interface into the ghost nodes. Then for each ghost node, the
¬‚uid and the normal component of the net stress in the normal direction closest point on the Lagrangian interface is determined. If the closest point
in the solid. However, this can be dangerous, for example, when the La- happens to be on the end of a linear segment, i.e., a Lagrangian grid node,
grangian material is in tension, since near-zero or negative stress might be then that velocity can be designated the closest interface velocity. Other-
calculated at the interface. While Lagrangian methods can be quite robust wise, the closest point is on an edge connecting two Lagrangian grid nodes,
under tension, Eulerian methods can su¬er a number of problems in treat- and the closest interface velocity is determined using linear interpolation
ing near-zero or negative pressures associated with rare¬ed or cavitated between those two nodes. The normal component of the interface velocity
¬‚uids. is combined with the tangential component of the extrapolated velocity
The one-dimensional interface is de¬ned by the location of the La- to determine the velocity at each ghost node. Once the Eulerian ghost
grangian boundary nodes that are adjacent to grid nodes of the Eulerian nodes have valid values for the extrapolated pressure, force boundary con-
mesh. This interface location is used to construct a signed distance function ditions can be determined at the Lagrangian interface. The midpoint of
in order to apply level set methods near the interface. Before de¬ning values each linear interface segment is de¬ned as a control point, and bilinear in-
in the Eulerian ghost nodes, a check is performed to see whether enough terpolation is used to determine the Eulerian mesh pressure at each of these
ghost nodes are present. That is, since the Lagrangian mesh is moving, one control points. Then this pressure is multiplied by both the length and the
needs to ensure that there is adequate overlap between the two meshes. inward-pointing normal of the line segment to determine the magnitude
This is done by examining the values of φ on the computational bound- and direction of the Eulerian pressure force on this segment. Finally, half
aries of the Eulerian mesh. If the computational boundary is an Eulerian of this Eulerian pressure force is added to each of the two nodes that make
ghost node, then the value of φ gives the distance to the interface and up the segment.
can be used to estimate the number of ghost nodes that exist between the Figure 17.1 shows the velocity ¬eld obtained as a shock wave propagates
interface and the computational boundary. Then the size of the Eulerian through an Eulerian gas sandwiched between two Lagrangian materials
mesh can be increased if there are not enough ghost nodes. with strength.
The Eulerian ghost nodes are de¬ned by ¬rst extrapolating S and p.
Then u at each ghost node is assigned the value of u at the nearest La-
grangian boundary node that lies on the interface between the Eulerian
and Lagrangian grids. Force boundary conditions are applied to the La-
grangian interface using the pressure from the Eulerian grid. First, the
pressure at the interface is determined using linear interpolation from the
17.3. Treating the Interface 207




18
Incompressible Flow


velocity field


0.7



0.6



0.5



0.4
18.1 Equations
0.3
Starting from conservation of mass, momentum, and energy, the equations
for incompressible ¬‚ow are derived using the divergence-free condition ·
0.2
V = 0, which implies that there is no compression or expansion in the ¬‚ow
¬eld. The equation for conservation of mass becomes
0.1

ρt + V · ρ = 0, (18.1)
0.1 0.2 0.5 1
0.3 0.7
0 0.6
0.4 0.9
0.8
indicating that the (possibly spatially varying) density is advected
along streamlines of the ¬‚ow. The Navier-Stokes equations for viscous
Figure 17.1. A shock wave propagating through a gas bounded on top and bottom
incompressible ¬‚ow are
by Lagrangian materials with strength. (See also color ¬gure, Plate 10.)
(2µux )x + (µ(uy + vx ))y + (µ(uz + wx ))z
px
ut + V · u+ = , (18.2)
ρ ρ
(µ(uy + vx ))x + (2µvy )y + (µ(vz + wy ))z
py
vt + V · v+ = + g,
ρ ρ
(18.3)
(µ(uz + wx ))x + (µ(vz + wy ))y + (2µwz )z
pz
wt + V · w+ = , (18.4)
ρ ρ
where µ is the viscosity and g is the acceleration of gravity. These equations
are more conveniently written in condensed notation as a row vector
T
· „)
p (
Vt + V · V+ = + g, (18.5)
ρ ρ
210 18. Incompressible Flow 18.2. MAC Grid 211

where “T ” is the transpose operator, g = 0, g, 0 , and „ is the viscous In order to update, u, v, and w on the appropriate cell faces, equations
stress tensor (18.2), (18.3) and (18.4) are written and evaluated on those appropriate
« 
cell faces. For example, in order to discretize the convective V · u term at
2ux u y + v x uz + wx
„ = µ  uy + v x v z + wy  , xi± 1 ,j,k we ¬rst use simple averaging to de¬ne V at xi± 1 ,j,k ; i.e.,
2vy (18.6) 2 2
u z + wx v z + wy 2wz
vi,j’ 1 ,k + vi,j+ 1 ,k + vi+1,j’ 1 ,k + vi+1,j+ 1 ,k
2 2 2 2
vi+ 1 ,j,k = (18.11)
which can be expressed in a more compact form as
4
2
«  « T
u u and
„ = µ v  + µ v. (18.7)
wi,j,k’ 1 + wi,j,k+ 1 + wi+1,j,k’ 1 + wi+1,j,k+ 1
w w 2 2 2 2
wi+ 1 ,j,k = (18.12)
4
2
The incompressible ¬‚ow equations model low-speed ¬‚uid events, includ-
de¬ne v and w, while u is already de¬ned there. Then the V · u term on
ing many interesting liquid and gas ¬‚ows observed in everyday life. A
the o¬set xi± 1 ,j,k grid is discretized in the same fashion as the V · ρ
number of classic texts have been written about both the analytical and nu- 2
term on the regular xi,j,k grid, for example with third order accurate
merical approaches to these equations, including Landau and Lifshitz [101]
Hamilton-Jacobi ENO. The convective terms in equations (18.3) and (18.4)
and Batchelor [11]. A rather inspiring collection of photographs obtained
are discretized similarly.
from various experiments was assembled by Van Dyke [169]. A standard in-
The viscous terms are discretized using standard central di¬erencing. For
troduction to numerical methods for the Navier-Stokes equations is Peyret
example,
and Taylor [133], which discusses both the arti¬cial compressibility method
of Chorin [46] and the projection method of Chorin [47]. Some of the most
ui+ 1 ,j,k ’ ui’ 1 ,j,k
popular modern-day numerical methods were introduced by Kim and Moin 2 2
(ux )i,j,k = , (18.13)
x
[96], Bell, Colella, and Glaz [13], and E and Liu [56]. A rather intriguing
ui+ 1 ,j+1,k ’ ui+ 1 ,j,k
look at all three of these methods was recently presented by Brown, Cortez, 2 2
(uy )i+ 1 ,j+ 1 ,k = , (18.14)
and Minion [19]. y
2 2



and
18.2 MAC Grid ui+ 1 ,j,k+1 ’ ui+ 1 ,j,k
2 2
(uz )i+ 1 ,j,k+ 1 = (18.15)
z
2 2
Harlow and Welch [79] proposed the use of a special grid for incompressible
are used to compute the ¬rst derivatives of u. The functions ρ and µ are
¬‚ow computations. This specially de¬ned grid decomposes the computa-
de¬ned only at the grid points, so simple averaging is used to de¬ne them
tional domain into cells with velocities de¬ned on the cell faces and scalars
elsewhere, e.g.,
de¬ned at cell centers. That is, pi,j,k , ρi,j,k , and µi,j,k are de¬ned at cell
centers, while ui± 1 ,j,k , vi,j± 1 ,k , and wi,j,k± 1 are de¬ned at the appropriate µi,j,k + µi+1,j,k
2 2 2
µi+ 1 ,j,k = (18.16)
cell faces.
2
2
Equation (18.1) is solved by ¬rst de¬ning the cell center velocities with
and
simple averaging:
µi,j,k + µi+1,j,k + µi,j+1,k + µi+1,j+1,k
ui’ 1 ,j,k + ui+ 1 ,j,k
µi+ 1 ,j+ 1 ,k = . (18.17)
2 2
ui,j,k = , (18.8) 4
2 2
2
vi,j’ 1 ,k + vi,j+ 1 ,k Then the second-derivative terms are calculated with central di¬erencing
2 2
vi,j,k = , (18.9) as well. For example, (µ(uy + vx ))y in equation (18.2) is discretized as
2
wi,j,k’ 1 + wi,j,k+ 1
µi+ 1 ,j+ 1 ,k (uy + vx )i+ 1 ,j+ 1 ,k ’ µi+ 1 ,j’ 1 ,k (uy + vx )i+ 1 ,j’ 1 ,k
2 2
wi,j,k = . (18.10) 2 2 2 2 2 2 2 2
2 (18.18)
y
Then the spatial derivatives are evaluated in a straightforward manner, for
example using third order accurate Hamilton-Jacobi ENO. at xi+ 1 ,j,k .
2
212 18. Incompressible Flow 18.4. Poisson Equation 213

18.3 Projection Method needs to be satis¬ed when the boundary condition on V is speci¬ed in
order to guarantee the existence of a solution. Here, “ represents the
The projection method is applied by ¬rst computing an intermediate boundary of the computational domain and N is the unit normal to that
velocity ¬eld V , ignoring the pressure term boundary. See, for example, [133] for more details.
T
V ’Vn · „)
(
+ V· V= + g, (18.19)
t ρ
18.4 Poisson Equation
n+1
and then computing a divergence-free velocity ¬eld V ,
In order to update the intermediate velocity V obtained with equation
V n+1 ’ V p
(18.19) to the divergence-free V n+1 using equation (18.20), we need to ¬rst
+ = 0, (18.20)
t ρ
¬nd the pressure by solving equation (18.21). This equation is elliptic, since
using the pressure as a correction. Note that combining equations (18.19) the incompressibility condition is equivalent to assuming an in¬nite speed
and (18.20) to eliminate V results in equation (18.5) exactly. of propagation for the sound waves. This elliptic treatment of the acoustic
Taking the divergence of equation (18.20) results in waves gives a CFL condition that depends only on the ¬‚uid velocity; i.e., the
sound waves do not play a role. Note that one should include the viscosity
·V
p
and gravity terms in the time-step restriction as well, unless, of course, the
· = (18.21)
ρ t viscosity is treated implicity.
The right-hand side of equation (18.21) can be evaluated at each cell
after we set · V n+1 to zero, i.e., after we assume that the new velocity
center using central di¬erencing; for example,
¬eld is divergence free. Equation (18.21) de¬nes the pressure in terms of the
value of t used in equation (18.19). De¬ning a scaled pressure of p = p t ui+ 1 ,j,k ’ ui’ 1 ,j,k
2 2
(ux )i,j,k = (18.27)
leads to
x
p
V n+1 ’ V + =0 (18.22) is used to compute ux . The components of the p/ρ term are computed
ρ
at cell faces, for example
and pi+1,j,k ’ pi,j,k
(px )i+ 1 ,j,k = (18.28)
p x
2
· ·V
= (18.23)
ρ
is used to compute the pressure derivative at xi+ 1 ,j,k . This discretization of
2
in place of equations (18.20) and (18.21), where p does not depend on t. the pressure derivatives is used both in equation (18.20) to update the ve-
When the density is spatially constant, we can de¬ne p = p t/ρ, leading
ˆ locity ¬eld and in equation (18.21) to ¬nd the pressure. In equation (18.21)
to a second derivative is needed as well. Once again, central di¬erencing is
used; for example,
V n+1 ’ V + p=0
ˆ (18.24)
(px /ρ)i+ 1 ,j,k ’ (px /ρ)i’ 1 ,j,k
and 2 2
((px /ρ)x )i,j,k = (18.29)
x
·V ,
p=

<< . .

 18
( 23)



. . >>

Copyright Design by: Sunlight webdesign