Nonperiodic Flow in the Numerical Integration of a Nonlinear Nonperiodic Flow in the Numerical Integration of a Nonlinear Differential Equation of Fluid Dynamics Differential Equation of Fluid Dynamics

Viscous incompressible fluid flow along a flat plate is modeled. by the Navier-Stokes equations with appropriate boundary conditions. A series solution is assumed and a set of three nonlinear ordinary differential equations is derived by truncating the series. The Reynolds number appears in these three equations as a parameter. These equations are solved by numerical integration. We show that these solutions exhibit qualitatively different behavior for different values of the Reynolds number of the fluid. The various modes include an asymptotic approach to a time-independent state, laminar (periodic) flow, and turbulence. We give several computer-generated pictures of the various modes.


I. INTRODUCTION
The Orr-Sommerfeld equation-provides the classic route for the determination of boundary layer stability of fiuid flow along a flat plate.' More recently, chaotic solutions to differential equations have become a tantaliz- ing possibility for the mathematical description of tur- bulent flow.
Here we report on some results of calcu- lations which combine aspects of both approaches.We derive a nonlinear partial differential equation to approxi- mate flow along a flat plate.From this equation we determine a set of three ordinary nonlinear differential equations, and demonstrate that solutions to these equa- tions exhibit both regular and chaotic behavior.We present also some details of this behavior.

II. PROCEDURE A. Derivation of equations
We begin with the Navier-Stokes equation for viscous flow.' Written in terms of the vorticity, these are dn, a'

Bx(
We further assume the flow field Uk to be a basic flow, U& (x2 ), with a perturbation of u ~(x ~,x2, t) and u2(x~, x2, t},i.e., Ul --Ul+u(, U3 --0. Consider Eq. ( 2) with k =3: Using Eq. (3c), U3 -0 and the second term on the right- hand side of Eq. ( 4) is zero.We evaluate 03 in terms of the basic flow and perturbative flow as where Qk is kth component of vorticity, t is time, p is ab- solute viscosity, xI is the lth space coordinate, Uk is the kth component of the flow velocity, 8 is the divergence of the flow velocity, and k is 1,2,3.We assume the fluid incompressible, so that the second term on the right-hand side of (1) is zero; also, 8, the divergence of the flow velocity, is zero and the fourth term on the right-hand side of (1) is zero.We retain Dn, With &3 «om Eq. ( 5), U3 --0, and with some rearrange- ment, Eq. ( 4) is a a a +(Ui+ui) .+u2 co3 at Bx) ax2 V'e -U V'q+U e"+~V4q Recall that Ui -= Ui(x2), so that derivatives of Ui with respect to t or x~are zero.We rearrange the resultant equation: We now introduce a stream function 'P=q'(x&, x2, t) to represent the pertufbative terms, such that BX V'%+ U, ", =0 .

P (10)
We wish to render Eq. ( 10) dimensionless.We make the assumption that fiow is restricted (in the y direction) to a boundary layer of thickness 5. ' We further assume that the fice-stream velocity is Uo.With these assump- tions we obtain the dimensionless equation at ax where R =p5Uo/p is identified as the Reynolds number.
(Note that by neglecting terms noriliriear in 4 and assum- ing U""~v anishes, Eq. ( 11) can be truncated.Assuming that 4 is of the form +=/(y)exp[ia(x ct)] in t-his trun- cated equation then yields the Orr-Sommerfeld equation.
Equation ( 11) is implicit in the standard derivation of the Orr-Sommerfeld equation, but does not explicitly appear in standard references on the subject.) toi -V %. --Using Eqs. (8)in Eq. ( 7) gives a'U a V2%, U a 72%, 1 a& (Sc) 1. Choosing the stream function 4 It is difficult to choose a stream function, %(x, y, t), which preserves the nonlinear character of Eq. ( 11).After To simplify the form of Eq. ( 9) we make the following substitutions: x for xi, y for xz, U for Ui, U~" for a Ui/ax2, U~~~f or a Ui/ax2, 4" for aV/axi, 4" for a+/ax 2. This~simplification gives where A =A (t), 8 =8(t), C =C(t), and k, l are positive real constants.This choice was suggested by a stream function used by Lorenz in his paper "Maximum Simpli- fication of the Dynamic Equations."' Using Eq. (12) in Eq. (11) gives

Choosing the flow fgnction U
The flow function U is a function of y only.The boundary conditions of flow along a flat plate suggest that U should be small for y near zero.Also, U should approach Uo as y approaches the boundary layer thickness, 5. (This latter re- striction can be achieved through the use of a multiplicative constant.) These restrictions still allow considerable free- dom in choosing U. We chose (and justify our choice in the Appendix) Uyy -sin(ly) (14) It is our intention to keep only those terms of Eq. ( 13) that are like the terms appearing in )Ii, Eq. ( 12), " namely, mul- tiples of cos(ly), cos(kx), and sin(ly)sin(kx).Therefore we make the trigonometric substitutions 1sin (kx) for cos (kx) and 1sin (ly) for cos (ly) in Eq. ( 13).We also substitute Fq. ( 14) in Eq. ( 13), noting that U and U""will give rise to terms that will not be kept.The resulting equation is In Eq. ( 15), the terms cos(ly), cos(kx), and sin(ly)sin(kx) are linearly independent, so that the coefficients of these terms must separately equal zero.' This yields B. Numerica1 integration of equations (Ref. 6).
Equations ( 17) are the equations we integrated numerical- ly.
The approximations leading to these equations have been severe, and much of the content of Eq. ( 11) may have been lost thereby.We do not expect Eqs. ( 17) to describe real flow to high precision.However, the nonlinearity and dependence on Reynolds number do remain, and we proceed to numerically integrate these equations, with an eye more toward qualitative than quantitative description of real flow.
To preserve the physical validity of the flow character- ized by Eqs. ( 17), some care must be taken in choosing values for the constants k, 1, and m.From Schlicting's discussion of a turbulent boundary layer we use the empirical result that the minimum wavelength of the dis- In addition, we impose the further restriction that k ~I so that f2&0 in Eqs. ( 17).As a result, analysis of the eigenvalues of Eqs. ( 17) is simplified.Lastly our choice, and justification, of the flow'function U require I «1 and n ~l.
We used double-precision advanced BASIC on the IBM Personal Computer for calculations and graphic results.
We also used double-precision FORTRAN on the Universi- ty of Rhode Island NAS 7000N mainframe computer for calculations and graphic results.
Typical integrations were carried out for 5000 time units, and we estimated the precision of the numerical calculations by studying the effect of reducing the time step At.After a time 5000 units, 8(5000) for bt =1 differed from 8(5000) for b, t= -, 'o, by about 0.1%.turbance be about 65, i.e. , the wavelength of the distur- bance is much larger than the thickness of the boundary layer.In our formulation, k and I play the part of wave numbers [as in the terms cos(ly), cos(kx), and sin(ly)sin(kx)] so that we have the restriction:

III. RESULTS
We observed a wide range of qualitative behavior of the solutions of Eqs. ( 19) corresponding to a range of values of R. For R &15.5914, we observed the solutions A, B,C to oscillate regularly (see Fig. 1), the amplitude of oscilla- tions decayed, and the solutions A, B,C approached a 0. fixed point.For 15.5914&8 &18, the solutions oscillate regularly [see Fig. 2(a)], the amplitude does not decay, and the oscillations appear to form a stable 4-cycle (see the following).For 18 &R &21, the solutions oscillate irregu- larly, the amplitude does not decay, and the oscillations appear to be changing from a 4-cycle to an 8-cycle.For 21&R&23, the solutions oscillate regularly [see Fig.

2(b)
], the amplitude does not decay, and the oscillations appear to form a stable 8-cycle (see the following).For R=23.4, the solutions oscillate regularly [see Fig. 2(c)], the amplitude does not decay, and the oscillations appear to form a stable 16-cycle.For R =40, the solutions oscil- late irregularly [see Fig. 2(d)], the amplitude does not de- cay, but varies over a wide range, and the oscillations ap- pear to be chaotic (see the following).
Our judgment as to whether the oscillations are stable n-cycles or chaotic derives from considerations of two types of information: the trajectory of the solution in A, B,C space, and the Poincare section of that trajectory with a plane, A =const.fixed point for R=10. Figure 3 (b) shows that the trajec- tory for R=16is stable and consists of two lobes.Figure 3(c) shows that the trajectory for R =22 is stable and con- sists of four lobes.Figure 3 (d) shows that the trajectory for R=23.4 is stable and consists of eight lobes.Figure 3(e) shows that the trajectory for R=40 is not stable and consists of many lobes of various sizes.
In the nonmenclature of mappings, a mapping is said to be an "n-cycle" if it takes n applications of the mapping to get back to the original point.For instance, a 4-cycle would consist of four points; point l would map into point 2, point 2 into point 3, point 3 into point 4, and point four back into point l.Thus, a 4-cycle takes four steps to return to the original point.Trajectories such as those traced by these solutions do not have such a simple structure.It is useful to take the intersection of this three-dimensional trajectory with a two-dimensional sur- face.The resulting set of intersection points is known as a Poincare section or Poincare map, and gives information about the periodicity of the trajectory."" For the two-dimensional surface we chose the plane 0.17 A =0.0116.As might be expected, each lobe of a trajec- tory strikes this plane in two points.For 8 =16, the tra- jectory of two lobes strikes the plane in four points; thus, R=16 is a 4-cycle [Fig. 4(a)].For R =22, the trajectory of four lobes strikes the plane in eight points; thus, R =22 is an 8-cycle [Fig.4(b)].For 23.4, the trajectory of eight lobes strikes the plane in sixteen points; thus, R =23.4 is a 16-cycle [Fig. 4(c)].For R =40, the many lobes strike the plane in many points; the trajectory and the intersection points do not repeat; thus, R =40 is nonperiodic, or chaotic [Fig. 4(d)].

IV. CONCLUSIONS
The simple mathematical model of a fluid dynamical system gives rise to equations whose numerical solutions behave in qualitatively different ways.This behavior in- cludes decay to constant values, regular oscillation, arid chaos.The authors believe these different behaviors resemble the flow of fluids under different circumstances, viz., fluid motion dissipating due to viscous forces, lami- nar fluid flow, and turbulent flow.
Moon et a/.have made a similar study of another equation of fluid dynamics, the Ginzburg-Landau equa- tion.They find periodic and chaotic regimes as the con- trol parameter is varied, just as we do.In addition, the transition to turbulence'in their case seems to proceed via the "three-frequency scenario" of Newhouse et al.In other work we will attempt to determine the predicted fluid velocities and their spectra from the functions A, B,C, and from there, to determine which, if any, of the FIG. ' 1. A, B,C, vs time for R =10.Time interval between points shown is 40 units.