ANALOG COMPUTING TECHNIQUES APPLIED TO ATMOSPHERIC DIFFUSION: CONTINUOUS POINT SOURCE by Fred. V. Brock E. Wendell Hewson Technical Report No, 3 ORA Project 03632 NATIONAL SCIENCE FOUNDATION GRANT G-11404 College of Engineering Department of Engineering Mechanics Meteorological Laboratories Institute of Science & Technology Great Lakes Research Division Special Report No. 12 THE UNIVERSITY OF MICHIGAN ANN ARBOR, MICHIGAN July 1961 ACKNOWLEDGMENT The authors wish to acknowledge the wholehearted cooperation and assistance of Mrs. Anne. Co Rivette, TABLE OF CONTENTS Page LIST OF ILLUSTRATIONS v ABSTRACT ix 1. INTRODUCTION 1 2. FINITE DIFFERENCE MODEL 5 3. ANALOG MODELS 4. COMPUTATION AND RESULTS 10 A. Isotropic Diffusion with an 84-Cell Model 10 B. Isotropic Diffusion with a 56-Cell Model 15 C. Nonisotropic Diffusion 19 D. Eddy Diffusivity Related to Plume Dimensions 22 E. Diffusion in Transitional States 23 F. Wind Direction Shear 32 G. Large-Scale Eddies 41 H. A Possible Combination of Two of the Models 50 5. SUMMARY AND CONCLUSIONS 51 APPENDIX 55 REFERENCES 56 iii LIST OF ILLUSTRATIONS Table Page I Cell Boundaries in the Horizontal and the Vertical for the First Quadrant of the 84-Cell Model 12 A-I Comparison of Plume Concentrations for the Meandering Plume and the Straight Plume 55 Figure 1 The 84-cell, finite difference model of isotropic diffusion from a point source. 6 2 The passive network analog of diffusion from a point source. 9 3 Computing circuit for the first quadrant of the 84-cell model. 15 4 Typical cell from the circuit of Fig. 3 showing the connections to the integrating amplifier. 14 5 Concentration plots for each cell as a function of distance downwind from a point source using the 84-cell model in isotropic conditions. 16 6 Concentration plots for each cell as a function of distance downwind from a point source using the 84-cell model in isotropic conditions. 17 7 Computing circuit for the first quadrant of the 56-cell model shown in (a)o An integrating amplifier for a typical cell is shown in (b). 18 8 Comparison of concentration plots for the source cell using the 84- and the 56-cell models with isotropic diffusion conditions. 20 9 Comparison of concentration plots for some cells using the 56-cell modelo 21 10 Computing circuit for the first quadrant of the 56-cell model which permits the inclusion of special forms of the coefficient of eddy diffusivity. 25 v LIST OF ILLUSTRATIONS (Continued) Figure Page 11 Computing circuit for the generation of the special forms of the coefficient of eddy diffusivity. 26 12 Plot of the modifying parameters of the coefficient of eddy diffuvisity. 28 13 Plot of the concentrations in the source cell as a function of the modifying parameter used. 29 14 Plot of the concentration in cell 1, 2 as a function of the modifying parameter used. 50 15 Plot of the concentration in cell 1, 3 as a function of the modifying parameter used. 51 16 The 25-cell, finite difference model of diffusion from a point source. 55 17 A schematic of the wind shear component V as a function of height in the 25-cell finite difference model. 54 18 Computing circuits for the 25-cell model for the case with wind shear. 36 19 Cross section of the plume at X = 1 showing lines of constant concentration in arbitrary units. 57 20 Comparison of the concentration plots for cells to the right and the left of the source cell in the presence of wind shear. 58 21 Concentration plots for some cells in the 25-cell model in the presence of wind shear. 59 22 Concentration plots for some cells in the 25-cell model in the presence of wind shear. 40 23 The 15-cell finite difference model of diffusion from a point source for use in the case of diffusion in the presence of large scale eddies. 42 24 Computing circuit for the 15-cell model including provision for variation of the horizontal wind component. 44 vi LIST OF ILLUSTRATIONS (Concluded) Figure Page 25 Computing circuit for the generation of the sinusoidal horizontal wind component. 45 26 Typical cell from the circuit of Fig. 23. 46 27 Plot of the concentrations for cells 0,2, 0,3, and 0,4 in the presence of large-scale eddies. 47 28 Plot of the concentration for cells -1,53 and +1,3 in the presence of large-scale eddies. 48 29 Plot of the concentrations for cells -1,2; +1,2; -1,4; and +1,4 in the presence of large-scale eddies. 49 30 Schematic of the fluctuating plume model showing the sinusoidal nature of the plume axis. 54 vii aR ABSTRACT The analog computer has been used to obtain solutions to the diffusion equation for the case of a continuous point source aloft. The basic model used is that of a source located between an inversion and the ground where the wind speed and eddy diffusivity are constant over the diffusing region. Solutions were also obtained with several special features added to the basic model: an increase of lateral diffusivity over the vertical; representation of the eddy diffusivity as a function of plume dimensions or as a function of distance downstream; presence of wind direction shear in the diffusing layers; and the presence of large-scale eddies which cause the meandering of the plume. / The solutions are presented as deviations from the reference solution, that solution derived from the basic model, to show graphically the effects of each parameter change. ix 1. INTRODUCTION Turbulent diffusion from an elevated point source is of primary concern in air-pollution meteorology. Although such diffusion is the central component of a common problem, it is by no means a well-understood phenomenon. This results in part from the lack of a sufficiently general mathematical model with a complete set of solutionso No matter how well diffusion is documented in the field, it cannot be completely and compactly expressed and understood except in reference to such a model. But atmospheric diffusion is such an intricate mechanism that no one model is ever likely to satisfy all damands. One model which is perhaps sufficiently general is the parabolic diffusion equation dx Kx~ ] +x fi 7 C x] where x.- aX ax - x ax dt +x y d dt and X = average concentration x,yz = space coordinates uvw = average wind components Ki = coefficient of eddy diffusivity in the i-direction t = time 1 This model has two main shortcomings: it yields only average concentration, and a complete set of solutions has not been found. In fact, such analytical solutions as are available apply only to special cases of limited application+1'2 Another deficiency is that the model implies an infinite diffusing velocity, but in cases where this is important, the hyperbolic diffusion equation can be used. In the case of steady-state diffusion, the velocity of diffusion is not important. The validity of the application of the parabolic diffusion equation to atmospheric diffusion has been questioned, but the lack of adequate solutions has made verification difficult. The use of electronic computing machines makes solutions available for any type of parameter variation but here another difficulty arises. Since analytical-solutions are not readily available, the functional dependence of some of the parameters has not been determined. For example Richardson has suggested that the coefficient of eddy diffusivity be a function of plume size, but the exact form of this function has not been indicated. The purpose of this paper is to show how the diffusion equation may be implemented on the electronic analog computer and to show the qualitative nature of the solutions. Means for obtaining improved solution accuracy are also indicateds The reasons for using the analog computer in this application are that it'is easy to program and that it presents the solution in graphical form. The problems chosen for solution-are based upon the case of diffusion in the steady-state from a point source located half-way between the surface 2 of the ground and the base of an inversion. The special cases considered here were selected from the many problems of interest described elsewhere. The methods outlined are perfectly general, however, and with an analog computer of sufficient capacity, a similar analysis may be made for nearly any meteorological conditions Added to the basic model are several functions of the average wind components and the coefficients of eddy diffusivity so as to consider increased lateral over vertical diffusivity, the case of a plume passing from one turbulent regime to another, the theoretical question of making the diffusivity a function of cloud size, the case of wind direction shear with height in the diffusing region, and the case of the meandering plume. The model is arranged so that the origin is at the source, the x-axis oriented downwind, the y-axis crosswind, and the z-axis vertical The basic model includes no provision for variations in the wind speed with height or in the eddy diffusivity in the y- or z-directions. This can be readily provided9 but the emphasis here is on special techniques. 2. FINITE DIFFERENCE MODEL The model presented in Eq. (1) can be simplified by deleting the terms that do not apply, to u ='K v (2) u a = y ^2 + Kz 2 " v ay (2) where the diffusion term in x was neglected in comparison to the transport term in x. The boundary conditions are 3 X -D co as x,yz 0oo LH H Kz > o as z - + for x > where H is the height above ground of the base of the inversion layer. These conditions state that the concentration is infinite at the source and that the ground and the inversion are impervious to the diffusing substance. The continuity condition is / 0 H H ii u x(x ylz) dz dy = (5). -00 2 for x > 0 where Q is the rate of emission from the source. It is convenient to state (2) in nondimensional form: Cs K a2s a2s as = = - - (4) ax Kz aY2 az2 )y where S = X/Xo Y = y/H Z = z/H x = H2 u x/K V = H V/Kz and the symbol X is reserved for later use. The constant Xo will be defined later. 4 The mathematical model to be used,. Eqo (4), must be put in a form suitable for programming on the analog computer. Since the computer can integrate only with respect to one independent variable, computer time, two of the independent variables must be replaced with finite differenceso The derivatives in Y' and. Z will be expressed as finite differences as follows: 2S 1 sm+l - Sm S m -l7 - Sm-1 2 (AY)m (AY) (AY)Im ml -Y 2(AY)mrsn+l - 2Sn + Sn-l Z i (AZ)n2 $n in The expressions for the second derivatives in Y and Z are different- because, in all the cases studied here, uniform differences are taken in Z but nonuniform differences are taken in. Y o If uniform differences were taken both ways, the two expressions would be identical. Equation (4) may now be stated as a set simultaneous, ordinary differential equations: dS Ky 1 S S.n n Sn- Smn d Kz (AY)m (AY)m+I (AY)ml 2 2 + - mjA j - ml n nm.n (6) (AZ)n2 2(AY)m where the subscripts m and n refer to discrete values in Y and Z respectively. The physical model corresponding to Eq. (6) is shown in Fig. 1. This is the particular model that will be used in the section on computation in connection with the isotropic diffusion caseo It consists of an 5 z ~ H A INVERSION.500 In N I - I.428 II I I I.357 4 -.286 i- - I.214 II 4.... I.143 I I I 4 4 0714... I I Xo Xo I II _ 01 xo xo.0714.1934 y H II c X,,, I, 1 i i Ti I I GROUND I 2 3 Fig. 1. The 84-cell, finite difference model of isotropic diffusion from a point source. 6 array of cells of rectangular cross section with width AY, height AZ, and with one end at x = 0 and the other at x = o.0 The cells are stacked vertically from the ground to the inversion base and horizontally from -- om to + X o The number of cells in Y is finite since the end cells have width AY = o The source cell is in the center of the array and is designated 0,0. All the properties of each cell are lumped in the center of the cross section so that the model is an array of rods with one end at x = 0 and the other at x = m o The source is no longer a point source but a vertical area source. In this model, the continuity Eq. (3) holds for all values of x, and by evaluating it at x = 0, we can define k o Q = H U S Xo dZ dY The integration at x = 0 is carried over the source cell only, and in a given cell all the parameters in the integrand are constant. Since the value of S at x = 0 in the source set at unity, Q = H2 u o (AY)o (AZ)o and therefore Xo -....'u H2 (AY)0 (AZ)o The solutions that will be obtained with this model will be the average concentration over the cross section of a cell as a function of its length, or distance downwind. 7 3. ANALOG MODELS In analog computation, it is helpful to realize that there exists a correspondence between the components of the circuit set up to solve Eq. (6) and the physical model shown in Fig. lo This can be shown by constructing on paper a passive network analog as shown in Fig. 2. In this circuit, there is a capacitor for each cell and a network of resistors connecting them. In this analog, the voltage at each node is a function of time and corresponds to the concentration as a function of distance downstream. The current into each capacitor is given by dV dX i = c- = u Ay Azdt dx and c = T Ay Az R l1 Ay R 1 Az R Ky A z z K- ^ A In this circuit the boundary conditions are implemented directly. Since there must be no flux through the ground or the inversion, the counterpart of flux, current, is made zero by simply making no connections across the boundaries. The initial concentration (voltage) in the source cell is provided by disconnecting the capacitor for that cell, charging it to arbitrary voltage, and switching it back into the network at x (time) equals zero. Since the end cells in y are infinite, the concentration (voltage) is always zero there so these nodes are connected to ground Provision for variation of the parameters K, u, Ay, and Az as a function of y and z is made by assigning appropriate values to the re8 INVERSION. I I I I I I I I I I I I I I I I I I I I I S-, I So, Si I 000r%* I I - ' II I I I I I ILZ I I. GROUND Fig. 2. The passive network analog of diffusion from a point source. The center cell is the source cell and its capacitor is initially changed to represent the source concentration. 9 sistors and capacitors involved. If any of these parameters are a function of x, then they must be made to change in time. This is not readily accomplished and constitutes one of the disadvantages of the passive network. In general, the passive network is much less flexible than the active circuit. An electronic analog circuit will be shown later for each case. They are constructed by assigning to each component a mathematical operation and these are connected so as to satisfy the governing set of differential equations. Identical circuits can be derived by simulating the passive network analog; thus each circuit has a dual nature. It is a series of mathematical operations designed to solve an equation and it is also an analogy of the system being studied where components and voltages represent physical parts of the atmospheric diffusion model. Both the active and the passive analogs are physical systems designed to behave like another, idealized physical system. The electronic analog circuits will be explored as they appear. 4. COMPUTATION AND RESULTS Since each of the cases considered is quite different, the form of Eq. (6) appropriate to each will be given along with the computing circuit and the results. A. ISOTROPIC DIFFUSION WITH AN 84-CELL MODEL The first case is diffusion from the point source in isotropic conditions, that is, the wind speed and eddy diffusivity are constant over the region. The 84-cell model used is shown in Fig. 1. It is symmetrical about 10 the coordinate axis so only the cells in the first quadrant need be simulated and of these, the ones in the outer column, m = 3, are always zero so that there are only 14 active cells in the simulation. In this case, Eq. (6) becomes dS 1 Sm+l.n Smn Smn - dX m,n 196 (AY)m Am+i (Y) -l +Smn+l - 2Sm,n + Sm,n-l 196 (AZ)2 where X = 196 x and the boundary conditions are ln = S-l,n Sml = Sm,-1 S1 (0) = 1.00 S.n = 0 Sm,7 = Sm,8 The factor 196 was inserted to make the coefficients more convenient, The cell boundaries given in Table I show that the spacing in the vertical is linear whereas the spacing in the horizontal is exponential and the source cell has a square cross section. 11 TABLE I CELL BOUNDARIES IN THE HORIZONTAL AND THE VERTICAL FOR THE FIRST QUADRANT OF THE 84-CELL MODEL Y = y/H Z = z/H 0 0 0-.0714 0.0714 0.1934 0.1428 00 0.2142 0.2856 0.3570 0.4284 0.5000 Since AZ.is a constant for each cell, Eq. (7) can be conveniently expressed as two sets of difference-differential equations for the columns used. dX = 08228(S2n - S,n) + (Sl,n+l - 2S,n + Sln-) l n dS dX 0.4937 Sn - 0.7410 S2 dX 2n l~n 2,n + (S2 n+l - 2S2,n + S2,nl) (8) The computer circuit for Eq. (8) which is shown in Fig. 3 has one amplifier for each of the 14 differential equations represented by (8). The top row of amplifiers simulates the first column of cells and the top equation in (8). The bottom equation is solved by the lower row of amplifiers. The S1 1 amplifier in the top row-represents the source cell and has an initial condition of S = 1.00 at computer time (distance in X) t = 0. The basic unit of computation is the integrating amplifier as shown in Fig. 4. for a typical cell. 12 +sl, 7 -J s +s - S 3, + -S 2, 2* 2 ' + 2, 45 2, Fig. 3. Computing circuit for the first quadrant of the 84-cell model. The integrating amplifier designated S1 1 represents the source cell and possess an initial condition of 1.00. 2, 7 -Sm-l,n +Sm,n -Sm,n-I -Sm, n+l I0 Fig. 4. Typical cell from the to the integrating amplifier. the coefficient potentiometers d-c AMPLIFIER -- COEFFICIENT POTENTIOMETER circuit of Fig. 3 showing the connections Nonunity equation coefficients are set on provided.. Solution curves for the first row of cells are shown in Figs. 5 and 6. B. ISOTROPIC DIFFUSION WITH A 56-CELL MODEL To facilitate the simulation, the number of cells can be reduced. In this case, only 4 cells in Y are to be used instead of six. The purpose of this reduction is that more involved problems to be treated are more readily handled with a simpler basic model. Not so much equipment is needed in this case but the reduction in the number of cells imposes a severe accuracy penalty. Before introducing the more interesting cases, the magnitude of the error introduced can be shown by solving the same problem employing the same techniques. Cell boundaries in Y for the 56-cell model are Y = 0 0.0714, o; while the cell boundaries in Z are as stated in Table I. In the first quadrant, there are still 7 active cells in Z but only 1 in Y; thus there is 1 active column and Eq. (8) simplifies to X = - 0.7220 S1,n + (S1n+ - 2S1 n + Si n-l) (9) X,n n where X = 196 x and the boundary conditions are S n = Siln S.1l = Sl_-l S1 1(0) = 1.00 S2n = 0 S1,7 = Si8 The computer circuit for Eq. (9) is shown in Fig. 7. As in Figo 35 there is 15 1.00 0.90 - = 196K x 0.80 - U/H2 0.70 - 0.60 z 0 o_ o 0.50 0.40 - 0.30 \ 0.20 -1 2 0.10 i1 - 0 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND. X Fig. 5. Concentration plots for each cell as a function of distance downwind from a point source using the 84-cell model in isotropic conditions. 16 0. 0.18 0.16 0.14.% I, 1,2 S= X/Xo 1 /96K X= uH X uH2,) O.KI C) I —: 0.10 z LIiJ 0 o 0.08 0.06 0,02 /,7 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND X Fig. 6. Concentration plots for each cell as a function of distance downwind from a point source using the 84-cell model in isotropic conditions. 17 (a) -S, +S, '7,..... ', 3 1.. 1.. OD (b) -S —, nl - -S — n-I Fig. 7. Computing circuit for the first quadrant of the 56-cell model shown in (a). An integrating amplifier for a typical cell in shown in (b). r one integrating amplifier for each active cell. The solutions for the source cell obtained from the 56-cell and the 84-cell models are compared in Fig. 8. Evidently the approximations are not bad for small values of X o For all the following work, the absolute error for values of X > 1 is large, the objective is to show deviations from the reference solutions in each case. C. NONISOTROPIC DIFFUSION Changing the ratio of the eddy diffusivity in the Y-direction to that in the Z-direction produces a very marked alteration in the solutions. As an example, set Ky = 10 Kz in the 56 —cell model. Equation (9) may now be written dS S l,n+l - 2S1,n + Sln-, dX = - 007220 S1 n+ 10 (10) l,n where X = 1960 x and the boundary conditions are S S S1 = S 1 l,n -ln S S 1,1 l,-l S1 (0) 1= o00 S2 0 2,n S =5 S1,7 = S1,8 The computer circuit for Eqo (10) is identical to the one shown in Fig. 7 except that the coefficient potentiometers are set to the values indicated in Eq. (lo). Solutions for the cases Ky = Kz and Ky = 10 Kz are plotted in Fig. 9. It shows that increasing the lateral diffusivity in this model 19 1.00 0.80O 1-.. 0.60 s = X/xo 0.40 196 K ' H -J\ -J o 0.20 --- —z 584-CECLL MODEL 0.10I23456 0.06 --- z DISTANCE DOWNWIND X Fig. 8. Comparison of concentration plots for the source cell using the 84- and the 56-cell models with isotropic diffusion conditions. 20 I. 1.10 Un z o 0.01 z z 0 0 -M 1 96 K X' Hm -X 0 0.1 0.2 0.3 0.4 0.5 0.6 07 DISTANCE DOWNWIND X Fig. 9. Comparison of concentration plots for some cells using the 56-cell model. This plot shows the changes produced in the solution by increasing the lateral diffusivity by 10 times over the vertical diffusivity. 21 has a pronounced effect. The magnitude of the effect is enhanced by the limited vertical diffusion due to the presence of the inversion. D. EDDY DIFFUSIVITY RELATED TO PLUME DIMENSIONS Plume dimensions may influence the magnitude of the effective eddy diffusivity acting on the plume since, as the cloud grows, larger eddies are able to operate to diffuse the materials. The eddy diffusivity might be made a function of the plume size and of the power spectrum of eddies present. In the usual parabolic diffusion equation, the plume has no easily defined boundaries except those imposed upon it, such as the surface of the ground. Therefore the diffusivity cannot be related directly to plume size but must be specified in some arbitrary manner. For example, it could be made proportional to the standard deviation of the concentration distribution taken perpendicular to the plume axis. If this distribution is normal, the standard deviation is inversely proportional to the central maxima or axial concentration, so that the diffusivity can be set inversely proportional to the axial plume concentration. Since the power spectra would never be continuous indefinitely, the diffusivity could not become infinite as the axial concentration went to zero but would reach some finite value and remain constant. In the finite difference model, the source is finite so that the diffusivity would increase from one small finite value to a larger one inversely as the concentration in the cell containing the source. Represent the diffusivity as 0(S1)K where K is the diffusivity for infinite plume dimensions and 0(S1) is a function of the axial concentration such that 0< (Si) < 1. 00. 22 For example, set 0(sl) = - 0.7 S,1 (11) since 1.00 1 S11 0 Setting Ky = Kz, Eqg (9) becomes dX 1 = (Si)(sn+l - 2Sl,n + Sl,n-l) - 07220 0(S1) Sln (12) with the same boundary conditions as before. The computing circuit, shown partially in Fig. 10, is similar to that of Fig. 7 in principle. This one uses 3 amplifiers, one integrator and two summers, per cell to isolate all the inputs to the integrating amplifier. With this configuration, it is only necessary to introduce the modifying parameter 0(S1) once per cell. Since 0(S1) is a function of the axial concentration, its form is determined by a feedback look defined by Eqs. (11) and (13). The function 0(S1) is generated as the solution proceeds, The computing circuits used in this problem are virtually identical with those of the next section. Therefore the description of the additional components required for the generation of 0(S1) and the presentation of the solutions will be deferred to that sectiono E, DIFFUSION IN TRANSITIONAL STATES A problem of special interest arises when the trajectory of the diffusing plume crosses the boundary from one turbulent regime to another. The physical model consists of a source located over land not far from the shore of a lake. The wind blows from the land out over the lake and the mechanical and perhaps 23 the thermal turbulence decrease markedly as it does so. The associated diffusion is referred to in this research as diffusion in a transitional state; the inversion breakup fumigation is the result of a corresponding time transition in diffusion. In the mathematical model, a change in the eddy diffusivity accompanies the transition from one turbulent regime to the other. This will be implemented by applying another modifying parameter 0(X) to the constant K o Then 0(X) will be unity in the first regime and fall to some lesser value in the second. For the present preliminary study, this function was simulated directly on the computer to conform to an intuitive impression of what it should be without using a mathematical expression. It is possible to contrive a mathematical description, which is: $(X) = 0.52 + o.48{l - h(X - X) + cos(X - X1) Fh(X - Xi) - h(X - X1 - r)]} (13) where h(a) is the unit step function which has the following properties: h(a) = 0, a < 0 = 1, a >0, and X1 is the value of X where the turbulence starts to decrease. The diffusion equation is just Eq. (12) with 0(S1) replaced with 0(X) The computing circuit is also that of Fig. 10 and the circuit for generation of 0(S1) and 0(X) is shown in Fig. 11. The plysical problems discussed in this section and in the previous one 24 +S, 2 +1 \ r. 2' +S, L -S- 2 + 1, 3 1,2 2.722 ~ SUMMING AMPLIFIER +S, 3, +Sg, 1. 2.722 ^ --- —-} +S, +SI. 1. L — -- - - - 1, SERVO 2.722 I --------------------— I (S I,.X ) 1.722 MECHANICAL LINKAGE TO FOLLOW-UP POTENTIOMETERS Fig. 10. Computing circuit for the first quadrant of the 56-cell model which permits the inclusion of special forms of the coefficient of eddy diffusivity. This configuration permitted the generation of solutions where the diffusion coefficients were a function of plume size and in the case of diffusion in transitional states. The dotted lines indicate mechanical connections. 25 DIODE LIMITERS ODIODE 0.500 AMPLIFI ER Fig. 11. Computing circuit for the generation of the special forms of the coefficient of eddy diffusivity. The diode limiters around amplifier 1 serve only to keep that amplifier from saturating and they do not contribute to the generation of the function. a are different despite a somewhat deceptive similarity in the governing equation used. The differences between the problems are indicated in the form of the modifying parameter 0 o It is only because of the similar way in which 0 is implemented that the problems are handled togethero It is possible to introduce both effects in the problemso In this case, a new parameter will be defined such that 0(S1,X) = 0(Sl) 0(X) o Figure 11 shows how each of three functions is generated. When switch 1 is on -100 and switch 2 is closed, the output of amplifier 4 is 0(S1) as defined by Eq. (11) and it passes through servomultiplier 1 unchanged to control servomultiplier 2. When switch 1 is on 0(X) and switch 2 is open, the function 0(X), generated by amplifiers 1, 2, and 35 controls servomultiplier 2. The function 0(S1,X) is generated when switch 1 is on 0(X) and switch 2 is closed. The functions 0(Sl), 0(X) and 0(S1,X) are plotted in Fig. 12 as they were generated by the computer. Solutions for Sll, S1,2 and Sl,3 are presented in Figs. 135 14, and 15, respectively, for 4 cases: 0 = 1, 0 = 0(Sl), 0 = 0(X) and 0 = 0(Sl,X). The case 0 = l is the reference solutiono Each modifying parameter acts to change the scale in X, that is, the new solutions can be found by translation in X from the reference solution. This is shown best in Figs. 14 and 15 where the maxima have been translated and expanded but the value of the maxima remains the same. The first modifying parameter 0(S1) is significantly different from unity only for small values of X so the greatest effect on the solutions should occur thereo Each solution shows large translation for small values of X o The net effect of coupling the eddy diffusivity to cloud dimensions is to increase the concentration values for small X 27 I (s/,) (X LL /.20 -a..20.10 -0 1 2 3 4 5 6 7 DISTANCE DOWNWIND X Fig. 12. Plot of the modifying parameters of the coefficient of eddy diffusivity. 0(S1) is a function of plume size while 0(X) represents transition from one turbulent regime to another. 0(S1,X) is a simple multiplicative combination of 0(X ) and 0(X). 28 1.00 z o.50 Fz Z \ \ 0.40.30.20 10 -0 I 2 3 4 5 6 7 DISTANCE DOWNWIND X Fig. 13. Plot of the concentrations in the source cell as a function of the-modifying parameter used. 29 .I.I.1.1 N _rr U) z 0 z w z.0.1 (S,, X).06.04.02 0 Fig. 14. modifying 5 6 7 DISTANCE DOWNWIND X Plot of the concentration in cell 1, 2 parameter used. as a function of the 30 .0 -z 0 h0' w z 0 o) (S,,X).0;.01 I l I_ 3 4 DISTANCE DOWNWIND X 5 6 7 Fig. 15. Plot of the concentration in cell 1, 3 as a function of the modifying parameter used. 31 along the plume axis and to delay the occurrence of the maxima off the axis. The solutions for 0(l) and 0(S1) do not converge for large X but reach a constant translation in X. But since the gradient is small for large X, the absolute value of the differences between the solutions is small. Since the function 0(X) can be determined independently, the solution could be found from the reference solution by plotting it on a new scale defined by rX -=X* 1r 0(X) dX The observed effect corresponds to this point of view in that 0(X) has no effect on the solution for X.- 1 but acts to translate the solution for X >1. F. WIND DIRECTION SHEAR The presence of wind direction shear in the diffusing layers causes the plume cross section to be skewed, and the amount of skewness can be determined by inserting terms in the diffusion equation for cross-wind components. Equation (6) contains the nondimensional parameter V = H v/K for these cross-wind components. Since the diffusion pattern is not symmetrical, it will not suffice to simulate only the first quadrant. In the interest of simplicity, a new 25-cell model is used as shown in Fig. 16. The source cell is 0,3 and there are 15 active cells. Figure 17 shows V(Z), which is an exponential function of height, in the finite difference model in relation to the downwind component H u/K. The exponential function was chosen as an approximation to the lower part of the Ekman spiral. Equation (6) may be tz H 1.00 -2,5 -1,5 0,5 +1,5 +2,5 ___________________.80___________________________ -2,4 -1,4 0,4 +1,4 +2,4 -2,3 -1,3 X0 0,3 +1,3 +2,3.40 -2,2 -1,2 0,2 +1,2 +2,2.................. -2,1 -1,1 0,1 +1,1 +2,1 - [ — --.1 -... ii | I.1 —1 -.60 -.10 0.10.60 y H Fig. 16. The 25-cell, finite difference model of diffusion from a point source. This model is used in the wind shear case. z 5 4 V(Z) 22 1.00 X H - Fig. 17. A schematic of the wind shear component V as a function of height in the 25-cell finite difference model. written dS 1 Sm+ln - mn Smn S Si nn Ixm n 25(AmY)mm (AY)mm+ (Y)-1 m,2 2 (Y)m ESm+l.n - Sm-ln + [Smn+i 2Sm n + Sm n-1] (14) where X = 25 x and the boundary conditions are -2,n = 0 S = 1 S2,n= 0 Smo = Sm 1 Sm,5 = Sm,6 So 3(0) = 1.00 The computing curcuit is shown in Fig. 18 and the solutions are in Figs. 20, 21, and 22. Each solution is compared to the reference solution, where V(Z) = 0. We can compare the concentrations in the cells for two cases at some distance downstream, say at X = 1.* The results, as shown in Fig. 19, indicate that the wind direction shear produces a marked skewness in the plume. The skewed plume is to be expected but a more interesting question is whether the average concentrations in the plume are reduced because of the skewness. Figures 20, 21, and 22 indicate that the concentrations. are reduced, *In making this comparison, note that changes in wind direction were provided by adding v components while holding u constant which increased the magnite of the total vector by an average 1.6.o Since X is a function of, we can compensate for this by evaluating the concentrations when V 0 at X = 1 and when V = 0 at X = 1016.. 35 S- 5 -.- +So 1.00 -Sl +S,2 -s +Sl -S,, +SOi -So,2 +S0o3 -So 4 +-So 5,-St I SI 2 -SI$ 1,4 -S1,$ S 6 -Fig. 18. Computing circuits for the 25-cell model 'or the case with wind shear. m _ -I 0 +1 n,.?..\\, 31 04 13 _ /79/ //A\\ at i i I \1 ill ~ i ii 2 (5 WH S (a) WITHOUT WIND SHEAR (b) WITH WIND DIRECTION SHEAR X =1.C00 Fig. 19. Cross section of the plume at X = 1 showing lines of constant concentration in arbitrary units. This figure shows the reference solution (a) and solution in the presence of wind direction shear with height (b). 37 0.020 0.018 0.016 0.014 00.012 \ a: 0.010 - 0 \ 00.008 -\ N 0.006 -,, 3 0.004 - ^.. 0.002 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND X Fig. 20. Comparison of the concentration plots for cells to the right and the left of the source cell in the presence of wind shear. The solid line is the reference solution for both cells. 38 0.020 0.018 0.016 0.014 tt 0.012 z - 0.010 - z 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND X Fig. 21. Concentration plots for some cells in the 25-cell model in the presence of wind shear. The solid lines are reference solutions. 0.006 N 0.004 0.002 0 2 4 6 8 10 12 14 DISTANCE — DOWNWIND X Fig. 21. Concentration plots for some cells in the 25-cell model in the presence of wind shear. The solid lines are reference solutions. 39 0.020 0.018 0.016 0.014 o 0.012 0.010- / // W "IX H / OI - I I N 1, 0.008 -4 Q002 -,7 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND X Fig. 22. Concentration plots for some cells in the 25-cell model in the presence of wind shear. The solid lines are reference solutions. 40 by about 4% in the center of the plume at X = 1. The integrated concentration over the plume at X = 1 was reduced about 2%. Thus it seems valid to conclude that wind direction shear does increase the diffusing properties of the atmosphere. It is planned to analyze the influence of wind direction shear in concentrations downwind from a continuous point source using an analog computer of much greater capacity. G. LARGE-SCALE EDDIES The presence of large-scale eddies or, since we are considering the steady state, a standing wave in the horizontal wind components can be represented as a sinusoid v = a sin bx. A stationary plume might be an example of this kind of process. The small-scale.eddies.which operate to diffuse the plume are represented by the coefficient of eddy diffusivity, whereas the large-scale eddies act only to transport the plume in the y-direction. Equipment limitations at the time of the research dictated the choice of a more coarse finite difference model as shown in Fig. 23. The mathematical model, based on Eq. (6), is dS mX= Sm+l n + Sm-ln + Sm n+l + mn-l - 4Smn m,n 18(AY) ESm+ln - Sm-ln] (15) 18(AY) Ha where V(X) = - sin BX and the boundary conditions are S2 n = S-2,n = so,3(0) = 1.00 41 to z! INVERSION It 1.00I 0 0 0 (9 -4 — ( I.6 r4) -- ~ ~ - - - - - - + 0.667i -- -- I -— 0 — - Y ro CJ a03331 I 0 0 0 I -.. -... - - - - -. - -- - - -| - - --- - - - -|I EARTH -0.667 -2 -I -0.333 0 +0.333 +1 +0.667 +2 Ii C: Fig. 23. The 15-cell finite differnce model of diffusion from a point source for use in the case of diffusion in the presence of large scale eddies. Sml = S, 2 Sm,5 = Sm,4 and AY = AZ for the active cells. If we set 5A = Ha(AZ)/2K, (15) becomes = [m,n-l - 4Sn + (1 - 5A sin BX) Sm+ln + (1 + 5A sin BX) Sm-ln (16) where X = 90 x. The computing circuit is shown in Figs. 24 and 25, and a typical cell is shown in Fig. 26. The solutions are plotted in Figs. 27, 28, and 29 along with the reference solution. This system exhibits a phenomena like resonance in a forced second-order system. It is not pure resonance but the sytem does respond markedly to changes in frequency of the cross-wind component. Changes in amplitude do not affect the qualitative nature of the solution. The plume meander lags behind the cross-wind components. The point at which the plume axis crosses the middle cell can be found by noting the point where the concentration of opposite pairs of cells is equal. This occurs for X = 1.506, 2.448, and 3.531, while the cross-wind goes to zero at X = 1.047, 2.094, and 3.141. The phase lage is not constant for the above three points but is 78.9 deg, 60.9 deg, and 67.1 deg, respectively. It is of interest to compare the concentrations in the plume, for sampling periods short enough to avoid plume shifting, when it meanders with when it does not. At an arbitrary point downwind, at X = 1.6, the total concentraion reported in the case with no meander was 16%o higher than in the case 45 +SI 4 20 + A sin BX - - - - - - - - -. - - - -J - - ---- -- -- - - -- I -j 20-A sin BX 0, 4 20 + A sin BX ----— 20- - - -A sin K 20-A sin BX S-1, 2 Fig. 24. Computing circuit for the 15-cell model including provision for variation of the horizontal wind component. +A sin BX \J1 SERVOMULTIPLIER / -100' MECHANICAL DRIVE I TO POTENTIOMETER I 20-A sin BX 20 + A sin BX Fig. 25. Computing circuit for the generation of the sinusoidal horizontal wind component. 20 +A sin BX S I m+i,n I I I SUMMING AMPLIFIER Sn mI Sm, n-i 0' sm-i, INTEGRATING AMPLIFIER Sm, n-+ i I I -I -! - --— I 20+A sin BX Fig. 26. Typical cell from the circuit of Fig. 23. The numbers inside the amplifier symbols represent the gain given each input. 1.00 ).10 0 n tr I — 0 0.01 0.001 V(X) = SIN 3X 0,2 \0,4."I% N. %... 014% q%% -4%, 0 0.5 1.0 1.5 DISTANCE Plot of the concentrations of large-scale eddies. The 2.0 2.5 DOWNWIND X 3.0 3.5 in the solution. Fig. 27. presence for cells 0,2, -0,3, and 0,4 solid line is the reference 47 V(X) = S/N 3X -/,3 /' "\ 0.10 \\ C \ W \/ z \ /. \ \ — Z - \ / \ 0.01 -- "._. 0.001 0 I I 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 DISTANCE DOWNWIND X Fig. 28. Plot of the concentration for cells -1,53 and +1,3 in the presence of large-scale eddies. The solid line is the reference solution. 48 1.00 V(X) = SIN 3X -/ 4 oi - -— ^ I — 0.10. z 0 z 0.01 0.001 I... I.I 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 DISTANCE DOWNWIND X Fig. 29. Plot of the concentrations for cells -1,2; +1,2; -1.4; and +1,4 in the presence of large-scale eddies. The solid line is the reference solution. 49 where the plume meanders. As shown in the appendix, this decrease in concentration in the meandering plume is due entirely to the longer path length followed by the plume. If one knew the total concentration across a straight plume as a function of distance downwind one could then compute the long-term concentration at a station subjected to a plume in meandering but otherwise identical conditions. The calculations may be so tedious as to be impractical since they would involve computing the length of the trajectories over the given interval, adjusting the concentration and integrating them. It would still be of interest to perform this calculation on the computer using random large-scale eddies, and it is intended that this topic will be the subject of a future report. H. A POSSIBLE COMBINATION OF TWO OF THE MODELS Two of the models presented above are related to some extent. The model in which eddy diffusivity increases with the plume dimensions can be related to that including large-scale eddies. In the first case, the complete spectrum of eddies can be included and the effect on the plume partially determined. The solution is partially complete in that, while the diffusing effects of all the eddies are included, the transport due to eddies much larger than the plume dimensions is neglected. The model simply ignores eddies of a given scale until the plume dimensions grow to that scale. The transport phenomenon is shown in the latter model where the large-scale eddies are treated as a wind component, but here they contribute nothing to the diffusion process no matter how large the cloud becomes. A combination of the two models might be an improvement. Let the wind components represent eddies larger than the plume while all smaller eddies act through the eddy diffusivity. The wind components would become smaller with distance downstream whereas the eddy diffusivity would increase. The wind components incorporated in the model are defined as averages which raises the question of how much they should be allowed to fluctuate or to possess time - or distance-variant properties. An averaged wind speed could be defined as the average or mean value obtained by the use of a sampling time just large enough to exclude all the turbulent eddies which contribute to the diffusion process. If we allow the assumption that the scale of eddies involved depends upon the scale of the plume, then the sampling time in the above definition would depend upon the plume size also. The wind-speed components could be represented as a Fourier series whose higher-frequency components are gradually eliminated as the plume grows, and conversely, the coefficient of eddy diffusivity would increase as more turbulent eddies were added, to its domain. This is a rather involved model even for analog computation and no attempt as yet has been made to set it up. 5o SUMMARY AND CONCLUSIONS As the model has become more complicated, it has been necessary to use fewer cells and to tolerate greater error. The amount of equipment needed is proportional to the number of cells and to the complexity of the model, and all the research described here was done on one machine with a very limited amount of computing equipment The direct way to reduce 'error, then, is to use a larger machine which would permit simulation of more cells. The errors are proportional to the fourth partial space derivative so that the error can be reduced by solving for a function related to the solution which possesses a smaller fourth derivative. Such a function could be found by perturbation techniques. The parabolic diffusion equation has been shown to be of general utility in the variety of problems covered. It is a model which is well adapted to analog computation. With sufficient equipment the actual programming seldom presents any difficulty even in the most involved problems discussed here. The trapping situation was used throughout this report to limit the scope of the research,and to facilitate handling the problems on a computer of given size. Extension of these methods to nontrapping situations can be readily done using the technique reported earlier. 52 APPENDIX In the case where the meandering plume was treated, there were two possible mechanisms which could affect the diffusion of the plume. The one noted in the main body of the report was the increased path length of the trajectory. The other is that large-scale eddies might change the shape of the plume, thereby influencing the diffusion process. An attempt was made to evaluate this mechanism by correcting for the first one. When there is a cross-wind which causes meandering, the magnitude of the total wind vector and the path length along the plume axis, as shown in Fig. 50, has been increased. To compensate for this the concentration for no cross-wind should be taken at a point further downstream in inverse proportion to the relative decrease in the magnitude of the wind vector. Taking the first zero crossing, at X1 = 1.506, the compensating factor is X2 = X1 VI7+ sin2 BX where 2 1 rxl 1 sin2 BX = sin BX dX = [BX - sin BX1 cox BXj = 0.4790 s 2BX0 - Then X2 = 1.506./.479 = 1.851. The following table gives the concentration in both cases. 55 Y / / / / / / / PL UME OUTLIN. ~ / / / / / / X AX/S \N —. — ^ -PLUME OUTLINE Fig. 30. Schematic of the fluctuating plume model showing the sinusoidal nature of the plume axis. TABLE A-I COMPARISON OF PLUME CONCENTRATIONS FOR THE MEANDERING PLUME AND THE STRAIGHT PLUME Cell 0,2 0,3 0,4 -1,2 -1,4 +1,2 +1,53 +1,4 Total concentration Straight Plume x = 1.831 0 0595 0.0596 0.0595 0.0412 0.0422 0.0412 0.0412 0.0422 0.0412 0.4278, Meandering Plume X = 1.506 0O.0599 0o0617 0.0599 0.o406 0.0420 0.o406 0.-o406 0.0420 O.0406 0.4279 - '' Evidently the plume shape is a little different since the individual cells do not report identical concentrations, but the sume of the concentrations for all the cells is not significantly different. 55 REFERENCES lo Sutton, 0. G., Micrometeorology, New York, McGraw-Hill Book Co., Inc., 19535. 2. Godson, W. L., "The diffusion of particulatematter from an elevated source," Arch. Meteorol., Geophys. und Bioklimatol., Serie A, 10, No. 4 (1958). 5. Goldstein, S., "On diffusion by discontinuous movements and on the telegraph equation," Quart. J. Mech. and Appl. Mech., 4:129-56 (1951). 4. Richardson, L. F., "Atmospheric diffusion shown on a distance-neighbour graph," Proc. R Soc., London, Series A, 110:709-737 (1926). 5. Gifford, Frank, Jr., "Statistical properties of a fluctuating plume dispersion model," in: Advances in Geophysics, Vol. 6, New York, Academic Press, pp. 117-138, 1958. 6. Hewson, E. Wendell, "Meteorological factors affecting causes and controls of air pollution," J. Air Pollution Control Assoc., 5:2-8 (1956). 7. Hewson, E. Wendell, "Atmospheric pollution," in: Compendium of Meteorology, Boston, Mass., The American Meteorological Society, pp. 1159-1157 (1951). 8. Hewson, E. Wendell, "Industrial air pollution meteorology," in: Vol. III, Industrial Hygiene and Toxicology, New York, Interscience Publishers, in press. 9. Brock, F. V., Analog Computing Techniques Applied to Atmospheric Diffusion: Infinite Line Source, Univ. of Micho ORA Report 05632-2-T, Ann Arbor, 1961. 10. Hewson, E. Wendell, "The meteorological control of atmospheric pollution by heavy industry," Quart. J. Rooy Meteor. Soc., 71:266-282 (1945). 56