Analysis of Laboratory Scale Aerated Basin Usi Computational Techniques By Robert Wayne Jenkinson B.A.Sc. Civi l Engineering, University of Waterloo, 1997 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R OF APPLIED SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES DEPARTMENT OF CIVIL ENGINEERING We accept this thesis as conforming to the required standard The University of British Columbia March, 2000 © Robert Wayne Jenkinson, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of ^vi- ^^'y/^e^^^ The University of British Columbia Vancouver, Canada Date DE-6 (2/88) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page ii Abstract This study endeavored to apply computational fluid dynamics (CFD) techniques to aerated stabilization basins (ASBs) on a laboratory scale. A laboratory study was set up to assess the hydraulic impacts of a small-scale mechanical surface aerator in a still-water tank. A small surface mechanical aerator was built that would circulate from 0.5 to 4 L/sec. The aerator was employed in a tank and the velocity patterns imported by the aerator were measured using an acoustic Doppler velocimeter (ADV). The velocity profiles were converted into dispersion parameters for use in a CFD model. A local turbulent time scale method was used to determine an approximate dispersion value as a function of distance from the aerator. The analysis examined different aerator power settings under otherwise identical conditions. Three successful runs were completed. The chaotic patterns within the tank induced by the aerator, produced turbulent structures that were not easily resolved using a time averaged statistical turbulence approach. Tracer studies were performed in the same basin, modified to provide a flow-through environment. The studies were performed with and without aeration and at different flow-through rates. Tracer study results were examined and some of the typical parameters used in the literature were extracted for comparative analysis. The trends were surprising as the dimensionless axial dispersion number was inversely related to the degree of mixing applied. A two-dimensional computational fluid dynamics (CFD) model was employed to model the fluid flow within the basin without aeration using a commercial package, FLUENT. The modeling results were compared with depth-averaged A D V velocity measurements taken within the basin. Many CFD modeling parameters were varied, including grid size and structure, turbulence modeling and boundary condition modeling. The model compared well with the A D V data. The results from the FLUENT modeling were applied to a second CFD model that solved the solute transport equations over a structured grid. The solute transport modeling only examined the transport equations with uniform dispersion over the basin and, although it was the original intent, no attempts were made to model aeration at this stage. The overall shape of the modeled tracer curves compared favourably to the tracer study data obtained in the lab. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page iii Table of Contents A B S T R A C T . . . . . II L I S T O F F I G U R E S V L I S T O F T A B L E S V I L I S T O F S Y M B O L S A N D A B B R E V I A T I O N S V I I 1 I N T R O D U C T I O N 1 1.1 SCOPE OF STUDY 2 1.2 STRUCTURE OF THESIS 2 1.3 ACKNOWLEDGEMENTS 3 2 L I T E R A T U R E R E V I E W 4 2.1 CHEMICAL REACTOR HYDRAULICS 5 2.2 APPLIED POND HYDRAULICS 11 2.3 TRACER STUDIES ON TREATMENT PONDS 20 2.4 COMPUTATIONAL FLUID DYNAMICS IN REACTOR HYDRAULICS AND ENVIRONMENTAL ENGINEERING 25 2.5 CONCLUSIONS 28 3 R E S E A R C H O B J E C T I V E S 30 4 A E R A T O R S T U D I E S 31 4.1 AERATOR AND BASIN CONSTRUCTION 31 4.2 VELOCITY MEASUREMENT INSTRUMENT 35 4.3 EXPERIMENTAL METHOD 38 4.3.1 Aerator Calibration 38 4.4 STILL-BASIN AERATOR TESTING 41 4.4.1 Determination of Sampling Time 43 4.4.2 Profile Data Analysis 44 4.4.3 Profile Results 46 4.5 DEPTH-INTEGRATED TRENDS 50 4.6 CONCLUSIONS 53 4.7 RECOMMENDATIONS 54 5 T R A C E R S T U D I E S 55 5.1 BASIN MODIFICATION 55 5.2 RHODAMINE CALIBRATION C URVE 56 5.2.1 Optimum Wavelength 57 5.2.2 Concentration Calibration Curve 57 5.3 TRACER ADDITION AND SAMPLING ROUTINE 58 5.4 TRACER DATA INTERPRETATION 61 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page iv 5.5 TRACER STUDY RESULTS 64 5.5.1 General Observations 64 5.5.2 Reproducibility 65 5.5.3 Parameter Trends 67 5.6 CONCLUSIONS .68 5.7 RECOMMENDATIONS 69 6 F L U E N T N A V I E R - S T O K E S M O D E L 70 6.1 FLUID FLOW THEORY AND EQUATIONS 70 6.1.1 Fundamental Equations 70 6.1.2 Turbulence Modeling 71 6.1.3 The k-e Model 73 6.1.4 Species Transport 74 6.1.5 Boundary Conditions for Laminar and Turbulent Flows 75 6.1.6 Initial Conditions and Solution Convergence 77 6.2 COMPARISON OF F L U E N T TO LABORATORY DATA 79 6.2.1 Grid Shape 80 6.2.2 Turbulence Modeling and Inlet Configurations 82 6.2.3 Grid Convergence 86 6.3 CONCLUSIONS 87 6.4 RECOMMENDATIONS 87 7 C F D T R A N S P O R T M O D E L 89 7.1 C F D TRANSPORT EQUATIONS AND THEORY 89 7.2 SPATIAL DISCRETIZATION APPROACHES AND ACCURACY 91 7.2.1 Dispersion Flux 91 7.2.2 Advective Flux 92 7.2.3 Semi-Discretized Equations 93 7.3 TIME DISCRETIZATION AND ACCURACY 94 7.4 NUMERICAL SOLUTION TECHNIQUES 95 7.4.1 Boundary Conditions 97 7.4.2 Velocity Boundary Conditions 97 7.4.3 Concentration Boundary Conditions: 98 7.4.4 Import of flow field from FL UENT 98 7.4.5 Initial Conditions for Tracer Study Simulation 99 7.4.6 Overshoot Management for Second Order Upwind 99 7.4.7 Programming Flow 101 7.5 MODEL VERIFICATION 102 7.5.1 Dispersion 103 7.5.2 Advection 103 7.6 R T D MODELING 104 7.6.1 Continuity Errors 105 7.6.2 Matching Laboratory RTD Profiles 109 1.1 CONCLUSIONS 110 7.8 RECOMMENDATIONS I l l 8 C O N C L U S I O N S A N D R E C O M M E N D A T I O N S 112 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page v 8.1 CONCLUSIONS 112 8.2 RECOMMENDATIONS 113 9 REFERENCES 115 APPENDIX A: AERATOR CALIBRATION 118 APPENDIX B: STREAM FUNCTION CALCULATION 119 APPENDIX C: CONTINUITY MACROS 120 APPENDIX D: AUTOCORRELATION CODE 125 APPENDIX E: TRACER STUDY RAW DATA 127 APPENDIX F: RTD PARAMETER MACRO 135 APPENDIX G: SCALAR TRANSPORT CODE 138 List of Figures Figure 2-1: Illustrative RTD Curve 21 Figure 4-1: Model Aerator Schematic 32 Figure 4-2: Basin Schematic 34 Figure 4-3: SonTek A D V Product Photo (A) and Mounted in Laboratory Basin (B) 35 Figure 4-4: Aerator Calibration Curves 41 Figure 4-5: Velocity Measurement Sampling Plane 43 Figure 4-6: Sampling Time Population Comparisons 44 Figure 4-7: U Velocity Profiles - Power 7, Q=3.1 L/sec 47 Figure 4-8: U Velocity Profiles - Power 6, Q= 2.7 L/sec 48 Figure 4-9: U Velocity Profiles - Power 5, Q=2.3 L/sec 49 Figure 4-10: Turbulent Kinetic Energy (k) Profiles 50 Figure 4-11: Sample Autocorrelation Output 52 Figure 4-12: Depth-Averaged Turbulence Values 53 Figure 5-1: Modified basin for Follow-through Conditions 56 Figure 5-2: Optimum Spectrophotometer Wavelength 57 Figure 5-3: Rhodamine-WT Calibration Curve 58 Figure 5-4: Aerator Flow Curve 61 Figure 5-5: Dispersion Number - Variance Curve 64 Figure 5-6: Tracer Study Reproducibility - No Mixing 66 Figure 5-7: Tracer Study Reproducibility - Maximum Mixing 66 Figure 5-9: RTD Parameters Based on Various Aerator Power Settings 67 Figure 5-10: RTD Parameters for Various Flow-through Rates without Aeration 68 Figure 6-1: Sample Convergence History 79 Figure 6-2: Velocity Sampling Points Within Basin 80 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page vi Figure 6-3: Bent 4 Grid Structure 81 Figure 6-4: Fluent Calibration - Grid Shape 82 Figure 6-5: Non-Convergent Laminar Residual History 83 Figure 6-6: Velocity Flow Fields - Comparison of Turbulence and Inlet Configurations 84 Figure 6-7: Flow Calibration - Comparison of Turbulence and Inlet Configurations 85 Figure 6-8: Flow Calibration - Grid Resolution 86 Figure 7-1: Comparison of Advective Flux Methods 100 Figure 7-2: TRANSPORT Code Flowchart 102 Figure 7-3: Dispersion Model Accuracy 103 Figure 7-4: Advection Accuracy 104 Figure 7-5: Sample Solute Transport Time Series - Pulse Tracer Input 105 Figure 7-6: Mass and Outlet Concentration Time Series - Pulse Tracer Input 106 Figure 7-7: Sample Solute Transport Time Series - Continuous Tracer Input 107 Figure 7-8: Mass and Outlet Concentration Time Series - Continuous Tracer Input 108 Figure 7-9: Normalized Continuity Error Map 109 Figure 7-10: Comparison of Model and Laboratory C-Curves 110 List of Tables Table 2-1 - Arceivala Dispersion Equations 15 Table 4-1: Exiting M i l l Volume to Flow Ratios 33 Table 4-2 : Modified Aerator Configuration Parameters 39 Table 4-3: Parameter Settings and Perturbations for Analysis 40 Table 4-4: Aerator Configuration for Profiling 45 Table 5-1: Flow Rate Characteristics 59 Table 5-2: Tracer Study Summary Table 62 Table 5-3: RTD Parameters 62 Table 6-1: Standard k-s Parameter Values 74 Table 6-2: Fluid Properties Used in FLUENT Modeling 75 Table 6-3: Velocity Sampling Points 80 Table 6-4: RMS Error by Mesh Type 82 Table 6-5: RMS Error by Turbulence Modeling and Inlet Configurations 85 Table 6-6: RMS Error by Grid Resolution 86 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page vii List of Symbols and Abbreviations Symbol Description a - approximate factorization left term /3 - approximate factorization center term volume correction factor s - turbulent energy dissipation rate y - approximate factorization right term Q - dimensionless time - basin theoretical hydraulic retention time Q - deviation from plug flow p f parameter p - fluid density K - Von Karman's constant (0.4) H - kinematic viscosity jueff - effective fluid dynamic viscosity /Ut - turbulent kinematic viscosity <Jk - k-s turbulent model constant aE - k-s turbulent model constant OQ - dimensionless time variance cr2 - variance TW - shear stress against wall v - fluid dynamic viscosity vt - turbulent dynamic viscosity yz - stream function A - aerator dimension elevation position calibration matrix for A D V signal conversion A D V - acoustic Doppler velocimeter ASB - aerated stabilization basin A S M - Algebraic Stress Model (turbulence) B - aerator dimension for propeller position C - C function concentration of reactant aerator dimension for draught tube position CIE - k-s turbulent model constant Cie - k-s turbulent model constant Symbol Description CFD - computational fluid dynamics Q concentration at sampling point i Co initial concentration cP concentration at peak of RTD Cr fraction of recovered tracer in tracer study CSTR - continuous flow stirred tank reactor Cu k-s turbulent model constant D diffusion constant d dimensionless dispersion number D water depth D* longitudinal dispersion coefficient based on turbulent diffusion D, longitudinal dispersion coefficient based on shear layer effects A dilution factor DT time increment between successive iterations E E function F F function G k turbulent production term H basin depth h height of water H height of top of aerator draught tube over water surface HRT - theoretical hydraulic retention time / depth of aerator propeller within draught tube / I function k 1 s t order reaction rate constant K constant of dispersion k turbulent kinetic energy I characteristic length L length of reactor L D V - laser Doppler velocimeter M mass (of tracer added to basin for tracer study) MI Morrill index N number of CSTRs in series n perpendicular distance away from wall Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page viii Symbol Description P fluid pressure P power setting of aerator rheostat P probability statistic Pe Peclet Number (1/d) PF plug flow reactor q backmix flow Q volumetric flow rate Qa active flow Qo total reactor flow QR recycle Flow R reaction group (kL/u) R() autocorrelation function RJ normalized update residual Rh hydraulic radius RTD - residence time distribution Sc Schmidt number SNR - signal-to-noise ratio t time t mean or average time tio time for 10% of tracer to pass in tracer study t90 time for 90% of tracer to pass in tracer study taa actual mean residence time tai ideal mean residence time U time of sampling point i tp time to peak concentration of RTD U velocity in x-direction u' turbulent fluctuation in u-velocity u* shear velocity u mean velocity V velocity in y-direction V volume VA active volume VR recycle volume vt beam-aligned A D V velocity vector V cartesian coordinate-aligned A D V velocity vector vd dead space parameter w velocity in z-direction w width of reactor WSP - waste stabilization pond z basin depth z dimensionless distance in reactor hydraulics Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 1 1 Introduction The aerated stabilization basin (ASB) is the most common secondary wastewater treatment system utilized by the forest products industry. ASBs are partially mixed, aerated reactors that degrade wastes biologically prior to discharge into the environment. These systems are characterized by hydraulic retention times of 3 - 10 days, are 1 - 5 m deep and occupy hectares of surface area. Structurally, these systems are very simple, being essentially shallow, lined, flow-through ponds, usually with mechanical surface aerators included to provide the necessary oxygen and mixing required for waste stabilization. It is known that the performance of ASBs and related systems is highly dependent on the flow conditions (Nameche and Vasel, 1998). The modeling of substrate decay typically uses an exponential decay equation based on the time the substrate spends in the reactor. The efficiency of the hydraulics of a system will have a direct impact on the amount of time the substrate will spend in the system, on average, and therefore will have a direct effect on the system performance. A variety of techniques have been developed to examine the hydraulics of reactors and specifically ASBs and waste stabilization ponds, but until recently, these were very simplistic mathematical models that attempted to model the reactor by simplifying the reactor to a representative model that reproduces the bulk fluid mixing response. Recently, more innovative techniques have arisen in the field of computational fluid dynamics (CFD) that examine the more detailed flow fields within reactors, allowing the hydraulic performance and mixing characteristics of the reactor to be extracted from the CFD solution. Some effort has been made to apply the CFD techniques to the problems of waste stabilization ponds and aerated stabilization basins. However, none have considered directly the hydraulic impacts of the aerators themselves. Literature unrelated to CFD analysis has shown the importance of aeration in determining the degree of mixing in an ASB (Nameche and Vasel, 1998). This study seeks to examine the detailed impact of aerators and the incorporation of that impact into a CFD model. It is an exploratory investigation into a single approach to the inclusion of the aerator into a CFD model. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 2 1.1 Scope of Study This study attempted to examine the impacts of aerators on ASBs, using a laboratory scale apparatus. When the apparatus was constructed, several factors that ordinarily may have an impact on ASBs were effectively neglected in this study. The lab apparatus utilized water as a medium rather than a waste stream substrate. No consideration was made for suspended solids that are ubiquitously present in ordinary ASBs. The analysis was not full-scale and the shape of the basin was restricted by the physical dimensions of the laboratory space, and as a consequence the applicability to a full-scale analysis was not addressed. Only a single aerator was included in the study, so the effects of multiple aerators or aerator interaction was not addressed. The aerator was always placed in the center of the tank, so that the effects of aerator placement were not addressed. The study did address several issues related to the inclusion of aeration in reactor tanks. The flow calibration of a model aerator was examined. The impact on the velocity fields in the surrounding waters of an aerator at various power settings was examined. The change in residence time distribution (RTD) was examined in a flow-through tank, with and without aeration, and the matching of flow field velocity measurements and RTD curves from CFD modeling to laboratory data was undertaken. 1.2 Structure of Thesis This thesis is divided into several sections: • Section 2 contains a literature review and outlines the progress made in the fields of ASB and waste stabilization pond (WSP) axial-dispersion analysis, surface mechanical aerator hydraulic studies, CFD techniques in ASBs and reactor design, tracer studies, and acoustic Doppler velocimetry. • Section 3 summarizes the objectives of the study and explains the studies motivation and scope. • Section 4 describes the laboratory setup, the model surface mechanical aerator and equipment used in the analysis of the hydraulic impacts of the aerator. Results of the flow Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 3 calibration of the aerator are presented. The results of the A D V data profile analysis for flow field structure and turbulence are presented. Continuity studies and results are also presented, as well as depth averaged turbulence and mixing profiles. • Section 5 outlines the tracer study procedure and includes analysis and results from the tracer studies performed at various flow rates and at various aerator power settings. • Section 6 describes the use of the commercial FLUENT CFD model package in predicting the flow fields in the flow-through configuration of the basin, including a grid convergence study and flow field verification with A D V laboratory data. A comparison of turbulent models in flow prediction is also presented. • Section 7 describes the development and implementation of a CFD mixing model using FLUENT flow data from Section 6 to predict the RTD profiles of the basin without aeration. • Section 8 summarizes the conclusions and recommendations of this study. 1.3 Acknowledgements The author would like to thank Professors Eric Hall and Gregory Lawrence for their supervision and support, the Natural Sciences and Engineering Research Council of Canada, the NSERC/COFI Industrial Research Chair in Forest Products Waste Management and the Sustainable Forest Management Network of Centres of Excellence for their financial support. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 4 2 Literature Review The science of reactor hydraulics has advanced significantly over the last 40 years. Simple hydraulic models that described ideal flow patterns, such as the plug flow reactor (PFR) and the continuously-stirred tank reactor (CSTR) gradually have given way to more complex models that characterize the degree of imperfection of the flow patterns within reactors. The earlier models were simple, one-dimensional equations with one or two extra parameters, evaluated experimentally, that described the mixing imperfections within the reactors. Difficulties immediately arose in how to evaluate these parameters and whole families of equations were developed to predict these parameters based primarily on reactor geometry, flow rates and mixing energy inputs. These parameters were predominantly evaluated from empirical data, not on mechanistic theory or from the physics of the mixing of the reactors. This fact made the transferability of these equations and the applicability of the one dimensional models, generally questionable. Additionally, the majority of reactors exhibit flow and mixing characteristics that are distinctly three-dimensional, and one dimensional hydraulic equations could not properly describe some hydraulic characteristics of importance, including dead-zones and short-circuiting patterns. The modeling of fluid flow in two or three dimensions can be accomplished by numerically solving the non-linear fluid flow equations over a representational computational domain. Techniques of this nature are often collectively called computational fluid dynamics, or CFD, and these techniques have been in existence since the 1970s. However, it has only been feasible to do small scale modeling at reasonable cost for the last 10 years, with the drop in the price of computational time and memory. This literature review traces the development of knowledge of reactor hydraulics with specific emphasis on the more revolutionary developments and those contributions that apply to aerated stabilization basins in some way. It begins with the general theory of chemical reactor engineering in which the foundations of the one-dimensional hydraulic modeling of reactors are based. The application of reactor theory to waste stabilization ponds, as well as ASBs, is presented. The techniques for performing and analyzing tracer studies are presented and also included are some reviews of recent applications of CFD techniques to fluid flow problems in Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 5 the water resources and pollution control industries. Where applicable, the associated theory is also discussed, with the exception of the CFD theory which is discussed in Sections 6 and 7 where the detailed modeling procedures are outlined. The literature pertaining to the interpretation and collection of velocity data using the acoustic Doppler velocimeter (ADV) is discussed in Section 4 where the laboratory experimental techniques are introduced. 2.1 Chemical Reactor Hydraulics The study of chemical reactor hydraulics begins with the seminal paper presented by Danckwerts (1953). Realizing that flow systems rarely conform to the ideal plug flow or completely mixed assumptions, Danckwerts derived a mathematical and analytical model for the one-dimensional mixing of fluid in a flow-through system as a family of dimensionless equations and diagrams. He was the first to describe the statistical distribution of the amount of time particles spend in a reactor, known as a residence time distribution, or RTD. The four age functions described by Danckwerts were the F-, I-, E- and C-functions. From an empirical point of view, the C-diagram corresponds to a conservative pulse tracer response at a reactor outlet, plotted as a concentration against time. The E-diagram has the identical shape as the C-diagram but is plotted as a dimensionless concentration against time. The area under the E-diagram should be unity and the area under the C diagram should the mass of the pulse tracer. The F-diagram is the integration of the E-diagram over time, so that it represents a cumulative distribution and monotonically approaches unity as time increases. The I-diagram is described as the internal age distribution and is the mirror image of the F-diagram with the I-function beginning at unity at time zero and monotonically approaching zero as time increases. The following equations outline the terminology defined by Danckwerts with a slightly modified and more current nomenclature; (2-1) l - F ( 0 = ^ / ( 0 (2-2) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 6 C(t) = —E(t) (2-3) where, F(t) Eft) C(t) I(t) F-function, E-function, C-function, I-function, V M t Q reactor volume, mass of tracer added, volumetric flow rate, and time. In this paper, Danckwerts also outlined some techniques for incorporating empirically-derived components for the generation of the C- and F- diagrams. For the analysis of water treatment ponds, the most significant was the inclusion of the diffusivity term, D, or its dimensionless dispersion coefficient d=D/LU, where L is a characteristic length of the pond or reactor and U is the flow velocity. Interestingly, the example used to illustrate this case was the flow of a conservative tracer through a packed column, comparing experimental data with predictions based on the F-diagram with a dispersion number. In this case, where the flow was slow, laminar and virtually one-dimensional, one would expect a Fickian-like dispersion to adequately represent the data, and it was shown to do so. The equation for the one-dimensional axial dispersion model is shown below in equation 2-4. dC d2C (2-4) dt dx2 where, x U C D t = axial velocity, = concentration of species, = dispersion number, = distance ordinate, and = time ordinate. Wehner and Wilhelm (1956) presented an analytical solution for the dispersed plug flow model with a decay equation incorporated. Much of the work of the paper examined the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 7 influence of the boundary condition assumptions of the Danckwerts flow model, including the mixing properties immediately before and after the reactor and the effects these assumptions have on the shape of the RTD. The contribution of greatest application to the study of waste treatment was the analytical solution of the dispersed plug flow model with first order decay. This equation has seen wide application in treatment system design (NCASI, 1983). C = go exp Pe-z (1 + a) exp a • Pe(\ - z) - (1 - a) exp a-Pe(z-l) (2-5) go = fa-Pe\ n .2 ( a-Pe^ -(I-a) exp (1 + a) 2 exp (2-6) J Pe = L • u D (2-7) a = ' 4R 1 + — v Pe (2-8) where, L z R k Pe u D C C0 = reactor length [L], = dimensionless distance normalized by L , = reaction group (kL/u), = first order reaction rate constant [T 1 ] , = dimensionless Peclet Number, = linear average velocity [LT 1 ] , = diffusion constant [L 2T"'], = concentration of reactant [ML" 3], and = initial concentration [ML" 3]. Octave Levenspiel contributed significantly to the development of chemical reaction engineering. His paper with Smith on the application of the diffusion model for longitudinal mixing in fluids (Levenspiel and Smith, 1957) outlined the requirements for application of the Danckwerts model to flow-through reactors. In this paper, Levenspiel and Smith outlined the conditions required for appropriate application of the one-dimensional advection-dispersion equations with the mixing parameter d=D/uL (defined above). The longitudinal dispersion must Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 8 be independent of position, the transverse dispersion must be very high with respect to the longitudinal dispersion and there should be no flow variations along the flow path. If any of these three criteria are not met, the application of the dispersed plug flow model is questionable. Particularly, the model does not hold for reactors with low Reynolds numbers or where the geometry of the reactor is such that uniform lateral concentrations are not possible. In addition to these observations, derivations were presented for the equations for RTD profiles along a reactor based on the dispersed plug flow model. Techniques were also presented for calculating dispersion numbers based on experimentally acquired concentrations through the calculation of the RTD variance or using a normalized probability graphs i f the dispersion is small. Elder (1959) presented a technique for the calculation of the dispersion number based on the turbulent shear flow profile typically found in pipe flow or open channel flow. Elder pointed out that the turbulent diffusion coefficient (D) can be related to the turbulence parameters of a fluid by the following equation: D = u"\~R(p)dp (2-9) where, u' - turbulent fluctuation in velocity, and R(p) = autocorrelation function . This equation applies i f the velocity fluctuations are uniform and isotropic and i f the distributions of dispersed particles would follow a Gaussian distribution. Also presented in this paper was the derivation of a longitudinal dispersion coefficient based on the logarithmic law of the wall profile. It was shown that the contribution of the dispersion due to the shear layer contributes much more to the longitudinal dispersion than the isotropic turbulence itself. Elder demonstrated that the longitudinal dispersion due to the shear layer depended on two factors: the height of the water and the shear velocity (which may be derived from the velocity profile). The following equations were presented, based on theoretical derivations. D. = 5.S6hu. (2-10) D = 0.07Au, where, Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 9 D* = longitudinal dispersion coefficient based on turbulent diffusion, h = height of water, and u* = shear velocity. It was shown by Elder that experimental results tended to agree with these dispersion numbers generally, except that the tails of the distributions were longer, indicating that some fluid was becoming trapped in the viscous sub-layer and released slowly. Cholette and Cloutier (1959) introduced the concept of combined models to model flow in reactors. The technique involved the networking of PF and CSTR reactors to manipulate the RTD profiles. Additionally, the authors presented a technique for analyzing a theoretical RTD profile assuming the deviation from the ideal (CSTR) profile could be accounted for by considering the active to theoretical volume ratio, and the active to theoretical flow ratio. Equation 2-11 shows the relationship between the active volume (Va) and the active flow (Qa)-The active flow is typically identical to the theoretical flow, Q. The following equation 2-12 shows the theoretical age distribution function based on the flow and volume ratios, and the third equation, 2-13, is the linearization of that equation that can be used to determine the Qa and Va parameters from experimental data. v.=Q,-i (2-H) Q I(6) = ^ e Q V ° (2-12) to 7(0) = In ^ - - ^ - ^ 0 (2-13) Q QVa Levenspiel (1962) integrated the science of chemical reactor engineering and presented the classical summary of the hydraulic characterization of reactors. With an emphasis on non-ideal flow, Levenspiel proposed a variety of models that predict or replicate RTD functions in the standard form originally presented by Danckwerts. He presented two models that could be used to characterize non-ideal flow systems: the dispersed plug flow and the equal-volume-CSTRs-in-series model, also relating these two systems mathematically. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 10 Levenspiel also outlined many forms of the combined models originally presented by Cholette and Cloutier (1959), using simple flow models in network-like structures to generate more complicated RTD functions. The elements of the network included plug flow reactors with and without dispersion, CSTRs, and dead space elements. The versatility of the combined models provided the ability to characterize the hydraulics of virtually any reactor system. However, very complex systems required a large number of characterization parameters, making the calibration of the system impractical, since the parameters were often not physically meaningful or easily measured. Levenspiel also presented techniques for tracer tests and the analysis of the resulting RTD data. When mapping the data to a dispersed plug flow model the only required parameter was the dispersion number (previously defined). The dispersion number is directly linked to the variance of the RTD distribution. The mean residence time could be calculated by determining the value of the first moment of the distribution. f t-Cdt Jo (.00 I Cdt Jo (2-14) cr = J " ( r - 0 2 Cdt /•CO f Cdt Jo (2-15) where, t = time, C = concentration of tracer at specific time and space within the reactor, i = mean fluid residence time, and a1 = variance. With the variance and mean known, the dispersion number can be calculated by solving the following equation (Levenspiel, 1962): o-' D , — = 2 - 2 t u L yU L j (l-e-uL/D) (2-16) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 11 Levenspiel and Turner (1970) developed techniques for the conduction and interpretation of RTD experiments by addressing the question of how a conservative tracer should be added to a reactor and how the sampling should be conducted. It was shown that the preferred method for accurate results is the addition of a tracer proportional to the velocity equivalent to a turbulent mixing and the sampling of the tracer using a mixed sampling cup reading. This was done by considering two possible approaches for both the addition of the tracer and the measuring of the tracer at the reactor outlet and, using a simplified illustrative example. A more general supplemental examination of this principle was illustrated by Levenspiel et al. (1970) using a laminar flow distortion model. Haddad and Wolf (1967) improved upon the tank-in-series-model by incorporating a back-mixing component allowing for counter-current movement between reactor segments. Therefore, two flow parameters were to be considered, the bulk flow (v) and the backmix flow (q). The q flow was coined as the degree of entrainment. The authors derived the equations for the F- and E-functions using Laplace transforms of the system's differential equations. This backmix-tank-in-series model has found useful application in its ability to model backward dispersion (NCASI, 1983). 2.2 Applied Pond Hydraulics The literature pertaining to studies of the hydraulics of aerated stabilization basins is sparse and fragmentary. Typically, research has focused on related treatment systems, usually waste stabilization ponds (WSPs) and sedimentation tanks. WSPs are hydraulically quite similar to ASBs, the only practical difference being the inclusion of mechanical aerators in ASBs. The hydraulic analysis of waste stabilization ponds borrows heavily from the seminal work of Danckwerts (1953) and Levenspiel (1962) and their techniques for residence time distributions for fluid flow systems. Thirumurthi (1969b, 1974) presented principles for the design of waste stabilization ponds based on the first order decay model. Thirumurthi rejected the notion of modeling WSPs as either PF or CSTR reactors, but rather, stressed the importance of using the axial dispersion model and the Werner-Wilhelm equation for BOD removal. Thirumurthi stated that the CSTR models for BOD removal should never be used in the design of WSPs. In the laboratory scale Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 12 systems analyzed, the axial dispersion model matched tracer data and effluent BOD data very well. It was suggested that the greatest degree of uncertainty lay in the determination of the dispersion parameter, d, and the biological decay coefficient, k. Thirumurthi presented a technique for the evaluation of k for WSPs that accounted for environment, temperature, organic and toxic corrections. The National Council of the Paper Industry for Air and Stream Improvement (NCASI, 1971) produced a technical bulletin which examined the mixing characteristics of aerated stabilization basins. The study examined, under field conditions, the degree of mixing that existed in aerated stabilization basins under variations in depth, geometry, mechanical energy input and other parameters. The study seems to be the first to have examined the hydraulic impacts of individual aerators. The mixing was determined by measuring the velocities in radial profiles at various distances away from a relatively isolated surface mechanical aerator. Velocity measurements were taken until the measured velocity imparted by the aerator could not be distinguished from the bulk flow patterns. It was found that three main factors affected the degree of mixing (or the "zone of influence") of an aerator. These were: (1) the type of aerator (high vs. low rotational speed); (2) the size or mechanical power of the aerator; and (3) the depth of the basin. Mangelson and Watters (1972) examined the effects of geometry and inlet-outlet configurations on the hydraulic performance of laboratory-scale waste stabilization ponds. A large variety of basin geometries was used and a large number of tracer studies were performed. The reproducibility of the tracer studies was found to be limited, especially around the RTD peak, where sensitivity to the hydraulic conditions was high. However, the bulk parameters were found to be more reproducible. It was found that a shallower water depth resulted in flow conditions more closely resembling plug flow. The flow velocity or Reynolds number was found to have no effect on the parameters. The inlet and outlet configurations were found to greatly effect the study results, and the length to width effects had a huge impact on the flow condition. As the length-to-width ratio (L/W) increased, the efficiency also increased, with a decrease in the overall dead space. It should be noted that the bulk parameters examined included the deviation from plug flow and the volume of dead space parameters and that the more standard parameters, namely the Morrill index and dimensionless dispersion number, were not examined. The definitions of the tracer study parameters will follow in Section 2.3. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 13 Murphy and Wilson (1974) presented an evaluation of the various techniques of hydraulic modeling of full-scale ASBs. They evaluated a variety of models as to their ability to predict RTD profiles and BOD degradation. The models included were the CSTR model, the PF model, the equal-volurne-tank-in-series model, the unequal-volume-tank-in-series model, the backmix-tank-in-series model and the axial dispersion model. The model equations were transformed from the time to the frequency domain using Fourier transforms. It was found that the axial dispersion models and the unequal-volume-tank-in-series model produced the most accurate results. First order kinetics were used in the evaluation of waste treatment. An attempt was made to predict the dispersion number, d, based on the basin flow rate, horsepower volume, and residence time, but no significant relationship could be found. Murphy and Wilson suggested a linear relationship between the dispersion number and the theoretical detention time over the square of the characteristic length, as described in equation 2-17 below: D V d = — = K (2-17) uL QL2 where, K = constant (2.88 xl0 3 m 3 /hr) , V = volume of basin, Q = flow rate through basin, L = axial length of basin, and d = dimensionless dispersion number. Ferrera and Harleman (1980,1981) examined the available modeling techniques to address the need of complex combined models or axial dispersion models. Some limitations to the combined and axial models addressed by the authors were the need for input parameters that cannot be measured, but which can only be determined by back-calculating from tracer studies. Further, the tracer study data would only apply to the unique conditions exhibited during the time over which the study was conducted and would not be applicable under different hydraulic conditions. The CSTR model was presented as the best of the traditional options due to its simplicity and performance in predicting waste treatment using first order BOD kinetics. The authors also suggested a new return-flow modeling approach in which the reactor was divided into active (A) and recycle (R) regions. The active zone would conduct flow directly to the outlet and be modeled as a dispersed plug flow reactor with the dispersion parameter (EA). The Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 14 return zone would be assumed to be completely mixed. Three parameters are required to model this type of reactor: Ds, the dilution factor relating the active and recycle flows; ft, relating the active to the recycle volumes; and Ea, the longitudinal dispersion equation mentioned above. Ds = Q± = Q r ~Q° (2-18) Q0 Qo P = ^ ~ (2-19) E A = 0-225 4 ^ (2"2°) where, Ds = dilution factor P = volume factor QA = active flow Qo = total reactor flow QR = recycle flow VA - active volume VR = recycle volume u* = shear velocity / = characteristic length K = Von Karman's constant (0.4) Rh = hydraulic radius Ferrera and Harleman calibrated this model to selected flow data of Mangelson and Waters (1972) and their three parameter model showed curve matches that could not be accomplished with an axial dispersion model alone. The calibration to the Mangelson and Waters data provided a range of Ds and ft parameters that the authors considered broadly representative of typical WSP hydraulic conditions. The dilution factor, Ds, was found to range from 1.4 to 4.0 with an average of about 2.0. The volume correction factor, /? was found to range from 0.24 to 0.76 with an average of about 0.5. Under a wide range of /?and with Ds equal to 1.5 and 2.0 it was seen that a CSTR model approximated the actual first order decay reasonably well and would only overestimate treatment with greater than 0.50 or so. It was demonstrated that plug flow modeling without dispersion would seriously overestimate treatment in all cases examined. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 15 In a discussion paper Arceivala (1981) responded to the presentations of Ferrera and Harleman (1981). The author proposed that the inclusion of baffles in a pond precluded the use of the CSTR model and presented one of the earlier and simpler empirical equations to predict the dispersion number of WSPs for the use of the axial dispersion model. The equations were based on studies of correlation between various geometric and flow parameters and the dispersion number. It was suggested that the dispersion number could be predicted sufficiently from the width of the basin alone, and also depended on the inclusion of baffles in the basin. Table 2-1 outlines these simple dispersion prediction equations: Table 2-1 - Arceivala Dispersion Equations Baffles No Baffles W>30m W < 10m where, D=33W D=16.7W D=11W" D=2WZ D = dispersion number (m2/hr), and W = width of WSP (m). Arceivala also mentioned that, from experience, the dispersion numbers of WSPs typically vary from less than 0.2 to greater than 4.0. Polprasert and Bhattarai (1985) examined various dispersion approximation formulae, including those presented by Ferrera and Harleman (1980) and a modified formula presented by Liu (1977). The authors presented an additional dispersion prediction formula for WSPs which relates the dispersion number to the theoretical HRT, the fluid viscosity and the pond geometry. OAS4[0v(W + 2Z)]0489 (W)xsn ~{LZy d = , , ^ 1 . 4 8 9 (2" 2 1) where, d = dimensionless dispersion number, W = basin width, Z = basin depth, L = basin length, 0 - basin theoretical hydraulic retention time, and v = fluid viscosity. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 16 Polprasert and Bhattarai compared their dispersion prediction equation to tracer studies performed on pilot scale ponds as well as the dispersion equations presented by Ferrera and Harleman (1981) and Liu (1977). It was found that Equation 2-21 produced the best results overall with the lowest standard error. It was also mentioned that the average dispersion values, d, are rarely greater than 1.0 because of low hydraulic loads and that values are expected to fall between 0.155 to 1.710. Sackellares and Barkley (1985) and Sackellares et al. (1987) examined several flow models for use in the analysis of aerated lagoons. The preferred model was to be selected and used in the D Y L A M O model developed for the Weyerhaeuser Company. As the D Y L A M O mechanistic model was using non-linear kinetics, the correct modeling of the flow path was essential. An examination of the following flow models was conducted: • dispersed plug flow model; • equal-volume-CSTRs-in-series model; • unequal-volume-CSTRs-in-series model; and • backflow cell model. The model selection was based on the lowest root-mean-square of the error comparing the estimated RTD functions to those generated from tracer studies. It was found that the equal-volume-CSTRs-in-series model performed the most poorly and that the remaining three were comparable. A decision was made to select the backflow cell model as it was believed that the ASBs assessed did not function hydraulically as CSTRs in series and that recycling was certainly occurring. It was also felt that an understanding of the backflow effects would be critical in determining the solutions to the non-linear biological reactions. Additionally, a CSTR model was required in order to capture the kinetics, as the analytical solution required for the dispersed plug flow would not be tractable. Some difficulties arose with this method in determining the appropriate number of cells and the degree of backmixing, as both are linked to the dispersion number by the following equation. ^ = d-N Qo (2-22) where, Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 17 QR - recycle flow rate, Q0 - bulk forward flow rate, d = dimensionless dispersion number, and N = number of CSTRs in series. This equation allows for an infinite number of recycle rates and N pairs that will produce identical dispersion numbers. However each will produce a distinctly different flow path and consequently different kinetic results. The evaluation of the recycle flow rate is difficult to measure and by the authors' own admission, the elucidation of the recycle flow requires more study. Marecos Do Monte and Mara (1987) reviewed problems with the predictive models developed in the literature. Two ponds were examined both in winter and summer, by means of tracer studies. Dispersion numbers were calculated using both the Levenspiel method (Equation 2-15) and the modified method used for non-uniform sampling frequencies described below in Section 2.3. The dispersion numbers calculated from the tracer data did not agree with either the model presented by Polprasert and Bhattarai (1985) or that of Arceivala (1981). More than anything else, this paper showed the limitations of the empirical techniques for the calculation of dispersion numbers. Preul and Wagner (1987) proposed a backflow model similar to that presented by Ferrara and Harleman (1981), but with a more detailed three-dimensional description and kinetic equations. This model contained an active zone and recycle zone as per the Ferrara and Harleman model but the active zone was modeled as a dispersed plug flow reactor followed by a CSTR reactor. The return zone was modeled as a CSTR including space in the reactor on both lateral sides of the active zone and then either above or below the active zone. That is, the model was conceptually split into two conditions: a "top flow" or a "bottom flow" condition specifying the location of the active zone with respect to the recycle zone. As the model had detailed kinetic equations including benthal layer feedback, the location of the active zone was important. This model suffered from the fact that a very large number of parameters were required even to accurately specify the hydraulics, particularly the size of the three component reactors, the recycle flow rate and the dispersion number in the first plug flow reactor. Tracer data would be required to fit this Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 18 model to an existing system and the top vs. bottom flow conditions would have to be determined with some velocity measurement device. Lyn and Rodi (1990) performed turbulent measurements on a model settling tank with different baffle arrangements using a two dimensional Laser Doppler Velocimeter (LDV) and performed tracer studies as well. Examined turbulence parameters were turbulent kinetic energy estimated using only two components as 0.75(u'2 + v' 2) and the Reynolds shear stress. The L D V was shown to be quite adequate in determining turbulence and flow patterns but was limited in that it could only perform in clear waters and that only two velocity components could be examined. The usefulness of various baffle configurations were compared by examining the impacts on hydraulic efficiency and the tradeoff with an increase in turbulence around the baffles and consequent resuspension of the materials. Some guidelines to modeling settling tanks were presented, stressing the importance of maintaining high Reynolds numbers so that turbulent flow is established in the lab model i f any comparisons are to be made with a full scale model. Agunwamba (1992) presented some cautionary notes regarding laboratory and pilot scale waste stabilization basin analysis and scalability. The ability to match the hydrodynamics, density differences, and shear stress were seen to make the overall interpretation of laboratory results difficult. Three general areas of concern were cited: pond geometry, flow conditions and environmental conditions. It has been found in the literature that both the L /W and W/H ratios play significant roles in predicting the dispersion number in a pond. The flow conditions are also significant in that the Reynolds number in smaller scale ponds is typically lower than those in full-scale systems and that the environmental conditions cannot be easily reproduced. Density gradients often exist in full-scale systems, leading to a degree of short-circuiting, and density stratification is also suspected to affect transverse mixing. The exchange of fluid over a stratification layer was found to vary with the square root of pond depth so that deeper ponds will exhibit a more stable stratification, i f present. Agunwamba et al.(1992) described a more mechanistic approach to the prediction of a dispersion number for one-dimensional longitudinal mixing equation. The model presented in this paper described the semi-empirical estimation of the dispersion number based on the longitudinal dispersion developed by a turbulent shear layer. The equations relied heavily on the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 19 theory developed by Fisher et al. (1979) in their work on diffusion in streams by applying these mixing phenomena to waste stabilization ponds. Some of the difficulties in this model lie in the use of the transverse mixing coefficient rather than the vertical mixing coefficient that influences the mixing in rivers. The transverse and lateral mixing coefficients were considered to be proportional to each other and the ratio of the two values was considered representable by a constant. The overall dispersion number is parameterized to three empirical coefficients. These coefficients relate to the geometry and hydraulic conditions of the WSP in question so that transferability is questionable. The basic structure of the equation is a slight modification of that from Polprasert and Bhattarai (1985). Pearson et al. (1995) examined the effects of different L /W ratios (from 1 to 6) on WSP performance and found there was little effect on performance and effluent quality. However, the hydraulic performance was not addressed directly. Rather, the performance of the ponds was related to the overall improvement in BOD, suspended solids and fecal coliform reduction. The authors suspected that density stratification was probably more important than pond geometry and that as long as the L /W ratios were of a reasonable nature (i.e. constructed to avoid severe short-circuiting) the performance should be adequate. The authors also noticed that the first order decay constants generally decreased over the length of a pond, which they attributed to recalcitrant organic chemicals. This effect would tend to mask the hydraulic improvements made by an increased L /W ratio. Nameche and Vasel (1998) examined the existing dispersion number prediction formulae for ASBs and WSPs and disagreed with Murphy and Wilson (1974), who had stated that aerator power was unimportant in determining the dispersion coefficient. Three new empirical equations were offered based on tracer studies performed at 30 full scale facilities relating the Peclet number (Pe) to the dimensional ratios and the input aerator power. Pe 0.52 + 0.14 — W (2-23) Pe 0.31 — W + 0.055-Z (2-24) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 20 Pe = 0.35 — + 0 .012-W Z (2-25) where, Pe L Z w Pi Peclet number (l/d), basin length (m), mean basin depth (m), basin width (m), and power input per unit volume (W/m ). The proposed equations statistically showed an improved performance over most other empirical models, including those presented by Polprasert and Bhattarai (1985), Liu (1977), Fisher et al. (1979), Arceivala (1981), and Murphy and Wilson (1974). Equation 2-23 applies to aerated lagoons, equation 2-24 applies to waste stabilization ponds and equation 2-25 applies to all basins. These equations suffer from the fact that they are purely statistical and that there is no physical basis for the relationships. The greatest contribution was to illustrate that the impact of aeration is significant in ASBs, an assertion that had not been refuted since Murphy and Wilson (1974). 2.3 Tracer Studies on Treatment Ponds The methods of conducting and interpreting tracer studies in reactors have evolved over the years through experience in application. A whole family of measured parameters have been extracted from the tracer study distribution curves in order to aid in the categorization and interpretation of them. Figure 2-1 below shows a sample RTD curve with some of the typical important statistical parameters indicated. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 21 U (0 o c o (0 u c o u Observed Recovery of Tracer Time tmin tlO tp t50 t t90 tmax Figure 2-1: Illustrative R T D Curve The key time parameters illustrated in Figure 2-1 are: tmin = time at which tracer is first detected, (max = time at which tracer is last detected, t„ = time at which n% of the added tracer has passed through the outlet, tp = time of maximum or peak tracer concentration, and t = time to reach centroid of the curve or actual mean detention time. Some other parameters have also been calculated from RTD data. One of the primary parameters that is extracted from tracer study data collected for ASBs and WSPs is the dimensionless dispersion number, d. This parameter is directly related to the variance of the distribution. However, the determination of the parameter depends on the inlet/outlet configurations assumed. The equation used for a system with a no-diffusion inlet or outlet is defined as follows (Levenspiel, 1962) al =2d-2d'[l-e-Vd] (2-26) where, Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 22 d - dimensionless dispersion number, and a2e = dimensionless distribution variance. However, i f the sampling data collected from the tracer study are not equally spaced in time, the following equation proposed by Marecos do Monte and Mara (1987) can be used to calculate the variance: .2 i=2 2 2 2 Qt 72 (2-27) where, Q t Q u = the mass of tracer recovered during the tracer study, = tm, the time to reach the centroid of the RTD curve, = the concentration at sampling point I, and = the time at sampling point i . Once the variance is known, d can be calculated from Equation 2-26 using trial and error methods. Thirumurthi (1969a) applied the theories of Danckwerts and Levenspiel to the analysis of settling tanks and performed a series of tracer studies on laboratory scale settling tanks. The series of 10 tests was used to define the hydraulic performance of these settling tanks and the relevance of various parameters to represent the hydraulics. It was shown that the dispersion number, d, was generally reproducible with a maximum error of 4.7%. However, other parameters were seen to be much less robust with ti and tp varying by as much as 84%. It was generally observed that d increased with an increase in the flow rate within the settling tanks. Tracer recovery ranged between 70-90% and this was seen as adequate for all calculations; no adjustment or correction of data was seen as necessary. Mangelson and Waters (1972) performed a large number of tracer studies on laboratory scale WSPs and used the hydraulic data to predict performance. The authors presented two parameters for the evaluation of the pond: the deviation from plug flow (6P/) and a dead space parameter Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 23 (Vd). The deviation from plug flow parameter was defined as the waste that leaves the detention pond before one detention time. As a reactor improves in efficiency it will approach plug flow conditions and 0pf wil l approach zero. The dead space parameter was defined by the amount of material that remained in the reactor after 2 theoretical hydraulic retention times. Additionally, the authors suggested that the parameters describing an RTD can be calculated using only the first 2 hydraulic retention times for sampling, because values beyond this point are typically at or near the limit of detectablility. The mathematical definitions of 9 f and Vd and #c,the dimensionless mean hydraulic retention time, are illustrated in equations 2-28, 2-29 and 2-30: \\\-9)(C/C0)d9 °Pr = i L 7 r — ( 2 - 2 8 ) l(C/C0)d9 Vd = \-9cF(9 = 2) (2-29) \29-(C/C0)d9 9C = i 5 _ (2-30) l(C/C0)d9 where, 9 = dimensionless time variable, normalized to the theoretical HRT, C = concentration recorded at the outlet, C0 = instantaneous concentration in the basin when tracer mass added assuming complete mixing, F(9) = tracer F-curve, evaluated at 9, 9pf = deviation from plug flow parameter, Vd = volume of dead space, and 9C = dimensionless mean time. Torres et al. (1977) performed two tracer studies on a pilot-scale, deep, wastewater stabilization pond for which thermal stratification was a concern. One test was performed in the summer, the other, in winter, to try to capture both stratified and unstratified conditions. Sulforhodamine B was used as the conservative tracer and it was applied such that a uniform Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 24 submerged in the inlet stream overnight so that no artificial density differences would be introduced due to temperature variations. Once the tracer was added, the pond was sampled at various locations and depths and at the outlet. The Cholette and Cloutier (1959) combined model was calibrated to the tracer data in order to determine the pond efficiency. The winter months exhibited a weak stratification with an average of a 70% active volume within the WSP. The summer test showed a strong stratification with the active volume dropping to only 22%. This study illustrated the serious impact stratification can have on hydraulic performance and the applicability of the Cholette and Cloutier modeling technique. The National Council of the Paper Industry for Air and Stream Improvement (NCASI, 1983) presented a technical bulletin outlining recommended procedures for conducting and interpreting tracer studies in ASBs. The bulletin examined techniques for both fluorescent dyes and lithium salt tracer studies, but only the former will be reviewed here. It was recommended that the peak measured tracer concentration should be approximately lOOx greater than the background detected concentration and that a conservative estimate is to assume that the tracer added is initially dispersed completely throughout the basin (i.e. CSTR). It is also recommended that the tracer be diluted into a medium that is similar to that of the basin to avoid unusual mixing or stratification. The point of introduction should result in complete mixing of the material into the influent flow or there should be a large degree of mixing immediately after the tracer enters the system. If the temperature difference between the tracer solution and the wastewater is significant, mixing becomes critical. Some additional mixing may be required at the inlet i f no natural mixing environment exists. The duration of the introduction of the tracer should be short relative to the basin theoretical HRT, so that the tracer represents a plug introduction as much as is possible. If the time for introduction is less than 1% of the basin HRT, then the modeling should be adequate. The effluent point should also be well mixed and sampling should be conducted for at least 2 HRTs with a minimum of 25 samples. Generally, sampling frequency should be higher in the first 25% of the testing period, with an increase in the interval between samples as the test progresses. Duplicate samples are recommended when possible. NCASI cited two of the more important analysis parameters as the efficiency (E) and the Morrill Index (MI), the equations for which are shown below. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 25 £ = ^ • 1 0 0 % (2-31) MI = ^ (2-32) 1^0 where, tna = actual mean residence time of the distribution tai = ideal mean residence time (plug flow) MI = Morrill index tgo = time for 90% of the tracer to leave the reactor tio = time for 10% of the tracer to leave the reactor Dorego and Leduc (1996) performed a number of tracer studies on a series of aerated basins with varying degrees of aeration, to examine the hydraulic performance of the systems. The lagoons were characterized by circulatory flow and short-circuiting, with long tracer tails indicating dead zones. These irregular flow patterns invalidate the use of the dispersed plug flow model. Rhodamine WT was used as a tracer and dosed such that complete mixing would yield a concentration of 2.3 rUg/L and the tracer sampling frequency was high at the beginning and reduced gradually toward the end of the tests. The predictive formulae proposed by Polprasert and Bhattarai (1985) in Equation 2-21, and Arceivala (1981) in Table 2-1 were evaluated. The dispersion numbers were calculated from the tracer data using both the Levenspiel (1962) method by Equation 2-15 and the Marcos do Monte and Mara (1987) method by Equation 2-27. It was found that, due to the staggering of the sampling, the Levenspiel equation provided unreliable results and the Marcos do Monte and Mara equation was used to calculate the observed dispersion numbers. The Arceivala equation tended to seriously overestimate the dispersion numbers. The Polprasert Bhattarai equation provided fairly good prediction of dispersion numbers on those basins with a regular geometry. 2.4 Computational Fluid Dynamics in Reactor Hydraulics and Environmental Engineering Recently, new computational techniques have emerged that promise to supercede the mathematically simplistic approaches historically employed in the analysis of the hydraulics of ASBs and other reactor systems. Computational fluid dynamics (CFD) is a group of techniques Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 26 for solving fluid flow equations over discretized computational space and time. The solutions gleaned from these techniques allow for approximate solutions of the Navier-Stokes equations, among others, for which analytical solutions are intractable. The solutions from these CFD techniques can provide spatial and temporal detail of the flows in the computational domain proportional to the degree of discretization resolution within the computational space. The detail provided by these techniques can be extremely valuable in the analysis of reactor hydraulics and can provide significant advantages over the bulk mixing analysis techniques discussed above. In the late 1990's researchers began to turn to computational fluid dynamics to predict the fluid flow characteristics in WSPs. One of the first papers to demonstrate the usefulness of this technique was that by Wood et al. (1995). The work by Wood et al. (1995) and Wood (1997) presented an application using a commercial CFD package, FID AT, and demonstrated its potential as a tool for flow prediction in WSPs. Wood examined primarily the influences of a variety of inlet and outlet configurations and compared both two-dimensional (2D) and three-dimensional (3D) models to empirical data. Some preliminary modeling of aerators was undertaken, but in cursory detail. Additionally, Wood examined the results of tracer studies performed on ponds and tested the ability of a CFD model to predict the RTD profiles, with reasonable success. The work of Scott et al. (1998a, 1998b) described the use of a 2D CFD model, designed using a commercial package, to predict improved hydraulic performance of ASBs through the strategic placement of baffles. The study showed distinct changes in the predicted RTD profiles with different baffle setups and illustrated the power and usefulness of this technique. In addition the work by Scott et al. did not consider the impacts of the surface mechanical aerators that were normally in operation in the ASBs. Armenante et al. (1997) outlined the development of a 3D CFD model using a commercial package (FLUENT) to predict flow properties in a closed, mixing vessel. The vessel contained only a single, liquid phase and was mixed by an axial propeller blade. The model results were verified with flow measurements taken from a laboratory system using a laser Doppler velocimeter (LDV). Data from the L D V were also used to provide velocity and turbulence boundary conditions at the mixing blade. The modeling incorporated both the k-e and the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 27 algebraic stress (ASM) turbulence models and examined both the accuracy of the flow fields within the model as well as the turbulent kinetic energy; The performance of the model was quite good, with the strongest predictions associated with the modeling with the A S M turbulence model over the k-e model. Celik et al. (1985) performed an analysis that was very similar to the present study, but concentrating on settling tanks. The authors examined the hydrodynamics of laboratory-scale rectangular settling tanks using 2D CFD modeling in an effort to improve hydraulic efficiency. The modeling employed the standard k-e model with the standard coefficients and used the law of the wall on the wall boundary conditions. The authors also included boundary conditions for inlet turbulent parameters and used a 3 part power law for the inlet velocity distribution. It was found that the inlet velocity profile had a much greater impact on the flow than the inlet turbulent parameters. The authors also performed tracer studies on the settling tanks to examine the effects of various configurations on bulk hydraulic performance and found very reasonable matches with the RTD profiles generated with the CFD model. The modeling results were very good, considering the coarseness of the grid (a 25x25 uniform mesh). For turbulent transport of the tracer the authors used a Schmidt number of unity. Ye and McCorquodale (1997) performed an interesting CFD model analysis on curved channels using the shallow water equations. These equations are modified versions of the Navier-Stokes equations but which include effects of bottom shear stress as a momentum sink. The authors also used the k-s model and no-slip wall boundaries. A non-uniform unstructured mesh was used to overcome the problems of modeling a curved computational space with an orthogonal grid. The model showed successful application in a variety of conditions including open-channel flows, supercritical flows and hydraulic jump and meandering channels. Ta and Brignal (1998) modeled storage reservoirs in much the same way that WSPs and ASBs are modeled using CFD techniques. The study examined the hydrodynamics of a 200,000 m 3 reservoir with a relatively steady flow rate of 10 ML/day. The residence time distributions were required for proposed inlet configurations to increase the mean hydraulic retention time and reduce short-circuiting. A three-dimensional model was employed using a commercial CFD package. The flow was modeled using the k-e turbulence model with various inlet and outlet Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 28 configurations. RTD curves were generated based on the flow fields calculated. No validation of the model results was performed, but the model did provide assistance in interpreting design options. He et al. (1999) examined the effects of a variable Schmidt number in the modeling of turbulent mixing in a jet-in-crossflow simulation CFD model. They employed a three-dimensional structured grid CFD approach with k-e modeling with species transport using a range of values for the Schmidt number. Although most modeling assumes a turbulent Schmidt number of 0.8 for this type of flow system, He et al. recommended a lower Schmidt number, with the best results with a number of 0.2. 2.5 Conclusions It can be seen from the 30 years of research into reactor hydraulics that progress has certainly been made along two distinct paths: bulk flow modeling using primarily axial dispersion theory; and a more recent evolution in the detailed flow modeling of fluids within reactors using CFD theory. The axial dispersion model benefits from a mathematical simplicity and an ease of implementation. The only parameters generally required for implementation are the reactor volume and some measure of the dispersion number, both of which can be estimated from empirical equations applicable for various physical conditions (Nameche and Vasel, 1998; Polprasert and Bhattarai, 1995; and Arceivala, 1981). If the underlying assumptions of dispersed plug flow can be met, the dispersion number can often be determined from physical or semi-empirical means. Despite its simplicity, the single reactor dispersed plug flow model has shown limitations in the shapes of the RTD profiles it can generate. Specifically, it is unable to model short-circuiting or recycling of flow (conditions not observed in true dispersed plug flow). To accommodate for this inherent inadequacy, researchers have developed combined models, subdividing the reactor into several distinct volumes and applying a dispersed plug flow model to each volume, each with different volumes and mixing parameters. These combined models have an infinite degree of versatility and could theoretically reproduce any RTD profile obtained with a tracer study. However, with the increased versatility comes an increase in the number of parameters to calibrate, raising questions of how to determine reactor volumes and dispersion numbers for each Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 29 component of the reactor model. The combined modeling approach marks a step away from the physics of mixing. The advent of computational fluid dynamics allowed for more detailed reactor modeling. CFD techniques are significantly more complex than the bulk flow modeling techniques described above, and require sophisticated fluid modeling code and skill in computer modeling and fluid dynamic theory. The tradeoff for the extra modeling effort is balanced by the detailed flow information generated for a reactor. There are a number of parameters (or boundary conditions) that have to approximated in the fluid modeling, but unlike the parameters in the combined models, these are purely physically-based and can be approximated with reasonable accuracy. In the future, for the hydraulic and chemical modeling of reactor performance one will likely observe a continuation of the trend away from the bulk flow modeling and toward the detailed flow modeling of CFD. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 30 3 Research Objectives Research has shown that CFD techniques can be successfully applied to WSP systems to model and predict hydrodynamic performance with reasonable success (Wood, 1997). However, unlike ASBs, the hydrodynamics of WSPs are dictated largely by the fundamental fluid flow equations and are not influenced by fluid mixing imparted by aeration devices. The influence of aeration in waste treatment ponds has not been directly addressed in regards to CFD modeling to date, and it was this study's main objective to explore the opportunities to introduce the effects of surface mechanical aerator mixing into a two-dimensional CFD model. The itemized list of objectives was as follows: 1. To determine, using a laboratory scale apparatus, the velocity and turbulence profiles induced around a small-scale surface mechanical aerator, whereby a determination of mixing imparted by the aerator could be calculated as a function of distance from the aerator and aerator flow rate. Detailed knowledge of flow patterns is critical for the incorporation of an aerator model into a CFD routine. 2. To conduct a series of tracer studies to provide detailed bulk-mixing relationships between degrees of mixing in the basin and the shape of the residence time distribution curves; and 3. To produce a two-dimensional CFD model that incorporates the mixing influence of a surface mechanical aerator and accurately reproduces the RTD curves of the of a basin with and without mixing. The scope and limitations of study for each objective are addressed with more detail within the chapters that follow. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 31 4 Aerator Studies In order to examine the effects of an aerator on the surrounding water body, a laboratory-scale apparatus was developed. A tank was built, as was a small scale surface mechanical aerator. Tests were performed without inflow into the tank to examine the hydraulic impacts without the influence of bulk fluid flow. This chapter reviews the following: • the construction of the basin and the aerator apparatus; • the Acoustic Doppler Velocimeter (ADV) and its use in flow measurement; • the flow sampling technique; and • turbulence and velocity data and analysis. 4.1 Aerator and Basin Construction A laboratory-scale tank and an aerator were designed and built to test hydraulic impacts of an aerator on its surrounding water. The aerator was designed to be similar to many of the types of surface mechanical aerators used by the wastewater treatment industry. A scale drawing of the aerator built for this study is provided in Figure 4-1. The flow through the aerator was driven by a 3.8 cm diameter impeller sleeved by a 20 cm diameter draught tube. The draught tube was removable and although several lengths of tube were available, only the 20 cm draught tube was employed in this study. The water drawn up through the tube was sprayed laterally by a splash plate. A 90-Watt DC motor, the power to which was controlled by a 10-point rheostat, powered the propeller. The entire apparatus was supported by stilts and rested on the floor of the study basin. The aerator was designed so that many of the important dimensions could be adjusted. The position of the motor could be moved relative to the splash plate, allowing for independent positioning of the propeller (shown as dimension " B " in Figure 4-1). The position of the draught tube could be moved up or down relative to the splash plate (dimension "C") and the entire apparatus could be moved up and down along the main support stilts allowing for accurate positioning within the water column (dimension "A") . The width and height dimensions of the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 32 aerator are shown in Figure 4-1 to illustrate the size of the aerator. It had a base diameter of 38.1 cm (15 inches) and a vertical support rod height of 83.8 cm (33 inches). f A ob-15" [0.38m] .MOTOR oo SPLASH .PLATE PROPELLER DRAUGHT TUBE Figure 4-1: Model Aerator Schematic The basin in which the testing of the aerator took place was 3.7 m wide, 5.5 m long and 0.61 m deep. The basin was configured for still water aerator studies (i.e. no ambient flow field). The basin was filled with water from the Vancouver water distribution network, and was seeded with silica particles of approximately 8 pm in diameter to increase the acoustic response of the velocity measurement device, described below. Approximately 2000 mL of the silica particles were used to seed the basin when full, but the amounts added were not precise and silica was Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 33 added only as required. Figure 4-2 shows a construction drawing of the basin, including the support struts along the sides. The rationale in determining an appropriate size for the test basin and aeration device involved some comparison of aerator flow rates to reactor volume ratios and attempting to size the aerator accordingly for this study. The size of the basin was essentially fixed by the available space in the laboratory at approximately 7.0 m 3 (at a water depth of 0.42 cm). The aerator was designed to produce a comparable flow ratio observed in the field. Communications with an aerator manufacturer (Aqua-Aerobic Systems Inc.) provided some approximate flow rates for 75 HP surface mechanical aerators. These aerators are employed at an ASB in Grande Prarie, Alberta. A n additional study provided in the literature (Aqua-Aerobic, 1998) sited a reactor with identical mechanical aerators in a mill aerated lagoon system in Boyle, Alberta. Table 4-1 below illustrates the approximate Volume to Flow ratios. Table 4-1: Exiting Mill Volume to Flow Ratios. Mill Volume Number of Volume / Flow Target Model Aerator (m3) Aerators Ratio (hrs) Flow Rate (L/sec) Grande Prarie, A B 836 000 30 6.9 0.30 Boyle, A B 100 000 21 1.2 1.9 Assuming that these were fairly typical values of aerator flow density, the aerator designed for the present study should have an approximate flow range of 0.3 L/sec to 1.9 L/sec. It was difficult to determine the flow of the aerator prior to actual construction but with the knowledge of the pitch, and diameter of the impeller blade and approximate rotations per minute of the motor known the flow rate of the aerator could be approximated. This was done by assuming the effect of gravity was small and that the water would be conducted through the draught tube very much like a solid, with the velocity equal to the rotational speed multiplied by the blade pitch. Allowances were made for efficiencies of 70% to 90% to provide a reasonable range of flow rates. The actual flow range after construction was higher than the proposed range with a lower limit of about 0.5 L/sec and an upper limit of about 4.0 L/sec. However, this higher flow rate was not seen as a real difficulty, as the scaling shown in Table 4-1 did not consider the increased contribution of viscosity at the lower scales, therefore a higher degree of mixing would be Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 34 desirable when linearly relating the flow scales. For mixing, one would like to find similarity with Reynolds numbers between the two systems in the surrounding waters, but with the flow fields around the full-scale aerators uninvestigated and the issues of splashing complicated the similarity analysis. It was decided that the aerator, once built with a reasonable flow and mixing rate, would exhibit a degree of mixing that could be measured and hopefully modeled. Success in modeling of the aerator in this tank would not necessarily scale directly to a full scale apparatus, but the techniques used in modeling this system could be applied to a full scale modeling endeavor. ] :0.30n If t 12f t [ . m] || i - T in i 1 [ 3,66m] j - - j - 2 f t [0.61n] t E f t [0.6 1 [0.30m] l i 1 • ion t [5.4 • l o r —— If-t [0.30m] 7 D J -1 7 / 1 6 f t [0.43m] - 2 f t [0.61m] 1 5 / 8 f t [0.50m]-Figure 4-2: Basin Schematic Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 35 4.2 Velocity Measurement Instrument Velocity measurements were made using an acoustic Doppler velocimeter (ADV) developed by SonTek, Inc. of San Diego, California. The A D V measures fluid velocities in three dimensions based on the Doppler phase-shift principle and is capable of resolving velocities as low as 0.1 cm/sec. The sampling rate of this instrument is 25 Hz, allowing for a degree of turbulence measurement. The A D V was mounted on a vertical shaft attached to a motorized bridge that moved along the length of the basin. This arrangement allowed for easy sampling anywhere within the basin. A picture of the mounted A D V is shown in Figure 4-3, below. Recently much work has been done by researchers examining the strengths and limitations of the SonTek A D V . Some of the relevant findings pertaining both to theory and to application are discussed below. Figure 4-3: SonTek A D V Product Photo (A) and Mounted in Laboratory Basin (B) Lohrmann et al. (1994) produced one of the first papers introducing the SonTek acoustic Doppler velocimeter or A D V as an alternative method for velocity measurement. The authors Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 36 described the science and structure of the A D V and outlined some of the strengths and limitations of the instrument for resolving mean velocity and turbulence parameters. The instrument consists of three parts - a 40 cm long probe, a conditioning module that attaches to the probe, and a processor that fits into the motherboard of a PC and is attached to the conditioning module via a 10 m cable. The probe itself consists of a long stem and three acoustic receiver elements positioned in 120° increments on a circle around the transmitter. Depending on the angle of the receivers, the sampling volume can be located at 5 or 10 cm from the probe, making the probe non-intrusive. The A D V operates on the acoustic Doppler principle, sending acoustic pulses from the transmitter that are reflected off particles in the water. The Doppler shift in the frequency of the acoustic pulse is translated into a velocity component in line with the reflected acoustic beam. It is important that the signal strength be significantly greater than the background Doppler noise in order to resolve the velocity. A variety of signal to noise ratio (SNR) limits have been recommended, and the SonTek manual suggests that a SNR should be greater than 15 dB i f any sort of turbulence measurement is to be resolved. The SNR can be increased with the inclusion of small particles in the water to increase the noise reflection, as mentioned above. The A D V does not need calibration for mean velocity magnitude or direction, as the calibration is only dependant on the shape of the probes and is performed by SonTek prior to distribution. The calibration consists only of a transformation matrix that translates the velocity signals collected by the reflected sound waves to Cartesian coordinates. The instrument can resolve velocities from 1 to 250 cm/s and laboratory studies by Lohrmann et al. (1994) showed excellent velocity resolution, even below 1 cm/s. The A D V is accurate to within 0.25% of the set velocity range i f well calibrated (i.e. probe not damaged), or ± 0.25 cm/s, whichever is greater. Precise physical placement of the probe is preferable but mild tilts of the probe (<5°) wil l not adversely effect the velocity measurement. Comparison with laser Doppler velocimeter (LDV) data in a laboratory wake channel showed excellent correlation (Kraus et al., 1994) Determining values of turbulence parameters is notably less straightforward than determining the mean velocity values. Fluctuations in a velocity signal obtained by the A D V can be due to one of two factors - a) actual velocity fluctuations; or b) Doppler noise (Lohrmann et al., 1994). If the Doppler noise and instrument noise can be eliminated from the beam signal, then the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 37 remaining turbulent values can reliably be considered. However, accomplishing this is complicated by the fact that the accompanying software for A D V signal interpretation provided by SonTek translates the beam velocities into Cartesian coordinates, effectively distributing the Doppler error into the translated signal. As a consequence, when turbulence measurements are taken they wil l consistently be biased high due to the additional Doppler noise. Due to the orientation of the receivers, the contribution of the Doppler noise was especially prevalent in the plane parallel to the circle of signal receivers. A n alternative to using the unmodified turbulence data is to use a low-pass filter and effectively remove the lower-frequency noise, and in this case the turbulent kinetic energy wil l be biased low (Lohrmann et al., 1995). Equations 4-1 to 4-4 below describe how the acoustic signal is processed and translated into velocity measurements in Cartesian coordinates. The original velocity vectors in the beam coordinate (Vt) are transformed by the calibration matrix (A), which is a 3x3 matrix, the values of its elements being determined in laboratory tests by SonTek prior to distribution (i.e. one-time calibration). The resulting vector (V) is the velocity signal in Cartesian coordinates. Equation 4-2 shows how the beam velocity signal can be broken up into a time averaged velocity ( V ) and two fluctuation components. The first component of the velocity is the actual velocity fluctuation (or turbulence) and the second is the Doppler noise. One can see in Equation 4-4 that the velocity fluctuations in the final Cartesian coordinate system include the sum of the actual fluctuations and the fluctuations due to the Doppler noise - both transformed by the calibration matrix (Lohrmann et al., 1994). V = A • V, (4-1) (4-2) V\+Vdi=0 (4-3) 1^ = 1 W)'2+IW,)2 (4-4) Voulgaris and Trowbridge (1998) sought to quantify the noise contribution to turbulence measurement by comparing the L D V and A D V on the same sampling volume in a laboratory Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 38 experiment. The nature of the error generated in the sampling by each instrument allowed for the determination of the actual or "ground truth" values of the turbulence parameters, as opposed to the comparison of two instruments without knowledge of the actual values of the parameters being measured. The authors presented a detailed and complicated analysis of the responses from each instrument, but in summary, it was shown that the A D V was very accurate in measuring velocities and the vertical turbulent fluctuations. In addition, the measurement of the Reynolds stresses were very accurate, except where the eddies were very small near boundaries. The A D V will overestimate the horizontal velocity fluctuations due to the Doppler noise addition mentioned above. The authors recommended that the removal of noise should involve the detection of the spectral noise floor sample and the integration of the spectra up to that point, shifting the turbulent bias from high to low. This present study did not filter the data, but used the data biased high. 4.3 Experimental Method The experimental method in this stage of the study involved two distinct phases: the calibration of the aerator and the velocity profile sampling of the aerator in a still basin. 4.3.1 Aerator Calibration The aerator was calibrated to determine how the physical configuration would influence the flow through it. Parameters that were varied were the power setting, the height of the top of the draught tube above the water surface (pumping head), the position of the propeller within the draught tube and the proximity of the splash plate to the top of the draught tube. There are five configuration parameters that can be adjusted on the aerator that were thought to have some impact on the flow field generated by the aerator itself. Three of these parameters describe the physical configuration of the aerator and were denoted as A, B and C and were illustrated previously in Figure 4-1. The power setting of the rheostat when the aerator was in operation regulated the rotational speed of the motor and hence the flow rate of water through the aerator. The incremental power setting as shown on the rheostat was called parameter P, for this study. The depth of the water was also believed to be an important factor that influenced the flow rate. As the water depth increases, the hydrostatic head the aerator must overcome is Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 39 reduced and the flow rate will increase through the aerator. The water depth was designated as D in this study. The aerator was calibrated by placing it in a smaller tank of water. A waterproof sheet was placed over the water and attached around the draught tube, producing a septum that separated the water that passed through the aerator from the remaining water. The sheet was loosely held in place, allowing it to sink with water that was pumped, thus allowing an accumulation of water above the septum without changing the hydrostatic head of the water in the tank. The aerator was allowed to run until the volume pumped was approximately 60 L. Runs were timed and then the water retained above the septum was siphoned off and the volume measured using graduated cylinders. The parameters A, B, C and D exhibit some physical redundancy, as they can be independently manipulated to produce similar physical situations with presumably identical flow patterns. In the interest of identifying the parameters that had a physical meaning, two substitute parameters were introduced to represent the possible physical configuration (I and H). The definitions of I and H are shown in Table 4-2 below. Table 4-2 : Modified Aerator Configuration Parameters Modified Parameter Formula Physical Description I 1 = 34.23- B - C The depth of the propeller within the draught [cm] tube H H = A - C - D [cm] The height of the top of the draught tube above the water surface. A set of typical or mean values for the parameters were established that seemed to represent the approximate physical means of the possible configurations. The parameters were varied about these typical values to determine the impacts of the variation. Table 4-3 below shows the typical values as well as the minimum and maximum values that were employed in the aerator calibration. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 40 Table 4-3: Parameter Settings and Perturbations for Analysis Aerator Configuration Typical Min Max Parameter Value Value Value A 55 cm 54 cm 62.5 cm B 20 cm 8 cm 23 cm C 5 cm 5 cm 8.5 cm D 48 cm 45 cm 48.2 cm P 3 2 8 H 2 cm 0.5 cm 6 cm I 11 cm 6.2 cm 17.7 cm A statistical multiple linear regression analysis was performed on the data set examining the influence of the last three parameters shown in Table 4-3 (P, H , and I) on the flow rate measured through the aerator. Appendix A contains the results from the regression analysis. The t-test values were compared to the two-sided t-test (Freund, 1962) for 59 degrees of freedom requiring an absolute t value of 2.045 for slope significance and only the P and H parameters were found to be significant. From the results of the statistical analysis, a series of flow curves was generated based on the variations that could be found when adhering to the typical parameter values, but adjusting the P and H values. The resulting calibration curves are presented in Figure 4-4, in which the flow rate through the aerator is expressed as a function of hydrostatic head (H) between the water elevation and the exit point of the draught tube for various rheostat power settings. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 41 4000 5 o 3000 — \ o CD £ 2000 to or 1000 o.o 1.0 2.0 3.0 Height (cm) 4.0 5.0 6.0 Figure 4-4: Aerator Calibration Curves Power Setting -e--e-- 0 -From visual observations it could be seen that the lower flow rates (around 0.5 L/sec) imparted virtually no mixing on the surrounding water, with the flow fields being very slow moving and virtually laminar. The highest flow rates exhibited a great degree of mixing with 4.0 L/sec being almost double the higher flow rate target previously set in Table 4-1. It was concluded that the aerator could provide quite a wide range of mixing conditions and would be adequate for agitating flow in the model basin. 4.4 Still-Basin Aerator Testing In order to examine the impacts of the aerator on the surrounding water body, the aerator was placed into the tank at the center and operated in water without any bulk flow. Prior to data collection, the aerator was activated and allowed to run for 30 to 40 minutes to allow the basin to reach steady state conditions. The ADV probe was set so that the sample volume was along the normal line that ran perpendicular to the cross bridge and the center of the aerator. Care was Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 42 taken to ensure that the A D V was properly oriented (i.e. the axes of the instrument corresponded to the three principal directions of the basin). Velocity measurements were obtained by the A D V at a frequency of 25 Hz for a specified period of time. Velocity measurements were made in a grid-like pattern in a plane extending radially from the aerator and over the vertical depth of the water. Figure 4-5 illustrates the location of the sampling plane with respect to the basin and the aerator. Samples were taken at 20 cm increments horizontally from the aerator from 50 cm distance to 190 cm distance and at 5 cm increments vertically from 4 cm from the surface of the water along the depth to a minimum depth of 7 cm. These limits on measurement were constrained by the physical shape of the tank and the aerator. A distance 50 cm from the aerator was only 30 cm from the edge of the aerator itself just beyond the splash zone. The A D V was found to not provide accurate velocity measurements when closer than 3 cm to the water surface or the basin floor. Velocity data were collected at each location for 7 minutes to ensure that the larger scale velocity fluctuations were captured in the sample time series. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 43 Plan l—Sampling P l a n e Profile £1 TBI •Sampling P l a n e z 4 \ ->x Figure 4-5: Velocity Measurement Sampling Plane 4.4.1 Determination of Sampling Time Some questions existed regarding the appropriate length of time required for sampling with the A D V . Long sampling runs were performed in order to establish the required sampling duration to capture the mean flow and turbulence parameters with a point sample. Samples were taken for fifteen minutes in duration and were collected 11 cm from the surface of the water at 35, 70 and 90 cm from the center of the aerator on the sampling plane. The velocity component considered was that of the x-direction, or the horizontal direction in line with the sampling plane. This component was chosen because it exhibited the largest velocity magnitudes at the sampling positions chosen. The samples were analyzed statistically using the t-test for mean similarity. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 44 The technique involved splitting each sample into two equal parts and examining the population distributions. The t-test will yield a probability statistic that indicates the probability that the two samples originated from the same population. As the sample sampling period increase one would expect the probability to increase as a greater number of the turbulent fluctuations would be captured. Figure 4-6, below, shows that this was not generally observed. The duration of the sample in Figure 4-6 represents the size of each velocity sample that was compared and the P statistic represents the probability that the two samples are of the same population. The target confidence limit threshold of P = 95% was never maintained with an increase in sampling time. There seemed to be a great deal of chaotic flow in the system on scales larger than those that could be captured with 15 minutes of sampling. It was determined that the largest reasonable sampling time for data acquisition was 7 minutes as larger sampling times would make obtaining a reasonably detailed grid of data difficult. This 7 minute duration was the sampling time used in subsequent analysis. 0 200 400 600 Duration of Each Population Sample (sec) Figure 4-6: Sampling Time Population Comparisons 4.4.2 Profile Data Analysis The data collected in the still basin tests were analyzed to capture the flow patterns in the basin and to assess how these flow patterns differed at various aerator flow conditions. The coordinate system employed here consisted of the following: x was defined as the horizontal direction, increasing away from the aerator; z was defined as the vertical direction, increasing Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 45 with increasing height; and y was defined as the horizontal direction perpendicular to both x and z in accordance with the 'right hand rule'. The orientations can be found in Figure 4-5, above. Velocities were defined as: u=dx/dt; v=dy/dt; w-dz/dt. Data analysis included the calculation of flow fields as velocities in three dimensions. Three separate runs were completed with the aerator configuration illustrated in Table 4-4, below. The configuration was identical to that of the typical configuration of the aerator calibration tests except that the depth of water was changed and the position A was adjusted to produce an H value identical to the standard configuration (Table 4-4). The power was varied for each of the runs. Table 4-4: Aerator Configuration for Profiling Aerator Configuration Value or Parameter Range A 49 cm B 20 cm C 5 cm D 42 cm P 5,6,7 H 2 cm I 11 cm Velocities were collected in time-averaged series and for analysis, these data series were separated into average and fluctuating, or turbulent components. Equation 4-5 below illustrates the separation of the velocity time series u into a time averaged and turbulent component. u{f) = u +u\t) (4-5) The velocity data were used to generate turbulence data, in the form of turbulent kinetic energy (k) defined by Equation 4-6 below (Rodi, 1980). k = (ua +v'2 +w'2)/2 (4-6) In addition to the analysis of the turbulent fluctuations in the flow, the velocity data were analyzed for average bulk flow. The mean velocity components were collected and analyzed to Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 46 see if any overall flow trends could be elucidated for each aerator configuration. Horizontal flow profiles (i.e. in the x-direction) were of particular interest and the patterns were expected to be similar to those presented by NCASI (1971) for full-scale aerators. The horizontal velocity profiles presented by NCASI exhibited high velocities away from the aerator at the water surface over a shallow depth with return flow to the aerator over a larger portion of the depth with lower velocities. The velocity magnitudes in both directions decreased with an increased distance from the aerators. 4.4.3 Profile Results Figure 4-7, Figure 4-8 and Figure 4-9 show three u-velocity plots at various distances from the center of the aerator. The plots show the in-plane horizontal velocity at the points measured against a dimensionless depth which is the height of the measuring point over the total water depth. The water depth was 42 cm in all runs. The configurations are identical except that the aerator flows were set to 3.1 L/sec, 2.7 L/sec or 2.3 L/sec. Although the profiles are quite different in form, a velocity field at the surface moving away from the aerator can be seen in all three cases with a recirculation zone of some shape beneath it. It can be seen that at some positions, all measured velocities move in a single direction over the depth (e.g. P=7 at distance of 50 cm). This observation seems to indicate that a flow mass-balance can not be resolved in the sampling plane using the sampling procedure employed. This could have been due to the presence of three-dimensional flows, inadequate sampling duration or too coarse a spatial sampling resolution. The inability to resolve flow fields may also have been a result of chaotic mixing patterns with the flow fields changing as the sampling progressed. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 47 P=7; 0=3.1 Usee 1.00 — I U Velocity (m/sec) 1.00 — | U Velocity (m/s) Figure 4-7: U Velocity Profiles - Power 7, Q=3.1 L/sec Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 48 U Velocity (m/s) Figure 4-8: U Velocity Profiles - Power 6, Q= 2.7 L/sec Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 49 P=5; 0=2.3 L/sec 1.00 - i 0.0 0.00 0.04 -0.04 0.00 0.04 -0.04 0.00 0.04 -0.04 0.00 0.0 x=130 x=150 cm x=170 x=190 U Velocity (mis) Figure 4-9: U Velocity Profiles - Power 5, Q=2.3 L/sec Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 50 Figure 4-10 shows a sample of turbulence data calculated from the same data presented in Figure 4-7, Figure 4-8, and Figure 4-9. These data show that the turbulence is predominantly a surface phenomenon and that the circulation zone matches the horizontal extent of the turbulent plume. A relationship between the size of the mixing zone and the aerator flow rate could not be clearly observed. Also, the intensity of turbulence in the flow field was very highly dependent on the flow through the aerator as seen in Figure 4-10. The profiles themselves seem to generally match those presented by NCASI (1971), in that the highest velocities are evidenced at the surface with an velocity magnitude decay with increasing distance and a bottom recirculation zone. P=7, Q=3.1 l/sec £ | 35.00 JS 30.00 •2 25.00 Z 20.00 re o tE > 35.00 30.00 25.00 20.00 60.00 80.00 P=6, 0=2.7 l/sec 100.00 120.00 140.00 160.00 180.00 60.00 80.00 100.00 120.00 140.00 160.00 180.00 200.00 P=5, Q=2.3 l/sec 60.00 80.00 100.00 120.00 140.00 160.00 180.00 200.00 Distance from Aerator (cm) Figure 4-10: Turbulent Kinetic Energy (k) Profiles K (m2/s2) 5.00E-004 4.50E-004 4.00E-004 3.50E-004 3.00E-004 2.50E-004 2.00E-004 1.50E-004 1.00E-004 5.00E-005 0.0OE+000 4.5 Depth-Integrated Trends In an effort to incorporate the aerator data into a depth-averaged 2D CFD model, the parameters of interest were examined by integrating over the depth and examining the trends. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 51 An attempt was made to approximate the dispersion number as a function of distance from the aerator. The technique used was one outlined by Fisher et al. (1979) for dispersion in isotropic turbulence. The technique requires the calculation of the turbulent time scale for the turbulent covariance in each of the principle directions. D = u'2^R(t)dT (2-9) where, 1 " — y^(u(t)-U~)(u(t + T)-u) R(T) = !LJ=± . (4-7) and, u' = turbulent fluctuation in velocity, R(r) = autocorrelation function, and T = autocorrelation time increment. The autocorrelation function was determined using an iterative technique illustrated in Visual Basic code presented in Appendix D. Figure 4-11 below illustrates a typical output of the autocorrelation routine. In the interest of computational efficiency, the autocorrelation function was only calculated 10% farther than the value of rthat resulted in a negative autocorrelation value. The area under the plot was calculated by incrementally summing the area under the plot until the negative autocorrelation point was reached. For example, in the case with the y-component in Figure 4-11, the area was calculated from rvalues from 0 to approximately 11 seconds. The area under the autocorrelation curve is known as the turbulent time scale (Fisher et al., 1979). Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 52 1.00 — i 0 10 20 30 40 Tau (sec) Figure 4-11: Sample Autocorrelation Output Prior to depth-averaging the values of the turbulent time scale and dispersion, the average of the three directional components was determined. The averaging of the components would provide a parameter of use in CFD modeling. It is this average value that was, in turn, averaged over the depth of the water column to produce the plots that follow in Figure 4-12. The results were difficult to interpret. The turbulent kinetic energy depth-averaged profiles seem to make reasonable sense, with the highest value profiles being associated with the greatest flow rate through the aerator and with the profile being generally higher near the aerator, than away from it. The turbulence time scale data are less straightforward. The 2.7 and 2.3 L/sec profiles show interesting trends with gradual increases in time step with a sudden drop and a further gradual increase. The higher velocity, 3.1 L/sec profile exhibits a much lower average turbulence time scale with a markedly different profile shape. The dispersion number profiles are products of the turbulent time scales and turbulent correlations and as a result, reflect the differences discussed in the previous two plots. As seen in Figure 4-12, the dispersion profile intensities are not clearly related to the aerator flow rate by this method. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 53 Depth Averaged Turbulence Values $ 2000 —i S 2.00E-4 — , C N 1.50E-4 — | ^ 1.00E-4 — \ cu 5.00E-5 0.00E+0 40.0 80.0 120.0 160.0 200.0 Distance From Aerator Centre (cm) Figure 4-12: Depth-Averaged Turbulence Values 240.0 4.6 Conclusions A laboratory scale basin was constructed, as was a small-scale mechanical aerator with the flow rate of the aerator approximating flow rates in existing ASB systems based on reactor volume. Exact similarity between the small scale aerator and one that is full size could not be accomplished due to the difficulty in accounting for splashing, but the mixing itself was found to be somewhat similar in the profiles that were formed, to those presented by NCASI (1971). It can be concluded from this aerator investigation that both the velocity magnitudes and the degree of turbulence (as measured by the kinetic energy of turbulence) are related to the flow rate through the aerator, as was expected. However, it has been shown that a radial sampling procedure seems insufficient to properly capture the flow regime (i.e. satisfy continuity). The flows seem highly three-dimensional, making the usefulness of a radial section relatively limited. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 54 The procedure for translation of the turbulent data profiles to depth-averaged dispersion profiles was unavailing. The expected trends in profile shape (higher dispersion near the aerator) were somewhat evidenced but a relationship to aerator power could not be rationalized. 4.7 Recommendations The incompleteness of the stream flow plots developed to date warrants an improvement in the experimental technique used. Although it may be difficult to acquire data any closer to the water surface, more data could be acquired closer to the bottom of the basin. Additionally, more data sets are required to fully explain and characterize the relationships between flow intensity and the resulting circulation and turbulent patterns. Specifically, one could perform detailed three-dimensional velocity sampling over a small volume to determine i f continuity could indeed be resolved in this basin using data from the A D V . The results presented only consider approximate dispersion number values based on the turbulent correlations. However, the flow shear is also expected to contribute to the mixing imparted by the aerator. The flow fields would have to be more carefully resolved (i.e. continuity accounted for) before the shear mixing could be evaluated. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 55 5 Tracer Studies Tracer studies were performed on the basin to determine the hydraulic characteristics and to see i f these characteristics could be modeled using CFD techniques. Several steps were required to prepare the basin for tracer studies and to conduct and interpret the studies themselves. This section is divided into several parts that address these issues in order: • basin modification; • rhodamine calibration curve; • tracer addition and measurement techniques; • results and data interpretation; and • recommendations. 5.1 Basin Modification In order to prepare the basin for tracer studies, the basin had to be adjusted from a still basin configuration, in which studies were performed on the aerators alone, to a flow-through configuration. In order to do this, two baffles were introduced to conduct the flow into and out of the basin. At the end of the outlet baffle section, a fixed weir was added that was 42 cm high from the floor of the basin. This inlet/outlet configuration was decided upon because of its potential to produce flow patterns that are not laterally uniform and at the same time to make use of the entire flow space and to minimize short-circuiting. This weir was used to control the height of the water in the basin. Beyond the weir a hole was cut in the side of the basin wall to allow water to drain out. Water was fed into the basin by a large hose that was attached to the city water supply. The flow rate was adjusted by a valve at the terminus of the hose. Flow rate was measured by using a 25 liter bucket that was calibrated to the liter and with a stopwatch to time the flow duration. Three measurements were taken and timed prior to the start of the tracer run and averaged to determine the flow rate. The issue of flow rate changes due to pressure fluctuations in the distribution system was a concern. Flow measurements were taken after the completion of some runs and it was determined that, i f flow rates changed at all, the changes were small and within the error of measurement. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 56 Plan 1 l/2f-t [0.46n]-l f t [0.30n] •15ft [4.57n]-•Water Source A e r a t o r Weir -•1 l /2 f t [0.46n] tz: I f t [0.30n] 12ft [3.6Sn] r 2 f t [0.61n] i Profile tt 1 3 / 8 f t [0.42m] Figure 5-1: Modified basin for Follow-through Conditions 5.2 Rhodamine Calibration Curve The tracer used in this study was Rhodamine WT. The dye was supplied in large barrels of Rhodamine mixed with methanol in a 60/40 Rhodamine/methanol ratio. It was unclear from the labeling i f the dye was mixed by weight or by volume so the true concentration of Rhodamine could never be established with certainty. It is for this reason that the concentrations of Rhodamine are referred to as dilutions from the original concentration. This does not effect the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 57 results of these dye studies as the absolute concentrations are not as important as the relative concentrations in tracer study analysis (Levenspiel, 1962). A set of seven Rhodamine dilutions was created to generate a concentration curve for use in interpreting the results from the tracer studies in the basin. The analysis of the Rhodamine concentrations was conducted with a H A C H DR/2010 spectrophotometer by transferring all samples into the same sampling tube, oriented the same way within the spectrophotometer. 5.2.1 Optimum Wavelength The optimum absorbency wavelength was determined to be 555 nm based on the absorbance response of a 40 000 times diluted sample and using the H A C H Spectrophotometer. The absorbance response based on wavelength can be found in Figure 5-2 below. A l l spectrophotometry was performed at 555 nm when interpreting the tracer study data. 490 500 510 520 530 540 550 560 570 580 590 Wavelength (nm) Figure 5-2: Optimum Spectrophotometer Wavelength 5.2.2 Concentration Calibration Curve In order to determine the concentrations of the Rhodamine for the tracer studies it was necessary to know the concentrations (or dilutions) of the Rhodamine as determined by the spectrophotometer. To determine the response of the spectrophotometer, a series of diluted samples was generated over the range of dilution from 1.6xl0"7 to 4x10"4. The response curve of the spectrophotometer was essentially linear over that range of concentrations with a R value of the linear fit of 0.998. The response curve is shown in Figure 5-3 below in a log-log plot. As the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 58 actual concentration was not important in the interpretation of concentration curves, and since the calibration response was almost perfectly linear, it was determined that all calculations could be performed in the more convenient units of absorbance as registered by the spectrophotometer. A l l samples were measured at room temperature. 2.0E+0 'E d) o c ro . a \ o W _Q <^ D ) O O.OE+0 -2.0E+0 -4.0E+0 -8.00E+0 -6.00E+0 Rhodamine WT-Log(Fraction of Original Concentration) Figure 5 - 3 : Rhodamine-WT Calibration Curve 4.00E+0 5.3 Tracer Addition and Sampling Routine The performance of the tracer studies was done with as much adherence to the NCASI procedures as was possible under the circumstances (NCASI, 1983). Flow rates through the basin were set to 2, 3, or 4 L/sec, because these flow rates produced hydraulic retention times that were reasonable for experimentation and were fast enough to induce turbulent flow in the basin. Flows in excess of 4 L/sec were very difficult to accurately measure. The mean hydraulic retention times and the entrance or exit Reynolds numbers are included in Table 5-1 below for each of the flow rates, using the width of the entrance as the characteristic length. As it was mentioned by Lyn and Rodi (1990), the need to maintain a degree of turbulence within the basin was met by the Reynolds numbers presented in Table 5-1. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 59 Table 5-1: Flow Rate Characteristics Flow Rate H R T Reynolds (L/sec) (min) Number 2.0 58 4800 3.0 39 7100 4.0 29 9500 Care was taken to ensure that, after the flow had been measured, the basin was allowed to reach hydraulic steady state. There was no straightforward way of measuring this, and previous literature reports were not clear as to how long hydraulic stabilization should take. It was decided that at least two hydraulic retention times would be allowed before any sampling was started. This was adhered to in each of the studies and often more time was allowed. Tracer was diluted and measured prior to addition to the basin. The NCASI guidelines suggest that the dye, when added, should be of similar consistency to that of the basin's fluid to avoid stratification. Additionally, Rhodamine WT, when undiluted, is highly viscous and tends to adhere to surfaces, making precise measurements of the amount of dye added very difficult. It was decided to dilute the Rhodamine 1:40 v/v. This degree of dilution provided a reasonably concentrated mixture for addition to the basin but also allowed the mixture to flow more easily and completely from the cylinder to the basin water. The dosing of the tracer was made such that the estimated peak of the concentration would be approximately 100 times the detection limit according to the NCASI guidelines. If the detection can be considered to be within the order of the spectrophotometer's incremental resolution (0.001 absorbance units) the approximate dosage of a 40 times dilution was about 3 liters assuming complete and immediate dispersion. This was found to be far in excess of the dosage required for the accurate generation of an RTD curve in this basin, as a good degree of short-circuiting was often observed resulting in relatively high peaks. The dosage amounts varied but were in the range of 100 mL to 500 mL of 1:40 Rhodamine-WT. The tracer was measured out using a graduated cylinder and was poured into the basin at about 10 cm upstream of the inlet to the basin. The flow fields were typically rather quiescent in this region, and certainly did not provide adequate mixing. To rectify this situation, immediately Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 60 after the dye was poured, the region containing the dye was quickly, mixed by plunging a mixing paddle vertically into the water and then out. This allowed for a reasonable distribution of the dye at the inlet. This sort of technique is in accordance with the NCASI recommendations that suggest mixing energy should be added i f there is none available at the inlet. Tracer sample collection was only performed at the outlet of the basin. The samples were collected with a syringe with a long tube attached. The tube was attached to a pole, one end of which was submerged into the outlet flow as the syringe was pulled. Attempts were made to extract water over the entire depth of the water column when the samples were taken in order to approximate accurately the 2-dimensional data required by the CFD model. A l l samples were collected at the very center of the 0.3m (1ft) wide outlet. Sampling intervals were staggered so that the greatest density of sampling was concentrated over the peak of the run, with a gradual decrease in frequency over time, as recommended by the NCASI guidelines. The flow through the aerator was determined using the calibration curve shown previously in Figure 4-4. A modified version of that graph was generated, as the aerator was always configured so that the head difference between the top of the draught tube and the water was 2 cm. Figure 5-4 illustrates the relationship between power level and flow through the aerator. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 61 5.4 Tracer Data Interpretation Tracer studies were performed in the basin with and without aeration at three different basin flow-through rates: 2, 3 and 4 L/s. The aerator power rates were also varied over the range allowed by the rheostat to determine the impact of the aerator on the tracer curves. Some runs were duplicated to examine the reproducibility of the tests. A l l raw data are tabulated and presented in Appendix E. A summary table is shown in Table 5-2 below. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 62 Table 5-2: Tracer Study Summary Table Run Basin Flow Aerator Flow Date Number Rate (L/sec) Rate (L/sec) 1 2.1 0.0 May-25-99 2 3.0 0.0 May-26-99 3 3.9 0.0 May-30-99 4 4.1 3.6 Jim-1-99 5 1.9 0.0 Jun-2-99 6 4.0 2.3 Jun-9-99 7 4.1 0.0 Jun-24-99 8 4.0 3.6 Jun-27-99 There is a series of parameters that can be derived from a tracer curve that were analyzed for comparison. The parameters calculated here are tabulated in Table 5-3 below. Table 5-3: RTD Parameters Parameter Description tbar mean time - centroid of the RTD curve tio time for 10% of tracer to pass t90 time for 90% of tracer to pass ho median time - time for 50% of tracer to pass cP the maximum concentration of the RTD tp the time to peak for the RTD MI Morrill Index d dimensionless Dispersion number Cr the fractional ratio of the dye that is recovered form the tracer study to the total dye added The dye recovery was calculated by summing up the product of the concentration at a sampling point by a time increment. The time increment was calculated based on equation 5-1 below. ' 2 2 2 The cumulative curve or F-curve was calculated incrementally as the dye recovery was calculated. The time-to-concentration values (tx) were calculated by examining the cumulative Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 63 curve for the two sampling points that straddled the target concentration. The time was determined by interpolating linearly between the two times. The calculation of the peak concentrations and times, Cp and tp, was done by simply determining the sampling point that exhibited the highest concentration and recording the concentration and the time (no smoothing or extrapolation was performed on the RTD data). The Morrill Index (MI) was calculated as the ratio tgo/tjo-The calculation of the dimensionless dispersion number was made according to the Marecos Do Monte and Mara method (1987) using the following equation with modified nomenclature. „ 2 1=2 (7 a -2 rc ,+c w ] 2 2 (2-27) e C t 1 lbar From the variance calculated with Equation 2-27, the dispersion number for a vessel with no dispersion at the inlet or outlet (i.e. closed reactor) can be calculated by trial and error by the Equation 5-2 provided by Levenspiel (1962). a] =2d-2{d)2 (l-e-Vd) (5-2) A reactor with an inlet and outlet open to dispersion (i.e. open reactor) has a relationship between the RTD variance and dispersion number defined by Equation 5-3, below (Levenspiel, 1962). a]=2d (5-3) The dispersion number can be calculated by trial and error or by using the curve shown in Figure 5-5, below. The solid line in Figure 5-5 represents Equation 5-2, or the reactor closed to dispersion at the inlet and outlet. The dashed line represents the dispersion relationship in Equation 5-3 for a reactor with an inlet and outlet open to dispersion, for comparison. It can be seen that at low dispersion the open and closed reactor equations are virtually identical. To remain consistent with the literature, the values of the parameters were calculated only for the first 2 HRT in the basin, as recommended by Mangelson and Waters (1972). This has a Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 64 significant impact on many parameters, especially those related to the variance of the RTD, as it tends to significantly reduce the variance. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Dimensionless Variance Figure 5-5: Dispersion Number - Variance Curve All parameters were calculated using a VBA macro developed in Microsoft EXCEL. This macro is shown in its entirety in Appendix F. Also, the results for each run, including the parameter calculation results are shown in Appendix E with the raw data. 5.5 Tracer Study Results 5.5.1 General Observations It was generally found that the shape of the tracer C-curves changed appropriately with increases in flow rate and aerator power, (i.e. increased mixing produced more CSTR-like curves). Visual observations made of the basin while the tracer studies were performed seemed to corroborate these findings. The aerators had a very significant impact on the flow fields, drawing tracer to the center of the aerator almost immediately. The slow moving regions of the tank tended to hold the tracer for longer periods of time, only gradually releasing the tracer into the faster moving flow path. Differences between various basin flow-through rates could not be Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 65 easily determined by visual inspection. The recirculation of tracer was very much evident, with the majority of the tracer being discharged after the first pass of the outlet, leaving the remainder to circle around and be discharged at the second pass. By the time the third complete circulation was reached, the tank usually appeared to be completely mixed. Some complications arose in the non-aerated runs, in that the lack of mixing often lead to streaks of dye that would stream through the outlet. Since the sampling was regular and samples were only taken at the prescribed time at the center of the outlet region, large amounts of dye could be missed entirely, or conversely, a streak could be captured, resulting in an unusually high concentration sample. The streaking of the dye was only observed at the time of the initial peak and was generally more pronounced in the slower moving, quiescent runs (i.e. 2 L/sec). This variation can be observed in the peaking in the runs without aeration. 5.5.2 Reproducibility In an attempt to assess the reproducibility of the tracer technique, several tracer studies were repeated under identical conditions. The results are presented here with the C-curves normalized for 100% mass recovery after 2 HRTs. For the case without aeration, shown in Figure 5-6, the reproducibility was quite good, with fairly good matching of peak times near 0.2. The peak values between the two runs did not match well with the Cp value for June 24 being almost twice that of the May 20 test. Additionally, the recirculation pattern indicated by the peak between 0=0.6 and 1.0 was captured in both runs. The reproducibility indicated by the tracer studies with aeration was also satisfactory and the results are shown in Figure 5-7. Both the June 1 and June 27 tracer profiles follow approximately the same shape. The peak dimensionless concentration values are both very near to 1.5 although the time to peak is quite different for each run. However, an examination of Figure 5-6 shows that, given the scatter in the RTD, the real peak location cannot be certainly known. Mangelson and Waters (1972) found that the specific timing of peaks and the appearance of the dye were not as reproducible as the parameters that represented the bulk shape of the RTD functions, and this observation seems to be supported by the present study. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 66 8.0 7.0 6.0 5.0 o O 4.0 o 3.0 2.0 1.0 0.0 Power=0, Flow=4 L/sec — | — June 24 ,1999 • © - - May 30, 1999 n r 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 t/HRT Figure 5-6: Tracer Study Reproducibility - No Mixing Power =8, Flow=4 Usee June 1, 1999 - © - - June 27, 1999 Figure 5-7: Tracer Study Reproducibility - Maximum Mixing Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 67 5.5.3 Parameter Trends Some interesting trends appeared in the analysis of the parameters derived from the tracer data collected. In the analysis of the effects of aerator power (measured as flow rate through the aerator) one can see interesting trends as seen in Figure 5-8. The peak concentration decreased with increasing aerator power; this was expected as mixing increased. Also, the mean and median times increased with an increase in aerator power. This phenomenon could be attributed to increased mixing in a system that would otherwise exhibit short-circuiting. A n interesting observation is that the dispersion numbers and variance of the distribution dropped with an increase in aerator power, as did the Morrill index, which is also a measure of degree of mixing. This is difficult to explain but may have much to do with the recirculation zone and the periodic re-occurrence of tracer in the C-curve. These re-occurrences of tracer will have a significant effect on the variance, as they occur some time after the first peak has already passed. This could explain the high dispersion numbers exhibited in those tracer studies with no aeration. Flow Through Rat* • 4 U s * o £ 0.2( E 0.40 — | 1 0.00 -0.0 1.0 2.0 30 4.0 00 1.0 2.0 3.0 Aerator R o w Rate (L/sec) Aerator Flow Rate (L/sec) Figure 5-8: R T D Parameters Based on Various Aerator Power Settings Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 68 No Aeration "2 10 -& a 5 -5, 0.11 # 5-0 -<3 0.60 -P K X £ , | 0.50 -s 0.40 -0.40 -P i E 0,30 -Flow Rate (L/sec) Figure 5-9: RTD Parameters for Various Flow-through Rates without Aeration The effect of varying flow-through rate on the parameters illustrates some interesting trends shown in Figure 5-9. The peak concentration dropped with an increase in flow rate. This seems to be rational as an increase in the Reynolds number should result in an increase in turbulence at the inlet and outlet regions. However the dispersion and variance parameters also decreased with an increase in flow rate, which seems to indicate a reduction in mixing. In direct contrast, the Morrill Index increased with an increase in flow rate. 5.6 Conclusions Tracer studies were successfully completed with and without mixing and at different flow rates. From the limited data set produced the reproducibility of the profiles seemed to be better without aeration than with aeration. A larger number of samples is required to produce more confidence in these observed trends. What is clear from this analysis is that the flow regimes, particularly those occurring under non-aerated conditions, do not follow the typical axial dispersion model and the recirculation periods add complexity when dealing with the parameters that typically are of great importance Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 69 in design, namely the dimensionless dispersion number d. It is illustrated here that when the underlying assumptions of dispersed plug flow are violated (cycling flow with non-uniform mixing) the interpretation of the RTD profile using this model becomes difficult. 5.7 Recommendations A l l the parameters here were calculated using the data only to 2 HRTs as recommended by Mangelson and Waters (1972). However, attempts should be made to develop techniques to extend data beyond that of the tracer curve by assuming exponential decay. Then a virtually complete data set could be obtained for each tracer and this may result in different estimated values for dispersion numbers. The number of tracer studies performed in this analysis was small, and some of the trends gleaned from these studies were weak and require more data for validation. A larger family of curves may reveal outliers in the existing set and new trends may emerge. A further suggestion is to obtain samples from within the basin during a tracer study test. Knowledge of the tracer concentrations within the basin would be very useful when calibrating the computational transport model to be discussed in Section 7. A series of time-lapsed photographs from above could also serve in comparing laboratory data to CFD output. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 70 6 FLUENT Navier-Stokes Model In order for the tracer model to be effective, it was necessary to know the flow patterns within the basin on a very detailed scale. A computational fluid dynamics model was employed to reproduce the fluid flow profiles exhibited in the previous section regarding tracer studies (Section 5). The commercial CFD model FLUENT, installed on a PC was used to calculate the fluid flow patterns. F L U E N T employs a finite-volume method for solving the fluid flow equations over structured meshes with a variety of turbulence models available. This was not a rigorous study of the effective CFD modeling of the basin but was only done to determine i f the flow could be accurately approximated for use in an RTD profile-generating routine. This section is divided into several sub-sections: • the fluid flow equations and theory; • the calibration of F L U E N T results to laboratory velocity data through the adjustment of grid size, shape and modeling parameters; and • recommendations for future modeling. 6.1 Fluid Flow Theory and Equations When reduced to its fundamental components, computational fluid dynamics is primarily a family of techniques for solving systems of differential equations with due care towards the physics these equations describe. 6.1.1 Fundamental Equations Computational fluid dynamics is a collection of techniques to solve fluid flow equations using numerical techniques. Most fluid problems involve the solution of the Navier-Stokes equations which are non-linear and which require numerical solution techniques to arrive at approximate solutions. The Navier-Stokes equations represent the conservation of both momentum and mass. The two-dimensional, incompressible, laminar flow equations are shown below. Equation 6-1 below illustrates conservation of momentum in the x-direction, equation 6-2 shows conservation of momentum in the y-direction and equation 6-3 represents conservation of mass in two dimensions. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 71 du du du 1 dP u{ d2u d2u — + u — + v — = + — - + -dt dx dy p dx p^dx by (6-1) dv dv dv vu hv — dt dx dy 1 DP fit d2v d2v (6-2) p dx p\dx2 dy2 (6-3) dx dy where, P u v velocity in x direction (m/s), velocity in y direction (m/s), fluid pressure (kg/m s2), fl = fluid viscosity (kg/m s), p = fluid density (kg/m ), and x,y,t = position and time ordinates (m, m, s). The unknown variables are the velocity components and the pressure, leaving a system of three equations and three unknowns. 6.1.2 Turbulence Modeling The flow equations provided in the previous section are applicable for laminar flow conditions. However, the flow conditions in this study are not laminar, as turbulence plays a significant role in the fluid flow as shown by the Reynolds number calculations in Section 5.3. In order to account for the turbulence, a traditional statistical approach was taken in which the primary variables were defined in terms of their time-averaged means and fluctuations from the mean as described in equations 6-4 and 6-5. (6-4) At (6-5) where, = instantaneous value, Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 72 </> = time-averaged mean component value, and (/>' = fluctuation from mean. The introduction of this substitution into equations 6-1 and 6-2 yields the following equations (Rodi, 1980). du _ du _ du \-U r V dt dx dy ]_dP_+d_ p dx dx p du p dx •uV J d + — dy p du P •uV (6-6) dv _ dv _ dv vu V V — dt dx dy }_dP_+d_ p dx dx p dv p dx •Vu' + • dy p dv — vu P 5y (6-7) It can be seen in equations 6-6 and 6-7 that new terms have been introduced by this transformation and these are the correlations between the fluctuations in flow in both directions. These correlations introduce complications into the solution of the continuity equations because the correlations between fluctuations can only be known i f the state of the turbulence itself is known. This problem can be overcome through the introduction of the turbulent or eddy-viscosity term (jut) and the development of equations that describe the behavior of the eddy-viscosity within the modeled space (Rodi, 1980). — M, •uv = — P du ^dv dy dx (6-8) P where, vt = turbulent kinematic viscosity (m2/s), and pt = turbulent dynamic viscosity (kg/m-s). The eddy-viscosity term is not a fluid property, but rather, represents an effective viscosity associated with the fluid state. The eddy-viscosity term is expected to vary in space and time, depending on the flow regime. The effective viscosity employed in the solution of the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 73 momentum equations is the effective viscosity represented as the sum of the turbulent and molecular viscosities as illustrated in equation 6-10 below. Jueff=M + Mt (6-10) where, jieff = effective dynamic viscosity (N-s/m2). 6.1.3 The k-e Model The turbulence model used in this study was the popular k-e model. This model introduces two equations that allow for the solution of the system of Equations 6-3, 6-6, and 6-7 by coupling the eddy-viscosity term to the velocity and pressure terms. It introduces two new parameters that must be solved for: the turbulent kinetic energy (k) and the dissipation rate of k (£•). It should be noted that the k-s model assumes that the eddy viscosities are identical for all Reynolds stresses, that is, the eddy-viscosity is isotropic (Rodi, 1980). These equations are shown below in Equations 6-11, 6-12, 6-13, and 6-14. M,=CU-p 8 (6-11) dk dk dk d h w — + V — = — dt dx dy dx ( pt dk} d —— + — p • crk dx J dy p-ak dy + Gk-e (6-12) de de de d VU h V = dt dx dy dx ( M, de^ d + — dy M, de yp-(T£ dyj k k (6-13) ML P f 2 ,fdv) 2 (dv du^ 2 \ 2 + 2 — + — + — V [dy) {dx) Kdx dy) ) (6-14) where, Gk k £ Cu, C/„ C2„ Ok, <JE = turbulent production term, = turbulent kinetic energy (m2/s2), = dissipation rate (m/s), and = empirical constants. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 74 The values chosen for the empirical constants in this study were the standard k-e parameter values shown below in Table 6-1 (Rodi, 1980). Table 6-1: Standard k-e Parameter Values Parameter Value CU 0.09 Ci. 1.44 C2E 1.92 1.0 1.3 6.1.4 Species Transport The modeling of species transport is done by assuming that the dispersion due to turbulence is proportional to the kinematic eddy viscosity adjusted linearly by the turbulent Schmidt number (He et al., 1999). The Schmidt number can be defined as the ratio between the diffusion of momentum and the diffusion of a scalar quantity within the medium. A variety of turbulent Schmidt numbers are employed in modeling. The fluent modeling system assumes a turbulent Schmidt number of 0.7. Lyn and Rodi (1990) used a Schmidt number of unity when modeling settling tanks. He et al. (1999) reported the common use of 0.8 for jet-in-crossflow modeling but also recommended a lower value of 0.2 for best results in this type of system. The governing equation for turbulent transport in two dimensions are shown in equation 6-16 below. 8C dC dC p, + u — + V- -dt dx dy p-Sc where, d2C d2C + dx dy 2 J (6-15) C — the concentration of a species at a point (kg/m3), p = fluid density (kg/m3), pt = turbulent viscosity (N-s/m2) u,v = velocity components in the x- and y-directions respectively (m/s), and Sc = the turbulent Schmidt number. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 75 6.1.5 Boundary Conditions for Laminar and Turbulent Flows The modeling of the flow using FLUENT involved three distinct boundary types: inlet, outlet and wall boundary conditions. The importance lies in determining how the boundary conditions (BCs) effect the primary state variables (u,v,P,k and e). Fluid Properties The density and molecular viscosity of the fluid were required for proper modeling of the fluid flow. The F L U E N T default water density and viscosity were used. The values are tabulated in Table 6-2, below. Table 6-2: Fluid Properties Used in FLUENT Modeling Fluid Property Value Density (p) 1000 kg/m 3 Molecular Viscosity (jti) 1.3 x 10"3 N-s/m 2 Inlet Boundary Conditions The inlets were modeled as velocity controlled inlets. In FLUENT the velocities are prescribed for the entrance either as a constant (uniform inlet velocity profile) or a polynomial equation (non-uniform inlet velocity profile). In this study both were options were examined. The pressure values are not enforced at the inlet. The turbulent kinetic energy (k) and dissipation rates (e) may also be set at the inlet i f values are known. In this study no values for the turbulent parameters were specified at the inlet. Outlet Boundary Conditions The outlet was modeled as a zero pressure outlet. This modeling set a point of reference for the development of relative pressure profiles within the solution space. A l l other parameters are not enforced at the boundary but are set to whatever the upstream value is prescribed to be by the solution. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 76 Wall Boundary Conditions The walls were modeled as no-slip boundaries. This boundary condition stipulated, in both laminar and turbulent flows, that all velocity components at the wall were zero. The pressure gradient was also assumed to be zero at the wall. However, the modeling of the near-wall flow varied significantly with laminar and turbulent flow modeling For laminar flow modeling at the boundary, the shear stress needed to be calculated in order to satisfy the momentum equations (FLUENT). The calculation of the shear stress is equivalent to the velocity gradient normal to the wall multiplied by the fluid viscosity as shown in equation 6-16 below. du on wall Au M— (6-16) An where, Tw = shear stress against the wall (N), u = velocity component parallel to the wall (m/s), and n = perpendicular distance away from wall (m). For turbulent flows, the modeling used the law-of-the-wall model to determine the wall shear stress and turbulent parameters based on the flow profile near the wall. The equation for wall shear stress is shown below in equations 6-17 and 6-18. pC0u25k°p5KUp f*= ,V +/ ( 6 - 17 ) and, y+ = H M " y" (6-18) J" where - 2 / 2 kp = turbulent kinetic energy at point P (m /s ), Up = mean velocity of the fluid at point P (m/s), E = empirical constant (9.81), yp = distance from point P to wall (m), Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 77 K = von Karman's constant (0.42), and CM = k-e turbulent model constant (0.09). Equations 6-17 and 6-18 are used i f the value ofy+ is greater than 11.225, otherwise the FLUENT solver uses the shear stress calculated from equation 6-16. The point considered is the centroid of the near wall cell, and FLUENT solves for the shear stress explicitly after every iteration. Boundary conditions were also applied to the two additional parameters introduced during turbulence modeling. The boundary condition for the turbulent kinetic energy was a zero flux across the wall. The boundary conditions for the dissipation rate are calculated based on the theory that the production of k and its dissipation rate were equal within the control volumes next to the wall. This assumption provided an equation for the calculation of the dissipation rates at the boundary cells shown below in equation 6-19. (6-19) where, £ p = dissipation rate at point P, kp = turbulent kinetic energy at point P, and yp = distance from point P to wall. 6.1.6 Initial Conditions and Solution Convergence The solutions were always initially set to a zero condition for all parameters throughout the solution space. The solution to the systems of linear equations that were employed based on Navier-Stokes and turbulent equations were solved using the Line Gauss-Seidel (LGS) iterative technique. The complete system of equations that describe the solution space is typically too cumbersome to solve simultaneously. The LGS method solves the equations along one column or row of the solution grid at a time. A solution along a single line of the solution space is much more easily generated than over the entire two-dimensional solution space. The solver performs the Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 78 numerical solution for each row and column in the solution space and then updates the solution. The change in a cell's value from one full LGS iteration to the next is considered the residual update for that iteration. In theory, as the solution converges, this update should become smaller. The normalized residual is a lumped parameter that represents the magnitude of the update over the entire solution space. The normalized residual used in FLUENT can be calculated using Equation 6-20 and is essentially an average of the update residual over the solution space. A l , j ] (6-20) i i m a x . / max where, the normalized update residual, the variable update residual for cell (i,j), the number of columns in the solution space, and the number of rows in the solution space. The solution was advanced iteratively until every solution residual fell below the conversion threshold, in this case set within FLUENT to be 10"3. After this point it was assumed that convergence had been reached and computation ceased for that flow scenario. Figure 6-1 below shows a typical convergence pattern illustrating the drop in the normalized update residual for the pressure and velocity terms. In Figure 6-1 the P, U and V variables represent the average update residual over the solution space for the pressure, x-velocity and y-velocity terms, respectively. Ri Imax Jmax Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 79 Bent-4 Grid (150x120) - 4 L/sec - k-epsilon turb and parabolic inlet 1.00E+0 - = » P U T3 <D N - — V co o 1.00E-3 - J 1.00E-4 0 200 400 600 Iteration Figure 6-1: Sample Convergence History 6.2 Comparison of FLUENT to Laboratory Data In order to determine i f the FLUENT model was producing accurate flow predictions, laboratory flow data were collected using the A D V to determine i f the calculated flow fields could be considered valid. The data were collected at various points within the basin as illustrated in Figure 6-2 below, at points marked A through J. At each point, the velocity was sampled at 4 depths. The velocities were averaged and the result was considered to be the representative velocity at that point. Velocities were converted to U and V velocity components that corresponded to the FLUENT flow directions and compared to the velocities found in FLUENT grid solutions. Table 6-3 below shows the coordinates of each of the points with the distances being measured from the bottom left corner of the active basin at the inlet. The sampling points were primarily selected to lie on paths of highest velocity (A through G). The remainder of the sampling points (H through J) were added afterward to resolve uncertainties in the flow path. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 80 y / \ >x Figure 6-2: Velocity Sampling Points Within Basin Table 6-3: Velocity Sampling Points Sample X(m) Y(m) Point A : 0 0.15 B: 1.11 0.26 C: 2.21 0.26 D: 3.21 0.81 E: 3.71 1.76 F: 4.26 3.46 G: 2.21 3.46 H : 4.26 1.76 I: 4.26 0.61 J: 0.36 1.76 6.2.1 Grid Shape Although initially the flow was modeled using a structured uniform mesh, some experimentation was performed on a structured, but slightly non-uniform mesh. One mesh examined was one that produced a more dense mesh near the inlet and outlet and this mesh was called "Bent 4". A n illustrative image of this grid structure is shown in Figure 6-3. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 81 Figure 6-3: Bent 4 Gr id Structure It was found that both mesh types produced reasonable results in terms of matching the flow, but the best results were obtained with the non-uniform mesh. The comparison was performed on the 2 L/sec model run with both meshes at a resolution of 120x96 grid elements. Table 6-4 below shows the values of the root mean square of the error between the predicted and the measured velocities, comparing grid shape. Three elements were compared, the velocities in the U direction, the velocities in the V direction and the velocity magnitude. Figure 6-4 below illustrates the associated errors with the change in grid shape for a mesh of 120x96. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 82 Table 6-4: RMS Error by Mesh Type Measurement Root Mean Square Error (m/s) 120x96, 2 L/sec, k-e, uniform inlet Uniform Mesh Non-Uniform Mesh U Velocity 0.00559 V Velocity 0.00112 Velocity Magnitude 0.00509 0.00505 0.00073 0.00468 125x100 Grid Size - 2 L/sec - U Velocity 125x100 Grid Size - 2 Usee - V Velocity + Uniform Grid o Non-Uniform Grid Exact 0.006 a 0 .004 E U 0 .002 0.00 0.01 U Velocity Measured (m/s) Figure 6-4: Fluent Calibration - Grid Shape -f- Uniform Grid O* Non-Uniform Grid Exact - i | i | i | i | i | -0 .002 0.000 0.002 0.004 0.006 0.00S V Velocity Measured (m/s) It can be seen from Figure 6-4 that the greatest relative deviation from the measured values was observed at points B and C. It seems that the model overestimated the U components of the velocities at these points. It is believed that this could be due to the lack of consideration of bottom shear stress in the model. The inclusion of bottom shear stress, perhaps with the shallow water equations (Ye and McCorquodale, 1997), would dissipate some of the kinetic energy in the model. Additionally, there were no turbulence input parameters at the inlet. It was assumed that the kinetic energy and the dissipation rate were zero at the inlet. The inclusion of some appropriate turbulence parameters may dissipate the kinetic energy to turbulent energy more quickly and reduce the velocities just downstream of the inlet. 6.2.2 Turbulence Modeling and Inlet Configurations Many parameters can be varied when modeling using CFD techniques. Two issues of concern were the inclusion or exclusion of turbulence modeling (k-e in this case) and the flow Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 83 inlet boundary conditions. The modeling was performed with and without the k-e model and with a uniform velocity inlet and a parabolic inlet velocity boundary condition. It was found that convergence could not be reached with the laminar flow modeling. Analysis of the residual histories indicated that the model seemed to approach convergence and then diverged to another flow condition and then oscillated about this solution. Figure 6-4 shows the convergence history of the U , V and P parameters as defined in Figure 6-1, above. Bent-4 Grid - 4 L/sec - Laminar with Uniform Inlet 1.00E+0 CO 13 •g co CD or CD CO T3 Q_ Z) T3 CD N la E 1.00E-1 1.00E-2 1.00E-3 - d 1.00E-4 V 400 800 1200 Iteration Figure 6-5: Non-Convergent Laminar Residual History Even though convergence was not reached with the laminar model, it was thought that perhaps the solution that the laminar flow model approached (but never reached) was a closer approximation to the measured velocities than that with turbulence modeling and proper convergence. Figure 6-6 below shows the various flow fields generated with either the laminar or turbulent models and with different inlet configurations. It can be seen that the flow fields are all generally similar with obvious differences observed at the inlets and in the shape of the recirculation patterns. The tighter circulation pattern observed for the parabolic inlet with laminar flow modeling (with two circulation zones) was of concern because it was difficult to Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 84 ascertain from observation i f the flow fields in the actual basin circulated with such a pattern. Comparison of the flow fields with additional data measured in the basin (points H, I and J) were used to examine this possibility. Table 6-5 below lists the RMS error values for the various flow components by turbulence model and inlet configuration. Figure 6-7 below shows the relative accuracy at these points, compared to laboratory A D V data. Contours of Flow Magnitude (m/s) - Bent4 (150x120) - 4 l/sec at inlet Uniform Inlet Parabolic Inlet Turbulent K-Eps 0.50 1.00 1 50 2.00 2.30 3 00 3.50 4.00 4.50 090 1 00 1 90 2 00 2 50 3 00 3 90 4 00 4 50 4.60E-O02 m/s 4.00E-002 m/s 3.60E-002 m/s 3.00E-002 m/s 2.50E-002 m/s 2.00E-002 m/s 1.50 E-002 m/s 1.00E-002 m/s 5.00E-003 m/S 0.0OE+OO0 m/s Figure 6-6: Velocity Flow Fields - Comparison of Turbulence and Inlet Configurations Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 85 Table 6-5: RMS Error by Turbulence Modeling and Inlet Configurations Root Mean Square Error (m/s) 150x120 Bent4 Measurement Laminar, Uniform Turb, Para Turb, Uniform Laminar, Para (non converge) (non converge) U Velocity 0.0076 0.0054 0.0057 0.0095 V Velocity 0.0043 0.0036 0.0037 0.0073 Velocity Magnitude 0.0080 0.0063 0.0066 0.0103 150x120 Grid Size, Bent4 - 4 Usee - U Velocity 150x120 Grid Size, BenM - 4 L/sec - V Velocity + o • o Laminar, Uniform inlet Turbulent, Parabolic inlet Turbulent, Uniform inlet Laminar, Parabolic inlet Exact 0.02 0.04 J Velocity Measured (m/s) 0.028 0.024 0.020 — | 0.016 0.012 0.008 0.004 0.000 -0.004 -0.008 "T~ Laminar, Uniform inlet O Turbulent, Parabolic inlet O Turbulent, Uniform inlet A Laminar, Parabolic inlet Exact I ' I 1 I 1 I 1 I 1 I 1 I ' I 1 I -0.012 -0.008 -0.004 0.000 0.004 0.008 0.012 0.016 0.020 0.024 0.028 V Velocity Measured (m/s) Figure 6-7: Flow Calibration - Comparison of Turbulence and Inlet Configurations A review of the information for specific calibration points revealed some interesting aspects of the modeling. Samples taken at point A (at the inlet) showed that the inlet is more like a parabolic velocity profile than a uniform velocity profile, but that the best representation was somewhere in between.. Most likely, the flow profile was logarithmic in nature. The results also distinctly show that the parabolic inlet more accurately modeled the outlet as well (point F). It can also be seen that the turbulence modeling improved the overestimation of the velocities at points B and C, due to the dissipation of the momentum into turbulent structures. Also, the tight circulation observed with the laminar flow regime with a parabolic inlet was refuted with the observations at point J, where the V velocity in that case indicates flow in the opposite direction. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 86 6.2.3 Grid Convergence A principle of CFD modeling is that, with an increase in resolution, the modeled solution will converge onto an actual solution. This principle was examined by looking at the results of four grid sizes while modeling the flow with k-e turbulence modeling and with a parabolic flow inlet. The error results are shown in Table 6-6 and are graphically displayed in Figure 6-8. Table 6-6: RMS Error by Grid Resolution Root Mean Square Error (m/s) by Grid Resolution Bent4 Grid, 4 L/s, k-e turbulence modeling, Parabolic velocity inlet. Measurement 73x50 120x96 125x100 150x120 U Velocity 0.008327 0.006993 0.006964 0.005400 V Velocity 0.006690 0.007186 0.007159 0.003610 Velocity 0.009948 0.009486 0.009464 0.006279 Magnitude UJ Bent4 Grid- 4 Usee - Turbulence modeling -Parabolic inlet - U Velocity Bent4 Grid - 4 L/sec -Turbulence modeling -Parabolic inlet - V Velocity + 73x50 o 120x96 • 125x100 o 150x120 Exact E w ui > > + 73x50 o 120x96 • 125x100 o 150x120 Exact U Velocity Measured (m/s) Figure 6-8: Flow Calibration - Grid Resolution 0.000 0.008 0.016 V Velocity Measured (m/s) It can be seen from Table 6-6 that the increase in grid resolution produced a fairly consistent decrease in the root mean square error values. It can also be seen in Figure 6-8 that the highest resolution grid generally produced the values that were closest to the exact match line. However, it must be noted that the actual correct answer, as judged by the laboratory velocity readings, may not correspond to the converged solution as reached by the computational model. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 87 6.3 Conclusions It can be seen that the flow modeling with FLUENT was quite successful in predicting the general flow patterns observed in the basin without aeration, when modeled as a turbulent flow (using the standard k-e model) with a parabolic inlet and using a non-uniform structured mesh. The model results generally matched velocity data collected within the basin quite closely. Some notable discrepancies included the samples downstream of the inlet for which the model over-estimated velocity magnitudes and directly at the inlet, where the measured inlet profile was found to lie somewhere between a uniform and parabolic inlet. The first discrepancy can likely be attributed to the fact that no bottom shear stress was accounted for in the two-dimensional model and that the turbulent parameters at the inlet were unspecified in the model (i.e. k and e set to zero). 6.4 Recommendations It is recommended that more experimentation be done with the grid structure for modeling the flow using FLUENT. Only two grids were examined in this study and a variety of grid techniques could be used. Although FLUENT 4.0 (the package in use) could not resolve unstructured meshes, other solvers do and an investigation of unstructured meshes might produce more accurate results. It is clear from the analysis that a two-dimensional approximation introduces some errors, especially regarding the basin floor shear stress. Incorporating the shallow water equations (i.e. including a momentum sink associated with the bottom shear stress) may yield results more in line with the velocity measurements in the tank. As was mentioned above, no initial inlet turbulence parameters were included in the modeling. Some effort could be made to measure these parameters using the A D V and perhaps improve the performance of the model. Additionally, it was shown that the inlet profile lay somewhere between a parabolic profile and a uniform inlet. The inlet could be more accurately modeled with some detailed velocity measurements at the inlet or by assuming some intermediate profile shape. The overall match between the depth-averaged velocity measurements and the velocities predicted in the FLUENT model for the highest resolution on the Bent 4 grid, with k-s Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 88 turbulence model and parabolic inlet, was quite good. It was this final flow regime that was used in the CFD Species Transport Model discussed in the next section. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 89 7 CFD Transport Model In order to calculate the RTD profile of the basin, it was necessary to incorporate a species transport model into the CFD modeling. The licensed version of the FLUENT modeling package that was employed for the Navier-Stokes equation modeling was not equipped with a species transport model. Therefore, it was necessary to develop a species transport model employing CFD techniques that could use the solutions provided by the FLUENT modeling to predict RTDs. The species transport code developed in this study was collectively named TRANSPORT. The model was developed using Visual Basic and incorporated CFD techniques that will be discussed below. This chapter addresses the following aspects of the CFD species transport modeling: • the theory of the development of a discretized model; • the spatial discretization techniques employed and the accuracy associated with them; • the time discretization techniques employed and the associated accuracy; • the numerical solution techniques; and • the modeling verification and results. 7.1 CFD Transport Equations and Theory In this section, a summary development of the finite volume approach for discretizing the species transport equations is presented. When examining the transport of a chemical species, the equation that is typically used is the advection-dispersion equation. Equation 7-1 below describes the contaminant transport equation that describes Fickian diffusion and advection of a species. DC dC 8C d2C d2C (7-1) hW + V dt dx dy dx2 dy2 where, Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 90 C = the concentration of a species at a point (kg/m3), D = the dispersion value in both x- and y-directions (m2/s), and u,v = velocity components in the x-and y-directions respectively (m/s). Equation 7-1 can be re-written in vector notation _ + V - v C = J D ( V - V C ) (7-2) dt Integrating Equation 7-2 over the control volume yields: f —dV+ \y-vCdV = D f V - V C W (7-3) tv Qt lev JCV where, dV = the incremental control volume, and CV = control volume. Using the Gauss' theorem (Spiegel, 1997) which states that f V-FdV = { F-ndA, (7-4) JCV Jd(CV) where, « = a unit vector normal to the control volume F = flux vector dA = an incremental area surrounding the control volume. Equation 7-3 becomes — V + { vCndA = D<{ VC-ndA. (7-5) dt Mcv) Mcv) Equation 7-5 can be applied to a two-dimensional uniform mesh, where each cell in the mesh becomes a control volume. The position of the center of the cell is denoted by the integer positions i and j, which denote a cell's position in the x- and ^ -directions respectively. This being the case, the four boundaries for cell ij will be located at half increment from the center (i.e. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 91 i±Vz,j±Vz). Applying this definition to the control volume equation (7-5) will yield the following equation for a finite-volume formulation for cell (i,j): dt Ax Ay + (uCy^ y Ay + (vQft! j ',\6x = D dc i+\12,j 8C 1,7 + 1/2 Ay + — Ax 8X M / 2 , 7 ^ , - , 7 - 1 / 2 (7-6) or dC dt 1 (uC-DdS. dx Ax V / i - l / 2, Av f vC-D d£ 1,7+1 / 2 J ! ' ,7-W 2 (7-7) 7.2 Spatial Discretization Approaches and Accuracy Continuous functions can be discretized using Taylor expansion techniques to establish the order of accuracy. A standard accepted level of accuracy in most CFD systems is 2 n d order. This implies that as the discretized interval in space or time is reduced by a certain factor, the errors associated with the discretization will be reduced by the square of that factor. The Taylor expansion equations are shown below in Equations 7-8 and 7-9 (Spiegel, 1997): dC C(x0+s) = C(x0) + s dx + • :2 d2C 2 dx2 + • £3 d3C 6 dx3 + . (7-8) C(x0-£) = C(x0)-£ 5C dx £2 d2c £3 d3c j 2 dx2 6 dx3 + . (7-9) where, C = the dependant variable; x = the independent variable; and £ = a small shift in location from x0. 7.2.1 Dispersion Flux In order to determine a reasonable approximation of the first derivative of C with respect to x we consider a case where xo= i+Vi and e= Ax/2. Subtracting Equation 7-8 from equation 7-9 in this case will provide the following: Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 92 CM-C,=2 V * J dC dx + -2(AxYd3C 61 2 J dx3 (7-10) Rearranging equation 7-10 will yield dx l + l / 2 C -C Ax + 0 ( A c 2 ) (7-11) The first term on the right hand side of equation 7-11 is a centered difference approximation of the derivative on the left hand side. The second term on the right hand side, 0(Ax 2), represents the errors associated with the approximation and that they are "of the order Ax ". In other words, the scheme is 2 n d order accurate for the first derivative of C with respect to x. It is this discretization that was employed throughout the modeling of dispersion. 7.2.2 Advective Flux The determination of the advective flux discretization was not as straightforward as that of the dispersion flux. In order to determine the advective flux across the control volume, the actual concentration at the boundary must be determined. Two different discretization schemes were examined in this study: first order upwind and second order upwind. First Order Upwind Perhaps the simplest approach to determining the concentration at the interface between two control volumes is to assume that the interface concentration is equivalent to the concentration in the control volume immediately upstream (in the opposite direction of the velocity vector normal at the interface). Simple manipulation of the Taylor expansion will yield: + 0 (Ax 2 ) (7-12) This approximation is clearly only 1 s t order accurate and will suffer from fairly serious accuracy errors. C i"+l / 2 „ Ax dC = C,. + 2 dx Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 93 Second Order Upwind This technique builds upon the principle of the first order upwind in that it looks to the direction of the flow at the interface and interprets the interface concentration from the upstream concentration data. This formulation includes an additional upstream control volume beyond the one adjacent to the interface in calculating the interfacial concentration. This approximation can be considered a linear extrapolation to the interface point using the two upstream control volume centers to construct the line. The formulation is second order accurate as can be illustrated through Taylor series analysis (Equation 7-13). C, 3C ;—C ;_, 3 , 2d2C i+1/2 ' " . - - A x " — T 8 dx2 + 0(Ax 3 ) (7-13) (+1/2 This formulation shows increased accuracy over the first order upwind scheme and increased stability over the centered difference scheme, but it is susceptible to overshoots and undershoots that can produce unreasonable and inaccurate results. The issues relating to overshoots will be discussed below. 7.2.3 Semi-Discretized Equations The incorporation of the spatial discretization into the original finite volume formulation results in a semi-discrete form, so called because the time derivatives remain continuous. The application of the two different advective fluxes and the single dispersion flux will result in two distinct semi-discreet forms which are shown below. The velocity values at the boundaries are calculated with the identical approximation technique as the concentration values (although the velocities are static and not updated). First Order Upwind This equation assumes that the velocity flow is positive in both the x and y directions. dCu _ 1 dt Ax f fC -2C +C ^ Ax _L_ 4v U- C — u- • ,C , — D i,J «.7 1.7 -1 ' . 7 ~ ! f C -2C +C ^ ° i , 7 + l T ^ i , 7 ' - l (7-14) A^ Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 94 Second Order Upwind This equation assumes that the velocity flow is positive in both the x and y directions. dC 1 dt Ax f 3-uUCU -4ui-iJCi-U +ui-i,jCi-i.j D (r -2C +C ^ Ax J _ Ay 3'UUCU -4uU-iCi,j-i + u i , j - 2 C U - 2 D f C - 2C + C ^ (7-15) v Ay J) 7.3 Time Discretization and Accuracy In order to ensure the time accuracy of the transport equations, a discretization of the left hand side of equations 7-14 and 7-15 of 2 n d order was required. The trapezoidal discretization scheme is second order accurate and the equation is supplied below, d C C"+l-C" = a(Cn+l) + a(C") dt ~ < At 2 ) where the a() function represents a linear equation, in this case the right hand side of equation 7-14 or 7-15. Essentially, the trapezoidal approximation evaluates the update in the concentration for a cell as the average of the calculated updates from the current and future time steps. We can rearrange Equation 7-16 into a form that is more convenient to solve using numerical methods through the introduction of the update variable, SC. X u = Clj1-CiJ (7-17) This can be introduced to Equation 7-17 to produce a system of equations that can be approached by discovering a solution for the concentration update. After each time step the update can be added to the original solution matrix. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 95 SCu a(SCj) _ a{C?j) At (7-18) To illustrate the fully discretized equation implementation, the first order upwind with trapezoidal time discretization is presented below. SCj 1 • + -At • 2Ax + • 2Ay UijSCj -ui_iJSC-ij-Dl w,jSCij j^SCij-i -D r SCMj-lSCij + dCt-xj^ Ax 5dj+\ - 25d,j + SCij-i 4y J) f l 2Ax u- C —u , C , — D (C -2C + C ^ (7-19) Ax J) 1 2Ay u C —u- ,C- • , — D •J i,J i, 7-1 i , 7-1 fC -2C +C ^ Ay J) The other two discretizations can be derived in the same manner as with equation 7-19. 7.4 Numerical Solution Techniques Equation 7-19 above can be written in more compact way that facilitates the computational techniques employed. (1 + &{3x:.. + Atprij + Atax:UE_h0 + Ata^.E^ + Atyx, >0 + AtyyJJE0A J5CU = - (A//?,,,. + Atax,dE_h0 + Atyx..jEh0 )c£ - (At/3y. . + A / a , , / , , + Atyy:iJE0<+l where, EabSCi,j = SCi+a,j+b (7-21) P x : , J 2Ax Ax2 (7-22) axi i = 2 Ax 2 Ax (7-23) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 96 r « , - f £ 07-24) In Equation 7-20 the a, /?, and y factors can be derived from the constants through the rearrangements of Equation 7-19. The form of Equation 7-20 will produce a sparse system of equations of size (n x m) where n and m are the horizontal and vertical dimensions of the solution space grid. The form of the equation can be changed by using the approximate factorization approach that will allow for a two-stage solution of two tri-diagonal matrix equations. Through the principle of approximate factorization, Equation 7-20 can be rewritten as: (1 + + A* «„• ;£_ 1 > 0 + toY^jE^) x (1 + Atj5y,. + Ata^E^ + Atyy:iJE0M )5Cij = - f a P r j j + ^ ocx:iJE_h0 + Aty^E^j-^tp^ + Atay,JE0<_1 + Aty^E^fcj This equation is second order accurate with respect to the time advance scheme and therefore does not reduce the accuracy of the trapezoidal time advance scheme (Pozrikidis,1997). The approximate factorization scheme is applied in two steps as follows, written now in vector notation. The solution involves solving for an intermediate update in the solution and using the update as the flux integral in the second step. (1 + At[/Jx ] + At[ax ]E_h0 + At[yx ]£+I>0 )8C= ( ? 2 6 ) - {At[px ] + At[ax ]E_ifi + At[yx ]EU0 )C" - {td[fiy ] + At[ay + At[yy }E0M p (1 + At[/3y ] + At[ay ]E0_, + At[yy ]E0M )5C = 8C (7-27) where, 5 C = the intermediate concentration in the approximate factorization routine. Each step of this system can be solved using the Thomas algorithm solution for tri-diagonal matrices and the update is added to the solution matrix for every time advance. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 97 7.4.1 Boundary Conditions The boundary conditions (BCs) were applied using the ghost cell technique. This technique expands the solution space by including an additional "ghost" cell beyond each boundary. The values of this cell are manipulated in such a way as to influence the conditions interpolated to the real boundary. For example, in a uniform structured mesh, i f one were to require a zero velocity at a boundary, one would set the ghost cell value equal to the adjacent real cell but opposite in magnitude. In this way the interpolated value at the boundary would be zero. 7.4.2 Velocity Boundary Conditions The velocity boundary conditions were required only because the imported structured grid from F L U E N T did not contain final ghost cell values associated with the converged solution. As a consequence the ghost cells employed here were calculated to be identical to those used in the FLUENT solution technique. Inlet For the inlet it was assumed that all flows were axial and that there was no crossflow. That is, all the flow was in the x-direction (u velocity) and there was no y-component (v velocity) to the flow. The boundary conditions for the inlet were set based on the preset values used in FLUENT. The x-component of the velocity would be known from the flow BCs in fluent and were incorporated. The equations are shown below with the index of zero denoting the velocity of the ghost cell control volume and an index of one being the first control volume inside the boundary. The index "wall" denotes the value specified for velocity at the boundary: uo = 2 ( 7 - 2 8 ) v 0 = - v , (7-29) Outlet It was assumed that the outlet velocities at the boundary would be equal to whatever the real cell values were adjacent to the outlet. The equations describing these follow. u0 = w, (7-30) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 98 v 0 =v , (7-31) Walls It was assumed that there was a no flow condition at the wall boundaries. The equations are identical to the no velocity condition seen in Equation 7-29 and follow below. (7-32) (7-33) 7.4.3 Concentration Boundary Conditions: The concentration ghost cells were employed to control the flux of concentration across the walls. The advective fluxes were controlled by the velocity boundary conditions but the diffusive fluxes at the boundaries were dictated by the concentration gradients at the boundaries. It was desirable in this case to set a zero diffusive flux at all wall boundaries as well as the inlet and outlet boundaries, which could be enforced by making the ghost cell value equal to the value of the cell immediately within the boundary. This principle is illustrated in Equation 7-34 below. C„ = Cx (7-34) It was found that the concentration boundary conditions for the advective flux conflicted with the boundary conditions for velocity when attempts were made to simulate a tracer study at the inlet. This is primarily because one concentration boundary condition cannot enforce a no flux condition for advection and diffusion simultaneously with the velocity boundary conditions set to realistic values. This issue is dealt with later in this section. 7.4.4 Import of flow field from FLUENT The flow field solutions that were determined to be the most accurate were calculated using a non-uniform structured mesh (see Section 6.2). This data structure was not directly employable in the scalar transport routine developed here, as the CFD code for conservative transport was designed for uniform, structured data. The data were converted to a uniform mesh by means of linear interpolation with two-dimensional triangulation. This technique generated triangular Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 99 surfaces that joined known data points and interpolated grid points by applying them to these interpolated surfaces. The computer package SURFER was used to implement this conversion. The results of the surfer conversion was a structured grid of data. Some concerns existed as the interpolation scheme considered the source values as points, when in fact they were control volume averages in the case of the FLUENT file. Additionally, the values that were output from the SURFER interpolation routine were point values and not the control volume averages that would have been ideal. A more accurate solution to this problem was not developed but it was suspected that i f the TRANSPORT computational domain was of a similar or lesser resolution than the F L U E N T input file, then the interpolation scheme should have been adequate. 7.4.5 Initial Conditions for Tracer Study Simulation Attempts were made to emulate a standard tracer study with the TRANSPORT model. The initial conditions within the tank were a zero tracer concentration throughout except where the addition of the tracer was concerned at the inlet. Two methods were employed to approximate this real life situation of tracer addition: 1. pulse input of tracer at the inlet, and 2. constant input of tracer into the inlet. The pulse input of tracer was approximated by uniformly distributing 100 units of tracer (arbitrary measure) across the column of cells just downstream of the inlet boundary at time zero. The simulation was allowed to run and the average concentration across the outlet was measured resulting in a characteristic C-curve for the simulation. The constant input of tracer was modeled by setting the ghost cell values just upstream of the inlet to a constant value throughout the simulation. This would simulate a constant input of tracer and the concentration at the output would generate a characteristic F-curve, from which a C-curve may be derived. 7.4.6 Overshoot Management for Second Order Upwind It was mentioned earlier that with the second order upwind technique for advection modeling there would be overshoots and undershoots. A technique called the Van Leer Total Variation Diminishing (TVD) Scheme was used to dampen the overshoots and undershoots (Van Leer, Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 100 1974). The advective flux calculation at the edge of the control volume can be changed slightly to: ^ , + . / 2 = C f + ^ ( C , - C H ) (7-35) in which can be seen that i f y/= 1 then the flux is a standard second-order upwind flux. The value of ^ i s adjusted based on the following equations to limit extrapolation overshoots and undershoots. r + \r\ ^ ( r ) = T77 (7-36) r s r (7-37) Figure 7-1 below illustrates the increased accuracy of 2 n d order schemes over 1 s t order schemes and the effect of the Van Leer TVD scheme to reduce overshoots. 4.0 6.0 \ 8.0 10.0 Distance ^ / Figure 7-1: Comparison of Advective Flux Methods Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 101 7.4.7 Programming Flow The code of the computer program that calculated the time-accurate transport of the conservative species is included in detail in Appendix G. The program was written in Visual Basic and the modules including the code that perform the calculations detailed in the previous sections are included. The control form used to input parameters and specify input file names is not included. Figure 7-2, below illustrates the flow chart that describes the basic computational structure of the solute transport routine. Prior to computation, a flow data file is required that specifies the grid size and the control volume velocity components. In all cases, except for debugging the flow data, a file was imported from a FLUENT solution data file into a form that could be interpreted by the solute transport routine. The flow data file and initial conditions were specified prior to the computation. After starting the routine, the solution solver would then iterate through the approximate factorization timesteps for a number of timesteps specified. The output interval specified at the user input stage prior to computation indicated when the RTD concentration data should be written in memory. The RTD data were stored in two different time series: as a mean concentration across the outlet boundary and as the total scalar mass in the basin. These data were output at the end of the run for analysis. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 102 User Input (Dispersion, Iterations) 1 Calc a, P/ Y ulate Factors No Set Boundary Conditions Approx Factor Columns Approx Factor Rows Update Solution Figure 7-2: T R A N S P O R T Code Flowchart 7.5 Model Verification In order to verify that the conservative tracer transport model was working satisfactorily, the model was compared to analytical solutions of the advection-dispersion equation. The model was verified by setting the flow data to a uniform velocity field and setting a known global dispersion parameter with a variety of spatial and temporal resolutions. Results could be compared directly. It was found that the model did indeed match the expected results and did converge on the theoretical exact solutions. Included below are some illustrative runs that were performed to establish confidence in the transport model. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 103 7.5.1 Dispersion Figure 7-3 below shows the precise modeling of dispersion flux with no advection. It was performed at a grid resolution of Ax = 0.123 and run for a total time of 10 with a dispersion value of D = 0.05. A mass in the center control volume (distance = 5) of 100 was established for the initial condition. The exact analytical solution is shown in Figure 7-3 as the thick line and represents the modeling target for accuracy. It can be seen that increased time resolution increased the accuracy of the routine. This figure shows a slice of the model in the x-axis but the y-axis verification was equally accurate. Tests for conservation of mass over the domain were also performed and the dispersion model showed no production or destruction of mass. 50 — | 0 2 4 6 8 10 Distance Figure 7-3: Dispersion Model Accuracy 7.5.2 Advection Figure 7-4 below illustrates the grid convergence of the 2 n d order advective flux procedure with T V D overshoot control. Shown here is the advection of a step initially positioned at 50 and advected in a positive direction with a velocity of 0.1 and a Ax of 1 with no dispersion. It can be seen that the increase in temporal resolution improved the solution as indicated by the change in the C F L number. The C F L number is a measure of the speed of the fluid medium to the "speed of information" transport and is shown in Equation 7-38 below (Abbott, 1989). It can be seen Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 104 that numerical dispersion does play a role but the results do converge to the exact solution. Tests for conservation of mass over the domain were also performed and the advection model showed no production or destruction of mass for the advection of a plug of tracer. CFL = v- At Ax (7-38) 2nd Order Upwind Advection Routine c o o O 40 60 Distance Legend Exact CFL = 0.1 CFL = 0.25 CFL = 0.5 Figure 7-4: Advection Accuracy 7.6 RTD Modeling The modeling of the advection-dispersion scalar transport equations in conjunction with the flow fields developed using the FLUENT model allowed for the determination of the concentration profiles within the basin and for the generation of RTD profiles. The modeling was performed without consideration of the aerator effects and the inclusion of dispersion number values was applied globally to the solution space, not considering turbulent transportation of the conservative tracer. It was found from visual inspection of the plots generated by the solute transport model that it seemed to be modeling with accuracy and the results seemed to match the results observed when the tracer studies themselves were performed on the basin. Of special interest was the location of the dead zones and the accumulation of dye at the top-left and bottom-right corners as well as in the center of the tank. These observed patterns were mimicked by the solute transport model. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 105 Figure 7-5 below shows four snapshot times for a model with a 4 L/s flow-through rate and a 100 mass unit addition of scalar across the inlet at time t = 0. Dx=Dy=0.0005 m2/s, 60x48, Q=4 L/s, DT=0.05 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 0.00 0.50 1.00 1.50 2.00 2.50 Figure 7-5: Sample Solute Transport Time Series - Pulse Tracer Input 7.6.1 Continuity Errors It was discovered after examining the time series of the mass of the scalar in the basin and the corresponding outlet concentration-based RTD curves, that some mass was being lost during the runs. It was found that the mass that was added to the basin quickly decreased long before the tracer reached the outlet of the basin. Figure 7-6 below illustrates the sum total of mass in the solution space and the concentration at the outlet time series for the same run illustrated previously in Figure 7-5. It can be seen that the mass decreased from the initial value of 100 to Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 106 about 63 after the first few iterations and this loss did not correspond to any concentration at the outlet. It is believed this initial dramatic drop was due to diffusion of mass out of the basin through the inlet. However, it was seen that there was a slow decrease in mass in the basin until about 250 seconds, when the scalar cloud actually reached the outlet and started to exit the system normally. This drop was believed to be caused by lack of continuity in the flow fields employed in calculation. Pulse Tracer Input at t=0 (C-Curve) 0=4 L/s, DT=0.05 sec, Dx=Dy=0.0005 m2/s, 60x48 grid size 100 — i Time (seconds) Figure 7-6: Mass and Outlet Concentration Time Series - Pulse Tracer Input In an attempt to assess these issues, two steps were taken. As a first step, a continuous input of tracer was assumed instead of a pulse input of tracer while modeling, thus reducing any negative flux of the tracer out of the basin at the inlet. The second step involved an analysis of the FLUENT flow field being employed to determine whether continuity was satisfied over the whole space. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 107 The switch from a pulse injection to a continuous feed of tracer was accomplished by setting the concentration at the inlet boundary to a specified value (100 concentration units). This allowed both advection and dispersion of the tracer into the basin, and conversely limited the amount of tracer dispersing out the inlet (because the tracer concentration at the inlet would be higher than that just within the basin). Figure 7-7 and Figure 7-8 below illustrate the snapshots of the concentration profiles and the mass and concentration time series respectively, for the continuous tracer input configuration. It can be seen by examining the outlet concentration time series in Figure 7-8 that mass was not conserved. With an inlet concentration of 100 units as specified earlier, one would expect the final outlet concentration to be 100 also. This was not the case and clearly mass was being lost. Dx=Dy=0.0005 m2/s, 60x48, 0=4 L/s, DT=0.05 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 Figure 7-7: Sample Solute Transport Time Series - Continuous Tracer Input Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 108 c to co <u 3 O o c <u o o o cn co i > < 3.0E+5 2.0E+5 —\ 1.0E+5 =§ 100.00 50.00 H 0.00 Continuous Tracer Input at t=0 (F-Curve) 0=4 Us, DT=0.05 sec, Dx=Dy=0.0005 m2/s, 60x48 grid size 0 4,000 8,000 12,000 16,000 20,000 Time (seconds) Figure 7-8: Mass and Outlet Concentration Time Series - Continuous Tracer Input For a further explanation of the mass errors associated with the RTD data, the flow field itself was looked to for errors in flow continuity. The basin flow field was analyzed using the same routine used to analyze the still basin profiles for continuity in Section 4.4. This routine was modified slightly in its output to produce a surface map of the magnitude of the continuity error, normalized to the magnitude of the local velocity gradient. The equation for this is shown below: I.J F • • + F • • * x:i,j * y.ij F . . + F • . i x-i,j y-ij' (7-39) where, Fx.iJ, Fy;iJ = the normalized continuity error magnitude at cell i,j, and = the velocity fluxes in the x and y directions at cell i j . Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 109 Equation 7-39 provides a measure of the error in continuity over the solution space normalized to a local flux magnitude. If continuity had been satisfied, then the plot would produce a surface of magnitude zero over the entire flow space. The worst local continuity error can be unity. The plot of Equation 7-39 over the solution space used in the RTD calculations is shown in Figure 7-9 below. It can be seen that the majority of the space exhibited good continuity satisfaction, but several areas violated the continuity criteria quite seriously, especially near the boundaries. It was believed that this continuity error was associated with the lost mass exhibited in the RTD routines. These errors were likely caused by the transformation routine that converted the non-uniform mesh to a uniform mesh for the tracer transport modeling. 150x120 - Turb, Para - Bent 4: Mapped to 60x48 Uniform Figure 7-9: Normalized Continuity Error Map 7.6.2 Matching Laboratory RTD Profiles Although issues have been presented regarding the evident loss of mass either by diffusion through the inlet or due to continuity errors in the flow fields, the RTD profiles generated by CFD can still be of value by considering the shape of the profiles in comparison to the laboratory tracer data discussed in Section 5. It was found that when the laboratory tracer studies and the generated RTD profiles were normalized for total mass for the first two hydraulic retention Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 110 times, the profiles matched reasonably well. Figure 7-10 shows the excellent match in RTD shape between the laboratory data and the model results for 4 L/sec and no aeration. The model even captured the recirculation period between 0.5 and 1.0 retention times, something that an axial dispersion model is incapable of doing. Additionally, the timing of the tracer appearance at the outlet, the time to peak and the time to recirculation are reproduced with satisfactory accuracy. Model Data Normalized - 4 L/sec o O O 8.00 4.00 0.00 6 » » o 0.00 0.50 - 0 — Laboratory Tracer Study (June 24, 99) CFD Output (D=0.0005 m2/s) 1.00 t/HRT 1.50 2.00 Figure 7-10: Comparison of Model and Laboratory C-Curves 7.7 Conclusions A two-dimensional conservative transport model was developed using finite-volume CFD techniques. The model performed well in matching the advection and dispersion model validations, and was accurate in both space and time. The program did exhibit some problems with mass conservation that was likely attributable to two effects: (1) the loss of flow continuity in the velocity grid as the non-uniform mesh was translated to a uniform mesh, and (2) the dispersion of material out of the basin system via the inlet. The numerical methods were likely not contributing to the error due to the good degree of mass conservation exhibited during the validation of the model. Notwithstanding the mass losses, the model still produced reasonable Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 111 results and the mass-corrected RTD profile matched the tracer profiles from Section 5 very closely. Of special interest was the ability of the model to resolve the recirculation zone influence in the RTD, something an axial dispersion model would be incapable of doing. 7.8 Recommendations The scalar transport aspect of the study is very much open to improvement. Clearly some deficiencies exist in the model where mass conservation is concerned. It seems that problems lie both in the unwarranted dispersion of the tracer out of the inlet and in the fact that the flow fields translated from the FLUENT data files exhibited continuity violations. It is likely that both of these issues contributed to the errors discussed above, but the mass transport needs to be examined more precisely at the inlets and at the areas of continuity violation before work can be done on the program. It may be possible to actually code the fluxes at the inlet to specified values but this would require a reworking of the solution routine. The introduction of the influence of aerators into the CFD modeling was not considered. This would require the mapping of the dispersion numbers calculated in Section 4.5 onto the computational domain, however more confidence in these dispersion profiles is required. The modeling did not take into account the turbulent Schmidt number for turbulent diffusion as discussed in Section 6.1. The model would have to arrange the import of the turbulent data generated by the F L U E N T modeling and calculate the dispersion numbers locally based on the eddy viscosity. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 112 8 Conclusions and Recommendations 8.1 Conclusions It was the objective of this study to examine the potential of modeling ASBs with computational fluid dynamics, using data generated with a laboratory scale apparatus. This analysis was primarily investigative, examining a variety of analytical and modeling approaches. The project was successful in that the apparatus was constructed and tested in a variety of ways, including still basin analysis for aerator impacts, tracer studies on a modified flow-through basin with a variety of flow rates and aeration levels, and reasonably successful fluid flow modeling and species transport. The construction of the aerator and the basin was completed with a good degree of care to ensure a range of flow rates and mixing that ensured reasonable similarity with actual ASBs. True similarity was difficult to accomplish due to the splashing and other scale issues, but the range of flow rates of the scale apparatus seemed to produce a good variety of mixing conditions. The trends presented regarding turbulence profiles in the still basin analysis were inconclusive and the difficulty in resolving flow continuity with the vertical plane sampling technique presented analytical difficulties. However, the analysis did provide insight into the extent of the complexity of the flow imparted by aeration. Additionally, it was found that the method employed to approximate the dispersion number profiles must be refined and tested, so that some mixing characteristics might be used in a CFD mixing model. The FLUENT flow modeling was a success in that it predicted the velocity measurements in the basin with a good degree of accuracy. The variety of techniques employed to refine the solution were reasonably successful, with a turbulent model with a parabolic inlet profile providing the most accurate results. The modeling might have benefited from the inclusion of more realistic governing equations, namely equations that accounted for the shear stress from the basin floor. The manufactured CFD code for species transport was also fairly successful. The model was developed successfully with an accurate replication of the analytical advection-dispersion cases. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 113 However, when combined with the modified flow results acquired from FLUENT, the results suffered from mass conservation losses. It is suspected that the mass losses are likely due to the continuity errors associated with the flow field and/or a steady-state dispersion out of the inlet. These problems could be addressed by converting the flow fields to a structured grid with a more complex routine and by coding in a strict flux restriction at the model inlet. The reproduction of the shape of the RTD from the tracer studies without aeration was quite good, although a mass correction was required to make such a comparison. 8.2 Recommendations The following list outlines some of the recommendations for further study: 1. It is recommended that more data be collected in the still basin study. With enough velocity profiles it might be possible to elucidate trends in turbulent time scale or dispersion numbers as related to the flow through the aerator and the distance of sampling point from the aerator. 2. The resolution flow continuity in the still basin data is also an issue of great importance. It could be worth investigating whether continuity can be resolved at all using the A D V within the basin, with a series of time-averaged point data. A detailed examination of a small area within the basin in three dimensions (rather than two) may yield more fruitful results. 3. The still basin aeration analysis might also benefit from examining the flow fields introduced by an aerator close to a wall boundary. If it could be found that the profiles of turbulent diffusion could be mapped to a 2D CFD model then knowledge of the impacts of boundaries would be valuable. 4. The tracer studies could be improved upon primarily through further data collection. However, the mixing impacts of other factors such as a change in water depth, the introduction of baffles or the placement of the aerator would also be of interest. 5. The FLUENT modeling could be improved upon by attempting to incorporate the shear stress at the bottom of the basin either through the introduction of the stress to the two-dimensional equations or through three-dimensional flow modeling. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 114 6. The CFD scalar transport model could be improved upon by examining the lack of mass conservation associated with the model and determining which phenomena are the causes. The model could also be improved by including turbulent transport into the dispersion model and including a technique to include space-variable dispersion into the model that could approximate the influence of aerators. Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 115 9 References A B B O T T , M I C H A E L B. ; B A S C O , D. R.; (1989), Computational Fluid Dynamics, Longman Group A G U N W A M B A , J. C. ; (1992) Field Pond Performance and Design Evaluation Using Physical Models, Wat. Research, 26 (10) pp.1403-1407 A G U N W A M B A , J. C. ; E G B U N I W E , N . ; A D E M I L U Y I , J. O.; (1992) Prediction of the Dispersion Number in Waste Stabilization Ponds, Wat. Research, 26 pp.85-89 A Q U A - A E R O B I C S Y S T E M S , I N C . Bulletin No. 950F 2/98 (1998) A R C E I V A L A , S O L I J A L ; (1981) Hydraulic Modeling for Waste Stabilization Ponds, J. Env. Eng. Div. ASCE, 109 pp.265-268 A J R M E N A N T E , P I E R O M . ; L U O , C H A N G G E N ; C H O U , C H U N - C H I A O ; F O R T , I V A N ; M E D E K , J A R O S L A V ; (1997) Velocity profiles in a closed, unbaffled vessel: Comparison between experimental L D V data and numerical CFD predictions, Chem. Eng. Sci. 52(20) pp. 3483-3492 C E L I K , L. ; R O D I , W O L F G A N G ; S T A M O U , A . I.; (1985) Prediction of Hydrodynamic Characteristics of Rectangular Settling Tanks, Turb Meas. And Flow Modeling, ASCE, pp.641-651 C H O L E T T E , A . ; C L O U T I E R , L . ; (1959) Mixing Efficiency Determinations for Continuous Flow Systems, Can. J. Chem. Eng., 37 pp. 105-D A N C K W E R T S , P. V . ; (1953) Continuous Flow Systems: Distribution of Residence Times, Chem. Eng. Sci., 2 (1) pp. 1-13 D O R E G O , N . C ; L E D U C , R O L A N D ; (1996) Characterization of Hydraulic Flow Patterns in Facultative Aerated Lagoons, Wat. Sci. Tech., 34 (11) pp.99-106 E L D E R , J. S.; (1959) The dispersion of marked fluid in turbulent shear flow, Fluid Mech., 5 pp.544-560 F E R R A R A , R A Y M O N D A . ; H A R L E M A N , D O N A L D R. F.; (1980) Dynamic Nutrient Cycle Model for Waste Stabilization Ponds, J. Env. Eng. Div. ASCE, pp.37-54 F E R R A R A , R A Y M O N D A . ; H A R L E M A N , D O N A L D R. F.; (1981) Hydraulic Modeling for Waste Stabilization Ponds, J. Env. Eng. Div. ASCE, 107 pp.817-830 F I S H E R , H U G O B. ; L I S T , E. J O H N ; K O H , R O B E R T C. Y . ; I M B E R G E R , J O R G ; B R O O K S , N O R M A N H. ; (1979), Mixing in Inland and Coastal Waters, Academic Press F R E U N D , J O H N E.; (1962), Mathematical Statistics, Prentice-Hall, Inc., Englewood Cliffs H A D D A D , A . H . ; W O L F , D.; (1967) Residence Time Distribution Function for Mulit-Stage Systems with Backmixing, Can. J. Civ. Eng., 45 pp. 100-104 H E , G U A N G B I N ; G U O , Y A H N U ; H S U , A N D R E W T.; (1999) The Effect of Schmidt Number on Turbulent Scalar Mixing in a Jet-in-Crossflow, Int. J. Heat Mass Transfer., 42 pp.3727-3738 K R A U S , N I C H O L A S C ; L O H R M A N N , A T L E ; C A B R E R A , R A M O N ; (1994) New Acoustic Meter for Measuring 3D Laboratory Flows, J. Hyd. Eng., 120 (3) pp.406-413 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 116 L E V E N S P I E L , O C T A V E ; (1962), Chemical Reaction Engineering, John Wiley & Sons L E V E N S P I E L , O C T A V E ; L A I , B. W.; C H A T L Y N N E , C . Y . ; (1970) Tracer Curves and the Residence Time Distribution, Chem. Eng. Sci., 25 pp.1611-1613 L E V E N S P I E L , O C T A V E ; S M I T H , W. K . ; (1957) Notes on the Diffusion-Type Model for the Longitudinal Mixing of Fluids in Flow, Chem. Eng. Sci., 6 pp.227-233 L E V E N S P I E L , O C T A V E ; T U R N E R , J. C. R.; (1970) The Interpretation of Residence-Time Experiments, Chem. Eng. Sci., 25 pp.1605-1609 L I U , H . ; (1977) Predicting Dispersion Coefficient of Streams", J. Env. Eng. Div. ASCE, 103 (EE1) pp. 59-69 L O H R M A N N , A T L E ; C A B R E R A , R A M O N ; G E L F E N B A U M , G U Y ; H A I N E S , J O H N ; (1995) Direct Measurements of Reynolds Stress with an Acoustic Doppler Velocimeter, Proc. IEEE Conf. Current Meas., pp.205-210 L O H R M A N N , A T L E ; C A B R E R A , R A M O N ; K R A U S , N I C H O L A S C ; (1994) Acoustic-Doppler Velocimeter (ADV) for Laboratory Use, Hyd. Meas. Exper., pp.351-365 L Y N , D. A . ; R O D I , W O L F G A N G ; (1990) Turbulence Measurements in Model Settling Tank, J. Hyd. Eng., 116 (1) pp.3-21 M A N G E L S O N , K E N N E T H A . ; W A T T E R S , G A R Y Z . ; (1972) Treatment Efficiency of Waste Stabilization Ponds, J. San. Eng. Div., ASCE, 2 pp.407-425 M A R E C O S D O M O N T E , M . H . F.; M A R A , D U N C A N D.; (1987) The Hydraulic Performance of Waste Stabilization Ponds in Portugal, Wat. Sci. Tech., 19 (12) pp.219-227 M U R P H Y , K E I T H L . ; W I L S O N , A L F R E D W A R R E N ; (1974) Characterization of Mixing in Aerated Lagoons, J. Env. Eng. Div. ASCE, 100 (10) pp.1105-1117 N A M E C H E , T H ; V A S E L , J. L. ; (1998) Hydrodynamic Studies and Modelization for Aerated Lagoons and Waste Stabilization Ponds, Wat. Research, 12 (10) pp.3039-3045 NCASI; (1971) A study of mixing characteristics of Aerated Stabilization Basins, NCASI Tech. Bull., No. 245 NCASI; (1983) A review of procedures for conducting conservative tracer studies in the hydraulic characterization of effluent treatment basins, NCASI Tech. Bull., No. 408 P E A R S O N , H . W.; M A R A , D U N C A N D.; A R R I D G E , H . A. ; (1995) Influence of pond geometry and configuration on facultative and maturation waste stabilisation pond performance and efficiency, Wat. Sci. Tech., 31 (12) pp.129-139 P O L P R A S E R T , C H O N G R A K ; B H A T T A R A I , K I R A N K . ; (1985) Dispersion Model for Waste Stabilization Ponds, J. Env. Eng. Div. ASCE, 111 pp.45-59 P O Z R I K I D I S , C ; (1997), Introduction to Theoretical and Computational Fluid Dynamics, Oxford University Press P R E U L , H . C ; W A G N E R , R. A . ; (1987) Waste Stabilization Pond Prediction Model, Wat. Sci. Tech., 19(12) pp.205-211 R O D I , W O L F G A N G ; (1980), Turbulence Models and their Application in Hydraulics, Int. Assoc. Hydraulic Research Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 117 S A C K E L L A R E S , R O B E R T W.; B A R K L E Y , W I L L I A M A . ; (1985) Selection and Verification of Hydraulic Flow Models for Lagoons, Forest Products, 81 pp.117-130 S A C K E L L A R E S , R O B E R T W.; B A R K L E Y , W I L L I A M A . ; H I L L , R. D.; (1987) Development of a Dynamic Aerated Lagoon Model, Wat. Env. Res., 57 (10) pp.877-883 S C O T T , J O N ; C E L I K K O L , B A R B A R O S ; B E A L , C. S C O T T ; S A N T A M A R I A , J O E ; (1998a) Improving the Performance of an Aeration Basin by Implementing the Results of a Computer Model of Both Sections of the Aeration Basin, TAPPIEnv. Conf. Proc, pp.1375-1387 S C O T T , J O N ; C E L I K K O L , B A R B A R O S ; S I E N E R , C A R L ; B U T A U D , B L A I N E ; S A N T A M A R I A , J O E ; (1998b) Restoration of Capacity to an Aeration Stabilization Basin Using Innovative Techniques, TAPPIEnv. Conf. Proc, pp.1388-1396 S P I E G E L , M U R R Y R.; (1997), Mathematical Handbook, McGraw-Hill T A , C. T.; B R I G N A L , W. J.; (1998) Application of Computational Fluid Dynamics Technique to Storage Reservoir Studies, Wat. Sci. Tech., 37 (2) pp.219-226 T H I R U M U R T H I , D H A N D A P A N I ; (1969a) A Break-Through in the tracer Studies of Sedimentation Tanks, J. AWWA, 41 (11) pp.405-418 T H I R U M U R T H I , D H A N D A P A N I ; (1969b) Design Principles of Waste Stabilization Ponds, J. San. Eng. Div., ASCE, 95 pp.311-330 T H I R U M U R T H I , D H A N D A P A N I ; (1974) Design Criteria for Waste Stabilization Ponds, J. Wat. Poll. ContorlFed., 46 (9) pp.2094-2104 T O R R E S , J. J.; S O L E R , A . ; S A E Z , J.; O R T U N O , J. F.; (1977) Hydraulic performance of a deep wastewater stabilization pond, Wat. Research, 31 (4) pp.679-688 V A N L E E R , B R A M ; (1974) Towards the Ultimate Conservative Difference Scheme. II. Monotonicity and Conservation Combined in a Second-Order Scheme, J. Compu. Phys., 14, pp.361-370 V O U L G A R I S , G E O R G E ; T R O W B R I D G E , J. H. ; (1998) Evaluation of the Acoustic Doppler Velocimeter (ADV) for Turbulence Measurements, J. Atmos. Ocean. Tech., 15(1) pp.272-289 W E H N E R , J. F.; W I L H E L M , R. H. ; (1956) Boundary Conditions of Flow Reactor, Chem. Eng. Sci., 6 pp.89-93 W O O D , M . G.; G R E E N F I E L D , P. F.; H O W E S , T.; J O H N S , M . R.; K E L L E R , J.; (1995) Computational Fluid Dynamic Modelling of Wastewater Ponds to Improve Design, Wat. Sci. Tech., 31 (12) pp.111-118 W O O D , M . G.; (1997) Development of Computational Fluid Dynamic Models for the Design of Waste Stabilization Ponds, Ph.D. Thesis, University of Queensland Y E , J I A N ; M C C O R Q U O D A L E , A L E X J.; (1997) Depth-Averaged Hydrodynamic Model in Curvilinear Collocated Grid, J. Hyd. Eng., 123 (5) pp.380-388 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 118 APPENDIX A: Aerator Calibration Statistical Analysis of the Aerator Calibration Data Examining the relative influence of the P, H , and I parameters. Statistical Output provided by the Microsoft Excel 97 Data Analysis tools. SUMMARY OUTPUT Regression Statistics Multiple R 0.985 R Square 0.969 Adjusted R Square 0.968 Standard Error 146.291 Observations 60.000 ANOVA df SS MS F Significance F Regression 3.000 38043142.436 12681047.479 592.545 0.000 Residual 56.000 1198454.997 21400.982 Total 59.000 39241597.433 Coefficients Standard Error f Stat P-value Lower 95% Upper 95% Intercept -248.190 94.329 -2.631 0.011 -437.153 -59.228 P 539.735 13.795 39.125 0.000 512.100 567.369 H -120.856 12.358 -9.780 0.000 -145.612 -96.101 I -9.913 5.464 -1.814 0.075 -20.859 1.032 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 119 APPENDIX B: Stream Function Calculation Subroutine for the generation of stream function plots based on two arrays of velocity data stored in SURFER grid file format. Sub CalcStreamArray(sU() As S i n g l e , sV() As S i n g l e , sStreamO As S i n g l e , sDX As S i n g l e , sDY As Si n g l e , iNx As Integer, iNy As Integer) 1 Where 1 sU, sV - V e l o c i t y arrays ' sStream -Stream flow output array 1 sDX, sDY - D i s c r e t i z a t i o n increments ' iNx,iNy -Array S i z e s Dim i As Integer Dim j As Integer For j = 2 To iNy sStream(l, j ) = sStream(l, j - 1) + s U ( l , j ) * sDY Next j For j = 1 To iNy For i = 2 To iNx sStream(i, j ) = sStream(i - 1, j) - s V ( i , j) * sDX Next i Next j End Sub Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 120 APPENDIX C: Continuity Macros V B A Macros - Cartesian Continuity Macro and Radial Continuity Macro Used to generate plots of continuity in two coordinate systems Sub CartesianContinuityMacro() 1 Macrol Macro 1 Macro recorded 04/05/99 by Wayne Jenkinson 1 This macro w i l l generate X, Z, du/dx, dw/dz column data 1 i n t o " C o n t i n u i t y " worksheet 1 One must have the data s t o r e d i n X,Z,U,V,W rows and run 1 the macro i n that sheet. The rows must be sorted by z, then x 1 C o n t i n u i t y du/dx + dv/dy + dw/dz = 0 ' assuming: dv/dy = 0 Dim rwlndex As Integer Dim c o l l n d e x As Integer Dim bolWithinDataSet As Boolean Dim sglFlow() As S i n g l e Dim sglDX() As S i n g l e Dim sglDYO As S i n g l e Dim sglFUp As S i n g l e Dim sglFDn As S i n g l e Dim sglDUp As S i n g l e Dim sglDDn As S i n g l e Dim strSizeMarker As S t r i n g Dim intNumX As Integer Dim intNumY As Integer Dim i As Integer Dim j As Integer •X,Z,U,V,W ReDim sglFlow(5, 150, 150) ReDim sglDX(150, 150) ReDim sglDY(150, 150) 'Run through the data t o get s i z e s intNumX = 1 ' f i n d r e p e a t i n g s e c t i o n l e n g t h i n z strSi z e M a r k e r = A c t i v e S h e e t . C e l l s ( 2 , 2).Text Do intNumX = intNumX + 1 Loop While ActiveSheet.Cells(intNumX + 1, 2).Text = strSizeMarker intNumX = intNumX - 1 1 F i n d t o t a l l e n g t h of record intNumY = 1 StrSizeMarker = "" Do intNumY = intNumY + 1 Loop U n t i l ActiveSheet.Cells(intNumY + 1, D.Text = strSizeMarker intNumY = intNumY - 1 intNumY = intNumY / intNumX ' C o l l e c t Data bolWithinDataSet = True rwlndex = 2 For j = 1 To intNumY Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 121 For i = 1 To intNumX For colIndex = 1 To 5 sgl F l o w ( c o l l n d e x , i , j ) = Acti v e S h e e t . C e l l s ( r w l n d e x , colIndex).Value Next c o l l n d e x rwlndex = rwlndex + 1 Next i Next j ' W i l l c a l c u l a t e the f l u x e s by f i n d i n g the average values between p o i n t s 1 and then i n t e r p o l a t i n g between those p o i n t s For j = 1 To intNumY For i = 1 To intNumX 1 c a l c dY I f i = 1 Then sglDUp = ( s g l F l o w ( l , i + 1, j) - s g l F l o w d , i , j ) ) s g l D X ( i , j) = (sglFlow(3, i + 1, j) - sglFlow(3, i , j ) ) / (sglDUp) E l s e l f i = intNumX Then sglDDn = ( s g l F l o w ( l , i , j) - s g l F l o w d , i - 1, j ) ) s g l D X ( i , j) = ( s g l F l o w d , i , j) - sglFlow(3, i - 1, j ) ) / (sglDDn) E l s e sglFUp = ( s g l F l o w d , i + 1, j ) + sglFlow(3, i , j ) ) / 2 sglFDn = (sglFlow(3, i - 1, j ) + sglFlow(3, i , j ) ) / 2 sglDUp = ( s g l F l o w d , i + 1, j ) - s g l F l o w d , i , j ) ) sglDDn = ( s g l F l o w d , i , j ) - s g l F l o w d , i - 1, j ) ) s g l D X f i , j) = (sglFUp - sglFDn) / (sglDUp + sglDDn) / 2 End I f 'calc dx I f j = 1 Then sglDUp = ( s g l F l o w d , i , j +1) - sglFlow(2, i , j ) ) s g l D Y ( i , j) = (sglFlow(5, i , j + 1) - sglFlow(5, i , j ) ) / (sglDUp) E l s e l f j = intNumY Then sglDDn = (sglFlow(2, i , j) - sglFlow(2, i , j - 1 ) ) s g l D Y ( i , j) = (sglFlow(5, i , j) - sglFlow(5, i , j - 1)) / (sglDDn) E l s e sglFUp = (sglFlow(5, i , j + 1) + sglFlow(5, i , j ) ) / sglFDn = (sglFlow(5, i , j) + sglFlow(5, i , j - 1 ) ) / sglDUp = ( s g l F l o w d , i , j +1) - sglFlow(2, i , j ) ) sglDDn = ( s g l F l o w d , i , j) - s g l F l o w d , i , j - 1)) s g l D Y ( i , j) (sglFUp - sglFDn) / (sglDUp + sglDDn) / 2 End I f Next i Next j Debug.Print Wo r k s h e e t s ( " C o n t i n u i t y " ) . A c t i v a t e A c t i v e S h e e t . C e l l s ( 1 , 1).Value = "X" A c t i v e S h e e t . C e l l s ( 1 , 2).Value = "Z" A c t i v e S h e e t . C e l l s ( 1 , 3).Value = "U" A c t i v e S h e e t . C e l l s ( 1 , 4).Value = "W" A c t i v e S h e e t . C e l l s ( 1 , 5).Value = "dU/dX" A c t i v e S h e e t . C e l l s ( 1 , 6).Value = "dw/dz" Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 122 rwlndex = 2 For j = 1 To intNumY For i = 1 To intNumX Ac t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , rwlndex = rwlndex + 1 Next i Next j .Value = .Value = .Value = .Value = .Value = .Value = sglFlow( sglFlow( sglFlowf sglFlow( s g l D X ( i , s g l D Y ( i , 1, i . 2, i , 3, i , 5, i , j) j) j ) j) j) j) End Sub Sub RadialContinuityMacro() 1 This macro i s as CartesianContinuityMacro except the r a d i a l coordinate ' c o n t i n u t y dtaa w i l l be generated. 1 This macro w i l l generate X, Z, du/dx, dv/dy column data ' i n t o " C o n t i n u i t y " worksheet 1 One must have the data s t o r e d i n X,Z,U,V,W rows and run 1 the macro i n that sheet. The rows must be sorted by z, then x 'coordinates: r , t h e t a , z 'Continuity: 1/r d ( r v ( s u b ) r ) / d r + 1/r d(v(sub)theta)/d t h e t a + d v(sub)z / dz = 0 'Assuming: 1/r d(v ( s u b ) t h e t a ) / d t h e t a = 0 Dim rwlndex As Integer Dim c o l l n d e x As Integer Dim bolWithinDataSet As Boolean Dim sglFlow() As S i n g l e 'Data from the spreadsheet, X,Z,U,V,W Dim sglDX() As S i n g l e Dim sglDYO As S i n g l e Dim sglFUp As S i n g l e Dim sglFDn As S i n g l e Dim sglDUp As S i n g l e Dim sglDDn As S i n g l e Dim str S i z e M a r k e r As S t r i n g Dim intNumX As Integer Dim intNumY As Integer Dim i As Integer Dim j As Integer 'X, Z,U, V,W ' (Let x = r) ReDim sglFlow(5, 150, 150) ReDim sglDX(150, 150) ReDim sglDY(150, 150) 'Run through the data t o get s i z e s intNumX = 1 StrSizeMarker = A c t i v e S h e e t . C e l l s ( 2 , 2).Text Do intNumX = intNumX + 1 Loop While ActiveSheet.Cells(intNumX + 1 , 2).Text = strSizeMarker intNumX = intNumX - 1 intNumY = 1 strSizeMarker = "" Do intNumY = intNumY + 1 Loop U n t i l ActiveSheet.Cells(intNumY + 1, 1).Text = strSizeMarker intNumY = intNumY - 1 intNumY = intNumY / intNumX Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 123 ' C o l l e c t Data bolWithinDataSet = True rwlndex = 2 For j = 1 To intNumY For i = 1 To intNumX For colIndex = 1 To 5 sgl F l o w ( c o l l n d e x , i , j) = ActiveSheet.Cells(rwlndex, c o l l n d e x ) . V a l u e Next c o l l n d e x rwlndex = rwlndex + 1 Next i Next j ' W i l l c a l c u l a t e the f l u x e s by f i n d i n g the average values between p o i n t s ' and then i n t e r p o l a t i n g between those p o i n t s For j = 1 To intNumY For i = 1 To intNumX 'cal c X I f i = 1 Then sglDUp = (sglFlow(1, i + 1, j) - sglFlow(1, i , j ) ) s g l D X ( i , j) = (sglFlow(3, i + 1, j) - sglFlow(3, i , j ) ) / (sglDUp) E l s e l f i = intNumX Then sglDDn = (sglFlow(1, i , j) - sglFlow(1, i - 1, j ) ) s g l D X ( i , j ) = (sglFlow(3, i , j) - sglFlow(3, i - 1, j ) ) / (sglDDn) E l s e sglFUp = (sglFlow(3, i + 1, j) + sglFlow(3, i , j ) ) / 2 sglFDn = (sglFlow(3, i - 1, j) + sglFlow(3, i , j ) ) / 2 sglDUp = ( s g l F l o w d , i + 1, j ) - s g l F l o w d , i , j ) ) sglDDn = ( s g l F l o w d , i , j ) - s g l F l o w d , i - 1, j ) ) s g l D X ( i , j) = (sglFUp - sglFDn) / (sglDUp + sglDDn) / 2 End I f ' Here i s the change f o r r a d i a l : ' 1/r d ( r u ) / d r = du/dr + u/r s g l D X ( i , j) = s g l D X ( i , j ) + sglFlow(3, i , j) / s g l F l o w d , i , j ) 'calc Y I f j = 1 Then sglDUp = (sglFlow(2, i , j +1) - sglFlow(2, i , j ) ) s g l D Y f i , j) = (sglFlow(5, i , j + 1) - sglFlow(5, i , j ) ) / (sglDUp) E l s e l f j = intNumY Then sglDDn = (sglFlow(2, i , j ) - sglFlow(2, i , j - 1)) s g l D Y ( i , j) = (sglFlow(5, i , j ) - sglFlow(5, i , j - 1)) / (sglDDn) E l s e sglFUp = (sglFlow(5, i , j + 1) + sglFlow(5, i , j ) ) / 2 sglFDn = (sglFlow(5, i , j) + sglFlow(5, i , j - 1 ) ) / 2 sglDUp = (sglFlow(2, i , j +1) - sglFlow(2, i , j ) ) sglDDn = (sglFlow(2, i , j) - sglFlow(2, i , j - 1 ) ) s g l D Y ( i , j ) = (sglFUp - sglFDn) / (sglDUp + sglDDn) / 2 End I f Next i Next j Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 124 Debug.Print Worksheets("RContinuity").Activate A c t i v e S h e e t . C e l l s ( 1 , 1).Value = "X" A c t i v e S h e e t . C e l l s ( 1 , 2).Value = "Z" A c t i v e S h e e t . C e l l s ( 1 , 3).Value = "U" A c t i v e S h e e t . C e l l s ( 1 , 4).Value = "W" A c t i v e S h e e t . C e l l s ( 1 , 5).Value = "dU/dR + A c t i v e S h e e t . C e l l s ( 1 , 6).Value = "dW/dZ" rwlndex = 2 U/R" For ] = 1 T o intNumY For i = 1 To intNumX Ac t i v e S h e e t . C e l l s ( r w l n d e x , 1).Value = s g l F l o w ( l , i , j) Act i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , A c t i v e S h e e t . C e l l s ( r w l n d e x , rwlndex = rwlndex + 1 Next i Next j 2) .Value = sglFlow(2, i , j) 3) .Value = sglFlow(3, i , j) 4) .Value = sglFlow(5, i , j) 5) .Value = s g l D X ( i , j) 6) .Value = s g l D Y ( i , j) End Sub Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 125 APPENDIX D: Autocorrelation Code Autocorrelation Code - executed when the A D V machine files are processed for velocity data. P u b l i c Sub Bat c h V e c t C a l c A u t o c o r r e l a t i o n O 'Subroutine t o c a l c u l a t e the a u t o c o r r e l a t i o n f u n c t i o n and time s c a l e of ' v e l o c i t y time s e r i e s 'Referes t o common block of v a r i a b l e s d e f i n e d as f o l l o w s 'intNumSamples -number of time samples 'sglDT -time increment between samples 'sglXRCorO -array h o l d i n g time and a u t o c o r r e l a t i o n f u n c t i o n data (as f o r y and z) 'sglXTimeScale -value f o r x - d i r e c t i o n time s c a l e (as f o r y and z) 'sglXv e l O - v e l o c i t y time s e r i e s (as f o r y and z) 'sglXvelBar -average v e l o c i t y over time s e r i e s (as f o r y and z) Dim i As Long Dim j As Long Dim sglSum As S i n g l e Dim sglDS As S i n g l e Dim bolNegativeFound As Boolean Dim intMaxIndex As Long ReDim sglXRcor(2, intNumSamples) ReDim sglYRcor(2, intNumSamples) ReDim sglZRcor(2, intNumSamples) 'Solve f o r A u t o c o r r e l a t i o n Function ' X bolNegativeFound = Fals e IngXRcorDim = intNumSamples For i = 1 To intNumSamples - 1 sglSum = 0 For j = 1 To intNumSamples - i - 1 sglSum = sglSum + ( s g l X v e l ( j ) - sglXvelBar) * ( s g l X v e l ( j + i ) - sglXvelBar) Next j s g l X R c o r d , i ) = i * sglDT 'Place S i n the f i r s t vect sglXRcor(2, i ) = sglSum / (intNumSamples - i ) / sglXX DoEvents I f s g l X R c o r d , i ) < 0 And bolNegativeFound = False Then IngXRcorDim = i bolNegativeFound = True End I f ' K i l l the loop once we've gone 10% f a r t h e r than the f i r s t negative ' value I f i > Int(IngXRcorDim * 1.1) Then E x i t For Next i 'Solve f o r the X Time Scale For i = 1 To IngXRcorDim sglXTimeScale = sglXTimeScale + sglDT * ( s g l X R c o r d , i ) + s g l X R c o r d , i + 1) ) / 2 Next i DoEvents Solve f o r A u t o c o r r e l a t i o n Function y Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 126 bolNegativeFound = False IngYRcorDim = intNumSamples For i = 1 To intNumSamples - 1 sglSum = 0 For j = 1 To intNumSamples - i - 1 sglSum = sglSum + ( s g l Y v e l ( j ) - sglyvelbar) * ( s g l Y v e l ( j + i ) - sglyve l b a r ) Next j s g l Y R c o r ( l , i ) = i * sglDT 'Place S i n the f i r s t vect sglYRcor(2, i ) = sglSum / (intNumSamples - i ) / sglYY DoEvents I f sglYRcor(2, i ) < 0 And bolNegativeFound = False Then IngYRcorDim = i bolNegativeFound = True End I f ' K i l l the loop once we've gone 10% f a r t h e r than the f i r s t negative ' value I f i > Int(IngYRcorDim * 1.1) Then E x i t For Next i 'Solve f o r the Y Time Scale For i = 1 To IngYRcorDim sglYTimeScale = sglYTimeScale + sglDT * (sglYRcor(2, i ) + sglYRcor(2, i + 1)) / 2 Next i DoEvent s 'Solve f o r A u t o c o r r e l a t i o n Function ' z bolNegativeFound = Fals e IngZRcorDim = intNumSamples For i = 1 To intNumSamples - 1 sglSum = 0 For j = 1 To intNumSamples - i - 1 sglSum = sglSum + ( s g l Z v e l ( j ) - sglZvelBar) * ( s g l Z v e l ( j + i ) - sglZvelBar) Next j s g l Z R c o r d , i ) = i * sglDT 'Place S i n the f i r s t vect sglZRcor(2, i ) = sglSum / (intNumSamples - i ) / sglZZ ' I F r e e s t y l e Debugging DoEvents I f sglZRcor(2, i ) < 0 And bolNegativeFound = False Then IngZRcorDim = i bolNegativeFound = True End I f ' K i l l the loop once we've gone 10% f a r t h e r than the f i r s t negative ' value I f i > Int(IngZRcorDim * 1.1) Then E x i t For Next i DoEvents 'Solve f o r the Y Time Scale For i = 1 To IngZRcorDim sglZTimeScale = sglZTimeScale + sglDT * (sglZRcor(2, i ) + sglZRcor(2, i + 1)) / 2 Next i End Sub CM CO (D (HOVH) ooucqjosqv • CN CD CD i- co m - .. • -r- CN CN CN O CN . - . . - . . j r - ^ - c o a o o o c Q — - m o o a S c N Q f O Q i ^ ^ r c o .__ _ __ioooo'-'^ -cncOT-or— : O O O O O O O o f - ; f - ; o o o o T - c O T f - ' ^ -<5 co m m cn in in x u . O) CM <P Q) CM m i : D ifl IO N Oi O ( I CM tO O CO CO ' : i- T- CM CM CM CO -! O O O O O O t I O) o co T- in ift > cn oo cn in m m . J <D ai N co •- ro >95754 524315 548482 S68255 D124 0.! B977 O.f 9698 O.f B289 O.f t cB m u> o o o eg (0 (o o o S CO r- <D d t*- *t T-. . _ CM in 3 co CM O^COCO'tCM'^CO'^-Tf r ^ c o r ^ i ^ - ' g - i n c D C M C Q C M MTTNOClcOCDrtcO ^ CO CO ^ CO O) O) ( 0 < 0 0 S o d o o o o o I c i & o s f f l i n o r t t N n p • C D S S N I f l C O T - O O C O N f ) 0 ' C 0 0 ' N N N C 0 0 ' C O O ) * ) o o o o o o o o > CO CM t. -• o i r o w s c n s s c . . . . imOPJi-NlfiNOOH o o m o o o ' 00 CM O O O 5 S S g o r o O O O S inininNcocom^cN'-cococoneoconnro^nnionrocoNNconn MSUBCOCCOOIVtDinfS'tONTfNVNNNOSSrtCOCpOOOO . . . . -J^ Olcn^ CfiCM^ CMin^ CnulCTO^ CMCM C^M O^lCMCNCQCNCNCOO COCOCOCOCOCOCD^COCDCM^tO^COT-l^lrtCO^W |OOOOOOOOo-OOC0O^CMOCMininoa>Ui^C0CMCMC0CMC0QCMCM^ CMCM^QQOOOO u i c o c o c o c o c o i n i n i c o 3 c 2 S c 3 t o c o c o c o c o a) E E O i -s t ^ c o n r t g s u i w o ) ( D t ^ ^ O t w c 9 < D C > i N r t o i i f i ^ c o ' t o ( D o i i n f f l o i D o c n ^ ^ s c j c s i i n c o r t o i i N i n Q i o ^ i n T-NrtlBoSCOOi-NnU5COKON^CV|OfiiBCO^ONrtCfllDOJa) OOOO 0-OO»-->-*-»-*-*-T-CMCMCM • "* "* -6 d 6 d d d d o o o o o d o o d o s p j c i m c o i - ' * c o n m M ' - s t o ScoSC*55r**CMh" M c i ( l ( 0 ^ i n i / ) c D i D s c o c o o ' - N n i n T ; o N i n o i n o d d d d d o d d o d d d ^ J ^ J * - ' ^ J ' « - : CN CM CN co cd ^ in ID s M _ CO CM CO O " • — r— in CM c • N >- n T- i 8 s f 1 i l = : • < s i l l s i - ° < h -8 1? £ SB ? C\| CD S C) O) I— i" _ . C O C M T - ^ ^ O O T - 0 ' > - 0 0 0 0 0 0 0 _ J O O O O O O O O O O O O O O O O O O O S CO CO CO CO O O O O O O O O O d d d d d > r - to TJ- m -t • O O O O O o ( " d d d d d d d d d d ; o o o o o o o o o o o o c o o o o o c i c i o o c i o c i t i i o n m « i o i f l U ) ( D m i - c N c o T f i n i n o m m m g i •^ -^ -^ -^^ CM C^M C^Oj i- CM CM ( > © m Q O © © ( r in in co r^ - co cn c > r-- r-- co -a- r-- c j S co T- T- m ro cn _ o c M c n c n i i i s o T - a> T-o o o o o o o o o o o o o ^- S S CO CI N CO w r \ i 3 u } ^ ( O i n c o ( \ i c 3 _ i n C O Q N ^ N O ' * l f i ( D S ' - o j N t N r O r o r i c o ^ Q ' t ' t ^ v v u j o o o o o c i d d d d d d o d d r ^ c c i n t o ^ S c M c o c ^ T - c ? c o c M ^ r > - i > - t D r ^ c o c M c o T - r ^ ^ c o < N O C r ) 3 ) c O f O S l f i Q C O N O l C O l O N C n Q S C O Q C O Q U I O O I O ) o d d d d d d co n N co in — TJ- cn r- co CO CO ^ CO tO CM 3 S) i n 0) O CO (D CD CM IQ CO CO O CO (D r - f \ l CD r -(O K CO CO Q CO ^ ^ d ^ °5 o" O O O O O CM O O O CM E O l | | OS T 3 : O h h P P S 5 2 S 5 CO CD CO S S 2 O CO < ^ co" co ko 3 = CO CD CD 43 „ to cq 3 E C 1= C/> • «Tl o o ro ro O f 2 E 5<S3! LQ CO f^ -CO CO CO CJ) CO CD S CM CM CO T-w P ^ d 5 0 0 CD 01 CO — r- 1-(D O r-- o . . ( D O O O d d o" 1 co CM i n co <S tn rC 8) ( ) 3 N i f t c i 5 c o ^ i D o u i c 6 » - T - ' i c o c N c o s c s i o r t N ' - l n o i n i . ^ c r i T j - c M ' ^ r - o i n c M i n i O O CM 1- CM CM CM 1 i n c o c o c o c o i n u i n c o c o c o i n f o i n i f ^ c v j o j o i o i n t D C N r o r t x r c n r r o ) l ( N i c j > o ) u ^ c M u i c M ( N ( O i o c n ( o a ) c o ; C M C M C M L O ^ ^ C O C N ^ - T - O O p O ; d o d d d b d d d o ' d d d d d d d r*- f -10 cn co to cS to 3 _. o^cr>c'KcT>Scr)ScBOT s i n Q i - m N u i c o ^ c n s w n T - m s o c o o i o i c n T - c o u ) T - c o u 3 ( N O ) u i N O ) t o « r ' " ' — • - •** O ' - r o i o i B c o o ' - r t i n t o c g o ' ^ ( N C M C N C N t N C M t O C O C O C O C O C O ^ l - ' ; in to r— co ai o i n ^ c o c o c o c M c o c T J u i t D T f r ^ ^ c o m c o r ^ r ^ c o c ^ O Q O O C O i n C C ) O r ^ O ^ C O C O U l U l ^ C O ' « - » - ' » - O Q O T - - r - * - 0 0 0 0 0 0 0 0 O O O Q - O O O Q - O Q - O O O O O O O O O O O O O O O O Q O O O O O O O o d d o d d d d d d d d d d d d d d d d d d d 0 0 0 * 0 0 0 0 • f a ' t o o > OI • 1 § i i s : <3°< I UO UO Tf CM CO LO i Tt- co r- uj co uo I N T- t CO CM O) I T - CO s o ci m ' S O) CM N ; o o o r - T- ^ i d o o d ° ° ° ° ° ° ° S CO CM CO i-TT CO LO LQ I Tf CM CM Tf CO : TT m TT r*- co 3° o o o d d m CM CJ) CD Tf Tf CO 00 S CO CO s CO O CO LO C M i n c o O T T f T r i n c M c o c o r ^ c N i n c o ^ c M C J i c D C M c o o i i r ) ^ -c o o c n i n c o ^ T r c o c o c o ^ c o T f Q c o c o c D o i n c M T - T -c o c M C n O T o m i n c j ) T f c D O f ^ c M c T ) ^ r ^ C M o c c S T f c o c M o i o n c o o i T - N i n c o T t c o o i i n ^ t r o n ' - u j T r s T -o c g t t o c o c o m c O ' - ^ - o i r o s T - i n c O ' j ' C O S ' - ^ c o — o d d d ci <^ d d cS d d o d o s ^ W T - T - N j N i o c o c o i n T - T - T - c o c D r i N i o s o J c o COCOi-^COCMCOintOCO^tOCDCDCOr^Tfr^CDOtOCMCO COSCOOIOllOOlOUJOCOrflO'-SniOtOIOOOOS r ^ c j j - ^ c o i n c o c M T f r ^ O T f c o x ^ c o T - c D C M c o c M i n c o o ) T-; v CM CM CM CM CO CO CO Tf Tf Tf U") in CD CD f1-- CO CO CO CO d d c i d o G c i c i c i c i O O O O O O O O O O O O O O O O O O CO CM CM LO f- Tf Tf CM O O i- o o CD co Tf o s o g CM CO CO CO CO LO CM CO CM CD CM S l i ) CI) Q O o —^ csi Tf (J) O co in cj) CO f - Tf 9 E .E? II d) H i - p ! 1 « I ^ ^ - O H W I cl <S S o o o E x ? ? ^ u Q. l ° 8 S o caOOO i m TT to to to T- • T- t- t- t o CD tO -IS l l e * x o E a TO Q-OI ^ £ " i •e s o > s s CD CD > co 25 i ^ CO CM I -» O CO I : E - o i CD g> i 5 £ ~ CO CO I CU E o m CM CM O CM (D CO CO O Tf T-d P "2 * o o o Ct ; > O O X CJ < m o o n i M ' - m o j N c o o M D i n m ^ m s ^ S ' -o c M C M s c n o ' - s c o T - c o c o o i T r r o c n c o s n i n n i O i -to 0) in m m co m - - - - i n o m CM co CM T- o r-T- CO T-o o o o o tocooicocMiocncMLOcocMincoT-Tf^r-S i - ^ c s i n c o c i ^ - T f ^ - i n i r i i n c D i o s s C M c o T f o c o t o c f t r M i n c o - < - T f r ^ o ( o c M - -m c o c o c o o ^ T f c M O ) c O ' * T - c o t o o i n c T ) (5 P O) S C) O) If) • Tf CD r - CM Tf ( Tf O t o i- CO • I CD Lf) CM CO UO • CM CM m CO CO o • to t o 35 j d d T < • o > CO CO CJ) ( > CM if) CM ( • N CM . to in co - co cn o > d d T-" ( o t o i n c g i n i o i n o ) ' - ' - Q l O i - C M f O C M — " i co to ii) CM oi n S CO —^ CO Tf Tf o m o to o i r d ^ r P 6 o o o o tD CO Tf CO t - S N l i ) CO N CO O CM to to in ^ . CO ^ -CM CO Tf O CO tO CM — — CO O f-- CM CM CO ^ r» O) O Tf t z = £ E to t 9-^ i o o o - o -o E co gs™ ! if fr.E ' ; N © E ! ^ CD in fl} j s 2 ~$ : Tf CM m - cS i n CO o CM cj> O) Tf CD CD CD CD CO CO CD Tf i n CO o d p p p O p o O p o o o o o o o d d d d d d d d o" d d d d d d d d -cMcoTfinintoini^incoincDinO' • c M c o T f i n i n o m m m o m o m o m o o o o o o ^T-T-T-^CM^CM^C0C0TfTfUlU)tDr^C0tJ)OCN \— —) i— \— l i . 'Btje. 13.1 m co i n CN co r-t o O fl ^ l O ' - i n C f l t N t ' i - f O S l N t I N O tO O I Si. CO ID _ _ . _____ _ . co j i o i o o o a i a i o i T - T - N ' o o p o ^ ° o o o 0 - o o ( o o o o o o o d o " • ( M T f - S . K C O C D C O C O O ) O J i ^ ^ c n c o _ s o q } i - t N O i . I ' - o i n s t c f t - i o c D M n s ^ » T j - S C O C N 1 ^ S C O C O ( N S a i C N ) _ - _ ( M m s r o m t i 9 0 T - - c o r - f M C N W c n T t i n Q t o c O N S N S . . . . . - . Tf CN O Ol s o - s o r o ^ - o i i o c o N T - i n t ^ c d Q C N i n c f i s c o o m c o a . o i a ) , - ; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 e _ 8° O O O O O 0 0 0 0 0 (DOlCD^ CNnTl-CD^ CD-2 ^ i . T - o 3 S T - r t i n i o N iu 1 1 u u i n n r. !*• iu O) N S N S IO (•) IO T - (D N s i n i o i n i a M O ' - . i D n N Q j Q q i o i O ' -OO r^t^ S^COCNS^OlOcOlSScSQ -^CNI OOOOOO^T-CNCN°'f^vU.inifi.tD_ ^ ._; O O O O O O O O O . . . ) S I D ^ N J O IT) CJ) O) fN CN CO 5 CM o m co (0 n N fl'iONCOfflCOCO'-'-CN oi o co i n o N co co in in <o 0 0 o o o o o ,_; o O O O O O O o ' O O ( N S ' - C O ^ S O r t M ( 0 ( 0 - - 0 ) S N T O O I i ^ C D C O ( N i n O ) t N O . ' * C N C 7 ! i r \ j i f i r o c o ? ' - 0 ) Q ' t Q ( 0 ' © o i n S o o o i - ^ . o i I O S T - I J J ^ ' - C N O J S C O C O i o co o o o T - T - ^ d d d E .!_> 11 or c • E « x ~ aj £ i3 ° " M • H H - H S S S 5 2 c o c N c o r t i n i n T j T j c o r - c o o ) h . c N i n c D c j ) n Q m c o _ _ _ 0 ) i o i i O S ( D r o p ) c J i S ^ \ t ^ C D C 0 T - - T j C N N ( D l O 5 - C 0 C O O . l D O O N N . C N - i O C S O S n O r j O O n r ; O O O C O V i f ix ' — • O CD I • i n ! m 1 N CJ) W CO h- P"- S- N-CN CO CO CO CD CO Tf CN • s c D ' - c o c n c D i n c T i - N N c o c p s i n i n i D c o l O a j . o c N T - t t N - m a j i n i n c o c N i o c N i f l l O C s i m c O C O C O i n c O O T O T O T C D C O C n C N O l C N C O i c o T - c o c o o T - i n x t O T - T - o o T - o t n o o • ™ " • • • • - 0 ^ & ^ . : o o o o 0 0 0 0 0 0 0 0 0 i n T - i n t N i n c o c o c o a > c O T - < i Q t O c p _ N Q C 7 ) i n c O T - T j i CD CO r- CO CD < S r ^ j i j N o i • C D T - _ T r - C J ) S m C N C O f 0 C O 1 , ( ' • " - c o ^ - T - s c o o i o r - ^ ' • ' CD S S O) - CM ' O if) Q Q T-l TJ- O) ift _ - 00 — E COCOCDCOT-ChcOCOTfcOr-OlCOCbTtQ O V O _ 0 ) O O r ( N C O T f ^ l O w N _ 0 ) 0 ) - n T t C J ) ( O S ^ l n O T d d o o o o d o o d d d o o o o 6 o o o 6 6 o d d c i d ' : ' - ' ' - V S co m w co co F , co m i n F -_. M t r ra c o 1- -3 u. o m < 3 ^ 1 co o • 2 a i TtCOCMCOCOTtTj-T-T-CO Q O O O O O O O Q O v , . O O O O O O O O O O O T - O C N C N Q d d d d d d d d d d d d d d d ^ T - s c o c o o j o i s - Q t ^ ^ ^ ^ w c o ^ c N o s c o c D ^ ^ c o i n M i n T t T - h - r ^ c o < N T f C M T - ^ c o i n r ^ c o T - T - o o c N O T - o o o o o Q o o o T - O C N ^ C N ^ c ^ C N ' - 0 0 0 0 0 0 0 0 0 ' 0 0 0 0 0 0 0 0 0 0 d d d d d d d d d d d d d d d d d d d d j i n c D i n r ^ i n c o i n o i i n o i n T - i n c N c o T f i o i n o i n i n i n g i n o w i r i t o 1^ co o i ^ d ^ r - ' ' - ^ ^ ^ s ^ c N N N " n , r , T i B u l 5 a : i 3 i o r j , f l 1 3 T~ T- T- CN CM 1-I O I ^> , bo s , K •2 a Ii 3 to 2 o > o 13 cn co r- m Tt co ( N U) r S T- T- 1 s t 01 S O (O < § § 8 5 I g i >54 0.179779 68 0.223452 !28 0.244391 81 0.267125 !02 0.280287 81 0.295842 18 0.31379 )53 0.338917 17 0.371223 45 0.401137 141 0.421478 33 0.457074 144 0.510918 0.109! 0.1361 0.148! 0.1627 0.170! 0.1 so; 0.191; 0.20C 0.226; 0.244< 0.256S 0.278! 0.311: 0 0 0 0 0 0 0 1 0 3 ( o o lotoroMoni-ooin j c M c o T r C M T i - t o t o c n T f ScoT-irJT-cgiQT-cNco ) T f C M C O C M C O C O » O C M ) ( O T - T f m o ) « - ( O N T -I N Q C M m f T - C J r-T-; d d d >- fg d d d 5 >. S 1: Q 2 ( to cu 1 A 8 o <} ^ 1 to o m w co c _ , , s c o i - o i f O i i n i O ( D i 3 C B ( 5 0 C I ) N T - T - N < V l t ' - N ? O M n O i n i O C O N O i f l f ' - C J8< i ( O N O C O S C O » - ' - S ^ C O C O C H i n r t l ._ _. . _ _ . . _. -iroionocnnnnfflnoJOifJcooic,.. ., ^r^NOT-nOQSINSN^Tft^^^ffl'O^W^NNinNNNCoSTrS^t S C N O C r C N T t C O O T C O l O C N T r C D ^ l O C O T r l Q C O c O C M M . 9 6 "7 H <H « ^ » d l f > 0 ' T , r > « « c ? o c ? d ^ d ^ d d d d o o ' - r i ' P r i P r i O O 0 * - T - T - T - T - 0 0 mmmmm m -~. " O O O O O IS -I cn 4^ E in I S 5 c Z. o 3 ^ Is ncBcoNU)0)Nioo)T-ioioi-ti*NTfT-iDNC«o>-snu)NraT-nioiovcon Nio^ nsT-(flQTfCTirtSNcc)ihnc\ifflcnT-r>i5tosO)N nosTt r>T-cBcoco5co^g to i r )NC»T -N inNcn^ " s ( D 0 3 S N i o c » c o N ( N i o o ( j i s i n T f o g « - 1 5 t T - c o c N t . - Tf CO i n CSI CO I Q . . _ . i C O C O O r t C O C O O T - C O l O C O , . . . . - O O T - T - T - T - C M C M C N C N C N C N C O t O r t i o d o d o o d o o d d o o o o d o ' o O O O O O O 1 - 1 —> CO Tt CO w m in co -1 CM T- co OCM CO T--» co 3 m 5 o co cn j= Tt cd cd > 2 *- o V CO g> C S. 2 • 4 o ! S I 3 a. =: 1 1 i l l l i s Sf8 I sis j : ^ 5 i n a> OT-T-<Ncn^ csii^ cncoincDCVi<ocsiT-r^ ini^ r^ cou^ QbQNCpiDIOCOSrtONCNONCSINT-T-T-r-^ i-t-^ O O O O O O O O O O O O O O Q O O O O O O O O O O O O O O O O O O O O O O O d o d o o d d d d d d d d o ' d d d d d d d d d d d d d d d d d d d o d o < t-si ~ O c B > T J g I 8 5 r - i n O) CO CN tO cn csi to m in i CJ o o d d = £ ~ -S s — O t - CO U 3 ^ C O C O L O C O l O N T t O L Q C O C O ^ C O C 7 > C n O > a ) r ^ C O O t O C Q C O ^ r ^ ^ C n i O ^ U j O O O O J O C O C O o t v ^ Q O « f f i ^ < ^ i n s ^ o c \ i T - T - c » ^ i 5 w ^ c o c p T £ C p C N C N i n C N r ' u l U 1 t O T r ^ C 0 r ^ Q T - ^ ^ C N r ^ j O O O O M" , J O C O C N t n O O O C O C N O O Q O ^ N c O C O ^ t O C O ' i - O O C N O O O O O O O O O O O ! o o ci ci oi ^ N r i d d d d c i r - ' c i d d i d d o d o o d d o d d 1 0 0 0 0 0 0 0 0 0 0 0 ' CO ID CO CO CO s s co ^ rt i- cn 10 , . C C Q ( Q 9 N ( O ^ ^ C O C V | S O ^ C D C N V W < Q ( O V C N [ N ( D < a , C N O ) P > •= O C O r t Q N l f l C Q W Q Q > S N ^ N » - C » © ^ ^ C O i n .2 oi f ^ i o c O Q T f c o N t o o i c o Q ^ i n c n c N ^ T r c N C ^ > c 3 Q 5 Q ( D S ^ C O U ) r t ! D O O I f l W Q ' * W S ' - ( D C \ l ' - i S o i n Q t 01 cu m c o o o c o i o ^ c j j o c y T r c N r ^ c n ^ c o c o c d ~ c r O O ^ ^ ^ Q T - N N N Q W N f O c O Q r t T t ^ ^ Q i n c B t o i D s s N o a o N n S £ •=: 0 d 0 0 0 d o d o 0 0 o 0 d o d o d 0 0 0 0 — — — — . • -^.10 i - . IO Q ^ ( M b -^ CO O (D tO j o w , • r " O " I" T -£ o * § P» c re D CO -r- CN E a ^ CO - m 4 2 - « CM 2 m j 5 S o K • O O I J csi csi C\i CO CO ' S C 0 C S I C N l 0 0 ) T - t 0 C N C 0 C Q t Q C 0 O ^ » * l O l O 0 3 O ^ T f C N C Q I ^ O o o o ^ t D c o T r T f o c o K o o c N r ^ o o o o ^ ^ N c o o ^ c s j ^ o o ^ o o o o o o o O O O O Q O C M C M ^ C N O O O O O O O O O O O O O O O O O O O O O O O O O O O O d o d o d d d d d o d d d o d d d d d d d d d d d d d d d d d d d d o o o co co to co 0 0 0 0 0 0 0 0 d o d o ' T - c M c o ^ i n i f l m < o w r » i o c o i n o > i o O ' - N n ' m s co cn o O T - CN CO lO l CN CM CSI CN CN t _ 5 _ R Q R •2 §.1 E 3 II - r ^ T - S S c ^ c o c S c M S c N c o ^ S . O - r ^ v o i D O i C N o p T f T - T - T - o O T - q . • t o i n c o c j j T - r t f M t D N W g ^ n i o c o c n c - ' - - ~ • o o o o o o o o o o o o o o • • • t « Q O - S C » T - f f l t O n < D C g - r - r t S Q Q Q i d d d 6 d d d d d d d d d d d d c i d d d d c i d c i c i c i c i c i d d d d d d S . I _ _ ° 8 8 l S o o o <3 _ s _ a % o _ o Si a R S >5 £ S T J S S> ro I I •t z Study 1-99 i l aj = o ^ o _• « _ £5 r i 1- « 1- 1-_ = 3 d -c £ " •5 4? » M 0) _ 3 o t O _ ' f 1 1 _ I t ^ c y c N T f o T j - t p c S c n o o e i - C T j o o ^ f f l ^ i n s c » s c o i t ) c o s _ _ i n i O i r . * T f T f ' j d d d d d d d d o d d d d d d d d d d d d o > o ift c 2 * " 9 t t " ' " < " N ' - - J - - - - - - - - - -d d d d d d d d o d d d d d d d o d d ' c - N n c j ) | g N < D - s N c o a c 5 i Q t t o ^ s n » t — o r ^ T f T - o i c g r ^ « T i - e N T - o i o o t n .2 'ffl ^ Q i n « ( D C 0 l _ ) C n Q T - 7 - C N f 5 3 > t= T f M c N c S t o i n o r ^ i n c N m <i) _ r t _ o t - r t - s e o o N n S s f t r t c s i * i g s c » T - 4 F E c o O T - T - T - ^ T - - c g N N c \ i c N i N 0 , e o r t ( o n n ^ ^ * i n i n i n Q p ° ° ° ° ° d d d d d d d d d c i d c i c i c i c i d c i c i c i c i c i c i d d d - r 1 ^ ^ - ^ c\i co to ( _ t N < o t o c o ( p c p c o < p r ^ i n T t o D C M t N c o g Q c s i o w c o f i n ^ n n i o f i f l M - c o M n w N N W N T - N O O Q O O C D O O O O O C D O O O O O O O O O O O O O Q O O O O O O O C o o o o o o o o o o o o o o o o o o o o o o O O O O ( ) - - - CO s - T - O o o o > 9 d 9 d 9 9 > o o o o - c M n u i ^ i / ) i n i n _ i O N i / ) e o ^ m i f i o w ^ t n c M c o ^ i / ) ( D N c o c j ) o i n i n i n o ^ co i r i co N: OO A I ' " ° ' " ^ T " ^ ' " ' " T * ' " ' " ' " ( N ^ R I ( ^ " " Q - l Q C O O O Q O O p Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 135 APPENDIX F: RTD Parameter Macro V B A Subroutine in Excel to Calculate the Dimensionless Cumulative Curve (F-Curve) from a dimensionless C curve. Also calculates the important parameters for RTD comparison. Sub DimensionlessStatsCalc() i * * * * * 'May, 1999 - RWJ 'Subroutine that c a l c u l a t e s the F-Curve from C-curve data 'C-curve data must be dimensionless and i n two adjacent columns 1 S e l e c t the columns before running the macro 'The output w i l l be generated i n columns t o the r i g h t i * * * * * Dim i As Integer Dim j As Integer Dim sglTimeO As S i n g l e Dim sglConcO As S i n g l e Dim sglCumulative() As S i n g l e Dim intNumRows As Integer Dim intActiveWriteRow As Integer Dim intWriteColumn As Integer Dim sglTracerMass As S i n g l e Dim sglDeltaT As S i n g l e Dim sglTIO As S i n g l e Dim sglT90 As S i n g l e Dim sglTPeak As S i n g l e Dim sglCPeak As S i n g l e Dim sglT50 As S i n g l e Dim sglFirstMoment Dim sglMean As S i n g l e 'Marecos do Mara V a r i a b l e s Dim sglQ As S i n g l e Dim sglTbar As S i n g l e Dim sglSigmaSquare As S i n g l e Dim rgRead As Range Set rgRead = S e l e c t i o n intWriteColumn = rgRead.Column + 3 intActiveWriteRow = rgRead.Row intNumRows = rgRead.Rows.Count ReDim sglTime(intNumRows) ReDim sglConc(intNumRows) ReDim sglCumulative(intNumRows) 'Populate Array sglTPeak = 0 sglCPeak = 0 For i = 1 To intNumRows sglTime(i) = ActiveSheet.Cells(rgRead.Row + i - 1, rgRead.Column) sglConc ( i ) = ActiveSheet.Cells(rgRead.Row + i - 1, rgRead.Column + 1) I f sglConc(i) > sglCPeak Then sglCPeak = sglConc(i) sglTPeak = sglTime(i) End I f Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 136 Next i 'Check Dye Recovery and Generate Cumulative D i s t ' n sglTracerMass = 0 sglFirstMoment = 0 For i = 1 To intNumRows I f i = 1 Then sglDeltaT = ( s g l T i m e d + 1) + sglTime(i)) / 2 sglC u m u l a t i v e ( i ) = sglDeltaT * sglConc(i) E l s e l f i = intNumRows Then sglDeltaT = sglTime(i) - (sglTime(i - 1) + sglTime(i)) / 2 El s e sglDeltaT = (s g l T i m e ( i +1) + sglTime(i)) / 2 - (sglT i m e ( i - 1) + sg l T i m e ( i ) ) / 2 End I f sglTracerMass = sglTracerMass + sglDeltaT * sglConc(i) I f i > 1 Then sgl C u m u l a t i v e ( i ) = s g l C u m u l a t i v e ( i - 1) + sglDeltaT * sglConc(i) sglFirstMoment = sglFirstMoment + sglDeltaT * sglConc(i) * sglTime(i) Next i sglMean = sglFirstMoment / sglTracerMass A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "Dye Recovery:" A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglTracerMass 'Write Cumulative Function A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w - 1, intWriteColumn + 2) = "Cumulative" A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w - 1, intWriteColumn + 3) = "Cumulative Corrected" For i = 1 To intNumRows A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w + i - 1, intWriteColumn + 2) = sglCumulative(i) A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w + i - 1, intWriteColumn + 3) = sglCumulative(i) / sglTracerMass Next i intActiveWriteRow = intActiveWriteRow + 1 'Check f o r t l O and t90 For i = 2 To intNumRows I f s g l C u m u l a t i v e ( i - 1) / sglTracerMass < 0.1 And sglCumulative(i) / sglTracerMass > 0.1 Then 'Found t l O - i n t e r p o l a t e t o f i n d time sglTIO = s g l T i m e ( i - 1) + (sglTime(i) - s g l T i m e ( i - 1 ) ) * (0.1 - s g l C u m u l a t i v e ( i -1)) / (sglCumulative(i) - s g l C u m u l a t i v e ( i - 1 ) ) E l s e l f s g l C u m u l a t i v e ( i - 1) / sglTracerMass < 0.5 And sglC u m u l a t i v e ( i ) / sglTracerMass > 0.5 Then 'Found t50 - i n t e r p o l a t e t o f i n d time sglT5 0 = s g l T i m e d - 1) + (sglTime(i) - s g l T i m e ( i - 1)) * (0.5 - sglCumulative ( i -1)) / (sglCumulative(i) - s g l C u m u l a t i v e ( i - 1)) E l s e l f s g l C u m u l a t i v e ( i - 1) / sglTracerMass < 0.9 And sglC u m u l a t i v e ( i ) / sglTracerMass > 0.9 Then 'Found t90 - i n t e r p o l a t e t o f i n d time sglT90 = s g l T i m e d - 1) + (sglTime(i) - s g l T i m e d - 1 ) ) * (0.9 - sglCumulative ( i -1)) / (sglCumulative(i) - s g l C u m u l a t i v e ( i - 1)) End I f Next i Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 137 'Calculate the variance based on the Marercos do Mara method 'There's a whole f a m i l y of v a r i a b l e s here but b a s i c a l l y tbar, Q and sigmasquare 'See the paper f o r d e t a i l s ' Calculate Q For i = 2 To intNumRows sglQ = sglQ + ((sglConc(i) + s g l C o n c d - 1 ) ) / 2 * (sglTime(i) - s g l T i m e d - 1))) Next i ' c a l c u l a t e Tbar For i = 2 To intNumRows sglTbar = sglTbar + ((s g l T i m e d ) + s g l T i m e d - 1 ) ) / 2 * (sglConc(i) + s g l C o n c d - 1 ) ) / 2 * (sglTimed) - s g l T i m e d - 1))) Next i sglTbar = sglTbar / sglQ ' c a l c u l a t e SigmaSquare For i = 2 To intNumRows sglSigmaSquare = s g l SigmaSquare + ( ( ( (sglTime (i) + s g l T i m e d - 1 ) ) / 2) * 2) * ( (sglConcd) + s g l C o n c d - 1 ) ) / 2) * (sglTime(i) - s g l T i m e d - 1))) Next i sglSigmaSquare = (sglSigmaSquare / (sglQ * (sglTbar A 2))) - 1 'Output A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "Mean" Ac t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglMean intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "CPeak" Ac t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglCPeak intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "TPeak" Ac t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglTPeak intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "TlO" A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglTIO intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "T50" Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglT50 intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "T90" Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglT90 intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = " M o r r i l l Index" A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglT90 / sglTIO intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "Max Calc Time" Ac t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglTime(intNumRows) intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "MOM Q" Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglQ intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "MOM Tbar" A c t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglTbar intActiveWriteRow = intActiveWriteRow + 1 Act i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn) = "MOM SigmaSquare" Ac t i v e S h e e t . C e l l s ( i n t A c t i v e W r i t e R o w , intWriteColumn + 1) = sglSigmaSquare intActiveWriteRow = intActiveWriteRow + 1 End Sub Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 138 APPENDIX G: Scalar TRANSPORT Code Scalar Transport Code Option E x p l i c i t 1 MAIN HOUSING CODE FOR TRANSPORT ROUTINE P u b l i c Sub ChemMixMainl(iNx As Integer, iNy As Integer, _ dbldt As Double, numiter As Long, _ d b l V e l O As Double, dblConc () As Double, Comments As S t r i n g , _ dblDX As Double, dblDY As Double, dblDeltaX As Double, dblDeltaY As Double, _ dblMaxX As Double, dblMaxY As Double, i n t E x p o r t l n t e r v a l As Integer, i n t T r a c e r l n t e r v a l As Integer, _ intOutBottom As Integer, intOutTop As Integer) Dim dblDeltaConc() As Double 'Update array sent to ApproxFactor Dim d b l F I O As Double Dim dblAxf) Dim dblAyO Dim dblBxO Dim dblByO Dim dblCxf) Dim dblCyO As Double As Double As Double As Double As Double As Double Dim i As Long Dim j As Integer Dim dtmLastTime As Date Dim dtmNewTime As Date dtmLastTime = Now() 'Dimensioned f o r dblVel(1,x,y):= U and dblVel(2,x,y) := V ReDim dblDeltaConc(iNx, iNy) ReDim d b l F I ( i N x , iNy) ReDim dblAxfiNx, iNy) ReDim dblAy(iNx, iNy) ReDim dblBx(iNx, iNy) ReDim dblBy(iNx, iNy) ReDim dblCx(iNx, iNy) ReDim dblCy(iNx, iNy) 'Because the problem i s steady s t a t e wrt V e l and Disp the ABC f a c t o r s ' w i l l not change from i t e r a t i o n t o i t e r a t i o n . ' Best to c a l c u l a t e them and store them f o r future 1 use. CalcABCFactors d b l V e l ( ) , dblDX, dblDY, dblDeltaX, dblDeltaY, d b l d t , iNx, iNy, dblAx () , d b l B x O , d b l C x f ) , d b l A y O , d b l B y O , dblCyO E x p o r t S u r f e r G r i d F i l e C o n c "Oconc.grd", iNx, iNy, dblConcO, dblMaxX, dblMaxY frmTransportl.MousePointer = 11 'Now do i t e r a t i v e l i n e solve For i = 1 To numiter 'Do approx f a c t o r ApproxFactor dblDeltaConc(), iNx, iNy, dblAx, dblAy, dblBx, dblBy, dblCx, dblCy, d b l F I 'Add update ma t r i x t o s o l u t i o n AddMatrix2 dblDeltaConc(), dblConcO, iNx, iNy 'Update Bndy Conditions Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 139 SetBoundaryConditions dblConcO, iNx, iNy 1 C a l c u l a t e f l u x i n t e g r a l C a l c u l a t e F l u x I n t e g r a l dblConcO, d b l V e l 0 , db l F I () , dblDeltaX, dblDeltaY, dblDX, dblDY, d b l d t , iNx, iNy I f i Mod i n t E x p o r t I n t e r v a l = 0 Then Exp o r t S u r f e r G r i d F i l e C o n c T r i m ( S t r ( i ) & "cone .grd") , iNx, iNy, dblConcO, dblMaxX, dblMaxY End I f I f i Mod i n t T r a c e r l n t e r v a l = 0 Then MassInTank "Tracer.csv", iNx, iNy, dblConc, i * dbldt ConcAtOutlet "TracerConc.csv", iNx, iNy, dblConc, i * d b l d t , intOutBottom, intOutTop End I f I f i Mod 20 = 0 Then dtmNewTime = Now() dtmLastTime = (dtmNewTime - dtmLastTime) / 10 frmTransportl.txtETA = Format((numiter - i ) * dtmLastTime + dtmNewTime, "Long Time") dtmLastTime = dtmNewTime End I f f r m T r a n s p o r t l . t x t l t e r = i DoEvents Next i frmTransportl.MousePointer = 0 End Sub 'Subroutine that sets up the flow up f o r the run button P u b l i c Sub FlowSetupRun(iNx As Integer, iNy As Integer, _ dbldt As Double, numiter As Long, _ InBottom As S i n g l e , InTop As S i n g l e , OutBottom As S i n g l e , _ OutTop As S i n g l e , strComments As S t r i n g , _ strlnputVelFileName As S t r i n g , dblDX As Double, dblDY As Double, _ dblMaxX As Double, dblMaxY As Double, i n t E x p o r t l n t e r v a l As Integer, i n t T r a c e r l n t e r v a l As Integer) Dim dblDeltaX As Double Dim dblDeltaY As Double Dim db l V e l ( ) As Double Dim dblConcO As Double Dim i n t l n B o t t o m As Integer Dim i n t l n T o p As Integer Dim intOutBottom As Integer Dim intOutTop As Integer dblDeltaX = dblMaxX / iNx dblDeltaY = dblMaxY / iNy ReDim d b l V e l ( 2 , iNx + 1, iNy + 1) ReDim dblConc(iNx + 1, iNy + 1) intlnBott o m = I n t ( i N y * InBottom / 100 +0.5) + 1 i n t l n T o p = I n t ( i N y * InTop / 100) intOutBottom = I n t ( i N y * OutBottom / 100 + 0.5) + 1 intOutTop = I n t ( i N y * OutTop / 100) 'Load the V e l o c i t y data f o r f u t u r e use... LoadVelData strlnputVelFileName, iNx, iNy, d b l V e l P l u g l n j e c t i o n O f T r a c e r D a t a iNx, iNy, dblConcO, in t l n T o p , intlnBottom 'NEED TO SET BOUNDARY VELOCITY AND CONC DATA FOR FLUX CALCS Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 140 SetVelocityBoundaryConditions d b l V e l , iNx, iNy, intlnBottom, i n t l n T o p , intOutBottom, intOutTop SetBoundaryConditions dblConc, iNx, iNy ChemMixMainl iNx, iNy, d b l d t , numiter, d b l V e l , dblConc, strCoraments, dblDX, dblDY, dblDeltaX, dblDeltaY, dblMaxX, dblMaxY, i n t E x p o r t l n t e r v a l , i n t T r a c e r l n t e r v a l , intOutBottom, intOutTop End Sub 'Subroutine that sets up the flow f o r the Checker Button ' This one I manipulate f o r output P u b l i c Sub FlowSetupRunChecker(iNx As Integer, iNy As Integer, _ dbldt As Double, numiter As Long, _ strComments As S t r i n g , dblDX As Double, dblDY As Double, _ dblMaxX As Double, dblMaxY As Double, i n t E x p o r t l n t e r v a l As Integer, _ dblVx As Double, dblVy As Double, i n t l t a r g e t As Integer, i n t J t a r g e t As Integer, dblMass As Double, i n t T r a c e r l n t e r v a l As Integer) Dim dblDeltaX As Double Dim dblDeltaY As Double Dim d b l V e l ( ) As Double Dim dblConcO As Double Dim i As Integer Dim j As Integer Dim i n t l n B o t t o m As Integer Dim i n t l n T o p As Integer Dim intOutBottom As Integer Dim intOutTop As Integer dblDeltaX = dblMaxX / iNx dblDeltaY = dblMaxY / iNy ReDim d b l V e l ( 2 , iNx + 1, iNy + 1) ReDim dblConc(iNx + 1, iNy + 1) intlnBott o m = 1 in t l n T o p = iNy intOutBottom = 1 intOutTop = iNy GenerateDummyVelData iNx, iNy, d b l V e l , dblVx, dblVy GenerateDummyConcData iNx, iNy, dblConc, i n t l t a r g e t , i n t J t a r g e t , dblMass 'NEED TO SET BOUNDARY VELOCITY AND CONC DATA FOR FLUX CALCS SetVelocityBoundaryConditions d b l V e l , iNx, iNy, intlnBottom, i n t l n T o p , intOutBottom, intOutTop SetBoundaryConditions dblConc, iNx, iNy ChemMixMainl iNx, iNy, d b l d t , numiter, d b l V e l , dblConc, strComments, dblDX, dblDY, dblDeltaX, dblDeltaY, dblMaxX, dblMaxY, i n t E x p o r t l n t e r v a l , i n t T r a c e r l n t e r v a l , 1, iNy ExportFlowDetailsCheck iNx, iNy, dblDX, dblDY, strComments, dblVx, dblVy, numiter, i n t E x p o r t l n t e r v a l , d b l d t , dblMaxX, dblMaxY, i n t l t a r g e t , i n t J t a r g e t , dblMass End Sub 'Subroutine to load the v e l o c i t y data i n t o a 3D matrix ' the f i r s t index of d b l V e l i d e n t i f i e s U (=1) or V (=2) P u b l i c Sub LoadVelData(strFilename As S t r i n g , iNx As Integer, iNy As Integer, d b l V e l 0 As Double) Dim i As Integer Dim j As Integer Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 141 Dim k As Integer Dim intLocateUV As Long Dim dblDummy As Double ' F i l e has already been found so no e r r o r checking required. Open strFilename For Binary Access Read As #1 'Must Seek the l o c a t i o n of the U and V data 1 i . e . s k i p over the P data dblDummy = iNx * iNy dblDummy = 1 2 9 + 1 2 8 + 8 * dblDummy intLocateUV = Int(dblDummy) Seek #1, intLocateUV For k = 1 To 2 For i = 1 To iNx For j = 1 To iNy Get #1, , d b l V e l ( k , i , j) Next j Next i Next k Close #1 End Sub 'Subroutine t o set the Ghost C e l l boundary c o n d i t i o n s P u b l i c Sub SetBoundaryConditions(dblConc() As Double, iNx As Integer, iNy As Integer) Dim i As Integer Dim j As Integer 'Do LHS ' Assume no d i f f u s i v e f l u x across i n l e t i = 0 For j = 1 To iNy db l C o n c ( i , j) = d b l C o n c ( i +1, j ) Next j 'Do RHS ' Assume no d i f f u s i v e f l u x across o u t l e t i = iNx + 1 For j = 1 To iNy db l C o n c ( i , j) = d b l C o n c ( i - 1, j) Next j 'Do Bottom j = 0 For i = 1 To iNx db l C o n c ( i , j) = d b l C o n c ( i , j + 1) Next i 'Do Top j = iNy + 1 For i = 1 To iNx db l C o n c ( i , j) = d b l C o n c ( i , j - 1) Next i End Sub P u b l i c Sub C a l c u l a t e F l u x I n t e g r a l ( d b l C o n c ( ) As Double, dblVel() As Double, _ d b l F I O As Double, dDeltaX As Double, dDeltaY As Double, _ dblDX As Double, dblDY As Double, dbldt As Double, iNx As Integer, iNy As Integer) Dim i As Integer Dim j As Integer Dim d b l F F l u x O As Double Dim dblGFluxO As Double ReDim d b l F F l u x ( i N x + 1, iNy + 1) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 142 ReDim dblGFlux(iNx + 1, iNy + 1) Ca l c u l a t e F l u x e s F i r s t U p w i n d True, dblConcO, d b l V e l 0 , d b l F F l u x O , dDeltaX, dDeltaY, dblDX, dblDY, d b l d t , iNx, iNy Ca l c u l a t e F l u x e s F i r s t U p w i n d False, dblConcO, d b l V e l 0 , d b l G F l u x ( ) , dDeltaX, dDeltaY, dblDX, dblDY, d b l d t , iNx, iNy ' CalculateFluxesSecondUpwind True, dblConcO, d b l V e l O , d b l F F l u x O , dDeltaX, dDeltaY, dblDX, dblDY, d b l d t , iNx, iNy 'CalculateFluxesSecondUpwind False, dblConcO, d b l V e l 0 , d b l G F l u x O , dDeltaX, dDeltaY, dblDX, dblDY, d b l d t , iNx, iNy 'CalculateFluxesSecondUpwindTVD True, dblConcO, d b l V e l 0 , d b l F F l u x O , dDeltaX, dDeltaY, dblDX, dblDY, d b l d t , iNx, iNy 'CalculateFluxesSecondUpwindTVD False, dblConcO, d b l V e l 0 , dblGFlux () , dDeltaX, dDeltaY, dblDX, dblDY, d b l d t , iNx, iNy For i = 1 To iNx For j = 1 To iNy d b l F I ( i , j) = - d b l F F l u x f i + 1, j) + d b l F F l u x ( i , j) - d b l G F l u x ( i , j + 1) + d b l G F l u x ( i , j) Next j Next i End Sub P u b l i c Sub C a l c u l a t e F l u x e s C e n t r e D i f f ( b o l V e r t i c a l As Boolean, dblConcO As Double, d b l V e l 0 As Double, _ d b l F l u x O As Double, dblDeltaX As Double, dblDeltaY As Double, dblDX As Double, dblDY As Double, dbldt As Double, iNx As Integer, iNy As Integer) 'Subroutine to c a l c u l a t e d i r e c t i o n a l f l u x e s 1 both h o r i z o n t a l (F) and v e r t i c a l (G) code i s in c l u d e d here 1 and the switch to choose the appropriate path i s ' b o l V e r t i c a l ' 'Generate Flux I n t e g r a l Matrices f o r n timestep 'The 1/2 increments are rounded up i n determining the i n t e g e r 'for the f l u x index. ' eg: F 2+1/2,4 => F(3,4) ' eg: G 3,1-1/2 => F( 3 , l ) 'Fluxes employ centred d i f f e r e n c e convection and 2nd order accurate d i f f u s i o n ' Flux accuracy confirmed on March 8, 1999 3:40pm Dim i As Integer Dim j As Integer Dim i U As Integer Dim i V As Integer i U = 1 i V = 2 I f b o l V e r t i c a l = True Then 'Calculate V e r t i c a l Flux M a t r i x 'NEED TO RUN THIS TO INX+1 AND INY+1 BUT NEED TO SET VELOCITY VALUES OUTSIDE THE STANDARD INX,INY MATRIX For j = 1 To iNy + 1 For i = 1 To iNx + 1 'Connective Flux d b l F l u x ( i , j ) = dbldt / dblDeltaX / 2 * ( d b l V e l ( i U , i , j) * d b l C o n c f i , j) + d b l V e l ( i U , i - 1, j) * dblConc(i - 1, j ) ) ' D i f f u x i v e Flux Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 143 d b l F l u x d , j ) = d b l F l u x d , j) - dbldt * dblDX / (dblDeltaX) * 2 * (dblConc(i, j) - d b l C o n c ( i - 1, j ) ) Next i Next j El s e 'Calculate H o r i z o n t a l f l u x M a t r i x For i = 1 To iNx + 1 For j = 1 To iNy + 1 'Convective Flux d b l F l u x d , j) = dbldt / dblDeltaY / 2 * ( d b l V e K i v , i , j ) * d b l C o n c ( i , j) + d b l V e K i v , i , j - 1) * db l C o n c ( i , j - 1) ) ' D i f f u x i v e Flux d b l F l u x d , j) = d b l F l u x d , j) - dbldt * dblDY / (dblDeltaY) " 2 * (dblConcd, j) - d b l C o n c d , j - 1)) Next j Next i End I f End Sub 'Subroutine to set the v e l o c i t y data i n the boundaries t o p r o p e r l y ' c a l c u l a t e the f l u x e s 'Subroutine checked q u i t e thoroughly - seems okay. P u b l i c Sub SetVelocityBoundaryConditions(dblVel() As Double, iNx As Integer, iNy As Integer, _ intlnBott o m As Integer, i n t l n T o p As Integer, intOutBottom As Integer, intOutTop As Integer) Dim i As Integer Dim j As Integer 'DO top j = iNy + 1 For i = 1 To iNx d b l V e l d , i , j ) = - d b l V e l d , i , j - 1) db l V e l ( 2 , i , j ) = -db l V e l ( 2 , i , j - 1) Next i 'Do bottom j = 0 For i = 1 To iNx d b l V e l d , i , j ) = - d b l V e l d , i , j + 1) d b l V e l ( 2 , i , j) = -dblVel(2, i , j + 1) Next i 'Do l e f t i = 0 For j = 1 To iNy '** Make i n l e t l i k e a w a l l so we have no f l u x 'across the w a l l at the i n l e t ' I f j >= intlnBot t o m And j <= intl n T o p Then ' d b l V e l d , i , j) = d b l V e l d , i + 1, j) ' d b l V e l ( 2 , i , j) = db l V e l ( 2 , i + 1, j) 'Else d b l V e l d , i , j) = - d b l V e l d , i + 1, j) db l V e l ( 2 , i , j) = -dblVel(2, i + 1, j ) 'End I f Next j 'Do r i g h t i = iNx + 1 For j = 1 To iNy I f j >= intOutBottom And j <= intOutTop Then d b l V e l d , i , j) = d b l V e l d , i - 1, j ) db l V e l ( 2 , i , j) = db l V e l ( 2 , i - 1, j) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 144 E l s e d b l V e l d , i , j ) db l V e l ( 2 , i , j ) End I f Next j = - d b l V e l d , i - 1, j) = -db l V e l ( 2 , i - 1, j) End Sub 'Subroutine that c a l c u l a t e s the A, B, and C f a c t o r s ' f o r the approximate f a c t o r i z a t i o n subroutine P u b l i c Sub CalcABCFactors(dblVel() As Double, dblDX As Double, dblDY As Double, _ dblDeltaX As Double, dblDeltaY As Double, dblDeltaT As Double, iNx As Integer, iNy As Integer, dblAx() As Double, dblBx() As Double, dblCx() As Double, _ dblAyO As Double, dblByO As Double, dblCyO As Double) Dim i u As Integer Dim i v As Integer Dim i As Integer Dim j As Integer i U = 1 i v = 2 For i = 1 To iNx For j = 1 To iNy 'Active M o d i f i c a t i o n - RWJ J u l y 5, 1999 d b l A x d , j d b l B x d , j d b l C x d , j d b l A y ( i , j d b l B y d , j d b l C y d , j = dblDeltaT * (-dblVeldU, i - 1, j ) / 4 / dblDeltaX - dblDX / 2 / (dblDeltaX) * 2) * dblDeltaT = (1 + dblDX * dblDeltaT / (dblDeltaX) * 2) = dblDeltaT * ( d b l V e l ( i U , i + 1, j) / 4 / dblDeltaX _ - dblDX / 2 / (dblDeltaX) " 2) * dblDeltaT = dblDeltaT * (-dblVel(iv, i , j - 1) / 4 / dblDeltaY - dblDY / 2 / (dblDeltaY) " 2) * dblDeltaT = (1 + dblDY * dblDeltaT / (dblDeltaY) " 2) = dblDeltaT * ( d b l V e l ( i v , i , j + 1) / 4 / dblDeltaY - dblDY / 2 / (dblDeltaY) " 2) * dblDeltaT Next ] Next i End Sub P u b l i c Sub ApproxFactor(dblDeltaConc() As Double, iNx As Integer, iNy As Integer, dblAxO As Double, dblAyO As Double, dblBxO As Double, dblByO As Double, dblCxO As Double, dblCyO As Double, dblRO As Double) Dim dblAvect() As Double Dim dblBvectO As Double Dim dblCvectO As Double Dim dblRvect() As Double 'This group of vectors i s used t o 'carry ABC i n f o t o the Thomas Algorithm ' l i n e s o l v e r Dim i As Integer Dim j As Integer 'Do l i n e s o l u t i o n s by columns(constant i and advancing upward) ReDim dblAvect(iNy) ReDim dblBvect(iNy) ReDim dblCvect(iNy) ReDim dblRvect(iNy) For i = 1 To iNx Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 145 'Assemble the Vectors For j = 1 To iNy db l A v e c t ( j ) = d b l A y ( i , j) db l B v e c t ( j ) = d b l B y ( i , j) db l C v e c t ( j ) = d b l C y t i , j) dbl R v e c t ( j ) = d b l R ( i , j) Next j ' C a l l Thomas Algorithm w i t h the Vectors ThomasAlgorithm iNy, dblAvect, dblBvect, dblCvect, dblRvect 'Insert the s o l u t i o n from thomas i n t o intermediate array For j = 1.To iNy db l D e l t a C o n c ( i , j ) = dbl R v e c t ( j ) Next j Next i 'Do l i n e s o l u t i o n s by columns(constant i and advancing upward) ReDim dblAvect(iNx) ReDim dblBvect(iNx) ReDim dblCvect(iNx) ReDim dblRvect(iNx) For j = 1 To iNy 'Assemble the Vectors For i = 1 To iNx d b l A v e c t ( i ) = d b l A x ( i , j) d b l B v e c t ( i ) = d b l B x ( i , j) d b l C v e c t ( i ) = d b l C x ( i , j) d b l R v e c t ( i ) = db l D e l t a C o n c ( i , j) Next i ' C a l l Thomas Algorithm ThomasAlgorithm iNx, dblAvect, dblBvect, dblCvect, dblRvect 'Insert the s o l u t i o n from thomas i n t o f i n a l s o l u t i o n array For i = 1 To iNx db l D e l t a C o n c ( i , j ) = d b l R v e c t ( i ) Next i Next j End Sub 'Subroutine f o r s o l v i n g t r i d i a g o n a l matrix system of l i n e a r eqns. P u b l i c Sub ThomasAlgorithm(intlmax As Integer, A() As Double, B() As Double, C O As Double, R() As Double) Dim i As Integer 'Note t o s e l f : I t h i n k we have t o solve over 0 to N+l ' Get t h i s s t r a i g h t before coding.... 'Linear combinations of rows t o e l i m i n a t e a l l A's 'Scale each row t o make B=l For i = 1 To intlmax - 1 C(i ) = C ( i ) / B ( i ) R(i) = R(i) / B ( i ) B( i ) = 1 B ( i + 1) = B ( i + 1) - C(i) * A ( i + 1) R ( i + 1) = R ( i + 1) - R(i) * A ( i + 1) A ( i + 1) = 0 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 146 Next i R(intlmax) = R(intlmax) / B(intlmax) B(intlmax) = 1 'Back S u b s t i t u t e e l i m i n a t i n g C's 'S o l u t i o n contained i n R's For i = intlmax - 1 To 1 Step -1 R(i) = R(i) - R ( i + 1) * C ( i ) C ( i ) = 0 Next i End Sub 'Sub to add one matr i x t o another i n 2 dimensions P u b l i c Sub AddMatrix2(Source() As Double, Target() As Double, i x As Integer, i y As Integer) Dim i As Integer Dim j As Integer For i = 1 To i x For j = 1 To i y T a r g e t ( i , j) = T a r g e t ( i , j) + S o u r c e ( i , j) Next j Next i End Sub P u b l i c Sub GenerateDummyVelData(iNx As Integer, iNy As Integer, d b l V e l ( ) As Double, dblVx As Double, dblVy As Double) Dim i As Integer Dim j As Integer For i = 1 To iNx For j = 1 To iNy d b l V e l ( 1 , i , j) = dblVx d b l V e l ( 2 , i , j) = dblVy Next j Next i End Sub P u b l i c Sub Pl u g l n j e c t i o n O f T r a c e r D a t a ( i N x As Integer, iNy As Integer, dblConcO As Double, i n t l n T o p As Integer, i n t l n B o t t o m As Integer) Dim i As Integer Dim j As Integer Dim d b l V a l As Integer 'Introduce a dimensionless mass of 100 db l V a l = 100 / (intlnTop - intlnBottom + 1) For i = 1 To iNx For j = 1 To iNy I f i = 1 And j <= intl n T o p Then db l C o n c ( i , j) = d b l V a l E l s e d b l C o n c ( i , j) = 0 End I f Next j Next i End Sub Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 147 P u b l i c Sub GenerateDummyConcData(iNx As Integer, iNy As Integer, _ dblConcO As Double, i n t l t a r g e t As Integer, i n t J t a r g e t As Integer, dblMass As Double) Dim i As Integer Dim j As Integer For i = 1 To iNx For j = 1 To iNy I f i = i n t l t a r g e t And j = i n t J t a r g e t Then d b l C o n c ( i , j ) = dblMass E l s e d b l C o n c d , j ) =0 End I f Next j Next i End Sub P u b l i c Sub Ex p o r t S u r f e r G r i d F i l e C o n c ( s t r F i l e n a m e As S t r i n g , iNx As Integer, iNy As Integer, _ dConcMatrix() As Double, dMaxX As Double, dMaxY As Double) 1iPUV - 1 i s P, 2 i s U, 3 i s V 'Surfer F i l e S t r u c t u r e '4 bytes of s t r i n g data "DSBB" - p o s l 'Integer (2bytes) inX - pos5 'Integer (2bytes) inY - pos7 '2*double (16 bytes) minx, MaxX -posl5,23 '2*double (16 bytes) minY, MaxY -pos31,39 '2*double (16 bytes) minZ, MaxZ -pos47,55 ' i n X * i n Y * s i n g l e --array of data i n rows of constant Y s t a r t i n g from Y=l t o Y=inY -pos56 Dim dMinY As Double Dim dMinX As Double Dim dMinZ As Double Dim dMaxZ As Double Dim s4 As S t r i n g * 4 Dim s V a r i a b l e As S i n g l e Dim s t r F i l e C h e c k As S t r i n g Dim i , j , k As Integer dMinX = 0 dMinY = 0 For i = 1 To iNx For j = 1 To iNy I f d C o n c M a t r i x d , j) > dMaxZ Then dMaxZ = dConcMatrix(i, j ) I f d C o n c M a t r i x d , j) < dMinZ Then dMinZ = dConcMatrix(i, j ) Next j Next i 'Check to see i f strFilename e x i s t s -- i f so .. k i l l i t . 'Perhaps I should s t i c k t h i s check outside the r o u t i n e ... w e ' l l see s t r F i l e C h e c k = Di r ( s t r F i l e n a m e ) I f ( s t r F i l e C h e c k <> "") Then K i l l strFilename End I f Open strFilename For Binary Access Write As #1 'Write Header Info s4 = "DSBB" Put #1, , s4 Put #1, , iNx Put #1, , iNy Put #1, , dMinX Put #1, , dMaxX Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 148 Put #1, , dMinY Put #1, , dMaxY Put #1, , dMinZ Put #1, , dMaxZ 'Write Values For j = 1 To iNy For i = 1 To iNx s V a r i a b l e = dConcMa t r i x d , j) Put #1, , s V a r i a b l e Next i Next j Close #1 End Sub P u b l i c Sub C a l c u l a t e F l u x e s F i r s t U p w i n d ( b o l V e r t i c a l As Boolean, dblConcO As Double, d b l V e l () As Double, _ d b l F l u x O As Double, dblDeltaX As Double, dblDeltaY As Double, _ dblDX As Double, dblDY As Double, dbldt As Double, iNx As Integer, iNy As Integer) 'Subroutine to c a l c u l a t e d i r e c t i o n a l f l u x e s ' both h o r i z o n t a l (F) and v e r t i c a l (G) code i s included here ' and the switch t o choose the appropriate path i s ' b o l V e r t i c a l ' 'Generate F l u x I n t e g r a l Matrices f o r n timestep 'The 1/2 increments are rounded up i n determining the i n t e g e r ' f o r the f l u x index. ' eg: F 2+1/2,4 '=> F(3,4) ' eg: G 3,1-1/2 => F ( 3 , l ) 'Fluxes employ centred d i f f e r e n c e convection and 2nd order accurate d i f f u s i o n ' Flux accuracy confirmed on March 8, 1999 3:40pm Dim i As Integer Dim j As Integer Dim i U As Integer Dim i v As Integer Dim d b l V e l D i r As Double i U = 1 iV = 2 I f b o l V e r t i c a l = True Then 'Calculate V e r t i c a l Flux M a t r i x 'NEED TO RUN THIS TO INX+1 AND INY+1 BUT NEED TO SET VELOCITY VALUES OUTSIDE THE STANDARD INX,INY MATRIX For j = 1 To iNy + 1 For i = 1 To iNx + 1 'Convective Flux d b l V e l D i r = ( d b l V e l ( i U , i , j) + d b l V e l ( i U , i - 1, j ) ) / 2 I f d b l V e l D i r >= 0 Then 'Flow moves i n forward x d i r e c t i o n d b l F l u x d , j) = dbldt / dblDeltaX * (dblVel (iU, i - 1, j) * d b l C o n c d - 1, j ) ) E l s e 'Flow moves i n negative x d i r e c t i o n d b l F l u x d , j) = dbldt / dblDeltaX * (dblVel (iU, i , j) * d b l C o n c d , j ) ) End I f ' D i f f u x i v e Flux Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 149 d b l F l u x f i , j) = d b l F l u x f i , j ) - dbldt * dblDX / (dblDeltaX) * 2 * (dblConc(i, j) - db l C o n c ( i - 1, j ) ) Next i Next j Else ' C a l c u l a t e H o r i z o n t a l f l u x M a t r i x For i = 1 To iNx + 1 For j = 1 To iNy + 1 'Convective Flux d b l V e l D i r = ( d b l V e l ( i V , i , j) + d b l V e l ( i v , i , j - 1 ) ) / 2 I f d b l V e l D i r >= 0 Then 'Flow moves i n forward y d i r e c t i o n ' d b l F l u x ( i , j ) = dbldt / dblDeltaY * ( d b l V e l ( i v , i , j - 1) * db l C o n c ( i , j -1) ) Else 'Flow moves i n negative y d i r e c t i o n d b l F l u x f i , j) = dbldt / dblDeltaY * ( d b l V e K i V , i , j) * d b l C o n c ( i , j ) ) End I f ' D i f f u x i v e Flux d b l F l u x f i , j) = d b l F l u x f i , j) - dbldt * dblDY / (dblDeltaY) * 2 * (dblConc(i, j) - d b l C o n c ( i , j - 1)) Next j Next i End I f End Sub P u b l i c Sub CalculateFluxesSecondUpwind ( b o l V e r t i c a l As Boolean, dblConcO As Double, d b l V e l () As Double,' _ d b l F l u x O As Double, dblDeltaX As Double, dblDeltaY As Double, dblDX As Double, dblDY As Double, dbldt As Double, iNx As Integer, iNy As Integer) 'Subroutine t o c a l c u l a t e d i r e c t i o n a l f l u x e s ' both h o r i z o n t a l (F) and v e r t i c a l (G) code i s included here ' and the switch to choose the appropriate path i s ' b o l V e r t i c a l ' 'Generate Flux I n t e g r a l Matrices f o r n timestep 'The 1/2 increments are rounded up i n determining the i n t e g e r 'for the f l u x index. 1 eg: F 2+1/2,4 => F(3,4) eg: G 3,1-1/2 => F( 3 , l ) 'Fluxes employ centred d i f f e r e n c e convection and 2nd order accurate d i f f u s i o n ' Flux accuracy confirmed on March 8, 1999 3:40pm Dim i As Integer Dim j As Integer Dim i U As Integer Dim i V As Integer Dim d b l V e l D i r As Double i U = 1 i V = 2 I f b o l V e r t i c a l = True Then 'Calculate V e r t i c a l Flux M a t r i x ' NEED TO RUN THIS TO INX+1 AND INY+1 BUT NEED TO SET VELOCITY VALUES OUTSIDE ' THE STANDARD INX,INY MATRIX Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 150 j ) ) For j = 1 To iNy + 1 For i = 1 To iNx + 1 'Convective Flux d b l V e l D i r = (dblVel d u , i , j) + d b l V e K i U , i - 1, j ) ) / 2 I f d b l V e l D i r >= 0 Then 'Flow moves i n forward x d i r e c t i o n I f i = 1 Then ' F i r s t Order d b l F l u x d , j) = dbldt / dblDeltaX * (dblVel (iU, i - 1, j ) * d b l C o n c d -El s e 'Second Order d b l F l u x d , j) = dbldt / dblDeltaX * (3 * d b l V e l d U , i - 1, j) * d b l C o n c d - 1, j ) - d b l V e l d U , i - 2, j ) * d b l C o n c d - 2, j ) ) / 2 End I f El s e 'Flow moves i n negative x d i r e c t i o n I f i = iNx + 1 Then ' F i r s t Order d b l F l u x d , j) = dbldt / dblDeltaX * ( d b l V e l d U , i , j ) * d b l C o n c d , j ) ) El s e 'Second Order d b l F l u x d , j) = dbldt / dblDeltaX * (3 * d b l V e l d U , i , j) * db l C o n c d , j) - d b l V e l d U , i + 1, j ) * d b l C o n c d + 1, j ) ) / 2 End I f End I f ' D i f f u x i v e Flux d b l F l u x d , j) = d b l F l u x d , j) - dbldt * dblDX / (dblDeltaX) * 2 * (dblConcd, j) - d b l C o n c ( i - 1, j ) ) Next i Next j El s e ' C a l c u l a t e H o r i z o n t a l f l u x M a t r i x For i = 1 To iNx + 1 For j = 1 To iNy + 1 'Convective Flux ' second order f l u x except at boundaries then f i r s t order. d b l V e l D i r = (dblVel ( i v , i , j) + dblVel ( i v , i , j - 1)) / 2 I f d b l V e l D i r >= 0 Then 'Flow moves i n forward y d i r e c t i o n I f j = 1 Then ' F i r s t Order d b l F l u x d , j) = dbldt / dblDeltaY * ( d b l V e K i v , i , j - 1) * d b l C o n c d , j - 1) ) E l s e 'second order d b l F l u x d , j) = dbldt / dblDeltaY * (3 * d b l V e l ( i v , i , j - 1) * db l C o n c d , j - 1) - d b l V e K i v , i , j - 2) * db l C o n c d , j - 2)) / 2 End I f E l s e 'Flow moves i n negative y d i r e c t i o n I f j = iNy + 1 Then ' F i r s t order d b l F l u x d , j) = dbldt / dblDeltaY * ( d b l V e K i v , i , j ) * d b l C o n c d , j ) ) El s e 'second order d b l F l u x d , j) = dbldt / dblDeltaY * (3 * d b l V e K i v , i , j) * d b l C o n c d , j) - d b l V e K i v , i , j +1) * d b l C o n c d , j +1)) / 2 Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 151 End I f End I f ' D i f f u x i v e Flux d b l F l u x f i , j) = d b l F l u x f i , j ) - dbldt * dblDY / (dblDeltaY) * 2 * (dblConc(i, j) - d b l C o n c ( i , j - 1 ) ) Next j Next i End I f End Sub P u b l i c Sub CalculateFluxesSecondUpwindTVD ( b o l V e r t i c a l As Boolean, dblConcO As Double, d b l V e l () As Double, _ d b l F l u x O As Double, dblDeltaX As Double, dblDeltaY As Double, dblDX As Double, dblDY As Double, dbldt As Double, iNx As Integer, iNy As Integer) 'Subroutine to c a l c u l a t e d i r e c t i o n a l f l u x e s ' both h o r i z o n t a l (F) and v e r t i c a l (G) code i s included here ' and the s w i t c h t o choose the appropriate path i s ' b o l V e r t i c a l ' 'Generate Flux I n t e g r a l Matrices f o r n timestep 'The 1/2 increments are rounded up i n determining the i n t e g e r 'for the f l u x index. ' eg: F 2+1/2,4 => F(3,4) eg: G 3,1-1/2 => F ( 3 , l ) 'Fluxes employ centred d i f f e r e n c e convection and 2nd order accurate d i f f u s i o n 1 Flux accuracy confirmed on March 8, 199 9 3:40pm Dim j As Integer Dim i U As Integer Dim i V As Integer Dim d b l V e l D i r As Double Dim dblR As Double Dim d b l P h i As Double i U = 1 i v = 2 Dim i As Integer I f b o l V e r t i c a l = True Then 'Calculate V e r t i c a l Flux M a t r i x ' NEED TO RUN THIS TO INX+1 AND INY+1 BUT NEED TO SET VELOCITY VALUES OUTSIDE ' THE STANDARD INX,INY MATRIX For j = 1 To iNy + 1 For i = 1 To iNx + 1 'Convective Flux d b l V e l D i r = ( d b l V e l ( i U , i , j) + d b l V e l ( i U , i - 1, j ) ) / 2 I f d b l V e l D i r >= 0 Then 'Flow moves i n forward x d i r e c t i o n I f i = 1 Then ' F i r s t Order d b l F l u x f i , j) = dbldt / dblDeltaX * ( d b l V e l ( i U , i - 1, j) * dblConc(i -1, j ) ) E l s e 'Second Order ' d b l F l u x f i , j ) = dblDT / dblDeltaX * (3 * d b l V e l ( i U , i - 1, j) * d b l C o n c ( i - 1, j ) - d b l V e l ( i U , i - 2, j ) * d b l C o n c ( i - 2, j ) ) / 2 I f ( d b l V e K i U , i - 1, j) * d b l C o n c f i - 1, j) + d b l V e K i U , i - 2, j) * d b l C o n c ( i - 2, j ) ) <> 0 Then dblR = ( ( d b l V e K i U , i , j) * dblConc ( i , j ) ) + (dblVel (iU, i - 1, j) * db l C o n c ( i - 1, j ) ) ) ' Flux R a t i o f o r TVD Scheme 'Flux C o r r e c t i o n f a c t o r f o r TVD Scheme Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 152 - 2 , j) * d b l C o n c d - 2, j) ) ) El s e / ( ( d b l V e l d U , i - 1, j) * d b l C o n c d - 1, j ) ) + ( d b l V e l d U , i dblR = 2 End I f 'Van Leer Scheme d b l P h i = (dblR + Abs(dblR)) / (1 + dblR) L i m i t e d E x t r a p o l a t i o n I f dblR <= 0 Then d b l P h i = 0 E l s e l f dblR <= 0.5 Then d b l P h i = 2 * dblR E l s e d b l P h i = 1 End I f d b l F l u x d , j) = dbldt / dblDeltaX * ( ( d b l V e l d U , i - 1, j ) * d b l C o n c d -1, j ) ) _ + d b l P h i / 2 * ( d b l V e l d U , i - 1, j) * d b l C o n c d - 1, j ) - d b l V e l d U , i - 2, j) * d b l C o n c d - 2, j ) ) ) End I f El s e 'Flow moves i n negative x d i r e c t i o n I f i = iNx + 1 Then ' F i r s t Order d b l F l u x d , j ) = dbldt / dblDeltaX * ( d b l V e l d U , i , j) * d b l C o n c d , j ) ) El s e 'Second Order ' d b l F l u x d , j ) = dblDT / dblDeltaX * (3 * d b l V e l d U , i , j ) * d b l C o n c d , j) - d b l V e l d U , i + 1, j) * d b l C o n c d + 1, j ) ) / 2 I f ( d b l V e l d U , i , j ) * d b l C o n c d , j ) + d b l V e l d U , i + 1, j ) * d b l C o n c d + 1, j ) ) <> 0 Then dblR = ( ( d b l V e l d U , i - 1, j) * d b l C o n c d - 1, j ) ) + ( d b l V e l d U , i , j) * d b l C o n c d , j) ) ) / ( ( d b l V e l d U , i , j) * d b l C o n c d , j ) ) + ( d b l V e l d U , i + 1, j) * d b l C o n c d + 1, j) ) ) El s e dblR = 1 End I f 'Van Leer Scheme d b l P h i = (dblR + Abs(dblR)) / (1 + dblR) L i m i t e d E x t r a p o l a t i o n I f dblR <= 0 Then d b l P h i = 0 E l s e l f dblR <= 0.5 Then d b l P h i = 2 * dblR E l s e d b l P h i = 1 End I f j) * db l C o n c ( i + 1 , j ) ) ) d b l F l u x d , j) = dbldt / dblDeltaX * ( ( d b l V e l d U , i , j) * d b l C o n c d , j ) ) + d b l P h i / 2 * ( d b l V e l d U , i , j) * d b l C o n c d , j ) - d b l V e l d U , i + 1," End I f End I f ' D i f f u x i v e Flux d b l F l u x d , j) = d b l F l u x d , j) - dbldt * dblDX / (dblDeltaX) * 2 * (dblConcd j) dblConc d - 1, j ) ) Next i Next j Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 153 E l s e 'Calculate H o r i z o n t a l f l u x M a t r i x 1) ) For 1 = 1 To iNx + 1 For j = 1 To iNy + 1 'Convective Flux 1 second order f l u x except at boundaries then f i r s t order. d b l V e l D i r = ( d b l V e l ( i v , i , j) + d b l V e l ( i V , i , j - 1)) / 2 I f d b l V e l D i r >= 0 Then 'Flow moves i n forward y d i r e c t i o n I f j = 1 Then ' F i r s t Order d b l F l u x ( i , j) = dbldt / dblDeltaY * ( d b l V e l ( i V , i , j - 1) * dbl C o n c ( i , j db l C o n c ( i , j - 1) E l s e 'second order ' d b l F l u x f i , j) = dblDT / dblDeltaY * (3 * d b l V e l ( i v , i , j d b l V e l ( i V , i , j - 2) * dbl C o n c ( i , j - 2 ) ) / 2 1) db l C o n c ( i , j db l C o n c ( i , j i . j I f ( d b l V e l ( i V , i , j - 1) * dbl C o n c ( i , j - 1) + d b l V e l ( i V , i , j - 2) * 2)) <> 0 Then dblR = ( ( d b l V e l ( i V , i , j) * db l C o n c ( i , j ) ) + ( d b l V e l ( i V , i , j - 1) - in: 2) * d b l C o n c ( i , j - 2))) El s e / ( (dblVel (iV, i , j - 1) * dblConc ( i , j - 1) ) + (dblVel (iV, dblR = 2 End I f 'Van Leer Scheme d b l P h i = (dblR + Abs(dblR)) / (1 + dblR) L i m i t e d E x t r a p o l a t i o n I f dblR <= 0 Then d b l P h i = 0 E l s e l f dblR <= 0.5 Then d b l P h i = 2 * dblR E l s e d b l P h i = 1 End I f 1) ) d b l F l u x f i , j) = dbldt / dblDeltaX * ( ( d b l V e l ( i V , i , j - 1) * dbl C o n c ( i , j + d b l P h i / 2 * ( d b l V e l ( i V , i , j - 1) * db l C o n c ( i , j - 1) - d b l V e l ( i V , 2) ) ) i , j - 2) * d b l C o n c ( i , j End I f El s e 'Flow moves i n negative y d i r e c t i o n I f j = iNy + 1 Then ' F i r s t order d b l F l u x f i , j) = dbldt / dblDeltaY * ( d b l V e l ( i V , i , j) * db l C o n c ( i , j ) ) El s e 'second order ' d b l F l u x ( i , j) = dblDT / dblDeltaY * (3 * d b l V e l ( i V , i , j) * db l C o n c ( i , j ) - d b l V e l ( i V , i , j + 1) * db l C o n c ( i , j +1)) / 2 j + 1)) <> 0 Then j) * d b l C o n c ( i , j ) ) ) * d b l C o n c ( i , j + 1))) I f ( d b l V e K i V , i , j) * db l C o n c ( i , j) + d b l V e K i V , i , j + 1) * dbl C o n c ( i , dblR = ( ( d b l V e K i V , i , j - 1) * dblConc ( i , j - 1 ) ) + ( d b l V e K i V , i , / ( ( d b l V e K i V , i , j) * dblConc ( i , j ) ) + ( d b l V e K i V , i , j + 1) El s e dblR = 1 End I f 'Van Leer Scheme d b l P h i = (dblR + Abs(dblR)) / (1 + dblR) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 154 'Limited E x t r a p o l a t i o n ' I f dblR <= 0 Then d b l P h i = 0 ' E l s e l f dblR <= 0.5 Then 1 d b l P h i = 2 * dblR ' E l s e 1 d b l P h i = 1 'End I f d b l F l u x d , j) = dbldt / dblDeltaX * ( (dblVel (iv, i , j) * d b l C o n c d , j ) ) + d b l P h i / 2 * (dblVel (iv, i , j) * d b l C o n c d , j) - d b l V e l (iv, i , j +' 1) * d b l C o n c ( i , j + 1))) End I f End I f ' D i f f u x i v e Flux d b l F l u x d , j) = d b l F l u x d , j) - dbldt * dblDY / (dblDeltaY) A 2 * (dblConcd, j) - db l C o n c ( i , j - 1)) Next j Next i End I f End Sub Sub ExportFlowDetailsCheck(iNx As Integer, iNy As Integer, _ dblDX As Double, dblDY As Double, strComments As S t r i n g , _ dblVx As Double, dblVy As Double, numiter As Long, _ i n t E x p o r t I n t e r v a l As Integer, dbldt As Double, dblMaxX As Double, _ dblMaxY As Double, i n t l t a r g e t As Integer, i n t J t a r g e t As Integer, dblMass As Double) Open "FlowDetails.dat" For Output As 1 P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, P r i n t #1, Close 1 "Date: " & vbTab & Now() "iNX: " & vbTab & iNx " iNY:" & vbTab & iNy "MaxX:" & vbTab & dblMaxX "MaxY:" & vbTab & dblMaxY "DispX:" St vbTab & dblDX "DispY:" & vbTab & dblDY "VelX:" & vbTab & dblVx "VelY:" & vbTab & dblVy " T a r g e t l : " & vbTab £. i n t l t a r g e t "TargetJ:" & vbTab & i n t J t a r g e t "Mass:" & vbTab & dblMass "Timestep:" & vbTab & dbldt " I t e r a t i o n s : " & vbTab & numiter "Output I n t e r v a l : " & vbTab & i n t E x p o r t l n t e r v a l "Comments:" & vbTab & strComments End Sub Sub MassInTank (strFilename As S t r i n g , iNx As Integer, iNy As Integer, dblConcd As Double, dblTime As Double) Dim dblMass As Double Dim i As Integer Dim j As Integer dblMass = 0 For i = 1 To iNx For j = 1 To iNy dblMass = dblMass + d b l C o n c d , j) Analysis of Laboratory Scale Aerated Basin Using Computational Techniques Page 155 Next j Next i Open strFilename For Append As 1 Write #1, dblTime, dblMass Close 1 End Sub Sub ConcAtOutlet(strFilename As S t r i n g , iNx As Integer, iNy As Integer, dblConc() As Double, dblTime As Double, intOutBottom As Integer, intOutTop As Integer) Dim dblConcSum As Double Dim i As Integer Dim j As Integer Dim intCVCount As Integer dblConcSum = 0 intCVCount = intOutTop - intOutBottom + 1 i = iNx For j = intOutBottom To intOutTop dblConcSum = dblConcSum + d b l C o n c d , j) Next j dblConcSum = dblConcSum / intCVCount Open strFilename For Append As 1 Write #1, dblTime, dblConcSum Close 1 End Sub
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis of laboratory scale aerated basin using computational...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis of laboratory scale aerated basin using computational techniques Jenkinson, Robert Wayne 2000
pdf
Page Metadata
Item Metadata
Title | Analysis of laboratory scale aerated basin using computational techniques |
Creator |
Jenkinson, Robert Wayne |
Date Issued | 2000 |
Description | This study endeavored to apply computational fluid dynamics (CFD) techniques to aerated stabilization basins (ASBs) on a laboratory scale. A laboratory study was set up to assess the hydraulic impacts of a small-scale mechanical surface aerator in a still-water tank. A small surface mechanical aerator was built that would circulate from 0.5 to 4 L/sec. The aerator was employed in a tank and the velocity patterns imported by the aerator were measured using an acoustic Doppler velocimeter (ADV). The velocity profiles were converted into dispersion parameters for use in a CFD model. A local turbulent time scale method was used to determine an approximate dispersion value as a function of distance from the aerator. The analysis examined different aerator power settings under otherwise identical conditions. Three successful runs were completed. The chaotic patterns within the tank induced by the aerator, produced turbulent structures that were not easily resolved using a time averaged statistical turbulence approach. Tracer studies were performed in the same basin, modified to provide a flow-through environment. The studies were performed with and without aeration and at different flowthrough rates. Tracer study results were examined and some of the typical parameters used in the literature were extracted for comparative analysis. The trends were surprising as the dimensionless axial dispersion number was inversely related to the degree of mixing applied. A two-dimensional computational fluid dynamics (CFD) model was employed to model the fluid flow within the basin without aeration using a commercial package, FLUENT. The modeling results were compared with depth-averaged ADV velocity measurements taken within the basin. Many CFD modeling parameters were varied, including grid size and structure, turbulence modeling and boundary condition modeling. The model compared well with the ADV data. The results from the FLUENT modeling were applied to a second CFD model that solved the solute transport equations over a structured grid. The solute transport modeling only examined the transport equations with uniform dispersion over the basin and, although it was the original intent, no attempts were made to model aeration at this stage. The overall shape of the modeled tracer curves compared favourably to the tracer study data obtained in the lab. |
Extent | 7935391 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0063992 |
URI | http://hdl.handle.net/2429/10260 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-0090.pdf [ 7.57MB ]
- Metadata
- JSON: 831-1.0063992.json
- JSON-LD: 831-1.0063992-ld.json
- RDF/XML (Pretty): 831-1.0063992-rdf.xml
- RDF/JSON: 831-1.0063992-rdf.json
- Turtle: 831-1.0063992-turtle.txt
- N-Triples: 831-1.0063992-rdf-ntriples.txt
- Original Record: 831-1.0063992-source.json
- Full Text
- 831-1.0063992-fulltext.txt
- Citation
- 831-1.0063992.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0063992/manifest