1. INTRODUCTION There exists a class of openchannel flow problems which can be adequately described in the context of depthaveraged, twodimensional mathematical models. Essentially, these models are capable of resolving flow currents in two horizontal directions with the fluid and flow properties assumed to be invariant along the vertical. Such simplified representations of a threedimensional phenomena are justified where turbulent mixing, due to bottom roughness, effectively generates a uniform velocity distribution over the flow depth. Thermal effluents from cooling facilities, contaminant discharges from chemical processing plants, and sediment movement by local scouring are all evidence of a progressive industrialization which has the potential to tax the natural cleansing action of waterways to the point where serious environmental consequences may result. Knowledge of the effect of proposed modifications is vital to the safe development and utilization of waterways. In this study, a fundamental investigation of the circulation phenomena in open channels is performed using a depthaveraged mathematical model. Circulation has been observed to have a primary influence on erosion and accretion processes, heat diffusion, and contaminant dispersion. Therefore, the clarification of the physical mechanisms leading to the generation of circulation in freesurface flow is a much needed research endeavor.
2. LITERATURE REVIEW Prior to the development or efficient, highspeed computer hardware, two alternatives were available to the researcher interested in modeling freesurface flow: (1) analytical solutions; and (2) physical hydraulic models. Despite providing the most accurate results, analytical solutions are extremely rare and apply only in highly simplified situations. Conversely, physical hydraulic models are used widely and quite successfully. Most of the drawbacks associated with physical models me of an economic nature owing to the high cost of construction and the time required to adequately simulate the prototype conditions. Lack of flexibility, however, is probably the largest deterrent to selecting a physical model since a change in configuration usually involves expensive modifications. The increased efficiency of digital computers has triggered a strong research effort in the refinement of mathematical models. While offering almost unlimited flexibility in the simulation of various alternatives, these models have the additional appeal of smaller development and operating costs. Certain specialized problems still remain within the realm of physical modeling, but the use of mathematical models continues to grow. Mathematical models of open channel flow exist at various levels of sophistication. The detailed description of flow phenomena is best accomplished in a threedimensional spatial framework; however, the complexity of a formulation in three dimensions often requires tremendous amounts of computational effort. Where valid, the simplification to a twodimensional representation can offer a considerable reduction in complexity and expense. This involves an integration of the threedimensional equations of fluid dynamics over the flow depth. The twodimensional depthintegrated system of equations is not without its pitfalls: there is a closure problem associated with the effective shear stresses which act tangentially on vertical sides of a fluid element. No rigorous relation between these stresses and the depthaveraged variables is currently available. Hansen (6) is credited for being the first to outline the depthaveraged twodimensional formulation. Later, several other researchers, most notably Leendertse (10, 11), followed Hansen in applying Iwodimensional modeling concepts to the study of estuarine and coastal hydrodynamics. In a detailed analysis of freesurface flow, Leendertse developed a computational model for long period wave propagation in wellmixed estuaries and coastal seas (10). Particular emphasis was given in Leendertse's work to the numerical properties of the twodimensional model in the form of a linear stability and convergence analysis following the von Neumann approach. Kuipers and Vreugdenhil (7) extended Leendertse's 1967 model to the realm of secondary flows. By imposing a steady condition al the open boundaries, they were able to use the unsteady character of Leendertse's model as an iterative technique to approach steady circulating flow under certain specified boundary conditions. A theoretical analysis of the vorticitygenerating mechanisms was performed in order to throw additional light onto the causes of circulating flow in depthaveraged twodimensional models. According to Kuipers and Vreugdenhil, if wind stresses are neglected, vorticity can be created either by the convective terms interacting with converging or diverging flow or through the effective shear stresses. Numerical experiments showed the necessity of including the convective inertia terms in order to model twodimensional circulating flow. The observed circulation in the model is attributed to the combined effect of the convective inertia terms and the artificial viscosity, since no explicit account of the effective stresses was made. In a companion report to that of Kuipers and Vreugdenhil, Flokstra (3) has directed attention to the importance of correctly modeling the effective shear stresses. Without actually resolving the closure problem associated with the modeling of these stresses, Flokstra made a detailed analysis of the relevant physical mechanisms contributing to the generation of circulation. According to Flokstra's vorticity balance, it is theoretically impossible to generate circulating flow without modeling the effective shear stress. In addition, Flokstra's analysis leads to the conclusion that a noslip velocity condition at closed boundaries is essential to the generation of eddy flow patterns. Abbott and Rasmussen (1) verified Kuipers and Vreugdenhil's conclusion that the convective inertia terms are necessary for the generation of circulation. Using physical reasoning, they attributed the circulating patterns to a direct consequence of the resistance effects dominating the inertial effects. Abbott and Rasmussen also concluded that "pseudocirculations," occurring strictly due to the truncation errors of firstorder difference schemes, were possible in depthaveraged twodimensional models. The occurrence of nonlinear instability in twodimensional numerical models was also studied by Flokstra (4). Three approaches were cited to handle this problem: (1) A spatial smoothing process similar to that used by Kuipers and Vreugdenhil; (2) the explicit introduction of an eddy viscosity term in the equations of motion; and (3) the use of a difference scheme affected with numerical viscosity. McGuirk and Rodi (12) developed a depthaveraged velocity and contaminant distribution model of openchannel flow which described the recirculation region immediately downstream of a side discharge into a flowing river. Considering the constant turbulent diffusion coefficient and nonexplicit representations of the turbulent structure too crude for the side jet phenomena, McGuirk and Rodi utilized an extension of Launder and Spalding's (8) twodimensional turbulence model. Lean and Weare (9) tested Flokstra's theoretically based conclusions using adepthaveraged circulation model of flows past a breakwater. The effective stresses were shown to have contributions from shear layer turbulence and turbulence generated at the bed. Criteria is presented to delimit the conditions under which the shear layer turbulence will predominate. An observation of numerical circulation (5) similar to that experienced by Abbott and Rasmussen but caused by a coarse computational grid is reported. At present, many uncertainties exist in the mathematical modeling of depthaveraged twodimensional openchannel flow. Several models which exhibit circulation and numerous theoretical reviews on the factors causing circulation are present in the literature. However, no comprehensive analysis of the various mechanisms leading to circulation has been attempted in conjunction with a numerical model. The objectives of this study are, thus, twofold: (1) the clarification of the physical processes which contribute to the generation of circulation; and (2) the identification of significant factors in the numerical specification of circulating flow problems. 3. GOVERNING EQUATIONS The derivation of the equations governing depthaveraged freesurface flow is accomplished by the successive simplification of the general threedimensional fluid flow equations. i.e., the NavierStokes equations. Four assumptions are basic to the equation set used in this study: (1) water is incompressible; (2) vertical velocities and accelerations are negligible; (3) wind stresses and geostrophic effects are negligible; and (4) average values are sufficient to describe properties which vary over the flow depth. Fluid density effects are not considered in this study; thus, incompressibility is a limitation to the applicability of the model. The second assumption is a consequence of the relatively large magnitude of the gravitational body force which, in flows of shallow depth, simplifies the zcomponent of the motion equation to a statement of hydrostatic pressure distribution in the vertical. For the type of flow under consideration, i.e., openchannel flow, the magnitude of the wind and geostrophic effects is insignificant compared to the driving forces found in the mean flow currents. These two terms can be easily incorporated into the model if desired, and their absence does not detract from the generality of the conclusions of this study. The most significant assumption is that an average value is capable of representing a property which normally varies over the flow depth. Until this assumption is made, the analysis is fundamentally threedimensional. Due to the depthaveraging of the equation set, information on the vertical distribution of velocity is partially lost. Fortunately, the shallow waters found in rivers and wellmixed estuaries usually do not require such detailed information. The resulting depthaveraged equations are:
in which ū, v̄ = depthaveraged velocities; t = time; x, y = coordinate directions; g = gravitational acceleration; η = water elevation (η = h + z_{b}); ρ = fluid density; τ_{bx}, τ_{by} = bottom shear stresses; and T_{xx}, T_{xy}, and T_{yy} = effective shear stresses defined as follows:
in which ν = kinematic viscosity; and u' and v' = turbulent velocity fluctuations. The foregoing equation set, although containing many approximations and simplifications, is not closed. Historically, turbulent flow theory has suffered from an incomplete physical representation of the turbulent momentum transfer (Reynolds stresses), i.e., those stresses due to the correlations of turbulent velocity fluctuations. Depthaveraging the formulation further complicates the problem by creating an additional stress due to the nonuniform velocity distribution in the vertical. These two stresses and the viscous shear stress combine into the terms previously identified as the effective shear stress. Representing the effective shear stresses in terms of the mean flow variables is by far the largest impediment to the proper depthaveraged modeling of circulating flow. Since this problem remains to be solved, the use of empirical parameters and calibration techniques is indicated. In this study, the procedure developed by Kuipers and Vreugdenhil (7) is adopted. This particular technique does not explicitly include the effective shear stresses in the equation set, but introduces them in a velocityaveraging routine which effectively simulates the contribution of the effective stresses. In this routine, an averaging procedure occurs after each set of new dependent variables has been generated, as follows:
in which ū_{j,k}* = spatially averaged ū_{j,k}; v̄_{j,k}* = spatially averaged v̄_{j,k}; α = weighting factor; and j,k = spatial indices. When these substitutions are made in the governing equations, closure terms appear such that
in which ∈ = (Δx)^{2} /(2Δt ); α = Δx = spatial increment; and Δt = temporal increment. The bottom shear stress, like the effective shear stress, has not been rigorously related to flow properties. However, years of experimentation has resulted in the development of several satisfactory empirical resistance equations. Any of the applicable resistance equations can be used to relate the bottom shear stress to the flow velocity, assuming the validity of a steady uniform flow roughness. The Chezy expression is preferred here for simplicity, due to the dimensionless friction factor f_{r} associated with it, as follows:
in which f_{r} = g /C^{ 2} and C = Chezy coefficient. 4. MATHEMATICAL MODEL
A finite difference approximation is used to represent the partial differential equation set in the numerical model. In order to visualize the computational structure, the reader must imagine levels of horizontal xy grids layered in the vertical time dimension. Conceptually, the four variables, u, v, η, and z_{b} should all be defined at every node location. However practical limitations in the computational procedure make it more convenient to define a separate grid system for each of these variables. These four grid systems are staggered in space in a form originally due to Platzman (13), es shown in Fig. 1. Normally, it is not necessary for Δx to be equal to Δy; however, the representation of the effective shear stresses used in this model does hinge on this assumption.
Among the several types of finite difference schemes available, the central difference approximations provide secondorder accuracy. Generally speaking, spatial derivatives can always be expressed in a central difference format. However, unless iterations are performed, temporal derivatives cannot be represented by central difference schemes because two entire time levels of unknowns would be present. The additional computational effort required by an iterative formulation is normally not warranted. Consequently, a less accurate but more expedient firstorder difference scheme is preferred for the temporal derivatives. Equations which are nonlinear with respect to the unknown variables present formidable difficulties for an efficient numerical solution. Although the governing equations contain nonlinear terms, a judicious specification of known and unknown values in the finite difference equations results in a linear representation of the unknown variables. This equation set is then solved by an efficient matrix inversion algorithm. The computational procedure used in this model is a multioperational solution mode based on the division of each time step Δt into two stages of halftime step each. Leendertse (10) modified the wellknown "alternatingdirection implicit" or ADI method, by including two explicit schemes in such a way that each stage comtained an implicit scheme followed by an explicit scheme. The advantage of the ADI method, in addition to those attributable to implicit schemes, lies in the solution procedure which solves the xmomentum equation separately from the ymomentum equation, permitting the twodimensional problem to be solved as a sequence of two onedimensional problems. After each implicit step, a single dependent variable remains unknown and can be directly solved for by an explicit method. The numerical model is based on the set of governing equations derived earlier. Each of the three equations has a difference scheme centered about a unique location on the grid system. The xmomentum equation is referenced to the node occupied by ū_{j,k} while the ymomentum equation is centered about node location v̄_{j,k}; η_{j,k} is the reference location for the continuity equation (for simplicity, overbars are omitted hereafter). In the mathematical model used in this study. the finite difference analogs of the various equations are, for the first stage:
Similar discretization techniques are applicable for the second stage: (1) ymomentum and continuity (implicit): and (2) xmomentum (explicit). See Ref. 14 for additional details. Two boundary types can be specified in the numerical model: closed and open boundaries. At closed boundaries, the most convenient specification is the condition of zero mass flux (i.e., zero velocity) in a direction perpendicular to the boundary. At open boundaries, either mean velocity or watersurface elevation may be specified, depending on the modeling needs. Spatial central finite differences such as those used on this model require information which lies outside the boundaries of the computational model. Although a zero velocity tangential to the wall is a realistic assumption (the noslip velocity condition), a zero water level at the boundary can be grossly inaccurate and may lead to numerical stability problems. A satisfactory alternative appears to be the relocation of interior values. A simple relocation technique was chosen in which exterior values were defined to be equal to those interior values adjacent to the boundaries. This is tantamount to a perfect slip condition, a reasonable assumption in turbulent flows in which the viscous effects have little influence on the horizontal velocity distribution. A pervasive problem in twodimensional mathematical modeling is the lack of adequate theoretical numerical stability criteria. Linear stability theory classifies the multioperational procedure as unconditionally stable. However, experience (15) has shown that it is only weakly stable. 5. RESULTS The objective of the numerical experimentation is the clarification of the circulation phenomena found in freesurface flow. Two hypothetical configurations are tested under a gamut of initial and boundary conditions. In addition to a sensitivity analysis, experiments are performed in which more than one element is varied from the baseline in order to determine the interaction between several components of the problem. The bulk of the testing program consists of experiments performed on a channelpool configuration. Channel dimensions are 4 m in width by 30 m in length, while the pool is a rectangle 14 m wide and 15 m long. The first series of tests are performed with a 0.5 m/s velocity specified in the channel as an initial condition and at the upstream end as a steady boundary condition. A water depth of 2.5 m serves as both the initial condition and the steady downstream boundary condition. There is no velocity in the pool initially. Closed boundaries are represented by zero velocities perpendicular to walls. Without a bed slope, the water in the system is driven solely by the upstream entrance velocity. Numerical parameters used in the baseline run are weighting factor α = 0.1, friction factor f_{r} = 0.0045, spatial increment Δx = Δy = 1 m, and temporal increment Δt = 1 sec. Figure 2 is a plot of the baseline velocity vectors after 200 sec of simulation. A welldeveloped circulation appears in the pool area with little or no divergence from the mainstream flow. Water levels in the pool range from 2.50 m to 2.52 m. Entrance levels. however, fall continuously from 2.50 m to 2.44 m.
The necessity of modeling the effective shear stresses where physical circulation is present has been explained theoretically by Flokstra (3). In the present model, the representation of the effective shear stresses is not explicit in the discretized equation set; rather, the action of the effective shear stresses is simulated by a velocityaveraging interphase which employs a weighting factor α to control the magnitude of the effect. The first experiment performed on the channelpool configuration omits the effective stress modeling by setting α = 0. Transfer of turbulent momentum from the mainstream to the pool is insignificant in this case. With no flow divergence into the pool, the initial conditions are virtually preserved. The convective inertia terms, u ∂u/∂x, v ∂u/∂y, u ∂v/∂x, and v ∂v/∂y, are often neglected in numerical modeling because the nonlinearity of these terms introduces additional difficulties. Several authors have commented on the necessary presence of the convective inertia terms in the equation set when secondary flow is to be resolved. To test this conclusion, all four convective inertia terms were set to zero. Figure 3 shows a plot at 100 sec of simulation time. Divergence of the flow from the channel into the pool is strong however, no circulation sets up. As water in the channel reaches the upstream end of the pool, flow immediately is directed into the pool along the boundaries, preventing circulation from occurring.
In the literature, the role of friction
in circulation is not clearly defined. To shed more light on fhis subject,
an experiment is performed wherein the friction terms,
The effect of an increase in flow depth
has been documented by Bengtsson (2) following experiments performed on a lake model.
His conclusion was that the influence of the horizontal turbulence terms was reduced with increasing depth.
In this study, the horizontal dispersion of momentum is simulated by the effective stress closure terms,
The earlier experiments
with friction appear to dismiss the importance attached to it by some authors. However, these tests were
all performed on problems of relatively small scale. To test the effect of friction in a large smile problem,
the size of the configuration was increased by setting
In earlier tests, a perfectslip velocity condition
had been implemented at all closed boundaries with satisfactory results. To test the validity of this
assumption, an experiment is performed in which the flow is subject to a noslip velocity condition as
closed physical boundaries. This condition can be approximated by setting all velocities located
outside of the physical boundaries to zero. For this experiment, a bed slope S_{o} = 0.0005 is introduced
into the channelpool configuration. Open boundary conditions at the channel end points are water levels
which correspond to the initial bed slope condition. Closed boundaries remain to be defined by zero
perpendicular velocities. All other parameters are the same as in the velocityspecified baseline testing,
i.e., α = 0.1, Δx = Δy = 1 m, Secondary flow in sudden expansions is of a slightly different nature than the channelpool system treated earlier. This is due in part to the bending of currents into the expansion and the increased exposure of secondary currents to mainstream effects. Abbott and Rasmussen (1) developed a channel expansion model from which results and conclusions were presented. An attempt is made here to verify these conclusions using a similar geometry. In this testing series, the configuration consists of an entrance channel 9 m wide and 7.5 m long while the expanded channel is 17 m wide and 22.5 m long. A fixed bed slope of S_{o} = 0.0005 in the direction of the entrance channel flow in specified for the entire configuration. The initial condition is a water surface slope parallel to the bed at a depth of 2.5 m. All velocities are set to zero at the beginning of the simulation. Open boundary conditions at both upstream and downstream ends are water levels which match the initial conditions. Closed boundaries are zero velocities perpendicular to the walls. Friction in this model is governed by a dimensionless friction factor f_{r} = 0.0045. The weighting factor is set to 0.1 while the space increment is 1 m and the time increment is 1 sec. Figure 4 is a plot of the channel expansion baseline results after 100 sec of simulation. A welldeveloped circulation forms a zone of separation in the corner of the expansion. Entrance velocities range from 0.72 m/s at the bottom wall to 0.77 m/s at the top wall while exit velocities vary from 0.44 m/s to 0.30 m/s, bottom to top. Water levels throughout the simulation are continuous and stable. Slightly lower elevations are found within the circulation pattern.
For the most part, the testing of the channel expansion produced results consistent with those of the channelpool configuration. A discrepancy was found, however, when the slopespecified models were tested without modeling the effective stress. In the velocityspecified channelpool model, the absence of the effective stresses resulted in no circulation being produced. Contrary to this result, both configurations, i.e., the channel expansion and the channelpool system, are able to generate a spiraling secondary flow without the presence of the effective stresses if a bed slope is the driving force. The question arose that the secondary currents which appear without the benefit of effective stress modeling may be of a strictly numerical nature. Theoretically, numerical effects should be minimized as the discretization of the problem domain is made exceedingly small. To this end, both space and time increments in the expansion model were reduced although the problem size remained constant. The values of these parameters for this test are Δx = Δy = 0.5 m and Δt = 0.5 sec. A horizontal water surface without velocity is specified as the initial condition. As the simulation begins, the downstream water level is lowered in 40 time steps to a depth of 2.5 m, the same depth as the upstream boundary condition. Thus, a line drawn through the and point water levels will he parallel to the bed slope. The weighting factor for this experiment is zero while the nondimensional friction factor is 0.0095. Figure 5 is a plot of the velocity vectors after 150 sec of simulation. A strong spiraling "circulation," with a character similar to those observed in earlier slopespecified experiments, sets up in the expansion corner.
6. EVALUATION The testing program for this study is fairly extensive; therefore, only a limited discussion can be presented here. However. the analysis and evaluation of the numerical experiments is based on all the tests performed. For a more detailed description, the reader in referred to Ref. l4. The use of specified water levels at the open boundaries proved to be a better boundary condition than the velocity specified in the early testing series. When a difference in water levels is the driving force for flow through the configuration, both velocity and water surface are steady and continuous, with few, if any, anomalies. This is in contrast to the velocityspecified model which requires a special initial condition before simulation can begin. In addition, oscillating water levels and velocities plague any of the velocityspecified simulations.
Specific differences in model behavior due
to a change in boundary conditions are found only in the testing of α, the
weighting factor used in the velocityaveraging routine. In the velocityspecified
model, circulation did not occur without the presence of the closure
terms,
While acknowledging the "noslip" velocity condition in nature, the testing of this model has been performed primarily with a "perfect slip" boundary specification. This is a reasonable assumption in models where the spatial resolution across a channel width is not very fine. In the experiment involving the specification of zero velocities outside the boundary geometry, overwhelming resistance effects resulted. Obviously, any consideration of a noslip boundary condition must be coordinated with an accompanying reduction in the bed resistance effects. Results from the testing program are reasonable and encouraging, despite the presence of two flaws traceable to the representation of the effective stresses. First, there is yet to be found a physical basis by which to choose the appropriate weighting factor, α, on the velocity averaging routine. There are instances in the literature where physical processes have been successfully replaced by numerical techniques in a general manner. Although the selection of the weighting factor is presently a manageable calibration task, a physical link to the turbulence process has yet to be advanced. The second drawback in the effective stress representation is somewhat more serious than the first. Much of the success is fulfilling the required profile of traits by the closure terms is due to the effect of numerical viscosity. As the weighting factor is increased, the model reacts as if the fluid is becoming more viscous, thus increasing the exchange of lateral momentum and increasing viscous damping. Essentially, viscosity is used to model a turbulence effect. Therefore, care must be exercised when using the velocity averaging routine. A balance must be struck between the simulation of the effective stresses and the associated change in fluid properties. This study has experimentally verified Flokstra's conclusion that true circulation, i.e., a flow pattern possessing a separation zone with closed circular streamlines, requires the modeling of the effective shear stresses. An orderofmagnitude analysis performed on the various terms in the equation of motion confirms the significance of the effective stresses in all instances where circulation occurs. Circulation requires a continuous exchange of turbulent energy across the shear layer to be maintained. Withdrawal of the effective stresses after the steady flow pattern has set up dissolves the circulation and ultimately leads to instability. The orderofmagnitude analysis of terms found in the momentum equation shows that convective inertia is significant in the channel and shear layer for all instances where circulation occurs. The presence of inertia allows the flow to retain a uniform structure for a distance beyond the location of a configuration change, enabling the development of a separation zone adjacent to the free extension of uniform flow. Ultimately, this separation of flow develops into a local zone of circulating flow. Merely including the convective inertia in the mathematical model will not ensure the generation of secondary flow. Even with the additional stipulation that the effective stresses and all other quantities be precisely described, circulation may yet be inhibited by frictional effects. In every case where secondary currents do not occur, the resistance term is significant and larger than the convective inertia terms. This study has found that bottom friction is the largest deterrent to the existence of circulating flow. A competition seems to exist between convective inertia and bed resistance forces. The fact that convective inertia terms contain spatial gradients of velocity makes them particularly sensitive to scale effects. In small scale problems, lateral differences in velocity occur over relatively small distances, creating large velocity gradients. The large magnitudes of the convective inertia terms render the bed resistance effects negligible. No reasonable friction factor can affect this mechanism. Conversely, largescale problems reduce the magnitude of both convective inertia and effective stresses to such a degree that the friction terms are promoted to the point where resistance effects totally inhibit the generation of secondary flow. This model displays at least two different instability mechanisms: Nonlinear and Courant. Nonlinear instability is theoretically described as the inability to resolve energy at scales smaller than twice the spatial grid increment. Physically, energy is cascaded to smaller and smaller scales by the action of the nonlinear terms in the governing equations. At the smallest scales, energy in dissipated by viscous effects. The discrete formulation of the numerical model interrupts the energy cascade at the resolution of the grid; thus, energy accumulates at this scale and eventually spoils the computations. Nonlinear instability is characterized by gradually developing water surface discontinuities accompanied by a similar behavior in the computed velocities. The process is usually slow, and can take as many as 100 time steps to develop. In this numerical model, nonlinear instability requires that the velocity averaging technique be present in all simulations although the weighting factor need not be large. Most runs were stable with α = 0.1. Physically speaking, the use of the velocityaveraging routine introduces a stronger viscous dissipation mechanism into the fluid. Energy no longer must be transferred to scales smaller than the grid resolution for viscous effects to act. Courant instability occurs when the ratio of physical celerity to numerical celerity (defined as the Courant number) exceeds a characteristic value. The physical celerity for this model is the maximum channel velocity, U_{max}, while the numerical celerity is defined as Δx/Δt = Δy/Δt. Thus, the criterion is of the form (U_{max} Δt /Δx) ≤ ε, in which ε is the characteristic limit. The two explicit stages contained in the multioperational procedure are presumably responsible for the Courant stability condition found in this model. Courant instability is characteristically a very rapid process. The simulation proceeds without noticeable difficulty until a sudden discontinuity appears in one of the dependent variables. Within a few time steps the entire calculation is spoiled. In this study the limiting value of the Courant number is dependent upon the weighting factor used in the velocity averaging routine. This is to be expected since large weighting factors can smooth instabilities that would otherwise occur if smaller weighting factors were used. For α = 0.1, the Courant number in this model must be less than or equal to 0.5. 7. CONCLUSIONS AND RECOMMENDATIONS Major conclusions from this study are as follows:
The following recommendations are offered for future research:
ACKNOWLEDGEMENTS This study was made possible by a National Science Foundation Grant No. CME7805458. The support of this institution is gratefully acknowledged. APENDIX I.  REFERENCES
APENDIX II. NOTATION The following symbols are used in this paper:

180422 21:50 
Documents in Portable Document Format (PDF) require Adobe Acrobat Reader 5.0 or higher to view; download Adobe Acrobat Reader. 