LINEBURG

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=

Copyright Design by: Sunlight webdesign