|
Items in cart [0] |
Questions? Call 888-624-8373 |
| Copyright © 2009. National Academy of Sciences. All rights reserved. Terms of Use and Privacy Statement |
Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 269
Pressure Transients in Transitional Boundary Layer
over a Solid Surface
Sin-T Cheng
Princeton University
Princeton, USA
I. INTRODUCTION
The Navier-Stokes system of equations of
an incompressible fluid, upon neglecting the
gravity and other extraneous forces, stands
as:
Continuity: ~ =
auf = o
axi
(1)
UP ~ .. ~ ui
Momentum: an] + Uj axj = - p axi T V ~X; ~X;
(2)
The continuity equation (1) is a limiting form
of
Bu; 1 dp 1 ~ ~
dXi ~ p dt = ~ p ( at + uj axj )P
(3)
for a compressible fluid, when dp/dt vanishes.
For a compressible fluid with a barotropic
relation p(p), equation (3) stands as
( ~3 + uj ~xa )P = PC2 ( anti ) (3a)
where c2 = dp/dp is the square of the speed of
sound in the medium for the specified
barotropic thermodynamic process, which speed
of sound becomes infinitely large for
incompressible fluids. Equations (2) and (3a)
define an initial value problem with the
specific choice of p(p) and the intial dataa
for the solution of the four scaler unknowns p
and ui, i = 1, 2, 3. Equations (1) and (2) is
a singular limit of (2) and (3a) when
d/dt p vanishes in the incompressibility
limit. There is no explicit means to advance
pressure in time.
269
In classical hydrodynamics(l), a velocity
potential ~ with Ui = vie = B/3xi ~ is
introduced for an irrotational (w = v x u = 02
flow, the continuity equation is reduced to v'
= 0 which can be solved without any
reference to pressure under suitable boundary
formulation. For steady state potential
flows, 2the Bernoullis relation
p + pu /2 ~ B as an integral of the momentum
equation enables us to find p where B can be
evaluated from the boundary data. For time
dependent flows, Be B¢/3t depends on both the
initial and the boundary data so that pressure
transients are not so easily determined from
the Bernoulli's integral even for potential
flows. The pressure p(t,x) and the velocity
u(t,x) should be solved simultaneously from
equations (1) and (2).
In the absence of a Bp/3t, (1) serves as
a subsidiary condition that the correct p(t,x)
must evolve simulataneously with u(t,x) so as
to keep the velocity field solenoidal at all
times and everywhere. In the computational
solution of some discrete form of (1) and (2),
we can only solve them iteratively. We hope
that the successive iterations would yield
better approximations to the pressure and the
solenoidal velocity fields and that the
sequence might eventually "converge" to the
true solution. Given a solenoidal initial data
ui(o) at t ~ to and the corresponding initial
pressure field p(o). We wish to find ui(l)
(to + fit) that is solenoidal at to + fit and
its corresponding P(~)(to + St). We do not
know the evolving pressure to advance Ui from
(2) throughout the interval St. Hence we take
it to be the known p(o) = p(to, x) or some
otherwise conceived pressure field P(o,1) for
evaluating the gradient Bp/3xi throughout x
and during the interval fit. Equations (2)
then give the estimate Ui(O,l) of Ui(l) This
timate of ui(o~l) will possess some
residual divergence A(o,1) over the field. We
wish to devise some format to estimate the
correction to P(o 1) such that P(0,2) = P(0,1)
+ [P(0 1) that will generate a new estimate of
Ui(O 2) with a smaller residual divergence
A(O 2); and its repeated application yielding
A(O,n) ~ O almost everywhere at some
sufficiently large n.
OCR for page 270
By taking the divergence of equation (2),
we obtain
( ad + Uj as - 2' aa2a ) ~ + ~ axe
1 32
=
p ~Xj ~Xj
(4)
which should stand as the evolution equation
for the divergence ~ of the velocity field.
Equation (4) is, however, often interpreted as
the Poisson equation defining the pressure
p(t,xi). In particular, for an incompressible
fluid with ~ = Bui/6xi = 0, we have
2 Bu Bu
v p = - P axj ~
_ = p (eij eii - wij Id) (4a)
so that the pressure field p(t,x) can be
determined from the solenoidal velocity field
ui(t,x) at any time t. It is often that
P(O 2) is taken as some discrete solution of
(4a) as a Poisson Equation with the source
terms estimated from the perturbed velocity
field Ui(O 1) or some weighted averages of
ui(o,O) and Ui(O,l). Alternatively, P(O 2)
may be taken as some weighted average of the
solution of this Poisson equation with P(o 1)
for further iteration. There are numerous
such alternatives. They give rise to a
pressure iterate P(o,2) with significant
spurious oscillations that quickly become
chaotic and apparently diverging with repeated
applications. The use of smaller fit fails to
help. The spurious oscillatory field has to
be filtered and smoothed to keep the results
bounded.(2) Filtering and smoothing has been
very popular and built into many robust
computational codes. The nature of such
filtering functions is obscure and its
desirable form highly problem-dependent.
Clearly, with A(o,1) = 8/3Xi Ui(O,l)
O. equation (4) should be used instead of
(4a). To solve (4) as a Poisson equation, we
need the estimate of B/8t ~ which may be
directly estimated from Ago 1) in a variety
of ways, such as +~(o,1)/6t, zero or -
A(o 1)/6t depending on how (4) is supposed to
achieve physically in the iterative
procedure. ) Such iterative solutions of p
from the Poisson Equation (4) also failed to
converge to some meaningful ui(l) and P(1) at
to + it, but led to chaos. Such chaotic
fields are numerical turbulence generated by
the iterative procedure for solving the
nonlinear system. They can hardly be
suppressed or damped by any reasonable amount
of artificial viscosity that might be
introduced into the computational algorithm.
Such computational difficulty of securing
detailed mass balance by varying pressure
prompts the development of computational
methods that avoid the determination of
pressure. Since the curl of the momentum
equations (2) eliminates the pressure
explicitly to give the vorticity transport
equations:
a ~ = v x ( u x a) + v v2 ~
Bt ~ ~ ~ ~
where
and
(5)
~ = v x u
~ ~ (6)
2 32
v u = axj~xj Ui = -v x ~ (7)
the vorticity ~ may be advanced by (5) and
the velocity u solved from (6) or (7),
apparently without any explicit reference to
pressure. This approach has been quite
successful for flows in two space dimensions
where the continuity relation can always be
satisfied with a stream function, and the
vorticity vector with only one nonzero
component vale is always solenoidal. The
vector equation (5) is simplified to a scalar
equation of vorticity transport. This stream
function ¢, and hence the velocity-vorticity
field, can be solved computationally.
Presumably the pressure p(t,x) can be obtained
by integrating the momentum equation
spatially. The pressure so evaluated at a
given point turns out to depend on different
paths of integration even if ~(t,x) has
converged to some well defined steady
state.(3) Thus the pressure transients cannot
be determined with such stream function
vorticity formulation even for flows in two
space dimensions.
For flows in three space dimensions, the
flow field evolution is generally visualized
according to equation (5) as the processes of
convection u vow interaction ~ vu and
dissipation u~ w. The effect of vorticity
interaction is often evaluated through the
Biot-Savart Law as a consequence of equation
(7). Such numerical solutions have often been
referred to as the solution of the Navier-
Stokes system(4). Without ingenious smoothing
and filtering, such computational solutions
quickly lead, however, to chaotic oscillations
even at Reynolds numbers significantly below
the critical Reynolds numbers of laminar flow
instability. Such chaotic oscillations are
clearly of numerical origin (or numerical
turbulences. The question of associated
pressure field is rarely mentioned. Why
should the computational solution through
vorticity formulations for the three-
dimensional flows be so much more difficult
than that of the two-dimensional flows?
270
OCR for page 271
It is simple to verify that a planar
field will rapidly become three-dimensional
and nonsolenoidal under its own convective
motion and their mutual vorticity interaction
as is evaluated by the Biot-Savart Law. A
nonsolenoidal vorticity field is
mathematically incompatible with its
definition equation (6). A nonsolenoidal
velocity field violates the continuity
relation (1). The derivation of equations (5)
from (2) implies solenoidal vorticity and the
derivation of equations (7) from (6) implies
solenoidal velocity. The popular approach of
considering vorticity transport and
interaction in discrete form without due
attention of retaining solenoidal velocity and
vorticity fields is fundamentally
questionable. It is often tolerated in the
construction of approximate solutions if the
residual divergences could be controlled to
remain small, and the global mass
conservation can be monitored and likewise
maintained. Thus, many forms of Mark and
Cell (MAC) and the Particle in Cell (PIC)
methods were developed as early as 1950's for
integrating (1) and (2) by introducing
Lagrangian particles to track the mass balance
of many sub-ensembles represented by such
particles. Alternatively, the Navier-Stokes
equations may be integrated in some mixed
Lagrangian-Eulerian formulations. The
calculated flow field can appear quite
reasonable and even be made to reproduce
selected details of global experimental data;
but adequate pressure field is yet to be
recovered from the computed velocity field
that is not quite solenoidal.
For flows in three space dimensions the
velocity vector can be represented as u = vat +
vxA where the scalar potential ~ is defined
through v2¢ = 0 and the vector potential A is
related to vorticity as ~ = v x v x A.. The
pressure and the velocity transients are
clearly not solely determined by the evolution
of vorticity ~ alone. For the solution of ~
and ~ (or v x A), the nonslip condition u = 0
is split into normal and tangential
components and applied separately. There is
no guarantee that the sum of v} and v x A will
satisfy either. It is not surprising that
computational solutions of such stream
function-vorticity formulation can hardly be
keep bounded without constant and repeated
"smoothing" and "filtering". While it is
possible to devise ingenuous "filtering"
schemes to reproduce some preconceived
solution, we hardly understand their physical
basis. Therefore we attempt to analyze the
situation and then develop a computational
algorithm free from such uncertainties.
II. Analysis of Artificial Compressibility
Approaches
The solution of (5) and (7) with
nontrivial boundary is not really independent
of pressure. On a solid boundary at rest, for
example, the nonslip conditions require
~ v x ~ = ~ Up (8)
under the solenoidal velocity and vorticity
fields. Any approximate formulation of the
vorticity boundary data in the solution of
(5), such as vex = 82/8xj~xj Ui evaluated
with the previous iterate, is actually some
approximation to the spatial gradient of the
pressure field on the solid surface. With the
appropriate pressure boundary data on the
solid surface and those at infinity (or other
closure surface), equation (4) or (4a) will
define the corresponding approximate pressure
fields. The iterative solution of the
vorticity transport equation (5) with
successive approximate boundary formulations
on the nonslip surface is equivalent to an
iterative process with the corresponding
sequence of approximate pressure field.
Accordingly, the iterative solution of (5) and
(7) through the intermediate variable
vorticity is simply an alternative form of the
pressure iterative solution of (1) and (2) in
terms of the primitive variables p and ui. It
suffices, therefore, to study the iterative
processes for the solution of (1) and (2) in
primitive variables, ui(t,xj) and p(t,xj), to
identify some property conducive to a
converging iterative process.
We shall show that the iterative
processes described in the previous sections
attempting to reduce the residual divergences
cannot converge for nontrivial initial
boundary value problems in three space
dimensions over a drag body (or surface).
For convenience, we shall refer to these
methods as "Artificial Compressibility" in
view of equation (3) that any residual
divergence ~ = Bui/3xi represents the
"condensation" or the rate of fractional
variation of fluid density, i. e.,
"compressibility. " The sequence of iterants
may be interpreted as describing the
succession of "compressible states" which, in
the limit of small residual divergence and
compressibility, was hoped to converge to the
limit of incompressible flow. (5) (6)
We describe the success ion of such
artificial compressible flows by the discrete
equivalent of the following set of partial
differential equations:
~ at=~~=~ axi (9)
271
OCR for page 272
Hi + uj aXj
- p ( axi+ Fi) + ~ a,.jaaX ui (10)
Here ~ is the function that describes how
~(O 1) ~ p(O,2) p(O,l) is evaluated from
O,1) in the iterative correction scheme to
the time step fit.
ran the ~PtAi 1 ~ of the
achieve ^(O,n) ~ 0 for
It clearly depends also ~
scheme of discretizing (9) and (10) for
computational solution.
For a physical medium with lop/8p ~ pc2 >
O and ~ p-1 [p/St, ~ is always positive and
vanishes as ~ ~ O in a physical limit process
of approaching an "incompressible" state
through successive "artificial compressible
states" of decreasing compressibility. A
residual divergence ~ in an incompressible
fluid is a volumetric source of the flow or
the mass source divided by fluid density.
This mass source will carry with it some
momentum which should be represented by some
external force Fi in the momentum equation
(10). Our attempt to construct a converging
computational algorithm for the solution of
(1) and (2) is thus recast into finding the
four scalar functions ~ and Fi such that a
converging solution of (9) and (10) through
some discrete equivalent system would exist
(and hopefully unique as a well posed problem)
in the limit of ~ ~ O at all times. It is
therefore necessary that:
km ~ (I) = 0
O
km Fi (I) = 0
O
(11)
so that (9) and (10) will reduce to (1) and
(2) respectively. The four functions c, Fi
are otherwise completely arbitrary, and may be
distinct for different discretization details
through finite difference, finite element,
spectral, or any other method. We attempt to
identify first the properties of ~ and Fi
that would facilitate, if not secure,
convergence.
By multiplying p into (9) and put into
(10) and integrating their sum over the volume
V of computation, the following temporal
variation of an energy integral E, is obtained
where the Stokes theorem has been used to
convert some volume integrals into surface
integrals over the bounding surface S of V.
Bt at J ( 2 u + 2 p2) dV
= ~~ TV ( axj dxj ) dv
+ ~ ui ~ P ui ~ -Fi ~ dv
+ ~ 1 am ( P 2 ) dS i + ~ - ( p + P 2 )
(12) 272
With ~ non-negative the quantity E = p u2/2 +
~ p2/2 can serve as the measure or the
magnitude of the set of computed results,
converging (or diverging) if BE/8t
O) at all later times.
The first volume integral on the right hand
side of (12) represents viscous dissipation
which is always negative and stabilizing. The
last two surface integrals depend exclusively
on the formulation of the computational
boundary. Then for a given boundary
formulation:
(i) The convergence behavior of an iterative
scheme depends only on the choices of Fi
(t,xj), i.e. the momenta carried by the mass
associated with the residual divergence A, not
the details ~ how the pressure correction dp
is evaluated from A. In other words, it is
futile to devise different iterative
algorithms ~ of advancing pressure from
residual divergence as is described in the
previous section, to promote convergence
O.
(ii) The boundary formulation that determines
the surface integrals on the left hand side of
(12) is important. Computational methods
highly successful for periodic boundary value
problems need not be as successful for
nontrivial problems as will be explained
below.
For periodic boundary values problems,
both the two surface integrals in (12) vanish.
The first volume integral is proportional to
viscosity and always negative. Therefore, if
we choose Fi = pui^/2 everywhere at every
iterative steps to render the second volume
integral zero, the iterative scheme will
always converge regardless of how the pressure
bp is evaluated from A. The convergence rate,
in terms of the weighted L2 or RMS norm
expressed as BE/8t globally, is proportional
to the fluid viscosity (or the dissipation).
The standard L2 norm convergence rate of
pressure corresponding to a converginig
velocity field will, however, vary as e-1/2
which can be troublesome as ~ ~ O.
Actual computational solutions of simple
periodic boundary value problems of (9) and
(10) with some preassigned small values of
and Fi =pui^/2 for different simple test
problems have rendered highly successful
approximations to the corresponding
incompressible flow velocity fields (5,6)
although without as satisfactory pressure
results. The schemes of discretization and
various other details apparently do not
matter. Such solutions of periodic, boundary
value problems can be much improved with
smaller ~ and through Richardson's
Extrapolation etc. The evaluation of the
pressure field is, however, more difficult as
~ ~ O as is suggested by the theoretical
result given above. Computational methods,
have often been developed with periodic
boundary value and applied to nontrivial
boundary value problems. The results have
OCR for page 273
been mixed. Various filtering and smoothing
schemes can be built into robust codes for
computing the velocity fields. They are not
solenoidal and cannot provide meaningful
solutions of the pressure field.
In the presence of a drag (or lift) body
(or surface) in the flow field, the boundary
formulation cannot be periodic. The last
surface integral in equation (12) representing
the net outflux of the stagnation pressure p +
pu2/2 over the external boundary, will remain
positive and non-zero even in the limit of
O. With both viscous contributions small at
sufficiently large Reynolds numbers and with
|~ Fib) = 0, equation, (12) showns that
BE/3t will eventually be dominated by the drag
contribution, (i.e. the positive surface
integral Ed = J - (P + 2 ) Ui dsi )
and remain positive as A-O. Thus BE/3t will
remain positive at sufficiently small ~ and
E will eventually diverge as A,O. Then these
iterative scheme of reducing the residual
divergence cannot converge.
It is of course possible to choose Fi
sufficiently large all the times so as to have
BE/3t < 0, and hence, an apparently converging
iterative process. This converged solution is,
however, not that of (1) and (2) since (10)
remains different from (2) with large Fi ~ O
We do not know how such a converged solution
might render a satisfactory approximation to
our problem. In any case we never compute
till ADO and have to stop at some finite ~ or
e. We could presume that Fi would, from there
on decrease toward the set of values pui^/2 or
the like such that the E would decrease to
some Emin before diverging under the influence
of Ed. Such an Emin could suggest some
asympototic approximation to the solution of
our problem.
There can be a variety of choices of Fi
including those ingeneous forms of smoothing
and filtering. Such choices may render good
approximations to the velocity field
especially if we have some preconceived ideas.
For such an approximate velocity field, the
asymptotic errors are concentrated in the
pressure field. In terms of the standard L2
norm, the r.m.s. error of the pressure field
will grow as c-l/2 which can be very large if
the velocity field should be nearly
solenoidal (e << l). Thus the pressure
solution has to be postponed. This is
reminiscent of the situation of the
computational solutions described in section
II. Our interest in the pressure transients,
however, are not well served by such methods
of artificial compressibility including those
with suitable filter(s).
273
A reasonable iterative solution for the
pressure transients should be constructed from
velocity iterants that are solenoidal
everywhere and at all times at least after
some fixed time to, i.e. iterations in
solenoidal subspace. The pressure field
associated with each solenoidal velocity
iterate can be consistently determined from
(4a) unique up to a constant. The evolving
solenoidal velocity field and its associated
pressure field will have to satisfy (2) or its
equivalent discrete set.
III. An Iterative Solution Method
in Solenoidal Subspace
The evolution of a solenoidal velocity
field ui(O) (to) under some prescribed smooth
pressure field p(°)(to) according to equation
(2) will not necessarily produce a solenoidal
velocity field. Indeed the ~ = Bui/3xi will
evolve according to equation (4). There is no
reason to expect the prescribed pressure field
p(o)(to) and ui(°) to satisfy 4(a) so as to
warrant the natural evolution into a
solenodial velocity field. In any case, the
flow system could be subjected to arbitrary
disturbance of pressure and/or velocity that
violates (4a) and generate residual divergence
in the course of time. Therefore we have to
face the question what is the physical meaning
of such a residual divergence although not
physically admissible to an incompressible
fluid.
H. Lamb(l) introduced the notion that a
residual divergence in an incompressible fluid
is "equivalent" to a set of associated
"impulsive forces". In otherwords, an equal
and opposite set of associated impulsive
forces would annihilate the local residual
divergence to leave a solenoidal velocity
field, being disturbed by the set of
associated impulsive forces. This concept has
been little developed; but gives equations (9)
and (10) a physical meaning different from the
artifical compressibility extension of the
equations system (l) and (2) for an
incompressible fluid. The residual divergence
~ is accompanied by an equal but opposite pair
of impulsive forces -Fi and Fi introduced or
developed at any point and any time. The -Fi
will remove the residual divergence ~ to
restore a solenoidal flow field. The Fi
remains as the "equivalent" set of impulsive
forces acting on the solenoidal flow. Thus
(9) can be replaced by (1) with ~ = 0; and
(10) by (lo) with Fi disturbances equivalent
to the nonsolenoidal velocity disturbance A.
For a given residual divergence A, this set of
equivalent impulsive forces Fi is not unique
and not in "dynamic equilibrium". An
impulsive force is a point source of
"discontinuous" wave, propagating into the
flow field and modifying the local velocity
and pressure as it passes. The wave will
reflect from boundaries and interact with
OCR for page 274
waves trom other sources to relax the
disturbed flow field to establish a new
"smooth" solenoidal velocity field with an
equilibrated smooth pressure field consistent
with (4a) and the physical boundary data.
This relaxation process is very fast in an
incompressible fluid where the wave speed (or
the speed of "sound") is "infinitely" large.
The dissipative forces of friction and heat
transfer have no time to act, so that this
relaxation process is essentially "potential".
This relaxed solenoidal disturbance velocity
field can be represented uniquely by the
gradient of a velocity potential ~ as vie
u(l) (to + fit) = u(l) (to + fit) + vat (1
3)
This potential function ~ is given by
V2¢ = -div u(l) = - ~ (14)
and can be solved with the appropriate
Dirichelet and/or Neumann type boundary
formulation over the field without altering
the physical boundary conditions of u(l). In
terms of the relaxed solenoidal velocity u(l)
from (13), we can solve for the associated
pressure field p(l) It + St) from (4a). This
p(l) (to + St) and uLl'(to + fit) is the first
physical iterate at the first time step to +
fit. It need not be consistent with the
discrete form of (2) in stepping from to to to
+ St. As such, they have to be iteratively
corrected. We designate this first iterate
for the first time step p(l,l) and u(l~l) etc
(i) If there is an external non-solenoidal
velocity disturbance, take the non-solenoidal
disturbed velocity field as u(l l). If there
is no external disturbance, advance the
initial solenoidal velocity u(O) to u(l~l) by
evaluating Bp/3x from the initial data p(°
over the time step St.
(ii) Obtain ¢(1,') according to (14) and
u(l~l) from (13)
(iii) Obtain p(l,l) according to (4a) and then
choose some average pressure
p(l,l)/= p(l~l)/2 + p(O)/2 or the like, for
the improved evaluation of the spatial
pressure gradient in (2).
(iv). Repeat ii) with the chosen average
pressure ~(1, ) and some averaged convecting
velocity u(l~l) to obtain a new estimate of
u(l,2), then 4(1,2), U(1,2) and p(l,2).
(v) Repeat the above until u(l~n , p(l,n)
remain essentially "unchanged" according to
some specified error norm. Take these
converged values as u(to + fit) and p (to +
it). They will serve as the initial data for
advancing the solution to the next time step
to + 2bt, etc.. As a standard initial value
problem without subsidairy conditions
iterative convergence can be expected with
many standard difference schemes (both
temporarily and spatially) for some finite
time interval T.
We can solve (5)-(7) iteratively in the
vorticity formulation in a similar manner. A
residual divergence of vorticity is equivalent
to an impulsive couple (i.e equal and opposite
pair ot impulsive forces). The relaxation of
such impulsive couples will likewise be
"potential" and the equilibrated vertical
field can be represented as
w(1) (to + St) = w(1) (to + fit) + v ~ (15)
with the potential function ~ governed by
v2 ~ = div w(1) (16)
to be solved with the appropriate mixed
boundary formulation for nontrivial problems.
If there is an external nonsolnenoidal
vorticity disturbance, take it to be the
w(1,1). If there is no external disturbance
advance the initial solenoidal vorticity w(O,
to w( ,1) with the vorticity transport
equation (S) under some approximate boundary
formulation. Then obtain ¢(1,1) from (16)
and w(1,1) from (IS) for future iterative
corrections to obtain the converged value of
w( , ). We find it much more convenient to
obtain the disturbed velocity field u as
solution of (7) rather than as the induced
Velocity evaluated from Blot & Savat Law.
u( ~ ) is then rendered solenoidal by v}
according to (13) and (14). If the updating
of the boundary data for the solution of (S)
involves pressure, we would solve (4a) to
obtain p(l). Otherwise, we can go direct to
equation(S) for up-dating w(1,1) to w(1,2) etc
with some average values of u and p until
(l,n) u(l,n) and p(l,n) converge before
proceeding to construct solution at to + 2bt.
In our earlier computations at rather
coarse meshes, we used the vorticity
formulation which avoids the differentiation
of the numerical velocity data to obtain
vorticity. In flow visualization experiments,
the flow fields are generally described in
terms of vorticity. We know of no
experimentally measured pressure transients
for checking or comparing with the presently
computed results. The most time consuming
part of the solution in vorticity formulation
is to solve the 6 Poisson equations (3 in u-p
formulation). The solution of Poisson
equation is tedious and prone to numerical
instability. We tried many "Fast Poisson
solvers" extended to 3 space; but adopted the
classical method with experimentally
determined relaxation parameter which we found
to be flexible, reliable and actually
considerably faster than most. We adopted
also the simplest difference scheme of
forward time and centered space for
discretization. Many discrete treatments of
the outer boundary in integrating (5) or (2)
have been tried to minimize the propagation of
boundary errors that limits the number of time
steps of computation before the flow field
development in the center of the computational
region is clearly affected. We tried many but
adopted the format of simple extrapolation
normal to the boundary.
274
OCR for page 275
We studied computationally the
development of disturbances in a uniform Shea
flow field between two parallel plates wtih
the top plate moving in its own plane with a
unit velocity relative to the stationary
plate. The flow Reynold number based on the
separation distance between the two plates is
3000, a typical value in the transition range.
Some computations were performed at Re ~ 300
to verify the stabilizing effect of the lower
flow Reynolds numbers. Artificial
disturbances were introduced at a cluster of
six mesh points mostly next to the stationary
plate. We began our study with the coarse
meshes (15-45) x 15 x 15 and At = 0.02 on the
IBM 3081 at Princeton University based on
vorticity formulation. We introduced
different types of artificial impulsive
disturbances of velocity and/or vorticity
components of wide range of magnitudes. They
are artificial or nonphysical since the
disturbed velocity (and/or vorticity) field is
not solenoidal and hence not physical. Each
such artificial disturbance field u(l~l)
and/or w(1,1), is accordingly equilibrated to
a solenoidal velocity disturbance field u(l)
(and or w(1)) and its associated pressure
disturbance p(l) to serve as a physically
valid initial disturbance. Both the
equilibrated u(l) and p(l) are highly
oscillatory over the entire field of
computation, with large peaks and valleys next
to where the impulses are applied. Away from
such peaks, the wide spread oscillations are
an order(s) of magnitude smaller. They are
highly variable in details. For convenience,
we shall refer to each case as the evolution
of the impulsive velocity (or vorticity etc.)
rather than the natural evolution of the
complicated pair of the physical initial data
in u( ), ~(1 and p(l)
For unit impulses of velocity and/or
vorticity components clustered in different
localities, the resulting pressure
disturbances vary considerably in peak
magnitudes. The magnitudes of their
solenoidal velocity disturbances remain 0(1)
while those of their peak pressures vary from
0(10-2) to 0(10-4). We computed the evolution
of dozens of such disturbances. The behavior
of their evolution falls into three
categories, each of which turns out to be
characterized by the initial peak magnitude of
the numerically small equilibrated pressure
disturbances.
(i) A small peak pressure disturbances
of 0(10-4) rapidly decays to restore uniform
shear flow within numerical noises,
commensurate with the available resolution.
(ii) An intermediate peak pressure
disturbance of < 10-3 spreads rapidly to form
some locally "turbulent" region but rapidly
decays to restore some disturbed laminar flow.
(iii) A large peak pressure disturbance
of 0(10-2) rapidly evolves into some
asymptotic pattern of an expanding local
"turbulent" region. The pattern propagates
and evolves slowly in details, reminiscent of
the circumstances revealed in many flow
visualization studies of "turbulent spots" in
transitional flows(7)(8)(9). The asympototic
stage of evolution of such large disturbances
is well formed at ~ 50 time steps in our
computations of 45x15x15 resolution in
IBM3081, well before the boundary errors
appear to distort appreciably the flow field
evolution in the center of the computational
field.
These qualitative features are confirmed
by finer mesh computations on CRAY 1 at NCAR
Colorado, with (31-45) x 20 x 45 meshes and At
= 0.01. We were able to observe some flow
field details of the evolution in these finer
mesh results. Further mesh refinements to
(45-61x20x45) have been carried out on CRAY 2
at Minnesota and C RAY XMP at Illinois for more
detailed studies of the evolution of selected
cases. We recomputed some cases with the
primitive variables, i.e. u-p formulation and
evaluated the vorticity field by numerical
differentiation of the computed velocity
field. These results of vorticity agree to at
least two significant figures everywhere and
at all times with those computed directly from
the vorticity formulation. Thus most of our
recent computations at fine meshes on C RAY 2
and XMP use the u-p formulation to half the
memory requirement and the computing time.
The results reported in the next section have
been computed with both formulations, which
our graphics can not distinguish.
The peak pressure magnitude of a given
initial cluster of impulsive disturbances
decreases rapidly as the cluster receeds from
the plate surface. Thus most of our computed
cases are for clusters in the proximity of the
stationary plate. For such clusters next to
the plate, unit impulses of all vorticity
-opponents turn out to be "small" while unit
impulses of transverse velocity components are
"intermediate" disturbances. Unit impulses of
streamwise velocity component is the only
"large disturbance". Thus, the magnitude of
the artificial impulse appears irrelevant as a
measure of disturbance "strength". The peak
magnitudes of the equilibrated pressure
disturbances correlate well with the eventual
course of evolution as is described above
despite that they are much less than those of
the associated solenoidal velocity
disturbances. This is physically
understandable from Newtonian Mechanics, that
the evolution of the velocity field d/dt u
(or B/8t ui) is primarily driven by the
unbalanced pressure forces, (at least at large
Reynolds numbers). We are therefore much
interested in the evolution equation of
pressure in an incompressible fluid.
275
OCR for page 276
By inverting approximately the Poisson
operator in (4a), we have
~ P + ~ ( Ui ~ Re v Ui) 2xi P
= ~ i i ax · axe P + ~ij k ~xi~x~ ark P ~
(17)
Here ~ ij and ,Bijk are the diffusivity and
dispersion tensors of the field, both as
complicated functions of the instantaneous
velocity field and the flow field boundary and
is a constant depending on the boundary
geometry only. If one dimensional model could
be any guide, the pressure evolution according
to (17) would be dissipative and heat like if
I I ~ij I 1/ | | pi k | | > > 1, and dispersive
and wavelike if 1 | eij | |/ | 1 13ijk 1 1 << 1 in
support of the three different courses of
evolution as depicted by our computed results.
We are reasonably confident in the
eventual convergence of our algorithm in the
limit of zero mesh sizes without smoothing and
filtering; but not so confident in the
quantitative aspects of our results of such
coarse mesh computations. Neither are
quantitative experimental results available to
check and guide our computations. The global
features of flow field evolution described
above are common to our results at different
mesh resolutions, and appear to be in
agreement with experimental observations. In
the next section, we report some flow field
details as have been observed in our finer
mesh results.
IV. Evolution of Large Disturbances in a
Uniform Shear Flow over a Wall
We focus our attention on the evolution
of unit impulses of streamwise velocity
components in the immediate vicinity of the
plate (j51) applied at a cluster of 6 mesh
points (I.J.K. ) = (20, 2, 14-16) . The
equilibrated pressure field over section J=2
is illustrated in Fig. 1 with two distinct
peaks of magnitudes ~ 3x102. The pressure
over the transverse section 1:20 inbetween
the two peaks is shown in Fig. 2 and that over
section K=14 (off from peaks in K=15) is given
as Fig. 3. The point (20, 2, 15) appears to
be a " Saddle " inbetween the two peaks . There
are two valleys in the transverse section. The
two distinct high pressure regions over the
plate surface appear to be two ne ighbouring
isolated globules or half domes encased in a
low pressure valley in the shape of some semi-
spherical caps in three space. This is
typical of "large" initial disturbances that
promptly evolve into propagating local
turbulent regions. When the magnitude of the
initial artificial u1 impulses is reduced to
10 ~ 1, it becomes marginally an " intermediate "
disturbance that fails to reach the asymptotic
state of development. Likewise a set of unit
impulses of up normal to and directed away
from the plate is " intermediate " or marginally
large. Unit U2 impulses directed toward the
plate or US impulses are " intermediate " and
marginally "small. " All unit vorticity
impulses are "small" disturbances that decay
monotonically. The detailed profiles of all
these disturbances, vary considerably but
without significant implications on the course
of development of such disturbances .
We take the evolution of the set of unit
streamwise velocity impulses u1 as the test
case for detailed study. The evolution
history of the peak pressure over section J=2
is illustrated in Fig. 4. The preciptous drop
of peak pressure in the first few time steps
is accompanied by a broadening of the high
pressure region and signficant changes in the
shape of the profile, including the
"disappearance" of one of the initial peaks
with the creation of many others around it.
An asympototic shape begins to form at about
20/\t, followed by some rise and fall in the
peak magnitude with significant broadening of
fluid volume under high pressure. The peak
pressure magnitude rises and falls but manages
to stays at ~ 10- 3 up to 200/`t, the longest we
have computed so far. Each time step (0 .01
H/U) at Re ~ 3000 corresponds to a phys ical
time ~ 10-5 sec. Our computation has thus far
covered only a couple of milk - second in
physical time during much of this period the
pressure transient in Section J=2 has been
developing in some complicated "asympototic"
or " quas i - steady" form .
The pressure elevation in Section J=2 at
time step 20^t is given as Fig. 5. The
downstream initial pressure peak has all but
disappeared while the upstream initial
pressure peak at I=l9 of magnitude 0. 39x10- 1
has become the only dominant peak at 1=20,
K=15 of magnitude 0. 96x10-2. There are now
four lower peaks surrounding the main peak
and many smaller ones further out forming
almost a ring of emerging peaks. The pattern
is evolving slowly, convecting downstream and
spreading out with additional "rings". Fig. 6
displays the pressure distribution at 20^t
over a vertical section perpendicular to the
stationary plate at K=14, i.e. one Liz off from
the peak in the J=2 plane to illustrate the
three dimensional nature of the many peaks and
valleys.
The pressure elevations over Section J=2
at the successive time steps, 40, 80, 120, 160
and 200 At are given as Fig. 7-11. The
primary peaks rise and fall as they propagate
outwards so that the location of the largest
pressure peak can suddenly shift in the plane
J=2. The secondary pressure peaks proliferate
and also very in magnitudes extensively
especially at later times. Their general
pattern as an expanding entity of high
pressure regions remains unchanged despite the
great variations in details and the rapid
increase of the number of secondary and/or
emerging pressure peaks. The pressure
276
OCR for page 277
bans lent has apparently reached some
asympototic stage of dynamic development. The
proliferation and the conspicuous variations
of the outer emerging peaks clearly suggest
actual wave dynamics ( rather than numerical
and/or graphical uncertainties ) . The pressure
elevations over Section K=14 remain little
changed through out the evolution in so far as
the variation in J are concerned. The dynamic
development of the pressure field appears to
lie largely in the streamwise direction.
Therefore, the diss ipative K. D . V . equation
at + P ~ = ~ ~ ~ ~ ~ (18)
with its well known asympototic behaviors),
at least in the nondissipative limit, may
serve as a useful model of ( 17 ) for studying
the global qualitative aspects of the
pressure evolution in a transitional flow
field as mentioned in the previous section.
To better appreciate the three
dimens tonal aspects of the pressure bans tents
next to the plate surface we give in Fig. 12
the colored 3D elevations of the pressure
field in J=2 from different view angles. The
color illustrates the magnitude corresponding
to the scale in the figure. We note that the
valleys appear as deep as the hill' s height.
The streamwise pair of the peaks are always
accompanied by a crossstrenm pair of valleys
to form a "quadruple" and the number of such
"quadruples" increases at later times. The
pressure elevations at planar sections with
somewhat larger values of J are s imilar
although with different magnitudes and
distributions. Thus the dynamics of pressure
transients appears to be describable by the
evolution of such "quadruples" across
localized half dome of high pressure region
encased in a shell of valley in three space.
As such it might be possible that the dynamics
of each quadruple may be governed by some
equation like the dissipative K. D.V. equation
(18) .
The evolution of the vector velocity and
the vector vorticity field appears much more
complex and extended much further out than
that of pressure. It is difficult to describe
with a few planar sections. We are unsure
that the movie in preparation can describe
fully the complexities of the evolution of the
contorted fields of u and w. The local
variations of velocity with time are chaotic,
similar to those "measured" hot wire
responses. There appears, however, some
"Order" out of the "Chaos" when viewed
globally. Such a global order is likely
associated with the quadruple structure of the
pressure disturbances. We describe first the
initial velocity impulse in Fig. 13 over the
transverse section I~20 inbetween the two
distinct pressure peaks at I '19 and 21
respectively (Fig. 1). The curved lines in
Fig. ,3 are formed by the projections of the
local instantaneous velocity vector at a point
in the Section 1=20 onto the section. It is
designated as "Streamlines " in the figure for
convenience, which describes the perturbation
velocity in the section over the uniform shear
flow u=y perpendicular to the plane. Note
that the initial artificial U1 impulses are
applied at j=2, while the equilibrated
solenoidal disturbed velocity field appears to
" converge " to a much higher pa int at
j ~ 7 to 8 . This "node " might be interpreted
as a " s ink" in the transversal surface and a
" source " for the longitudinal ( s treamwise )
flow. It is also apparent in Fig. 13 that the
flow next to the surface below this node is
too complex for our computation to resolve. A
similar projection of the velocity vector for
surface J=2 (Fig. 14) suggests the presence of
counter rotating spiral fluid motion at this
"node " .
This node in the transferal plane I~20
shifts toward the surface and reaches halfway
at But. By 14At the transversal flow appears
to collapse toward the surface as is
illustrated in Fig. 15. This collapse is
strengthened and widened at 20At suggesting
the formation of a counter-rotating vortex
pair next to the wall (Fig. 16) that becomes
conspicuous at 120At (Fig. 17). This vertical
pair disappeared from the surface at 140At
(Fig. 18). It is apparently lifted from the
wall, displaced outward and replaced by a jet
of fluid toward the surface as in Fig. 15.
Such a transverse flow pattern is evident in
Fig. 19 at 160^t and Fig. 20 at 200At, with
the emergence of a new set of rapidly
expanding counter - rotating vortices of
opposite sign midstream. A strong transversal
flow away from the surface ( i . e . the
"ejection") is accompanied by the fluid flow
down toward the plate surface, displacing the
existing complex flow pattern next to the
surface outwards. The successive events
appear to result from the arrival of a local
high pressure globule with its associated
vortex pattern in the transferal section
I .20. Under the "convective" motion, such
high pressure "domes" or "quadruples" carry
with them the spiral flow pattern generated
the large local pressure gradients. The
successive patterns at 1=20 reflect the
simultaneous flow patterns over the various
neighbouring transverse sections as a high
pressure "dome" or "quadruples " passes by
according to the dispersion dynamics of
pressure described by equation (17).
Fluid elements, "collapsing" or
"converging" toward a high pressure "dome" are
deflected sidewise to fall into a pressure
valley on the side, temporarily trapped in the
valley. They may escape and then trapped
again further downstream into the next spiral
of a neighorbouring oncoming pressure dome.
Or they may escape collectively as a vigorous
j et away from the plate in the absence of an
immediate oncoming high pressure dome, i . e .
277
OCR for page 278
"Ejection". Different fluid elements will
undergo "spiral" motions of different
dimensions and intensities from "Collapse" to
"Ejection". Such an event will generate
intense chaotic fluid motions further out.
There results then a turbulent or chaotic flow
region centered around an expanding array of
spiral vortices with repeated collapses and
ejections next to the solid surface. The
phenomena of "collapse" and "ejection" and the
formation of counter rotating axial vortex
pairs have been widely observed in flow
visualization experiments). We are much
encouraged to have reproduced all these
without numerical artifices. It appears that
the "Burst" into local flow turbulence at a
point is associated with the passage of a high
pressure dome or quadruples over this
observation point.
The motion of the pattern of dispersive
pressure waves described by equation 17 is
much different from the local instantaneous
velocity ui. Where the local pressure
gradient is large with sufficiently large v2
Ui or v x a, the motion of the high pressure
region can differ much from the local fluid
velocity ui. Thus fluids are "drawn" toward
the pressure quadruples, "stirred up" through
the spirals and "ejected" into neighbouring
flow region as chaotically moving fluids. The
quadruple nature of the pressure field
transients appears to be crucial in driving
the laminar to turbulent transition. The
substantial velocity differential of the fluid
and the quadruple motion is responsible for
the spread of the "chaos" or "turbulence"
beyond the confine of the propagating pressure
quadruples.
When the water (or other liquid)
saturated with dissolved air (or other gases)
is swept over by such pressure quadruples, the
rapid decompression accompanying the arrival
of a pressure quadruple will lead to the
formation of "air bubbles". Such air bottles
moving with the fluid will escape the
quadruples and remain in the fluid after the
passage of the quadruple. Accordingly the
proliferation of such quadruples in a
propagating local turbulent region becomes a
powerful source of "Cavitation". The
pressure elevation on the solid surface (j=l)
around a pressure quadruple is well
represented in Fig. 12. Air bubbles are
generated in the pressure valleys; they need
not "collapse" there. The cavitating bubbles
can escape into the surrounding turbulent
region, be drawn together (Fig. 14, 15) by
some oncoming quadruples, and "coalesce" into
large bubbles before their eventual
"collapse". The collapse of a reasonably sized
bubble can be a "large" pressure disturbance
of magnitudes comparable to 10-2 - 10-3 pU2
i.e the local Froude Number based on the
bubble diameter U2/gD < 102 or 103. Such a
collapse generates a new local turbulent
region to perpetuate the propagation of
278
chaotic flow as flow transition and
turbulence. With the magnitudes of the
pressure peaks and valleys generally some
fraction (> 10-3) of the dynamic head (pU2) of
the relative motion of the solid surface with
the "outer free stream", pitting would
naturally be more serious nearer to the tip of
a propeller or for propellers of higher
speeds.
A pressure quadruple in the proximity of
the plate surface around a given point gives
rise to a surface pressure pattern like Fig.
12. This surface pressure pattern will set
the plate into vibration in some
characteristic manner. The multiple of such
convecting pressure quadruples in a local
turbulent region is thus a powerful
"acoustic" source, with identifiable
characteristics, directionally, spectrally
and/or in selected correlation functions to
distinguish itself from the prevailing noisy
environment. Better understanding of the
quadruple structure of the propagating local
turbulent region in a transitional flow field
could have far reaching implications. The
absolute magnitudes of the "large"
disturbances that will generate propagating
pressure quadruples ~ 10-3 pU2 at Re ~ 3000
are of the order of 10~1 Kg/m2 (or 2 x 10~1
lbs/ft2). They are comparable to the pressure
differentials produced by a rather small or
even insignificant free surface waves in open
sea. Thus at Re > 3000, turbulent flow
conditions will likely prevail over a
submerged solid surface. A lOOdb noise in
air will generate such quadruples at Re ~ 3000
to cause flow transition.
IV. Concluding Remarks
We developed an algorithm to compute the
evolution of the pressure and the velocity
disturbances in a flow field of an
incompressible fluid with the Navier-Stokes
equations system through iteration in
solenoidal subspace(s). The discrete
equivalent of the initial value problem can be
obtained through any simple, consistent and
stable schemes with At satisfying the
condition of zone of dependence without any
numerical aritifices of smoothing, filtering
and/or damping. We computed the evolution of
many artificial impulsive disturbances of
different types and magnitudes in a uniform
shear flow between two parallel plates at Re =
3000. Computations have been repeated at
successive mesh refinements and both in terms
of the primitive variables (u,p) and with the
help of the derived variable vorticity w.
The eventual course of evolution depends
primarily on the peak magnitude of the
equilibrated pressure disturbance, not so much
on the profile details, nor on the nature of
the artificial disturbances generating it.
Small disturbances are dissipative and stable.
OCR for page 279
Intermediate disturbances are diffusive,
spreading out and decay rapidly.
Sufficiently large disturbances will propagate
and proliferate to become asymptotically a
multitudes of propagating local high pressure
regions each surrounded by valleys. Over the
planar section next to the stationary plate
they appear as quadruples with spiral vortices
attached to the sides and fluid ejection and
collapse fore and aft. These convecting
quadruples generate "chaotic" motions of the
surroundings fluids. The flow is in
Transition from the "Laminar" to the
"Turbulent" state. Such asympototic "Order"
persists, nevertheless, within the apparent
tl Chaos".
The evolution of the pressure field can
be described by a three dimensional analog of
the dissipative KDV equation possessing the
different limiting properties described above.
As such the propagating large pressure
disturbances are likely dispersive waves whose
"convective" velocity can be much different
from the local instantaneous velocity of the
fluid. If the "Solitone" like character of
the solution of KDV equation should prevail,
"Turbulent" flows would retain such "orderly"
transitional structure, even if not as the
strict superposition of such "quadruples".
Such convecting quadruples of pressure are
powerful acoustic sources that may possess
distinct characteristics for recognition or
identification. For water "saturated" with
air, a pressure quadruple is a source of
cavitating bubbles to help spreading
turbulence through the collapse of coalesced
bubbles.
The author wishes to acknowledge the
support of the National Science Foundation
under grant MAE 80-10876 and MSM 83-12094, and
to thank Dr. Sylvian Roy and Ms. Min Chen and
Mr. Jerry Syms in carrying out the
computations and data analysis reported in
their theses.
279
OCR for page 280
a
5
~ ~~ 1 1
MPXlMUM PRESSURE = 0.Z474E-05
LIE Ha ~ M~MUH PR~dURE
ON [~^12 Si]CE j at- 2
~ -,- --_ - - . -. -.-, .
MnX l MU E ~ 0. 3511E - C11
.. ..
ICE STEP
Fig. 4
E = ~~
Eli 3 280
Fig. 6
OCR for page 281
~ ~~ ~ ~ z
~ ~~ ~ ~ z
~1~ ~~ ~ 0.~1~-~
Fig. 7
~~ ~ Z
TREBLE Fit CnSE 51
MnXIMUM PRESSURE ~ O.ZQ53E-OZ
Fig. 8
Z
= ~
~1~ ~ = 0.~1~-~
Fig. 9
~1~ = ~ 1
Fig. 10
PRESSuRE
~ ~ ~ ~ =
~~ ~ Z
~1~ - ~ - ~
Fig. 11
281
OCR for page 282
282
OCR for page 283
Representative terms from entire chapter:
residual divergence
sTREaMLlNEs RT 1= 20
THE TIME STEP IS - I
T0R3.J C~T SP2T ~E\--R~T12\ F:R CPC- S;
—
Fig. 13
STREqML 'NES RT J- 2
TURBULENT SRL4T GtNERqT10N F0R CRSE Sl
THE TIME STEP IS - I
Fig. 14
STREAMLINES RT 1= 20
THE TIME STEP IS - 11
TU230!-\l SP2T GEMER-Tl2N FZ? CRSF S
l
\\~\lll/ /W
Fig. 15
STRERt
REFERENCES
Fig. 19
1. Lamb, Hydrodynamics, 6th Edition, Dover
Publication. A general reference and
specifically Art .1, p. 10. Art. 19. pp. 161-2
and Art. 152, pp. 214-216 (1932).
2. Nichols, B. & Hirt. C., "Calculation of 3D
Free Surface Flows in the vicinity of
submerged and exposed structures". Jour. of
Comp. Phy. Vol. 12, no. 173, pp. 234-246.
3. Rimon Y. & Cheng, S. I., "Numerical
Solution of a Uniform Flow over a Sphere at
Intermediate Reynolds Numbers", The Physics of
Fluid, Vol. 12, No. 5 (1969).
4. Leonard, A., "Vortex Stimulation of 3D
Spotlike Disturbances in a Laminar Boundary
Layer", Second Symposium of Turbulent Shear
Flows", London (1979) or NASA TM 78579.
a. Temam, R., "Navier Stokes Equation", North
STREAMLINES RT I- 20 Holland Elseviers, (1978).
TU2BJ!ENT SPIT GENERRT 1 IN FeR CASE 51
Fig. 20
6. Chorin, A. J., "Numerical Solutions of the
Navier Stokes Equations", Math. Comp. pp. 745-
762, (1968).
7. Hinze, J. O., "Turbulence", McGraw Hill,
Second Edition (1975). A general reference on
Turbulence and specifically on Transition and
Turbulence Structure, pp. 605-614, Fig. 7.7,
p. 611, pp. 639-641, p. 647, p. 659-666 and p.
724-742.
8. Cantwell, B. J., "Organized Motion in
Turbulent Flow" Section in Ann. Review Fluid
Mech. 13, pp. 457-515, (1981), Annual Reviews
Inc.
9. Kline, S. J., Reynolds, W. C., Schraub, F.
A., Runstadler, P. W., "The Structure of
Turbulent Boundary Layers", J. Fluid
Mechanics, 30, pp. 741-773, (1967).
10. Lax, P. & Levermore, C. D., "The Small
Dispersion Limit of the Korteweg-deVries
Equation I", Comm. on Pure and Applied Math.,
Vol. XXXVI, pp. 253-290 (1983).
284
