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
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., . 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.,  and the ref-
problems: One can treat both the solid and the п¬‚uid with Eulerian numer-
erences therein, speciп¬Ѓcally Hancock , Noh , and McMaster ,
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 , 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.  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. .
While there are a number of sophisticated Lagrangian programs, we par-
ticularly like the approach of Caramana et al.  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 , , and . 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 , , With boundary conditions speciп¬Ѓed on both the Eulerian and Lagrangian
and . 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  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
(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 
Hamilton-Jacobi ENO. The convective terms in equations (18.3) and (18.4)
and Batchelor . A rather inspiring collection of photographs obtained
are discretized similarly.
from various experiments was assembled by Van Dyke . 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 , which discusses both the artiп¬Ѓcial compressibility method
of Chorin  and the projection method of Chorin . 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
, Bell, Colella, and Glaz , and E and Liu . 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 . 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  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,  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)
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)ОГЛАВЛЕНИЕ След. стр. >>