Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
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.
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
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
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 < 0 (or > 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
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
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
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
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
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
"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.
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
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
~ ~~ ~ ~ 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
282
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<L ~ NFS PT I - 20 THE T IME STEP IS - 20 C~"l ENT S~Zi GENC2mTl0N FZ? C~'E 51 Fiz. 16 STREnMt.INES RT 1= 20 THE T I ME STEP I S - 120 R~cNT 5P2T E_~-~, i2\ F22 CqS- Fig. 17 STRERML I NEs RT I - 20 THE TIME STEP IS - 150 TtJ?a-, ENT S'3T GFNERPTiZ~ F22 CR~E S~ - - j~-,-'~ =-~//~= Fig. 18 283
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
