RANDOM WALK MODEL APPLICABLE TO RIVERS By Jason Vine B.A.Sc. (Civil Engineering) University of British Columbia 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 CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A August 1996 © Jason Vine, 1996 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 CC,v, I F^,n,er,ry The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract The disposal of waste effluent from industrial and municipal sources into rivers for the purpose of dilution is a common practice. Dilution occurs in the near-field due to specifics of the diffuser design, and in the mid-field due to parameters of the ambient environment. The ability to predict downstream dilution is of interest to regulatory officials and engineers who wish to design such disposal systems. A random walk model has been developed which is geared towards use in rivers. The model only performs mid-field calculations, however the potential exists to accept the results from a near-field model as input. The model attempts to account for the major velocity structures that occur in an open channel, including transverse and vertical shear, secondary currents due to flow through a bend, and turbulent fluctuations. Equations describing such phenomena are obtained from the literature. Laboratory experiments are conducted to verify the validity of some of these equations. As expected from theory, the research confirms that random walk modelling predicts the same results as Fickian diffusion. Arguments are presented as to the choice of a suitable timestep for use with a random walk model. Several key parameters of the random walk are identified, and their effect on the stability and accuracy of the model results are examined. The model is applied to the disposal of effluent from a pulp and paper mill located in Prince George, B.C. Results agree acceptably with a popular diffusion model as well as with digital imagery of the plume. Plots of plume width and dilution are produced for various river flow rates which may be expected at Prince George. i i i Table of Contents Abstract ii List of Tables v List of Figures vi List of Symbols vii Acknowledgments ix 1. Introduction 1 1.1 Disposal of Effluent 1 1.2 Random Walk Modell ing 2 1.3 Objectives and Scope of Research 3 1.4 Application Site 3 2. Literature Review 4 2.1 Frames of Reference 4 2.2 Advection-Diffusion Equations 5 2.3 The Random Walk 8 2.4 Lagrangian Particle Tracking Models 9 2.5 Shear Velocity 10 2.6 Transverse M i x i n g Coefficient 11 2.7 Vertical M i x i n g Coefficient 12 2.8 M i x i n g Length Scales 13 2.9 Mean Longitudinal Velocity Distribution 13 2.10 Secondary Currents 16 2.11 Turbulence Intensity 20 3. Experiments 28 3.1 Purpose 28 3.2 Equipment 28 iv 3.3 Experiments 29 4. Model Development - PATRAC 42 4.1 General Specifications 42 4.2 Geometry 43 4.3 Flowlines 44 4.4 Mean Longitudinal Velocity Distribution 44 4.5 Mean Vertical And Transverse Velocity Distributions 44 4.6 Shear Velocity 45 4.7 Secondary Currents 45 4.8 Turbulent Fluctuations 46 4.9 Transverse Calibration 47 4.10 Vertical Calibration 47 4.11 Model Output 48 4.12 Choice of Suitable Random Walk Parameters 48 5. Application to the Fraser River 53 5.1 Validation of the Northwood Site Simulation 53 5.1.1 Comparison of P A T R A C With Cormix 54 5.1.2 Comparison of P A T R A C with Field Data 55 5.2 Extreme Flow Case Simulations 57 5.3 Calculation of Downstream Dilution 58 Chapter 6 Conclusions and Recommendations 67 6.1 Summary and Conclusions 67 6.2 Recommendations 69 6.3 Further Research 69 References 71 APPENDIX A - PATRAC USERS MANUAL 73 APPENDIX B - PATRAC LISTINGS 85 List of Tables 2.1 Length Scales for Complete Transverse Mixing 13 2.2 Open Channel Flow Regions 14 2.3 Comparison of Longitudinal Velocity Profiles 16 5.1 Summary of Simulations 54 A . l P A T R A C Software Components 73 List of Figures 2.1 Definition sketch of channel geometry and velocity profiles 23 2.2 Various Frames of reference 24 2.3 Open channel flow regions 25 2.4 Comparison of velocity laws 26 2.5 Plot showing how secondary current surface magnitude normalized by the mean velocity varies with the ratio of depth to radius of curvature 27 3.1 Time series of velocity fluctuations taken at 25 Hz for deep experiment 33 3.2 Time series of velocity fluctuations taken at 25 Hz for shallow experiment 34 3.3 Power spectral density plots for deep experiment 35 3.4 Power spectral density plots for shallow experiment 36 3.5 Filtered low frequency power spectral density plots for deep experiment 37 3.6 Longitudinal velocity distributions showing predicted low velocity law and experimental values 38 3.7 Reynolds stress distributions 39 3.8 Standard deviation of velocity signals for deep experiment 40 3.9 Standard deviation of velocity signals for shallow experiment 41 4.1 Flowchart of the random walk engine of P A T R A C 50 4.2 Transverse calibration of a single plume 51 4.3 Sensitivity/stability analysis of random walk simulation 52 5.1 Digital image of the Northwood site on the Fraser River 60 5.2 Comparison of P A T R A C and CORMIX using a straight, rectangular channel 61 5.3 Plume width as predicted by the Northwood site simulation with P A T R A C 62 5.4 Overhead snapshot of the Northwood site at 805 m^/s flow 63 5.5 Variation of parameters affecting transverse mixing given various flow rates and distance downstream 64 5.6 Contour maps of dilution for various river flow rates 65 5.7 Plume centreline and bank dilution plot 66 vii List of Symbols a dimensionless transverse mixing number = £t/du* as acceleration along the s direction A constant of integration b channel width C concentration VC spatial gradient of concentration Q generic constant d depth of water D ,D ,D longitudinal, transverse, and vertical molecular diffusion coefficients transverse, vertical turbulent diffusion coefficients 8 gravitational acceleration H0 total energy at a point K von Karman constant m friction factor M mass of pollutant M mass/unit time P pressure Pobs' Psim observed, simulated width of waste plume Py probability of particle being located between y and y+Ay P density of water q solute mass flux R radius of curvature of channel R. Reynolds number = du,* 1 v K hydraulic radius s channel slope t time At random walk time step bed shear stress arc angle at which secondary current circulation is developed u advective velocity vector u,v,w longitudinal, transverse, and vertical velocities u',v',w' longitudinal, transverse, and vertical turbulence intensities u,v,w longitudinal, transverse, and vertical mean velocities random normal deviates viii u max free stream velocity U+ dimensionless velocity = u/u* U max dimensionless free stream velocity (surface velocity) W» shear velocity Vs transverse secondary current surface velocity Vsc>Wsc transverse and vertical secondary current velocities V kinematic viscosity wf wake function Wfall fall velocity of a particle ™s maximum vertical secondary current velocity x,y,z longitudinal, transverse, and vertical directions Ax,Ay,Az longitudinal, transverse, and vertical random walk distances z+ dimensionless depth = zw, / v c standard normal deviate Acknowledgments The author would like to acknowledge his supervisor and friend Dr. G. A . Lawrence, without whom this thesis would not be possible. Thanks to the Fraser River Action Plan which provided funding for this research. The digital image of the Fraser River was acquired and processed by Borstad Associates Ltd. of Sidney, British Columbia with partial funding of Department of Environment, Department of Fisheries and Oceans and the Environmental Innovations Program (EIP). Many thanks to the shop guys (especially Dick Postgate and Kurt Nielson) for helping me with so many lab problems. My thanks to the Civi l Engineering faculty (especially Dennis Russell, Michael Quick, and Bi l l Caselton) at UBC, who provided excellent guidance and support throughout my seven year academic career at UBC. More thanks to the office personnel for all of their help with my many requests. Finally, thanks to all my friends that I made at U B C , and especially to Bonnie Marks, who is a great friend and the finest coffee buddy I could have found. Chapter 1 Introduction 1.1 Disposal of Eff luent The disposal of liquid waste from industrial and municipal sources into receiving waters for the purpose of dilution is a common practice. Indeed, many large scale industrial processes such as pulp and paper milling are located near a body of water to provide both a source of water as well as a method of disposal. The effluent which is disposed of may contain many different substances including heavy metals, organic and inorganic chemicals, and suspended solids. Regulation of the permissible concentrations of these substances in the receiving waters is undertaken by governmental environment ministries. In British Columbia, these concentrations are outlined by the British Columbia Ministry of Environment (1994). Acceptable concentrations are usually based on acute and chronic effects to wildlife and humans. Given that maximum allowable concentrations of pollutant are specified by law, it is up to the discharger to ensure that acceptable dilution of the waste plume is achieved. The dilution that is obtained is a function of both the method of disposal and the ambient environment into which the effluent is disposed. The engineer has control over the disposal mechanism, but not over the ambient environment. Disposal of the effluent in the near-field is achieved through a diffuser, which can be of either single or multi-port design. The sizing, angles and number of ports are some of the variables which may be modified in order to maximize the amount of near-field dilution that is obtained. Once the effect of near-field dilution has subsided, the remaining dilution that occurs is due to factors in the environment. Here the effects of ambient velocity structure and turbulence dilute and transport the waste plume. In a river, the velocity profiles can become quite complex, and may include boundary effects, rapid changes in channel geometry, and additional velocity structures such as circulation due to bends in the channel. It has also been suggested that the interaction of effluent with sediment particles may be of concern (Evans 1996). The ability to predict the dilution achieved by both near and far-field processes is of great importance. Government employees may wish to study the effect of the addition of a new pollution source. Engineers want to ascertain if a particular design for an effluent discharge will achieve the desired results. Company or municipal employees, anxious to avoid costly fines associated with violation of disposal laws, may wish to predict what downstream effects will occur due to changes in an industrial process. Various methods, equations and programs exist to predict the dilution of a plume (see Fischer et al. (1979) and Wood et al. (1993)). All are simplifications of an extremely complex problem which may never be completely solved due to the stochastic nature of turbulence. 1.2 Random Walk Modelling The random walk or "drunkards" walk, which is described in Fischer et al. (1979), has been used for years to aid instructors in explaining the concept of Fickian diffusion. Practical application of the method to pollution dispersion has only recently become practical due to the large computational effort required. The development of a pollutant transport model using this approach is an interesting exercise and provides a useful tool from which further work on the interaction between sediment and pollutant may be undertaken. Interest in this technique should only increase in coming years with more powerful computing equipment being made available. 1.3 Objectives and Scope of Research This thesis describes the development of a particle tracking computer model which is applicable to modelling effluent dilution in rivers. The objectives of this research are to: • identify the major mechanisms responsible for transport and dispersion of pollutant in a river channel. • conduct research into random walk modelling, which is potentially very powerful, but only now becoming practical due to computational requirements. • develop a dilution prediction model which is capable of being used to simulate different sites. • apply the model to a specific site in British Columbia. 1.4 Application Site Funding for this thesis was provided by the Fraser River Action Plan, which has a specific interest in pulp and paper mill effluent disposal into the Fraser River. The computer model which was developed is specifically geared towards the Northwood site located in Prince George, B.C., however the model is capable of being reprogrammed to accommodate a different location. The reader may wish to consult Environment Canada (1994) for further information on research/reports conducted in the Fraser River basin. Chapter 2 Literature Review This chapter concentrates on the background information used to produce a particle tracking model suitable for use in rivers. Some previous work on particle tracking models is also reviewed. To aid the reader, figure 2.1 shows a sketch of the coordinate system and important variables. The first section aims to help the reader understand a fundamental difference between a random walk model and one which is based on Fickian diffusion. Following this, sections are presented on advective-diffusion equations, random walk, and Lagrangian particle tracking. The remainder of the chapter is devoted to describing the various mechanisms which are responsible for turbulent diffusion in a river. 2.1 Frames of Reference A major difference between Fickian diffusion and random walk modelling is the frame of reference that each uses. Fickian diffusion describes concentration in an Eulerian frame, whereas random walk models use the Lagrangian frame of reference. 1) Eulerian - The Eulerian frame of reference does not move, and specifies quantities at all points in space at a given time, providing a "snapshot". The Eulerian frame may be used to describe unsteady conditions, such as the motion of a cloud of pollutant as shown in figure 2.2(a). If steady-state conditions exist, such as the case of a diffuser which has been operating continuously for some time, the time variable (t) may be replaced with downstream distance divided by velocity (x/u). The result is an equation for concentration which is only a function of location and mean velocity. The use of the Eulerian frame sometimes makes the determination of any given quantity (such as concentration of a pollutant) as simple as using a formula such as described in the next section. These analytic solutions give the concentration as a simple function of location and various other parameters. 2) Lagrangian - Contrast Fickian diffusion and the Eulerian frame with random walk modelling, which uses the Lagrangian frame. The Lagrangian frame shown in figure 2.2(b), specifies the variation of a particle or packet of fluid as a function of time. Rather than specifying a quantity of interest (concentration) as a function of space and time, random walk modelling and the Lagrangian frame specify a particle position as a function of time. The output from a Lagrangian particle tracking model does not usually translate into a quantity of interest. Additional analysis must be performed to convert these particle "tracks" into scalar quantities such as concentration. This is described in more detail later. 3) Instrumental - There is one additional frame of reference which is valid when spatially fixed sensors are used in experimental setups. This Instrumental frame, shown in figure 2.2(c), describes the variation of a quantity at a fixed point in space as a function of time. Ideally, experimental sensors would be able to determine the quantity of interest at all locations and at all times. This is not possible with typical experimental sensors, however there is a way to relate the Instrumental frame to Eulerian-space. Taylor's hypothesis suggests that if a turbulent field sweeps past the sensor much faster than it changes, then the Instrumental frame of reference is equivalent to the Eulerian-space frame. The acoustic Doppler velocimeter, described in chapter 3, uses this frame of reference. 2.2 Advection-Diffusion Equations The traditional approach for determining contaminant dispersion in a river is to use advection-diffusion equations developed from Fickian diffusion theory. There are several general works on this process, which include Fischer et al. (1979), Csanady 6 (1972), and Lam et al. (1984). A l l of these references include sections on Fick's law, which can be stated mathematically as : where: q = solute mass flux D = diffusion coefficient V C = spacial gradient of concentration In rivers, it is desirable to add an advective term to this equation resulting in : where: C= concentration U = three dimensional advective velocity vector Fischer (1979) uses Green's theorem to write conservation of mass as : dC ~dt- = - V q Equations 2.2 and 2.3 may be combined as follows : q=-DVC (2.1) q = UC + (-DVC) (2.2) (2.4) If U is assumed to be constant, than equation 2.4 becomes : dt (2.5) Equation 2.5 may be written out in Cartesian coordinates as dC dC dC dC — + U— + V— + W — dt dx dy dz dx dx dy\ • dy _ + — + • dz\ z dz (2.6) where: Dx,Dy,Dz = diffusion coefficients in the three principal directions u,v,w = longitudinal, transverse, and vertical mean velocities Solution of this equation now depends on any simplifying assumptions the modeller makes and the specifics of the application. Solutions can be accomplished by developing analytical equations or by numerical methods. The above mentioned references provide various analytical solutions to the advective-diffusion equation. One of the most fundamental is the solution for an instantaneous point source in three dimensions, which is given by Fischer et al. (1979) and Lam et al. (1984) as : C(x,y,z,t) = - M (4ntf2(DxDxDz) I 7 Te x P -1 4t (2.7) where: M = Mass of pollutant Of more direct use for modelling effluent discharge is the steady state solution combined with the assumption that mixing occurs very quickly over the depth : C(x,y) = - M ud-^AnDx I u exp In. 4Dx (2.8) where: M = mass / unit time which both Fischer et al. (1979) and Csanady (1972) provide as the two dimensional steady state solution to the advective-diffusion equation. 2.3 The Random Walk Fischer et al. (1979) provides a detailed description of the random walk process, however a quick summary of the major points is helpful. We begin by defining what a tracer particle is for the purposes of random walk modelling. It is desirable that a tracer particle follow the flow pattern exactly, without the effects of gravity, buoyancy or inertia. These conditions are difficult to obtain in the laboratory, but simple to impart in a computer simulation. Fischer et al. (1979) describes a particle as "...any sufficiently small, but identifiable, neutrally buoyant entity whose size is such that its dynamical behavior is essentially indistinguishable from that of the fluid". To define the random walk, begin by assuming that a tracer particle is allowed to follow a series of random steps of size Ay and time At, with an equal probability of moving forward or backward. On average the particle will not move, but if a series of particles are allowed to follow such motion, a concentration distribution begins to develop. Fischer et al. (1979) describes the use of the central limit theorem to develop a relationship for this concentration distribution. The following describes the probability that a particle will be between the points y and y+dy : py(y,t)dy = f V 44KD1 ) \4Dt dy (2.9) Fischer et al. (1979) finishes the description by stating that if a large group of particles are released at time zero at the origin, the concentration distribution of particles can be described by: 9 C(y,t) = •sjAnDt M exp - (2.10) This result is the same as the one dimensional analytic solution to the diffusion equation, indicating that the random walk process produces the same prediction as Fickian diffusion. Another useful equation relates the dispersion coefficient to the values of Ay and At used in the random walk : The above equation will be used in chapter 4 to aid in the choice of a suitable time step for the particle tracking models described in the next section. 2.4 Lagrangian Particle Tracking Models Rather than trying to numerically solve an advection-diffusion equation or use an analytic solution to one, a Lagrangian particle tracking model follows as many as tens of thousands of particles through a simulated flow field. The number of particles used is typically a function of the accuracy required and available computing power. Allen (1991) provides a comparison of simulation run-time on various computing platforms. Simulations using as many as 100,000 particles were possible using parallel processor systems. The path of any one particle provides little information, however the average of a large number of particles produces useful results. The term "particle" can have different meanings. If the particle is a tracer (which is affected solely by the flow field in the model), it is possible to model concentrations of pollutant released from a source such as a diffuser. A particle, if it is not considered a tracer, can also represent a unit of sediment or a bio-solids floe if assigned suitable dynamic properties such as fall velocity and drag coefficient. These particles are acted upon not only by the surrounding flow field, but also by gravitational forces. Any combination of particles can be introduced in the model, allowing for simulation of multiple sediment sizes as well as tracer particles. Heslop & Allen (1993) describe a two-dimensional random-walk model used to simulate a section of the River Severn. The model was initially run with laboratory-derived turbulence parameters, which were later replaced with data obtained from field measurements. The assumption was made that particles which cross over a boundary are reflected back into the water column an equal distance. Results obtained were encouraging, although the authors acknowledge that the model did not perform well in sections of the river which had large, tight meanders. The paper suggests that this was due to the model being two-dimensional, and that a three-dimensional model might solve some of the discrepancies. This hints to the importance of secondary currents, which are described in more detail in section 2.10. 2.5 Shear Velocity The shear stress ( T o ) on the bottom of an open channel flowing with water is often defined for convenience in terms of shear velocity (u t ) . Shear velocity is an important quantity which must be determined when modelling open channel flow. Among other uses, the equations presented in the following sections describing velocity distribution, turbulence intensity, and secondary currents all do so in terms of the shear velocity. For river work, there are three methods for determining this value. Nezu & Nakagawa (1993) describe these methods, which are summarized below. 1) Use of the following relationship allows determination of the shear velocity when actual flow measurements are not available: where: p The wide channel assumption is made that y~Rh, resulting in: (2.13) 2) If velocity measurements are available in the wall region (z/d<0.2), the shear velocity can be determined from the log-law: Assuming a value of K=0AI, u* can be determined so as to best fit the measured data. 3) If measurements of Reynolds stress are available, the shear velocity can be determined from the theoretical Reynolds stress distribution: 2.6 Transverse Mixing Coefficient The turbulent transverse mixing coefficient, et, is often substituted in place of the molecular diffusion coefficient D to simulate dispersion due to turbulent fluid motion. When used with such equations, the value of et specifies the rate of transverse mixing. On the other hand, the value of et can be inferred from the output of a particle (2.14) dz \ d) (2.15) tracking model and be used to confirm that the model is producing reasonable results. Several textbooks on the matter including Fischer et al. (1979) and Rutherford (1994) report this value as the dimensionless group et/dw*. For convenience, the following definition is made : (2.16) where: a = dimensionless transverse mixing number The reported values of this group vary somewhat, but both Fischer et al. (1979) and Rutherford (1994) present similar results. For an experimental flume, the reported value of a = 0.13 - 0.15. In slow meandering streams, the value increases so that a = 0.6 ± 50%. The increased value of a in real rivers is due to effects of meandering and sidewall irregularities. Section 2.10 describes how this value can be even higher in channel bends. 2.7 Vertical Mixing Coefficient Similar to the turbulent transverse mixing coefficient, the turbulent vertical mixing coefficient, e v , can also be defined. Both Fischer et al. (1979) and Rutherford (1994) provide an approximate value of £ v/dw, =0.067. This simplification comes from the following equation which describes the vertical mixing coefficient for momentum : (2.17) £„ = Kdut\ — .d Fischer et al. (1979) references an experimental study by Jobson and Sayre (1970) which confirms the validity of this equation for describing the vertical mixing of mass as well. When equation 2.17 is averaged over the depth, the value of ev/dw* =0.067 is obtained. 2.8 Mixing Length Scales Formulae that predict the distance downstream at which complete mixing occurs are useful, although the definition of complete mixing varies depending on the literature. Fischer et al. (1979) bases this on the concentration being within 5% of the mean value everywhere on the cross section, while Rutherford (1994) uses 2% of the mean value. Table 2.1 summarizes the simplified formulae. Table 2.1 - Length Scales for Complete Transverse Mixing Case Fischer (1979) (95% of Mean Value) Rutherford (1994) (98% of Mean Value) Transverse/Bank Source x = 0.4ub21 et x = Q.536ub21 e, Transverse/Centreline x = 0. lub21 e, x = QA34ub2 /£, The above table is a summary of constants derived for use with the following equation: x = q— (2.18) More information on the derivation of values for C, may be obtained from Fischer et al. (1979). 2.9 Mean Longitudinal Velocity Distribution The distribution of longitudinal velocity with depth is described in many textbooks including Henderson (1966), Dougherty et al. (1985), Fischer et al. (1979), and Rutherford (1994). A l l of these describe Prandtl's mixing length theory and/or the derivation of the well-known velocity distributions. Open channel flow consists of several different regions. Nezu (1993) provides a summary of the equations that describe each region of the flow. Table 2.2 provides a summary of these regions and the various formulae describing each. Table 2.2 - Open Channel Flow Regions Region Criterion Formula Viscous Sublayer z+ « 26 Buffer Layer 5 < z+ < 30 u+ =-3.05 + 5.00 ln(z +) Wall Region 26 < y+ < 0.2 R* z/d < 0.2 u+ = -ln(z+) + A K v ' Outer Region z/d > 0.2 u+ = — ln(z +) + A + wf In the above table, the dimensionless groups are + zut + u dut z =—, u =—, Rt = V U* V These quantities are the dimensionless depth, dimensionless velocity, and Reynolds number based on shear velocity. A typical velocity profile is shown in figure 2.3(a). The formulae given in Table 2.2 are plotted in figure 2.3 in dimensionless coordinates and on log scales, emphasizing the viscous sublayer and buffer layer. In practice, these layers are very small (< 1 cm) and are generally neglected. Traditionally, the formula associated with the wall region (known as the log-law) is applied to the entire depth of flow as a result of recommendations made by Keulegan (1938) and cited in Coleman (1981). The formula can be best-fit to the data for a particular river by adjusting the von Karman constant (K) and the constant of integration (A). A simpler approach is to use the standard values of 0.4 and 5 respectively for pipe flow. Recent work has shown, however, that the log-law may only be valid in the wall region, and that the law may be adjusted by the addition of a "wake" function, wf, and not by adjusting K and A. The resulting law is known in the literature as the "log-wake law" or "velocity defect law". Nezu & Rodi (1986) have performed detailed measurements using a laser Doppler anemometer confirming the modification of flow in the outer region. As a result, the von Karman constant (K) and the constant of integration (A) are truly "constant", with accepted values of K = 0.41 and A=5.29. One of the most widely accepted wake functions is that of Coles (1956), given as : The term EI, known as Cole's wake strength parameter, was found by Nezu and Rodi (1986) to be a function of the Reynolds number R,. For values of i^>2000, found in typical large rivers, the value of II has a constant value of 0.2. Combining the wake function with the log-wake law produces the following velocity defect law: £ / + m a x is specified as the free stream velocity, which is assumed to be the surface velocity. The above relations are more useful when specified in terms of the mean velocity rather than the surface velocity. The log law described in Table 2.2 gives the following relation, which is then substituted back into the above defect laws. (2.20) 'max = U H K (2.21) Table 2.3 gives a summary of the final velocity profile laws specified using the mean velocity. Table 2.3 - Comparison of Longitudinal Velocity Profiles Log Law _ ut U = U H K In Vd) Velocity Defect Law (Coles) _ w* U = U H K In f - l +1 - 2TICOS2 K2dJ_ The two equations are plotted in figure 2.4 using typical values for the Fraser River (if = 1.25 m/s, u*=0.1 m/s, n=0.2). 2.10 Secondary Currents Secondary currents in river bends result in the flow taking on a helical pattern. Secondary currents are discussed in great detail in Rozovskii (1961), while a good summary is available in Henderson (1966). Assume the total energy at any point in the river is given by: d + ^- = H0 (2.22) 2g Taking the derivative: dd udu . — + - — = 0 (2.23) dn g dn In the above equation, n is the direction outwards from the centre of the curve. In addition to the total energy, Henderson (1966) describes the Euler equation which is essentially a balance between the forces exerted on a unit volume of fluid. These forces are caused by a pressure gradient, as well as the weight of the fluid, resulting in the following which is the Euler equation: where: —{p+Yz) + pas=0 (2.24) p = pressure as = acceleration along the s direction If the pressure distribution is assumed to be hydrostatic, then p + yz = yd and the Euler equation becomes: dd u2 — = — (2.25) dn gr Combining 2.20 with 2.22 yields: du u _ —+ - = 0 (2.26) dn r The above equations shows that u decreases and d increases as one moves from the inner to the outer bank. In addition, the log velocity law shows that u increases as one moves from the bed to the surface. The resulting pressure gradient on the bed dominates and produces an inward flow, with a corresponding outward flow at the surface. The simplest way to model secondary currents in bends is to assume constant radius of curvature and fully established flow within the bend. The distribution of transverse secondary currents with depth are described by many researchers, including Rozovskii (1961), Engelund (1974), Leschziner and Rodi (1979), Kalkwijk and De Vriend (1980) and Kikkawa et al. (1976). Some of these works present analytic solutions while others present the results of numerical experiments. In any case, the calculation procedures are somewhat lengthy. Several recent papers including Odgaard (1986) and Yeh & Kennedy (1993) have assumed that the transverse velocity profile is linear with depth. Aside from some boundary effects at the surface and the bed, this assumption seems to match with experimental data fairly well. The vertical velocity profile is assumed also to be linear, with values at the side walls set to ensure that continuity is maintained. To use the assumption of linear profiles, the value of the transverse surface velocity must be determined. Several researchers present analytic expressions, including Odgaard (1986): where: 2K2m u m = K— Jin & Steffler (1993) propose: vt = N—u R (2.27) (2.28) Rozovskii (1961) suggests: (2.29) where: Fj = 1.261 at the surface F 2 =.436 at the surface C = Chezy friction factor Finally, Kikkawa et al. (1976) use: -d v=u F, R K c a (2.30) where: FA = 3.33 at the surface F B = 1.11 at the surface Figure 2.5 shows the various formulae for secondary current surface magnitude plotted as functions of the ratio d/R. A l l of the formulae presented are linear functions of this ratio. Field data from Anwar (1986), Bathurst et al. (1979) and Rozovskii (1961) are also plotted. As an example, typical values for the river at Prince George (B=300 m, d=2 m, u =1.25 m/s, w*=0.1 m/s, R=450 m), used with equation 2.28 give a value of vs=0.058 m/s, or about 5% of u. For the Fraser River at Northwood, the above equations all tend to estimate that the maximum transverse velocity is in the range of 3 to 5% of the mean velocity. In the Fraser River, the value of the assumed vertical secondary currents are much smaller than those in the transverse direction. When the assumed linear profiles are applied, the resultant maximum vertical velocity is only 0.00039 m/s, two orders of magnitude smaller than the transverse velocity. This is due to the width to depth ratio of the river being so large. In reality, and contrary to the simplified model, secondary currents do not instantly develop the moment a bend begins. Rather, they gradually build up to a fully developed state as flow through the bend progresses. Additionally, the currents do not stop immediately, but die off gradually upon leaving a bend. Rozovskii (1961) presents the following formula for determining how far into a bend the flow must travel before becoming fully developed: 0,.„ = ^ | (2.31) where: 0 ] i m = Arc angle at which circulation is developed Rozovskii (1961) also provides an estimation of the downstream distance required for the circulation to drop to 10% of the fully developed value: 2 3C x~-^d (2.32) As discussed in chapter 4, the distances required for the circulation to develop and decay are quite small for the Fraser River, validating the assumption that secondary currents may be switched on and off immediately upon a particle entering and leaving a bend. One more formula is available which attempts to predict the value of a in the bend of a river. Rutherford (1994) suggests that in curved channels, the value is such that 1 < a < 3. Yotsukura and Sayre (1976) produced the following relationship from experiments performed on the Missouri River: a = - ^ - = 0.4j du* yU*R j (2.33) 2.11 Turbulence Intensity Nezu & Nakagawa (1993) provide the following semi-theoretical relations which specify turbulence intensity as a function of depth: (2.34) (2.35) (2.36) These equations are only appropriate away from the bed where turbulent energy production is approximately equal to viscous dissipation. An additional term can be added to the above formulae for use very near the boundary, but this term is neglected as the viscous sublayer in a typical river is negligible. The constants were obtained experimentally in a laboratory flume. These constants indicate that the turbulence intensities are such that u' > v' > w'. Nezu & Nakagawa (1993) state "Because smaller-scale eddies characterizing the fine structure of turbulence are nearly isotropic, they behave as random motions that can be described by the Gaussian probability density function." Thus velocity fluctuations in any one direction can be described by a normal distribution with a standard deviation equal to the turbulence intensity in that direction. Generation of a random normal variate can be accomplished using the following formula derived from Box and Muller (1958) and cited by Kottegoda (1979): (2.37) where: £ = standard normal deviate (ji = 0, a = 1) uvu2 = random uniform deviates (0 < w, < 1) A typical random walk step can then be generated by first determining turbulence intensity and then calculating the fluctuations as: AJC = u+u'£ Ay = v + v ' £ Az = w + w'£ (2.38) (2.39) (2.40) This chapter has provided the background and equations required to program a random walk model. The next chapter is concerned with checking the validity of some of these equations. Figure 2.1 - Definition sketch of channel geometry and velocity profiles, (a) Looking upstream; (b) Side view. Figure 2.2 - Various Frames of Reference, (a) Eulerian; (b) Lagrangian; (c) Instrumental 25 Figure 2.3 - Open Channel Flow Regions, (a) typical longitudinal velocity profile; (b) non-dimensional plot showing the various regions and their associated velocity laws. T3 \ IS 0.0 0.5 1.0 1.5 u (m/s) Figure 2.4 - Comparison of Velocity Laws. Values of u =1.25 m/s and u*=0.1 m/s were used, representing typical values of the Fraser River near Northwood. 0.40 I 1 1 1 1 I1 1 i I I I M i i i i i i i j i i i i i i i i i I i i i i i.'i i i i I i i i i i II i 0.00 I i i i i Ii i i i i I i i i i i i i i i I i i i i I i i i i i i i i i 1 i i i i I 0.000 0.010 0.020 0.030 0.040 0.050 d/R Figure 2.5 - Plot showing how secondary current surface magnitude normalized by the mean velocity varies with the ratio of depth to radius of curvature. Lines indicate theoretical equations, symbols indicate field data. Chapter 3 Laboratory Experiments 3.1 Purpose Laboratory experiments were carried out to collect turbulent velocity time series and to compare some of this data with equations for velocity distribution, shear velocity and turbulence intensity. Additionally, spectral plots are presented allowing for a comparison between estimated and observed maximum eddy size. 3.2 Equipment The experiments were performed in a flume in the Civi l Engineering hydraulics lab at U.B.C. The flume is 0.5 m wide, 1 m deep, and 12 m long. The flume has two pumps with a combined total recirculating capacity of 0.34 m-Vs, however only 0.17 m-Vs was available at the time of the experiments. Velocity measurements were accomplished with a Sontek 3-D acoustic Doppler velocimeter. Details on the theory of operation are available in Kraus et al. (1994). The probe works by sending out acoustic pulses into a small cylindrical measuring volume. The Doppler shift caused by particles in the flow is then used to calculate the velocity components in three directions. The measuring volume is 5 cm away from the probe head, making the system non intrusive. Measurements can be taken at up to 25 Hz. The probe also has the capability of determining the location of the sampling volume relative to the bed (or surface), allowing for very precise positioning of the probe. This is particularly useful when working near the bed. 3.3 Experiments Two different cases were studied with the depth of water set to 9 cm (known as the shallow experiment) and 50 cm (the deep experiment). The flume was left in a horizontal position. Two different probe heads were used. One of these heads permits measurements in a downward looking direction, the other in an upward direction. Combined use of the probes allowed for a complete velocity distribution to be obtained for the deep flow. The shallow flow combined with the probe head geometry prohibits measurements of the upper region, so only the wall region was studied in this case. A l l measurements were taken in the middle of flume, away from the side walls. Work in progress by Gomm (personal communication) shows that this flume does not produce large boundary layers on the sides, so measurements taken in the middle of the flume should be unaffected by the side walls. The probe was suspended using an aluminum mounting rack which was specially designed for this flume. Each experiment consisted of the establishment of steady flow in the flume followed by systematic measurements taken at various vertical locations. The probe was moved manually, while the acquisition software displayed vertical location to sub-centimeter accuracy. The probe was then secured in place for acquisition of the velocity signal. Record lengths of just under three minutes were collected at 25Hz, resulting in 4096 data point time series suitable for analysis with an FFT routine. The assumption is made when taking these measurements that turbulence is stationary (i.e. it is assumed that the statistical characteristics of turbulence do not change with time). The collection of the data was carried out using Sontek's proprietary software, while processing of the data was completed using Visual Numerics P V - W A V E . Collection and analysis were performed on PC platforms. Sample output from the velocity probe is shown in figures 3.1 and 3.2. Figure 3.1 shows data from z/d=0.2 for the deep case, while figure 3.2 shows the output from the shallow case. Shown are the fluctuations only (the mean values have been removed). The flow is clearly turbulent, with a Reynolds number (based on mean velocity) for the deep flow equal to 2.8x10^. Results from the shallow case are also turbulent, having a Reynolds number equal to 3.5xl()4. Both time series contain a large amount of high frequency information, however one can also clearly see larger scale fluctuations in the deep case. The shallow case also contains some low frequency information, but it is much less pronounced. To quantify these results, figures 3.3 and 3.4 show the power spectral density (PSD) plots for each experiment. In the case of the deep experiment, figure 3.3 shows a significant amount of energy in the lower frequencies. Low frequency energy is seen to be one to two orders of magnitude larger than the high frequency information. This tends to support the idea that large, low frequency eddies are responsible for the majority of turbulent transport. Figure 3.4 shows that there is much less energy in the lower frequencies of the shallow flow. This is expected as the maximum size of an eddy in the shallow flow is much smaller than that of the deep flow. To examine the viability of using d I ut as an order of magnitude estimate of the period of the largest eddy, figure 3.5 shows the same spectral information except that the spectra have been smoothed and plotted so as to concentrate on the low frequencies. Figure 3.5 shows that for the deep experiment, the peak value of the PSD corresponds roughly with the predicted largest eddy period of d I u*. The figure also shows that the predicted period for the shallow experiment does not correspond as well to the measured spectra. This is not necessarily incorrect, as d I w, is merely an estimate of the largest eddy. Additional effects such as low frequency variations in the flow rate of the flume, geometric irregularities, and errors in the determination of ut may all contribute to differences between the observed and maximum eddy size. The mean velocities from each time series are assembled into composite longitudinal velocity profiles, as shown in Figure 3.6. In the lower portion of the flow (z/d<0.2) the velocities in both experiments follow the log law quite closely. When equation 2.14 is applied, shear velocities of 1.43 cm/s and 1.2 cm/s are calculated for the deep and shallow experiments respectively. The upper portion of the flow in the deep experiment appears to be fairly uniform, with some curving back of the velocity as surface effects become important. It should be noted that the reduction of velocity near the surface is not handled by the wake functions discussed in chapter 2. This is due to the assumption that the maximum velocity occurs at the surface, which is not always true. This is not of concern for comparison reasons as even when the log law is used to describe the velocity over the entire depth, a maximum error of only 10% occurs at the surface. Similar to the composite velocity profiles presented, Reynolds stress distributions can be plotted as functions of depth such as those shown in figure 3.7. In the lower portion of the flow where turbulent kinetic energy is being generated due to shear, the Reynolds stress takes on a positive value which allows for determination of a shear velocity near the bed. At the peak value, the shear velocity is calculated by equation 2.15 as 1.64 cm/s for the deep flow, and 1.14 cm/s for the shallow flow. These values are similar to those calculated from the log velocity method (1.43 cm/s for the deep flow, 1.2 cm/s for the shallow flow). The standard deviation of the velocity signals as functions of depth may be described by the turbulence intensity equations (2.34-2.36). These equations may be plotted along with the experimental data to show the amount of agreement. Figure 3.8 shows that in the upper region of the deep flow, where there is little shear, all three of the distributions seem to show no trend. In the lower region where shear exists, however, trends are clearly visible. According to Nezu & Nakagawa (1993), turbulence intensity should increase with depth. Agreement with the deep experiment is mixed. The predicted trends correspond with longitudinal and transverse measurements, although the vertical measurements seem to contradict Nezu & Nakagawa (1993). In the shallow experiment, shown in figure 3.9, much better agreement with the predictions is found. In both experiments, however, the constants given in equations 2.34-2.36 seem to under predict the longitudinal and transverse measurements and over predict the vertical measurements. This chapter has attempted to verify the validity of several equations presented in chapter 2. The equations for longitudinal velocity distribution (table 2.3), turbulence intensity (equations 2.34-2.36), largest eddy period (d/w») and shear velocity for open channel flow (equations 2.14 and 2.15) have all shown reasonable agreement with the laboratory data. Unfortunately, facilities are unavailable to confirm the validity of equations for secondary currents. Section 2.10 attempted to remedy this by comparing the equations presented in chapter 2 with some observed field data. Chapter 4 will now describe the computer model which was developed to apply these equations to simulate the dispersion of a pollutant in a river. Figure 3.1 - Time series of velocity fluctuations taken at 25 Hz for deep experiment, (a) longitudinal direction; (b) transverse direction; (c) vertical direction. 34 201 ~i 1 1 r -20 _i i i_ ~i 1 1 1 1 r~ (a) _l I L_ 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 Time (sec) -10 h-~i 1 r (b) _j i i i i i i i_ 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 Time (sec) 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 Time (sec) Figure 3.2 - Time series of velocity fluctuations taken at 25 Hz for shallow experiment, (a) longitudinal direction; (b) transverse direction; (c) vertical direction. 35 Frequency (Hz) Figure 3.3 - Power spectral density plots for deep experiment, (a) longitudinal direction; (b) transverse direction; (c) vertical direction. io -I , . , , • , , , , , , .11 1 . . 1 0 0 0.5 1.0 1.5 2.0 Frequency (Hz) Frequency (Hz) Figure 3.4 - Power spectral density plots for shallow experiment. (a) longitudinal direction; (b) transverse direction; (c) vertical direction. 37 1.0000 0.1000 ac < < E o.oioo 0.0010 b-0.0001 0.00 < 1.0000 0.1000 E 0.0100 a a! 0.0010 b-0.0001 0.00 1.0000 N 0.1000 < m \ < E 0.0100 0.0010 0.0001 0.00 Deep Exper iment Shallow Exper iment (a) 0.20 0.20 0.20 Figure 3.5 - Filtered low frequency power spectral density plots, (a) longitudinal direction; (b) transverse direction; (c) vertical direction. The solid vertical line indicates the value of d/u* for the deep experiment, while the dashed vertical line indicates the value of d/u* for the shallow experiment. T3 0 10 20 30 40 50 60 u (cm/s) Figure 3.6 - Longitudinal velocity distributions showing predicted log velocity law and experimental values. - 3 - 2 - 1 0 1 2 3 -u'w' ( c m -2 / s -2 ) Figure 3.7 - Reynolds stress distributions showing experimental values. (a) Standard Deviation (u) (cm/s) 1.0 0 1 2 3 Standard Deviation (v) (cm/s) 0 1 2 3 Standard Deviation (w) (cm/s) Figure 3.8 - Standard deviation of velocity signals for deep experiment, (a) longitudinal direction; (b) transverse direction; (c) vertical direction 41 T3 0.30 - i—i—i—.—i—i—i—i—i—i—i—.—i—i—i—i—.—i—|—,—,—,—,—,—,—,—,—, 0.25 — \ \ (a) : 0.20 \ V 0.15 0.10 <^ > F lume Exper iment \ \ ~~ 0.05 — Nezu & Nakagawa Predic t ion \ <d> — 0.00 ' ~ Standard Deviation (u) (cm/s) 0.30 0.25 0.20 — 0.15 0.10 0.05 — 0.00 F lume Exper iment — Nezu & Nakagawa Predic t ion o (b) Standard Deviation (v) (cm/s) 0.30 - (c) : 0.25 \ 4 0.20 \ . • -E 0.15 \ H 0.10 \ <>-^> Flume Exper iment _z 0.05 ^ \ Nezu 8c Nakagawa P r e d i c t i o n s 0.00 - , 13 \ El Standard Deviation (w) (cm/s) Figure 3.9 - Standard deviation of velocity signals for shallow experiment, (a) longitudinal direction; (b) transverse direction; (c) vertical direction Chapter 4 Model Development - PATRAC This chapter is intended to provide a description of the computer model. It is not intended as a detailed operations manual, which is provided in Appendix A. A flow chart is provided as figure 4.1, however, to enable the reader to better understand the sequence of operations performed by the model. 4.1 General Specifications P A T R A C (short for PArticle TRACking model) is a three dimensional Lagrangiah particle tracking model written in Visual Basic. The general equations that P A T R A C uses to follow particles through the simulated flow field are as follows : Ax = u + u (4.1) At Ay -21 = v + v' + v j c (4.2) Az At where: = w + w' + wsc + wfaIl (4.3) vsc,wsc are the secondary current velocities wfall is the fall velocity of the particle The model is capable of simulating either an instantaneous or continuous release of particles. The duration of the continuous release that the model can sustain depends on the total number of particles in the model as well as how many are introduced per time step. The maximum number of particles that P A T R A C can handle is 30000. Increasing the number of particles greatly increases the run-time of the simulation, while having too few particles does not afford enough resolution when attempting to determine concentrations. A balance must be found between acceptable run-time, resolution, and desired duration of release, as described in section 4.12. P A T R A C is very flexible with regards to specification of initial conditions. Each particle can be given a separate initial starting position, allowing for complex release conditions. For tracer particles, these release conditions can either be the port locations of a diffuser, or they can be the coordinates and concentration distributions obtained from a near field mixing model. 4.2 Geometry P A T R A C assumes a rectangular channel. This assumption is made mostly for computational convenience, but is also due to a relative lack of site specific geometric information. The width to depth ratio of the Fraser River is large, so it is assumed that errors introduced due to the vertical side walls of the rectangular channel are small. Many rivers including the Fraser may have a compound structure consisting of a deep fast moving channel on one side offset by a shallow slow moving channel on the other side. If more site specific information of this type becomes available, the model could conceivably be reprograrnmed to use such information. Information on the width of the channel, determined from 1:50000 topographic maps, is stored in a separate file. This file contains the width of the river digitized every 200 metres downstream from an initial starting point. P A T R A C performs a linear interpolation between these points to give the width of the channel at any point as a function of distance downstream. The radius of curvature of the channel, also obtained from topographic maps and stored at 200 metre intervals, is stored in another file. This information is used to help determine the magnitude of secondary currents generated by flow through the bends in the channel. Given the flow rate and width of the channel, the model calculates average depth of flow and mean velocity by using continuity and the wide channel approximation of the Manning equation in a rearranged form : 4.3 Flowlines To correct for the changing width and depth of the channel, correction velocities are calculated and added to each particle during each time step. If the model is run in laminar mode with no turbulent fluctuations, the resultant paths of the particles represent laminar flowlines. When viewed from above, particles in the centre of the channel follow a straight course, while particles on the banks of the channel follow the banks. When viewed from the side, particles on the surface follow the surface, while particles on the bottom follow the changing depth of the channel. 4.4 Mean Longitudinal Velocity Distribution P A T R A C can use either of the longitudinal velocity distributions described in the literature review. Unless detailed river measurements show a need otherwise, the log law (table 2.3) is used because of its simplicity. 4.5 Mean Vertical And Transverse Velocity Distributions The assumption is made that the mean components of velocity in the vertical and transverse directions are zero, although they need not be. Secondary currents, described in section 4.7, make these terms non-zero in the bends of the river. However in straight sections, without detailed river measurements, these components are assumed to be zero. (4.4) 4.6 Shear Velocity Recall from chapter 2 that there are three methods for determining shear velocity for open channel flow. The simplest is equation 2.13, which describes shear velocity as a function of depth and slope of the channel. The other equations (2.14 and 2.15) allow determination of the shear velocity from a known distribution of velocity or Reynolds stress. Nezu & Rodi (1986) found that, in general, the values obtained by equation 2.13 did not vary by more than 5% from equations 2.14 or 2.15 in controlled laboratory flume experiments. This is fortunate, as equation 2.13 is a quick method of determining shear velocity which lends itself well to the very repetitive calculations performed by P A T R A C . 4.7 Secondary Currents Recall that figure 2.5 shows the various formulae for secondary current surface magnitude plotted as functions of the ratio d/R. A l l of the formulae presented are linear functions of this ratio. Field data from various sources of literature are also plotted. The data seems to best fit the formula provided by Jin and Steffler (1993), although at the d/R ratio for the Fraser River there seems to be little difference between the formulae. In view of this, equation 2.28 (Jin & Steffler 1993) is chosen for use in the model. This equation is simple to apply and requires no information which is not already stored in the model for other purposes. This makes the equation computationally efficient and a good choice for use in P A T R A C . The model uses the value of vs predicted by equation 2.28 along with the assumed linear distributions described in the literature review for calculation of the secondary currents in bends. For simplicity, secondary currents are activated immediately upon a particle entering a bend, and are shut off immediately upon a particle leaving a bend. Justification for this assumption comes from equations 2.31 and 2.32. For the conditions in the Fraser described in section 5, the value of 0 l i m =12° , which translates into an arclength of approximately 95 metres. Additionally, the approximate distance for the circulation to drop to 10% of the fully developed value is again about 95 metres. At the Northwood site, the simulation ends before the bend does, so we need only be concerned with the bend entrance. The length of reach in which this transition would be important is only about 7% of the total simulation length, and thus the effect is not included in the simulation. Rozovskii (1961) provides formulae for calculating the magnitude of the circulation in this transition region. These equations could be included in the program code at the expense of additional computing time. 4.8 Turbulent Fluctuations Turbulent fluctuations are generated from a Gaussian random number generator as described by equation 2.37. Equations 2.38-2.40 are used to describe the magnitudes of the fluctuations. The question of how long to apply these fluctuations can be answered by using equation 2.11 along with some consideration of the scales involved. If one assumes that the largest eddy size is on the order of depth d, and that the characteristic velocity of the eddy is the shear velocity ut, then the ratio of d/u* gives the time scale of the eddy. To check this assumption, equation 2.11 can be applied. If Ay is set to d, and At is set to d I w*, then by algebraic manipulation one can produce the following result: a = -J- = — - — = -T- = - = 0.5 (4.5) da* 2dutAt 2du — ^ u* Recall from chapter 2 that from field observations a = 0.6 ± 50%. Equation 4.5 shows that the assumptions of eddy size and timestep produce a random walk growth rate which is well within the observed values of a. Thus the time step chosen for use in the model is the average value of d I ut for the given stretch of river being simulated. 4.9 Transverse Calibration The transverse calibration of the model involves simulating a single point source in a rectangular channel which has dimensions and flow characteristics approximate to those of the Fraser River at Prince George. The channel is made wide enough so that the boundaries will not affect the growth rate of the plume. The growth of the plume can be seen on figure 4.2(a). The width of the plume is defined as enclosing 95% of the particles. Also plotted on figure 4.2(a) is the curve for 8t/d«»=0.5, along with the 50% error bounds associated with it. The plume predicted by the model is clearly within the error bounds. Figure 4.2(b) shows an overhead view of the simulation results frozen after 30 minutes of simulation time (90 timesteps of 20 seconds each). 4.10 Vertical Calibration The assumption that the largest eddy period is equal to d I ut combined with the long time step associated with this eddy results in a coarse resolution in the vertical direction. This is because after one time step, most of the particles will have already mixed across the depth. The complete growth of the plume in the vertical direction after only one time step precludes precise checking of the vertical mixing rate. As a simple check, simulation 1 (described in chapter 5) is tested with the equation from Fischer et al. (1979) in table 2.1. This equation predicts that complete vertical mixing will occur within 39 metres of release. One time step for this simulation is 20 seconds, resulting in complete vertical mixing at 25 metres downstream. The model over predicts vertical mixing, but the result is just slightly more than 50% larger than the Fischer et al. (1979) equation. In order to obtain a better vertical resolution, a model would need to break each timestep down into a series of substeps, which would allow for a better vertical resolution and hence the ability to calibrate the model more accurately. Such a model would require the number of calculations to be at least one order of magnitude larger than the current model. This small improvement in the model is not justified given the large amount of additional computing time associated with it. 4.11 Model Output Before beginning simulation, the user may select how often the model generates an output file. These output files consist of the x, y, and z coordinates of each particle at the instant the output file is generated. Any one individual file is of little value, however when several files are averaged the results stabilize and produce useful information. Section 4.12 discusses this in greater detail. The output files can be interpreted in several ways, depending on what the user wishes to look at. If the coordinates of each particle are plotted on a scatter plot such as figure 4.2(b), the resulting image provides a good indication of the size and shape of the plume. P A T R A C can calculate the location of the plume edges, based on a user definable criterion. This is typically the 95th percentile of particle concentration. Another useful result can be obtained by processing the data files with the histogram routine of P A T R A C , which generates 2 dimensional histograms. These histograms can be used to determine downstream concentrations with the help of contour maps, as demonstrated in chapter 5. 4.12 Choice of Suitable Random Walk Parameters In general, the random walk process has three main parameters which may be specified by the user. These parameters are the timestep, the number of particles released per time step, and the number of timesteps over which an average is taken. The timestep is not available for user selection in PATRAC, as this is based on a physical process. This leaves two parameters which must be chosen so as to produce stable and accurate answers. The results of several simulations attempted, varying the number of particles released per step as well as the number of timesteps averaged are shown in figure 4.3. The error shown represents the difference between the simulation and the theoretical results predicted by Fickian diffusion. It is clear that given a sufficient number of particles, the error stabilizes quickly so that only 20 or so time steps need be averaged to obtain a stable solution. Any increase in the number of timesteps used produces no benefit in error reduction. Combining this fact with the added run time required to simulate additional time steps makes the optimal value clearly visible on figure 4.3. The effects of the number of particles released per time step is also shown in figure 4.3. From the figure, and also intuitively, more is better. It is clear that even after the simulation has achieved stable results, the absolute value of the error is affected by the density of particles. The value of 80 particles/time step produces negligible error, while halving the number of particles produces some error, but the results are still quite acceptable. The value of 20 particles/time step yields unacceptable results. It is clear that the effect of the number of particles on the accuracy of the simulation is non linear. The easiest way to quantify this effect is with a graph such as figure 4.3. It is recommended that such a graph be prepared whenever a random walk model is used to ensure accurate and stable answers. For all simulations completed in this thesis, 80 particles/time step were used. It is important that the simulation is run over enough time steps to produce a stable result. Once stability is achieved, additional time steps are wasted effort, and the computational effort should be focused on the number of particles used. It is recommended that the number of particles released per step be on the high side to produce the lowest error possible. In addition to reducing the absolute error in plume width detection, having more particles will produce better resolution when attempting to determine downstream concentrations. 50 Initialize variables ' (End)-* Yes Release particles for this timestep Calculate width, depth, and mean velocity I Calculate shear velocity I Calculate longitudinal velocity from log law I Calculate secondary current velocities T Calculate turbulent fluctuations I Calculate flowline correction velocities I Move particle to new location I Check boundary conditions and correct Generate output file Figure 4.1 - Flowchart of the random walk engine of P A T R A C . Figure 4.2 - Transverse calibration of a single plume, (a) Plume width plotted as a function of downstream distance and compared to Fickian theory (a=0.5); (b) Overhead snapshot of plume after 30 simulation minutes. 0 20 40 60 80 Number of Timesteps Averaged Figure 4.3 - Sensitivity/stability analysis of random walk simulation. Error (defined as £ (Pobs ~ Psimf i s plotted against the number of timesteps averaged for various particle release rates. Chapter 5 Application To The Fraser River P A T R A C is used to simulate conditions downstream of the Northwood Pulp and Paper mill, located in Prince George, B:C. (see figure 5.1). The Northwood plant is located just upstream of a large bend in the river, which has a radius of curvature of approximately 450 m. The plant disposes of effluent into the river via a multiport diffuser; details of which are available along with an experimental study of the near-field discharge in Marks (1996). The simulation is chosen to end where the channel suddenly widens, approximately 1400 metres downstream of the diffuser. This sudden extreme change in geometry is accompanied by the division of the river into two channels, which is beyond the scope of the current study. To model the 1400 metre section of the river downstream of the diffuser, width and curvature measurements were taken from Energy, Mines and Resources Canada topographic map 93 G/15, Prince George, edition 4. The scale of the map is 1:50000. Validation of the model using the popular software package C O R M I X is performed, along with comparison of the model predictions to field data. The model is then used to predict plume width and downstream dilution for several river flow rates. 5.1 Validation of the Northwood Site Simulation In this scenario, all of the particles used are designated as tracer particles. These particles will exhibit no inertial effects and will travel perfectly with the simulated fluid. Flow rates for the simulations were taken from a compilation of the entire record of flows from the Shelly hydrometric station, located upstream of the site. Annual hydrographs from the years 1951-1991 were averaged to produce one hydrograph from which flows could be selected. This data is available from Water Survey of Canada on a CD entitled H Y D A T . The bed slope in the area is assumed to be 0.001, as taken from Northwest Hydraulic Consultants (1994). Manning's n is assumed to be 0.03 (Henderson (1966)). Following recommendations from section 4.12, the number of particles used in the simulation is 30,000. This is approaching the upper limit on an array size in Visual Basic, and thus sets the number of particles available. Table 5.1 summarizes the simulations which were performed. Table 5.1 - Summary of Simulations Sim # Simulated Date Simulated Flow (m3/s) Description 1 Sept. 15, 1993 805 Constant width channel (300 m) intended as a simplification of Northwood site for CORMIX comparison 2 Sept. 15, 1993 805 Northwood site with varying width, depth, and secondary currents 3 Low Flow 150 Low flow case 4 High Flow 2370 High flow case 5.1.1 Comparison of PATRAC With Cormix For validation, C O R M I X was used to simulate the Sept. 15, 1993 flow. Field data is also available for this date, as described in the next section. CORMIX, short for CORnell MLXing model, is a popular mixing model marketed by Gerhard Jirka of Cornell University, New York. A description of the program is available in Jirka and Akar(1991). Figure 5.2(a) shows the width of the plume as a function of distance downstream. The width of the plume is obtained using the method described in chapter 4. Also plotted is simulation 1 from P A T R A C , which is intended to be directly comparable to the C O R M I X output. Two interpretations of the CORMIX output are shown. C O R M I X does not predict that bank attachment occurs, and thus predicts a much smaller plume which does not touch the right bank within the 1400 metres of simulated reach. This contradicts the results of both P A T R A C and the field data which clearly show that bank attachment does occur. If it is assumed that the CORMIX plume attaches to the right bank immediately, the results between CORMIX and P A T R A C are similar and in quite good agreement. Both sets of output also compare acceptably with the field data provided by Borstad and Associates (discussed in the next section). In any case, the disagreement involves the width of the plume, and not the growth rate which is very similar for both CORMIX and PATRAC. Differences in the width of the plume may be partially accounted for by near-field conditions, which CORMIX attempts to simulate but P A T R A C does not. It is possible that slightly different conditions would cause C O R M I X to predict bank attachment. The overhead view of the simulation, frozen after 30 minutes of simulated time, is shown in figure 5.2(b). The image is transparent in depth, so that all particles regardless of their vertical location are shown. The plume is seen to grow with downstream distance. The plume visibly begins to touch the right bank of the river at a distance of 600 to 800 metres downstream from the diffuser. At this point, the growth of the plume is hampered on one side, as can be seen in figure 5.2(a) which clearly shows the plume growth rate slowing. 5.1.2 Comparison of P A T R A C with Field Data Digital information is available to correspond with the September 15, 1993 simulation from Borstad and Associates. The image was obtained by a Compact Airborne Spectral Imager (CASI), which was flown over Prince George in a small plane. The image was provided as a TIFF file, which was then processed with Image 1.57 for the Macintosh. Figure 5.1 shows the digital image. Digital measurements of width of river and width of plume were taken every 100 m downstream using the image analysis program. The growth of the plume as predicted by the model for the September 15, 1993 simulation of the Northwood site is shown in figure 5.3(a). Also plotted is the data taken from the CASI. Figure 5.3(a) shows that the field data contains a large amount of scatter. This scatter is due to uncertainty in the width measurements taken with the image analysis software. The uncertainty is partially due to a fairly course pixel resolution, and partially due to the uncertainty in detecting the plume edge. The most disagreement occurs during the first 500 metres, where the field data indicates that the plume width actually decreases. One possibility for this discrepancy is that the thalweg of the river may cut across from the left to the right bank of the river in this area. It is possible that this could influence the plume to move towards the right bank, resulting in a decrease of the plume width in this region. Effects of the location of the thalweg are not accounted for in P A T R A C . Another explanation is that the CASI may not "see" into the entire depth of the river. If it is only detecting the plume in the upper layer, than the true width of the plume can not be determined from the digital data. In general, however, the simulation does predict the growth to within a reasonable accuracy, as shown by the error bars on the field data. A snapshot of the plume which has been frozen in time is shown in figure 5.4. To aid in visualization, the image has been remapped onto a curvilinear coordinate system, even though the model does not use curvilinear coordinates. The transformation results in slight distortion of the plume, but the visualization benefit is quite useful. Comparing figure 5.1 with 5.4, good qualitative agreement is achieved both in the geometric simplification of the model as well as in the plume location. 5.2 Extreme Flow Case Simulations Given the extensive run times involved with the random walk model, it is advisable to carefully choose the desired cases for simulation. Two obvious choices include the low and high flow rates expected in the river. Figure 5.3(b) shows the plume width results from simulations 2, 3 and 4. Generally, results are as expected, with the higher flow rates producing more dispersion. The high flow case shows several effects which are visible but not as clearly in the other cases. Firstly, at a distance of 500 metres downstream the growth rate of the plume slightly increases, corresponding to the point where the bend begins in the simulation. A more significant increase occurs when the channel contracts sharply, which begins at x=900 metres. The growth rate of the plume increases sharply in the high flow case, while only mildly affecting the medium and low flow cases. Figure 5.5 helps to highlight why this occurs. Figure 5.5(a) shows the value of the secondary current surface velocity normalized by the mean velocity, while figure 5.5(b) shows the shear velocity. Both are plotted as functions of downstream distance. Figure 5.5(a) indicates that the secondary currents in the low flow case are minimal throughout the entire simulation. The high flow, and to a lesser extent the medium flow, show a substantial increase in the magnitude of the secondary currents at the contraction. This increase in secondary current velocity increases the transverse shear associated with the secondary current velocity profiles. This results in increased dispersion, and is part of the reason why the plume width is seen to jump sharply at the constriction as shown in figure 5.3. Additional dispersion is caused due to the change in shear velocity, as shown in figure 5.5(b). Equations 2.34-2.36 link the shear velocity to the turbulence intensities, so it is logical to assume that an increase in the shear velocity is accompanied by an increase in the turbulence intensities. Since the random walk step is a function of the turbulence intensities, the conclusion is that increasing the shear velocity results in an increase in the plume growth rate. 5.3 Calculation of Downstream Dilution Contour maps of dilution, calculated assuming an effluent flow rate of 1.5 m^/s (see Marks 1996) are shown in figure 5.6. The dilutions are calculated by the routine described in Appendix A. The maps are not plotted in a curvilinear system as it is desired to avoid distortion with the contour lines. As one would expect, increasing the flow rate clearly increases the amount of dilution. The maps give clear results everywhere except near the right bank, where the contouring is difficult to interpret due to the rapid contraction of the channel. In addition, the contouring package used does not perform well near a boundary. To supplement the results, figure 5.7 shows centreline and bank dilution as functions of downstream distance. The results are the same as presented in the contour maps, only with added clarity near the boundary. A curious oddity occurs in figure 5.6(c), with a resulting effect on figure 5.7. The high flow of 2370 m-Vs produces a large time step of 30 seconds, as determined from the method described in section 4.8. As a result, the overhead view of the plume no longer appears to produce a continuous release of particles, but rather appears to "puff" particles once every 30 seconds. This is visible on figure 5.6(c) near the diffuser as rings of low dilution surrounded by areas of high dilution. The puffs eventually catch up with each other due to longitudinal shear in the simulation. This result carries over to figure 5.7, where a sharp rise and dip in dilution occurs on the high flow lines. The predicted bank dilution never becomes stable, whereas the low and medium flows clearly show that bank dilution is decreasing with downstream distance. It is clear that the simplicity of a single time step simulation is beginning to show at the higher flow rates, which could be refined if a "substep" model such as mentioned in section 4.10 were developed. This chapter has shown that a random walk model can produce results which compare reasonably with CORMIX and available imagery. The model has been used to predict low, medium, and high river flow dilutions and plume widths at the Northwood site in Prince George. Plots of plume width versus downstream distance show the growth of the effluent plume. Contour maps of dilution show dilution spatially except near the bank where they are unclear. To clarify this area, a graph which shows centreline and bank dilutions versus downstream distance has also been provided. Figure 5.1 - Digital image of the Northwood site on the Fraser River. The plume from diffuser is clearly visible. The river flows counterclockwise through the bend. Image provided by Borstad Associates Ltd. (a): o • A Borstad Northwood PATRAC Straight Channel CORMIX Straight Channel(No Bank Attachment) CORMIX (Assuming Instant Bank Attachment) 200 400 600 800 1000 Downstream Distance (m) 1200 1400 Figure 5.2 - Comparison of PATRAC and CORMIX using a straight, rectangular channel, (a) Plot of plume width vs. distance downstream as predicted by P A T R A C and CORMIX, compared against data from Borstad and Associates, (b) Overhead snapshot of the plume. Figure 5.3 - Plume width as predicted by the Northwood site simulation with P A T R A C . , (a) P A T R A C predicted plume width compared with Borstad and Assoicates data; (b) Different flow rates compared : High flow = 2370 m3/s, Intermediate flow = 805 m3/s, Low flow =150 m3/s. Figure 5.4 - Overhead snapshot of the Northwood site at 805 m3/s flow. The orthogonal coordinates used by P A T R A C have been remapped on to a curvilinear system, resulting in slight distortion of the plume. 64 0.10 0.08 h 0.06 o xn > 0.04 0.02 0.00 High Flow Intermediate Flow <y—0 $ e — o-Low Flow • B B B B B- - B -(a) * * o -B B B 200 400 600 800 1000 Distance Downstream (m) 1200 1400 u CO xn D 0.20 0.15 0.10 0.05 \ -0.00 High Flow * e — * * -Intermediate Flow * * * * * * Low Flow • B B B B B B-(b) * * * * -B B B 200 400 600 800 1000 Distance Downstream (m) 1200 1400 Figure 5.5 - Variation of parameters affecting transverse mixing given various flow rates and distance downstream. High flow = 2370 m3/s, Intermediate flow = 805 m3/s, Low flow = 150 m3/,s. (a) Variation of secondary current surface velocity normalized by mean velocity; (b) Variation of shear velocity. Figure 5.6 - Contour maps of dilution for various river flow rates, (a) Low flow = 150 m3/s; (b) Intermediate flow = 805 m3/s; (c) High flow = 2370 m3/s. 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 Distance Downstream (m) a o -4-> Q 500 400 300 200 100 01 C.L. C.L. C.L. (b) High Flow Inf. Flow Low Flow / | \ / [ \ / T \ /TS / [ \ / | \ 3t 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 Distance Downstream (m) Figure 5.7 - Dilution plots for the Fraser River at Northwood for various flow rates, (a) Dilution at the right bank of the river, (b) Dilution in the centre of the channel. Chapter 6 Conclusions and Recommendations 6.1 Summary and Conclusions From the work presented, conclusions can be made with regards to various aspects of the research: Velocity Field Experimental investigation in a laboratory flume shows that in general the equations to describe the velocity field presented in chapter 2 perform as expected. Specifically: • the log velocity law performs well in the lower 20% of the flow. Although not intended to, the law is commonly applied over the entire depth of flow. Experiments conducted show a maximum error of approximately 10% when this approach is used. It is still desirable to obtain field measurements to confirm the velocity structure. • the three methods presented to determine the shear velocity are found by Nezu and Nakagawa (1993) to produce results within 5% of each other in controlled experiments. Experiments performed with an acoustic Doppler velocimeter confirm that the log velocity method and Reynolds stress method of calculating shear velocity produce similar results. Collection of field data to determine the shear velocity may prove difficult or too expensive, but a good estimate is still available if the slope of the channel and depth of flow are known. • numerous equations are available to describe the flow patterns generated by secondary currents in the bend of a river. Their effect on transverse mixing is not well understood (Rutherford 1994). The diffusion model developed fails to predict the large transverse mixing coefficient in a bend that is suggested by some researchers. • the turbulence intensity equations presented by Nezu and Nakagawa (1993) are acceptable estimates of how the turbulence intensity in open channel flow varies with depth. The constants do not quite agree with results from these experiments, but a much larger suite of experiments would need to be conducted to confirm this. Both experiments and equations agree that u'>v'>w'. Random Walk Modelling of Pollutant Dispersion • the random walk process can be used successfully to model pollutant dispersion in an open channel flow. • the choice of the random walk step size (either spatial or temporal) is critical to producing reasonable estimates of dispersion. The step size must not be chosen to satisfy computational efficiency, but rather to emulate the important low frequency eddies which are responsible for the majority of pollutant transport. • the stability of the results from a random walk model is found to be a function of the number of timesteps over which the results are averaged. The absolute error of the results is also affected by the number of particles used in the simulation. • the computational cost associated with this method detracts from the practicality of using this model to simulate large numbers of cases. Faster equipment, especially parallel processor computers, are desirable to make application of this method more practical from an engineering point of view. Application to the Northwood Site on the Fraser River • acceptable agreement with field data and CORMIX indicate that the model produces reasonable results when applied to the Northwood site. Estimates of plume width and dilution are provided for low, intermediate, and high flow rates which are expected at Prince George. • the model suggests that secondary currents have a moderate to strong effect on mixing where the channel contracts. • from the model predictions and available field data, the effluent plume from Northwood seems to be sharply confined to the right bank of the river until the sudden channel expansion approximately 1400 metres downstream of the diffuser. • this method could be applied to other areas of the Fraser, such as the CanFor mill located just downstream of the Northwood site. 6.2 Recommendations • P A T R A C was designed to provide an estimate of diffusion using a minimal amount of field data. Formulae obtained from the literature are used to specify the flow field and turbulence intensities, and a simple rectangular shaped channel is assumed. Most of these equations were derived using data collected in laboratory flumes and may not be as appropriate as equations derived from actual river measurements. P A T R A C , as with any model, can only benefit from the collection of field data. The geometry of the channel and the velocity structure are the most valuable. The straightforward programming used to develop the random walk model allows any of the formulae to be replaced, thus the user may wish to accommodate any available field data in this way. • the stability and accuracy of a simulation will vary with the site chosen and the conditions that the model is required to predict. Err on the side of caution and always produce a stability/sensitivity chart to ensure the results are stable and accurate. 6.3 Further Research Many of the subjects in this thesis, both directly and indirectly related to random walk modelling, have not been fully examined. More research could be done to: • obtain detailed bathymetry and velocity measurements at the Northwood site. • perform a mid-field die study to delineate the plume width and determine dilution. • study secondary currents in the field, including their effect on transverse mixing coefficients. • develop a numerical finite difference model to produce a more accurate velocity field from which the random walk model could operate. • develop a multiple time scale model, which would remove some of the difficulties associated with assuming a single time step equal to d/u*. • use the existing model to develop graphs of how various quantities such as dilution and plume width vary as functions of input variables (such as flow). REFERENCES Anwar, H.O., 1986. "Turbulent Structure in a River Bend", Journal of Hydraulic Engineering, 112(8), pp. 657-669. Coleman, N.L. , 1981. "Velocity Profiles with Suspended Sediment", Journal of Hydraulic Research, 19(3), pp. 211-229. Coles, D., 1956. "The Law of the Wake in the Turbulent Boundary Layer", Journal of Fluid Mechanics, 1, pp. 191-226. Csanady, G.T., 1972. Turbulent Diffusion in the Environment. D. Reidel Pub. Co., Boston. Daugherty, R.L, Franzini, J.B., Finnemore, E.J., 1985. Fluid Mechanics with Engineering Applications. McGraw-Hill, New York. Engelund, F., 1974. "Flow and Bed Topography in Channel Bends", Journal of the Hydraulics Division ASCE, 100(HY11), pp. 1631-1648. Evans, W., 1996. "Flocculation of Pulp M i l l Effluents with Fraser River Sediment", Thesis submitted in partial fulfillment of the degree Master of Applied Science, University of British Columbia, Vancouver, Canada. Fischer, H.B., List, E.J., Koh, R.C.Y., Imberger, J., Brooks, N.H. , 1979. Mixing in Inland and Coastal Waters. Academic Press, New York. Henderson, F .M. , 1966. Open Channel Flow. Macmillan, New York. Heslop, S.E., Allen, C M . , 1993. "Modelling Contaminant Dispersion in the River Severn using a Random-Walk Model", Journal of Hydraulic Research, 31(3), pp. 323-Hinze, J.O., 1987. Turbulence. McGraw-Hill, New York. Jirka, G., Akar, J., 1991. "Hydrodynamic Classification of Submerged Multiport-Diffuser Discharges", Journal of Hydraulic Engineering, 117(9), pp. 1113-1128. Johannesson, H. , Parker, G., 1989. "Secondary Flow in Mildly Sinuous Channel", Journal of Hydraulic Engineering, 115(3), pp. 289-306. Kikkawa, H. , Ikeda, S., Kitagawa, A. , 1976. "Flow and Bed Topography in Curved Open Channels", Journal of the Hydraulics Divison ASCE, 102(HY9), pp. 1327-1342. Kottegoda, N.T., 1979. Stochastic Water Resources Technology. John Wiley & Sons, New York. Kraus, N . C , Lohrmann, A. , Cabrera, R., 1994. "New Acoustic Meter for Measuring 3D Laboratory Flows", Journal of Hydraulic Engineering, 120(3), pp. 406-412. Lam, D.C.L. , Murthy, C.R., Simpson, R.B., 1984. Effluent Transport and Diffusion Models for the Coastal Zone. Springer-Verlag, New York. Marks, B., 1996. "Initial Dilution of a Horizontal Jet in a Strong Current", Thesis submitted in partial fulfillment of the degree Master of Applied Science, University of British Columbia, Vancouver, Canada. Nezu, I., Nakagawa, H. , 1993. Turbulence in Open-Channel Flow. A . A . Balkema, Rotterdam. Nezu, I., Rodi, W., 1986. "Open-Channel Flow Measurements with a Laser Doppler Anemometer", Journal of Hydraulic Engineering, 112(5), pp. 335-355. Northwest Hydraulic Consultants, 1994. "Determination of Sediment Deposition Zones Fraser River Basin", Department of Environment Report No. 1994-32. Odgaard, A.J . , 1986. "Meander Flow Model. I: Development", Journal of Hydraulic Engineering, 112(3), pp. 1117-1136. Rozovskii, I.L., 1961. Flow of Water in Bends of Open Channels. The Israel Program for Scientific Translations, Jerusalem. Rutherford, J .C , 1994. River Mixing. Wiley, New York. Steffler, P .M. , Jin, Y.C. , 1992. "Predicting Flow in Curved Open Channels by Depth-Averaged Method", Jounral of Hydraulic Engineering, 119(1), pp. 109-124. Yeh, K.C. , Kennedy, J.F., 1993. "Moment Model of Nonuniform Channel-Bend Flow. LFixed Beds", Journal of Hydraulic Engineering, 119(7), pp. 776-795. Yotsukura, N . , Sayre, W.W., 1976. "Transverse mixing in natural channels", Water Resources Research, 12, pp. 695-704. Appendix A PATRAC USER MANUAL Introduction This manual is intended to provide guidance and information about P A T R A C and associated programs. Theory of the program is provided in the main body of the thesis. The main program is simple enough that it may be used most who are competent with the Windows operating environment. More advanced topics, such as reprogramming the model to a new location, may only be accomplished by someone familiar with Visual Basic. A working knowledge of P V - W A V E is also useful if the W A V E analysis tools are to be used. This program is provided as is for the purpose of academic study only, and the author takes no responsibility for the use/misuse of the software. The program is O N L Y FOR A C A D E M I C INTEREST, not for commercial/regulatory use. Program Structure The complete set of routines is written in two languages. Most of the simulation software is written in Visual Basic, while most of the graphics software is written in PV-W A V E . A l l of the components were used to produce the results in the thesis. The following table lists the components : Table A.l - PATRAC Software Components Software Component Language Purpose P A T R A C . F R M Visual Basic This is the main body of the program, and contains the following routines: - Initialization of the variables - Random walk engine - Production of output files (movie files) P A T R A C E . F R M Visual Basic This component is used to process the movie files into 2D histograms. These histograms may then be used to generate contour plots such as those produced by JCON.PRO E D G E . F R M Visual Basic This component is used to perform edge detection of the plume. The program uses movie files produced by PATRAC.FRM. The edge detection files are then used by P30EDGE.PRO to produce a steady-state image of the plume width. DILUTION.FRM Visual Basic This component is used to generate the centreline and bank dilutions which are then used by JDIL.PRO to plot dilution as a function of downstream distance. PATRAC.BAS Visual Basic Variables which are required by the Visual Basic routines are defined in this section. P A T R A C . M A K Visual Basic The main file which houses all of the above Visual Basic routines under one interface. P30BIRD.PRO PV-WAVE This program is used to generate the overhead view of the plume in Cartesian coordinates. P30EDGE.PRO PV-WAVE This program is used to average the results from EDGE.FRM into one file which can then be plotted to show plume width versus downstream distance. JCON.PRO PV-Wave This program is used to produce contour maps of dilution from the 2D histograms produced by P A T R A C E . F R M . JDIL.PRO PV-Wave This program is used to produce a plot of dilution versus downstream distance from the results produced by DILUTION.FRM. REMAP.PRO PV-Wave Remaps one movie file into a curvilinear system, producing a more realistic overhead view of the plume. Using PATRAC This section will describe the procedure for using PATRAC and the subprograms to produce results similar to those presented in the thesis. It is assumed that the model is being used to simulate the Northwood site. Modifying the programs to relocate the site should only be undertaken by an experienced programmer. Although most of the variables which need to be described may be done so from the main control screens, some programming changes may be needed as well. Where this is the case, the changes required are also outlined. The relevant section of code is included, and areas which need to be changed are bold faced. Part 1 - Running the main program 1) When in the Visual Basic environment, load P A T R A C . M A K . This will load a project which contains all of the Visual Basic components required. 2) Modify code as required. Specifically, code must be modified as follows: a) to modify the diffuser design, or to incorporate the results of a near field mixing model, the statements under "Diffuser Release Characteristics" must be changed. This section contains loop statements which specify the initial starting positions of each particle. The current routine assigns particles evenly to three locations, representing the three diffuser risers at Northwood. If northwoodoption.value=true then yi = -150 for counter=l to nop if yi >= -110 then yi = -150 x(counter)=1800:y(counter) = yi: z(counter) = initiald / 2 : status(counter) = 2: dprev(counter) = initiald yi = yi + 15 next counter end if b) change the DOS paths to the geometry files. The three files, BNORTH.WID, BNORTH.RAD, and BNORTH.CUR describe the width, radius of curvature, and direction of curvature of the site at 200 meter intervals. The 200 metre spacing is not an option which may be easily changed, as almost every simulation and graphing routine assumes 200 metre spacing. The user modifies (as well as uses) the program at his/her own risk! open "c:\users\jason\frap\patrac30\bnorth.wid" for input as #1 for counter = 0 to 63 input#l,bf bf = val(bf) bfrommap (counter) = bf next counter close #1 ...similarly for the radius of curvature and directions files... c) modify boundary conditions. If necessary, the boundary conditions may be modified under "CHECK B O U N D A R Y CONDITIONS". Modification of these statements is not recommended. 3) Start the program. You will be presented with the main screen. Many options are available here: PARAMETERS The following parameters and options must be set before using the main random walk engine : General Parameters : Flow Rate - of the river (m^/s) Manning n - Manning friction factor n # of Particles - The number of particles to be used. Maximum is 30000 (recommended) Slope - Slope of the channel Secondary Currents Either set to E N A B L E D or DISABLED. Provides a convenient way of shutting off the secondary currents routine. Release Options Either set to STRAIGHT C H A N N E L or NORTHWOOD SITE. If set to STRAIGHT C H A N N E L , the user may specify a starting X , Y , and Z location for the particles. If set to NORTHWOOD SITE, these starting locations are ignored. Release/Step - the total number of particles released per time step. Movies Either DISABLED or E N A B L E D . If E N A B L E D , the user specifies the interval at which output files are to be generated as well as a DOS path for their location. PROGRAM BUTTONS START - after describing all of the input parameters, this starts the random walk engine. HISTOGRAM - initiates the PATRACE routine for calculating 2D histograms. EDGE DETECT - initiates the EDGE routine for determining plume width. DILUTIONS - initiates the DILUTION routine for determining centreline and side dilutions 4) When all of the parameters have been set, press START. The program will begin to simulate, and status will be shown on the status line. Depending on the simulation, run-times of up to 6 hours may be expected on a Pentium 90. The movie files have the extension ".MOV". The program may be broken and then restarted usually without loss of results. The only way to end the program is to select E N D from the Visual Basic menu. Part 2 - Running PATRACE When the main body has finished running, you will be left with a set of files which contain the X , Y , and Z locations of each particle during each snapshot. P A T R A C E may be used to produce 2D histograms, which may then be used with JCON.PRO to produce a contour map of dilution. 1) From the main control screen, press HISTOGRAM. 2) Specify the path and filename of the movie files to be processed. 3) Specify the starting and ending file numbers. Note that a number of the files at the beginning of the simulation will not be of use as steady state conditions will not have developed. Similarly, some files at the end will also not be useful as steady state conditions will have ended. Use P30BIRD to examine the files to ensure that the images are suitable before processing them. 4) Press S C A N FILES. The program will search each file and determine the minimum and maximum X and Y values. These will serve as a guideline for choosing the parameters in the next step. 5) Specify the output specifications. The program requires that you specify the minimum and maximum X and Y values for the histogram, as well as the bin sizes. Use the results from step 4 as guide. 6) Press GO. The program will generate a series of 2D histograms which may then be used with JCON.PRO to produce a contour map of dilution. These files have the extension ".HIS". Part 3 - Running EDGE Another analysis option is to run EDGE, which determines the plume width and centreline of each movie file. The results may then be averaged with P30EDGE.PRO. 1) Modify code as required. a) The size and step of the box which is used to locate the plume edge are coded directly into this routine. The box is also programmed to only scan from the centreline of the channel to the right bank, as the Northwood plume never crossed the centreline during simulation runs. This increases the computational efficiency of the routine by a factor of 2. If simulations produce a plume which travels across the centreline of the river, then these values must be recoded. ...for i = 1 to particles if ((edgex(i) > (xi-100)) and (edgex(i) < (xi+100))) then slice(si) = edgey(i) si = si +1 end if ...for scany = 0 1 to -170 2 step - 1 J ng = 0 for i = 1 to si if slice (i) > scany then ng=ng+l next i ng = ng / si * 100 diff = abs (2.5 4 - ng) if diff < mindiff then mindiff = diff: lb = scany next scany ...for scany = 0 1 to -170 2 step -l-> ng = 0 for i = 1 to si if slice (i) > scany then ng=ng+l nexti ng = ng / si * 100 diff = abs (50 4 - ng) if diff < mindiff then mindiff = diff: cl = scany next scany ...for scany = 0 1 to -170 2 step - l 3 ng = 0 for i = 1 to si if slice (i) > scany then ng=ng+l nexti ng = ng / si * 100 diff = abs (97.5 4 - ng) if diff < mindiff then mindiff = diff: rb = scany next scany In the above, the superscripted numbers indicate the following : 0 - The width of the slice taken at each X location 1 - Maximum value of y to begin scan 2 - Minimum value of y to end scan 3 - Step size of y to scan 4 - Edge detection value (97.5-2.5)=95%, with 50 equaling 50% or the middle of the plume 2) From the main control screen, press "EDGE DETECT". 3) Specify the number of particles, the minimum and maximum values of X , and the X step size. Supply a path and filename, along with the starting and ending file numbers. Press START. The routine will produce one edge detection file for each movie file. The edge detection files have the extension ".TRA". Part 4 - Running D I L U T I O N 1) From the main control screen, press "DILUTIONS". 2) Specify the number of particles, a path and filename, and the starting and ending file numbers. You must also supply the name of the EDGE95 file, which is the output from P30EDGE (i.e. P30EDGE.PRO must be run before this routine may be used.) 3) Press START. The program will scan the movie files and produce a series of files with a ".DIL" extension. Each of these files contains a ID histogram along the centreline and right bank of the channel. This histogram is then converted into dilution with JCON.PRO. Part 5 - Running P30BIRD P30BIRD is a P V - W A V E program which is used to visualize the plume from overhead. This routine and the other P V - W A V E routines are V E R Y specific to the Northwood site. Extensive recoding would be required to move these routines to another location. 1) Modify the code as required : a) Specify the path, filename and file number under "SPECIFY THE FILE N U M B E R S " . for i = 14,14 do begin fnam='c:\users\jason\movies\runO\nor'+strcompress(i,/remove_all)+'.mov' b) Comment out the appropriate lines depending on whether a straight channel or the Northwood site was used in the simulation. These areas are clearly marked in the code. 2) Run the program at the Wave console by typing .RUN P30BIRD. A graphics screen will open showing the plume, along with the channel boundaries. Part 6 - Running P30EDGE P30EDGE takes the ".TRA" files produced by EDGE and averages them into one file which specifies an average of the plume. The file "EDGE95.CSV" contains the left bank, centerline, right bank, and width of the plume, all as functions of downstream distance. 1) Modify the code as required : a) specify the path, filename and file numbers of the files to be averaged, for i=20,40 do begin fnam='c:\users\jason\movies\runO\nor'+strcompress(i,/remove_all)+'.tra' b) specify the number of files being averaged. widavg=widavg/21.0 c) specify the path for the output file. status=dc_write_free('c:\users\jason\movies\run0\edge95.csv',x,lbavg, clavg,widavg,/column) 2) Run the program at the Wave console by typing .RUN P30EDGE. The program will run and produce a single ".csv" file as the output. This file can then be plotted by any standard graphing package. Part 7 - Running JCON JCON.PRO produces a contour map of dilution from the ".HIS" files produced by PATRACE. 1) Modify the code as required : a) specify the path, filename and file numbers of the ".HIS" files to be averaged. for i=20,40 do begin fnam='c:\uses\jason\movies\run3\nor'+strcompress(i,/remove_all)+'.his' b) specify the number of files to be averaged. his=his/21.0 c) specify the flow rate, slope, and Manning n that were used in the simulation. q=805. s= 0005 n=.03 d) specify the dimensions of the 2D arrays. for i=0,52 do begin his(*,i)=his(*,i)/d(i) endfor for i=0,52 do begin for j=0,34 do begin... e) specify the number of particles per time step and effluent flow rate, along with the bin sizes of the 2D array. dilution=(12.1/1.52)/(his/(5.3*25.4)) In the above, the superscripted numbers indicate the following : 1 - Number of particles per time step used 2 - Effluent flow rate (m3/s) 3 - Y Bin size 4 - X Bin size f) specify the contour levels. (Knowledge of the PV-Wave contour command is essential for this step. levels=[20,40,80,100,200] 2) Run the program at the Wave console by typing .RUN JCON. The program will average the files specified, open a graphics windows, and plot a contour map of dilution. Part 8 - Running JDIL JDIL.PRO produces a plot of centreline and bank dilution. 1) Modify the code as required : a) Specify the path, filename and file numbers of the ".DLL" files to average. for i=20,40 do begin fnam='c:\users\jason\movies\run3\nor'+strcompress(i,/remove_all)+'.dil' b) Specify the number of files being averaged. bdavg=bdavg/21.0 clavg=clavg/21.0 c) Specify the flow rate, slope and Manning n which were used in the simulation. q=805. s=.0005 n=.03 d) specify the number of particles per time step and effluent flow rate dilution=(12.1/1.52)/(bdavg/( 10.* 10.)) In the above, the superscripted numbers indicate the following : 1 - Number of particles per time step used 2 - Effluent flow rate (m3/s) 2) Run the program at the Wave console by typing .RUN JDIL. The program will average the files specified, open a graphics window, and plot centreline and bank dilution as a function of downstream distance. Part 9 - Running R E M A P REMAP.PRO takes one movie file and remaps it onto a curvilinear coordinate system, suitable for use at the Northwood site. The results are output to a graphics screen. 1) Modify the code as required : a) specify the path and filename of the movie file to remap. status=dc_read_free('c:\users\jason\movies\run4\nor71.mov',x,dl,y,d2,z, /column) 2) Run the program at the Wave console by typing .RUN R E M A P . The program will load the file and redisplay it on a curvilinear coordinate system. The Code The following pages contain printouts of the code in Visual Basic as well as PV-W A V E . The code is documented with comment lines which indicate the beginning of each major step. This should make identifying major sections of code fairly simple. The most likely reason to change the code would be to try new formulae for components such as secondary currents. Replacement of these equations should be a simple matter. The flow chart presented in chapter 4 will aid the reader in following the PATRAC code. The remaining modules are fairly simple, containing small loops and conditional statements which are relatively easy to understand. Appendix B PATRAC Listings PATRAC.FRM - This is the main body of the program which performs the variable initialization, the random walk engine, and the production of the output files. Sub Command l_Click () P A T R A C 3.5 - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ' THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y 1 T H E A U T H O R ASSUMES NO RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ' P R O G R A M ' Read the input values from the input screen and ' set up initial values statustext.Text = "Setting up arrays..." q = Val(qtext.Text) n = Val(ntext.Text) slope = Val(slopetext.Text) nop = Val(noptext.Text) initialx = Val(initialxtext.Text) initialy = Val(initialytext.Text) straightb = Val(straightbtext.Text) rpt = Val(rpttext.Text) movieinterval = Val(movieintervaltext.Text) moviename = movienametext.Text Randomize Timer t = 0 timestep = timesteptext.Text kappa = .41 nopground = 0 diam = diam /1000 framenumber = 0 movienumber = 0 hp = 0 ' Preload arrays with initial values If straightoption.Value = True Then b = straightb If northwoodoption.Value = True Then b = 340 initiald = (n * q / (b * slope A .5)) A (3 / 5) For counter = 1 To nop x(counter) = initialx: y(counter) - initialy: z(counter) = initiald / 2: status(counter) = 2: dprev(counter) = initiald Next counter If northwoodoption.Value = True Then 86 yi =-150 " For counter = 1 To nop Ifyi >= -110 Then yi = -150 x(counter) = 1800: y(counter) = yi: z(counter) = initiald / 2: status(counter) = 2: dprev(counter) = initiald yi = yi + 15 Next counter End If Open "c:\users\jason\frap\patrac30\bnorth.wid" For Input As #1 For counter = 0 To 63 Input #l,bf bf=Val(bf) bfrommap(counter) = bf Next counter Close #1 Open "c:\users\jason\frap\patrac30\bnorth.rad" For Input As #1 For counter = 0 To 63 Input #1, rf rf=Val(rf) rfrommap(counter) = rf Next counter Close #1 Open "c:\users\jason\frap\patrac30\bnorth.cur" For Input As #1 For counter = 0 To 63 Input #1, cf cfrommap(counter) = cf Next counter Close #1 'Begin particle tracking statustext.Text = "Particle tracking..." While nopground < nop 'Update screen display timetext.Text = t nopgroundtext.Text = nopground 'Release rpt number of particles If hp <= nop Then For c = hp To (hp + (rpt - 1)) status(c) = 0 Next c hp = hp + rpt End If For counter = 1 To nop If status(counter) = 0 Then ' Calculate width, depth, and velocity bi = Int(x(counter) / 200) xc - x(counter) - bi * 200 If northwoodoption.Value = True Then b = bfrommap(bi) + (bfrommap(bi + 1) - bfrommap(bi)) / 200 * xc Else b = straightb End If d = (n * q / (b * slope A .5)) A (3 / 5) ubar = q / (b * d) ' Calculate shear velocity ustar = Sqr(9.81 * d * slope) ' Determine u from the log velocity law (no wake functions) If z(counter) > d Then z(counter) = d - .01 u = ubar + ustar / kappa * (Log((d - z(counter)) / d) + 1) ' Determine if curved section exists then calculate secondary currents If secondary.Value = True Then If cfrommap(bi) = "s" Then vs = 0: ws = 0: wsc = 0: vsc •- 0 Else If cfrommap(bi) = "ccw" Then m = .41 * (ubar / ustar) vs = (ubar / .41 A 2) * (d / rfrommap(bi)) * (1.3 + 2.4 / m) ws = vs * (d / b) vsc = -vs * (2 * z(counter) / d. - 1) wsc = ws * (y(counter) / (b / 2)) Else If cfrommap(bi) = "cw" Then m = .41 * (ubar / ustar) vs = (ubar / .41 A 2) * (d / rfrommap(bi)) * (1.3 + 2.4 / m) ws - vs * (d / b) vsc - vs * (2 * z(counter) / d - 1) wsc = -ws * (y(counter) / (b / 2)) End If End If End If End If ' Run for Straight Channel Only If secondary .Value = True Then ' m = .41 * (ubar / ustar) vs = (ubar / .41 A 2) * (d / 450) * (1.3 + 2.4 / m) ' vs = .2 * ubar ws = vs * (d / b) vsc = -vs * (2 * z(counter) / d - 1) wsc = ws * (y(counter) / (b / 2)) End If ' Calculate fluctuations uprime = ustar * 2.3 * Exp(-(d - z(counter)) / d) vprime = ustar * 1.63 * Exp(-(d - z(counter)) / d) wprime = ustar * 1.27 * Exp(-(d - z(counter)) / d) rudl =Rnd(l) Ifrudl =0Then rudl = Rnd(l) rud2 = Rnd(l) If rud2 = 0 Then rud2 = Rnd(l) stnd = (-2 * Log(rudl)) A .5 * Cos(2 * 3.14 * rud2) uhf = uprime * stnd rudl = Rnd(l) Ifrudl = 0 Then rudl = Rnd(l) rud2 = Rnd(l) If rud2 = 0 Then rudl = Rnd(l) stnd = (-2 * Log(rudl)) A .5 * Cos(2 * 3.14 * rud2) vhf = vprime * stnd rudl = Rnd(l) Ifrudl =0 Then rudl = Rnd(l) rud2 = Rnd(l) If rud2 = 0 Then rudl = Rnd(l) stnd = (-2 * Log(rudl)) A .5 * Cos(2 * 3.14 * rud2) whf = wprime * stnd ufluc = u + uhf vfluc = vsc + vhf wfluc = wsc + whf ' Calculate deltax, deltay, and deltaz then move particles to new positions ' Including Stream Tube Corrections bi = Int(x(counter) / 200) xc = x(counter) - bi * 200 If (ufluc * timestep) < (200 - xc) Then deltax = ufluc * timestep x(counter) = x(counter) + deltax bi = Int(x(counter) / 200) xc = x(counter) - bi * 200 If northwoodoption.Value = True Then b = bfrommap(bi) + (bfrommap(bi + 1) - bfrommap(bi)) / 200 * xc Else b = straightb End If d = (n * q / (b * slope A .5)) A (3 / 5) If northwoodoption. Value = True Then If xc> 0 Then vcorrmax = ufluc * (((b / 2) - (bfrommap(bi) / 2)) / xc) vcorr - ((2 * vcorrmax) / b) * y(counter) vfluc = vfluc + vcorr End If deltay = vfluc * timestep y(counter) = y(counter) + deltay Else smallstep = ((200 - xc) / (ufluc * timestep)) deltax = ufluc * smallstep * timestep x(counter) = x(counter) + deltax - .01 bi = Int(x(counter) / 200) xc = x(counter) - bi * 200 If northwoodoption.Value = True Then b = bfrommap(bi) + (bfrommap(bi + 1) - bfrommap(bi)) / 200 * xc Else b = straightb End If d = (n * q / (b * slope A .5)) A (3 / 5) If northwoodoption.Value = True Then If xc> 0 Then vcorrmax = ufluc * (((b 12) - (bfrommap(bi) / 2)) / xc) vcorr = ((2 * vcorrmax) / b) * y(counter) vfluc = vfluc + vcorr End If deltay = vfluc * smallstep * timestep y(counter) = y(counter) + deltay deltax = ufluc * (1 - smallstep) * timestep x(counter) = x(counter) + deltax bi = Int(x(counter) / 200) xc = x(counter) - bi * 200 If northwoodoption.Value = True Then b = bfrommap(bi) + (bfrommap(bi + 1) - bfrommap(bi)) / 200 * xc Else b = straightb End If d = (n * q / (b * slope A .5)) A (3 / 5) If northwoodoption.Value = True Then If xc> 0 Then vcorrmax = ufluc * (((b 12)- (bfrommap(bi) / 2)) / xc) vcorr = ((2 * vcorrmax) / b) * y(counter) deltay = (vsc + vhf + vcorr) * (1 - smallstep) * timestep y(counter) = y(counter) + deltay End If End If ' Calculate streamtube correction velocity If northwoodoption.Value = True Then wcorr = (z(counter) / dprev(counter)) * ((d - dprev(counter)) / timestep) wfluc = wfluc + wcorr End If deltaz = wfluc * timestep z(counter) = z(counter) + deltaz Check boundary conditions If z(counter) > d Then z(counter) = d - Abs(z(counter) - d) If z(counter) < 0 Then z(counter) = 0 + Abs(z(counter) - 0) If z(counter) > d Then z(counter) = d - Abs(z(counter) - d) If z(counter) < 0 Then z(counter) = 0 + Abs(z(counter) - 0) If y(counter) < -b / 2 Then y(counter) = (-b 12) + Abs((y(counter) + (b / 2))) If y(counter) > b / 2 Then y(counter) = (b / 2) - Abs((y(counter) - (b / 2))) If x(counter) > 3200 Then status(counter) = 1: nopground = nopground + 1 If x(counter) < 0 Then x(counter) = 0 End If dprev(counter) = d Next counter t = t + timestep framenumber = framenumber + timestep ' Check if it is time to produce a movie file If (framenumber = movieinterval And movieoption.Value = True) Then framenumber = 0 moviefilename = moviename + Format(movienumber) + ".mov" movienumber = movienumber + 1 statustext.Text = "Last Movie File :" + moviefilename Open moviefilename For Output As #1 For moviecounter = 1 To nop Print #1, x(moviecounter);y(moviecounter);z(moviecounter) Next moviecounter Close #1 End If Wend End Sub Sub Command2_Click () patrac.Visible = False gradatio.Visible = True End Sub Sub Command3_Click () patrac.Visible = False edge.Visible = True 'edgelmovienametext.Text = patrac!movienametext.Text 'edgelfinishframetext.Text = movienumber - 1 End Sub Sub Command4_Click () patrace.Visible = True patrac.Visible = False patracelmovienametext.Text = patraclmovienametext.Text patracelfinishframetext.Text = movienumber - 1 End Sub Sub Command5_Click () patrac.Visible = False form 1. Visible = True End Sub Sub gotoanalysis_Click () End Sub Sub gotofileoptions_Click () End Sub Sub Timerl_Timer () If proi = 1 Then statustext.Text If proi = 2 Then statustext.Text If proi = 3 Then statustext.Text If proi = 4 Then statustext.Text If proi = 5 Then proi = 0 proi = proi + 1 End Sub "PATRAC 3.5 (C) 1996 Jason Vine" "No Unauthorized Use without Author Consent" "White Indicates User Entry Field" "Green Indicates Program Output" PATRACE.FRM - This component is used to process the movie files into 2D histograms. These histograms may then be used to generate contour plots such as those produced by JCON.PRO Sub Command2_Click () ' P A T R A C E 3.5 - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ' THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ' T H E A U T H O R ASSUMES N O RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ' P R O G R A M 'Read the input values from the input screen and setup 'initial values moviename = movienametext.Text startframe = Val(startframetext.Text) finishframe = Val(finishframetext.Text) minx = 100000: maxx = -100000: miny = 100000: maxy = -100000 particles = Val(particlestext.Text) 'Begin scanning the files for extreme values statustext.Text = "Scanning Files for Extreme Values..." For frame = startframe To finishframe 'Open the relevant movie file statustext.Text = "Scanning File # " + frame moviefilename = moviename + Format(frame) + ".mov" Open moviefilename For Input As #1 'Scan the file for the maximum values For counter = 1 To particles Input #1, xin, yin, zin If xin < minx Then minx = xin If xin > maxx Then maxx = xin If yin < miny Then miny = yin If yin > maxy Then maxy = yin Next counter Close #1 Next frame 'Display the maximum values on the screen minxtext.Text = minx maxxtext.Text = maxx minytext.Text = miny maxytext.Text = maxy statustext.Text = "Extremes Processed..." 93 End Sub Sub Command3_Click () 'Read the input variables from the input screen 'and setup initial values xmin = Val(minxouttext.Text) xmax = Val(maxxouttext.Text) xbin = Val(xbinouttext.Text) ymin = Val(minyouttext.Text) ymax = Val(maxyouttext.Text) ybin = Val(ybinouttext.Text) moviename = movienametext.Text startframe = Val(startframetext.Text) finishframe = Val(finishframetext.Text) particles = Val(particlestext.Text) For frame = startframe To finishframe statustext.Text = "Processing File #" + frame ' Blank out the Histo array For counter = 0 To 250 For subcounter = -100 To 100 histo(counter, subcounter) = 0 Next subcounter Next counter ' Load Histo array with y values index = Int(ymin / ybin) yindexmin = index For counter - ymin To ymax Step ybin histo(0, index) = counter index = index + 1 Next counter yindexmax = index ' Load Histo array with x values index = 1 For counter = xmin To xmax Step xbin histo(index, yindexmin - 1) = counter index = index + 1 Next counter xindexmax = index ' Open the relevant file moviefilename = moviename + Format(frame) + ".mov" Open moviefilename For Input As #1 ' Scan the file and begin counting For counter = 1 To particles Input #1, xin, yin, zin xbinno = Int((xin - minx) / xbin) + 1 ybinno = Int(yin / ybin) histo(xbinno, ybinno) = histo(xbinno, ybinno) + 1 Next counter Close #1 ' Write the Histo array to disk histofilename = moviename + Format(frame) + ".his" Open histofilename For Output As #1 For counter = 1 To (xindexmax - 1) For subcounter = (yindexmin) To (yindexmax - 1) Print #1, Format(histo(counter, subcounter));"," Next subcounter Print #1, Next counter Close #1 Next frame statustext.Text = "Finished..." End Sub EDGE.FRM - This component is used to perform edge detection on the plume. The program uses movie files produced by PATRAC.FRM. The edge detection files are then used by P30EDGE.PRO to produce a steady-state image of the plume width. Sub Command2_Click () ' E D G E 3.5 - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ' THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ' T H E A U T H O R ASSUMES N O RESPONSIBILITY FOR T H E USE/MISUSE O F THIS ' P R O G R A M ' Read the input values from the input screen and ' set up initial values moviename - movienametext.Text particles = Val(particlestext.Text) startframe = Val(startframetext.Text) finishframe = Val(finishframetext.Text) ep = eptext.Text minx = minxtext.Text maxx = maxxtext.Text xstep = xsteptext.Text 'Load the relevant file For frame = startframe To finishframe statusbox.Print "Loading File #"; frame moviefilename = moviename + Format(frame) + ".mov" Open moviefilename For Input As #1 For counter = 1 To particles Input #1, edgex(counter), edgey(counter), edgez(counter) Next counter Close #1 'Slice at given x intervals edgefilename = moviename + Format(frame) + ".tra" Open edgefilename For Output As #1 For xi = minx To maxx Step xstep statusbox.Print "Slicing at x="; xi; "m" For si = 1 To 3000 slice(si) = 0 Next si 'Count the number of particles in the slice si= 1 For i = 1 To particles If ((edgex(i) > (xi - 10)) And (edgex(i) < (xi + 10))) Then slice(si) = edgey(i) si = si + 1 End If Next i cl = 0 mindiff = 100 96 'find left plume edge For scany = 0 To -170 Step -1 ng = 0 For i = 1 To si If slice(i) > scany Then ng = ng + 1 Next i ng = ng / si * 100 diff=Abs(5-ng) If diff < mindiff Then mindiff = diff: lb = scany Next scany 'find centreline mindiff =100 For scany = 0 To -170 Step -1 ng = 0 For i = 1 To si If slice(i) > scany Then ng = ng + 1 Next i ng = ng / si * 100 diff=Abs(50-ng) If diff < mindiff Then mindiff = diff: cl = scany Next scany 'find right plume edge mindiff =100 For scany = 0 To-170 Step-1 ng = 0 For i = 1 To si If slice(i) > scany Then ng = ng + 1 Next i ng = ng / si * 100 diff=Abs(95-ng) If diff < mindiff Then mindiff = diff: rb = scany Next scany Print #l,xi;",";lb;","; cl; ",";rb Next xi statusbox.Cls Close #1 Next frame End Sub DILUTION.FRM - This component is used to generate the centreline and bank dilutions which are then used by JDIL.PRO to plot dilution as a function of downstream distance. Sub dilution_Click () ' DILUTION - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ' THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ' T H E A U T H O R ASSUMES NO RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ' P R O G R A M ' Read the input values from the input screen and ' set up initial values moviename = movienametext.Text particles = Val(particlestext.Text) startframe = Val(startframetext.Text) finishframe = Val(finishframetext.Text) edgename = edgenametext.Text ' Load the width file used with this routine Open "c:\users\jason\frap\patrac35\dilution.csv" For Input As #1 For counter = 1 To 13 Input #1, dilb(counter) Next counter Close #1 nib = 0 ncl = 0 'Load the relevant movie file For frame = startframe To finishframe statusbox.Print "Loading File #"; frame moviefilename = moviename + Format(frame) + ".mov" Open moviefilename For Input As #1 For counter = 1 To particles Input #1, edgex(counter), edgey(counter), edgez(counter) Next counter Close #1 movieedgename = moviename + Format(frame) + ".tra" 'Load the E D G E file Open movieedgename For Input As #1 For counter = 1 To 13 Input #1, xtra(counter), lbtra(counter), cltra(counter), rbtra(counter) Next counter Close #1 ' Open the output file dilfilename = moviename + Format(frame) + ".dil" Open dilfilename For Output As #2 ' Scan the movie file For scani = 1 To 13 scanx = 1900 + (scani - 1) * 100 ' Count the number of particles on the bank and centreline For counter = 1 To particles If (edgex(counter) > (scanx - 5) And edgex(counter) < (scanx + 5) And edgey(counter) > ((-dilb(scani) / 2# + 10) - 5) And edgey(counter) < ((-dilb(scani) / 2# + 10) + 5)) Then nib = nib + 1 If (edgex(counter) > (scanx - 5) And edgex(counter) < (scanx + 5) And edgey(counter) > (cltra(scani) - 5) And edgey(counter) < (cltra(scani) + 5)) Then ncl = ncl + 1 Next counter statusbox.Print; scanx; nib; ncl ' Output the results to the file Print #2, scanx;","; nib;","; ncl nib = 0 ncl = 0 Next scani Close #2 statusbox.Cls Next frame statusbox.Print "Finished..." End Sub 99 P30BIRD.PRO - This program is used to generate the overhead view of the plume in Cartesian coordinates. ; P30BIRD - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ; THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ; T H E A U T H O R ASSUMES N O RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ; P R O G R A M ; Set up initial parameters set_plot,"win32" ;set_plot,"ps" ;device,filename="c:\users\jason\wave 1 .ps",/landscape !p.multi=[0,0,0,0,0] !x.style=l !x.range=[1800,3100] !x.ticks=13 !y.range=[-200,200] !x.tickformat='(F5.0)' Ix.title-Distance Downstream (m)' !y.title='Y (m)' x=dblarr(30000) y=dblarr(30000) z=dblarr(30000) xp=dblarr(13) lb=dblarr(13) cl=dblarr(13) rb=dblarr(13) wid=dblarr(13) ;mov=intarr(500,250,50) i=0 ;LOADS T H E WIDTH FILE FOR N O R T H W O O D INTO A N A R R A Y SO T H A T ;IT C A N O V E R P L O T T E D O N TOP OF T H E P L U M E bx=findgen(64)*200 by=dblarr(64) status=dc_read_free('c:\users\jason\frap\patrac30\bnorth.wid',by) ;mov(*,*,0)=tvrd(0,0,500,250,order=l) ;SPECIFY T H E FILE N U M B E R S (MULTIPLE FOR M O V I E MAKING) for i=14,14 do begin fnam='c:\users\jason\movies\runO\nor'+strcompress(i,/remove_all)+'.mov' status=dc_read_free(fnam,x,d 1 ,y ,d2,z,/column) ; mov(*,*,i)=tvrd(0,0,500,250,order=l) plot,x,y,psym=3 ; A C T I V A T E T H E S E LINES W H E N USING N O R T H W O O D A N D FRASER RIVER ; oplot,bx,by/2 ; oplot,bx,-by/2 ; A C T I V A T E T H E S E LINES W H E N USING A STRAIGHT C H A N N E L oplot, [ 1800,3300], [-150,-150] oplot,[1800,3300],[150,150] ; A C T I V A T E T H E S E LINES W H E N USING N O R T H W O O D T O P L O T T H E ;95 A N D 99 % LINES ; status=dc_read_free('c:\users\jasori\movies\run2\edge99.csv',xp,wid,/column) ; patb=[320.,300.,307.5,315.,327.5,340.,320.,300.,237.5,175.,180.,185.,185.] ; oplot,xp,-((patb/2.0)-wid) ; status=dc_read_free('c:\users\jason\movies\run2\edge95.csv',xp,wid,/column) ; patb=[320.,300.,307.5,315.,327.5,340.,320.,300.,237.5,175.,180.,185.,185.] ; oplot,xp,-((patb/2.0)-wid) j A C T I V A T E T H E S E LINES W H E N USING A STRAIGHT C H A N N E L T O P L O T ;THE P L U M E E D G E (EITHER 95 OR 99%) ; status=dc_read_free('c:\users\jason\movies\run3\edge95.csv',xp,lb,cl,rb,wid,/column) ; oplot,xp,lb ; oplot,xp,rb endfor ;device,/close end P30EDGE.PRO - This program is used to average the results from EDGE.FRM into one file which can then be plotted to show plume width versis downstream distance. ; P30EDGE - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ; THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ; T H E A U T H O R ASSUMES NO RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ; P R O G R A M ;Setup initial parameters x=dblarr(13) cl=dblarr(13) rb=dblarr(13) lb=dblarr(13) rbavg=dblarr(13) lbavg=dblarr(13) wid=dblarr(13) widavg=dblarr( 13) clavg=dblarr(13) ;Specify the file numbers and file path/name for i=20,40 do begin fnam='c:\users\jason\movies\runO\nor'+strcompress(i>/remove_all)+'.tra' ; A C T I V A T E THIS LINE W H E N USING T H E U N B O U N D E D ROUTINES TN P A T R A C status=dc_read_free(fnam,x,dl,lb,d2,cl,d3,rb,/column) ; A C T I V A T E THIS LINE W H E N USING T H E B O U N D E D ROUTINES IN P A T R A C ; status=dc_read_free(fnam,x,d 1 ,wid,/column) ;Calculate the average value of the quantities of interest wid=abs(rb-lb) rbavg=rbavg+rb lbavg=lbavg+lb clavg=clavg+cl widavg=widavg+wid print,i endfor widavg=widavg/21.0 ; A C T I V A T E T H E S E LINES W H E N USING T H E U N B O U N D E D ROUTINES IN P A T R A C clavg=clavg/21.0 rbavg=rbavg/21.0 lbavg=lbavg/21.0 status=dc_write_free('c:\users\jason\movies\run0\edge95.csv',x,lbavg,clavg,rbavg,widavg,/column) ; A C T I V A T E THIS LINE W H E N USING T H E B O U N D E D ROUTINE IN P A T R A C ;status=dc_write_free('c:\users\jason\movies\ded\run2\edge95.csv',x,widavg,/column) end JCON.PRO -This program is used to produce contour maps of dilution from the 2D histograms produced by PATRACE.FRM. ; JCON - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ; THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ; T H E A U T H O R ASSUMES NO RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ; P R O G R A M ;Setup initial parameters set_plot,'win32' set_plot,'ps' device,filename='c:\users\jason\wave.ps',/landscape ;IGNORE THIS SECTION (OLD C O D E - M A Y H A V E S O M E USE) ;tek_color ;x=dblarr(30000) ;y=dblarr(30000) ;z=dblarr(30000) ;status=dc_read_free('c:\users\jason\movies\run3\nor30.mov',x,dl.y,d2,z,/column) ;plot,x,y,psym=3,xrange=[l 800,3100],yrange=[-170,0],$ ;xsty le= 1 ,ystyle= 1 ,color=4,$ ;xtitle='Distance Downstream (m)',ytitle='Y (m)' ;Load the relevant files and calculate the average value his=dblarr(35,53) temp=dblarr(35,53) for i=20,40 do begin fnam='c:\users\jason\movies\mn3\nor'+strcompress(i,/remove_all)V.his' status=dc_read_free(fnam,temp) his=his+temp print.i endfor his=his/21.0 x=1800+findgen(53)*25 y=0+findgen(35)*(-5.) b=dblarr(53) status=dc_read_free('c:\users\jason\frap\patrac35\contour.csv',b,/column) ;Turn the histograms into dilution q=805. s=.0005 n=.03 d=((n*q)/(b*sA.5))A.6 for i=0,52 do begin his(*,i)=his(*,i)/d(i) endfor for i=0,52 do begin for j=0,34 do begin if his(j,D eq 0 then his(j,i)=.0000000001 endfor endfor dilution=( 12./1.5)/(his/(5. *25.)) ;shade_surf,dilution,y,x,shades=bytscl(dilution,min=0,max=400),ax=90,az=0 dilplt=rotate(dilution, 1) dilplt=reverse(dilplt) dilplt=reverse(dilplt,2) ;Plot the results plot,x,-b/2,psym=-3,thick=3.0,$ xrange=[ 1800,3100] ,yrange=[-170,0] ,ysty le= 1 ,xstyle= 1 contour,dilplt,x,y,/noerase,$ ;levels=[300,400,500,1000] ,$ ;c_labels=[ 1,1,1,1,1 ],/follow,$ ;levels=[20,40,80,100,200],$ ;c_labels=[ 1,1,1,1,1,1 ],/follow,$ levels=[20,40,60,80,100,150,200],$ cJabels=[l,l,l,l,l,l,l],/follow,$ yrange=[-170,0] ,xrange=[ 1800,3100] ,xstyle= 1 ,ystyle= 1 ,$ xtitle='Distance Downstream (m)',ytitle='Y (m)' xcl=dblarr(13) lb=dblarr(13) cl=dblarr(13) rb=dblarr(13) wid=dblarr(13) status=dc_read_free('c:\users\jason\movies\run3\edge95.csv',xcl,lb,cl,rb,wid,/column) oplot,xcl,cl device,/close end JDIL.PRO - This program is used to produce a plot of dilution versus downstream distance from the results produced by DILUTION.FRM. » JDIL - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ; THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ; T H E A U T H O R ASSUMES NO RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ; P R O G R A M ;Setup initial parameters set_plot,'win32' ;set_plot,'ps' ;device,filename='c:\users\jason\wave.ps',/landscape x=dblarr(13) bd=dblarr(13) bdavg=dblarr(13) cld=dblarr(13) cldavg=dblarr(13) b=[320.,300.,307.5,315.,327.5,340.,320.,300.,237.5,175.,180.,185.,185.] ;Load the relevant files and calculate the average value for i=20,40 do begin fnam='c:\users\jason\movies\run3\nor'+strcompress(i,/remove_all)+'.dil' status=dc_read_free(fnam,x,d 1 ,bd,d2,cld,/column) bdavg=bdavg+bd cldavg=cldavg+cld print,i endfor ;Convert the average value into dilution on the centreline and bank bdavg=bdavg/21.0 cldavg=cldavg/21.0 for i=0,12 do begin if bdavg(i) eq 0 then bdavg(i)=.00000001 if cldavg(i) eq 0 then cldavg(i) =.00000001 endfor q=805. s=.0005 n=.03 d=((n*q)/(b*sA.5))A.6 bdavg=bdavg/d cldavg=cldavg/d dilutionb=( 12./1.5)/(bdavg/( 10. * 10.)) dilutioncl=( 12./1.5)/(cldavg/( 10. * 10.)) ;Plot the results 105 plot,x,dilutionb,xrange=[ 1900,3100] ,y range=[0,400] ,xsty le= 1 ,ysty le= 1 ,$ xtitle='Distance Downstream (m)',ytitle='Dilution',psym=-4 oplot,x,dilutioncl,psym=-2 oplot, [2800,2850] ,[40,40] ,psym=-4 oplot, [2800,2850] ,[30,30] ,psym=-2 xyouts,2900,39,'Bank' xyouts,2900,29,'C.L.' ;device,/close end REMAP.PRO -Remaps one movie file into a curvilinear system, producing a more realistic overhead view of the plume. ; R E M A P - JUNE 25/96 - ORIGINAL P R O G R A M B Y JASON VINE ; THIS P R O G R A M IS FOR A C A D E M I C INTEREST O N L Y ; T H E A U T H O R ASSUMES NO RESPONSIBILITY FOR T H E USE/MISUSE OF THIS ; P R O G R A M ;Set up initial parameters set_plot,'win32' set_plot,'ps' device,filename='c:\users\jason\wave.ps',/inches,ysize=7.0,xsize=7.0,yoffset=0.75 x=dblarr(30000) y=dblarr(30000) z=dblarr(30000) xp=dblarr(30000) yp=dblarr(30000) ;Load the relevant file status=dc_read_free('c:\users\jason\movies\run4\nor71.mov',x,dl,y,d2,z,/column) ;Transform the coordinate system x=x-1800 for i=0,29999 do begin if x(i) gt 1400 then x(i)=1400 if x(i) It 200 then xp(i)=x(i) else xp(i)=200+(450-y(i))*sin((x(i)-200)/450) if x(i) It 200 then yp(i)=y(i) else yp(i)=450-(450-y(i))*cos((x(i)-200)/450) endfor ;window,0,xsize=500,ysize=500 plot,xp,yp,psy m=3 ,xrange= [0,1000] bxp=dblarr(15) byp=dblarr(15) bx^fO.,100.,200.,300.,400.,500.,600.,700.,800.,900.,1000.4100.,1200.,1300.,1400.] b=[340.,320.,300.,307.5,315.,327.5,340.,320.,300.,237.5,175.,180.,185.,185.,185.] lb=b/2.0 rb=-b/2.0 for i=0,14 do begin if bx(i) It 200 then bxp(i)=bx(i) else bxp(i)=200+(450-lb(i))*sin((bx(i)-200)/450) if bx(i) It 200 then byp(i)=lb(i) else byp(i)=450-(450-lb(i))*cos((bx(i)-200)/450) endfor oplot,bxp,byp,psym=-3,thick=3.0 for i=0,14 do begin if bx(i) It 200 then bxp(i)=bx(i) else bxp(i)=200+(450-rb(i))*sin((bx(i)-200)/450) ifbx(i) lt 200 then byp(i)=rb(i) else byp(i)=450-(450-rb(i))*cos((bx(i)-200)/450) endfor oplot,bxp,byp,psy m=-3 ,thick=3.0 device,/close end
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Random walk model applicable to rivers
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Random walk model applicable to rivers Vine, Jason 1996
pdf
Page Metadata
Item Metadata
Title | Random walk model applicable to rivers |
Creator |
Vine, Jason |
Date Issued | 1996 |
Description | The disposal of waste effluent from industrial and municipal sources into rivers for the purpose of dilution is a common practice. Dilution occurs in the near-field due to specifics of the diffuser design, and in the mid-field due to parameters of the ambient environment. The ability to predict downstream dilution is of interest to regulatory officials and engineers who wish to design such disposal systems. A random walk model has been developed which is geared towards use in rivers. The model only performs mid-field calculations, however the potential exists to accept the results from a near-field model as input. The model attempts to account for the major velocity structures that occur in an open channel, including transverse and vertical shear, secondary currents due to flow through a bend, and turbulent fluctuations. Equations describing such phenomena are obtained from the literature. Laboratory experiments are conducted to verify the validity of some of these equations. As expected from theory, the research confirms that random walk modelling predicts the same results as Fickian diffusion. Arguments are presented as to the choice of a suitable timestep for use with a random walk model. Several key parameters of the random walk are identified, and their effect on the stability and accuracy of the model results are examined. The model is applied to the disposal of effluent from a pulp and paper mill located in Prince George, B.C. Results agree acceptably with a popular diffusion model as well as with digital imagery of the plume. Plots of plume width and dilution are produced for various river flow rates which may be expected at Prince George. |
Extent | 4676124 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-14 |
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.0050372 |
URI | http://hdl.handle.net/2429/6045 |
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 | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-0500.pdf [ 4.46MB ]
- Metadata
- JSON: 831-1.0050372.json
- JSON-LD: 831-1.0050372-ld.json
- RDF/XML (Pretty): 831-1.0050372-rdf.xml
- RDF/JSON: 831-1.0050372-rdf.json
- Turtle: 831-1.0050372-turtle.txt
- N-Triples: 831-1.0050372-rdf-ntriples.txt
- Original Record: 831-1.0050372-source.json
- Full Text
- 831-1.0050372-fulltext.txt
- Citation
- 831-1.0050372.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050372/manifest