A detailed description is given of a method for obtaining numerical solutions to nonsteady incompressible fluid flow problems. Included are discussions of stability properties of the finite difference equations and boundary and initial conditions appropriate to current applications of the numerical program. Solutions are given for flows about obstacles in. a channel at various Reynolds numbers, with emphasis given to tbe process of development of the Karman vortex street.
The author is indebted to Professor A Thom, not only for the experimental. photographs of Fig 12 and Fig 13, but also for his pioneering work in numerical methods which has been drawn upon in the development of the method here used. I wish to thank Dr A M Lippisch for having performed the experiments which have provided us with the results shown in Fig 29. I am also much indebted to my superior, Dr F H Harlow, for his encouragement and assistance in many aspects relating to the develapment of the numerical method.
Fig. 1 Layout of the lattice of points, indicating where variables are defined relative to the mesh and the geometry of the channel.
Fig. 2 Vorticity patterns of unstable growths of a perturbation using the time-centered equation.
Fig. 3 Samples of initial solutions displayed in terms of stream and vorticity functions.
Fig. 4 X component velocity profiles and vorticity contours for a fluid initially at rest with walls moving in its own plane.
Fig. 5 Velocity profiles, vorticity contours, and streamline patterns after fluid interaction with obstacle.
Fig. 6 Comparison of analytic and calculated vorticity profiles in a fluid adjacent to a wall moving in its own plane.
Fig. 7 Vorticity and streamline plots for a fixed obstacle in a moving fluid.
Fig. 8 First part of streamline sequence for R=100. d/h = 1/6.
Fig. 9 Second part of streamline sequence for R=100.
Fig. 10 First part of streakline sequence for R=100.
Fig. 11 Second part of streakline sequence for R=100.
Fig. 12 Qualitative comparison between examples of experimental and calculated streakline patterns.
Fig. 13 Qualitative comparison between examples of experimental and calculated streakline patterns.
Fig. 14 First part of streamline sequence for R=200. d/h = 1/6.
Fig. 15 Second part of streamline sequence for R=200.
Fig. 16 First part of streakline sequence for R=200.
Fig. 17 Second part of streakline sequence for R=200.
Fig. 18 First part of streamline sequence for R=300. d/h = 1/6.
Fig. 19 Second part of streamline sequence for R=300.
Fig. 20 Streamline sequence for R=300.
Fig. 21 Transformed streamline sequence for R=300.
Fig. 22 First part of streamline sequence for R=1000. d/h = 1/6.
Fig. 23 Second part of streamline sequence for R=1000.
Fig. 24 Streakline sequence for R=1000.
Fig. 25 Transformed streamline sequence for R=1000 and R=6000.
Fig. 26 First part of streamline sequence for R=6000.
Fig. 27 Second part of streamline sequence for R=6000.
Fig. 28 Streakline sequence for R=6000.
Fig. 29 Experimental streakline patterns for R=6000.
Fig. 30 Streamline sequence for R=100. d/h = 1/12.
Fig. 31 Streamline sequence for R=200. d/h = 1/12.
Fig. 32 First part of streakline sequence for R=100.
Fig. 33 Second part of streakline sequence for R=100.
Fig. 34 First part of streakline sequence for R=200.
Fig. 35 Second part of streakline sequence for R=200.
Fig. 36 Transformed streamline sequence for R=200.
Fig. 37 Transformed streamline sequence for R=300. d/h=1/12.
Fig. 38 First part of streakline sequence for R=300.
Fig. 39 Second part of streakline sequence for R=300.
Fig. 40 Streamline sequence for R=300.
Fig. 41 Streamline sequence for R=1000. d/h = 1/12.
Fig. 42 First part of streakline sequence for R=1000.
Fig. 43 Second part of streakline sequence for R=1000.
Fig. 44 First part of a streamline sequence of flow about a projection from a fixed wall.
Fig. 45 Second part of a streamline sequence of flow about a projection from a fixed wall.
Fig. 46 First part of a streakline sequence of flow about a projection from a fixed wall.
Fig. 47 Second part of a streakline sequence of flow about a projection from a fixed wall.
Fig. 48 Various representations of the flow past an obstacle for R-200. d/h=1/6.
Fig. 49 Various representations of the flow about a projection from a fixed wall.
Fig. 50 Patterns of unstable growths of a perturbation using a nontime-centered vorticity equation.
With the advent of the high-speed computer, it has become poesibl.e to develop methods for studying theoretically many of the nonsteady incompressible fluid flow problems which previously had been hopelessly complicated for analysis. The purpose of this report is to describe in detail. one such numerical method, and to illustrate its applicability through a munber of examples. Though this new computer method has depended upon analytical and experimental results to guide its early development, it now provides a powerful tool for attacking many problems which cannot be approached by analysis nor understood by experimental observations alone.
The numerical method can fill a gap betveen the limiting cases that have been studied analytically and the more practical results obtained by the experimentalist, Hopefully it will provide new insight to the theoretician. It can serve as a tool not only in terms of end results, but also to test the validity of exact methods. To engineers it can become an aid in planning experiments and in interpreting results. Also the numerical method may serve to unify much of the diverse empirical data which currently exists in the literature.
A finite difference method is used. The derivatives of the Helmholtz vorticity equation are replaced With central differences both in space and time. The resulting algebraic equation is used to obtain advanced values or vorticity at all points in a lattice in which initial values are known. These values then constitute a source term in Poisson's equation, which is satisfied by the stream function. A successive approximation method is used to simultaneously evaluate the stream function at all points of the lattice. The advanced values of the stream function in turn lead to velocities which may again be used to advance vorticity values to a new time.
We consider two-dimensional flow of a fluid relative to a rectangular frame of reference. The mesh is composed of square cells whose junction points are points at which the stream function and vorticity are evaluated. An interlaced system of points is used to prescribe the velocity components. Problems considered for examples in this report are those of a fluid flowing betwen parallel walls with periodic end boundaries. Obstacles or projections from walls affect the flow.
The computer used is the IBM Electronic Data Processing Machine Type-7090, with 32K memory.
This is not the first application of a numerical method to the study of incompressible flow. Several authors have applied numerical methods to steady flow problems. In most of these cases desk calculators were used (see, for example, Thom [1], Allen and Southwell [2], and Jenson [3].) Recently Payne [4] developed a numerical method for nonsteady flows using the Manchester University Mark I Electronic Computer. Payne's method differs from the method described here in that velocities of the flow field are obtained directly from the vorticity distribution by numerical integration.
A group at Los Alamos [5] has employed a method suggested by von Neumann for computing the time-dependent flow at an interface between two inviscid fluids. Also, Thomas [6] has formulated a method in which a truncated Fourier expansion of the field variables is used to represent approximately the fluid configuration.
The finite-difference mesh used in the examples discussed is composed of 1425 lattice points. The program uses input data to determine the configuration of walls and obstacles. Every point of the mesh requires a storage space in the memory which serves to indicate the type of treatment the variables at that point are to be given. The input program generates indicator numbers which form correspondence with the configurat1on being studied. In the main body of the program the indicator numbers are tested, and the implied treatment of the points is hence determined.
Output data are plotted by use of the Stromberg-Carlson SC-4020 Microfilm recorder. Plot data are placed on magnetic which is processed by the recorder. In all plots, the boundary and obstacle configurations are first established, and then contour patterns are superimposed upon the plot. Contour plots of any of the variables can be made. Plots are made at equally spaced time intervals, using some multiple of the timee step increment. Values ot the functionals of the motion and other quantities of interest may simultaneously be printed on film with corresponding plots.
Independent plots are made of motions of test particles. These particles are identified only by their coordinates, which are influenced by the velocities at each step ot the computation. The familiar streak-line patterns are produced in this way.
In the examples given the patterns are those obtained through the normal processing procedure. In no case have the plots been retouched or altered in any way. Computation of a given problem may take as much as four hours, about 10%, of this time being used for data output. Since in most cases computation can continue indefinitely, problems are terminated when nothing more is to be learned.
This report may be divided into three somewhat distinct areas. First, a major part o:f the discussion concerns details of the calculational method. This discussion involves Sections 1 through 5 for the basic calculations with incidentals to the calculation appearing in Section 7-A and Sections 8-A and 8-B. Additional difference forms of the vorticity equation are discussed in Section 9. Second, a summary ot equations and procedures is given in Section 6 to provide a programmer with the essentials of the calculational method, Third, discussions are given of the solutions obtained. These are found in Section 7 and Section 8-C. Although it has not been possible to divorce these discussions from all aspects of the numerical method, it is hoped that some benefit may be derived from study of the solutions and examples without a detailed review of all the preceding material.
The Eularian form of the equations of motion of a viscous fluid in full generality are
ρ(∂/∂t + u⃗ ⋅ ∇) u⃗ = ρg⃗ - ∇P + ∇(λ∇⋅u⃗) + 2(∇ ⋅ μ∇)u⃗ + ∇ ×(μ∇ × u⃗ ) (1-1)
See Harlow [7]. Here u⃗, ρ, P and t are, respectively, the velocity, density, pressure, and time.The quantities μ and λ are the coefficients of first and second viscosity, respectively, and g⃗ is a body force. We take μ and λ as constants, in which case (1-1) becomes
ρ(∂/∂t + u⃗ ⋅ ∇) u⃗ = ρg⃗ - ∇P + (λ+μ)∇(∇⋅u⃗) + μ∇2 u⃗ (1-2) The continuity equation is ∂p/∂t + ∇⋅(ρu⃗) = 0 (1-3) or dp/dt + ρ(∇⋅u⃗) = 0 (1-4) where dp/dt = (∂/∂t + u⃗ ⋅ ∇) p
If the fluid is incompressible, then
dρ/dt = 0 or (∇⋅u⃗) = 0 (1-5) In such a case (1-2) becomes ρ(∂/∂t + u⃗ ⋅ ∇)u⃗ = ρg⃗ = ∇P + μ ∇2u⃗ (1-6) and if we include the equation (∂/∂t + u⃗ ⋅ ∇)p - 0 (1-7) then (1-5), (1-6), and (1-7) form a complete set. Now consider (1-6) in the form ∂u⃗/∂t + ½∇u⃗2 - u⃗×(∇×u⃗) = - ∇P/ρ + ν∇2u⃗ + g⃗ (1-8) where ν=μ/ρ is the kinematic viscosity. If we take the curl of (1-8) and define velocity as ω⃗ = ∇ × u⃗ (1-9) we obtain ∂ω⃗/∂t + ∇×ω⃗×u⃗ = - ∇ (1/p)×∇P + ν∇2ω⃗ (1-10)
We have taken g⃗ as a conservative force and hence it does not appear in (1-10).
If the density is everywhere the same, (1-10) becomes
∂ω⃗/∂t + ∇ × ω⃗ × u⃗ = ν∇2ω⃗ (1-11) This may be written in terms of dyads as ∂ω⃗/∂t + ∇ ⋅ (u⃗ ω⃗ - ω⃗ u⃗ - ν∇ ω⃗) = 0 (1-12) and is hence seen to be a conservative equation. If we now take the divergence of (1-6) we obtain for an incompressible fluid ∇ ⋅ (u⃗ ⋅ ∇) u⃗ = - ∇ 1/ρ ⋅ ∇P - 1/ρ ∇2p + ∇ ⋅ g⃗ (1-13) and making use of the identity ∇ ⋅[∇ ⋅(u⃗ u⃗)] = (u⃗ ⋅ ∇)u⃗ + u⃗ (∇ ⋅ ω⃗) we can write ∇ ⋅[∇ ⋅(u⃗u⃗)] = - ∇ 1/ρ ⋅ ∇P - 1/ρ ∇2P + ∇ ⋅ g⃗ (1-14) If both the density and the body force are everywhere the same (1-14) becomes ∇2P = - ρ ∇ ⋅ [∇ ⋅(u⃗u⃗)] (1-15)
In the applications to follow we shall restrict the discussion to a one-material, homogeneous fluid. The basic equations to be considered are (1-5), (1-9) and (1-12).
Since ∇ ⋅ u⃗ = 0, we may define a vector A⃗ ≡ ξ i⃗ + η j⃗ + ψ k⃗ such that u⃗ = ∇ × A⃗ Considering a case where u⃗ has the components u in the x-direction and v in the y direction, we may write u = ∂ψ/∂y (1-16) v = - ∂ψ/∂x (1-17) where ψ is the stream function. In two dimensions the definition (1-9) becomes ω = - ∂u/∂y + ∂ν/∂x (1-18) or in terms of the stream function ∂2ψ/∂x2 + ∂2ψ/∂y2 = - ω (1-19) Also the time-dependent vorticity equation (1-12) becomes ∂ω/ωt + ∂/∂x (uω - ν∂ω/∂x) + ∂/∂x (νω - ν∂ω/∂y) = 0 (1-20)
We could use (1-16) and (1-17) to write (1-20) in terms of ψ and ω only. In such a case the basic set to be solved would be (1-19) and (1-20). We shall, however, make use of (1-16) and (1-17) throughout the discussion to follow.
In the reference frame selected we note that (1-15) becomes
∂2P/∂x2 + ∂2P/∂y2 = -ρ Q (1-21) where Q is defined as Q = ∂2u2/∂x2 + 2 ∂uv/∂x∂y + ∂2v2/∂y2 (1-22)
In this report the pressure will be considered as an auxiliary variable to be evaluated after the basic solution has been obtained.
To solve an initial-condition boundary value problem governed by the equations given in Section 1-B, we approximate the derivatives with finite differences. We prescribe a square mesh of cells through which the fluid flows. The individual cell may be considered to be a volume element of unit thickness in the z direction, in relation to which the fluid must obey certain laws. Since the fluid is incompressible we must require that the net inflow into the volume element be balanced with the net outflow. Consider a cell whose lower left corner carries the index i,j-1 (see Fig. 1). We consider an average inflow velocity on the left face of the cell end identify it as ui,j-½. With similar prescriptions on the other faces, we may write an equation which balances the inflow and outflow in the form
ui+1,j-½ - ui,j-½ + vi+½,j - vi+½,j-1 = 0 (2-1)
This equation is a difference form of the continuity equation. We guarantee that (2-1) will hold and that mass will hence be conserved by prescribing
ui,j-½ = (ψi,j - ψi,j-1) / a (2-2) and vi+½,j = (ψi,j - ψi+1,j)/a (2-3)
where a is the length of the side of the square cell. Note that this prescription makes (2-1) an identity.
The ψ's are positioned at the intersection of the cells. Equations (2-2) and (2-3) are difference forms of (1-16) and (1-17), respectively.
If we now consider the curl of the velocity relative to the point i,j, then, in terms of the nearest velocity values, we may write
ωi,j = (ui,j-½ -ui,j+½ + vi+½,j - vi-½,j) / a (2-4)
where ω, the vorticity, is taken to be positive in a counterclockwise sense. Equation (2-4) is a difference form of the defining equation of ω, Eq. (1-18). Using (2-2) and (2-3), Eq.(2-4) may be written
ψi,j = (ψi+1,j + ψi-1,j + ψi,j+1 + ψi,j-1 + a2ωi,j)/4 (2-5)
This is a difference approximation of Poisson's equation. We shall use (2-5) to obtain ψ values and then use (2-2) and (2-3) to obtain velocities. Equation (2-5) is applied in a successive approximation technique described in Section 3. Note that the variables in the above equations are defined. at points as indicated in Fig. 1.
Now the vorticity values in the fluid vary in time, in that vorticity may be transported through fluid motions or may be diffused tbrough the dissipative action of viscosity. Vorticities may arise as a result of interaction of the fluid with solid objects; otherwise the properties of the vorticity are simply a function of the initial conditions and the transport and diffusion properties. Equation (2-4) is used to establish the properties of the vorticity as related to fluid interactions with solid objects (see Sections 4-B and 4-C).
Time variation of the vorticity is taken into account by a finite difference approximation of (1-20). In the development of the calculational procedure, experimentation was carried out on several difference forms of (1-20) (see Section 9). We give the form currently in use, the form which has been used in the applications of Section 7. We write
ωi,jn+1 = {1 / [1+(4νδt)/(a2)]} {ωi,jn-1 + 2δt/a [ (uω)i-½,jn - (uω)i+½,jn + (νω)i,j-½n
-(νω)i,j+½n ] + 2νδt/a2 (ωi+1,jn + ωi-1,jn + ωi,j+1n + ωi,j-1n - 2ωi,jn-1)} (2-6)
Here we have replaced the time derivative by
(ωi,jn+1 - ωi,jn-1)/2δt
where n is the time index such that ωn = ω(nδt). Note that (2-6) makes use of a centered time difference. This property is esential for sufficient stability for small viscosities (see Section 9).
To obtain the values required on the right of (2-6), which are not positioned at points of definition, we use averages of the given values. Thus
ui+½,j = (ui+1,j+½ + ui,j+½ + ui,j-½ + ui+1,j-½)/4 (2-7) νi,j+½ = (νi+½,j+1 + νi-½,j+1 + νi-½,j + νi+½,j)/4 (2-8) ωi+½,j = (ωi+1,j + ωi,j)/2 (2-9) ωi,j+½ = (ωi,j+1 + ωi,j)/2 (2-10)
Any additional values required may be obtained by permutation of the indices of these equations.
In Fig. 1 we have indicated the presence of an obstacle occupying the space of a single cell. In the applications considered obstacles occupy an integral number of complete cells. By keeping the geometry of obstacles simple, in relation to the mesh, we avoid the concern that (2-1) will not hold.
We shall consider the stability properties or (2-6) in terms of very small perturbations. While such an assumption cannot lead to conclusive results, some insight is to be gained. We shall further siimplify the discussion by considering (2-6) in reduced forms, first with the inertia properties absent and then again with zero viscosity. Finally, the properties of the complete equation are inferred from a series of stability tests.
With the nonlinear terms absent, (2-6) may be written in the form
ωi,jn+1 - ωi,jn-1 = 2νδt/a2 (ωi+1,jn + ωi-1,jn + ωi,j+1 + ωi,j-1n - 2ωi,jn+1 - 2ωi,jn-1 (2-11)
This form of central time differencing was first given by Dufort and Frankel [8]. The stability properties for the one-dimensional form have been discussed by Richtmyer [9]. We here give an alternate discussion and show that the properties are the same in two dimensions.
We consider a typical Fourier component solution to (2-11) in the form
ωi,jn = ω0 eik1i eik2j rn (2-12)
where the i's before the k's are the imaginary number √-1 k1 and k2 are wave numbers of the component variations in the x and y direction, respectively, and ω0 is the reference level of the component. rn is a growth factor of the variations, and we note that the condition for boundedness of (2-12) is
|r| ≤ 1 (2-13)
Substitution of (2-12) into (2-11) leads to
r - 1/r = 2νδt/a2 ( eik1 + e-ik1 + eik2 + e-ik2 - 2r - 2/r) or r2 = 4α/(4α+1) (cos k1 + cos k2 )r + (4α-1)/(4α+1) = 0 (2-14) where α = νδt/a2 (2-15) Thus r = 1/(4α+1) (γ ± √(γ2 - 16α 2 +1) (2-16) where γ = 2α(cos k1 + cos k2) We consider two cases, first γ2 + 1 - 16α2 > 0 for which (2-16) is real and γ2 + 1 - 16α2 < 0 for which (2-16) is imaginary. Note that -4α ≤ γ ≤ 4α
hence in the real case the right-hand aide of (2-16) has its maximum value for γ=4α with a plus sign for the radical, and r ≤1 must be satisfied. With these substitutions we obtain r = 1 or, for this case, a condition of neutral stability.
The right-hand side ot (2-16) is a minimum if γ=4α and the minus sign is used at the radlical. Here r≥-1 must be satisfied. We obtain r=-1 or a condition of alternating stability.
In the imaginary case we must require (rr‾) ≤ 1. Note that (2-16) becomes (rr‾)=(4α-1)/(4α+1) and hence the condition (rr‾)≤1 is satisfied for any α.
With the transport terms absent then, we note that small perturbations will not grow when a central time difference is used as in (2-11).
We now consider (2-6) with zero viscosity. That is
(ωi,jn+1 - ωi,jn-1) / 2δt = 1/a [ (uω)i-½,jn - (uω)i+½,jn + (νω)i,j-½n - (νω)i,j+½n] (2-17)
We assume that the velocity components vary only slightly from steady values u0 and ω0 and write
ui,jn = u0 (1 + δi,jn) δ ≪ 1 (2-18)
and
νi,jn = ν0 (1 + ξi,jn) ξ ≪ 1
If we now make use of (2-4) in terms of (2-18) it will be found that to first order (2-17) becomes
ωi,jn+1 - ωi,jn-1 = 2δt/a [u0(ωi-½,jn - ωi+½,jn) + ν0(ωi,j-½n - ωi,j+½n) ] (2-19)
We may now consider the component solution (2-12) for Eq. (2-19) as was done previously. This leads to
r-1/r = -i4δt/a (u0 sin (k1/2) + ν0 sin (k2/2) )
or
r2 + i4δt/a (u0 sin(k1/2) + ν0sin (k2/2)r -1 = 0 (2-20)
Thus
r = -iξ±√(1-ξ2)
where
ξ = 2δt/a {u0sin(k1/2) + ν0 sin(k2/2) }
Note that if ξ2≤1
rr‾ = ξ2 + 1 - ξ2 = 1
But if we assume ξ>1, for example, then r becomes imaginary and we note that the condition
-ξ -√(ξ2-1) ≥ -1
must be satisfied for stability. This leads to the contradiction
ξ ≤ 1
We conclude then that the region of stability is such that
ξ2 ≤ 1
This same conclusion results from assuming ξ < -1.
We note that the stable region is characterized as one of neutral stability, i.e., it is a region in which small perturbation will neither grow nor damp.
We have noted from the above analyses of the reduced forms of (2-6) that there is a restriction on δt for the ν=0 case, but not for the case in which the transport terms are absent. While the overall analysis of (2-6) has not been carried out, experiments have been performed to check whether any instabilities exist. The experiments do indeed indicate that instabilities will result if the condition
uδt/a < 1
is not satisfied in its two-dimensional form.
Even if the complete (2-6) were unconditionally stable, experience suggests that for accuracy we would have to use a δt such that
β = ( { |u0 + |ν0|}/a) δt ≤ 1 (2-21)
and also
α = νδt/a2 ≤ 1/4 (2-22)
(See Harlow [10].) We hence use these criteria as a basis for experimentation. Equation (2-22) is the two-dimensional stability criterion for the diffusion equation without the central. time difference (see Section 9).
In the experiments a perturbation was applied to a steady state flow u=4, μ=0, and ω=0 everywhere. Periodic end boundaries were used while containing walls moved with the fluid (see Section 4 for a discussion of boundary conditions). The perturbation ω=102, was inserted at t=0 at a single point 1/7 of the field length from the left boundary and centered vertically in the field. Its growth or decay was observed as time progressed.
A series of experiments for β=0.8 and α varying within the range specified by (2-22) showed a decay in the perturbation at t = 100 Δt. Growth of the perturbation was, however, observed for β>0.8. Some typical examples of the form of the growth are displayed in Fig. 2. Here the left-hand frames are vorticity patterns for cases where β=1.
The upper of these is for α=¼, the equals case of (2-22). The center frame is for α=0.15 and the lower frame for α=0. The upper and center frames are records taken at 40δt, while the lower was taken at 170δt. In the example of δ=0 growth was indicated at t=100δt but appeared to be very slow. The growth was four cells in wavelength with no growth at alternate mesh points.
The central frame of this left group shows a more rapid growth than for either α=0 or α=¼. The plot interval in all of these cases was Δω=0.05. Both α=¼ and α=0.15 cases have growth components of four cells in wavelength as does the α=0 case.
The series on the right of Fig. 2 are for a constant α corresponding to that of the center left frame α=0.15. The center right frame is a later time record of the center left frame, being taken at 80δt. Here Δω=3.2. The upper right frame is for β=0.9 and shows a much slower growth since Δω=0.05 and the plot was made at 90δt. By contrast, the lower right frame shows a very rapid growth of the perturbation by 10δt for β=1.6. In this case Δω= 12.8.
In the application of the computer method then, we must require β≤0.8. It has been found, however, that there is no real necessity to make δt as large as possible. This stems from the fact that the iteration procedure consumes the major part of the computing time (see Section 3). It has been noted that a reduction of δt by a factor or two permits convergence of the successive approximation solution to occur in almost half the time, because the new solution is displaced less from the old solution which is always used as a first guess. We here use (2-21) and (2-22) with a factor 0.45 on the right to govern δt for the calculation.
Judgments about the method to be used in obtaining the ψ solution, using difference approximations to Poisson's equation, are many and varied. The view taken here is that this part of the calculation is to be handled in the simplest possible way. One may argue that convergence rates of more subtle methods would be a vast improvement. This is, of course, entirely possible. On the other hand, if the calculations involve a large number and variety of decisions at each step of the calculation, more rapid convergence may not necessarily be an advantage in terms of computer techniques. For example, if the decision-makng instructions are vast and complex, tests may consume computing time that could otherwise be used to advantage in making more sweeps of the field. Also, the decision-makng apparatus will usually require modification if one is to handle different boundary conditions and configurations.
The method used here is based on the simple procedure that has come to be kow as Liebmann's Method. Liebmann used the method in connection with Laplace's equation in the form
ψi,jk+1 = ¼Σ of the kth iterate at neighboring points
A simple improvement applied by many workers, which is compatible with computer methods, is to use any available advanced values of neighboring points. Thus if the field of values is swept from left to right as text is read on a printed page, then (2-5) may be written
ψi,jk+1 = (ψi+1,jk + ψi-1,jk+1 + ψi,j+1k+1 + ψi,j-1k + a2ωi,j)/4 (3-1)
Use of this procedure requires only a single storage section for the ψ field in the computer.
Equation (3-1) is applied in successive sweeps of the field, at each time step, until the field of values has settled. Tests of convergence, to be discussed, are made periodically to determine whether additional sweeps will improve the solution. If, for some reason, convergence to within the desired accuracy does not take place after 500 sweeps of the field, the problem is terminated with printouts of all the variables involved. In the present formulation, the 500 sweeps of the field or 1425 points take less than three minutes.
Thom and Anelt [11] give a number of large molecule formulas which may be applied in the same manner as (3-1). Of these formulas only the 20 formula was tried, since the still larger molecule formulas require special tests for off-boundary points. While the 20 formula showed convergence to occur more rapidly in terms of the number of iterations, the more complex formula used more machine time per iteration. Also, with (3-1) the correspondence with the continuity equation is direct in terms of the velocity components, whereas the interpretation of the large molecule formulas in terms of velocities is not meaningful.
The attributes of a convergence test that are of particular interest are:
While the criteria used are empirical in character, they have served their intended purpose in the applications considered.
We note that changes in a given ψi,j between successive sweeps are at first large and then approach zero asymptotically. If we assume that the ψ's will everywhere change in such a manner that the magnitude of the change is always less than the preceding change, then we may apply the criterion
( | ψi,jk+1 - ψi,jk | )max / ψr < 0.0002 (3-2)
where ψr (a positive number) is the reference value prescribed at some point in the field. The number on the right was experimentally determined, a smaller number producing no appreciable change in the results.
We choose
ψr ≈ a2 |ω|
such that for significant values of ω the terms on the right of (3-1) are of the same order of magnitude. To further insure that velocities of significant size are accounted for, a second criterion is applied after (3-2) is satisfied. We assume that the values of ω that are of a magnitude of 1% of the maximum value or greater are significant, and we require for the totality ot such points that
| { ∑i ∑j |ωi,j| } / { a½∑i∑j | ψi+1,j + ψi-1,j + ψi,j+1 + ψi,j-1 -4 ψi,j| } -1 | < 0.0002 (3-3)
The number on the right in (3-3) was again determined through experimentation. Normally (3-2) will be satisfied with fewer sweeps than (3-3). Equation (3-3) is not applicable to initial solutions which satisfy Laplace's equation (see Section 5).
The normal procedure in obtaining the solution is to perform 20 sweeps of the field and then test (3-2). If (3-2) is not satisfied, ten more sweeps are made prior to testing again. When (3-2) is satisfied, then (3-3) is also tested. Further groups of ten sweeps are made until both criteria are satisfied. In the applications thus far considered, all cases in which the criteria were properly satisfied have given a solution which was the appropriate one.
One form of treatment of the inflow-outflow boundaries is to prescribe them to be periodic. In the finite difference calculation, periodicity is accomplished readily by computing left and right boundary points in such a way that they involve the same quantities. Thus in computing the vorticity, only one of the sets of boundary points needs to be calculated using (2-6). The opposite boundary set may then be prescribed equal to these. If the left boundary is the computed boundary, then appropriate values as required by the calculation are taken at points just inside the right boundary. While functional values of the u velocity component are required on end boundaries, no special treatment of (2-2) is needed.
In performing the ψ iteration as discussed in Section 3, it may be considered worthwhile to make the sweeps in a random fashion to reduce any residual effect of asymmetry that might arise from constant sweepig in a single direction. In such a case it has been found that the periodic boundaries may be calculated separately without causing any difficulty in terms of convergence. In this treatment the advanced values of (3-1) that are used may differ on the left and right boundaries, resulting in a slight difference in the final boundary values. If a sufficiently stringent convergence criterion is used, these differences in the right and left values become insignificant. On the other hand, if the convergence criterion is sufficiently stringent, a one-directional sweep will usually give good results without asymmetries.
We consider cases in which the wall is either stationary or moving with constant velocity relative to the fixed reference frame. A stationary wall is treated in an identical manner to obstacles in the flow (see Section 4-C). For a moving wall we use boundary conditions derived by assuming no influencing factors in the flow except those due to the motion of the wall itself. Corrections are then applied where disturbances in the flow affect these wall conditions.
Let the plane y=0 define a wall adjacent to a fluid occupying the upper half space y>0. If the wall is moving in its own plane, fluid motion obeys the equation
∂u/∂t = ν ∂2u/∂y2 (4-1) or if we define ω = - ∂u/∂y (4-2) and take the partial derivative of (4-1) with respect to y, we obtain ∂ω/∂t = ∂ ∂2ω/∂y2 (4-3) Now ∫0∞ ∂ω/∂t dy = ν ∫0∞ ∂2ω/∂y2 or d/dt ∫0∞ω dy = ν |∂ω/∂y|0∞ = - ν (∂ω/∂y)0
since sufficiently far from the origin ω is henceforth zero.
Using (4-2) and then (4-1) we obtain for a constant wall velocity
d/dt ∫0∞ ω dy = (∂u/∂t)0 = 0 (4-4)
or the integral of the vorticity is conserved.
Thus we note that any initial vorticity present at the wall will prescribe the total vorticity that exists in the field of interest for all time. We conclude, then, that, in terms of the finite difference calculation the treatment of such a boundary condition requires that there be no loss of vorticity as a consequence of transport or diffusion at the wall. We are guaranteed that no extraneous vorticity will creep into the system in connection with interior points by the application of conservative difference forms (see Section 9-C). This conservative property must be extended also to the application of (2-6) at wall boundaries.
In regard to the diffusion part of the equation, one can readily prevent the gain or loss of vorticity by prescribing an ω value equal to the boundary value at a hypothetical point one cell outside the region of interest. Since diffusion will not occur between equal values, this prescription will satisfy the requirement.
We note that if we sum the vortices along the positive y axis, we may write
ω0 + ω+1 +ω+2 + ... + ωj = constant (4-5) for any time t, provided ωj+1 = o. But since ωj = (uj-½ - uj+½)/a (4-6) (4-5) requires that u-½ = constant (4-7)
for all time, u-½ is a hypothetical value of the velocity ½ cell distance outside the region of interest. Simple reasoning indicates that u-½=g(u0) only, where u0 is the wall velocity. Indeed experiments show that u-½ = u0 (see Section 7-A).
We may now also prescribe a ψ value, ψ-1 outside the wall, such that
ψ0 - ψ-1 = u0 a
Since the wall must constitute a streamline, ψ=ψ0 everywhere on the wall, then ψ-1 = constant at all points one cell outside the wall. Thus we can define v-1 for use in (2.6). Th!s prescription prevents transport or vorticity through the wall. This argument !s admittedly heuristic since no transport is involved in the problem discussed. We are, however, assuming that the full Eq. (2-6) is being applied to this simplified example.
In regard to transport parallel to the wall, it is reasonable to substitute u0 for the average of the velocities as prescribed by (2-7). Hovever, in Section 7-A the above problem is discussed and comparisons are made with the analytical solution. It is concluded from these experiments that in terms of the diffusion equation alone the results of the calculation imply that in the finite difference calculation the effective position of the wall is ½ cell below the region of interest. We have previously noted that u-½=u0, and for consistency With the diffusion solution we shall always substitute u0 for any value of velocity, called for in the application of (2-6), which requires an x component of velocity outside the region of interest.
With these considerations we write for an upper wall
ωi,j0n+1 = [1/(1+4νδt/a2)] {ωi,j0n-1 + (2δt/a)[(uω)i-½,j0n - (uω)i+½,j0n + (νω)i,j0-½n ]
+ [2νδt/a2] (ωi+1,j0n + ωi-1,j0n + ωi,j0n + ωi,j0-1n - 2 ωi,j0n-1) } (4-8)
where j0 is the y index at the upper wall. The corresponding equation for the lower wall is
ωi,0n+1 = [1/(1+4νδt/a2)] {ωi,0n-1 + (2δt/a)[(uω)i-½,0n - (uω)i+½,0n - (νω)i,j0+½n ]
+ [2νδt/a2] (ωi+1,0n + ωi-1,0n + ωi,+1n + ωi,0n - 2 ωi,0n-1) } (4-9)
Again note that averages are used as obtained by (2-7) through (2-10) with permutations of the indices as required. In the averages the respective wall velocities are subsituted for uj0+½ and u-½.
The treatment at the lower wall is derived from the fact that
u-½ = u0 We note that ψ-1 = ψ0 - u0a But also at the lower wall 2ψ0 = ψ-1 + ψ+1 + a2ω0 Hence we may write ψ0 = ψ+1 -u0a + a2ω0 (4-10)
While (4-10) is to be satisfied it does not insure that the wall will constitute a streamline. The current method used to insure that the walls are streamlines is to use the proper form of (4-10) for each point; then average these values for the full length of the wall, prescribing the average value so obtained as the ψ value for the wall. Further, values so obtained are iterated values; they are obtained prior to each sweep or interior points of the iteration process (see Section 3).
The equations for the upper and lower values are, respectively,
ψj0 - u0a = 1/(i0 +1) ∑i=0i=i0 ( a2ωi,j0 + ψi,j0-1 ) (4-11) ψ0 - u0a = 1/(i0 +1) ∑i=0i=i0 ( a2ωi,0 + ψi,+1 ) (4-12)
where i0 is the x index for the right boundary.
The above conditions on ω and ψ suffice for the simple example of a wall moving adjacent to an infinite region of fluid. But with disturbing influences such as an obstacle in a channel between two such moving walls, constricted flow modifies wall conditions in that vorticity is induced into the wall where constrictions occur. It will be noted that the process used to prescribe ψ permits this variable to adjust to interior conditions but does not guarantee that ψ is consistent with ω at the wall in the local sense or (4-10), for example. The procedure, therefore, is to invert (4-10) after the ψ solution has been obtained and re-evaluate ψ at the wall, to be everyvhere consistent with the ψ solution. This correction applied to the moving wall vorticity values suggests a treatment applicable at obstacles.
In the treatment of obstacles the first requirement is that the periphery of the obstacle must be a streamline of the flow. This means that ψ must be held fixed on obstacles during the iteration process or, in some fashion, the final results must be such that ψ is everywhere the same on the periphery. At first this appears to be a dilemma since in the normal calculation we have vorticities in the field which act as source functions in the iteration of Poisson's equation. Thus if ψ were held fixed on the obstacle, the vorticities that existed there would not be taken into account.
This apparent dilemma turns out to be a blessing since we need not use (2-6) in connection with obstacle-fluid interfaces. We advance the time on ω using (2-6) everywhere except at the obstacle. We then let the obstacle take on the reference ψ value, ψ = ψr, and iterate, holding the obstacle ψ fixed. We proceed then by obtaining the component velocities, and finally the vorticities are brought up to date using (2-4). Since (2-4) and (2-5) are equivalent we note that if we applied the resulting vorticities as source functions in (2-5), the end result would be the same, even though ψ were not held fixed on the obstacle. Thus we see that even though the character of the equations requires a given order of calculation, the fundamental requirement that the interface be a streamline requires changing this order.
The vorticities prescribed on the obstacle in the above-described manner are affected by the transport and diffusion properties of ω in the sense that transport and diffusion may take place between these points and adjacent points in the fluid. Ths occurs through the normal time advancement of ω, using (2-6), at adjacent points.
In general, the integral of the vorticity does not remain a fixed value when an obstacle is in the fluid stream. This is not surprising when it is noted that, unlike with the uncorrected wall conditions considered, there are accelerations occurring in the vicinity or the interface of the obstacle and the fluid.
The application of (2-4) at the obstacle is simple if the required velocities on and interior to the obstacle are prescribed (zero). Note that the obstacle should be at least one cell in thickness and/or height, since otherwise one would have to store the inferred vorticity values separately with regard to the side from which the obstacle is observed. Further, tests would have to be applied in evaluating these vorticities, using (2-4), such that proper velocity values would be used. Velocity values on opposite sides of the obstacle certainly cannot be used for a given computation.
Note that obstacle-fluid interfaces are such that vertical elements pass through u points, and horizontal elements through ν points. This prescription is the natural one for satisfying the property that the normal component of velocity be zero, or that the periphery of the obstacle is a streamline. Some efforts were made with reverse placement of the obstacle interface such that horizontal elements passed through u points and vertical elements passed through ν points. These attempts were not successful. Such a method may possibly be made to work, but indications are that if it is to be successful it will involve an elaborate testing program.
In the examples given in Section 7 the inflow-outflow boundaries are periodic. To more closely represent the characteristics of an idealized wind tunnel, we consider the possibility of prescribing uniform inflow conditions on the left and exrapolated conditions on the right. By uniform inflow we mean thet u=u0, a constant, and ν=0 on the left boundary. Here ω is zero on the inflow boundary while ψ is a linear function of y, both remaining fixed throughout the calculation, ψ is also a fixed constant on the walls. The limitations of space surrounding an obstacle for the current mesh used are such, however, that an inflow condition of this type can only be interpreted to represent a jet, since obstacles cannot be placed sufficiently far downstream.
Perhaps a better approximation, considering the lack of space, would be to let u, at the inflow boundary, adjust to interior conditions of the flow but maintain ν=0 at this boundary. In such a case ω is given by
ω0,j = (u0,j-½ - u0,j+½)/a (4-13) and 2ψ0,j = ψ0,j+1 + ψ0,j-1 + a2 ω0,j (4-14)
Here (4-14) is the treatment given ψ on the left boundary in the iteration process, while (4-13) may be calculated after convergence in the same manner as vorticities on an obstacle.
In regard to the continuative outflow conditions we shall simply indicate a method which has given qualitative results in some test cases.
We take
ψi0,j = 2 ψi0-1,j- ψi0-2,j (4-15)
in the iteration process, ignoring the possible existence of vorticity on the right boundary. Equation (4-15) simply states that the y component of the velocity does not change in the space of the last two cells on the right.
For ω we require that it take on the same value as the mesh point just interior to the boundary. This prescription affects the calculation only in diffusion and transport at the mesh points just inside the right boundary. We assume that the vorticity values are of insignificant size in this region.
We have given the above discussion not for purposes of actual application, but only as a basis for further experimentation with boundary conditions of this type.
In the applications to be discussed. in some detail in Section 7, the initial solutions are all approximate solutions to Laplace's equation, i.e., the initial solutions are obtained using the equation
ψi,j = ( ψi+1,j + ψi-1,j + ψi,j+1 + ψi,j-1 )/4 (5-1)
ψ values on walls and obstacles are held fixed during the iteration process using (5-1). While convergence is tested using (3-2). In the examples discussed, the end boundaries may be considered to be periodic, although some other possible prescriptions would lead to essentially the same results for these initial solutions.
We present the initial solutions separately since a given initial solution forms the basis for a whole series of time-dependent problems Thus irrotational flows take on an aspect of usefulness here even in cases where such solutions are simple enough to be prescribed without any calculation.
As a first example we consider a case where no iteration is required. Let the fluid, occupying the space between horizontal walls, be stationary while the walls move at a velocity u=um0. The initial, non-viscous solution is, in this case, ψ = constant everywhere, i.e., the frictionless motion of the walls does not affect the fluid in any way. The complete description in terms of the variables used is v = 0 everywhere, u=ω=0 at all interior points, u = u0 at the walls, ω=-∞ at the upper wall, and ω=+∞ at the lower wall. The latter two prescriptions come from the definition of ω (1-18). We, however, modify the prescriptions of ω at the wall for the calculation, since ω is given in terms of finite differences. For example, at the lower wall the vorticity is given by
ωi,0 = (ui,-½ - ui,+½)/a (5-2)
Here ui,-½ is a hypothetical velocity defined at a point ½ cell below the region of interest (outside the lower wall). We consider two interpretations of ui,-½. These are
ui,-½ = u0 (5-3) and ui,-½ = 2u0 - ui,+½ (5-4)
To Section 4, however, we obtained the result ui,-½ = constant for all time; hence (5-4) cannot be correct because ui,+½ varies in time, Thus with the prescription (5-3) and
U1,j0+½ = u0 (5-5)
for the upper wall, we may obtain initial vorticity values using (5-2) and
ωi,j0 = (ui,j0-½ - ui,j0+½)/a (5-6)
In Section 7-A, it is shown that the prescriptions (5-3) and (5-5) indeed lead to good agreement with analysis. By contrast {5-4) prevented convergence after the first time step.
The top two configurations of Fig. 3 represent the initial condtions of the problem we have been considering. An obstacle has been placed in the stationary fluid, but it takes no part in the initial solution. All left frames of Fig. 3 are examples of streamline patterns, while the right frames are the corresponding vorticity contours. In the case we are considering there are no streamlines because ψ=constant everywhere. The vorticity contours are observed only at the cells adjacent to the walls, because ω is represented to vary linearly, through the plot program, from its value at the wall to zero at one cell distance from the wall. A continuation of this problem with viscosity is given in Section 7-A, and plots for succeeding times are given in Fig. 4 and Fig. 5.
As a second example of an initial solution we consider a case where the walls and the fluid are in motion but a stationary obstacle is present in the fluid stream. The center frames of Fig. 3 represent the initial solution for this case. The initial ψ solution is first guessed and then iterated to obtain the resultant ψ field. Vorticity is present only at the obstacle-fluid interface. These initial vorticities are indicated in the contour patterns on the center right.
The guessed solution is obtained by prescribing ψ= ψr, the reference value on the obstacle and on a horizontal line through the center of the obstacle. Along some line joining the walls, but not including the obstacle, we may apply (2-2), taking ui,j+½=u0 >, to obtain guessed ψ values for all points on the line. Now ψ's at all points on horizontal lines are set equal except for those on the obstacle. With this guessed solution the iteration program is entered, leading to the final solution after a few iterations.
After the ψ solution has been obtained, Eqs. (2-2) and (2-3) are applied to obtain interior velocity values consistent with the ψ solution. The velocity values in turn are used in (2-4) to obtain the vorticity values. The vorticity values so obtained are essentially zero everywhere except at the upper and lower surfaces of the obstacle, provided that the ψ solution has converged sufficiently. This example of an initial solution forms the basis for the results discussed in Section 7-B.
We next consider an example in which the fluid and the upper wall are in motion but the lower wall is stationary. In such a case we may treat the lower wall as an obstacle, prescribing ui,-½=0. The initial ψ solution is again ψ=constant everywhere. With v=0 and u=u0 the upper wall and at all interior points, we may apply (2-4) to obtain
ωi,0 = (ui,-½ - u0)/a = -u0/a (5-7)
at the lower wall. This solution can form the basis for the study of the boundary layer on a flat plate.
As a final example of an initial solution we simply attach an obstacle to the stationary lower wall of the previous example. we take ψ=ψr on the lower wall and obstacle, establish a guessed solution as described above, and then iterate the ψ solution. Again this is followed by re-evaluating the u and v fields using (2-2) and (2-3). The vorticity values for the initial solution are then obtained using (2-4). The stream function plot for this example is the lower left frame of Fig.3, while the corresponding vorticity contours are displayed in the lower right frame.
To set up the program for calculation, input data, containing parameters which specify boundary placement and types of boundary conditions to be used, are supplied to the computer. An initial program converts these data into a predetermined set ot indicator bits which are stored in an area containing a space for every point of the mesh. The indicator bits are tested throughout the code; thereby the proper formulas applicable at each stage of the calculation for each mesh point are selected.
Constants for the calculation are based on a minimum number of parameters in terms of input data. Basic required parameters are the viscosity, the reference ψ value, velocities of the walls, and an indication as to whether the fluid is to be initially stationary or moving, other parameters such as the print and plot intervals may be supplied or may be made a basic part of the program. The δt used in the calculation is based upon criteria discussed in Section 2-C. Provisions exist for modifying δt downward should conditions arise in which these criteria may be violated. There is no need for upward revision of δt since its size does not strongly influence the speed of calculation (see Section 3).
Provisions exist in the setup part of the program to permit the use of restart data for continuation of a problem. The restart data is obtained in routine fashion when a requested amount of computer time has elapsed.
The diagram designated Simplified Control Block indicates the order of calculation. In Section 5 we discussed the initial solution. This part of the program often requires revision to permit an efficient initial calculation. While the procedure for iterating the ψ solution is the same as that to be discussed in Section 6-D, a good guessed solution should be provided, unless, of course, the initial solution is one that can be immediately prescribed. While the guessedd. solution generally involves only simple calculations as has been discussed in Section 5, we have thus far not succeeded in making the program a general one for a large variety or configurations.
The guessed solution forms the basis from which the iterated solution is obtained. When all the variables have been consistently established (see Section 5), initial solution plots and prints are obtained to provide a permanent record of the initial conditions.
Time is now advanced and ω is calculated for this new time at all mesh points except those on an obstacle-fluid interface. The designated ωi,j calculation is discussed in detail in Section 6-C. The ψ calculation, which uses advanced ω values, is discussed in Section 6-D. Velocity components corresponding to new ψ values are calculated, and then used to obtain corrected vorticity values at walls and new time, vorticity values at obstacles (see Section 6-D).
At this stage the criterion
{(|u| + |v| )max δt} /a < 0.45
may be tested to determine any trend which may indicate that a smaller δt should be used. If such is indicated the calculation can be repeated using a smaller δt; otherwise the fields of storage may be rotated in preparation for a new advancement in time.
The fluid momentum is next calculated for later use in connection with fuctionals of the motion (see Section 7-A).
Test particle coordinates sre modified in increments obtained by using a two-dimensional linear interpolation of the four nearest u velocity values and the four nearest v velocity values.
Finally a test is performed on the time to determine if the calculational results currently obtained are to be recorded. If not, the calculation proceeds to a new time. If the calculational results are to be recorded, the desired functionals may be calculated and included in the plotted and printed results.
The essentials of the vorticity calculation are outlined in the flow diagram ωi,j Calculation. The index, which is initially set and later revised, references storage spaces of the ωi,j's and the corresponding indicator storage spaces. An indicator storage space is tested for each mesh point to determine the treatment to be given the vorticity at the point under consideration. We require every point to have an indicator specification; otherwise the calculation is terminated. Special labels are carried in the indicator storage spaces to discriminate between extreme points in the sense of the lower and upper wall or the left and right boundary points. While the program has a large number of tests to perform which consume calculation time, it is advantageous to have the flexibility to add other conditions that may become of interest. Further, this technique permits an almost unlimited variety of configurations to be prescribed through the input data.
The following formulas are those specified by the flow diagram. (Note that points on obstacles are bypassed in this part of the calculat1on.)
A. ωi,jn+1 = {1/(1+4νδt/a2)} {ωi,jn-1 + (2δt/a) [(uω)i-½,jn - (uω)i+½,jn + (uω)i,j-½n
-(νωi,j+½n ] + (2νδt/a2) (ωi+1,jn + ωi-1,jn + ωi,j+1n + ωi,j-1n -2ωi,jn-1) }
B. ωi,j0n+1 = {1/(1+4νδt/a2)} {ωi,j0n-1 + (2δt/a) [(uω)i-½,j0n - (uω)i+½,j0n + (uω)i,j0-½n]
+ (2νδt/a2) (ωi+1,j0n + ωi-1,j0n + ωi,j0n + ωi,j0-1n -2ωi,j0n-1) }
C. ωi,0n+1 = {1/(1+4νδt/a2)} {ωi,0n-1 + (2δt/a) [(νω)i-½,0n - (νω)i+½,0n - (νω)i,0n
+ (2νδt/a2) (ωi+1,0n + ωi-1,0n + ωi,+1n + ωi,0n -2ωi,jn-1) }
D. ωi0,jn+1 = {1/(1+4νδt/a2)} {ωi0,jn-1 + (2δt/a) [(νω)i0-½,jn - (νω)+½,jn + (νω)i0,j-½n
-(νωi0,j+½n ] + (2νδt/a2) (ω+1,jn + ωi0-1,jn + ωi0,j+1n + ωi0,j-1n -2ωi0,jn-1) }
E. ω0,jn+1 = ωi0,jn+1
In A and D the average values as required in the transport terms are obtained using (2-7) through (2-10) and variations thereof. Similar forms are used for B and C except that the wall velocities are used in place of any u velocities required outside the walls.
We assume that right-hand boundary points are computed first, hence leading to the simplification E.
We have not indicated the treatment which must be given to the four corner points of the mesh. These must in general carry a dual prescription. That is, for example, they may be both periodic boundary points and wall points. Hence the applicable equations are of such a form that both properties are satisfied.
The stream function values are calculated in a simultaneous fashion through the iteration procedure (Section 3). A record is kept of the number of sweeps which are made of the ψ field in the iteration process. If the number of sweeps becomes unreasonably large, the caculation is terminated on the assumption that some error is preventing convergence.
We enter the ψi,j Calculation program by first testing the time to see if we are dealing with the initial solution. In the initial solution the wall >ψ values are held fixed at values established through the guessed solution (see Section 5). For a moving wall succeeding time values are calculated prior to sweeping the field of interior points. All points carry an indicator specification, in the same manner as previously discussed in connection with ω, to identify the ψ formula to be used. Note that obstacle points are bypassed since ψ=ψr, a constant there. Containing walls may be part of an obstacle and in such a case they are treated the same as an obstacle. Moving wall points are bypassed since they have already been calculated.
The indicated formulas are:
F. ψi,j0 = u0a + 1/(i0+1) ∑i=0i0 (ψi,j0-1 + a2ωi,j0) G. ψi,0 = u0a + 1/(i0+1) ∑i=0i0 (ψi,+1 + a2ωi,0) H. ψi,j = (ψi+1,j + ψi-1,j + ψi,j+1 + ψi,j-1 +a2ωi,j) / 4 I. ψi0,j = (ψ+1,j + ψi0-1,j + ψi,j+1 + ψi,j-1 +a2ωi,j) / 4 J. ψ0,j = ψi0,j
Equation F is for the upper wall, and G is for the lower wall. The wall velocities are here specified to be the same (u0).
We have not specified the iterate index in the above formulas. The reader is referred to Section 5 where the iteration method is discussed.
No indicator specification is necessary for the mesh points in connection with the velocity component calculation, since the velocities are everywhere related to ψ in the same way. We have
ui,j+½ = (ψi,j+1 - ψi,j)/a and vi+½,j = (ψi,j - ψi+1,j)/a
everywhere in the field.
The final calculations required to obtain a complete solution for time t are those of the vorticity at obstacles and walls. For obstacles not on the periphery of the mesh we use
ωi,j = (ui,j-½ - ui,j+½ + vi+½,j - vi-½,j) / a
A stationary wall is treated like an obstacle with all velocities taken to be zero except that u value just inside the wall.
If the wall is moving, then an adjustment is made of those vorticities at the wall which had earlier been advanced in time. For example, if the upper wall velocity is u, we correct the vorticity using
ωi,j = (ui,j0-½-u0) / a Similarly, if the lower wall velocity is u0 we use ωi,j = (u0 - u1,+½) / a
In Section 5 the initial solution was given for the simple example of a fluid at rest between horizontal walls which moved with a velocity u = u00 (see Fig. 3). The initial solution was the nonviscous flow solution, and hence the walls slipped without affecting the fluid. Vorticity, however, was analytically implied to exist at the walls and was prescribed through a finite difference mechnism. Now with the advancement of time using (2-6) the initial vorticity at the walls begins to diffuse out into the fluid, the rate of diffusion being, of course, governed by the size of ν prescribed. In Fig. 4 and the upper frame of Fig. 5 this diffusion is observed to progress as indicated by u velocity profiles and vorticity contours at some selected times. In Fig. 5 the fluid is observed to be starting to interact with the obstacle. The lower frame of Fig. 5 is a streamline pattern at a time corresponding to that of the upper frame of the same figure. The u profile is taken at a cross section coincident With the u=0 line as it appears above the vorticity plot. The plot above the ψ pattern is for a ν velocity profile (ν=0 everywhere on the cross section). Printouts are shown above the vorticity plots these are, respectively, t, Δt, ν, m1Δω, m2Δω, Δω, the numbers of extremes in interior ω values, and the number of iterations required to achieve the given solution; similarly above the ψ plot we have t, CD, CL, ψr+m3Δψ, ψr-m4Δψ, Δψ, the number or extremes in interior ψ values, and again the number of iterations for the given solution. Δω and Δψ are the plot intervals for the corresponding plots while m1, m2, m3 and m4, are positive integers. The plot intervals in conjunction with the extreme values which exist in the fields (m1Δω>ωmax, -m2Δω<ωmin, ψr+m3Δψ>ψmax, and ψr-m4Δψ<ψmin) and the number of extremes are all parameters required for the contour plots. CD and CL are the drag and lift coefficients, respectively, of the obstacle.
Our first interest is that the diffusion of vorticity, as calculated by the finite difference method, dose indeed follow closely analytic results. An analytic solution, which may be used for comparison with the results of this calculation, was first given by Stokes [12]. The analytic solution is that of a wall moving in its own plane through a fluid. The vorticity profile taken at points along a line normal to the wall is given by
ω = - du/dt = u0/(√(πνt)e - y2/4νt (7-1)
where u0 is the constant velocity of the wall.
Four successive time plots of this analytic result are given in Fig. 6. We have included the first two calculated profiles (taken for only the lower half of the channel) corresponding to the times prescribed for the analytic results. The discrepancies between calculated and analytic values are reconciled when it is recalled that we had to prescribe a hypothetical value u-½=u0 outside the wall (see Section 4 and Section 5). Effectively, in relation to the analytic result, the wall for the calculated result lies ½ cell below the line y=0. Translating the calculated values by ½ cell makes the results compare favorably with the analytic solution. The indicated points represent the calculated values after the translation. Note that the agreement is good even for very early times when only a few mesh points are involved in the calculation.
The force exerted by the walls on the fluid is proportional to the vorticity at the walls. Using the experimental equation defining the viscosity coefficient μ, namely,
Ft = μ du/dy A
we may write for the tangential force for the two-dimensional case
Ft = ρ ν du/dy L
where A is the area over which the force acts. With unit distance assumed in the z direction, we replace A by L, the length of the wall. Using (1-18) ve may further write
Ft/ρ= -νωL (7-2)
For the example of Figs. 4 and 5 we may compute the total force exerted by both walls in the form
Ft/ρ = νa ∑i=0i0-1 (ωi+½,0 - ωi+½,j0) (7-3)
where use is made of (2-9).
We may also compute the momentum increase of the fluid which has resulted from the forces at the walls. For this purpose we consider the momentum change of a given cell and sum over all cells. In consistency with the central time difference used in the vorticity equation, we here calculate the momentum change over a double time step. Again assuming unit distance in the z direction, we may write
1/ρ (δMx/δt)totaln = (a2/2δt) ∑i=0i0-1 ∑j=0j0-1 (ui+ ½,j+½n+1 - ui+ ½,j+½n-1) (7-4) where ui+ ½,j+½ = (ui,j+½ + ui+1,j+½) / 2 (7-5)
Prior to the fluid interacting with the obstacle, (7-4) and (7-3) should give the same answer since the wall forces are the only forces acting. This has indeed proven to be the case, and we make use of this balance of force and momentum change to evaluate the force exerted on the obstacle. Thus the drag force on the obstacle becomes
FD/ρ = 1/ρ (∂Mx/∂t)total - Ft/ρ (7-6)
Now, using the usual definition of the drag coefficent, we may tabulate Cd. That is,
CD = FD / (½ρuf2d) (7-7)
where d is the vertical dimension of the obstacle. The velocity ur is not well-defined for the problem considered here, since it is the free stream velocity for a case of an obstacle in an infinite fluid. For this preliminary problem we have simply used u0, the wall velocity. This prescription leads to the indicated drag coefficient. We shall discuss the lift coefficient in Section 8.
The results here discussed are those which follow from the second initial solution example given in Fig. 5. We consider only two obstacle shapes, one composed of four cells stacked vertically, the other of only two cells. A variety of viscosity values are considered for each of these obstacles. We refrain from presenting a detailed discussion of quantitative results since these are given by Fromm and Har1ow [13], In that paper, drag coefficients for the two obstacle shapes are given, as are Strouhal numbers and other quantitative features which are compared with experimental results. Our principal aim here is to make comparisons of flow patterns in a variety of representations. Hopefully this will lend an insight into the processes that take place.
We begin with considerations of the familiar symmetric wake for low Reynolds numbers (see Fig. 7). The Reynolds number for the flows is defined as R=u0d/ν where d is the vertical dimension of the obstacle. The cell length is taken as 1/8 unit while the wall velocity is four units, leading to the time scale used in the examples.
In Section 5 we discussed the technique of prescribing the initial vorticity values on the obstacle. Now through the advancement of time, the vorticity present initially only on the obstacle becomes both diffused and transported from the obstacle surface into the fluid. These processes occur solely in connection with calculations made at points inside the fluid since obstacle mesh points are bypassed in this part of the calculation ( see Section 4). We may omit modifying the vorticity on the obstacle through the time-dependent equation because such changes mean nothing in terms of the ψ iteration (ψ = constant on the obstacle). We therefore postpone revision of ω values on the obstacle and make them later conform to the new ψ field. This latter process introduces all additional phenomenon which allows the obstacle vorticity to be replenished in accord with the fluid motions adjacent to the obstacle, Thus three effects are present which modify the vorticity contours as displayed in Fig. 7. These effects are not readily separable. We know that vorticity diffuses more rapidly at lower Reynolds numbers, but the calculations also indicate that vorticity is less rapidly replenished at the obstacle. For example, the integrals of the magnitude of vorticity in the field at the specified time are in ratios of about 20:23:26:36, in the order of Reynolds numbers given in Fig. 7. In all cases the integrated vorticity magnitude grew from that initially present, having been the same for all cases. Indications are that, at lower Reynolds numbers, the diffused vorticity tends to decelerate more the flows in the immediate vicinity of the obstacle, hence slowing the replenishing process. The decelerations of flow in the neighborhood of the obstacle are a consequence of interaction of the mainstream flow and a flow that has its origin in the diffused vorticity. This interpretation of greater decelerations occurring for lower Reynolds numbers is consistent with experimental results which show that greater drags occur in such cases.
In all cases presented in Fig. 7., a symmetric eddy pair exists in the wake at the time for which the plots were made. These do not appear for R = 15 and 30, but they would if a sufficiently small Δψ were used in the stream function plots. The origin of the eddy pair comes about mathematically in that the vorticity transported into the wake region, acting as a source term in (2-5), produces closed contours of streamlines, i.e., produces a circular motion of fluid whose direction of rotation is a function of the sign of ω. At R = 100, we note that there are distinct separated extremes of vorticity in the wake. This effect may be indicative or this flow later forming into a vortex street. The separated extremes of vorticity are of the same sign as the adjacent vorticities on the obstacle surface, having been derived from these originally. We assume that once circular flow has become established, vorticity tends to become trapped into regions of extreme values near the center of fluid rotation, since less transport of vorticity can occur from such points. This would happen more readily, of course, if the diffusion rate were slow. In all cases observed the detached extremes in vorticity diffused avay, forming less peaked extremes. Indications are that there is no regenerative action that would tend to make the extremes grow. The Δω of the plots is 6.4, and hence the detached extreme is of a magnitude of about 19.2. This is only about a quarter of the magnitude of the extreme at the leading corner of the obstacle.
The indicated extremes in vorticity appearing ahead of the obstacle are a calculational instability related to the transport terms of the vorticity equation. At low Reynolds numbers the diffusion rate is sufficient to eliminate these vorticity extremes which alternate in sign. The disturbance becomes increasingly more severe at higher Reynolds numbers, and periodicity permits the disturbance to enter at the right and even work back into the wake region. The reader will note these fluctuations throughout the examples, and we emphasize their existence here so that no attempt is made to read a physical significance into these purely calculational effects.
Recently this instability has been reduced by about a factor of 2 by preventing a contribution to transport which occurs calculationally in a direction counter to the flow. This has been done only in connection with points adjacent to the obstacle since the instability has its origin there. Tests are being performed to see what further progress can be made toward complete elimination of the instability.
In the examples of Fig. 7 the channel-width to obstacle-height rat1o d/h = 1/6. The channel is 24 cells in height.
We proceed to a consideration of a flow sequence (Fig. 8 and Fig. 9) of a vortex street formation for the case R=100. The numerical program has been shown to persist in giving a symmetric flow pattern unless a perturbation is introduced. We introduce a perturbation by augmenting the vorticity on the center line ahead of the obstacle. This is done by adding, at one instant of time, a positive vorticity of about 10% of the maximum vorticity value in the field to three adjacent points. The perturbation does not manifest itself in the wake for some time following its introduct1on, but has in all cases produced a vortex street at Reynolds numbers for which asymmetric flow has been shown to occur experimentally. At lower Reynolds numbers the perturbation diffuses away and produces no visible asymmetry in the wake.
With the boundary conditions used for the problems considered, the fluid momentum continually decreases, as time progresses, due to the drag of the obstacle. Also because of periodicity and the short channel length, the flow at late times becomes less rapid in the central region of the channel. These properties require adjustment of the free stream velocity, and hence affect both the drag coefficient and the Reynolds number of the flow. These effects are discussed in detail by Fromm and Harlow[13]. The effect is manifested in Fig. 8 and Fig. 9 in the overall spreadiog or the streamlines for the general decrease in momentum, and a more noticeable spreading at late times in the flow in the central region due to the periodicity.
In Fig. 10 and Fig. 11 the streakllne sequence for the same problem (R = 100) is recorded. The streaklines are produced by inserting massless test particles on the left boundary at eight points. The particle positions are adjusted at each time step, being influenced by the four nearest u and four nearest ν velocity points. New particles are added at a rate which would allow a maximum of 178 particles, for each line, to cover the field if the velocity were everywhere uniform at the wall velocity. The rate of introduction is adjusted to keep the spacing beween particles nearly the same when the velocity of the field has decreased. This is done by making the rate proportional to the wall velocity modified by the ratio of the current momentum to the initial momentum. If more particles are required on the left after all storage spaces are used, the first ones introduced are removed from the display. This process is observed in the sequence and should not be misconstrued as an effect of the flow. The removal of points from the tails tends to leave the immediate region behind the obstacle deficient in particles, some of which would be caught up in the eddies at early times. The final effect, however, is surprisingly like the corresponding flow representation results obtained experimentally (see Fig. 12 and Fig. 13).
Figures 12 and 13 are included only to show a qualitative comparison because the experimental results obtained by Thom [1] are for a circular cylinder, and the ratios of the obstacle-height to the channel-width differ. In Fig. 12, Thom's results are for d/h=1/10 rather than 1/6. Also, in Fig. 12 the Reynolds number for Thom's result is R = 71, whereas the calculated result is from the R = 100 series for a time t = 3.516, at which time R has dropped to about half its original value. Fig. 12 has the experimental R = 54 with d/h=1/12, whereas the calculated example is taken from the sequence R = 100, d/h = 1/12 (Fig. 32 and Fig. 33) at a time t= 2.461.
Proceeding to a higher Reyno1ds number we include streamline patterns for R=200 Fig. 14 and Fig. 15. Here the mechanism of the street formation is displayed somewhat better than in the R=100 case, since not only are the eddies in more rapid rotation, but we have also reduced the Δψ of the display. Using these resu1ts and the results for R=300 (Fig. 18 and Fig. 19), we consider the process by which a vortex street is initiated.
Some of the distinctive features in the development of the vortex street are as follows:
While the process of the vortex street development can, through the examples, be observed and interpreted. by the reader, we include a preliminary interpretation that has grown out of the results up to this time.
The effect ot the perturbation introduced ahead of the obstacle is to allow, at some later time, an imbalance to occur in the eddy pair. The first indication of this imbalance is that the streamlines show one eddy to have a more expanded interior contour. The imbalance is evident in the R=200 case at t=0.703. For the R=300 case the imbalance is evident at t= 2.285. The big difference in these tmes is that in the R=300 case the perturbation was introduced at a later time.
We assume that the imblance indicates that one eddy involves somewhat more fluid than the other. Some portion of the fluid may, in fact, pass around this eddy and then move along the obstacle to the opposite aide in an S-shaped pattern as is observed at later times. It appears that this eddy interacts more with the obstacle and therefore shifts downstream relative to its partner. With the eddy centers shifted relative to one another, the eddy nearer the obstacle becomes larger while the downstreem eddy diminishes. The downstream eddy diminishes because of an interaction of the reverse flow with the mainstream flow from the opposite side of the obstacle. The upstream eddy grows because its interaction with the obstacle and the downstream eddy forces it further into the mainstream flow,causing its rotational velocity to be agented by the mainstream. Continuation of this process leads either to another shift in the relative positiou of the eddies or to actual shedding.
Note that the S-shaped flow pattern at late stages, prior to shedding, follows the contour of the larger eddy in the reverse flow region and forms the smaller eddy into an elongated region. In this elongated region the reverse flow and the mainstream flow come into very close proximity. This effect leads to the shedding process since the revere and mainstream flows eventually interact to isolate the downstream vortex and form a third vortex near the obstacle.
At this late stage, prior to shedding, the downstream eddy has become extremely week. Its complete return to a position nearest the obstacle hss been blocked by the larger eddy. Once the third eddy appears, the shedding process has begun. The new eddy, being nearer the obstacle, experiences a thrust into the mainstream and hence grows rapidly. The thrust into the mainstream again is a result of interaction with the obstacle in a confined area between as large downstream eddy and the obstacle The downstream eddy again loses rotational velocity because of its interaction with the mainstream, but it is shed prior to becoming as weak as the first eddy. Apparently the shedding occurs more easily since now the relative positions of the centers of rotation and the size of the eddies are such that the entire breadth of the wake is occupied by a predominant eddy.
In this discussion we have said nothing of the vorticity and its influence on what occurs. It is, of course, the mathematical basis for formation ot the eddies, but the motion of the fluid, and hence the eddies, partially determines one distribution of the vorticity. Trapped extremes in vorticity are the mathematical property that provides the calculated rotational character of the eddies as they pass dowstream. The eddy in nearer proximity to the obstacle grows because it is in a field of larger magnitude vorticity, and shifts upstream in its center of rotation becase the larger vorticities are upstream. Finally, the rotation of the eddy is instrumental in detaching an extreme of vorticity from the obstacle, and this trapped vorticity extreme sustains the eddy's rotation as it moves in the direction of the mainstream flow.
There is a somewhat different process in the displayed formation of the vortex street as given for R=300 (Fig. 18 and Fig. 19). This is attributable to the very late introduction of the perturbation, after the eddy pair has developed considerably compared to the R=200 case.
In Fig. 16 and Fig. 17 a streakline sequence is given for the R=200 case. The times are not, in every plot, the same as those for the streamline plots, having been chosen in such a way as to permit the changes to be followed with a minimum mumber of frames. Fig. 20 gives the streakline sequence for R=300.
In Fig. 21 we include late time plots of streamlines for the R=300 case in which a transformation has been made such that the walls are stationary and the obstacle is moving at a uniform velocity to the left. Thus
ψ' = ψ - u0y
This form of display is often found in the literature, since it is useful in obtaining the wave length of the vortex street, its velocity of motion relative to the obstacle, and, hence, the frequency of the oscillations.
Figure 22 and Figure 23 display the streamline sequence for R=1000, whereas Fig. 26 and Fig. 27 show a corresponding sequence tor R=6000. We observe a somewhat more violent oscillation of the vortex street as R is increased, possibly because the trapped vorticity diffuses less rapidly. At late times the fluid in the wake 1s disturbed to a much greater extent than at lower Reyolds numbers, even very near the obstacle.
The disturbed character of the flow at late times is partly due to the calculational instability discussed earlier. Streakline plots for R=1OOO (Fig. 24) and for R=6000 (Fig. 28), however, show oscillations of short wave length, at early times, before the calculational instability affects the flow in the wake region. The streaklines, which are less susceptible to the effects of the calculational instability, also show a very disturbed character of the flow in the wake. In Fig. 29 the experimental results obtained by Lippisch, for a dynamically similar case, show like properties. (See Lippisch [14] for a discussion of the flow representation technique used.)
In Fig. 25 a comparison is made in the plots for R = 1000 and R=6000. It is of interest, calculationally, that a difference in the flows exists even at these very small viscosities. One would expect to reach a stage at which error terms related to transport in the finite difference approximation would be of a magnitude comparable to the diffusion terms rendering any further decrease in ρ undetectable. Also of interest calculationally is that convergence of the iteration process continues to be the same in spite of the more rapid fluctuation in the solution. In general, the convergence rate is more rapid for smaller viscosities, and this continues to be true for these extended cases. It Is believed, however, that the R=6000 case is about as far as the calculational method can be extended in its current form, since here the instability discussed earlier tends to confuse the display patterns.
We proceed to considerations of a smaller obstacle, an example of which was already given in Fig. 13. We include streamline sequences for R=100 (Fig. 30), R=200 (Fig. 31), R=300 (Fig. 40), and R=1000 (Fig. 41) with d/h=1/12. It is somewhat surprising that detailed information is obtainable with an obstacle only two cells in height, yet good quantitative agreement has been obtained in terms of drag coefficients and Strouhal mumbers for these flows (see Fromm and Harlow [13]).
We include more extended streakline sequences for these smaller obstacle examples since little can be observed in the corresponding streamline sequences. The streakline sequences are Fig. 32 and Fig. 33 for R=100, F!gs. and 35 for R= 20O, Figs. 8 and 39 for R = 500, e:od FIga. 2 and 43 for R = 1000, Fig. 34 and Fig. 35 for R=200, Fig. 38 and Fig. 39 for R=300, and Fig. 42 and Fig. 43 for R=1000.
Transformed streamline plots are included for R=200 (Fig. 36) and for R=300 (Fig. 37). These are late time plots after shedding has commenced. The reader is cautioned that times do not always coincide where plots such as Fig. 36 and Fig. 37 are placed together for comparison. While it was possible in most cases to use corresponding times, f1aws in the film sometimes prevented this.
We consider a time-dependent flow based upon the last example of an initial solution given in Fig. 3. We consider a case were R=400 (for the projection height). Again, periodic end boundaries are used, and again, there is no force to sustain the fluid motion except that small amount due to the motion of the upper wall.
At early times the flow becames asymmetric about the obstacle with the formation of an eddy very near the top and on the leeward side of the projection (see Fig. 44 and Fig. 45). This eddy grows, and its center of rotation moves tovard the lower wall and downstream. The reverse motion of the eddy adjacent to the lower wall first induces at the lower wall a vorticity of opposite sign (counter clockwise) to that which was instrumental in producing the eddy itself (see Fig. 49). Shortly following this, the reverse flow transports a portion of the induced vorticity into the re-entrant corner, thus providing the mechanism for the development of a second eddy whose rotation is counter to that of the original eddy.
The growth of the first eddy presumably is not because of more rapid rotation, but simply involves more fluid. The extreme of vorticity associated with the eddy does not grow in magnitude; it spreads slowly owing to diffusion, while the integrated vorticity associated with the eddy increases through the addition of vorticity transported from the obstacle. The vorticity transported from the obstacle is added to that associated with the eddy in the form of isolated extremes that are manifest as ripples near the outer edge of the eddy. These ripples are most visible in the streakline plots of Fig. 46 and Fig. 47, particularly at time t= 0.703. It is uncertain whether this phenomenon is real or simply a consequence of the calculational instability previously discussed. Presumbly the addition of vorticity to the region of the eddy could occur in a continuous manner. Ripples of the same type were involved for the higher Reynolds numbers of the problems discussed in the previous section (see, for example, Fig. 24 and Fig. 28).
The second eddy produced in the leeward re-entrant corner grows (involves more fluid in its rotation) as a consequence of more vorticity being transported from the lower wall into the region of the corner by the reverse flow. Eventually the second eddy's growth causes the flow to form a figure-eight pattern of which one part constitutes a reverse flow in the proximity or the mainstream flow. At this stage an interaction in these flows takes place which prevents the initial eddy from receiving any more vorticity from the obstacle, and its remaining properties simply move downstream with diminishing strength. The isolation of the initial eddy immediately permits the formation of a third eddy as a consequence of the vorticity transported from the top of the projection. By this time, however, because of the nature of the boundary conditions used, the average flow rate has gone down and the new eddy cannot again develop to the degree of the first eddy. Undoubtedly the character of the third eddy would be somewhat different from the first, but the slowing also has some effect here, In any event, the process constitutes the development of a type of vortex street which requires no perturbation to be initiated. The interaction with the lower wall is such as to produce the periodic wake.
It is interesting to note that the second eddy also induces vorticity into the lower wall; a characteristic feature of this type of flow that the lower wall vorticity alternates in sign as it is observed along the leeward side of the obstacle ( see Fig. 49 vorticity plots). Presumably with sufficiently rapid fluid motions of alternating eddies along the lower wall, vorticity which is induced, as in the case of the second eddy, can also be transported into the flow region, causing further formations of eddies.
We include a method which has been used successfully to obtain relative pressures in the flow fields. Since P/ρ satisfies Poison's equation (1-21) with a source term given by (1-22), we may use the same iteration technique as that for ψ to calculate the pressure. We make use of the finite difference form
(P/ρ)i,j = 1/4 [ (P/ρ)i+1,j + (P/ρ)i-1,j + (P/ρ)i,j+1 + (P/ρ)i,j-1 + a2Qi,j] (8-1)
where a2Q is obtained from
a2Qi,j = (ui+1,j2 + ui-1,j2 -2ui,j2 + νi,j+12 + νi,j-12 - 2νi,j2)
+2[(uν)i+½,j+½ - (uν)i-½,j+½ + (uν)i-½,j-½ - (uν)i+½,j-½] (8-2)
In (8-2) the velocity components are averages of the two nearest values.
For example,
ui+1,j2 = ([ui+1,j+½ + ui+1,j-½]/2)2 (8-3)
and
(uν)i+½,j+½ = 1/4(ui+1,j+½ + ui,j+½) (νi+½,j+1 + νi+½,j) (8-4)
Evaluating pressure by the iteration technique is not a very efficient method for obtaining an auxiliary variable, particularly since the iteration process of the basic calculation is the most time-consuming element. There is compensation, however, in that the method leads to a complete description of the pressures, tbroughout the field, with a minimum of complexity in programming. This is in contrast to a system in which pressures are calculated using difference forms of the time-dependent velocity equations. In such a case, calculations would have to be made using the u component equation for calculations made along lines parallel to the walls and the ν component equation for calculations made along lines normal to the walls. If an obstacle were present in the fluid, a number of complicated judgements would have to be made by the computer to permit the calculation to proceed about the obstacle. The possibility of formulating a sufficiently general decision-making program to handle a variety of configurations by such a method seems extremely remote.
In the present program the pressures are calculated only at the plot-time interval rather than at every time step. This is not as good a time-saving device as might be expected, since the old solution becomes a relatively poor guess for the solution to be calculated. Here, however, the rather involved calculation of a2Q is not made as often as for the case where the pressure is always kept current with the other variables.
Convergence criteria analogous to (3-2) and (3-3) may be used.
Here use has been made of only the simpler form
[ |(P/ρ)i,jk+1 - (P/ρ)i,jk|]max / (P/ρ)r < 0.0002 (8-5)
The implied reference value (P/ρ)r in (8-5) has meaning only in regard to this inequality. It is a number of the order of magnitude or a2|Q|, analogous to the prescription for ψr, but it does not specify a fixed value somewhere in the field. With the currently used (periodic) boundary conditions, there is no way of knowing where the pressure contours will lie. This makes it impossible to prescribe a fixed pressure at any point, since if Q has a nonzero value at a point which is prescribed as fixed in pressure during the iteration process, the Q value there will have no effect on the calculation. The procedure has been simply to use the initial solution as a base. Surprisingly this approach does not cause the numbers to drift during the iteration processes.
No program has been set up to establish a guess for the initial solution, which is iterated from the base, P/ρ=0 everywhere. Equation (8-5) is usually satisfied after about 200 to 300 sweeps of the field leading to the initial solution. Subsequent solutions will usually take considerably fewer sweeps since guess solutions are available.
Since the velocity component values are given everywhere in the field and on the boundaries, no difficulty arises in the calculation of Q. In Section 4 we discussed the use of hypothetical velocity values outside the usual boundaries. Here we simply make use of this interpretation in obtaining Q on the boundaries. In connection with walls, the ν component values required in (8-2) that are on or outside the wall are taken to be zero. Similarly the u component values required outside the wall are prescribed to be the wall velocity value.
In regard to periodic boundaries, the treatments of Q and P/ρ are completely analogous to those of ω and ψ, respectively.
At walls and obstacles the treatment of P/ρ must be such that the condition of zero pressure gradient normal to the interface is satisfied. In the former use this is done by simply using the kth iterate value at the wall for any value required outside the wall, in obtaining the (κ + 1st) iterate. For obstacles, however, the treatment becomes much more involved. This is mainly because at obstacle corners the normal gradient is not defined. For nonreentrant corners the current procedure is simply to treat the point as being interior to the fluid. At re-entrant corners, however, a zero normal gradient is taken in both x and y directions.
It is seen in the above that to properly treat the pressure at an obstacle requires an extensive testing program to indicate the relation of the point of interest to its neighbors. We were spared this complexity in the ψ calculation since the obstacle-fluid interface had to be a streamline. In contrast the pressure must take on values as dictated by Q values and neighboring pressures.
We do not here go into detail in the discussion of the computer testing program. It is of primary importance, however, that the implied or real values of P/ρ always be consistent when used in connection with computations at all adjacent points. If this requirement is not met, convergence will not be achieved.
In application of (8-1) and (8-2) it is assumed that the obstacle is at least a cell distance in width and height. In such a case Q may be calculated on the obstacle in the same manner as interior points, provided the interface velocities are prescribed. Also this requirement avoids multiple values of the variables at a given mesh point.
In application of (8-1) and (8-2) it is assumed that the obstacle is at least a cell distance in width and height. In such a case Q may be calculated on the obstacle in the same manner as interior points, provided the interface velocities are prescribed. Also this requirement avoids multiple values of the variables at a given mesh point.
One of the benefits of knowing pressures is that calculation of lift forces on an obstacle can be made. In Section 7-A, we discussed the method of obtaining drag forces from vorticity values. Knowing pressures should further allow separation of drag and lift into components associated with pressure forces and with frictional forces. Such separation hinges on the possibility of obtaining forces directly at the obstacle. Thus, for example, we consider the obstacle problem used to study formation of the vortex street ( see Fig. 48). Pressure force contributing to the drag may be obtained by an equation of the form
FDP = a ( ∑Pfront face - ∑Pback face (8-6)
while frictional forces contributing to drag are give by an equation of the form
FDF = νa ( ∑ωbottom face - ∑ωtop face) (8-7) Thus FD = FDP + FDF (8-8)
Lift force should be separable in a similar way. Unfortunately, attempts to make such calculations have failed. In the drag force case the force, while steady and consistent, is always low. We attribute the error in these results to the coarseness of the mesh which is used to describe the obstacle.
Fortunately, the balance of forces technique discussed in Section 7 gives results for drag coefficients that compare favorably with experiment. Apparently this is so because gross properties of the flow are correctly represented and can be accounted for by measurements of overall momentum change and forces at the walls. We expect then that lift coefficients can likewise be calculated correctly by a balance of forces. Thus if
Fu/ρ = a ∑i=0i0-1 (P/ρ)i+½,j0 (8-9) and Fl/ρ = -a ∑i=0i0-1 (P/ρ)i+½,0 (8-10)
are the forces exerted by the fluid on the upper and lower walls, respectively, then the lift force is given by
FL/ρ = - Fl/ρ - Fu/ρ (8-11) We may then calculate the lift coefficient using CL = FL / (½ρuf2d) (8-12)
If the wake behind a symmetric obstacle becomes unstable, a lift force is experienced by the obstacle. This has been observed to occur for the example of Fig. 48. Indications are that the calculated lift force oscillates with the period of the vortex street. The lift forces are also influenced by the fluid momentum decrease; hence, quantitative studies of comparisons with experiment have not yet been made.
In Fig. 48 we portray a variety of representations for the solution of the R=200, a/n=1/6 case. This was a rerun of the same problem given earlier (Figs. 14-17), differing only in that the introduction of the perturbation was at a slightly different time. The positions of vorticity extremes have been of some interest for purposes of discerning properties of the street which may be compared with analytic models [15]. We note from the vorticity and pressure displays that the extremes in vorticity may be correlated with pressure minimums. We also note that centers of ψ' contours are at a greater distance from the street center line than are extremes in vorticity, in agreement with the discussion given by Shaefer and Eskinazi [15].
In Fig. 48 we note that a pair of low and high pressure is associated with each eddy that is shed. The pressure gradient is greater in the low pressure areas than in the high. With the exception of the first eddy that is shed, the associated low pressure regions are connected by saddle points.
At leading corners of the obstacle there are two low pressure areas, while at the center of the obstacle there is a high pressure region. It is in connection with this property that pressure drag cannot be accurately calculated at the obstacle. In reality the pressure is high over the entire front face of the obstacle with a sharp drop at the corners. The finite difference method cannot resolve this steep gradient with the present coarse mesh representation of the obstacle. It is expected that a calculation in which the obstacle is of the order of ten or more cells in height will permit accurate pressure forces to be obtained.
In Fig. 49 we again show a variety of representations, for the given times, of flow about a projection attached to a fixed wall. Here also it is evident that the pressure minimum may be correlated with the vorticity extreme. In the final pressure display only one low pressure area and one high pressure area are displayed, these being associated with the single eddy that has been shed.
We again note that the ψ' extreme is displaced outward from the fixed wall (or the center line) relative to the vorticity extreme. Because of the weakness of the eddies after the first, only the first is well portrayed on the ψ' plot.
In the vorticity plots we observe, at the earlier of the two times given, the induced vorticity in the lower wall. Some transport of this positive vorticity, by reverse flow, has taken place into the leeward side of the obstacle. At the second time, we note that positive vorticity has formed into a closed contour through the action of the counterclockwise rotating eddy. The presence of an induced negative vorticity is here observable at the lower wall. Of interest in both Fig. 48 and Fig. 49 is the similarity between streakline plots and vorticity plots. This follows from a theorem due to Helmholtz, related to inviscid flow, which states that a vortex is composed of the same fluid particles and therefore moves with the fluid. The similarity in the patterns hence shows to some degree how useful the theorem is for cases of small viscosity.
Throughout the early development of the calculational method the placement of variables relative to the mesh has been the same (Fig. 1). The only differences from the current method were in the form of the vorticity equation. The first form did not make use of a central time difference. This form was
(ωi,jn+1 - ωi,jn) / δt = 1/2a [ (uω)i-1,jn - (uω)i+1,jn + (νω)i,j-1n - (νω)i,j+1n]
+ν/a2[ωi+1,jn + ωi-1,jn+ ωi,j+1n + ωi,j-1n -4ωi,jn] (9-1)
where
ui-1,j = (ui-1,j+½ + ui-1,j-½)/2 (9-2)
and
νi,j-1 = (νi+½,j-1 + νi-½,j-1)/2 (9-3)
Here again permutation of the indices in (9-2) and (9-3) gave the other velocity values required in (9-1).
Considerations of the stability properties of (9-1) led to its rejection, since they tended to limit the range of applicability of the method. We here follow the arguments involved and show those further considerations which led ultimately to the current form.
The stability criterion for (9-1) witbout transport terms is
νδt/a2 < 1/4 (9-4)
This result may be arrived at by a treatment analogous to that given to the corresponding one-dimensional diffusion equation (see, for example, Harlow and Meixner [16]). For the case of transport but no diffusion we assume, as in Section 2, that the velocity components vary only slightly, at some point in the field, from some steady state values. See Eqs. (2-18). Considering only first order variations, we obtain in a fashion analogous to that in Section 2
ωi,jn+1 - ωi,jn = (δt/2a) [ u0(ωi-1,jn - ωi+1,jn ) + ν0(ωi,j-1n - ωi,j+1n)] (9-5)
Substitution of the typical Fourier component solution (2-12) into (9-5) leads to
r = 1 - (iδt/a)(u0 sin k1 + ν0 sin k2) (9-6) The condition or boundness requires that r r‾ ≤ 1 (9-7) where r‾ is the complex conjugate of r. Hence, we must have (δt2 / a2)(u0 sin k1 + ν0 sin k2)2 ≤ 0 (9-8)
for all possible values of k1, and k2. This condition obviously camot be met, and we conclude that (9-1) is unconditionally unstable if ν = 0.
We may infer from the foregoing that (9-1) may be stable for ν>0, in which case the stability must result from the dissipative effects of the viscosity. Also it is expected that greater stability will occur for larger viscosities if the condition (9-4) is not violated. We will see that this is indeed the case.
For analysis of the complete equation we obtain an analogous relation to (9-6) including the viscous effects in the form
r = 1 + (2νδt/a2) (cos k1 + cos k2 -2) - (iδt/a)(u0 sin k1 + ν0 sin k2) (9-9) To simplify the analysis we take K1 = k2 = k; hence (9-7) requies that r r‾ = [1 + 4α(cos κ - 1)]2 + >β2 (1 - cos2κ) ≤ 1 (9-10) where α = νδt / a2 (9-11) and β = ( |u0| + |ν0| ) δt/a (9-12)
We have chosen the form (9-12) for β because if u0 = -ν0, β will be zero without the absolute values. We know that the stability properties in such a case would have to be the same as for u0=ν0; hence, if the steady state components are of opposite sign, then we would have to let κ1=-κ2 and β>4α and the condition for boundness, Eq.(9-7), implies that we must have
(β2 - 4α)2 ≤ 0
But this condition cannot be met except from β2=4α. We hence conclude from the theory that, for stability, we must first require
β ≤ 4α (9-13)
Though no zero slope occurs over the range of values of rr‾ if (9-13) is satisfied, we must still test the end values, cos κ=±1. The case cos κ=1 implies a disturbance of one cell in wave length or an infinite wave length, neither of which can grow. This is indicated by substitution of cos κ= 1 into (9-10), leading to rr‾=1 for all α and β.
For cos κ=-1 we again obtain the condition (9-4) or the stability property without transport.
We note then that the analysis gives the stability properties
α ≤ ¼ and β ≤ 4α (9-14)
for the vorticity difference equation form (9-1). It has, however, been found through experimentation that a stable region exists for which
β2 ≤ 4α (9-15)
We portray the region of stability as indicated analytically by the following diagram and also that added region which has been found to exist experimentally. We account for the added region of stability in that growth can take place at only certain discrete wave lengths, none of which occurs in this region.
We have also included a representation of the stable regions in terms of a δt,ν plot for which we have taken u0=4, ν0=0. With a= 1/8, the stable regions are as portrayed. Check marks are included, in both diagrams, of a selected sampling of experimental test points whose unstable growths are portrayed in Fig. 50.
Without the added region of stability we would surely have to reject the difference form (9-1), since for a given velocity and mesh size small viscosity problems could not be calculated at all. The added region of stability permits calculation at small viscosity, but here δt must be reduced to a smaller and smaller increment. By contrast we have observed that convergence of the iteration procedure is more rapid for smaller ν, indicating that larger time steps would be permissible in order to reduce computation time. It is for this reason that other possible difference forms have been sought.
Prior to consideration of other possible difference forms we consider the manner in which unstable growth is manifested. In Fig. 50 the upper three plots of both columns are vorticity contours, while the two lower frames are streamline plots. The latter correspond to the cases displayed just above them. In the upper left plot we have α=0 or ν=0 and hence no diffusion is taking place. It is for this reason that the perturbation growth is only on the center line where the perturbation was introduced in the same manner as discussed in Section 2. From the plots we cannot discern between positive or negative vorticity, but it is known from numerical listings that the predominant wave length of the growth is four cells. The apparent existence of smaller wave lengths is because of interference between wave lengths of four, five, and six cell length.
With the introduction of viscosity growth of the perturbation is accompanied by an outward diffusion. In the second frame on the left of Fig. 50 we have α=0.1, β=0.7. Here the growth is associated with transport only, and as in the case of ν=0 the disturbance contains components of wave length of four, five, and six cells. The growth involves only the x component, the spread in the y direction resulting from a normal diffusion process. The same is true of the next example in which the diffusion, while free of unstable growth, is more rapid. We see in the last frame (below) the streamline plot for this case.
In he first example of the right-hand column we include a case which is, according to the analysis, in a stable region. We account for the growth in that the perturbation is not sufficiently small for the analysis to be valid in this near proximity to the unstable region. We note that the growth is now a combination of transport and diffusion instabilities in that now a growth is evident in the y component also. This becomes more meaningful through comparison with the next frame where the stability property associated with diffusion is, in fact, violated. Here the y component growth is even more evident, and a similar growth of the x component is noted. With the x component growth, which is here predominantly of the diffusion type, there is still evidence of the presence of a small amount of transport instability.
In the third case we have a pure diffusion instability with a characteristic patchwork pattern in which both x and y component growths are of two cell wave length. Here the diamond-shaped contours are not all displayed because of certain restrictions inherent in the plot program. In the number listings, this diffusion instability is manifested in a given point vorticity value, being surrounded by vorticity values of opposite sign at the four nearest neighbor points.
Finally for contrast with the characteristic appearance of the transport instability we include also the streamline pattern for a purely diffusion instability.
It may be concluded from the foregoing that the stringent stability properties of the first form of the vorticity equation hinge on the fact that the ν= 0 case is completely unstable. If this case could be made stable, then the stability properties would be improved for all values of ν. We therefore look for a difference form of the vorticity equation which will be stable even though the viscosity is zero. This consideration has led. to the use of the central time difference. As has already been noted in Section 2, the central time difference will give neutral stability when ν= 0, provided that (2-21) is satisfied.
If we now, however, consider the two-dimensional diffusion equation with a central time difference, we discover that the diffusion term differences must be modified, since the central time difference form, with the differencing as in (9-1), is unconditionally unstable (see Richtmyer [9], page 94, example 7). Thus if we add a diffusion part that is always unstable to a transport part that has neutral stability, the combined effect can only produce unconditional instability.
We note from example 8 of Richtmyer that the method proposed by Dufort and Frnnkel [8] can give stability with a centered time difference. Use of this technique leads to the current difference method discussed in Section 2.
The other forms of the vorticity equation that have the same stability properties (to first order) as those currently used are
ωi,jn+1 = [1/(1+4νt/a2)] {ωi,jn-1 + δt/a [ (uω)i-1,jn - (uω)i+1.jn + (νω)i,j-1n
-(νω)i,j+1n ] + 2νδt/a2) [ ωi+1,jn + ωi-1,jn + ωi,j+1n + ωi,j-1n -2ωi,jn-1 ] } (9-16)
and
ωi,jn+1 = [1/(1+4νt/a2)] {ωi,jn-1 + δt/a [ui,jn (ωi-1,jn - ωi+1,jn) + νi,jn (ωi,j-1n - ωi,j+1n ) ]
+ 2νδt/a2 [ωi+1,jn + ωi-1,jn + ωi,j+1n+ ωi,j-1n -2 ωi,jn-1 ] } (9-17)
These equations make use of averages of the velocities at defined points as specified by (9-2) and (9-3). There is no concrete evidence that (9-16) and (9-17) are not of equal merit with (2-6). Equation (9-17), however, is a nonconservative form, and experience indicates that such a form can result in errors which are related to the fact that exact cancellation does not occur when the terms are summed over the mesh (see Fromm [17]). This property has not been shown to cause difficulty here, but should be examined prior to making use of (9-17).
1. Thom, A, Proc. Roy. Soc. (London), Ser. A, 141, 651 (1933).
2. Allen, D N, de G, Southwell, R V, Quart J, .Mech. and Applied Math., 8, 129 (1955).
3. Jenson, v G, Proc. Roy. Soc. (London), Ser. A, 249, 346 (1959).
4. Payne, H B, Aeronautical Research Council (London) ARC-R and M 3047 (1956).
5. Blair, A, Metropolis, N, von Neumamn, J, Taube, A H, Tsingou, M, A Study of a Numerical Solution to a Two-Dimensional Hydrodynamical Problem, Los Alamos Scientific Laboratory Report LA-2165 (1957).
6. Thomas, L H, IBM Corporation. Unpublished documents received April 1967.
7. Harlow, F H, Dynamics of Compressible Fluids, Los Alamos Scientific Laboratory Report LA-2412. (April, 1960), p. 120.
8. Dufort, E C, Frankel, S P, Math Tables Aids Comp., 1, 135 (1953).
9. Richtmyer, R D,, Difference Methods for Initial-Value Problems, Interscience Publishers, Inc., New York (1957).
10. Barlow, F H, Stability of Difference Equations. Selected Topics, Los Alamos Scientific Laboratory Report IAXS-2452 (July 1960).
11. Thom, A, Anelt, C J, Field Computations in Engineering and Physics, D. Van Nostrand Company, Ltd., London (1961).
12. Stokes, G C, Cambr. Phil. Trans., 2, 8 (1851}; Math. and Phys. Papers, 3, 1 (1901).
13. Fromm, J E, Harlow, F H, Numerical Solution of the Problem of Vortex Street Development Phys. Fluids, to be issued.
14. Lippisch, A M, Aeronaut, Eng. Rev, 17, No. 2, 24 (1958).
15. Shaefer, J W, Eskinazi, S J, Fluid Mech., 6, 241 (1959).
16. Barlow, F H, Meixner, B D, The Diffusion of Radiation, Los Alamos Scientific Laboratory Report LA-2196 (April 1958).
17. Fromm, J E., Lagrangian Difference Approximations for Fluid Dynamics, Los Alamos Scientific Laboratory Report LA-2535 (March 1961).