Modelling a recirculating density-driven turbulent flow

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS I N FLUIDS, VOL. 6,241-253 (1986) zyxwvu zyxwv zyxwvutsrqp MODELLING A RECIRCULATING DENSITY-DRIVEN TURBULENT FLOW BRUCE A. DEVANTIER Department of Civil Engineering and Mechanics, Southern Illinois University, Carbondale, I L 62901, U.S.A. AND BRUCE E. LAROCK Department of Civil Engineering, University of California, Davis, C A 9561 6, U.S.A SUMMARY zyxwvu A mathematical model of turbulent density-driven flows is presented and is solved numerically. A form of the k--E turbulence model is used to characterize the turbulent transport, and both this non-linear model and a sediment transport equation are coupled with the mean-flow fluid motion equations. A partitioned, NewtonRaphson-based solution scheme is used to effect a solution. The model is applied to the study offlow through a circular secondary sedimentation basin. KEY WORDS Finite-element model Turbulence Density flow INTRODUCTION Fluid motions which are initiated or significantly influenced by gradients in fluid density are called density currents. On a large scale both the earth’s gravitational field and its rotation influence fluid motions when less dense fluid is at a lower potential level than more dense fluid. Fortunately the rotational effect becomes insignificant on smaller scales. Fluid density differences may be caused either by differences in temperature or by the presence of contaminants, such as salt or suspended solids. The prediction of flow patterns in such density-driven flows must often simultaneously cope with the facts that the flow is predominantly turbulent and that the structure of the turbulence is itself influenced by the non-constant density field. In confined regions these effects may lead to the creation of zones where the flow turns back on itself or recirculates. The computation of the behaviour of a recirculating density-driven flow is usually complicated further because the flow is driven by spatially continuous density gradients and is not amenable to modelling with the assumption of the presence of distinct density layers, as can happen when the flow field is infinite. Finally, such flows are dramatic examples of a situation where a small deviation from uniformity (in density) causes not a small but a large behavioural change in related phenomena (the flow field), contrary to one’s usual expectation. This paper considers the modelling of what is essentially a turbulent density current driven by the presence of a distribution of suspended solids of higher density than the carrier fluid (water); the settling of the solids creates spatial gradients in the bulk fluid density (solids plus fluid). A higherorder turbulence closure model, the k--E model, is used to model the turbulent transport effects. The use of simpler models would result in a much simpler mathematical form for the final equation set, zyx zyx zyxwvu zyxwvuts 0271-2091/86/040241-13$06.50 0 1986 by John Wiley & Sons, Ltd. Received 8 February 1985 Revised 10 October 1985 242 zyxwvutsrq zyxwvut 9. A. DEVANTIER A N D 9. E. LAROCK but those models would require a priori knowledge of the nature of the turbulent transport. For strong density currents the turbulent transport is intimately linked to a flow pattern which is determined primarily by the density distribution. Consequently a more general turbulent flow closure model is preferred even though it is computationally more complex. The finite element method is used to solve approximately the resulting set of coupled partial differential equations, owing to the ease with which boundary conditions can be implemented, and also to the ease with which the computational grid spacing can be modified. The model equations are therefore presented directly in integral form. No numerical upwinding is employed in the solution of the equations, so the proper choice of grid spacing is important. A practical application of the model to the study of flow through a circular secondary clarifier, which is a common waste-water treatment process, will also be presented. In this way the capabilities of the model, its operational features and (good and bad) characteristics can be displayed. MODEL EQUATIONS zyxw zyx A set of equations that suitably model a sediment laden, density-affected turbulent flow has recently been developed and presented by DeVantier and Larock.' The model is based on laws of mass, volume and momentum conservation for the transporting fluid (usually water) and the transported, suspended solid particles. The use of local volume averaging and the usual Reynolds decomposition of variables in turbulent flow produces a continuum model that is appropriate for application of multidimensional flows. For ready reference the partial differential equations which make up the model are listed in Appendix I. For steady two-dimensional radial flow the element-level residuals that result from application of the Galerkin finite element method to the mean flow transport equations are F, 1. Mi( &(rU) = jQeNi[r( U g + W aZ ~ + a ar p)_w]dR- F,= F, + = IQe aw + W-aw + ap + qAAVEg aZ aZ zyxwvu zyxwvut Nir( Uar F,=JQeNir(Ug+ + j Q e r ( E ~ + aZ W ~ ) d R 1 Q. - W ~ ) d R + ~ Q ~ r A ( l - aNi h)V,dn aZ aZ rE, (aAaNi -~ +--aAaNi) dR ar ar aZ aZ Following consistency arguments advanced by King2 the density of the bulk fluid (the weighted zyxwvut zyxwvutsr zyxwvutsrqponmlkjihg zyxwvu zyx zyxwvut 243 RECIRCULATING DENSITY-DRIVEN TURBULENT F L O W density of the fluid and the transported sediment) must have a single value for each element. Hence a single average value of A, A A v E , appears in the body force term v ] l \ A v & of the vertical momentum equation (2). The factor v] = (p, - pf)/pf is the normalized density difference between the solid of density p, and fluid of density pf, so the body force term represents only the body force associated with the sediment fraction. The modified pressure P at a point is the difference between the actual pressure and the clearwater hydrostatic pressure at that point. The two co-ordinates are the horizontal radial direction rand the vertical direction z in cylindrical co-ordinates. F,, F,, F , and F A are the element level residuals for the modified pressure P, radial velocity U , vertical velocity W , and sediment volume fraction A, respectively. Integration is performed over an element domain Re having an element boundary re,and I, and I, are the components of the unit outer normal from the boundary. A mixed weighting of the residuals has been chosen by analogy to the usual treatment of the Navier-Stokes equations. Hence N i is the eight-node biquadratic function that is used to interpolate all variables except P and weight the momentum and mass transport equations; M i is the bilinear function that is used to interpolate P and weight the continuity equation. As equation (4)shows, the transport of the sediment is influenced not only by convection by the fluid mean velocity components U and W but also by two parameters characteristic of the sediment motion itself; they are the sediment settling velocity V, and the turbulent diffusion coefficient E,. The other notable transport terms in equations ( 2 ) and ( 3 ) are the turbulent momentum flux or velocity correlation or ‘stress’ terms, the overbarred terms; they are also called Reynolds stresses (per unit mass) and must be suitably represented in terms of other variables to close the equation set. The viscous contributions to the fluid stress are assumed to be negligible in this turbulent flow. The terms in these equations are all non-dimensional, based on a characteristic velocity U,, height h and the clearwater kinematic viscosity v and density pf. Many turbulence models have been proposed; they range from the simple to the very complex. The simpler models assume local equilibrium between turbulence generation and dissipation and neglect history and transport effects. The more advanced models also consider the latter effects which are believed to be significant in the present study. The turbulent stress terms have therefore been modelled via the k--E turbulence model proposed by Rodi3 for density-affected flows. The model is built on the kinematic eddy viscosity constitutive relation, which in cylindrical co-ordinates may be written zyxw zyxwvu - uu = $k au 2~,--, ar U uv = #k - 2 ~ , - - , r .- + - aw - ww = 3k - 2 ~ , - - , aZ IC\ IJI - Here k = (G+ uu ww)/2 is the turbulence kinetic energy per unit mass. The kinematic eddy viscosity v, is expressable by dimensional analysis in terms of k and its rate of dissipation E as V , = C,k2/&, (6) where c, is a model constant. Final closure of the equation set is accomplished by developing transport equations for k and 8. These equations are listed in Appendix I. The Galerkin residuals of the resulting equations are ak ak U-+W--Pr+E-B ar aZ ) dQ+ ,,(akaNi akaNi) r- -----+-dQ J- ok ar ar aZ aZ 244 and zyxwvutsrq zyxwvu zyxwvutsrqpon zyxwvuts zyxwvuts zyxwvu B. A. DEVANTIER A N D B. E. L A R O C K aE a& E U p + W--cc,,-[Pr+(l ar aZ k -cE3)B] The model constants cZ1= 1.44, cE2= 1.92, ce3= 0.8, cc,= 1.3, ok= 1.0, and c, = 0.09 are given by RodL3 In addition, E , = yv, with y = 1.2 is chosen.’ In these equations Pr and B are the production of turbulence kinetic energy by the mean shear rate and by buoyancy and take the forms Pr = v , [ 2( g)2+ (g+ + 2( 2) : g) + 2( r>] and Several of the boundary conditions that are required to complete the model require some explanation. For boundaries where solid walls confine the flow domain, the velocity can be determined once the wall shear stress zw is known. The wall shear is related to the near-wall velocity U , by a logarithmic law for wall-bounded turbulent flows where u* = ( ~ ~ / p is) ”the~ shear velocity, ti is the von Karman constant, and v is the kinematic viscosity. The wall shear stress, which is determined by solving equation (11) for , ,z is then used for the boundary integral forms of shear stress in equations (2) and (3). The boundary of the computational region is displaced a small distance a from the wall to ensure that flow in the computational region is indeed turbulent so that the k--E model is valid. The constant K is 9.0 for hydrodynamically smooth walls. The standard k--E model then assigns wall values to the turbulence parameters as and zyxwvuts zyxwv El, 4 = -. tia Fluid free surface boundaries are treated as planes of symmetry for velocities and the turbulence variables k and E. Along a free surface there is no sediment flux, so + + yv,( gl, El,) A(1 - A)Vs12= 0. A similar relation is applied along wall boundaries, where a fraction the wall is resuspended by the local turbulence; thus yV, (g + 1, 1,) + pA (1 - A) V,1, = 0. p of the sediment reaching RECIRCULATING DENSITY-DRIVEN TURBULENT FLOW EQUATION SOLUTION TECHNIQUES zyxw zy 245 It is not difficult to see that the model governing equations are strongly coupled and non-linear. A Newton-Raphson scheme for solution of the equations is chosen because of its potentially quadratic convergence rate. The method is applied in the form zyxwvuts zyxwvu zyxwvuts zy zyxwvut where F iis a residual and q j is an array of nodal values of the unknowns. These residuals F i are global nodal residuals that are assembled by adding together all of the individual element contributions in the usual finite element way. The full incremental change Aqj for the nodal variables is not necessarily imposed; instead the unknowns are updated via qyew= qoy + aAqj. (17) The choice for a in equation (17) should be made so that the nodal unknowns converge to a solution in some optimal fashion. Following Fletcher: a line search was implemented to ensure a decrease in an objective function which is representative of overall solution convergence. The extreme non-linearity of the equations did not allow the use of more powerful search techniques, so a simple interval halving scheme was employed in conjunction with an objective function defined as the sum of the absolute values of the nodal residuals. A further limitation upon CL was imposed because the variables k, E and A must by definition be non-negative. The limitation that precludes negative values of certain variables is known as a realizability constraint. Thus the realizability constraint was first satisfied before a was determined by optimization, and the upper limit for a during the optimization was defined by the realizability limit on a. The earliest attempts to obtain solutions of the equations applied the Newton-Raphson scheme directly to the full equation set, but these attempts were not successful. Consequently a partitioned matrix approach similar to that proposed by Schamber and Larock' was developed and used. The method sacrifices the convergence rate of the full Newton-Raphson method in return for a more stable convergence by partitioning the equations into two sets and applying the Newton-Raphson scheme separately to each while holding the other set constant. This approach was deemed necessary because the values of the turbulence parameters k and E could not reasonably be estimated without prior knowledge of the flow field. The structure of the flow field is in turn strongly influenced by density differences that are affected by the sediment transport equation. The partitioned matrix approach turned out to be a method that is less sensitive to relatively poor initial estimates for k and E. The best choice for the partitioning of variables was not immediately obvious. Partitioning the variables into the subsets ( U j ,W j ,P j ) and ( k j ,c j , A j ) leads to subsets of equations of roughly equal size and was initially attempted. This grouping was not acceptable because of the strong dependence of P upon A. By operating on these variables in separate subsets, the pressure P (above clearwater hydrostatic pressure) could not react to changes in excess body force caused by new estimates for A, and the solution scheme diverged. The successful partitioning of the variables is therefore ( U j ,W j ,P j , A j ) and ( k j ,cj). These two groups then alternately become the qj array in equations (16) and (17), with one set held constant while the other is operated upon by the Newton-Raphson scheme. The first set of iterations is performed upon the mean flow variables ( U j ,W j ,P i , Aj)while holding the turbulence variables ( k j ,E ~ )constant. This step requires the prescription of an initial distribution of v, without a knowledge of k, and so the 2k/3 contribution to the normal stresses iiii and WW were set to zero. Subsequent sets of iterations on the mean flow variables included the 2k/3 contributions, once the turbulence variables had been operated upon for a set of iterations. zyxwvut zy 246 zyxwvutsrq zyxwvut zyxwvut zyx zyxwv zyxwvuts zyx B. A. DEVANTIER A N D B. E. LAROCK The first sets of iterations on the ( k j ,c j ) subset were computed with the use of modified forms of the Newton coefficient matrix dFi/aqj to avoid divergent results caused by the strong non-linearities in the two governing equations. There are two major ways in which the non-linearities affect the convergence of the Newton estimates: 1. The first cause is the behaviour of the sink term involving E2/k in the E transport equation. This functional dependence causes the Newton estimates to diverge rapidly from the solution whenever the E estimates are incorrect and especially when the E estimates are poor. To avoid this behaviour, the Newton coefficient matrix contributions that are associated with the E equation are formed while assuming for the sink term that &/kis a locally specified constant. 2. The second cause is the k2/&dependence of v,. The turbulent viscosity v, appears throughout both the k and E equations in both the terms containing Pr and the diffusive terms. When v, is assumed to be a locally specified constant in forming the Newton coefficient matrix, the iterates converge more smoothly although more slowly. With these assumptions the Newton matrix contributions for F , are dF,MODI akj =- 1. & rNiNjc,,&Pr +(1 - c E 3 ) BdQ, ] (18) The boundary integrals that appear in F , in equation (8) do not appear either because the boundary conditions cause the integrals to be zero owing to flux considerations or because k and E are specified on the boundary, and the entire nodal residual equation is omitted from the set. Note that no simplification is made for the source term c,,(&/k)Pr because the factor v, in Pr causes the overall term to depend linearly on k. The superscript MOD1 denotes that equations (18) and (19) are the first modified forms of the coefficients. They are used until the E estimates are sufficiently well behaved that the full non-linear nature of the c,,c2/k sink term can be tolerated. The second modified form, MODII, is then employed in the forms and aF ,MOD" - aF,MOD' aEj aEj + & r N i NjcEzidQ. Finally, the full Newton form which considers the non-linearity of vt is used; for example, a F y L - aF,MoD" aEj aEj The use of these modified equation forms will be described further in the next section. MODEL RESULTS Radial flow circular sedimentation basins are used extensively throughout the world for both primary and secondary sewage treatment. In some primary treatment units the solids or sediment 247 RECIRCULATING DENSITY-DRIVEN TURBULENT FLOW zyxwvuts zyxwvutsr zyx zyxwvut concentrations are sufficiently low to allow the sediment transport to be considered independently of the determination of the flow field.5In most practical applications this is not the case, however, and this is particularly true of secondary clarifiers. Indeed the current mathematical model was developed so that it would aid in the understanding of secondary clarifiers. Figure 1 is a radial-section schematic diagram of a centre feed radial flow secondary clarifier. The sediment-bearing fluid enters the centre of the basin, is clarified by gravity sedimentation as it flows outward, and most of the fluid exits over the eMuent weir. The sediment, which is deposited on the basin floor, is periodically removed by a slowly rotating suction apparatus (not shown) which also removes some of the fluid. A circular basin was chosen for a number of reasons. (1) These basins are by far the most common design configuration used in the United States. (2) The flow ought to be more nearly two-dimensional than the flow in rectangular basins; in rectangular basins both corners and side walls encourage the development of secondary currents, and these geometric effects are absent here. (3) The velocities at the point sink model of the eMuent weir are not nearly so large as those at the exit from a rectangular basin, owing to radial deceleration, and thus they should affect the flow field only locally. It is regrettable that no suitable data exist on solids-induced density currents in circular basins,6 but qualitative comparisons can still be made. Currently there is considerable interest in the measurement of the details of flows in settling basins, but practical problems such as probe fouling, the avoidance of sludge removal mechanisms, and the very low velocities themselves make the collection of such data very difficult. The authors are aware of one measurement programme in progress at the University of California, Davis, which seeks to collect mean flow and turbulence data. It is hoped that this paper might encourage additional investigators to obtain detailed flow measurements in operating circular basins. This model was used to determine flow patterns and sediment concentration patterns in such clarifiers. The final computational mesh for the problem is shown in Figure 2. The line r = 0 is z t PHYSICAL DOMAIN Figure 1. Schematic diagram of radial flow clarifier cross-section zyxw 1.0 0.0 1.o 2.0 3.0 4.0 Figure 2. Computational mesh. Vertical and horizontal lengths in basin depths (125 ft, 3.8 m) 248 zyxwvutsrq zyxwvutsrq zyxwvutsrqp zyxwvu zyx zyxwvu zyx B. A. DEVANTIER A N D B. E. LAROCK the axis of symmetry for the basin, and the computational region begins 0.65 basin depths away from r = 0 just outside the baffle. The baffle creates a nearly uniform inlet profile at the edge of the computational region. The lines z = 0.0 and r = 4.0 coincide with wall boundaries, and thus computation begins a small distance inward from the walls. Outflow over the effluent weir at the upper right corner ( r = 4.0, z = 1.0) is locally modelled as potential flow towards a sink in the corner. The strength of the sink then defines the outflow velocities near that corner. In the presumed absence of wind across the basin it is assumed that the vertical cross-section is representative of the basin, and thus that the flow may be modelled as two-dimensional. The grid shown in Figure 2 was the result of a number of grid refinements. The most important of these refinements was above the inlet region near z = 1.0, a free surface. The solution shows that this is a region where a return current takes a sharp downturn, and k and E values change drastically radially here. Owing to an insufficiently fine mesh, the k k iterates locally diverged at a few nodes near the free surface. The divergence remained localized because one or two free surface nodal values of k were being forced toward negative values so strongly that the realizability constraint allowed only very small relaxation factors, and thus the iteration algorithm ground to a halt. The refined grid of Figure 2 contains 260 elements and 853 nodes with 2662 U - W - P - A equations and 1530 k--E equations. The proper initialization of variables (as in any Newton scheme) was very important for this problem. The values of k and E were initialized as decreasing linearly in the r direction. This produces a linearly decreasing v, which is what might be expected as a consequence of the linear deceleration caused by the radially increasing cross-section. Initialization of the velocity components to zero caused no algorithmic difficulties, but the values of A and P must be compatible to avoid immediate divergence of the algorithm. This agreement is required because the gradients in P which drive the flow are three orders of magnitude smaller than most of the values of P itself. Thus the deviation in P from the clearwater hydrostatic pressure must primarily balance the excess density body force VAg but also secondarily provide the driving force for the mean flow via its gradient. With this in mind, a linearly decreasing profile in z was chosen in initializing A; then the equation was integrated to obtain P, and the constant of integration was chosen to agree with the chosen reference value for P. A vector representation of a velocity field is shown in Figures 3 and 4. The inflow velocity to the computational region is 0.396cm/s, and the inlet concentration of solids is 1400mg/l. The low velocity is typical of sedimentation basins. However, the flow is turbulent in spite of the low velocity because the length scale is so large (basin depth is 3.8m). The fluid withdrawal velocity along the bottom of the basin is not noticeable because it is 4.6 per cent of the inflow velocity, and yet this results in the removal of 33 per cent of the inflow from the basin. The density current created by the sediment concentration gradients is clearly evident in the Figure. There is a strong current along the bottom and a return current near the top of the basin. This flow pattern is similar to those measured by Larsen’ in rectangular basins. The return current cannot be discerned in the outer portion of Figure 3 because the velocities are very low in comparison with the bottom current and thus appear as x s. Figure 4 presents a magnified view of the velocity field in the outer region. A stagnation point that divides the return flow and the exit flow is also apparent near the upper right corner. The sediment concentration distribution which primarily drives the flow is shown in Figure 5. Clearly the inlet region is subject to relatively strong density gradients. The 500mg/l contour Velocity scale zyxw zyxwvutsr zyx RECIRCULATING DENSITY-DRIVEN TURBULENT FLOW 249 zyxwv C,, = 1400mg/L Lwy 0.0 Olft/s (003 m/s) -~ ~ - * + / + + * / + + / ~ zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA : : ,“. x L A 4 + 4 4 t + + + 4 A X X x x x x X x x x x 4 x x x x 4 x x x x x x x x zyxw zyxwvutsrq f X x X * + + * x + + x * + + + x x * * * * Figure 3. Plot of predicted velocity vector field (small velocity scale). C , , = 1400mg/l Velocity scale , I 0.0 I I I I 0.05 ft/s (0.016m/s) zyxwvu zyxwv Figure 4. Plot of predicted velocity vector field (expanded velocity scale). 130 200 Contours in mg/L Figure 5. Contour lines of predicted suspended solids concentration in mg/l islands near the end wall are an anomaly that is the result of an instability in the A-field prediction near the corner stagnation point. Steep concentration gradients near the corner cannot be properly resolved by this particular grid. Further mesh refinement was not attempted owing to computer core storage limitations and also because the instability disappeared in subsequent solutions for flow with lower inflow sediment concentrations and accordingly lower underflow current velocities. The behaviour is localized, and the overall effect on concentration prediction is minimal. 250 zyxwvutsrq zyxwvu zyxwvutsr zyxw B. A. DEVANTIER A N D B. E. LAROCK Solutions for the k- and &-fieldswere only obtained after an important modification was made in the transport equations given in equations (7) and (8). When the buoyancy term B was included in these equations, no solutions could be obtained for k and E for an initially specified velocity and sediment concentration field. This occurred because the basin was stably stratified and B became a net sink of k. Its value was so large that it swamped the shear production Pr and drove values of k and E toward negative values. Turner' notes that it is a physical impossibility for ( P r + B ) to act as a net sink. Turner also cites a number of sources which assert that -B/Pr, the flux Richardson number, has some theoretical limit less than unity, at which turbulence collapses. He concludes from these sources that the maximum Richardson number should lie in the range 0.1 to 0.3. Therefore buoyancy damping should have only a secondary effect upon values of k and E, and thus upon v,. Larsen,' as well as others, has consistently found strong turbulence in secondary sedimentation basins in normal operation. Consequently, the present numerical solutions have been obtained with B arbitrarily set to zero, with the thought that such a solution should still capture the basic features of the turbulence. Although the initial solution to the problem was difficult to obtain, the development of solutions to subsequent related problems was substantially easier and less time consuming. In those cases the solution to the first problem could be used as the first estimate to start the iterative solution procedure for a subsequent problem. Solutions for problems involving lower inlet solids concentrations were thus bootstrapped successively from one to the next. Each subsequent solution could be viewed as a solution obtained by continuation in A-space or as an incremental unloading. The inflow concentration of sediment was the primary factor determining the strength of the bottom current, so that each lower inlet concentration solution converged slightly faster than the previous one. A comparison of computation schedules for the initial solution and a subsequent solution (see Appendix 11) reveals that the bootstrap solution required less than one third the number of iterative cycles (14 iterations vs. 52). Both solutions required more than twice as many k - E iterations as U - W - P - A iterations, because the k--E equation set contained the strongest non-linearities. In general, a k--E iteration consumed approximately one third the CPU time of a U - W - P - A iteration owing to the smaller number of equations and simpler boundary conditions. The iteration schedules shown in Appendix IT are only in part determined by the computational algorithm. Within each set of iterations (i.e. U - W-P-A, k--E, MODI, etc.) the line search routine chose c1, but the choice of the number of iterations within a set and the timing of the switches from one k--E form to another were controlled by the operator. The results of each set were examined before subsequent sets of iterations were performed so that convergence rates could be monitored. The logic behind the choice of the iteration schedule was necessarily often intuitive, and several subsets of iterations had to be rerun with different k--E modes or increased iterations to avoid divergence. At this point it would not be a simple task, indeed it may not be possible, to automate the iteration schedule completely. Further study is definitely needed in this area. zyxw zyxwvuts zyxwvutsr SUMMARY AND CONCLUSION A two-dimensional model that captures the dominant features of a density-driven turbulent flow has been presented. The Galerkin finite element method has been applied to this model; the resulting non-linear equation set has been solved by partitioning the equations into two subsets and employing the Newton-Raphson method (and two modified forms thereof). The solution scheme has been aided by use of a line search technique. It has been further modified by the need to satisfy a non-negativity constraint on k , E and A. The application of the model to flow through a secondary sedimentation basin illustrates the features of the solution techniques. RECIRCULATING DENSITY-DRIVEN TURBULENT FLOW zyxw 25 1 From this work it is clear that work remains to be done in several areas. In the model itself the role of the buoyancy production terms in the turbulence closure is not fully understood, and additional research to resolve this behaviour is desirable. Further developments that contribute to the reliability or robustness of the solution technique are desired; implementation of realizability constraints and a line search method in a finite element solution algorithm are also features that have received relatively little attention heretofore. It is felt that these techniques will become more important as more research investigators tackle turbulent flow problems where the distribution of a major driving force for the flow, here the creation of density gradients, is itself a part of the solution (and a determinant of the turbulence properties) and is therefore unknown a priori. zyx ACKNOWLEDGEMENT The major portion of this work was completed with the support of U.S. National Science Foundation Grant CME-7914762, which the authors greatly appreciate. zyxwvutsr zyxwvu zyxwvutsr APPENDIX I zyxwvut The partial differential equations which describe the mean flow model and the turbulence model are listed here in cylindrical co-ordinates. The individual variables are defined and explained in the main text. The mean flow transport equations are a -(rU) ar + Y- aw aZ = 0 (Z r U- (continuity), + W-Z +r + -(ruu) ) i r - - - uu a+ r-(uw) aZ =0 (r momentum), + r]Ag = 0 ( z momentum), (sediment conservation). (27) The k--E turbulence model transport equations are + r(- Pr+E-B)=O E2 E +rc,,--rc,,-[Pr+(l k k ( k transport), (28) -c,JB]=O (E transport). (29) Equations (5),(6), (9) and (10) and the accompanying main text define these turbulence variables further. zy 252 zyxwvu zyxw zyxwvut zyxwv zyxwv zyxwvu B. A. DEVANTIER A N D B. E. LAROCK APPENDIX 11. COMPUTATION SCHEDULES Original solution, inflow concentration = 1400 mg/l Iteration numbers 1, 2 3-6 7-13 14, 15 16, 17 18, 19 20, 21 22, 23 24, 25 26, 27 28-30 31, 32 33,34 35, 36 37, 38 39, 40 41,42 43 44-46 47 48,49 50 51, 52 Mode U - W-P-A v, = 009 U-W-P-A V , = 0.09 -+ 0.00 1 1 k-& MODI k-& MODII k-& Full U-W-P-A k-& MODI k-6 MOD11 k-& Full U - W-P-A k-& MODII k-& Full U - W-P-A k-& MODII k-& Full U - W-P-A k-& Full U - W-P-A k-& Full U-W-P-A k-& Full U - W-P-A k-& Full Relaxation factors 0.6, 1.00 1.0, 1.0, 1.0, 1.0 0.9, 1.0, 0.9, 1.0, 0.9, 0.5, 0.8 1.0, 1.0 1.0, 1.0 1.0, 1.0 1.0, 1.0 0.9, 1.0 0.6, 1.0 1.0, 1.0 1.0, 1.0, 1.0 05, 1.0 1.0, 1.0 1.0, 1.0 1.0, 1.0 1.0, 1.0 1.0, 1.0 zyxwvutsrq zyxwvutsrq 1.o 1.0, 1.0, 1.0 1.o 1.0, 1.0 1.o 1.0, 1.0 Bootstrap solution, inflow concentration = 11 20 mgll Iteration numbers 1 2, 3 4, 5 6 7, 8 9 10, 11 12 13 14 Mode Relaxation factors U - W-P-A MOD11 k-& Full U - W-P-A k-8 Full U - W-P-A k-& Full U-W-P-A k-& Full U - W-P-A 1 .o 1.0, 1.0 1.0, 1.0 1.o 1.0, 1.0 1.o 1.0, 1.0 1 .o 1.o 1.o k-8 zyxw zyxwvu zyxwv zyx zyxwvu zyxwvutsrqp RECIRCULATING DENSITY-DRIVEN TURBULENT FLOW 253 REFERENCES 1. B. A. DeVantier and B. E. Larock, ‘Sediment transport in stratified turbulent flow’, J . Hyd. Engineering, A S C E , 109, (HY12), 1622-1635 (1983). 2. 1. P. King, ‘A linite element model for three dimensional flow’, Report 9150 to U S . Army Corps of Engineers. Vicksburg, Mississippi; RMA Inc., Lafayette, California, 1982. 3. W. Rodi, Turbulence Models and Their Application in Hydraulics, IAHR Publication, Delft, Netherlands, 1980. 4. R. Fletcher, Practical Methods of Optimization: Unconstrained Optimization, vol. I , Wiley, New York, 1980. 5. D. R. Schamber and B. E. Larock, ‘Numerical analysis of flow in sedimentation basins’, J . Hyd. Div., A S C E , 107, (HYS), 575-591 (1981). 6. A. I. Stamou and W. Rodi, ‘Review of experimental studies on sedimentation tanks’, SFB 210/E/2, University of Karlsruhe, F.R. Germany, August 1984. 7. P. Larsen, ‘On the hydraulics of rectangular sedimentation basins: experimental and theoretical studies’, Report N o . 1001, Dept. Water Resources Engr., Lund Inst. Tech., University of Lund, Sweden, 1977. 8. J. S. Turner, Buoyancy effects in Fluids, Cambridge University Press, Cambridge, U.K., 1973, pp. 38-164.