- Email: [email protected]
Recent extensions to Eulerian methods for numerical fluid dynamics
RECENT EXTENSIONS TO EULERIAN METHODS FOR NUMERICAL FLUID DYNAMICS
F. H.
HARLOW, C . W . IITRT (Los Alamos, USA)
The paper is a review of Eulerian methods for numerical fluid dynamics carried out at the Los Alamos Scientific Laboratory (USA) under the direction of well-known scientists . A new solution technique for the difference equations is suggested . There are considered flows at very low Reynolds numbers as well as the problems of atmospheric dynamics, turbulence models in numerical fluid flow calculations and so on . § 1 . Introduction Eulerian techniques for numerical fluid dynamics are characterized by a fixed mesh of computational cells . Relative to these, the fluid motions are calculated by means of finite-difference statements of the principles of mass, momentum and energy conservation . In contrast, Lagrangian methods apply those same principles to a computational mesh of cells that move with the fluid. For each viewpoint there are one or several basic versions, together with numerous extensions that have been designed to enhance applicability and improve accuracy [`] . Our purpose in this paper is to describe some of the most recent of these that have been developed at the Los Alamos Scientific Laboratory for the Eulerian viewpoint ., and in some cases to compare the capabilities of these extensions with alternative versions developed elsewhere . We assume that the basic features of Eulerian methodology are well known to the reader in the forms exemplified by the Fluid-in-Cell (FLIC) method for high-speed flows [ 2 ] and by the Marker-and-Cell (MAC) method for incompressible flows ['] . § 2 . A new solution technique for the difference equations In the MAC technique for two-dimensional, incompressible fluid flows, we require the simultaneous solution of three difference equations for pressure, p, x-direction velocity, u, and y-direction velocity, v . In their simplest form, these are
* Zh . vychisi . Mat . mat. Fiz ., 12, 3, 656-i72, 1972 .
123
F . H . Harlow and C . W. Hirt
1 24
istl
nti
n+l
u~+'!'l .l C U,-%,i
1
Vi.i+Yc - V,,
(1)
ay
T
nti n (u$)1 .1-(u')i+I .i + ui+'¢ .j = "++`S i +atI ax
(2)
+ (uv) :+k .i-''.,-(nv)i+'l: .i+1= + Pi .j - Pi+ia + pax by 2u ;, .'1a ,f 1 U,+'laj+ ui-%,i - 2 u=+'6d + uj+Ya .1+i + u,+'h.i-i + v
ay`
ax` n vf.l+'(s
V
+ at{(,
+
(UV
ax
+ (VI ) i .i(y"') :,j+l + Pio- Pi.i +1 +
ay
pay 2vi .j+'&
vi+I .i, Yj +
+
v
6x 2
;3+'!,+ vl,j-ii, - 2vi v by` +
j+'G
(3)
The new procedure, which has been given a preliminary description in several places [' _], will be generalized in a forthcoming paper by Hirt [ s ] . For the case at hand, the ideas can be illustrated by noting that large parts of the two momentum equations depend only on quantities available at the beginning of each computational cycle . Gathering these together into terms that are labeled with superscript tildes, u, l,j and v,, we can then write Eqs. (2) and (3) as follows +1
St pOx
(4)
u,+Yr .j=
(5)
Vf,j+'/. _ giJ+'IA +-(Pf" -Pi .i-I)
.+1
St
pay
If we define p~ Du = (6)
j - ul_Y,j
V1,j+% - vi .l-Y:
+
ax by Then the combination of Eqs . (1), (4) and (5) shows that the pressures can be obtained by solving St
St
ay (P"j+I+R P (7) Pf+a ;+Pf-le - 2ps>i)+ , per
I -
2pi,i)=D-i •
The simplicity of this approach is enhanced if one notes from Eqs . (4) and (5) that vanishing of the tilde normal velocity at a rigid wall results in
Recent extensions to Rulerimt methods
125
art especially simple homogeneous boundary condition for the pressure there . In the original description of this simplified procedure ['], it was noted that the tilde velocities carry the correct vorticity for the end of the cycle (except at vertices lying on boundaries) but that they fail to satisfy the conservation of mass condition (vanishing divergence) . From that viewpoint, it was sufficient to utilize a potential function for the conversion of tilde velocities to final velocities, the object being to create a divergence-free field while simultaneously preserving vorticity . This is now recognized to be equivalent to the procedure described her ,,, this present version having the advantages, of a direct interpretation for the potential function (namely p / p) and of direct generalization to a large class of similar problems, for example, to situations in which it is desired to have advanced time quantities in the convection of viscous terms, as described by Hirt [°] The essence of Hirt's generalization is the use of a specialty modified Newton's method for tlu' iteralive solution . The system if equations is split into those parts that enter into the solution procedure and those parts (like the tilde velocities above) that are determined front data available at the beginning of the cycle . Also, the velocities and pressures are simultaneously iterated, as in reference [`] . This will receive a detailed description, with examples, in the forthcoming paper by Hirt [ a] . § 3 . Flows at very low Reynolds number The methodology developed for this purpose by Pracht ['] is illustrated by its application to the dynamics of the Earth's mantle material where Reynolds numbers are on the order of 1D` . The purpose of this investigation was to develop a numerical tool for the utilization of available geophysical data to formulate a realistic model of the Earth's interior, The numerical problems in such an investigation are several-fold . The Navier-Stokes equations for an incompressible fluid require the following features of generalization : 1 . A heat transport equation to describe the effects of convection, diffusion, radioactivity sources, and dissipation of kinetic energy . (This last, usually negligible for incompressible flows, is appreciable because of the extremely large viscosity of the mantle material .) 2. A viscosity coefficient that varies appreciably as a function of both temperature and pressure, exhibiting a minimum (corresponding to maximum fluidity or mobility) at a depth of around five hundred kilometers . 3 . A viscosity density that depends strongly on depth (through the compressional effects of pressure) and on temperature (through both volumetric expansion and the effects of phase transitions) . Accordingly, the Boussinesq approximation for the buoyancy effects must be generalized considerably, the neglect of inertial terms being allowed only for the perturbations from an
F. H. Harlow and C . W. Hirt
1 26
equilibrium density distribution, rather than from merely a constant density . The numerical aspects of the solution procedure incorporate a variety of new features . One of these is the requirement for an implicit formulation of the viscous terms, such that all velocities in the strain-rate tensor are at the advanced time . Since the Eulerian changes with time and the advection terms are negligible in the momentum equations, the problem reduces to that of maintaining a perpetual stress balance among the pressure gradient, the viscous shear, and the buoyancy . This last varies because of the transient properties of the heat flow, in which the Eulerian time variations and the convection terms are of crucial importance . Another matter is the potential numerical instabilit introduced by the large spatial variations in viscosity . Difficulties arise i achieving convergence of the iterative solution for the velocities whenever ax O v
v 8x
or
ay av
v
ay
anywhere exceeds a value near unity . To mitigate this difficulty has required the introduction of an under relaxation in the iterative solution of the equations . With the successful incorporation of these features, it has been possible to investigate a variety of geophysical models . A typical computational region contains a slice of Earth along a two-dimensional latitudinal cut through central South America, several thousand kilometers long and nearly oae thousand kilometers deep . The detailed methodology, description, and results will be presented elsewhere [e] . § 4 . Atmospheric dynamics Global modes for meteorological studies have utilized a variety of Eulerian techniques . Some have been purely two-dimensional, with coordinates along latitude and longitude, the vertical effects being modeled through an equilibrium assumption . Others employ a small number of vertical layers, but also usually utilize an equilibrium assumption to simplify the procedure . Our goal, however, has been to calculate a variety of much smaller scale meteorological processes, in which the lateral scale is comparable with the depth of the atmosphere . Accordingly, for two-dimensional investigations we utilize a plane or cylindrically-symmetric section in which the vertical zoning is as fine as the horizontal . For the investigation of localized storms (thunderheads or tornados, for example) it is necessary to extend the calculational region to as high as ten kilometers, over which the equilibrium density varies by approximately a factor of three . To model the effects of intense explosions
Recent extensions to Euieriar. methods
127
in the atmosphere may require extension of the calculation region to as high as one-hundred kilometers or more, and moreover introduces the necessity to account for the effects of high-speed compressibility . For this last, we utilize the newly-developed ICE technique described below . For investigations of most meteorological problems, even including the most violent storms, the air flows are sufficiently subsonic as to he classed as cincompressibleti . Nevertheless, the variations in atmospheric density must be fully taken into account, including those produced by advection, temperature changes and water vapor variations . This means that the MAC method, for example, requires some significant modifications for applicability . In addition, an azimuthal velocity must be included for such cyclonic flows as are found in tornados, water spouts or dust devils . More complicated, however, than any of the extensions described so far, is the requirement for realistic representations of turbulence effects . Our version of this is described below . Aside from turbulence, the principal modifications, those of heat transport, buoyancy, and density layering, are accomplished in much the same way as in the geophysical calculations . The general Navier-Stokes equations are specialized to circumstances in which the density varies only slightly from that which is in gravitational equilibrium with an initial temperature distribution . Inertial effects of the perturbation are neglected, only the buoyancy effects are included in a manner that extends slightly the usual lloussinesq assumption . Our calculations have examined the concentration of vorticity during tornado formation and the distortion of a thunderhead plume rising from heated ground in a cross wind ['] . § 5 . An implicit, continuous- fluid, culerian (ICE) extension to span all flow speeds fit many problems, it is necessary to allow for both low and high speed flows, either simultaneously in adjacent spatial regions or sequentially in the same region. Closely related, and perhaps even more frequently encountered, are fluid flow problems in which the flow speeds are supersonic in one direction and far subsonic in an orthogonal direction . For any of these and similar circumstances, it is necessary to have a single formulation of numerical techniques capable of spanning the entire range of Mach number, M, from zero (the #incompressible* limit) to much greater than unity (the hypersonic limit) . For this purpose we have developed a n implicit formulation of the equations that reduces to the MAC method for M =0 and becomes a variant of the FLIC technique for M % 1 .0 . The essence of the ICE method is that it isolates those parts of the equations related to density variations and treats them purely implicitly . While
F . H. Harlow and C . W . Hirt
128
most of our applications have been to multidimensional problems, ['0] it is convenient to illustrate the technique for a simple one-dimensional configuration, for which the differential equations are apu 2 ax
apu Op apu ax at +3t 89E+BauE='
a [(p+
a(P+q) ax
For the viscosity function, q, we choose for illustration au ax it is convenient to define,
x
R,+I = 6
[n++'1=(P+u+-'1
-P++n++ •, ) + qc - qi++],
so that the difference equations for mass and momentum become ,+1
n
ht
8 = Ox
[(Pu}+-'=-(Pn) ;+v ]~
1
+
8x0
[(pu)Ij+-(PU)+ w]
and (P(pu)+v, 6t ~
1-¢ „ 0 (P"-P,+t)+R++Y+ • _((pPi;-Pi++)~ 6x
The hybrid function, p, is formed from the equation-of-state pressure function, p(p, I), in a manner similar to that used by Janenko and Neuvazaev ["] for a somewhat different purpose . P< =
pi ,
+ ei n (P
,n+1-
Pi'),
where c ;" - (Op l dp) i" is related to the speed of sound. (With a constant specific heat, ci is the square of the isothermal sound speed .) Thus Pi singles out from the equation of state its principal dependence upon the density and for it uses the advanced-time density, while the rest remains at the beginning-of-cycle values . Without this procedure, i.e., with c;" = 0 and ; or with 0 = ¢ = 0, the procedure reduces to an explicit Eulerian treatment, in which pressure depends only on the previous-cycle values of p and 1 . Note the occurrence of weighting constants, 0 and ~, by which the technique can be varied from purely explicit (8 = ~ = 0) to partially time cen-
Recent
extensions to Eulerian methods
129
tered (0 = b = 0.5) to advanced-time implicit (0 = ,~ = 1 .0) . To solve( ' these equations, we eliminate the advanced-time mass fluxes, which ate the same as the momentum densities, and thereby derive an equation for the hybrid pressure function, p . To exhibit this, it is convenient to rather together all of the beginning-of-cycle quantities by defining the function G, as follows G; = P' C'"
+ bt dx
+ 0 (1- 0) 0t2
[(Pu)1n 9.-(Pu) ;+61+
n
Sx'
estz
6x
- An+vJ
in terms of which the equation to be solved is 46v (Pi-, + pi+, - 2P) -{- G ; . 6X' This, then, describes the essential features of the ICE technique . The steps for each computational cycle can be summarized as follows : 1 . Calculate the required beginning-of-cycle terms, R; + G, ; 2. Solve for pi for every cell. For some circumstances a direct solution method can be used ; for others it may be necessary to employ a relaxation technique ; 3. Solve for pi+t using the definition of p, ; 4. Solve for the new values of the velocities ; 5 . Calculate the new energy for every cell, using an appropriate finitedifference approximation to the energy equation . There are various wags that the energy and momentum convection fluxes can be written, the precise form being irrelevant to the conceptual formulation of the ICE technique, but nevertheless crucial to the assurance of accuracy and numerical stability of the calculations . At least four types of the flux expressions are appropriate for consideration, each having advantages and disadvantages . Consider, for example, the momentum flux term, (pa-), . While the density is already stored for the position at which the flux is to be evaluated (that is, at the center of cell i) the velocity positions are at i± t/2, so that some type of interpolation is needed . Four flux expressions that could be considered are Centered : (Pu 2 )i- pi[(ui-'h + ui+ ;.) /21'' . ZIP : (pa 2 )i -' Piui-I u,-'G~
F . H. Harlow and C. W . Hirt
1 30
Partial Donor: Ott)
ut+%) /2, { piui-'h (u,-k + p .u=+%(u--'i,+ui+S)/2 .
(U'_'/' + u,+ss) > 0, (u_v.+ut+v,)<0.
Complete Donor. )t (pa')
z .,'h,
(ui-y --u++'+ h) > 0, (u+-'ti -f' ut+'G) < 0 .
In addition there is an interpolated donor form, not considered here. For the energy-equation convection terms there are analogous flux expressions for two of these forms : Centered : '/s (p+ + p ;+a) u++v, (Et + Ei+, ), (PuE) i+ Donor : p;u;+9,Et, ut+a > 0, (P nE ) ;+u, ui+ { P ;+,u++wE++t,
< 0.
Some properties of each of these are the following : 1 . Centered - This form eliminates one order of Sx truncation error, compared with the donor cell flux, but this advantage is overshadowed by the tendency to numerical instability that is produced by not time centering the convective flux . This tendency to instability is usually mitigated by the addition of a counterbalancing diffusive flux (either real or artificial) and car) also be tolerated in circumstances in which the growth of the instability is bounded by non-linear dissipative effects . 2. ZIP-This form is also one order in 8x more eaccurate* than the donor cell flux, but in addition has the advantages of eliminating a nonlinear contribution to instability ['2] . For these reasons, the ZIP flux is often preferred for the momentum convection . The ZIP flux also needs a counterbalancing diffusive flux when it is not centered in time . 3 . Donor-Neither the partial nor the complete donor cell flux is space centered, and accordingly both have low order (in Ox) truncation errors . These contribute a positive diffusive effect and tend to automatically stabilize the numerical calculations . The magnitude, however, can be excessive, thereby masking the true diffusive effects and leading to erroneous interpretations (for example, for flows at high Reynolds numbers) . The once-popular donor-cell flux representation is being abandoned by most invertigators except for very special circumstances . Our recommendation for ICE calculations is to use ZIP fluxes or centered fluxes, the instability of these mitigated by the controlled use of artificial diffusion with a coefficient that just slightly exceeds the effective value pre-
Recent extensions to Eulerian methods
1 31
dieted by truncation error analysis . The only circumstance in which a donorcell flux appears desirable is in the early stages of problems with especially violent initiation (for example, in high-speed projectile impact), the later stages of the calculation automatically shifting to the centered or ZIP flux . To show some of the ICE -method properties, we may expand the difference equations in a form that reveals the lowest order diffusional truncation effects . With ZIP differencing of the momentum convection terms, for example, the result is Op
8pu
r)t
OX
8u l
au
6x2 0 5 2p 8t =1 (20-1)(c +u`) 2 4 dx
I ap
+a-+0x p (2¢-1)c~
dx 3uE ]
6t 2
usx= ap 8p
Ox
~.
+ p 1
a° u OX
We have used c ° to denote either of several partial derivatives of p with respect to p ; the symbol as used in this section is only approximately related to the square of the sound speed, being progressively more so in the limit as C'- W .
Several classes of effects can be noted from these equations . For r2 -. 00 (the incompressible limit), the effective diffusion coefficient becomes large and negative if 0 and/or ~ is less than 0 .5, corresponding to intolerable mimerical instability . In contrast, for 0 and ¢ both greater than 0 .5 the diffusion positive (hence stabilizing), and the equations suggest that this could become intolerable as cz -> 00 . This excessive diffusion, however, is purely longitudinal, so that two-dimensional incompressible flow calculations with the ICE method ordinarily will not suffer from the effect, since longitudinal signal propagation is extremely rapid, anyway . Despite the implicit nature of the ICE method, there are still some restrictions that must be observed in order to assure numerical stability . These are difficult to describe precisely ; the matter is discussed more fully in ['°] . § 6 . Turbulence models in numerican fluid flow calculations There is considerable reason for believing that the most important effects of turbulence depend only on some of the collective effects of the fluctuations in a way that is independent of the detailed structure of the turbulent flow . Because of the widespread occurrence and profound effects of turbulence there has been a strong effort to discover how to represent these col-
F . H. Harlow and C. W . Hirt
132
lective effects in a useful way for predicting turbulence and for analyzing its numerous important effects . One of the results of these investigations has been the relatively new idea of the transport-equation model . Several versions of this have been proposed and proof tested at Los Alamos, with results that are highly encouraging . These results, together with a bibliography of earlier work are given in references P ." "] To illustrate the central problems involved in these investigations, consider the following equation
0R,;
OR„; 0k
-- -1- U,
ox,
(it
1a,(' - ( ax, +
On,
R*A+ ax,,
au,'¢ \'
on, 1j
X,
aa ;
r3x, I +
8 dx,
Rla-t-(ui'uJ au ;
( dx ; + 3 %x,
) +
z
" - 4D;,} -
v ~-
T"
ax„'
(g;ui T' + g,ui T') .
This equation for the Reynolds stress tensor, R,;, is a rigorous consequence of the Nr avier- Stokes equations for an incompressible fluid, To derive it, the velocity of a fluid element, U, is considered to be the sum of a mean-flow part . a,, and a fluctuating part, a,' . Likewise, the pressure, ~, and the temperature, T, are split into the mean part (designated in multi-factor terms by an overbar) and a fluctuating part (designated by a prime) . With a formulation normalized to unit density, the Reynolds stress is R;; = u ;'a ;' . The other quantities in the equation are g„ the gravitational acceleration ; To, the reciprocal bulk expansion coefficient in the Boussinesq approximation to the buoyancy term ; v, the molecular kinematic viscosity coefficient ; and
1
a u;' 8u ;
2
ax„
ax,
which describes the decay rate of Reynolds stress . It is likewise possible to derive a rigorous transport equation for D, ; from the Navier - Stokes equation ['"] . Although this is a fundamental part of the theory, it is omitted here, not being necessary for illustration of the principles that are, to be demonstrated in the present paper . In the form presented above, the Reynolds stress equation cannot be solved because of the presence of various moment terms (the overbar terms and D ;) that are not related to the available variables . The central problem in turbulence modeling, therefore, in to discover appropriate ways for representing these terms as functionals of a set of variables that is no more numerous than
Recent extensions to Eulerian methods
133
the set of equations that can be found to desci ibe their variations . One form for such an equation, derived in["] . introduces a scale function, s, together with various universal, dimensionless constants and functions . a, 0, r, to, 0, 4 in such a way that oR i;
du;
oR,
`+UAW-
fdl
+
ox,
UX,f
1~ ~ ;~,
d 8
d'R,, ~ , r)X,
[
du, ox,
`R,A--'-"R ;A-+
s-q
\'
OR,A
r)x,',vi\
2v4 s-
dx.A / 'US'
-R1
tj
dx,
-(g,R,A
(3TA .N
d s R, 'R„ a-v,
-J 11,
I
ox,
q o R,A A A' .A
dx .~,
/1
oT ± g;A ; ;J-,
where ~3 is the thermmuetric condur-tivit.v, and (D„ = - ~,~ (R„ - 1 ,i s' aU
RA
,b n l -
UIE ; }
I oU A
x„ ,
R ' `a,r,
tz~R,4~dx ;
du r t~ x,.~~
if the scale function is specified or calculated, then this tensor equation for the Reynolds stress is complete, and its simultaneous solution with the mean-flow equations determines the flow . To specify the scale function, ona could resort to simple mixing-length theories, or else utilize a second transport equation . Reference [ "] discusses this latter, and shows how this can be accomplished by means of moment approximations in the transport equation for D,,. Several properties of a flow must he present in order for a transport-equation model to be valid . In a division of the flow into mean and fluctuating parts, the spectrum of the fluctuations must not strongly deviate from some universal, self-similar, equilibrium spectrum . There simply are not enough degrees of freedom in a low-order moment representation to describe the effects of a complicated spectrum that varies by large amounts from case to case . Indeed, there must be some cases in which strong departure from equilibrium removes all validity from the low-order moment equations . At present, however, we have very little information as to how far from equilibrium a turbulent flow can be and still be capable of description by a transport model, For many kinds of investigations it will be necessary to treat both the turbulence and the mean flow as time varying, the whole being so complicated that to get solutions will require numerical methods with high speed computers . At present, there are investigations in progress to determine just how extensively the turbulence transport models by themselves can he applied to flows with superimposed large-scale fluctuations, so that the matter is still far from settled . What is clear is that for truly one scale turbulence problems, the
1 34
F . H . Harlow and C. W . Hirt
present models for Reynolds stress transport are capable of producing excellent agreement with experiments . Whenever there is a need for simultaneous coupling of the mean-flow equations to the turbulence transport calculations, we may expect one of the high-Reynolds-number difficulties usually associated with finite-difference calculations to be mitigated by the diffusive properties of the effective eddy viscosity . Numerical stability considerations usually restrict the Reynolds number to magnitudes less than the square of the number of computational cells required to resolve a typical linear dimension of the flow field . For laminar flows at high Reynolds number, it is accordingly necessary to introduce an artificial viscosity (either explicitly or else as a consequence of the differencing technique) whereby the effective Reynolds number is reduced to the required magnitude . In many cases the results are nevertheless satisfactory because of insensitivity of the laminar flow field to Reynolds number variations when the Reynolds number is sufficiently large . For fully turbulent flows, however, this decrease in effective Reynolds number, resulting from the eddy viscosity can be sufficient to preclude numerical instability, even when the resolution is sufficiently modest for practical computation .
§ 7 . Curved surfaces In Eulerian fluid-dynamics calculations, the occurrence of curved boundaries requires special treatment . A variety of examples come to mind : 1. Free surfaces ["] . 2. Moveable interfaces between materials ["] . 3. Curved rigid walls 4. Moveable walls ["] . In every case, the principal difficulty arises from the partial cells created by the boundary . For free surfaces, it is necessary to balance both the normal and tangential stresses . Experience shows that the original version of the MAC method ['] is not sufficiently accurate in resolution of the surface configuration for problems with large viscosity and / or low amplitude waves . To alleviate the difficulty, it has been found useful to introduce surface marker particles, and to apply the stress conditions at the actual surface position rather than at the cell center ["] . For moveable interfaces between materials, a similar marker-particle set has also been useful as a means for determining how to better resolve the continuity of normal and tangential stresses ["] . Curved rigid walls present a somewhat different challenge that has been solved by Viecelli ["] in an intriguing manner . He introduces a pressure boundary condition that accelerates the fluid in exactly that manner that is re-
Recent extensions to Eulerian methods
135
quired to keep the fluid velocity tangential to the wall . For moveable walls, he uses a similar technique [f 8 ], in which the normal velocity is adjusted to the prescribed magnitude, rather than vanishing . A somewhat different approach to the problem of curved surfaces (fixed or moving) has been taken by Hirt and Amsden [10 ] . They have devised an Arbitary-Lagrangian-Eulerian (ALE) technique, in which the computer mesh can vary from completely fixed (but not necessarily with rectangular cells) to moveable with the fluid velocity . With this, they can resolve curved and /or moving wall configurations with the Lagrangian capability, and at tho same time resolve the strong internal distortions of the fluid with the Eulerian capability .
Fig. 1
1 36
F . H. Harlow and C . W. Hirt
Fig . 2
§ 8. Theree space dimensions Time-varying problems in three space dimensions are easily solved with Eulerian techniques, especially when the computer program is based on the efficient new iteration method of Hirt [°] . Recent calculations by Hirt and Cook ['°] have shown both efficiency and accuracy in the investigation of at • mospheric wind flows over variable terrain and around buildings . Extensions that include heat transport, buoyancy effects and smoke transport have demonstrated a capability for the realistic analysis of atmospheric pollution studies, and the techniques are readily extendable to underwater pollution and sediment transport and deposition investigations . § 9 . Sample calculations A variety of sample calculations are described here to illustrate some of the techniques described in the preceding sections .
S 2TJ
N
I PV
£ t?3
al
spoil aw
U
Aajng o7 suo]suapxa lua ag
138
F . H . Harlow and C . W. Hirt
Fig . 6
Fig . 1 illustrates the evolution of a block of viscous material such as far, resting on a table and sagging tinder the influence of gravity . This calculation was performed using the technique described in § 3 and ['] . These
Recent extensions to Eulerion methods
139
pictures show the configuration of the marker particles used in the calculation to delineate regions occupied by fluid . The first frame a is of the initial block shape while
b
through
d
show subsequent stages in the time
history of the block . Fig . 2 illustrates a type of problem that can be studied in detail with the improved Marker-And-Cell techniques mentioned in § 7 . The example illustrates what happens when a spherical drop of fluid falls into a container of fluid . The picture sequence follows the splash from the moment of initial impact a, through the point of maximum cavity formation c, and into the secondary splash resulting from cavity collapse
f.
Another new technique mentioned in Section 7 is the Arbitrary-Lagrangian-Eulerian method . The unique feature of this method is its ability to resolvo flexible boundary walls . This capability has been utilized to calculate the pulsating flow of blood through an elastic artery . Fig . 3 presents three kinds of data produced directly by the computer . In a the calculational mesh for a cross-section of the tube is shown with a two pulse disturbance visible on its right edge . The left edge of the mesh is an axis of symmetry . The mesh is continuously rezoned to maintain the radial grid lines in fixed positions, and the axial grid lines were required to intersect the radial lines with equal spacing between intersection points . This choice of rezoning maintains a uniform mesh pattern, which is necessary for accurate results . In b and c are shown velocity vectors and contour ones of constant pressure, both at the same lime as the grid configuration shown in a . It will be observed that maximum pressures occur under the surface bulges and minimum pressures under the surface indentations . corresponding to the application of elastic stresses in the confining tube . A variety of three-dimensional calculations have been carried out with a new version of the Marker-And-Cell technique . Fig . 4-6 illustrate the results of several calculations involving the incompressible flow of air about large rectangular structures . In Fig . 4 a uniform flow is entering the computational region (large rectangular box) front the left silo and, after flowing around a building (smaller rectangular box), is allowed to flow out the right side of the computational region . In a the velocity vectors in a layer of cells passing over the top of the building are shown in perspective, while b shows the velocities in a layer of cells near the bottom of the mesh . A large stagnation region behind the building is clearly visible in b . Fig . 5 shows the results of a similar calculation performed with a second, smaller, structure located behind the first one . Finall, Fig . B illustrates a distribution of particulate matter (air pollution) calculated, for the One and two building problems . In a there is a strong concentration of pollution in the, wake region behind the building,
140
F . H . Harlow and C. W . Hirt
while in b no such concentrated region exists, because the smaller structure deflects the flow across the near face of the large building and prevents the formation of a large recirculation region .
Translated by D . E . Brown
REFERENCES 1.
HARLOW, F . H. Numerical methods for fluid dynamics, an annotated bibliography . Los Alamos Scient. Lab. Rept. LA-4281, 1969 . GENTRY, R . A ., MARTIN, R . E . and DALY . B . J . An Eulerian differencing method for unsteady compressible flow problems . J . Comput . Phys ., 1, 87, 1966 .
3.
HARLOW, F . H . and WELCH, J . E . Numerical calculation of time-dependent viscous incompressible flow of luid with free surface, Phys . Fluids, 8, 2182 . 1965 ; WELCH, J . E ., HARLOW, F . H ., SHANNON, J . P . and DALY, B . J . The MAC method . Los Alamos Scient . Lab . Rept . LA-3425, 1966 . CHORIN . J . Numerical solution of the Navier - Stokes equations . Math . Comput, 22, 745, 1968 .
5.
AMSDEN, A . A . and HARLOW, F . H . The SMAC method . Los Alamos Scient . Lab. Rept . LA-4370, 1970 .
6.
HIRT, C . W . A general algorithm for the solution of implicit difference equations . Los Alamos Scient . Lab ., Group T-3, unpublished notes, 1971 .
7.
PRACHT, W . E . Implicit solution of creeping flows, with application to continental drift. Proc . See. Internat . Conference on Numer. Methods in Fluid Dynamics, Univ . California, Berkeley, California, Sept . 15-19, 1970 ; PRACHT, W . E . A numerical method for calculating transiet creep flows, J . Comput. Phys ., 7, 46, 1971 .
8.
PRACHT, W . E . and HARLOW, F . H . A numerical solution technique for the investigation of planetary interior dynamics . Los Alamos Scient. Lab ., Group T-3 manuscript, 1971 ; HARLOW, F . H . and PRACHT, W . E ., A theoretical study of geothermal energy extraction, J . Geophys . Res ., submitted for publication .
9.
GENTRY, R . A . and NAKAYAMA, P . 1 . Loa Alamos Scient Lab., unpublished manuscript, 1971 .
10 .
HARLOW, F . H . and AMSDEN, A . A ., A numerical fluid dynamics calculations method for all flow speeds, J . Comput . Phys ., 8, 197, 1971 .
11 .
YANENKO, N . N . and NEUVAZHAEV, V . E . A method of computing gas-dynamical movements involving non-linear heatconductions,Tr . Matem . in-ta Akad . Nauk
12 . HIRT, C . W . Heuristic stability theory for finite difference equations, J . Comput, Phys ., 2, 339, 1968 .
Recent extensions to Eulerian methods
141
13 .
DALY, B . J . and HARLOW, F . H . Transport equations in turbulence . Phys . Fluids, 13, 2634, 1970 .
14-
HARLOW, F . H . Turbulence models . HARLOW, F . H . and AMSDEN, A . A . Fluid dinamies, an introductory text . Sec . ed . Univ. California, Los Alamos Scient . Lab, Monograph, LA-4700, 1971 .
15 . CHAN, R . K .C, STREET, R . L. and STRELKOFF, T . Computer studies of finiteamplitude water waves . Stanford Univ . Techn . Rept. 104, 1969 ; CHAN, R . KC and STREET, R . L . A computer study of finite amplitude waves, J . Comput . Phys ., 6, 68, 1970 ; NICHOLS, B . D . and HIRT, C . W . Improved free surface boundary conditions for numerical incompressible-flom calculations, J . Comput . Phys ., 8 . 434, 1971 . 16 .
DALY, B . J . Phys . Fluids . 10, 297, 1967 ; DALY, B. J . and PRACHT, W . E . Numerical study of density-current surges, Plevs . Fluids, 11, 1, 1968 ; DALY, H . J . A technique for including surface tension effects in hydrodynamic calculations, J. Comput . Phys ., 1, 4, 1969 ; DALY, B . J, Numerical study of the effect of surface tension on interface instability, Phys . Fluids, 12, 1, 1969 .
17 .
VIECELLI, J . A . A method for including arbitrary external boundaries inn the MAC incompressible fluid computing technique, J . Comput . Phys ., 4, 543, 1969 .
18 . VIECELLI, J . A . A computing method for incompressible flows bounded by moving walls, ;J . Comput . Phvs ., 8, 119, 1971 . 19 .
HIRT, C . W . and AMSDEN, A . A ., An arbitrary Lagrangian _ Eulerian computing technique for all flow speeds, J . Comput Phys ., to be published.
20 . HIRT, C . W . and COOK, J . L. The calculation of three-dimensional flows around structures and over rough terrain, J . Comput . Phys ., accepted for publication .
F. H.
HARLOW, C . W . IITRT (Los Alamos, USA)
The paper is a review of Eulerian methods for numerical fluid dynamics carried out at the Los Alamos Scientific Laboratory (USA) under the direction of well-known scientists . A new solution technique for the difference equations is suggested . There are considered flows at very low Reynolds numbers as well as the problems of atmospheric dynamics, turbulence models in numerical fluid flow calculations and so on . § 1 . Introduction Eulerian techniques for numerical fluid dynamics are characterized by a fixed mesh of computational cells . Relative to these, the fluid motions are calculated by means of finite-difference statements of the principles of mass, momentum and energy conservation . In contrast, Lagrangian methods apply those same principles to a computational mesh of cells that move with the fluid. For each viewpoint there are one or several basic versions, together with numerous extensions that have been designed to enhance applicability and improve accuracy [`] . Our purpose in this paper is to describe some of the most recent of these that have been developed at the Los Alamos Scientific Laboratory for the Eulerian viewpoint ., and in some cases to compare the capabilities of these extensions with alternative versions developed elsewhere . We assume that the basic features of Eulerian methodology are well known to the reader in the forms exemplified by the Fluid-in-Cell (FLIC) method for high-speed flows [ 2 ] and by the Marker-and-Cell (MAC) method for incompressible flows ['] . § 2 . A new solution technique for the difference equations In the MAC technique for two-dimensional, incompressible fluid flows, we require the simultaneous solution of three difference equations for pressure, p, x-direction velocity, u, and y-direction velocity, v . In their simplest form, these are
* Zh . vychisi . Mat . mat. Fiz ., 12, 3, 656-i72, 1972 .
123
F . H . Harlow and C . W. Hirt
1 24
istl
nti
n+l
u~+'!'l .l C U,-%,i
1
Vi.i+Yc - V,,
(1)
ay
T
nti n (u$)1 .1-(u')i+I .i + ui+'¢ .j = "++`S i +atI ax
(2)
+ (uv) :+k .i-''.,-(nv)i+'l: .i+1= + Pi .j - Pi+ia + pax by 2u ;, .'1a ,f 1 U,+'laj+ ui-%,i - 2 u=+'6d + uj+Ya .1+i + u,+'h.i-i + v
ay`
ax` n vf.l+'(s
V
+ at{(,
+
(UV
ax
+ (VI ) i .i(y"') :,j+l + Pio- Pi.i +1 +
ay
pay 2vi .j+'&
vi+I .i, Yj +
+
v
6x 2
;3+'!,+ vl,j-ii, - 2vi v by` +
j+'G
(3)
The new procedure, which has been given a preliminary description in several places [' _], will be generalized in a forthcoming paper by Hirt [ s ] . For the case at hand, the ideas can be illustrated by noting that large parts of the two momentum equations depend only on quantities available at the beginning of each computational cycle . Gathering these together into terms that are labeled with superscript tildes, u, l,j and v,, we can then write Eqs. (2) and (3) as follows +1
St pOx
(4)
u,+Yr .j=
(5)
Vf,j+'/. _ giJ+'IA +-(Pf" -Pi .i-I)
.+1
St
pay
If we define p~ Du = (6)
j - ul_Y,j
V1,j+% - vi .l-Y:
+
ax by Then the combination of Eqs . (1), (4) and (5) shows that the pressures can be obtained by solving St
St
ay (P"j+I+R P (7) Pf+a ;+Pf-le - 2ps>i)+ , per
I -
2pi,i)=D-i •
The simplicity of this approach is enhanced if one notes from Eqs . (4) and (5) that vanishing of the tilde normal velocity at a rigid wall results in
Recent extensions to Rulerimt methods
125
art especially simple homogeneous boundary condition for the pressure there . In the original description of this simplified procedure ['], it was noted that the tilde velocities carry the correct vorticity for the end of the cycle (except at vertices lying on boundaries) but that they fail to satisfy the conservation of mass condition (vanishing divergence) . From that viewpoint, it was sufficient to utilize a potential function for the conversion of tilde velocities to final velocities, the object being to create a divergence-free field while simultaneously preserving vorticity . This is now recognized to be equivalent to the procedure described her ,,, this present version having the advantages, of a direct interpretation for the potential function (namely p / p) and of direct generalization to a large class of similar problems, for example, to situations in which it is desired to have advanced time quantities in the convection of viscous terms, as described by Hirt [°] The essence of Hirt's generalization is the use of a specialty modified Newton's method for tlu' iteralive solution . The system if equations is split into those parts that enter into the solution procedure and those parts (like the tilde velocities above) that are determined front data available at the beginning of the cycle . Also, the velocities and pressures are simultaneously iterated, as in reference [`] . This will receive a detailed description, with examples, in the forthcoming paper by Hirt [ a] . § 3 . Flows at very low Reynolds number The methodology developed for this purpose by Pracht ['] is illustrated by its application to the dynamics of the Earth's mantle material where Reynolds numbers are on the order of 1D` . The purpose of this investigation was to develop a numerical tool for the utilization of available geophysical data to formulate a realistic model of the Earth's interior, The numerical problems in such an investigation are several-fold . The Navier-Stokes equations for an incompressible fluid require the following features of generalization : 1 . A heat transport equation to describe the effects of convection, diffusion, radioactivity sources, and dissipation of kinetic energy . (This last, usually negligible for incompressible flows, is appreciable because of the extremely large viscosity of the mantle material .) 2. A viscosity coefficient that varies appreciably as a function of both temperature and pressure, exhibiting a minimum (corresponding to maximum fluidity or mobility) at a depth of around five hundred kilometers . 3 . A viscosity density that depends strongly on depth (through the compressional effects of pressure) and on temperature (through both volumetric expansion and the effects of phase transitions) . Accordingly, the Boussinesq approximation for the buoyancy effects must be generalized considerably, the neglect of inertial terms being allowed only for the perturbations from an
F. H. Harlow and C . W. Hirt
1 26
equilibrium density distribution, rather than from merely a constant density . The numerical aspects of the solution procedure incorporate a variety of new features . One of these is the requirement for an implicit formulation of the viscous terms, such that all velocities in the strain-rate tensor are at the advanced time . Since the Eulerian changes with time and the advection terms are negligible in the momentum equations, the problem reduces to that of maintaining a perpetual stress balance among the pressure gradient, the viscous shear, and the buoyancy . This last varies because of the transient properties of the heat flow, in which the Eulerian time variations and the convection terms are of crucial importance . Another matter is the potential numerical instabilit introduced by the large spatial variations in viscosity . Difficulties arise i achieving convergence of the iterative solution for the velocities whenever ax O v
v 8x
or
ay av
v
ay
anywhere exceeds a value near unity . To mitigate this difficulty has required the introduction of an under relaxation in the iterative solution of the equations . With the successful incorporation of these features, it has been possible to investigate a variety of geophysical models . A typical computational region contains a slice of Earth along a two-dimensional latitudinal cut through central South America, several thousand kilometers long and nearly oae thousand kilometers deep . The detailed methodology, description, and results will be presented elsewhere [e] . § 4 . Atmospheric dynamics Global modes for meteorological studies have utilized a variety of Eulerian techniques . Some have been purely two-dimensional, with coordinates along latitude and longitude, the vertical effects being modeled through an equilibrium assumption . Others employ a small number of vertical layers, but also usually utilize an equilibrium assumption to simplify the procedure . Our goal, however, has been to calculate a variety of much smaller scale meteorological processes, in which the lateral scale is comparable with the depth of the atmosphere . Accordingly, for two-dimensional investigations we utilize a plane or cylindrically-symmetric section in which the vertical zoning is as fine as the horizontal . For the investigation of localized storms (thunderheads or tornados, for example) it is necessary to extend the calculational region to as high as ten kilometers, over which the equilibrium density varies by approximately a factor of three . To model the effects of intense explosions
Recent extensions to Euieriar. methods
127
in the atmosphere may require extension of the calculation region to as high as one-hundred kilometers or more, and moreover introduces the necessity to account for the effects of high-speed compressibility . For this last, we utilize the newly-developed ICE technique described below . For investigations of most meteorological problems, even including the most violent storms, the air flows are sufficiently subsonic as to he classed as cincompressibleti . Nevertheless, the variations in atmospheric density must be fully taken into account, including those produced by advection, temperature changes and water vapor variations . This means that the MAC method, for example, requires some significant modifications for applicability . In addition, an azimuthal velocity must be included for such cyclonic flows as are found in tornados, water spouts or dust devils . More complicated, however, than any of the extensions described so far, is the requirement for realistic representations of turbulence effects . Our version of this is described below . Aside from turbulence, the principal modifications, those of heat transport, buoyancy, and density layering, are accomplished in much the same way as in the geophysical calculations . The general Navier-Stokes equations are specialized to circumstances in which the density varies only slightly from that which is in gravitational equilibrium with an initial temperature distribution . Inertial effects of the perturbation are neglected, only the buoyancy effects are included in a manner that extends slightly the usual lloussinesq assumption . Our calculations have examined the concentration of vorticity during tornado formation and the distortion of a thunderhead plume rising from heated ground in a cross wind ['] . § 5 . An implicit, continuous- fluid, culerian (ICE) extension to span all flow speeds fit many problems, it is necessary to allow for both low and high speed flows, either simultaneously in adjacent spatial regions or sequentially in the same region. Closely related, and perhaps even more frequently encountered, are fluid flow problems in which the flow speeds are supersonic in one direction and far subsonic in an orthogonal direction . For any of these and similar circumstances, it is necessary to have a single formulation of numerical techniques capable of spanning the entire range of Mach number, M, from zero (the #incompressible* limit) to much greater than unity (the hypersonic limit) . For this purpose we have developed a n implicit formulation of the equations that reduces to the MAC method for M =0 and becomes a variant of the FLIC technique for M % 1 .0 . The essence of the ICE method is that it isolates those parts of the equations related to density variations and treats them purely implicitly . While
F . H. Harlow and C . W . Hirt
128
most of our applications have been to multidimensional problems, ['0] it is convenient to illustrate the technique for a simple one-dimensional configuration, for which the differential equations are apu 2 ax
apu Op apu ax at +3t 89E+BauE='
a [(p+
a(P+q) ax
For the viscosity function, q, we choose for illustration au ax it is convenient to define,
x
R,+I = 6
[n++'1=(P+u+-'1
-P++n++ •, ) + qc - qi++],
so that the difference equations for mass and momentum become ,+1
n
ht
8 = Ox
[(Pu}+-'=-(Pn) ;+v ]~
1
+
8x0
[(pu)Ij+-(PU)+ w]
and (P(pu)+v, 6t ~
1-¢ „ 0 (P"-P,+t)+R++Y+ • _((pPi;-Pi++)~ 6x
The hybrid function, p, is formed from the equation-of-state pressure function, p(p, I), in a manner similar to that used by Janenko and Neuvazaev ["] for a somewhat different purpose . P< =
pi ,
+ ei n (P
,n+1-
Pi'),
where c ;" - (Op l dp) i" is related to the speed of sound. (With a constant specific heat, ci is the square of the isothermal sound speed .) Thus Pi singles out from the equation of state its principal dependence upon the density and for it uses the advanced-time density, while the rest remains at the beginning-of-cycle values . Without this procedure, i.e., with c;" = 0 and ; or with 0 = ¢ = 0, the procedure reduces to an explicit Eulerian treatment, in which pressure depends only on the previous-cycle values of p and 1 . Note the occurrence of weighting constants, 0 and ~, by which the technique can be varied from purely explicit (8 = ~ = 0) to partially time cen-
Recent
extensions to Eulerian methods
129
tered (0 = b = 0.5) to advanced-time implicit (0 = ,~ = 1 .0) . To solve( ' these equations, we eliminate the advanced-time mass fluxes, which ate the same as the momentum densities, and thereby derive an equation for the hybrid pressure function, p . To exhibit this, it is convenient to rather together all of the beginning-of-cycle quantities by defining the function G, as follows G; = P' C'"
+ bt dx
+ 0 (1- 0) 0t2
[(Pu)1n 9.-(Pu) ;+61+
n
Sx'
estz
6x
- An+vJ
in terms of which the equation to be solved is 46v (Pi-, + pi+, - 2P) -{- G ; . 6X' This, then, describes the essential features of the ICE technique . The steps for each computational cycle can be summarized as follows : 1 . Calculate the required beginning-of-cycle terms, R; + G, ; 2. Solve for pi for every cell. For some circumstances a direct solution method can be used ; for others it may be necessary to employ a relaxation technique ; 3. Solve for pi+t using the definition of p, ; 4. Solve for the new values of the velocities ; 5 . Calculate the new energy for every cell, using an appropriate finitedifference approximation to the energy equation . There are various wags that the energy and momentum convection fluxes can be written, the precise form being irrelevant to the conceptual formulation of the ICE technique, but nevertheless crucial to the assurance of accuracy and numerical stability of the calculations . At least four types of the flux expressions are appropriate for consideration, each having advantages and disadvantages . Consider, for example, the momentum flux term, (pa-), . While the density is already stored for the position at which the flux is to be evaluated (that is, at the center of cell i) the velocity positions are at i± t/2, so that some type of interpolation is needed . Four flux expressions that could be considered are Centered : (Pu 2 )i- pi[(ui-'h + ui+ ;.) /21'' . ZIP : (pa 2 )i -' Piui-I u,-'G~
F . H. Harlow and C. W . Hirt
1 30
Partial Donor: Ott)
ut+%) /2, { piui-'h (u,-k + p .u=+%(u--'i,+ui+S)/2 .
(U'_'/' + u,+ss) > 0, (u_v.+ut+v,)<0.
Complete Donor. )t (pa')
z .,'h,
(ui-y --u++'+ h) > 0, (u+-'ti -f' ut+'G) < 0 .
In addition there is an interpolated donor form, not considered here. For the energy-equation convection terms there are analogous flux expressions for two of these forms : Centered : '/s (p+ + p ;+a) u++v, (Et + Ei+, ), (PuE) i+ Donor : p;u;+9,Et, ut+a > 0, (P nE ) ;+u, ui+ { P ;+,u++wE++t,
< 0.
Some properties of each of these are the following : 1 . Centered - This form eliminates one order of Sx truncation error, compared with the donor cell flux, but this advantage is overshadowed by the tendency to numerical instability that is produced by not time centering the convective flux . This tendency to instability is usually mitigated by the addition of a counterbalancing diffusive flux (either real or artificial) and car) also be tolerated in circumstances in which the growth of the instability is bounded by non-linear dissipative effects . 2. ZIP-This form is also one order in 8x more eaccurate* than the donor cell flux, but in addition has the advantages of eliminating a nonlinear contribution to instability ['2] . For these reasons, the ZIP flux is often preferred for the momentum convection . The ZIP flux also needs a counterbalancing diffusive flux when it is not centered in time . 3 . Donor-Neither the partial nor the complete donor cell flux is space centered, and accordingly both have low order (in Ox) truncation errors . These contribute a positive diffusive effect and tend to automatically stabilize the numerical calculations . The magnitude, however, can be excessive, thereby masking the true diffusive effects and leading to erroneous interpretations (for example, for flows at high Reynolds numbers) . The once-popular donor-cell flux representation is being abandoned by most invertigators except for very special circumstances . Our recommendation for ICE calculations is to use ZIP fluxes or centered fluxes, the instability of these mitigated by the controlled use of artificial diffusion with a coefficient that just slightly exceeds the effective value pre-
Recent extensions to Eulerian methods
1 31
dieted by truncation error analysis . The only circumstance in which a donorcell flux appears desirable is in the early stages of problems with especially violent initiation (for example, in high-speed projectile impact), the later stages of the calculation automatically shifting to the centered or ZIP flux . To show some of the ICE -method properties, we may expand the difference equations in a form that reveals the lowest order diffusional truncation effects . With ZIP differencing of the momentum convection terms, for example, the result is Op
8pu
r)t
OX
8u l
au
6x2 0 5 2p 8t =1 (20-1)(c +u`) 2 4 dx
I ap
+a-+0x p (2¢-1)c~
dx 3uE ]
6t 2
usx= ap 8p
Ox
~.
+ p 1
a° u OX
We have used c ° to denote either of several partial derivatives of p with respect to p ; the symbol as used in this section is only approximately related to the square of the sound speed, being progressively more so in the limit as C'- W .
Several classes of effects can be noted from these equations . For r2 -. 00 (the incompressible limit), the effective diffusion coefficient becomes large and negative if 0 and/or ~ is less than 0 .5, corresponding to intolerable mimerical instability . In contrast, for 0 and ¢ both greater than 0 .5 the diffusion positive (hence stabilizing), and the equations suggest that this could become intolerable as cz -> 00 . This excessive diffusion, however, is purely longitudinal, so that two-dimensional incompressible flow calculations with the ICE method ordinarily will not suffer from the effect, since longitudinal signal propagation is extremely rapid, anyway . Despite the implicit nature of the ICE method, there are still some restrictions that must be observed in order to assure numerical stability . These are difficult to describe precisely ; the matter is discussed more fully in ['°] . § 6 . Turbulence models in numerican fluid flow calculations There is considerable reason for believing that the most important effects of turbulence depend only on some of the collective effects of the fluctuations in a way that is independent of the detailed structure of the turbulent flow . Because of the widespread occurrence and profound effects of turbulence there has been a strong effort to discover how to represent these col-
F . H. Harlow and C. W . Hirt
132
lective effects in a useful way for predicting turbulence and for analyzing its numerous important effects . One of the results of these investigations has been the relatively new idea of the transport-equation model . Several versions of this have been proposed and proof tested at Los Alamos, with results that are highly encouraging . These results, together with a bibliography of earlier work are given in references P ." "] To illustrate the central problems involved in these investigations, consider the following equation
0R,;
OR„; 0k
-- -1- U,
ox,
(it
1a,(' - ( ax, +
On,
R*A+ ax,,
au,'¢ \'
on, 1j
X,
aa ;
r3x, I +
8 dx,
Rla-t-(ui'uJ au ;
( dx ; + 3 %x,
) +
z
" - 4D;,} -
v ~-
T"
ax„'
(g;ui T' + g,ui T') .
This equation for the Reynolds stress tensor, R,;, is a rigorous consequence of the Nr avier- Stokes equations for an incompressible fluid, To derive it, the velocity of a fluid element, U, is considered to be the sum of a mean-flow part . a,, and a fluctuating part, a,' . Likewise, the pressure, ~, and the temperature, T, are split into the mean part (designated in multi-factor terms by an overbar) and a fluctuating part (designated by a prime) . With a formulation normalized to unit density, the Reynolds stress is R;; = u ;'a ;' . The other quantities in the equation are g„ the gravitational acceleration ; To, the reciprocal bulk expansion coefficient in the Boussinesq approximation to the buoyancy term ; v, the molecular kinematic viscosity coefficient ; and
1
a u;' 8u ;
2
ax„
ax,
which describes the decay rate of Reynolds stress . It is likewise possible to derive a rigorous transport equation for D, ; from the Navier - Stokes equation ['"] . Although this is a fundamental part of the theory, it is omitted here, not being necessary for illustration of the principles that are, to be demonstrated in the present paper . In the form presented above, the Reynolds stress equation cannot be solved because of the presence of various moment terms (the overbar terms and D ;) that are not related to the available variables . The central problem in turbulence modeling, therefore, in to discover appropriate ways for representing these terms as functionals of a set of variables that is no more numerous than
Recent extensions to Eulerian methods
133
the set of equations that can be found to desci ibe their variations . One form for such an equation, derived in["] . introduces a scale function, s, together with various universal, dimensionless constants and functions . a, 0, r, to, 0, 4 in such a way that oR i;
du;
oR,
`+UAW-
fdl
+
ox,
UX,f
1~ ~ ;~,
d 8
d'R,, ~ , r)X,
[
du, ox,
`R,A--'-"R ;A-+
s-q
\'
OR,A
r)x,',vi\
2v4 s-
dx.A / 'US'
-R1
tj
dx,
-(g,R,A
(3TA .N
d s R, 'R„ a-v,
-J 11,
I
ox,
q o R,A A A' .A
dx .~,
/1
oT ± g;A ; ;J-,
where ~3 is the thermmuetric condur-tivit.v, and (D„ = - ~,~ (R„ - 1 ,i s' aU
RA
,b n l -
UIE ; }
I oU A
x„ ,
R ' `a,r,
tz~R,4~dx ;
du r t~ x,.~~
if the scale function is specified or calculated, then this tensor equation for the Reynolds stress is complete, and its simultaneous solution with the mean-flow equations determines the flow . To specify the scale function, ona could resort to simple mixing-length theories, or else utilize a second transport equation . Reference [ "] discusses this latter, and shows how this can be accomplished by means of moment approximations in the transport equation for D,,. Several properties of a flow must he present in order for a transport-equation model to be valid . In a division of the flow into mean and fluctuating parts, the spectrum of the fluctuations must not strongly deviate from some universal, self-similar, equilibrium spectrum . There simply are not enough degrees of freedom in a low-order moment representation to describe the effects of a complicated spectrum that varies by large amounts from case to case . Indeed, there must be some cases in which strong departure from equilibrium removes all validity from the low-order moment equations . At present, however, we have very little information as to how far from equilibrium a turbulent flow can be and still be capable of description by a transport model, For many kinds of investigations it will be necessary to treat both the turbulence and the mean flow as time varying, the whole being so complicated that to get solutions will require numerical methods with high speed computers . At present, there are investigations in progress to determine just how extensively the turbulence transport models by themselves can he applied to flows with superimposed large-scale fluctuations, so that the matter is still far from settled . What is clear is that for truly one scale turbulence problems, the
1 34
F . H . Harlow and C. W . Hirt
present models for Reynolds stress transport are capable of producing excellent agreement with experiments . Whenever there is a need for simultaneous coupling of the mean-flow equations to the turbulence transport calculations, we may expect one of the high-Reynolds-number difficulties usually associated with finite-difference calculations to be mitigated by the diffusive properties of the effective eddy viscosity . Numerical stability considerations usually restrict the Reynolds number to magnitudes less than the square of the number of computational cells required to resolve a typical linear dimension of the flow field . For laminar flows at high Reynolds number, it is accordingly necessary to introduce an artificial viscosity (either explicitly or else as a consequence of the differencing technique) whereby the effective Reynolds number is reduced to the required magnitude . In many cases the results are nevertheless satisfactory because of insensitivity of the laminar flow field to Reynolds number variations when the Reynolds number is sufficiently large . For fully turbulent flows, however, this decrease in effective Reynolds number, resulting from the eddy viscosity can be sufficient to preclude numerical instability, even when the resolution is sufficiently modest for practical computation .
§ 7 . Curved surfaces In Eulerian fluid-dynamics calculations, the occurrence of curved boundaries requires special treatment . A variety of examples come to mind : 1. Free surfaces ["] . 2. Moveable interfaces between materials ["] . 3. Curved rigid walls 4. Moveable walls ["] . In every case, the principal difficulty arises from the partial cells created by the boundary . For free surfaces, it is necessary to balance both the normal and tangential stresses . Experience shows that the original version of the MAC method ['] is not sufficiently accurate in resolution of the surface configuration for problems with large viscosity and / or low amplitude waves . To alleviate the difficulty, it has been found useful to introduce surface marker particles, and to apply the stress conditions at the actual surface position rather than at the cell center ["] . For moveable interfaces between materials, a similar marker-particle set has also been useful as a means for determining how to better resolve the continuity of normal and tangential stresses ["] . Curved rigid walls present a somewhat different challenge that has been solved by Viecelli ["] in an intriguing manner . He introduces a pressure boundary condition that accelerates the fluid in exactly that manner that is re-
Recent extensions to Eulerian methods
135
quired to keep the fluid velocity tangential to the wall . For moveable walls, he uses a similar technique [f 8 ], in which the normal velocity is adjusted to the prescribed magnitude, rather than vanishing . A somewhat different approach to the problem of curved surfaces (fixed or moving) has been taken by Hirt and Amsden [10 ] . They have devised an Arbitary-Lagrangian-Eulerian (ALE) technique, in which the computer mesh can vary from completely fixed (but not necessarily with rectangular cells) to moveable with the fluid velocity . With this, they can resolve curved and /or moving wall configurations with the Lagrangian capability, and at tho same time resolve the strong internal distortions of the fluid with the Eulerian capability .
Fig. 1
1 36
F . H. Harlow and C . W. Hirt
Fig . 2
§ 8. Theree space dimensions Time-varying problems in three space dimensions are easily solved with Eulerian techniques, especially when the computer program is based on the efficient new iteration method of Hirt [°] . Recent calculations by Hirt and Cook ['°] have shown both efficiency and accuracy in the investigation of at • mospheric wind flows over variable terrain and around buildings . Extensions that include heat transport, buoyancy effects and smoke transport have demonstrated a capability for the realistic analysis of atmospheric pollution studies, and the techniques are readily extendable to underwater pollution and sediment transport and deposition investigations . § 9 . Sample calculations A variety of sample calculations are described here to illustrate some of the techniques described in the preceding sections .
S 2TJ
N
I PV
£ t?3
al
spoil aw
U
Aajng o7 suo]suapxa lua ag
138
F . H . Harlow and C . W. Hirt
Fig . 6
Fig . 1 illustrates the evolution of a block of viscous material such as far, resting on a table and sagging tinder the influence of gravity . This calculation was performed using the technique described in § 3 and ['] . These
Recent extensions to Eulerion methods
139
pictures show the configuration of the marker particles used in the calculation to delineate regions occupied by fluid . The first frame a is of the initial block shape while
b
through
d
show subsequent stages in the time
history of the block . Fig . 2 illustrates a type of problem that can be studied in detail with the improved Marker-And-Cell techniques mentioned in § 7 . The example illustrates what happens when a spherical drop of fluid falls into a container of fluid . The picture sequence follows the splash from the moment of initial impact a, through the point of maximum cavity formation c, and into the secondary splash resulting from cavity collapse
f.
Another new technique mentioned in Section 7 is the Arbitrary-Lagrangian-Eulerian method . The unique feature of this method is its ability to resolvo flexible boundary walls . This capability has been utilized to calculate the pulsating flow of blood through an elastic artery . Fig . 3 presents three kinds of data produced directly by the computer . In a the calculational mesh for a cross-section of the tube is shown with a two pulse disturbance visible on its right edge . The left edge of the mesh is an axis of symmetry . The mesh is continuously rezoned to maintain the radial grid lines in fixed positions, and the axial grid lines were required to intersect the radial lines with equal spacing between intersection points . This choice of rezoning maintains a uniform mesh pattern, which is necessary for accurate results . In b and c are shown velocity vectors and contour ones of constant pressure, both at the same lime as the grid configuration shown in a . It will be observed that maximum pressures occur under the surface bulges and minimum pressures under the surface indentations . corresponding to the application of elastic stresses in the confining tube . A variety of three-dimensional calculations have been carried out with a new version of the Marker-And-Cell technique . Fig . 4-6 illustrate the results of several calculations involving the incompressible flow of air about large rectangular structures . In Fig . 4 a uniform flow is entering the computational region (large rectangular box) front the left silo and, after flowing around a building (smaller rectangular box), is allowed to flow out the right side of the computational region . In a the velocity vectors in a layer of cells passing over the top of the building are shown in perspective, while b shows the velocities in a layer of cells near the bottom of the mesh . A large stagnation region behind the building is clearly visible in b . Fig . 5 shows the results of a similar calculation performed with a second, smaller, structure located behind the first one . Finall, Fig . B illustrates a distribution of particulate matter (air pollution) calculated, for the One and two building problems . In a there is a strong concentration of pollution in the, wake region behind the building,
140
F . H . Harlow and C. W . Hirt
while in b no such concentrated region exists, because the smaller structure deflects the flow across the near face of the large building and prevents the formation of a large recirculation region .
Translated by D . E . Brown
REFERENCES 1.
HARLOW, F . H. Numerical methods for fluid dynamics, an annotated bibliography . Los Alamos Scient. Lab. Rept. LA-4281, 1969 . GENTRY, R . A ., MARTIN, R . E . and DALY . B . J . An Eulerian differencing method for unsteady compressible flow problems . J . Comput . Phys ., 1, 87, 1966 .
3.
HARLOW, F . H . and WELCH, J . E . Numerical calculation of time-dependent viscous incompressible flow of luid with free surface, Phys . Fluids, 8, 2182 . 1965 ; WELCH, J . E ., HARLOW, F . H ., SHANNON, J . P . and DALY, B . J . The MAC method . Los Alamos Scient . Lab . Rept . LA-3425, 1966 . CHORIN . J . Numerical solution of the Navier - Stokes equations . Math . Comput, 22, 745, 1968 .
5.
AMSDEN, A . A . and HARLOW, F . H . The SMAC method . Los Alamos Scient . Lab. Rept . LA-4370, 1970 .
6.
HIRT, C . W . A general algorithm for the solution of implicit difference equations . Los Alamos Scient . Lab ., Group T-3, unpublished notes, 1971 .
7.
PRACHT, W . E . Implicit solution of creeping flows, with application to continental drift. Proc . See. Internat . Conference on Numer. Methods in Fluid Dynamics, Univ . California, Berkeley, California, Sept . 15-19, 1970 ; PRACHT, W . E . A numerical method for calculating transiet creep flows, J . Comput. Phys ., 7, 46, 1971 .
8.
PRACHT, W . E . and HARLOW, F . H . A numerical solution technique for the investigation of planetary interior dynamics . Los Alamos Scient. Lab ., Group T-3 manuscript, 1971 ; HARLOW, F . H . and PRACHT, W . E ., A theoretical study of geothermal energy extraction, J . Geophys . Res ., submitted for publication .
9.
GENTRY, R . A . and NAKAYAMA, P . 1 . Loa Alamos Scient Lab., unpublished manuscript, 1971 .
10 .
HARLOW, F . H . and AMSDEN, A . A ., A numerical fluid dynamics calculations method for all flow speeds, J . Comput . Phys ., 8, 197, 1971 .
11 .
YANENKO, N . N . and NEUVAZHAEV, V . E . A method of computing gas-dynamical movements involving non-linear heatconductions,Tr . Matem . in-ta Akad . Nauk
12 . HIRT, C . W . Heuristic stability theory for finite difference equations, J . Comput, Phys ., 2, 339, 1968 .
Recent extensions to Eulerian methods
141
13 .
DALY, B . J . and HARLOW, F . H . Transport equations in turbulence . Phys . Fluids, 13, 2634, 1970 .
14-
HARLOW, F . H . Turbulence models . HARLOW, F . H . and AMSDEN, A . A . Fluid dinamies, an introductory text . Sec . ed . Univ. California, Los Alamos Scient . Lab, Monograph, LA-4700, 1971 .
15 . CHAN, R . K .C, STREET, R . L. and STRELKOFF, T . Computer studies of finiteamplitude water waves . Stanford Univ . Techn . Rept. 104, 1969 ; CHAN, R . KC and STREET, R . L . A computer study of finite amplitude waves, J . Comput . Phys ., 6, 68, 1970 ; NICHOLS, B . D . and HIRT, C . W . Improved free surface boundary conditions for numerical incompressible-flom calculations, J . Comput . Phys ., 8 . 434, 1971 . 16 .
DALY, B . J . Phys . Fluids . 10, 297, 1967 ; DALY, B. J . and PRACHT, W . E . Numerical study of density-current surges, Plevs . Fluids, 11, 1, 1968 ; DALY, H . J . A technique for including surface tension effects in hydrodynamic calculations, J. Comput . Phys ., 1, 4, 1969 ; DALY, B . J, Numerical study of the effect of surface tension on interface instability, Phys . Fluids, 12, 1, 1969 .
17 .
VIECELLI, J . A . A method for including arbitrary external boundaries inn the MAC incompressible fluid computing technique, J . Comput . Phys ., 4, 543, 1969 .
18 . VIECELLI, J . A . A computing method for incompressible flows bounded by moving walls, ;J . Comput . Phvs ., 8, 119, 1971 . 19 .
HIRT, C . W . and AMSDEN, A . A ., An arbitrary Lagrangian _ Eulerian computing technique for all flow speeds, J . Comput Phys ., to be published.
20 . HIRT, C . W . and COOK, J . L. The calculation of three-dimensional flows around structures and over rough terrain, J . Comput . Phys ., accepted for publication .