,1 LIBRADY mchlgen State University ‘ L 1 PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. | DATE DUE DATE DUE DATE DUE MSU Is An Affirmative Action/Equal Opportunity Institution ‘ cadmium“: NUMERICAL SIMULATION OF THE AIR AND FUEL MIXING PROCESS IN FUEL-INJECTED ENGINES USING DISCRETE VORTEX DYNAMICS By Asuquo Bassey Ebiana A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1990 (045-61996 ABSTRACT NUMERICAL SIMULATION OF THE AIR AND FUEL MIXING PROCESS IN FUEL-INJECTED ENGINES USING DISCRETE VORTEX DYNAMICS By Asuquo Bassey Ebiana The numerical modeling of turbulent combustion problems continues to invite the application of novel analytical techniques. Contemporary numerical models seek to utilize the observable eddy structures of the turbulent flow to simulate the physical systems directly without recourse to conventional finite difference formulation. Vortex methods represent such a class of numerical techniques. Vortex methods are particularly appealing because they are suitable for use in high Reynolds number flows without the problems of spatial resolution, stability and numerical viscosity associated with finite difference schemes. The discrete vortex method is used in this study to simulate the mixing process prior to chemical reaction in the fuel-injected cylinder of a piston engine. Mixing descriptions are restricted to a two—dimensional, incompressible, inviscid, formulation of the bulk diffusion process. The analysis employs the fundamental conservation equations of mass, momentum and scalar species with appropriate simplifying assumptions. Although emphasis is limited to inviscid calculations, the the theoretical and numerical considerations relating viscous and molecular diffusion effects are fully addressed. The incorporation of viscous and molecular diffusion effects will be the focus of a future study. In the present study, the initial vorticity distribution inside the cylinder is approximated by discrete fluid vortices each of quasi Gaussian-shaped distribution. The initialization scheme makes use of prescribed velocity data and an imaging system for the vortices found by modifying the circle theorem of Milne-Thomson. Fuel injection is simulated using a discrete jet model analogous to Ashurst’s scheme. At each time step, a grid—insensitive Vortex-In-Cell (VIC) solution technique is used to approximate the velocity field inside the cylinder. Sub grid accuracy is achieved through quadratic spline interpolation using Couet’s numerical filter. Each vortex is advanced every time step using Heun’s second order accurate time integration scheme. Predictions of the numerical studies show excellent agreement with theoretical results obtained using the method of images. Plausible, qualitative flow patterns consistent with theory and intuition are illustrated. The primary conclusion of this study is that optimum condition is reached for a 60° injection angle, swirl rate proportional to 0.125 m2/s and 1.26 m/s injection speed. The artificial complete mixing time for this optimum condition is estimated to be 0.159 sec. The implication that there is an optimum swirl strength for minimum complete mixing time could be of significance in engine design. Dedicated to the memory of my father, Bassey Esang Ebiana, who instilled in me the value of education. ACKNOWLEGEMEN TS I would first like to express my deepest gratitude to Professor Merle C. Potter, my dissertation advisor, for the unselfish aid he has given me during the research and writing of this dissertation. Without his invaluable insight, new ideas, direction and guidance, much of the intricacies of this project would have been difficulty to unravel. ' I am very much appreciative of the generous disposition of his time. Thanks are also due Professors James V. Beck, Mahlon C. Smith and David H. Y. Yen for serving on my committee. Their comments, corrections and suggestions are well noted and appreciated. I am deeply indebted to the Mathematics Department for its generous offer of graduate assistantships without which I would not have been able to sustain myself and family financially this far. Special appreciation is extended to Barbara Miller, Kathy Trebilcott, Dr. Elizabeth Phillips, Professor Charles Seebeck and Professor Jacob Plotkin for their sensitivity to my financial plight. I am also deeply appreciative of the department’s recognition of my efforts by way of the Excellence-In-Teaching citation award accorded me. I wish to thank my department, the Mechanical Engineering Department, too for its financial support via laboratory teaching opportunities. In particular, the Fluid Dynamics Laboratory teaching experience was quite rewarding and I am grateful to Professor John F. Foss for refurbishing my ‘ undergraduate background with the fundamentals of experimental fluid dynamics and for his guidance and direction. Special thanks are extended to Dr. Kenneth Ebert of the Office for International Students and Scholars for his kind disposition and professional understanding of foreign students’ concerns. Appreciations are extended too to the Stafi of the Case Center for Computer Aided Engineering who acquainted me with the VAX/VMS, SUN and HP computer systems and helped me with system related problems. In particular, Mike Mcpherson, Margaret Wilke, Fred Hall, Greg Fell, Dan Mcknown and George Lupanoff deserve to be specially acknowledged. I would also like to thank Mr. Tom Volkening, the chief librarian for his help in searching for literature materials for me. Thanks too to my colleagues Dr. Bashar AbdulNour and Dr. Syed Ali for their friendship throughout the years. I wish to thank my mother, Mrs. Ukpong Bassey Ebiana for all she has done for me. My parents-in-law Elder and Mrs. Ene Okon Iban and family have been great sources of financial and emotional support and for this I am indeed very sincerely thankful. My late sister-in-law Carol Okpok Ene Iban (may God bless her soul) deserves special acknowledgement for her relentless efforts to encourage and spur me on. My cousins Mr. Nweme Eyo Eyo, Mr. Asuquo Eyo Eyo and Mr. Abang Okon Eyo and their families have also been emotionally very supportive and I thank them for their support. My sincere thanks and deep appreciation to a very special fiiend Dr. Abiodun Oyewale, for his encouragement and incredible financial sacrifice in my behalf. I am especially and deeply indebted to my wife, Bassey and daughter Victoria, for their encouragement in the face of great odds, emotional support, sacrifice, patience and understanding throughout the duration of my study. Above all I thank God Almighty for His care, protection and guidance always. TABLE OF CONTENTS Page LIST OF TABLES ....................................................................................... viii LIST OF FIGURES ....................................................................................... ix NOMENCLATURE .................................................................................... x111 INTRODUCTION .......................................................................................... 1 CHAPTERS 1. PREVIOUS INVESTIGATIONS ............................................................ 8 1.1 Flow Field .................................................................................... 8 1.2 Mixing Models ........................................................................... 11 2. DESCRIPTION OF PROBLEM ......................................................... 14 3. THE DISCRETE VORTEX METHOD ............................................... 19 3.1 Introduction ................................................................................ 19 3.2 The Flow Field ........................................................................... 20 3.2.1 The interior region .......................................................... 23 3.2.2 Boundary layer region ..................................................... 28 3.2.3 Intake flow ...................................................................... 34 3.2.4 Summary ........................................................ 35 3.3 Exact Solution ............................................................................ 36 iv 3.4 Mixing ........................................................................................ 39 3.5 Choice of Numerical Parameters and Functions ........................ 46 NUMERICAL FILTERS ..................................................................... 47 4.1 Introduction ................................................................................ 47 4.2 Description of Filters ................................................................. 47 4.3 Analysis of Filter Performance .................................................. 51 4.4 Summary ..................................................................................... 54 NUMERICAL METHOD .................................................................... 56 5.1 Conformal Transformation .......................................................... 56 5.1.1 Grid transformation ......................................................... 56 5.1.2 Velocity transformation .................................................. 59 5.2 Point Vortex Initialization Scheme .............................................. 61 5.3 Vortex-In-Cell (VIC) Algorithm ................................................. 63 5.4 Vortex Sheet Algorithm ............................................................. 66 5.5 Advection Scheme ..................................................................... 68 5.6 Sorting Algorithm ....................................................................... 71 5.7 Mixing Algorithm ....................................................................... 72 NUMERICAL EXPERIMENTS AND DISCUSSION OF RESULTS .................................................... 74 6.1 Introduction .................................................................................. 74 6.2 Optimization of Numerical Parameters ....................................... 75 6.2.1 Mesh size ......................................................................... 75 6.2.2 Injector nozzle diameter .................................................. 77 6.2.3 Vortex core radius ........................................................... 78 6.2.4 Time step .......................................................................... 81 6.3 Validation of Numerical Codes ................................................. 83 6.3.1 Conformal transformation ................................................. 83 6.3.1.1 Grid transformation ........................................... 83 6.3.1.2 Velocity transformation .................................... 84 6.3.1.3 Corner singularities ............................................ 85 6.3.2 Point vortex initialization .................................................. 86 6.3.3 Velocity field and vortex trajectories ............................... 87 6.3.3.1 Velocity field .................................................... 87 6.3.3.2 Vortex trajectories ............................................. 88 6.4 Numerical Simulations .............................................................. 91 6.4.1 Swirling motion with no jet flow ...................................... 95 6.4.2 Jet flow with no swirling motion ...................................... 96 6.4.3 Swirling motion with jet flow ........................................... 97 SUMMARY AND CONCLUSION .................................................... 99 RECOMMENDATIONS .................................................................... 103 TABLES ...................................................................................................... 106 FIGURES ...................................................................................................... 1 13 APPENDICES Numerical Filters ............................................................................... 168 Al Derivation of Christiansen’s filter as an illustrative example 168 Functions ........................................................................................... 173 8.1 Smoothing function Gij for convergence method [39] ............ 173 vi B.II Derivation of the streamfunction ‘1’ inside a finite rectangular domain using method of images .................. 174 C. Coordinate Transformation ............................................................... 179 CI Square (Cm) -> Circle (x,y) ...................................................... 179 CH Circle (x,y) -> Square (Cm) ....................................................... 181 D. Evaluating the Strengths of the Point Vortices .................................. 187 E. List of Computer Programs ................................................................. 189 BIBLIOGRAPHY ........................................................................................ 286 Table 1. Table 2. Table 3. Table 4. Table 5. Table 6. Table 7. Table 8. LIST OF TABLES One-dimensional filter models. Dyer’s experimental data. Comparison of the absolute value of the relative error (%) between the modified exact and numerical tangential velocity values at grid points. Error statistics for grid point coordinate and velocity transformations. Velocity field initialization data. L, 0- optimization data (swirl strength S0 = 12.5 cm-m/s). Nc optimization data (swirl strength S9 = 12.5 cm°m/s). Optimization data for swirl strength 39' Figure 1. Figure 2. Rm 3. Figure 4. Figure 5. Figure 6. Figlue 7. Figure 8. Figure 9. Figure 10. Figure 11. LIST OF FIGURES Three phases of heat release. Sketch of intake flow. Numerical filters for vorticity deposition on the grid and for interpolation of local effects of the velocity field from the grid to the Lagrangian points. Vorticity distribution in a grid, with shading showing assignment of vorticity to grid point for VIC method. ISO-vorticity contours as represented by the various filters. Large tick marks represent a half-grid length. Fourier transform of the various filters. Inverse Fourier transform of the various filters. Error estimates between original filter and the inverse of its transform. Grid transformation. Quadratic spline vorticity distribution of three typical vortices on their nearby grid points. Zone of influence of a vortex sheet. Figme 12. Figure 13. Figure 14. Figure 15. Figure 16. Figure 17. Figure 18. Figure 19. Figure 20. Figure 21. Figure 22. Figure 23. Vortex sorting. Relative error between exact and numerical streamfunction values in the square domain using a single vortex and a system of 105 vortices. Relative error between the initial and final positions of a point vortex after one complete revolution. Radial distribution of tangential velocity component. Geometric distribution of the system of 115 discrete vortices (10:1m 1113). Relative errors in Vr and V6 velocity components as a function of time. Trajectories of systems of discrete vortices. Tangential velocity profiles at discrete time steps representative of short and long time errors. Radial velocity profiles at discrete time steps representative of short and long time errors. Uniform grid squares superimposed on circular domain. Complete mixing time as a function of the number of cells (Swirl strength 89:12.5 cm-m/s, Number of vortices per cell Nc=2). Complete mixing time as a function of the inlet jet velocity (Swirlstrength 89:12.5 cm-m/s, Number of vortices per cell N c=2). Figure 24. Figure 25. Figure 26. Figure 27. Figure 28. Figure 29. Figure 30. Figure 31. Figure 32. Figure 33. Average complete mixing time as a function of the number of cells (Swirl strength Se=12.5 cm°m/s. Number of vortices per cell Nc=2). Complete mixing time as a function of thenumber of vortices per cell (Swirl strength 80:12.5 cm-m/s. Number of cells L=32). Complete mixing time as a function of inlet jet velocity (Swirl strength 89:12.5 cm-m/s. Number of cells L=32). Complete mixing time as a function of the swirl strength (Number of cells L=32. Number of vortices per cell Nc=2). Complete mixing time as a function of the inlet jet velocity (Number of cells L=32. Number of vortices per cell Nc=2). Geometric distribution of the system of 115 discrete vortices and velocity vector representation of the induced swirling motion. Fuel particle trajectory and vorticity profile at y=0 as a function time. Injector at 60 degrees (counterclockwise from horizontal). Um=210.0 cm/s Streamlines and velocity profiles at y=0 as a function of time (without swirl). Injector at 60 degrees (counterclockwise from horizontal). Um=210.0 cm/s. Air and Fuel particle trajectory and vorticity profile at y=0 as a function of time. Injector at 60 degrees (counterclockwise from horizontal). Um=210 cm/s. Streamlines and velocity profiles at y=0 as a function of time (swirl). Injector at 60 degrees (counterclockwise from horizontal). Um=210.0 cm/s. ’ Figure 34. Figure 35. Figure 36. Figure 37. Figure 38. Figure 39. Figure 40. Figure 41. Figure 42. Figure 43. Fuel particle trajectory and vorticity profile at y=0 as a function time. Injector at 90 degrees (counterclockwise from horizontal). Um=210.0 cm/s. Streamlines and velocity profiles at y=0 as a function of time (without swirl). Injector at 90 degrees (counterclockwise from horizontal). Uin=210°0 cm/s. Air and Fuel particle trajectory and vorticity profile at y=0 as a function of time. Injector at 90 degrees (counterclockwise from horizontal). Um==210 cm/s. Streamlines and velocity profiles at y=0 as a function of time (swirl). Injector at 90 degrees (counterclockwise from horizontal). Um=210.0 cm/s. Fuel particle trajectory and vorticity profile at y=0 as a function time. Injector at 120 degrees (counterclockwise from horizontal). Uin=210.0 Cm/S. Streamlines and velocity profiles at y=0 as a function of time (without swirl). Injector at 120 degrees (counterclockwise from horizontal). Um=210.0 cm/s. Air and Fuel particle trajectory and vorticity profile at y=0 as a function of time. Injector at 120 degrees (counterclockwise from horizontal). Um=210 cm/s. Streamlines and velocity profiles at y=0 as a function of time (swirl). Injector at 120 degrees (counterclockwise from horizontal). Um=210.0 cm/s. Vorticity distribution in a grid, with shading showing assignment of vorticity to grid point for VIC method. Method of Images. NOMENCLATURE English Letter Symbols a - Radius of circle in Milne-Thomson’s Circle Theorem. C,CI,CZ - Constants. c,s,d - Jacobian elliptic functions with real arguments (also c1,sl,d1). cn,sn,dn - Jacobian elliptic functions with complex arguments. D - Finite circular domain representing the combustion chamber. 1). - Total derivative. Dt d - Distance between mesh points. dAiJ - Area element of vorticity assigned to the grid point (i,j). dj - Injector nozzle diameter. 911,512 - Elements of a 2 x 2 matrix. F - Matrix of geometric factors. f(z) - Complex potential function. f(¢,t) - Distribution function. Gij - Green’s function or vortex core smoothing function h - Boundary discretization length. Also, length of vortex sheet segment. xiii Grid index in the x (or C) direction. Grid index in the y (or 1]) direction. The complete elliptic integral of the first kind. Also, dimension of square domain. Modulus. Also, wave number. Number of square meshes or cells within the circular domain. Initial number of vortex sheet segments created at wall test points. Vortex sheet integer tags. Also number of dimensions. Number of grid points in either C or 1] directions. Number of vortex elements or fluid particles. Number of vortex elements required per cell. Number of indices. Unit normal vector on S. Averaged joint probability distribution function. Pressure of fluid. Also, instantaneous joint probability distribution function and parameter of interest. Radius of combustion chamber. Reynolds number. Discrete vortex core radius. Dimensionless coordinate r/R. Matrix of discrete point vortex strengths. Boundary of cireular domain. ‘ Swirl strength. xiv SI W(z) =(X.y) (x,y) z=x + iy Symbols Boundary of square domain. Unit tangent vector on S. Time. Artificial time characterizing complete mixing. Inlet jet velocity. Velocity at the outer edge of the numerical shear layer (i.e., velocity at"infinity"). Velocity component in the x-direction. Velocity component in the y-dircction. Velocity vector. Radial velocity component Tangential velocity component. Average tangential velocity. Matrix of tangential velocity components. Complex potential function. Position vector. Cartesian coordinates in the circular domain. Space coordinate in the complex circular plane. Dummy variable. Mixing fiequency. Also average distance between vortices. XV Length of time step. Half-widths for the vorticity distribution in Men g/Thomson’s filter. Numerical shear layer thickness. Dirac delta function. The del operator. The Laplacian operator. The coefficient of molecular diffusion. Turbulent diffusion coefficient. Discrete point vortex circulation. Vortex sheet circulation. Maximum allowable vortex circulation. Vector with Gaussianly distributed random variables. x—component random walk due to vortex j. y-component random walk due to vortex j. Of the order of. Discrete point vortex strength. Reduced wave number. Mean. Kinematic viscosity of fluid. Angle of injection. xvi .5 a, (C31) Subscripts Space coordinate in the complex square plane. Set of o-dimensional scalars. Velocity potential. Also, scalar function. Filter function. Density of fluid. Number of independent scalar species. Also, standard deviation (variance = oz). Set of o-dimensional probability spaces corresponding to (D. Dummy variable. Stream function. Also probability space. z-component vorticity. Non-dimensional cartesian coordinates in the square domain. Initial parameter. For constants and functions. Variable parameter, a = 1,2,...,o. Valve lip region. Valve seat region. Index in the x(or C) direction. Index in the y(or 11) direction. Also, jet. xvii p - Discrete point vortex. pt. - Piston. r - Radial component. 3 - Vortex sheet. Also, shear layer, and source. t - Tln'bulent. e — Tangential. C - Vorticity. Superscripts n - Time step number. 6 - Number of scalar species. t - Turbulent — - Mean. - Fluctuating component. Also, transformed boundary. Abbreviations CA. - Crank Angle. CID - Coalescence and Dispersion. CIC - Cloud-In-Cell. CPU - Central Processor Unit. ESCIMO - Engulfinent, Stretching, Coherent, Inter-diffusion, Moving, Observer. xviii FISHPAK IMSL NCAR PDF PLOTIT TDC 2D VAXIVMS VIC Fast Fourier Transform. A portable elliptic boundary value problem solver. Hewlett Packard. International Mathematical and Statistical Libraries. Laser-Doppler Velocimetry. National Center for Atmospheric Research (Graphics package). Probability Density Function. Interactive Graphics and Statistics. Top Dead Center. Two dimensions. Virtual Architecture Extended/Virtual Memory System. Vortex-In-Cell. INTRODUCTION In recent years, growing concern over automobile emissions and fuel economy has spurred great demand for further improvement of the internal combustion engine. Since efficiency and emissions are influenced by the details of the combustion process, it is logical that much of the recent engine research has focused on this process. The important details of combustion include processes such as mixing, chemical kinetics, vaporization, and so on. These processes are aided and abetted by turbulent motions which produce large concentration gradients and large surface areas of interface between fuel and air. To model these processes more accurately, it is necessary to incorporate new ideas and approaches as our understanding of fluid motion, turbulent structures, and their effects on combustion in an engine improves. The current knowledge of the above processes is limited, and so a precise prediction of these events in an engine has yet to be realized. Turbulent motion of fluids pose problems of great mathematical and numerical complexity. The theoretical analysis of turbulent motion has traditionally been by statistical methods. Because the mechanism of turbulent mixing is intricately linked to the dynamics of turbulent motion, the statistical theory of turbulent mixing has been developed parallel to the statistical theory of the turbulent motion problem. In either case, the choice of a statistical approach leads to a closure problem which 2 requires arbitrary postulates of the unknown terms. In turbulent motion, the closure is a direct consequence of the non-linearity of the Navier-Stokcs equations. However, in turbulent mixing it is caused by the non-linearity in the stochastic variables even though the scalar conservation equation is essentially linear. Because randomness and non-linearity combine to make the equations of turbulence nearly intractable, turbulence theory suffers from the absence of sufficiently powerful mathematical methods. In the absence of a completely formal mathematical theory, experimental measurements and computer modeling of the actual behavior in a firing engine have proven to be invaluable tools in the engine design process. Description of experimental methods to measure and observe turbulent fluid motion in an engine cylinder abound in the literature. These methods fall into either of two broad classifications: direct measurements or Schlieren photography. Direct measurements include the use of hot-wire anemometers for time-resolved temperature distributions, continuous-wave Raman spectroscopy for fuel and air concentrations, pulsed laser Raman spectroscopy for density fluctuations, Laser- Doppler Velocimetry (LDV) for the velocity and turbulence intensity and more conventional measuring devices for the cylinder pressure. The single point measurement and directional ambiguity of the hot-wire anemometer, the large pressure and temperature corrections that must be applied to hot-wire measurements and the fact that hot wires do not survive the combustion process present significant difficulties. The use of Laser-Doppler anemometry, which in principle has none of the above disadvantages, evolved from the recognition that engine design must be firmly based on the techniques of analysis and fundamental understanding of the in-cylinder processes. Schlieren photographic methods have been used to provide a highly 3 informative qualitative assessment of the fluid motion inside the combustion chamber. When coupled with other diagnostic efforts, the visualization allows the detailed behavior at a point to be placed in context with other phenomena within the chamber. Thus, the value of each diagnostic is enhanced by the application of the other. For example, the understanding of the flow structure would alleviate the directional ambiguity of the hot-wire anemometer and also enable investigators to choose more strategic locations for point measurements. A broad range of predictive models and mechanistic theories also abound in the literatrue. ESCIMO, "population balance", "eddy breakup", coalescence/dispersion and probability models are but a few of the prominent models proposed by researchers in the field of turbulent combustion. The most widely applicable and popular numerical implementation of the predictive models usually consists of finite difference formulations of the appropriate equations. Although such techniques work well in many cases, the need for a fine grid in the boundary layer where sharp gradients exist places a computational upper bound on the size of the Reynolds number that can be effectively modeled. Analysis implies that at least several grid points must fall within the boundary layer whose thickness is O(Re'm) and one always has a finite amount of space available on the computer. This problem is especially important in cases dealing with flows in wakes or separated flows when the initial boundary layer may not be visible to a coarse finite difference grid. Another drawback of finite difi‘erences is that in areas near the boundaries sharp gradients may give rise to large truncation errors which swamp the original approximation. This problem is further complicated by the truncation errors which may cause a numerical viscosity to form which is greater than the viscosity of interest, [32]. Also, strongly 4 interacting transient fuel sprays are difficult to model using finite difference techniques. The difficulties stem from the fact that the spray contains a wide range of droplet sizes, length scales and time scales. Furthermore, the inclusion of detailed combustion chemistry in this approach is far beyond the capability of modern computers, and cannot be seriously considered as a technique for modeling the combustion processes in practical combustion devices. The past two decades have brought enormous advances in computers and computational techniques on the one hand, and measurements and data processing on the other hand. The impact of such capabilities has led to a revolution both in the understanding of the vortical structures which are important in the mechanics of turbulence as well as in predictive methods for application in technology. The existence of identifiable, large scale coherent structures (eddies) in turbulent flows, as observed by Brown and Roshko [l2], Winant and Browand [72], Roshko [61] and Dimotakis and Brown [27], has led to the growing belief that turbulence or vorticity fluctuations are not so random as was commonly thought. This development has prompted several researchers to propose turbulent mixing and combustion models based on the creation, evolution, interaction and decay of these eddies. Vortex methods represent a class of this new generation of numerical models. The discrete vortex method was proposed by Chorin [16,17,18] for modeling fluid flow. In Chorin’s scheme, the vorticity in the fluid and in the boundary layers is subdivided into vorticity functions with finite cores (discrete vortices), and the equations of motion are solved by following the discrete vortices throughout the fluid. The fundamental advantage of this method is that it is suitable for use in high we“; link», ."w be». -\1“ .v— “as. it: S Reynolds number flows without the problems of spatial resolution, stability and numerical viscosity associated with conventional finite difference schemes. Some reliable results have been obtained using this scheme [2,16,64]. However, for Green’s function formulation [45], the dynamics of the vorticity field translates into a computational effort which scales as the square of the number of vortex elements. This is due to the fact that if N discrete vortices are present in the fluid, then O(N2) interactions must be computed per time step, since each vortex influences all the others. Furthermore, the discrete vortices move according to two components, one of which is a random displacement. This displacement leads to a partially random distribution of vorticity which becomes more random as the number of discrete vortices present in the fluid decreases. Hence to get accurate results, one must use many vortices or average over ensembles. In this study, air and fuel mixing processes are examined within the fuel- injected cylinder of a piston engine. The mixing process is important because it controls combustion which in turn controls the pollutant formation mechanisms. A detailed understanding of the mixing phenomena allows advanced engine concepts to be explored in search of higher efficiency, lower emissions and greater fuel flexibility. The idea of mixing is best described perhaps by the combination of bulk, eddy and molecular diffusion effects. Bulk diffusion includes diffusion other than eddy and molecular, e.g., diffusion caused by the bulk effect of fluid motion. In nu'bulent flows, the bulk motion of large groups of molecules (eddies) give rise to eddy diffusion. Molecular diffusion is a product of relative molecular motion across the boundaries of fluid elements. In any fluid system. molecular diffusion is always occurring and may 6 not always be neglected. Superimposed on this could be either bulk or eddy diffusion or both. The attainment of sub-microsc0pic homogeneity, where molecules are uniformly distributed over the field, is the ultimate measure of the effectiveness of mixing in any mixing process. Chemical reaction occurs only where the fluids are molecularly mixed. Fuel-injected internal combustion engines typically employ injectors producing finely atomized sprays with complex distribution patterns and time-varying injection rates. The interaction of the injected fuel jet with the primary air motion inside the engine cylinder compounds the complexity of the fluid motion. The ensuing chemical kinetics of the combustion process makes the complexity of the fluid motion even more staggering. Given this consideration, attention is limited to the mixing- dominated phase prior to chemical reaction. Specifically, the bulk effect of the inviscid fluid motion on the scalar field is numerically simulated. The two—dimensional, inviscid calculations employ the discrete vortex method developed by Chorin in conjunction with the Vortex-In-Cell (VIC) solution technique to approximate the fuel intake and the mean flow generated by the primary air swirl inside the combustion chamber. The VIC technique combines the best features of the vortex tracking concept of the discrete vortices with the stream function-vorticity finite-difference method. It employs a numerical filter to spread the vorticity onto a grid for use by a fast Fourier transform routine. The VIC technique is conceptually more efficient than Green’s function or Biot-Savart integration techniques, in part, because the work load increases linearly with the number of vortex elements. Although the numerical simulation of the viscous and eddy diffusion effects and n.” ‘ AU vi? A {We 7 molecular activity are neglected for simplicity, the theoretical fiamework relating these effects are fully addressed in this study. The basis for the simulation is discussed and tested against theoretical results for the velocity field obtained by the method of images. The numerical predictions of the effect of bulk diffusion on the scalar field are qualitatively assessed and shown to be consistent with theory and intuition. Finally, a conclusion is drawn with regard to the optimum swirl condition for minimum mixing time. The plan of the rest of this dissertation is as follows: Chapter 1 provides a brief historical perspective on the evolution of the discrete vortex method. The direction of earlier efforts, difficulties encountered and the work of previous investigators are outlined. A discussion of some of the mixing models consistent with the ideas and techniques of the discrete vortex method is also presented. Chapter 2 provides a clearer focus on the problem under investigation. The problem parameters, assumptions and - governing equations are outlined. A more in-depth discussion of the theoretical basis of the discrete vortex method and its numerical implementation are presented in Chapters 3 and 5. Chapter 4 provides a description of the available numerical filters in the literature for use in the VIC method and a discussion of the considerations leading to the choice of an optimal filter. In Chapter 6 numerical experiments and discussion of results are presented. A summary and concluding remarks regarding the study are provided in Chapter 7 and in Chapter 8 suggestions for future studies are outlined. L1 CHAPTER 1 PREVIOUS INVESTIGATIONS 1.1 Flow Field Calculations using vortices to simulate the gross features of an incompressible two—dimensional fluid motion have a long history of evolution from Rosenhead’s [60] hand calculations in 1931 using a few line vortices to recent computer calculations involving hundreds of vortices. Rosenhead was the first to introduce a method of analysis in which the interface between two uniform parallel streams of an inviscid, incompressible fluid is approximated by an array of line vortices. Using similar calculations, Abernathy and Kronauer [l] simulated the von Karman vortex street and Westwater [71] studied the rolling-up of a vortex sheet of finite breadth. In this classical line vortex method, it is assumed that the vorticity is concentrated at a number of points. This corresponds to approximating the vorticity by a sum of delta- functions. Each line vortex is moved by the velocity field induced by the other vortices. I . “‘3 l we iqtpw on“ v in V U! 'I’ if z»... “a i l ‘5 A... mL‘ “‘5‘ V‘s-h '- t“ W. L ,3 I In ’ '5. v' It; pr) 9 With the advent of computers, Birkhoff and Fisher [10] reexamined and refined Rosenhead’s calculation by following more line vortices with more accurate integration time schemes. They showed that some of the vortices actually move about irregularly in contrast to Rosenhead’s result. It was then apparent that the velocity field becomes unbounded near a line vortex, and this leads to a spurious interaction of neighboring vortices. This effect is not present in the original calculations by Rosenhead and Westwater, possibly because of the small number of vortices used or the limited accuracy of their calculations. Takami [69] and Moore [52], using a large number of vortices, indicate that the classical line vortex method is unreliable. Most of the earlier efi‘orts were directed toward understanding the time evolution of finite-area vorticity regions or initial break up of the shear layer in an inviscid, incompressible fluid. Due to the difficulty with the velocity divergence near the line vortex center, numerical solutions using line vortices have had some poor convergence results. In the early part of the last two decades, Chorin and Bernard [19] initiated a resurgence of interest in the use of vortex methods, for modeling real high Reynolds number flows of practical significance, when they suggested the use of discrete vortices (modified line vortices with finite cores) for smoothing the singularity associated with line vortices. The finite core radius is analogous to the introduction of a small viscosity which allows the vorticity in the line vortex to diffuse only a small distance. The discrete vortices are assumed to retain the same shape for all time. Chorin improved the convergence of flow past a circular cylinder [16] by an appropriate choice of the core or "cutoff" radius. The inviscid component of the real flow is solved via the discrete vortices while the viscous effects between vortices are 10 incorporated using a random walk for each vortex. Chorin [17] further improved the slow rate of convergence near solid boundaries when he introduced a vortex sheet algorithm to generate finite vortex sheets only as needed to . satisfy the no- slip boundary conditions in the surrounding potential flow. The rate of convergence for the various forms of the vortex core smoothing function have been discussed by various authors including Hald and Del Prete [39], Milinazzo and Saffman [50] and Beal and Madja [8] Other investigators have performed discrete vortex calculations on a variety of flow problems. Ashurst simulated turbulent mixing-layers [2,3], fluid motion in a form-stroke piston cylinder [4] and various combustion problems [5]. Shestakov [64] derived a hybrid method for the cavity flow problem. Men g and Thomson [49] combined the discrete vortex method with the VIC technique to study a class of non- linear hydrodynamic problems which includes the transport of aircraft trailing vortices in a wind shear, and the gravity current. A literature search shows that the majority of the calculations using discrete vortices have been performed for external flows. The method has not been extensively applied to problems involving internal flows. A comprehensive review of the general class of vortex methods which should serve as a good entry point into the literature is given by Leonard [45]. Other notable reviews appearing elsewhere include Clements and Mull [21] and Zabusky [74,75]. C31: 1 A. e “a... U bu en. mi 1 3 £1:- E J to}. iii-5‘ £5 1 1 1.2 Mixing Models Contemporary mixing models which utilize the turbulent structure of the flow in engines is perhaps best illustrated by Spalding [66] in his ESCIMO theory. The acronym stands for Engulfment, Stretching, Coherent, Inter-diffusion, Moving Observer. In this model, Spalding describes the engulfment process of the eddies as well as the details of the internal processes within an eddy. Although the model appears to conceptually embody much of the physics of turbulent mixing, validation of the model is still at an early stage. The Coalescence and Dispersion (CID) model, which admits finite rate micromixing, was first suggested by Curl [26] when he proposed a rather complicated integro-differential equation for which no analytical solution is known. Curl’s CID model was originally intended to describe drop interactions in a two-liquid phase chemical reactor. It was later extrapolated to gaseous turbulent systems and utilized in the fields of chemical engineering and turbulent combustion [35,36,57]. The model suggests a simple mechanism for the exchange of scalar concentration between initially segregated fluid particles that are fed continuously into a perfectly-stirred reactor. The feed rate was set to a predetermined rate so that a constant particle population was maintained. The fluid particles are pictured as undergoing independent pair collisions, representative of the molecular diffusion time scale. Levenspiel and Spielman [46] first used the Monte Carlo numerical technique to solve the integral term in Curl’s equation. The Monte Carlo method is a stochastic fluid particle interaction model of the turbulent mixing process. It is conceptually 12 simple and especially suited for mixing in discrete systems. The flexibility, clarity and practicality of this formulation has resulted in its subsequent widespread use. Flagan and Appleton [36] give a lucid formulation of the Monte Carlo technique for the solution of Curl’s equation. The works of Shreider [65] and Hammersley and Handscomb [40] describe various applications of the Monte Carlo method. The probability density function (PDF) model plays an increasingly significant role in the theoretical treatment of turbulent flows and turbulent mixing of scalar fields. Because transport of the PDF is accounted for in both position and probability spaces, these equations give a more complete description of the system than their moment counterpart and the physics of the system can be more exhaustively explored. The PDF formalism is particularly attractive for reacting flows, since the effects of reaction appear in closed form irrespective of the complexity and non- linearity of the reaction scheme. These favorable attributes have led several authors to base turbulence closures on transport equations for the PDF’s. Lundgren [47] devised a method to derive a hierachy of transport equations for N-point probability density functions for fluctuating variables in turbulent incompressible flow. The Lundgren approach has been applied to turbulent flames by Pope [55] and J anicka, et. al. [42]. Fox [37] and Dopazo and O’Brien [30,31] extended this method to compressible and reacting flows. The attraction of the PDF approach appears to diminish when solving the transport equation. Analytical solutions have been obtained in a few simple cases, but in general, numerical methods are required. Such numerical methods (e.g., the Monte Carlo method) have to cope with the integro-differential nature of the an 45‘ 1 about“ I not w" - SG..3C 3 1.; A v $45 A 13 equations and with the large dimensionality of the PDF when many species are involved. An order of magnitude analysis [56] shows that the finite-difference solution of the PDF transport equation is impractical for a number of species greater than 3. The computational expense is found to rise exponentially with the number of species. CHAPTER 2 DESCRIPTION OF PROBLEM The numerical simulation of air-swirl generated inside a flat disc-shaped cylinder representing the compression volume of a piston en gine combustion chamber is described in this chapter. Under the effect of a high pressure differential, a fuel jet is introduced at an adjustable angle into the cylinder through a nozzle orifice located on the cylinder wall and the problem of the evolution of the velocity and vorticity fields due to bulk diffusion is addressed. For liquid fuel injection into the cylinder of a piston engine, an aerodynamic force will be encountered due to the swirling air in the cylinder. The fuel, internally ttu'bulent and non-uniform at the exit plane of the nozzle, will exhibit a tendency to disintegrate that ultimately will lead to the breakup of the fuel into droplets of different sizes and concentrations in the spray. The distribution of droplet sizes will depend on the relative magnitudes of the aerodynamic, viscous, inertial, and surface forces present. The atomized fuel spray will be deflected by the air swirl and will evaporate and mix with air at the fringe of the jet at a rate determined by the primary air swirl, turbulent motion, and molecular viscosity. At some time during fuel atomization and mixing, combustion would occur with the ensuing flame spreading l4 15 across the chamber. The time between injection and attainment of chemical reaction conditions is the physical delay period. Lyn [48] observed that the combustion process can be divided into three distinguishable heat release phases, as shown in Figure 1. The first phase is associated with a period of rapid combustion and rapid cylinder pressure rise and lasts for only about 3° .crank angle (CA). During this phase, fuel that has evaporated and mixed with air during the physical delay period burns with a non- luminous flame. The burning rate decays in the second phase which lasts about 40° CA and the jet burns with a luminous, highly turbulent diffusion flame. In the final phase, the remaining .fuel burns at a very low rate extending through the entire expansion stroke. The present study is concerned with the mixing—dominated physical delay period prior to chemical reaction. The addition of combustion chemisn-y will be the topic of a future study. Experimental studies of the fluid flow in a firing engine cylinder prior to and during injection, reveals a complex combination of fluid motions. The nature of the motion in the engine cylinder can be broadly classified into three types : swirl, turbulence and squish. Swirl denotes a large scale rotary motion of the charge in the combustion chamber about the chamber axis, more or less. It is created during induction and amplified during compression by the combustion chamber configuration. Turbulence is the fluctuating velocity component superimposed on the mean velocity 16 of the fluid motion. Squish denotes the organized radial motion created by the piston movement near t0p dead center (TDC), that is, near the end of compression. Because it is anticipated that squish effects will significantly alter the flow patterns at the end of compression, a right cylindrical disc-shaped cylinder with no squish zones is used to model the combustion chamber near TDC. The formulation of a mixing model, which neglects bulk fluid motion typified by the air swirl inside the cylinder cannot hope to realistically predict design criteria. The theoretical formulation will consider a system of equations composed of the conservation of mass, momentum and scalar species. The problem of evolution toward a uniform final state of the fuel undergoing inviscid bulk diffusion at constant pressure and temperature will be numerically investigated using discrete vortex dynamics. At high Reynolds number, the aerodynamics of the flow are dominated by inertial forces except at the immediate vicinity of the solid surface where viscosity plays a dominant role. In this study, the effects of molecular viscosity, turbulence and molecular mixing are neglected for simplicity. The closure of the turbulent transport term and a satisfactory model for the molecular mixing term in the scalar conservation equation are particularly challenging. A future study will attempt to incorporate these effects. A discrete jet model analogous to Ashurst’s scheme [4] and a Vortex-In- Cell solution technique will be used to approximate the fuel injection and the fluid motion inside the cylinder, respectively. Because of the three-dimensional and turbulent nature of the phenomena occurring within the engine combustion system, it is not advisable to solve the fundamental equations even with the use of the largest computers. To keep the l7 mathematical difficulties involved within reasonable bounds, assumptions must be introduced to allow the simplification of the equations and the evolution of a tractable mixing theory: (i) It will be assumed that the controlling features of the problem can be described by a two-dimensional, constant volume, mean flow formulation. The flat, disc-shaped configuration of the cylinder and the fact that combustion in an internal combustion engine occurs near TDC where the volume does not change appreciably, justify this assumption. (ii) For two-dimensional flows, the stretching and rotation term in the vorticity transport equation is zero. The only nonzero component of the vorticity vector will be normal to both the velocity vector and the direction of the piston motion. (iii) Under engine operating conditions, the chamber conditions are near or above the thermodynamic critical point of the fuel. Test results by Shahed, et. al. [63] show that a major portion of the spray appears in vapor form. Droplets appear only near the jet nozzle exit within a few orifice diameters. This finding agrees with the conclusions of several other workers [13,48,53]. For simplicity, this study represents the fuel jet as air jet. The consideration of a single phase eliminates complications caused by fuel droplets. 18 (iv) The change in temperature and pressure of the vaporized mixture in the physical delay period is minimal. This assumption can be justified from energy conservation analysis. Thus the kinematic viscosity is assumed constant since the pressure is assumed constant prior to combustion. (v) We will also assume that, to a first approximation, the atomized fuel spray, idealized as discrete vortices, moves with the entrained flow. This is a fair assumption because parametric estimates [59] of the aerodynamic forces on typical droplets in a fuel spray suggest that the relative velocities of droplets in air rapidly decay so that droplets move with the entrained air in the jet flow. (vi) An incompressible fluid is assumed since there is no volume source or sink in the flow situation considered. Volume displacement effects are neglected since the volume occupied by the fuel is small compared to the volume occupied by the swirling air flow. (vii) Gravity is not a factor. Centrifugal convection is due primarily to the centrifugal body force which is accounted for in the total derivative. CHAPTER 3 THE DISCRETE VORTEX METHOD 3.1 Introduction Modern experimental investigations identifying the existence of large scale coherent structtues (eddies) in turbulent flows have rekindled much interest in the numerical modeling of turbulent combustion problems. Contemporary numerical models seek to utilize the observable eddy structures of the turbulent flow to simulate the physical systems directly without recourse to conventional finite difference formulation. Vortex methods represent such a class of numerical techniques. The central theme of Chorin’s discrete vortex method is based on the realization that it is the interaction of the eddies that produces transport of mass, momentum and energy. The basic concept of the method for flow simulation is to raprescnt the observable eddy structures in the vorticity-containing regions of the flow by discrete vortices and then track this discretization in a Lagrangian reference frame. The essentially grid-flee formulation permits the simulation of the growth and decay of eddies encountered in real fluids. The method is particularly appealing 19 20 because it is specifically designed to deal with problems of increased spatial resolution, stability and numerical viscosity associated with finite difference formulation. 3.2 The Flow Field At any point in a non-reactive mixing flow system, the state of the fluid can be characterized by the velocity V, the pressure p and a set of scalars = _a_(1. 833).) _ .02_ + vi' (3,47) The averaged joint PDF PC?) can be defined as the ensemble average of the instantaneous joint PDF p(\I'), a random function due to the dependence of the distribution of 4’01 upon the initial conditions ca (x,0), and the random forcing turbulence field V(x,t). Averaging over all possible realizations of initial conditions and V fields yields a a PC?) = H < P(Va;x.t) > = H < 6 (Va - ¢a(xit)) > (3.48) a=l a=1 where the angled brackets denote an ensemble averaged quantity and 8 is a Dirac delta function. The justification for defining a probability density as the ensemble average of a collection of Dirac spikes can be found in Stratonovich [67]. 41 It is interesting to observe the correspondences between the physically recognizable terms in equation (3.45) for the fluctuating scalar quantity ¢a and equation (3.46) for the averaged joint PDF PCP). The convective term retains its structure in the POP-equation but the turbulent convective contribution is added. The molecular transport term appears in the PDF-equation twofold : as a transport term in position space of the same structure as in the (pa-equation, and as a transport term in probability space. Molecular transport in position space represents molecular diffusion while molecular transport in probability space represents mixing by molecular action. It is shown in [54] that the effect of the mixing term is to reduce the variance of PC?) without affecting the mean. Omitting the molecular diffusion terms, which are negligible at high Reynolds number, equation (3.46) reduces to [ML 4. a—PSZL] =- 3p (vi p(‘I’) > - .32— (3.49) p at 1 axi axi a\l’2 aXi aXi The transport equation in this form is not immediately soluble, because, whereas the term on the left is exact, the terms on the right, representing turbulent transport and molecular mixing, are unknown and need to be modeled in order to close the transport equation for POP). 42 The closure of the turbulent transport term is the fundamental stumbling block to progress in probability formulation of turbulent scalar mixing. This is because turbulent mass transport or eddy diffusion is a complex process that is dependent on the scale and intensity of turbulence in the flow. Various models, which tend to make the mathematics easier but at best are a crude representation of reality, have been suggested in the literature [28,42,54] to describe the phenomenon. Pope’s [55] simple gradient diffusion model, at) _a_ 3P9?) axi - - Bxi I‘t axi ’ (3°50) introducing the turbulent diffusion coefficient I‘t(x,t), seems promising and better behaved for inert scalars and for large departures from Gaussiarrity and homogeneity. As with the turbulent transport term, a number of models for the molecular mixing term have been proposed in the literature [26,31,42,55] but a completely satisfactory model is yet to be demonstrated. The fluid particle coalescence and dispersion model, proposed by Curl [26], is the only model applicable to an arbitrary number of scalar quantities and can be expected to produce qualitatively correct results. The model makes no attempt to describe the detailed dynamics of the turbulent flow but, instead, relies on a probability density function of concentration, itself dependent on empirically determined mixing frequency, to describe the non- uniformities in the flow. Flagan and Appleton [36], Evangelista et. al., [35], and Pratt [57] have modeled the mixing term in the form 43 - 3f? 0. A numerical solution for POP) can be attempted by treating the simultaneous action of the advection, diffusion and mixing processes sequentially using the technique of operator splitting mentioned in Section 3.2. The component equations are Advection : ayatfli + (Vi)? = 0 (3.53) Edd an 2 8P9?) Diffleionz p at = 32;; (rt 88 ) (354) miti‘g'i” p A1313: = - pptv) + 263 f P(‘P+‘P’)P(‘P-‘I”)d‘1” (3.55) The advection of POP) by the mean flow (i.e., bulk diffusion) is obtained using the solutions fi'om the analysis of the flow field discussed in Section 3.2. It is assumed that, to a first approximation, the motion of the injected scalar particle is identical to the motion of the Lagrangian point of fluid that would occupy the position of the particle if it were not there. That is, the conserved, passive, scalar field is advected by the fluid. Turbulent difiusion of PCP) can be approximated by adding turbulent velocity components to the droplet velocities. This requires that the turbulence intensity and length scales be specified. This point remains a formidable closure problem. The use of random walk in modeling turbulent diffusion was suggested by Chorin [15] as a generalization of the random walk model used in his vortex sheet method for viscous 45 diffusion. Molecular mixing of PC?) can be approximated using the Monte Carlo method, mentioned in Chapter 1 for the simulation of the coalescence-dispersion term. An essential feature of the Monte Carlo computional technique is that at any location in the combustion chamber, the joint PDF is represented by an ensemble of fluid particles which are identified in terms of their individual thermodynamic states. The turbulent mixing process is described as a sequence of interactions between randomly chosen groups of n representative fluid particles which mix (with frequency [3) completely with one another to form it new particles having an identical thermodynamic state. The simulation is allowed to continue until the ensemble statistics (mean and variance of properties) are observed to stabilize. Finally, the ensemble average and other statistics are obtained to represent the properties of interest. The Monte Carlo method is insensitive to the dimensionality of the PDF and the computational expense rises only linearly with 6, which is the best that can be achieved by any algorithm. The combination of the solutions for the advection, diffusion and mixing equations produces an approximation to the joint PDF solution for the turbulent mixing problem. However, because of problems posed by closure approximations, the present study focuses primary attention on the numerical solution of the advection or bulk difi‘usion component equation. Eddy diffusion and molecular effects will be incorporated in a future study. 45 3.5 Choice of Numerical Parameters and Functions Specifying threshold values for the following nine control parameters : (i) d, mesh size, (ii) dj, injector nozzle diameter, (iii) r0, discrete vortex core radius, (iv) h, boundary discretization length, (v) At, length of time step, (vi) 58, numerical shear layer thickness, (vii) 7m, maximum allowable vortex circulation. (viii) 1, initial number of "smaller" vortex sheet segments created at wall test points, (ix) [3, mixing frequency, and the proper choice of a vortex core smoothing function Gij and a numerical filter function (p will produce a discrete vortex calculation which describes a real flow. If the evolution of the important features of the real flow are to be predicted, then an appreciation for the interplay between these control parameters must be developed. Details of the numerical implementation of the foregoing theoretical analysis are presented in Chapter 5. The numerical experiments germane to the present study and discussion of results are presented in Chapter 6. CHAPTER 4 NUMERICAL FILTERS 4.1 Introduction q Numerical filters serve a dual purpose in Vortex-In—Cell (VIC) formulation. They are used for vorticity deposition on the grid and for interpolation of local effects of the velocity field from the grid to the Lagrangian points. These interpolation procedures are subject to errors. The extent to which these errors are minimized is dependent upon the design of the filter. The choice of a carefully constructed filter for use in the VIC algorithm is therefore a critical factor in obtaining consistent numerical results. This choice is dictated by several considerations including accuracy, grid independence, a reduction of undesirable grid effects or aliassing, Fourier transformability and computational efficiency. 4.2 Description of Filters A survey of the literature shows that there have been essentially five methods used to derive the proper numerical filter. Some have assumed a "shape" for the vortex core in real space while others have generated the shape in wave vector space to obtain some desirable features during the Fast Fourier Transform (Fr-T). A sixth 47 48 filter developed by Bartholomew [7] is included in this description. An important consideration is that the remainder of the fluid not sense that an interpolation has occurred, hence, the vorticity carried by each vortex should be spread to as few mesh points as possible. The Cloud-In-Cell (CIC) method, first reported by Christiansen [20] and later used by Baker [6], assumes each point vortex to possess a uniform or "top hat" vorticity distribution within a square cell. The distribution involves the nearest two grid points in a given coordinate direction. Meng and Thomson [49] used a slightly modified "top hat" profile of Christiansen. They assumed a vorticity distribution over a rectangle with dimensions x and y of one cell according to §(x.y) «4 1 - 1 (4.1) 3(x-x-) 2 3y-y-) _.L _.L [ 1 + Ax ] I: l + Ay :(2 where Ax/3 and Ay/3 are the half-widths characterizing the spreading and xj, yj are th the coordinates of the j vortex. Like Christiansen’s filter, the distribution involves the nearest two grid points in a given coordinate direction. To improve the numerics of the . CIC method, Hockney et. al. [41] proposed the use of a truncated-Gaussian core having a larger cross-section (2d x 2d) for the 49 vorticity distribution function : I m .2 50%) .. exp {- EL } lxil s d =1 2 2r0 (4.2) =0 Ixi|>d where 1 s i S m, r0 = 0.455d, d is the distance between mesh points, xi is the coordinate in the ith dimension relative to the vortex and m is the number of dimensions. This distribution involves the nearest three grid points in a given coordinate direction. That is, the Gaussian distribution is truncated beyond 1.5 grid lengths. Couet [23,24] and Wang [70] developed their Gaussian distributions for the vorticity in wave vector space and then transformed them into real space. The results gave quadratic and cubic splines, respectively, and involve the nearest three and four grid points, respectively, in a given coordinate direction. Bartholomew [7] assumed a polynomial expression to describe the fraction of the total circulation credited to the nearest three grid points in a given coordinate direction. These distribution functions can be put into a common format so that in two dimensions the filter is the product of two one-dimensional filters. These one- dimensional filter models, the scheme used in generating them, the number of nearest 50 grid points involved for each scheme and their authors are listed in Table 1. It is from Bartholomew’s perspective of using the definition of circulation that the compact forms of the filters, with the exception of Couet’s and Wang’s, were independently derived. Profiles of the six filters are illustrated in Figure 3. It is observed that the profile of Meng/Thomson’s filter is in fact a slight modification of Christiansen’s triangular profile. The others appear to assume a quasi-Gaussian profile which is brought smoothly to zero at the extent of each filter. A generalization of the formulations listed in Table 1 can be made as follows : A uniform rectangular grid system having coordinate directions x and y is first defined as shown in Figure 4. Next, an assumed vorticity distribution function §(x,y) = f(x,y). (4.3) with the point vortex coordinates (xp,yp) taken to lie at the center of this distribution, is normalized, C(Cm), where C=x/d and n=y/d, and a local grid point which serves as the origin for the coordinates of the discrete vortex is chosen. The normalized filter tp(C,t|) which determines the fiaction (71,)”; of the total circulation 71, that is credited to the surrounding grid points of interest is now obtained by integrating the normalized vorticity distribution over the overlapping area in the cell enclosing the grid point : (Cm), involves an elaborate four—step procedure: (i) The functions s(C,k) and 31(11 ,k) are pretabulated. (ii) Equation (5.1) is reduced to another system of non-linear, algebraic equations of the form (see Appendix C) f(k,s,sl,X) = O (5.5) f(krssSPY) = 0 (5.6) (iii) Equations (5.5) and (5.6) are then solved simultaneously for s and SI, for every x and y input, using a quasi Newton-Raphson numerical method. 59 (iv) The coordinates (Cm) are recovered floor the pretabulated functions s(C,k) and sl(n,k) via linear interpolation. 5.1.2 Velocity transformation Let a fluid motion inside the square region of the (u-plane be given by the complex potential, W(m) = q, + hp (5.7) where e is the velocity potential and \y is the stream function. Then at corresponding points z and 0), given by Eq. (5.1), W and therefore it and ‘1’ take the same values. Now, 8’ is a boundary and so a streamline, and therefore 1]] = C, a constant, at all points of 8’. Since S corresponds point by point with S’, u! = C at all points of S. Therefore S is a streamline, in the motion given by (5.7) and (5.1) together, in the z-plane. The actual form of the complex potential in terms of 2 would be obtained by eliminating (0 between (5.1) and (5.7), but it is often preferable to look on 2 as a parameter and forgo the elimination. Thus, to find the fluid velocity at (x,y) in the z- plane corresponding with (Cm) in the (D-plane, we have d_W_ d_VY.d_w dz ‘ d0) dz (5'8) 60 The equivalent matrix form of Eq. (5.8) can be written as [3] _ c11 e12 11] (x,y) ‘ [W ”11 " (5'9) (Cm) where en, cm are elements of the 2 x 2 matrix. The velocity components in the z- plane can be transformed to give velocity components in the cylindrical coordinate system as : V9 = vcosO - usin0 (5.10) V = -vsin0 - ucosO (5.11) Conversely, the fluid velocity at (Cm) in the (tr-plane corresponding with the fluid velocity at (x,y) in the z-plane is given by dW dW dz m=arrs (m) 01’, 1°] [in “1 = 5.13 V (Cm) ’12 ell " (x,y) ( ) Notice that the matrix system in Eq. (5.13) is the inverse of that in Eq. (5.9). 61 Representation of the grid velocity field in the circular domain is done via conformal transformation as discussed above. Off- grid point velocity calculation in either domain can be accomplished using an interpolation routine relevant to each grid system. 5.2 Point Vortex Initialization Scheme The irrotational vorticity field at the initial time can be approximated using discrete vortices in two ways. The first scheme involves analyzing the problem from the point in time when the mechanism that is responsible for the production of vorticity within the boundary layer first introduced vorticity into the interior domain to the point where the flow field matches the desired initial conditions. This method is computationally expensive. It would be computationally more advantageous to develop a method that generates the initial distribution and strengths of the vortices from a prescribed velocity field directly, or alternately, from the derivative of the velocity field. In either case an iterative method is implied. The initialization scheme employed in this study makes use of the tangential velocity field prescribed by the experimental data of Dyer [33]; see Table 2. The data is classified according to a time to which is the time from the closing of the shrouded intake valve to the time of measurement. Dyer’s data were used to represent different initial swirl conditions for this study. The strengths and optimum distribution of the discrete vortices are determined as follows : 62 (i) Discrete vortices are arbitrarily but uniformly placed in concentric circles about the origin of the circular domain. (ii) The modified circle theorem of Milne-Thomson is applied to the system of vortices to compute geometric factors and using Dyer’s velocity data, a matrix inversion routine supplied by IMSL is used to determine the strengths of the vortices. (iii) The optimum distribution is found when the sum of the circulation due to each vortex is equal to the global circulation 21tRVe—, i.e. N N 2 (719,- =21tR_V9_ or 2 K,- = Rig (5.14) i=1 i=1 where N is the number of discrete point vortices. (7p)j is the circulation of vortex j, 1% is the strength of vortex j, R is the radius of the combustion chamber and-VB is the average tangential velocity or slip velocity obtained by numerically integrating Dyer’s experimental velocity data. To obtain the strengths of the discrete point vortices, consider an arbitrary distribution of the vortices as noted above. The velocity at any point zk inside the circular domain is the sum of the velocities due to all vortices located at 2]- given by 63 the modified circle theorem: N [u - i v ]k = 21 -itq [ 1/(zk - zj)-1/(zk - aZ/z‘jn (5.15) 1... Equation (5.15) is composed of three factors : the desired velocity (u—iv)k, the discrete point vortex strength Kj and a geometric factor relating the influence of vortex j at point 21‘. This information can be represented in matrix form as where (V e)k’ ij, Sj are matrices for the prescribed tangential velocities of the swirl at the points zk, geometric factors and vortex strengths respectively. Solving for the strengths by inverting the ij matrix, we get _ -1 Sj — ij (Ve)k (5.17) See the Appendix D for the development of the matrices. 5.3 Vortex-In-Cell (VIC) Algorithm A VIC technique is essentially a formulation that combines the Lagrangian vortex tracking concept with the stream function-vorticity finite difference 64 formulation. A distinguishing feature of the VIC method is the fact that vorticity information is contained in the Lagrangian vortex element rather than on the Eulerian mesh. The Eulerian mesh is introduced for fast Fourier transform and convenient representation of results, but is not explicitly used to advance the solution in time as would be done in the finite difference method. As a first step towards solving Poisson’s equation (Eq. 3.8) on the square of fixed mesh size d, we must create, at the beginning of each time step, a mesh record of the vorticity field starting fi'om the Lagrangian representation. This is done by an interpolation scheme using Couet’s filter 2 .thlnl] ; é—slnlsg- (Pm) ={r} ~n2 ; In |s-2L} (5.18) 0 ; Otherwise where 11 = y/d or x/d, to distribute the vorticity from each vortex over the nearest three grid points in a given coordinate direction. Figure 10 illustrates an application of the above one-dimensional model of Couet’s filter. Three quadratic spline distributions of vorticity with different total amplitudes are shown to have their centers located at positions a, b, and c. The three nearest grid points to each of the centers will share the vorticity distribution according to the spline function weighting on them. Grid points N-2 to N will share the vortex "a" with grid points N-l receiving the largest weight. Similarly, points N—l to N+2 will share the vortices "b" and "c". 65 After the vorticity information is represented on the uniform grid, a grid stream function can be generated by using FFI‘ techniques to invert Poisson’s equation (3.8). The stream function is centrally difl‘erenced on the grid to produce a velocity field which is conformally transformed onto an equivalent grid system in the circular domain. The velocity at the vortex position is now obtained by interpolation. The vortices are then convected to new positions and the procedure is repeated each time step until the flow evolution is completed. The convenience of FFI‘ is available only when boundary conditions are imposed over the square domain. In order to satisfy the no slip normal boundary condition, the value ‘1’ = 0 is specified on the grid boundary for this study. The grid used is a modest 60 x 60 but extension to larger dimensions is feasible with the use of mainframe computers. Swartztrauber, et. al. [68] have developed an elliptic boundary value problem solver (FISHPAK) which allows for general boundary conditions and is portable, i.e., it runs on a Cray, Cyber, IBM, or Vax computer with no modifications. This software is used in this study to accomplish the Poisson inversion. In general, if the grid is M x M, this is an order leogzM calculation. In a typical 2D calculation, the number of grid squares M x M is of the same order of the total number of vortices N and thus FFI‘ requires on the order of NzlogzN operations per time step. The interpolations take of the order N operations per time step. When N is large, NzlogzN is considerably smaller than N2. Thus the VIC technique offers a resolution of the problem associated with excessive computation time. 66 5.4 Vortex Sheet Algorithm The velocity field in the boundary layer region can be constructed using vortex sheets as described in Section 3.2.2. The u and v component velocities are given by Eqs. (3.26) and (3.27). At the initial time, the tangential velocity V6 along the wall is fixed by the potential flow solution of Eq. (3.8). Theoretically, the tangential no slip boundary condition, Ve=0, can be satisfied by generating a continuous sheet of vorticity along the wall with circulation 78 = — V9 per unit length of the wall. However, it has been observed computationally [2,18] that 78 = —V9 will satisfy the tangential no-slip boundary condition only in the limit as At —> 0. In order to satisfy the boundary condition exactly (except possibly at a finite set of values for t), Chorin [18], algorithmically, extended the velocity field across the wall (y=0) at the beginning of each time step by anti-symmetry u(x,-y) = -u(x,y). As a result C(x,-y) = C(x,y) for y at 0, and a continuous vortex sheet of circulation Vs = -2 V0 per unit length of the wall appears at y = 0 if the tangential velocity V6 does not vanish at the wall. This sheet is segmented into elements of finite length h, centered at test points along the wall, specifying the spatial resolution in the x- direction. In order to have a reasonable approximation of the diffusion equation (3.7), the resolution in the ydirection can be improved by replacing each vortex segment by 1 number of "smaller" segments with circulation (ys)j=ys/l such that |('Ys)j I S l¥nax' 67 where gnu is a predetermined value. Vortex elements from each segment h diffuse normal to the wall via random walks and induce a velocity field which advects them downstream At subsequent time steps, V9 is calculated using the obvious specialization of Eq.(3.26). The correction term (ys)j/2 in Eq. (3.26) was added by Chorin [18] to account for the effect of the image of the discrete point vortex the sheet can generate (see Figure 12). Unlike the "elliptic" flow modeled by point vortices, where the effect of each point vortex extends throughout the field, the zone of influence of a vortex sheet is restricted to the "shadow" below it, a consequence of Eq.(3.24). Thus, Eqs. (3.26) and (3.27) imply that the sum is over all vortex sheets whose "shadow" on the wall engulfs the sheet whose center is located at (xj,yj). This is illustrated in Figure 11. The darker the shadow, the smaller u(x,y) becomes. Whatever fluid enters a shadow region from the left and cannot leave on the right must leave upwards. To reduce the total number of sheets retained each time step, the vortex sheet algorithm can be implemented such that exactly half of them disappear. This is accomplished as follows : at each wall test point an even number of sheets are created and their intensity is adjusted such that if I (Yslil < IXnaxl no sheets are created. A rejection technique [18] is then used to ensure that any two successive 112’s used at the wall will have differing signs. This rejection technique can be used only at the wall, or else it would destroy the independence of the successive "I? ’s in the interior and thus fail to describe the diffusion process correctly. 68 A tagging and variance reduction technique can be used to effect a substantial reduction in the statistical error. Since diffusion takes place only in the y-direction, the random numbers (n2)j used in (Eq. 3.33) need be independent of each other only when they are used with vortex sheet segments whose centers lie in a narrow strip perpendicular to the wall. As vortex sheets are created, they are assigned integer th vortex sheet segment. All elements with tags, mj being the tag assigned to the j the same 112 are assigned the same tags. The effect of tagging is to piece together the elements created at the several boundary points into coherent vortex sheets with the elements of each sheet identified by a common tag. 5.5 Advection Scheme Vortices can be advected in only two ways, by convection or by diffusion. Once the local velocity is known, the convection step is accomplished by using an appropriate time integration scheme. For vortices in the interior region, Heun’s time integration method is used to evaluate the integrals in Eq.(3.l7) and Eq.(3.l8), namely n+1 _ n Atl n n Xj — Xj + 2 11(Xj ,Yj ) + u(xj*’yj*) ] (5.19) .n .n .n =1- * t3? I °* = on on on 0* = on on on where xJ xJ + AH: u(xJ ,yJ ) ] and yJ yJ + A2L[ V(xJ .YJ ) ] (5.21) 69 Heun’s technique is second order accurate in time. Such accuracy lessens the tendency for discrete vortices to spread out in spinning motion as they come close together. A first-order scheme merely takes the tangent at a given point and constructs a straight line approximation to the trajectory, in effect spreading out the vortices. Since Heun’s method requires two evaluations of the velocity field per time step, this can prove expensive when using the direct O(N2) method of calculating vortex interactions. The FFI‘ technique is advantageous in this respect. Near the solid boundary, an explicit Euler integration which is first-order accurate can be used to evaluate the integrals in Eq. (3.32) and (3.33) in order to accomplish the convection step, namely n+1 ___ on .n xJ xJ + uJ At (5.22) n+1 _ .n .n yJ - yJ + vJ At (5.23) The use of a second-order time integration scheme for the vortex sheets is not necessary for two reasons. First, the convection field within the boundary layer is essentially a one-dimensional flow in a direction parallel to the wall and thus poses no great problem for a first-order scheme. The main component driving the sheets away from the wall is the diffusion term. Second, the implementation of a second- order scheme is a very complicated procedure, since sheets may turn into point vortices during the predictor step and thus must be treated differently. 70 It can be assured that the vortices are not convected beyond one mesh width (for discrete point vortices) or sheet length (for vortex sheets) as long as the Courant condition IVmax IA} 5 dorh (5.24) is satisfied. The Courant condition sets an upper bound on At for numerical stability. Diffusion as discussed by Chorin [18] is equivalent to a random walk which is governed by statistical laws derived from the solutions of the diffusion equation (3.7). The integration of this equation yields a Gaussian distribution for the flux of vorticity [40]. Tschebysheff’s theorem [44] states that the Gaussianly distributed random variable with mean zero will rarely be greater in magnitude than three standard deviations. During a time step of duration At, vorticity will diffuse with a standard deviation of \l 2vAt . Hence, in a time step, interior vortices will not likely diffuse beyond one mesh width if 3% S d . In some physical sense, statistical simulation of diffusion steps back from the concept of a material continuum, and employs a sort of kinetic theory of giant and imaginary particles which wander randomly through the region of interest carrying macroscopic elements of the diffusing field. The outcome of a calculation is a list of particle positions and the elements of the concentration field they are employed to transport. 7 1 5.6 Sorting Algorithm Vortex sheets and discrete vortices are objects of a very similar nature; they are both determined by position and circulation. A computational element (xj, yj, Cj) can be treated as either a sheet or a discrete point vortex, depending on the circumstances. A sheet of negative circulation casts a shadow which slows the fluid under it; by the equation of continuity, this creates an upward flow to the left and a downward flow to the right, just as if the sheet were a point vortex. Conservation of circulation arguments show that a sheet segment of circulation Vs should transform into a point vortex of circulation yp = ysh, where h is the length of the sheet. These facts can be used to create a transition between the discrete vortices and the vortex sheets. Following Chorin’s argument [l8], vorticity is transferred from the boundary layer to the interior by allowing those vortex sheets located a distance greater than 58 from the wall to become point vortices. Any vortex which finds itself less than 58 from the solid boundary (inside or outside) becomes a sheet or is reflected accordingly (see Figure 12 ). Reflection of any sheet which crosses the wall back into the fluid imposes the antisymmetry condition. By taking 58 large relative to the variance of the steps, it is unlikely that a discrete point vortex will move further outside the domain than 53 in one time step; if this does occur, it is removed. 72 5.7 Mixing Algorithm The physics of the evolution of scalar field undergoing turbulent and molecular mixing is embodied in the PDF transport equation (3.47). The Lagrangian view of molecular mixing can be numerically simulated using the Monte Carlo coalescence and dispersion (CID) scheme. The initial bimodal joint PDF, P(\|!;X,0), can be characterized by the random distribution inside the combustion chamber of N equal-mass, chemically-inert fluid particles with values 0 designating the swirling air particles and 1 designating the injected air particles. The fluid particles are idealized turbulent eddies (i.e., finite- cored vortices) assumed to be larger than molecular dimensions, but sufficiently small with respect to the turbulence microscale. The thermodynamic state of each particle (i.e., pressure, mass fraction, temperature, etc.) is assumed to be uniform throughout its volume at any instant of time, but the number of particles in a given state may vary with time as a result of mixing between particles. The Monte Carlo C/D scheme is a stochastic fluid particle interaction model. Each fluid particle is assigned a number index, 1 S i S N, and a pseudo-random number generator is used to generate it different indices, i1, i2, . . ., in (ij at ik) thereby identifying the particles to be involved in the first interaction. In the simplest binary C/D mixing model (n=2) proposed by Curl [26] pairs of fluid particles are allowed to mix completely by adopting the average values of their properties and then separate. The mixing interactions are computed sequentially; the first interaction representing a 73 forward step in time equal to _ __1_ At(mix) - MN (5.25) where [3(t) is the empirically determined, time dependent mixing frequency and N is the fluid particle population in the combustion chamber. A new group of pairs of fluid particles is next chosen at random, from the entire ensemble of N fluid particles, to be involved in the next mixing interaction; this way the calculation proceeds forward in time inside the combustion chamber. At any time during the mixing process, the mean composition and other mean properties of the fluid may be evaluated by taking an ensemble average over N fluid particles. Furthermore, the shape of the evolving distribution function f(tion field induced by a single discrete point vortex of strength K=21t, placed cellh‘ally in the computational grid. The second examined the effect of a system of 105 ClisCI'ete point vortices (0.0 2.0d intmum M ooooooooooooo ooooooooooooo +++++++++++++ mmmmmmmmmmmmm ooooooooooooo .ooooooooooooo. Efiooooooooooooo, ooooooooooooo ooooooooooooo ooooooooooooo~ ooooooooooooo 0000000000000 Maximum HHHHHHHHHHHHH ooooooooooooo +++++++++++++ mmmmmmmmmmmmm ommmNthomumm mnmmmoabmmmmo EfiNmMNNmHmmvmom: MMHHOm¢MMFPPm omvmovmrmoamm HHommHmHmNmNo MMMNNNHNNMMVW mmmmmmmmmmmmm l intmum M 0000000000000 0000000000000 ooooooooooooo ooooooooooooo sfiooooooooooooo:fl ooooooooooooo ooooooooooooo ooooooooooooo ooooooooooooo 0000000000000 Cutoff Value oamm¢mmnomoam .Hmmmmmmmmmm¢¢”.u ;Hooooooooooooognl 0000000000000 4 Maximum ooooooooooooo ooooooooooooo +++++++++++++ mmmmmmmmmmmmm mammmflm¢mvmmo mommHHmnNmovH Efiflmmammaaommmm: MOOmemHWNmHH mmmmm¢mmmmwmm hmammmmmamabv mmmmmmhbmmmmo Q'VV‘Q‘V‘VQ‘Q‘Q‘Q‘Q‘Q‘ID For All Grid Points > l.5d intmum M ooooooooooooo ooooooooooooo +++++++++++++ mmmmmmmmmmmmm ooooooooooooo ooooooooooooo Efiooooooooooooo: ooooooooooooo ooooooooooooo ooooooooooooo ooooooooooooo 0000000000000 POINT VORTEX AT CELL CENTER For All Grid Points Maximum HHHHHHHHHHHHH ooooooooooooo +++++++++++++ mmmmmmmmmmmmm mm¢mmmombmmma mmmmmmmmmmbwm ffimmmmmmmmmmavv: NH¢mommh¢mmmm havommmhvmwmo mmmmmv¢ormwoa mbmmvvvmwmavm NNNNNNNNNNMMM A intmum M 0000000000000 0000000000000 +++++++++++++ mmmmmmmmmmmmm 0000000000000 0000000000000 gfiooooooooooooogH 0000000000000 0000000000000 0000000000000 0000000000000 0000000000000 Cutoff Value mmwmmbmmoamm¢ vvvvvc¢vmmmmm _ gfioooooooooooooih 0000000000000 A 109 Table 4. Error statistics for grid point coordinate and velocity transformations. Coordinate Transformation Velocity Transformation C 11 u v Avg. absolute error 2.0397014507 1.4282386507 1.5206969501 1.4264177501 your. absolute error 1.3470206505 6.3128371506 3.4926014E+00 4.0869408E+00 lMax. relative. error (76)) 1.9993558502 1.0198522502 4.09980205+07 9.25844255+05 Root mean square error 6.8885259E-O7 35780748307 3.3185267E-01 3.0976245E-01 Grid point location (2,1), (2,61), (1,3), (159), (2,1), (2,61), (1,3), (159), °f m‘ mu” m (60,1), (60,61) (61,59), (61,3) (60,1), (60,61) (61,59), (61,3) Note : Average absolute error = Relative error (%) Root mean square error In M [‘2 [2(1’5 ' PN)2/M2T i=1 j=1 where p is the parameter of interest. (6.6) (6.7) (6.8) 110 55.5 55.5 55.5 555 55.5 555 5 55 55 55 5 5 5 5 555 55.5 55.5 55.5 55.5 55.5 555 5 55 55 55 55 5 5 5 555 55.5 55.5 55.5 55.5 55.55 555 5 55 55 55 55 5 5 5 555 55.5 55.5 55.5 55.5 55.55 555 5 55 55 55 5 5 5 5. 555 55.5 55.5 55.5 55.5 55.55 555 5 55 55 55 5 5 5 5 555 15%. .1 5::— .55__ £5:— .5¥_ .55. N 2 5.2533550 58:55 5505— A 0 ”805° 0 " one a c I e o o o 0 555.55 Maw ".55 mm 5 55 55 T55 _ 55 _ 55 1 55 _ 55 55 255 £8.55 om.m_ n .58 N 3555—55 2...: £5.85 mm. 5 u .55: 5 a} £8.85 5 u 5:. ~ .555 555 v 5p v 55 8: 5.55 85555555 5.55 558.; .355 6.55.5. 55.5 55.5 55.5 55.5 55.5.5 .555 5 .55 55 55 5 5 5 5 55.. 55.5 55.5 55.5 55.5 55.55 555 5 55 55 55 5. 5 5 5 55 55.5 55.5 55.5 55.5 55.55 5.55 5 55 55 55 5. 5 5 5 55 55.5 55.5 55.5 55.5 5.5.5.5 555 5 55 55 55 5. 5 5 5 55 55.5 55.5 55.5 55.5 55.55 555 5 55 55 55 5. 5 5 5 5 55.5 55.5 55.5 55.5 5.5.55 555 5. 55 55 55 5. 5 5 5 5 55.5 55.5 55.5 55.5 55.55 5.55 5 55 55 55 5. 5 5 5 5. 55.5 55.5 55.5 55.5 55.55 555 5 55 55 55 5 5 5 5 5 :5 55.5 55.5 55.5 55.55 555 5 55 55 55 5 5 5 5 5 55.5 55.5 55.5 55.5 55.55 555 55 55 55 55 55 55 55 5 5 55.5 555 55.5 55.5 55:55 555 55 55 55 55 55 55 55 5 5 55.5 55.5 55.5 55.5 55.55 555 55 55 55 55 55 55 55 5 5 55.5 55.5 55.5 55.5 55.55 555 55 55 55 55 55 55 55 5 5 15555. 1 5515; 5515; .55. 5 2 8555555 55> .555 .O 51 -AENV u 555 5E 55.55 n .P .5 5.5 u 5 5.5 5.5 _ 5.5 _ 5.5 _ 5.5 _ 5.5 5.5 5.5 an. . . _. $5.5": 55. 555555u55vr55u5o58535 .55 .555 85 u 55.5 555 555525 5555555555 5...... 558.55 .35 6.55.5. Table 6. L, 6. optimization data (Swirl strength S6 = 12.5 cm-m/s). 111 Nc = 2 ej = 30° Grid Size U“ 6 x 6 8 x 8 10 x 10 (cm/s) (L=32) (L=52) (L=80) t’ (1118) 112.3 173 185 344 125.7 149 229 327 152.6 112 213 319 210.0 124 184 250 NC = 2 Bi = 60° Grid Size Uin 6 x 6 8 x 8 10 x 10 (cm/s) (L=32) (L=52) (L=80) t' (m8) 112.3 185 207 343 125.7 159 198 307 152.6 130 213 304 210.0 99 198 343 Nc = 2 Bi = 900 Grid Size Uin 6 x 6 I 8 x 8 I 10 x 10 (cm/s) (L=32) (L=52) (L=80) 6" (ms) 112.3 215 267 356 125.7 191 236 371 152.6 138 268 373 210.0 108 167 330 112 Table 7. Ne optimization data (Swirl strength S9 = 12.5 cm-m/s). L = 32 9j = 60° N Uin c 1 2 3 (cm/s) 1‘ (ms) 112.3 179 185 185 125.7 128 159 173 152.6 113 130 167 210.0 82 99 152 Table 8. Optimization data for swirl strength Se. N,= 2 Swirl strength 89(cm-m/s) Uin 4.1 7.9 12.5 17.5 . 27.8 (cm/s) t‘ (1118) 112.3 175 178 185 152 146 125.7 210 135 159 150 128 152.6 141 164 130 123 149 210.0 116 88 99 100 118 Rate of Heat Release 113 Figure l. Crank Angle Three phases of heat release. 114 >X . Fuel Oj Angle of injection Figure 2. Sketch of intake flow 115 5.505. cmmwcfiws on. o. 2% on. 5:85 So: 55:00—53 2: .50 5.06550 .82 mo coca—55825 .58 55:5 2% 65 5 55:62.25 55:35:; 58 53:51. Hutu—55:2 .m Semi 5823..» 0.558 552.5: .5 oi 0.» ad 0.. 06 red 5 rad 1:. w ( r 1.6 5 :3 1 1 11 1 11 .111 I111 155 5823.55 055.935.5528 .5 0... ad ad 0.. 0.0 b 1% L I-l h L1 1? ".0. 106 1a.: .3 w, ( 1 1.6 W 1.6 1 11111 a. 82:28; 65 855.83 2830:5525: 5.5 6.» ed a; o... L L «d1 '31 mo 11--- 15.. 825.8; 65. 53255.2:an .5 0.9 ad 0d a. - ed 10.6 10.0 («)4 rod 111111111 1111 o.— - 1: - 5535598 85225 62x8: 5. ad 0.» ad 9.— ad 1 I? L L («)4 j (0)6 116 11111111 o 0.0.. o o o / '0’... 1 I 0.0.0. . . C 0.0... 0.0.0. .mo 0 Figure 4. Vorticity distribution in a grid, with shading showing assignment of vorticity to grid point for VIC method. is assigned to grid point (i,j) @ is assigned to grid point (i+1.j) 111111 is assigned to end point (i+1.j+1) g is assigned to grid point (i,j+1) 117 © WW (POLYNOINL) WSW/BAKER (AREA wacumc) COUEY (OUIDRATIC SPURS) RENO/THOMSON (MODIFIED AREA WGHHNG) CZ3© WANG (WC W5) 30 SPREADING flLTER Figure 5. Iso-vorticity contours as represented by the various filters. Large tick marks represent a half-grid length. 118 .825 got? 2: Ho 5535: “organ .0 Bani 523mm 03:8 oz; 82:59 6.2 858:. 233536: 2.336 63.52: fizxoo: u} a} n) 3 0.0 6.? Gd 0.“ O. a 0.0 3 3 3 O.” a." O. - Gd 0.. 0d 6.9 3 O.“ O. — Gd LyllulllrtrlLl I? L n L P Lullbl If"? TLIIIPILIIIrII—ILI‘rIlPLILI..-NdI ILIILI lb: I—u LIOIPVL... In L L h Ndl 10.6 r! .0d 10.6 Y V -«d «.6 find 7 + rod W, I. W. :3 w, ( ( ( 10.0 #96 T06 r r v 36 .3 r3 Fltltll t ll lilnl -llllau liltl 3 IO.— ItlIDIII-lsl|tolll;l .Inlulloillnl- [04 (ill! I a-ll .- . I? i-| -5!!! IO.— Aflzaam 35938 5:8 82:39. 93 Exa>mmzém§6 .5 S 5 ad 3 0.0 ad ad 0.. ad ad ad 3 3 o.“ o. — od 3 3 3 a.” a c. — 3 LIrLLPiPLILILIILILILIIFN-onl b b hLDLIb L I-LLbLl Nu? LLILL blbfLFbihLl "3°!- 1 6.0 t 0.0 r 0.6 1N6 12 rad rvd W tvd 1W.” z. W mod .od 1.6 . W.56 and 0.9 119 .328 Set? 05 «o .5838“ .330”— 382: .5 uuzwfi 522m 05:8 oz; .. oi 9n 6d 0.. 6.0 L p TIL In L» b g t36 rllalllll‘lllln'l’llllll l'lsllllrlou. 522% 95.953 558 c 6.6 6." 6..“ 0.— 3 LI |—[ r? b! L L b (w («M 8256.“; 65 8583 28332: 64 6.6 6 6d [I llhlIIlrl. Ir.I|LII-L1uIILIllJ 82:28; 65. Evanxfimzsgu 1.? rod :2 -od 1.6 r60 0 Ea.— -3. :2 7... rod r66 0 7| (up as» 6.0 0.0 T I .lbll. LIIIIPIEL- 253.3 83228 65.8: e 6.” 6d 6.- 6.6 .l cLll'LfaulLll 0. I1 filfiuoll as» Aggy 3:95:56 6.” 6d 6.— 6.6 115-Iii»?! ilrtllLIlttrlt LII-Pliltfldl °a° (by tvd 6.6 v|!| Ir 1|: ID..!D..1'!- I.'I.Vt.l.|alv! 0.! . rue. 120 . 6.50285 m: 00 832: 2: new 8:: 3538 5033 88:58 Bum .w 2:3... 522%. 9.53 92; e 0.? 0.n 0n“ 0.. 0.0 I- .LI lLI Ir IPI IILII. 1.. .IL.. iwl0dfll Amuznam 25.933 5:8 9 .v 66 6d 6.. jiblillhl [P lb! [F I. i» o 6 6 q 0 .. 8 8 I I I (8)4083 (20m 82580; 55 BESS 2829505: e 0.9 0." 0." 0.. 0.0 #1 - IL .- Ir -Ilr-.I.I.—I I. ,IPI-I.L|IIIILII-&|T°.°NI v v mod? v V I wad? V v v wed-I v V v 0 won! 7 n 0 v v v .166 V v V wad. V v v wad. JV I V "0.0“ v V Fl--. . .III . I..!II.. ....III II III! IIIIIII, [cog 6 6.0 0.» 6d 6.. 6.6 II LI lb! b [bl L LII'IIPIIIII TO“. M .v6.6NI mod? m "0.0.! moan v . v o 100 v v (20m (shuts 25838 5.6225 55.8: 06 0.» II . LI .. IIIFIIIILII I. .u 0.." 0.. 0.0 L! IL.-TI.- .IIrI L100“! : W I H II 1.....- m Mod? m3... No.2: Ma..- (2)-M3 121 casein 63358255 35 d oSmE A305 A8323...“ ocoala 122 mg... 10. Quadratic spline vorticity distribution of three typical vortices on their nearby grid points. 123 1"‘11211—1‘1 hdilh t —l Hat-1 , Figure 11. Zone of influence of a vortex sheet. 124 ..s [I \II 0.8, \ I \\.6 . w . . . \NKstKNKN 7 . / .lfl. Eaten .825 .3 2:3”. / mm m I. Nut... J. RN N \NKF\K\NKN% x0 7 _ _ _ 1 . . “ xxxkx «\N\\\Ku\N~huh _ . n _ _ Me a 125 .332; no. mo 829$ a can it? Ewan a mafia :BEoc 825m 2: E 83? 8.823883 .8589: 95 830 5033 echo 03:23. .2 oSwE 6.. v u. v 0.9 8.. u 5 83:? no. go 803% a mafia 3.8 a»... 5:3 295 a mafia 83 amok 4.5 x .53 Ha Qmo $8 x .E3mum 9% 66336. 6936 6966 6VX6¢ 6NX6N 663.66. 6986 6986 6¢X6¢ 6NX6N . . . _ . . . . 11 1»l|L\plIIIAt0.O /|\. a ,2 m 3 a a O .\.0 I0.N \.I z (\ 1 11111.15. 611113.1111II, -2. .66 16.6. I 6.6N .66... I 6.6.. r3... (:4) 30333 mums [RELATIVE ERRORI (7a) 126 40.0 q r‘-o.o 30.0- 20.0- 10.0- r‘-o.a r‘-o.7 o.o- , . I . 0.0 4.0 8.0 12.0 16.0 At x 1.0E—03 sec. Figure 14. Relative error between the initial and final positions of a point vortex after one complete revolutl' 1011. 20.0 127 40.0 30.0 - Vs (m/S) 300 mg. . b 500 me. GOO o 60 Ct >00 0 oooooo oooooooo 00° 00 -o.sa -1 .C - - - -1.0 -O.5 0.0 0.5 1.0 Figure 16. Geometric distribution of the system of 105 discrete vortices (10:100 ms.) 128 30.0 — Vg—component . --- V,-component 25.01 " ‘1 \ \s f I V V \ \ 1‘ \ ' ' A I J g? 20-0- P ‘ “i ' 1"(i J» t I‘\ ’\ ‘T.’ ’v ‘\ " Ii 0 ‘ l L. l t. 3 1 5.0 " I .2. l -0-’ .I I 2 a) Q: 1 0.0 -4 I 5.0- 4 0.0 ' T ' I ' I T I ' I ' l ' 1 ' I ' l ' T ' 0.0 0.1 0.2 0.3 0.4- 0.5 0.6 0.7 0.8 0.9 1.0 1.1 Time (sec) Figure 17. Relative errors in V1. and V9 component velocities as a function of time. EXACT 129 NUMERICAL .0 -11 1.0 fi1!1 .1 T Y rfir r—1 1 r r V r r 0.8 i 0.5 r 0.5 0.0 r rft T r T r Y rf‘rrr'r r 2 CUFF .1.‘ - a; - at; r T I V f Y I I 0 1 1 1 +14 1414. 1.1—141414 4H,.414I4 44143414141MHJH 0T4 14141113111414 14L 0 u o. ”w .. .. u u .4. m .. u u a m 1.0 1 1.0 1 I 1 1 1 1 1 ‘00 1 1 1 1 | I 1 I 1 I 1.0 frj'f" I" 0.5 0.5 frrfTrrfir ' I 0.0 0.0 0.0 lfirrrfllrrrrl} rYY—I'Irrr rrrv l .0 41.5 'vrl'rrrrrfij r I 1 r #14 .4 J J 4 1J.J.J-4.J-I4.-4t41414.414.JL .111414141414..414.JIJ14.41414 -4L 114.14.4J.J14141414.4|41414IJ|41411VM c m n... M w M..- u M m .4... u u «m. Figure 18. Trajectories of systems of discrete vortices EXACT 130 NUMERKIAL ‘9 3 J8 .1.. r.) u/ z\ \\ \m\\ J/ .t } / ./ (5 3 1.0"" .IJJJJJ ilgLvl‘ YJI].J.414 11.41114141411414 gr. . _ ._ u w 5 e _ .. 4 I. T4114 J 4-1.4.41- 1.114113 a u u M . n 1 0.: on -o.s I: (con’td) Figure 18. 131 NUMERICAL EXACT 1»\ '1 1 1/ 0.54 / 0.0-1 \( Figure 18. (con’td) 132 f r :5 I5 n— n. n- O T .I. a m a t . I 1 . I I M 5 m , a .lo [a m . m . u . I r us 3 t J w.wiwio mime. m mimic Wm» . _ T35 33! .3593 333 #8.? 3259.3 O 0. 1| ll 15 rd 0 O 1 In R a D V to t A m A I .l m 4 m 4 0 CI 0.. nu. t n A! d 1 I Wlfij ‘-4 1 ‘ 1 I ‘ d 1 a a a {34. T35 goo—u? 33:55:. we: Tax-.5 330.2. 35:055.. x/R QR 1 £31“!!! was.» d #4 93.5 58.2.. .3... '93... Figure 19. Tangential velocity profiles at discrete time steps representative of short and long time errors. 133 0. 5 In. rm m .. IO M . - v t o. fiaaua id... 1- .. .. ,... w w m. 385 98.2, 363.. n _ :5 n a a rd. 0 f I. m a. a. I v Q. n m w m c m 4 4. 3:5 332.5 .23”. x/R ~--” 1.0 In. n. e A IO . d T q I u q 1 «.40.: 325 £83. 3.3m ‘ fl 1 f T 0.0 0.3 -O.5 d d d .1 - awn -1-... 3E. 33.2 3..qu ‘ x/R x/ R In: . a . . v . [0 . a m 5 m .. 6. . . I. a a ,. . b.- o 1 1 a a q Am ‘ n-u ‘ a 1 4.. w w m a. a m ?\E. 33.... .33. 0. . to. a 5. 1.0. s I .I u. .t m .4 = f o.- o 1 d 1 d 4 m 1 1 11 44 1.1. a w a a. a a 5.5. are: .33: Figure 20. Radial velocity profiles at discrete time steps representative of short and long time errors. 134 \ 6x6 (L = 32) \J/ // \\ / \ / \ f ‘ 8 x 8 \ l (L = 52) X / \\\ /// \.___.m 10 x 10 (L = 80) Figure 21. Uniform grid squares superimposed on circular domain. 1. (ms) 1‘ (ms) t. (ms) 135 500.0 400.0“ M U1n-1 12.3 cm/s +4- 03,-12557 em/n ‘ o-o UMIISZJ urn/a ‘1'30" H 0210.0 cm I 0.0 2C .o ' 06.0 fdo # 36.0 V 06.0 ' 76.0fi 06.0 Number of cells. L V 90.0 500.0 400.0% 300.0-1 H UM-‘HZJ cm/s +4- U-m-uSJ em/o ‘ o—o u,,,-132.e an/o 0,-600 lb! - 10.0 an o .o - 36.0 i r fl 1 r v V 1 V r 40.0 30.0 60.0 70.0 80.0 Number of calls. L 20.0 ' ' I ' r ' I ' I V 40.0 50.0 00.0 70.0 Number of cells. L I 30.0 , 1 00.0 90.0 Figure 22. Complete mixing time as a function of the number of cells. (Swirl strength 89:12.5 cm-m/s. Number of vortices per cell Nc=2) 1. (ms) 1‘ (ms) t. (m a) 136 400.0 100.0“ “0 HOFOO" L-32 .G-OOOW .0 7 100.0 100.0 200.0 Inlet jet velocity. Uin (cm/s) 250.0 400.0 I 300.0-1 200.0- 100.0“ M 0 Jo' 4—1- 01:80° L-52 .0 ' 100.0 f 100.0 2000 Inlet Jet velocity. ”in (cm/s) 25 3.0 500.0 400.0 "* L-00 .0 ' 100.0 v 150.0 200.0 Inlet jet velocity. ”in (cm/s) f 250.0 Figure 23. Complete mixing time as a function of the inlet jet velocity. (Swirl strength 89:12.5 cm-m/s. Number of vortices per cell Nc=2) 137 500.0 W 0":30‘3 +4- 9j=60° ‘ (3'0 91=90° 400.01 0'0 ' l ' I F l ' I ' I ‘ T ' 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Number of cells, L Figure 24. Average complete mixing time asa function of the number of cells. (Swirl strength 89:12.5 cmom/s. Number of vortices per cell Nc=2) 138 300.0 8-8 U1n=112.3 cm/s -l-—-l- Uin=125.7 cm/s e—e Uin=152.8 cm/s 31"600 . BK-éfi MIL-210.0 cm/s 200.04 A U] E v *-H 100.0“ 000 I I I 0 1 2 3 Number of vortices per cell, NC Figure 25. Complete mixing time as a function of the number of vortices per cell (Swirl strength Se=12.5 cm-m/s. Number of cells L=32). 1* (ms) 139 300.0 M N¢=1 .1—1- Nc=2 BJ—BOO e—e NL=3 200.0-1 .1 100.0- 0.0 w I ' I ‘ I 50.0 100.0 150.0 200.0 Inlet jet velocity, Uin (cm/s) Figure 26. Complete mixing time as a function of the inlet jet velocity (Swirl strength 80:12.5 cm-m/s. Number of cells L=32). 250.0 140 400.0 B—B U1n=112..‘5 cm/s +-l- U3n=125.7 cm/s J G-O U1n=152.6 0111/: 316-36 Inn-210.0 cm/e 300.0- 0.-80° J t* (ms) N O P O l 0.0 ' I I 0.0 10.0 20.0 Swirl strength 56 (cm m/s) Figure 27. Complete mixing time as a function of the swirl strength. (Number of cells L=32. Number of vortices per cell Nc=2) 50 .0 t* (ms) 141 300.0 8-8 59: 4.1 cm. m/s 4—l— Sg= 7.9 cm. m/s 91-600 G—O Sg=12.5 cm. m/s - H 89-175 cm. m/a H Sa‘27.8 cm. m[s 200.0- 100.0d 090 r I r I ' I 50.0 100.0 150.0 200.0 Inlet jet velocity. Uin (cm/S) Figure 28. Complete mixing time as a function of the inlet jet velocity. (Number of cells L=32. Number of vortices per cell Nc=2) 250.0 142 \\\\\O\“‘Ol9l"”,,” \\\\\‘9\‘OIOII'OI”,,/ \ \ K x \ t i I r I I I I \\\\‘0\OloI///// \ \\\\cI//I \ I”-\\\\ ’(ll.\‘\\ ’l”l.‘\\\\ \“~“fl”’ f’po-o—Q-sKN \\~“a” {(l'nau-fiu\ \ s s s . . . z I I I I I l l t l l s \ \ \ \ / ’lp/lllclclo“\\\\\ \ \ / l / I I I i \ \ \ \ \ \ x \ ’ll‘"l“\“ VR ”-.--‘§~\' ”"p*—~~~\K\\ ‘u‘In‘l-l'", \.\\\.\9\||9|ol””cI/// \ \\\\\.\\|‘o|olol”cI//// I \ \\ \\o\o\1tO\olnlaI”///// \\\\.II.IIII/ \\\ ‘t"”’ s\\‘.:/// \ x x x x s e \\\\”” I z z z x x M . s \ x x x \ x \ \ x k \ _ . . u I I I / I / Ill/IIIIIv—~~~\\\\\ / 2 z x I a M ~ e x x x x x \ l!!//’oltolo|.\\\\\\\ Ill/I/Ill'i\\\\\\\ ’ r I I I I x z / I I 1 I I I l l 1 l l i \ \ \ \ \ I I ’lllolololv“\\ “““‘\ Q0 x/R .4fi "-.--~‘\‘ 'l’p'p--~.~~\\\ I/(//""~uh\\\\\\\ \\\\\\\‘li§§1§l//// ““""’ \\\\|\‘ |""”/ \ \\ \\ \\ \i \i II II I! I! II I \“II|'"’/ ciao. \\ \ I l o \\ I\ o o a o a .o.. u u . -\\ . . . . \ {Ivhr I y ‘ - .0 ‘\‘§..— ’\\“\ \‘~\~-’a” ’10- \ \ \ \ \ IIII“IQ\I\\\\\ \\\\ ~~~-.-.M’/// I I ’v““- I [flyo—o-o—u“\ \ III“"*" I I l I I i i \ \ \ \ \ \ x \ III/lll“‘\\\\\\\\ s} non. discrete vortices duced swirling mo \ 0.0+ V 0.0-1 2: IQ . s > 4-5‘ -z.s~ J -‘o: -5.3 . Y -1.0 -I.0 -0.5 0.0 0.5 1.0 x/R t = 0 ms. 1.: $0.: 1 0.54 23.0-4 '3‘ 9 I m > \ 0.0+ V 0.0-4 .E‘ o E a > -o.s- 45.04 -‘-= '59.: f -1.0 4.0 -0.5 0.0 0.5 1.0 x/R 1.- teen: 0-3" saw “. 8 O \ a: C \ Uri m >. 3‘ :5 E > 4.4 ammo-4 -‘.‘: -zem.: r r r -1 4.0 48 01 0.3 1.0 I/R Figure 30. Fuel particle trajectory and vorticity profile at y=0 as a function of time. Injector at 60 degrees (cormterclockwise from haimntal). Uin’ 210.0 m/s. v/R 145 1311.0 -mu Vodidly (1 [609) —w Vorticity (t/uc.) Vodiaty (t/uc.) Figure 30 (con’ul). ‘ ‘4_l oi: olo x/R -l .0 0.34 -0.5'1 t=0rns. ’0 I“ . 4 fl 0 3.5 :83, oi: 43.: 0.0 x/R -‘I° “ I t=10ms. (without swirl). 210 0 m/s. time Uin= Sqeamflneeandvebcitypmfilesatysoasafimctionof Injector at 60 degrees (counterclockwise from haizontal). Figure 31. t=29tns. t a 271 ms. Figure 31 (con’td). 148 -1'.o 43.: 0.0 of: 1.0 x/R t = 0 ms. 0.5% v/R 0.0-l .osJ -“ -l.0 Vorticity (l/sec ) thuculy (l/Icc.) Vorticity (l/soc ) t=20 ms. Injecurat 60mm“ Airand Fuelparticlenajecmandvcrticityprofileatyw counterclockwise x000..." scam and \//\ 400.04 \ HON: . - I .0 -0.5 0.0 0.5 1.0 x/R VOWL‘T— 543.04 om -SN.O‘ "M:% Y .I .0 -o.3 '10 0.5 1.0 r/R 500:: * ,/\~ ., \ V, ‘l l -w.0-4 \ \ \ - noon-4 4500:: , . -v .0 «3.3 0.0 0.5 In x/R as a function of time. from horimntal). Uin = 210.0 m/s. 149 1.0-1 500.3 0.3‘ ’1‘ 0.04 o 3 E C v § .4 Q *3 i ; q 4.3-1 -$.: ~10“: 1 r . -1.0 -1.0 -o.: 0.0 0.: 1.0 x/R t a 30 ms. 500:: "5 0.0-1 0 3 \ 5 .5." .2 E W4 -‘m 1 fi' .' -|.0 “10 0.0 0.5 1.0 x/R t a 32 ms. ‘- soon: [22:] - .. -. . e . . 0“ 0.0 ‘ . ‘ a 0.: :“o: "' “6' 8' J” on .5 53' .- -~. 0 o ’ g ' ' o A no 0 . O . . .': '0 , . .1 I ' a . ' . * 8 . g 0 o to flu . 5.. < K 0 o. . . o O o‘ .. ' 0.. a \> :0. .3. :3 0., . 3:.:'.c: ; ' o .0 ".. '0 0 s 000 3 O. 0 .0. . ~ . . E J . . ‘ o .0 v. .3 .. O". > 0.0-1 0 0. .. . 'O'.$ . . 4 ~ '- ."-°1\.°. .‘-.:, o4 . o v . . .. 0‘ . . O . 0.. l e 5’ . ‘0 .. .~ 0 .O‘ o o. - d 1m 43.3 0.0 0.: no "'° '°" M u 1.0 I]! V” t=198 ms. Figure 32 (con’td). 150 Velocity (tn/0) 8 45.04 -JOO: T I V. -1 .0 -01 0.0 0.8 i .0 "/R t=0ms. Velodly (m/c) «.0 «5.: do 0.: 1.0 x/R t=lOrns. Velocity (tn/o) t=20ms. Figure 33. Streamlines and velocity profiles at y=0 as a function of time (with swirl). Injector at 60 degrees (counterclockwise from horizontal). Uin= 210.0 m/s. 151 - «QC-Ill ~00 1T \ E v g 0 > Velocity (m/e) -l.0 -65 0.0 0;: 1.0 VoIocRy (m/e) t=l98ms. Figure 33 (con’td). 152 -‘ C: \/ -|.0 é. $3 0.5 1.0 0.5-1 A 0.0-{ ml 4M -1 d u n \ 3 .4“ .2 E > t=0ms. Vorticity (l/aec.) Vorticity (I /:cc ) 1.3 0.3% 0.0-1 4 4.54 q -‘ u: r r r - 1.0 -0.0 0.0 0.3 1.0 r/R 1m; son-l 0.0-1 -50.0'l um; . , r -l.0 -0.0 0.0 0.0 1.0 x/R ”5 250.0% 0.01 .mm 410.: fi—i v ' -0 ‘03 0.0 0.0 l 0 x/R Figure 34. Fuel particle trajectory and vorticity profile at y=Q as a function of time. Injector at 90 degrees (counterclockwise from horizontal). Um: 210.0 m/s. 153 olo r/R 5“ 25.04 0 30.x 3 b_o..to> 0‘ . u....¢.: . ~65 015 W -l.O 0.54 a} -1.‘.‘ x/R t=30ms -0.0 0.01 “mm duo-x .v bitg «coon-4 m l/R x/R t=53ms. 1:00.01 a. m M. a A.uo-\ .v 3535) -1 $0.04 r/R r/R taflfinm 1mmm34amwmx L __J r . v Velocity (tn/s) Velocity (tn/e) Velocity (tn/e) t=20tn8. Figure 35. Streamlines and veloci profiles at y=0 a a function of time (without swirl). Injector at 90 degrees counterclockwbe from haizontal). Uin" 210.0 m/s. Velocity (m/e) Velocity (tn/e) Velocity (In/e) t-256ms. Figure 35 (MM). 156 tom: Vorticity (I/uc.) -‘m‘ r v v -10 -0.5 0-0 0.5 10 Vorticity (t/uc.) m4 1/8 Vorticity (l/eec.) l/R Figure 36. AirandFuelpmicleuajectoryandvaficitypmfileatyaOasafimctionoffime. Injecttr at 90 degrees (counterclockwise from horizontal). Uin’ 210.0 m/s. 0.5 0.0 v/R ~13 . -10 0.5 v/R 2 -1.: . 157 Vorticity (I /eoc) Vorticity (I /uc.) Vorticity (i/eoc.) 1M3 is 010 0.0-l «mom 4.5 010 Figure 36 (con’td). 43.3 Velocity (m/e) 1 t=0ms. t-lOms. Figure 37. Streemlinee and velocity profiles at y=0 as a function of time (with swirl). Injector at 90 degrees (counterclockwise from haizonml). Um: 210.0 m/s. 1? \ E v 3‘ § > Veloeky (In/e) Velocity (tn/e) t-257ms. Figure 37 (con’td). 0.54 00-4 v/8 fl. y/R Figure38. Vorticity (l/eoc ) t=0ms. Vertlclty (t leee.) t=lOms. thlculy (l/eoc) 104 0.0‘ .o -o.e 010 0.: 1.9 x /R 23.04 0.0-1 - e «is ole 03 1.0 If" ace: m val-4 tonne-4 new w «new J -: .e -o.s ale 0.: m I/R Fuel particle trajectory and vorticity profile at y=0 as a function of time. Injector at 120 degrees (countacloclrwise from horizontal).U m: 210.0 mls 161 m 1m ‘ A.uvu\ .v 3.093) 0.04 flu: 0.0 0.0 '03 x/n x/I t=30tns. l.0 ob m4- Ade.\.v .93...) w 004 4w: l.0 r/R -6.e an: t-36ms 4 w J1 m - ..uee\.v been.) -l.0 x/R tszfium. figmnw(umm& Velocity (m/s) 4.54 1: E V 3. m 010 ole .o x/I ‘ 1 W A i ‘ /"‘\ é ./ 3 '/ 4 V/ -W .9; .0 _a‘. 0.. 0:0 “ t-20ms. Figure 39. Streamlines and velocity profiles at y=0 as a function of time (without swirl). Injector at 120 deuce: (counterclockwise from horizontal). um: 210.0 m/s. 163 Figure 39 (con’td). 164 10m; scum 0 $0.04 0 e: eee eeeeeeee 091 . Vorthlty (l /sec.) seem Vorticity (l/eoc.) 4”: v/R Vorticity (l/eec.) x/R Figure 40. Air and Fuel Barrie le trajectory and vaticity profile at y=0 as a function of time. lnjccttr at l deucee (counterclockwise from horizontal). Uin’ 210.0 m/s. 165 1m; ”.04 00‘ Vorfiaty (I lac.) 40mm d”; r/R Vorticity (t/eec.) l an Vorticity (I [sec] «mead - «can: . , . m i " .0 ‘05 0.0 0.5 Lo ~05 ‘/R “an 40 (con’td). 166 Velocity (m/e) Velocity (m/e) t-20rns. Velody (tn/e) ”its: tan-4 0.: 45.04 J - 4.0 43.6 on 0;: 1.0 n/R olu'l -1.o «is 03 ole 1.0 Figure 41. Streamlines and velocity profiles at yaO as a function of tune (with swirl). Injector at 120 degrees (counterclockwise from horizontal). ”in“ 210.0 m/s. 167 Velocity (tn/e) 0.5 1.0 1W 3‘ . \ 5 g 0.0- > 45.01 -l .e 4'» is of u Vcbdty (tn/a) t=332rns. Figure 41 (con'td). APPENDICES APPENDIX A APPENDIX A NUMERICAL FILTERS AI. Derivation of Christiansen’s filter as an illustrative example YA r l r | I l —fl-——+-——P- I I l = 1 : #1: I I I d--| .———I—— I l J " l : i+l,j ' I I l -—|—-—+-—-|—-d I l I I 3 1 l >X Figure 42. Vorticity distribution in a grid, with shading showing assignment of vorticity to grid point for VIC method is assigned to grid point (i,j) is assigned to grid point (i+l.j) is assigned to grid point (i+1J+l) W§Q§ is assigned to grid point (i,j+1) 168 Notes : (1) (2) (3) (4) (5) (6) 169 Each point vortex is assumed to possess uniform vorticity within a square cell with dimensions (d x d). Unit circulation (yp=l.0) is assumed for each vortex. The point vortex coordinates (xp,yp) lie at the center of the square constant-vorticity core and must fall between four grid points. The origin for (xpyp) is the lower left hand grid point location. The filter which credits the circulation to the surrounding four grid points is obtained through a bilinear interpolation scheme (area weighting). AL)" (719i J, are the area and vortex circulation assigned to the grid point (i,j) respectively; §(x,y) is the vorticity and (xg,yg) are the grid point coordinates. §A=1.0 f §(x.y)dA Wit = f anymu = EAIJ Dividing equation (A.2) by (A.1) : We (7pm i 'Yp A Aid. (A.1) (A2) (A3) (A.4) :45? Ai+1,j 170 [(xg+dl2)-(xp-d/2)l-[(yg+d/2Hyp-d/2)] [xg+d12-xp+d/2]°[yg+d/2-yp+d/2] [xs'xpfifl'ws'ypm] (A.5) [xg'xpw] . U8- Yp+dl (1 ° (1 [(xg-xp/d+1]-[(yg-yp)/d+1] [1-(xp-xg)/d]-[l-(yp-ygvd] [l-Tl(X)]°[1-11(y)] (A6) [(Xp+dl2)-(Xg+d/2)]°[(Yg‘I-d/Z)-(Yp-df2)] [xp+d/2-xg-d/2]°[yg+d/2-yp+d/2] pr-ng-Iyg-ypwl (A7) [xp-xgl-[yg-YPM] d- d [(xp-xgydlttyg-ypxdnl [(xp-xgydI-Il-(yp-ygydl [Tl(X)]°[1-Tl(y)] (A8) 171 Am” = [(xp+d/2>-I-[C I I L__@.f3__I l J #5 I I I L.______I I l I I I ~ #7 I I---Q---I Figure 43. Method of Images 175 Notes : (i) Sign Convention : Counter clockwise is positive, clockwise is negative. (ii) to = §+in ; 0) = C-i'n (B-5) «)0 = co+ino . «00 = co-ino (13.6) Consider the ‘n-direction, #1 #2 #3 W((D) = 'iK{1n[m-(Co+iflo)I-1n[(0-I§0+i(2K-flo)II-IIIIOJ-(CO-iflo)]+ #4 #5 1n[a)-{ §0+i(2K+110) ] ]+ln[to-{I;O-i(2K—n0) H- #6 #7 ln[o)-{§0+i(4K-no)}]—1n[co-{§O-i(2K-I-n0)}]+... } = -iIc{ln(u)—coo)-ln(to—EO-iZIQ-In(a)-E'O)+ln(co-mo-i2K)+ 1n(o)—coo+i2K)-ln(0)-coo-i4K)-ln(0)-'cb‘0-i2K)+...} = -iI<{[ln((o—mo)+ln(0)—mo-i2K)+ln(tn—coo+i2K)+...]- [ln((I)—(To)-ln(co-50-i2K)-ln(0)-50+i2K)+...] } nz-oo = -iK {Eln[m—mo+i2nK]-2In[w—6fo+i2nK]} n=-oo = f(w) + f(_m) (13.7) 176 where f(to) = 41:23, [Manx] (B.8) [1.1-co f(co) = ixifh [to—50+i2nK] (13.9) n=-oo Now, consider f((o) : f(m) = - iKln{(to-coo)[(to—mo)+i2K)][(m—wo)-i2K][(m—mo)+i4K]+...} = - iK1n[(m—coo)fi£(0)-too)2+(2nK)2]] n = = “II "—‘§’§°’- ”n‘ nll+r(2£—n'£301121(2nK)2] = -ix{ln[m: HIM-m nK)}2]] + m[2,th:ll(2 nK)2]I I = -ix[1n[3inh{“(—2HK01}] + 1n[2K “321111)le = -ixln[smh{"—(“;'—I;°i1] + Cl (13.10) where Shun—(92%} = [fi—ZmI—(ioflfigl+r—(%mfl}2]] (3.11) (13.12) D II 5' l—'—I ”I; 5 .__Z 177 Similarly, it can be shown that : E25) = ixm[smh{’i(—°;-_me)—I] — C1 (B.13) Substituting (3.10) and (B.13) in (B.7) : 2K )II} “((03) = -iK{1n[Sinh{n—(%_iw-Q)-}] — 1n[3inh{7£(_‘L“—b_ ' (“(20.1%) - -ncln 1t( —) Sinh{ 2K } . SinhICm+inml 1 = -1K1n —,—-—-,—— Smh[§m+111p] ) , Sinhtcm)Cos( nm)+iCosh+Sin2< _____n_m)] 2 Sinh2(§m)+Sin2( 11p) 12 1n[Sinh2(§m)+Sin2( no) 2 Sinh2(§m)+Sin2( nm) (B.17) APPENDIX C CI. where APPENDIX C COORDINATE TRANSFORMATION Square (Cm) -) Circle (x,y) We begin with the transformation (5.1) (See Bowman [11]): z = sn(co)dn(to) cn(0)) z =x+iy a) = §+in sn(0)) = [(sd1+icdslc1)/(1-d2512)] dn(o)) = [(dcldl-ikzscs1)/(1—dzslz)] cn((o) = [(ccl-isdsldl)/(l-dzslz)] c = cn(c.k) cl = cn(n.k’) d = dn(c.k) d1 = dn(n,k’) s = sn(C.k) 81 = sn(11.k’) k’ = m 179 (CI) (C2) (C3) (C4) (C5) (C6) (C7) (C.8) (C9) (C10) (CH) (C. 12) (C. 13) 180 Substituting equations (C2), (C4), (C5) and (CG) in (CI) we get +. _ [(sd1+icdslc1)/(l-d2812)] [(dcldl-ikzscsl)/(142312)] (C14) x 1" " [(ccl-isdsldlyu-dzslzn ' (sdl-I-icdslclxdcldl-ikzscsl)(cc1+isdsldl) (ccl-isdsldlxl-dzslzch1+isdsld1) [sdlzdcl-ikzszdlcsl+icdzslc12dl+k2c2dslzcls)(ccl+isdsldl) [clds(d12+k2c2512)+icdlsl(clzdz-k232)](ccl+isdsldl) (02c12+d2d1252312X1-dzslz) _ (cclzds+ic1d2d1szsl)(d12+k2czslz)+(ic2cldlsl-cdlzdSSIZXClzdz-kzsz) _ (czc12+d2dlzszslz)(l-d2312) ic ldl sl [d252(d12+k2czs 12)+c2(c 12d2,k252)] = cds[c12d12(1-d2812)+k2312(c2c124-d1252)] + icldls,[d2(d12s2+c2c12)+k2c2s2(d2s1-1)] (c2c12+d2d12s2512)(l-dzslz) _ [cds[c12d12(1-d2s12)+k2512(c2012+d1232)] ] (c2c,2+d2d12s2s,2)(1-d2s12) . I:cldlsl[d2(d1252+c2c12)+k2c252(d251-1)] ] l (C2C12+d2d1252512)(1412512) 18 1 Thus, x = = f(Canvk) (c2c12+d2d1252812X1—dzslz) c d s [d2(d 282+c2c 2)+k2c252(dzs -1)] =[111 izziz 22l :I=g(c’“’k) (czcl +d (11 s 51 )(l-d 51 ) II. Circle (x,y) --> Square (Cm) X-Coordinate Squaring both sides of equation (C. 15) we get 2 [cds[c12d12(1-d2312)+k2312(c2c12+d1252)] T x = (c2c12w2d1232312)(1-d2812) (C. 15) (C. 16) (C. 17) Now, c2, d2 can be expressed in terms of 32; and of, d12 in terms of 512 (see [11]) as c2 =l-s2 d2=1-k232 012:1 ' $12 (112: 1 ‘ k’zslz (C. 18) (C. 19) (C.20) (C.21) 182 Substituting equations (C.18)—-9 (C21) in equation (C.22) gives (1-s2)(1-k2s2)s2[(1-s12)(1- 2312){1-(l-k232)312]+ x2 = 1:2.o.12{(1-s.2)(1-s12)+(1.1831232I]2 (C22) [ (1—s2)(1-512)+(l-k252)(1-k’2512)s2312]2[1-(1-k252)512]2 an252+an4s4+an656+an858+an10510 (C .23) P0+P282+P4s4+P686+P888+P1Os10+P12812 where an2’ an4, an6' an8' an10' P0, P2, P4, P6, P8, P10, P12 are polynomial functions of 51: Pm = [1-2(3-2k2)s12+(15-18k2+4k4)s144(5-3k2)(1-k2)s16+ (c 24) (15-13k2)(1-k2)s13-6(1-k2)2s110+(1-k2)2s112] ' Pm = -[(1+k2)-2{3+2(1-k2)k2}812+(15+7k2-22k4+4k6)sl4- 4(5+7k2-5k4)(1-k2)s15+(15+22k2-29k4)(1-k2)313- (c.25) 2(3+8k2)2(1-k2)25110+(1+3k2)2(1-k2)2s112] ans = [l-2(4—k2)s12+(5+k2)(5-4k2)s14-4(5-k2)(2+k2)(l-k2)s16+ (7+5k2)(5-4k2)(l-k2)s18-2(8+7k2)(1-k2)s110+ (C.26) 3(1+k2)2(1-k2)23112]k2 ans = -[-2812+(11-7k2)sl4-4(6-k2)(l-k2)315+2(13-7k2-2k4)(1-k2)318 ‘2(7+2k2)( 1-k2)281 10+(3+k2)(1-k2)281 12] [(4 (C.27) an10= [$14-4(1-k2)s15+2(3-2k2)(1-k2)s18-4(1-k2)s110+(1-k2)23112]k6 (C.28) 183 [1-4812+6Sl4'4516+818] -2[l-(5+k2)slz+2(5+k2)sl4-10s15+(5-2kz)s18-(1-k2)23110] 1-6(1+k2)812+(15+22kz-k4)s14-4(5+7k2-2k4)s15+ 3{5+4(1-k2)k2)}$13-2(3+2kz)(1-k2)sl10+(1,k2)2sl12 -4[-s12+(5+k2)sl4-(10+k2-k4)sl6+2(10-3k2-2k4)318- 5(1-k2)2s11°+(1-k2)23112]k2 [5314-4(6-k2)s16+(36-20k2-k4)s13-4(6.1;2)(”(2)8110+ 6(l-k2)23112]k4 -2[-2315+3(2-k2)s18-(6-k2)(1-k2)sl1°+2(1-k2)251121165 [SIS-2(1-k2)sl10+(1-k2)5112]k3 Equation (C.23) can be written as where 2 4 6 8 10 12 PXO+PX28 +Px4s +PX68 +Px88 +leos +PX128 = O P x0 = ~P0x2 sz - P -P x2 " nx2 2 P14 = an4'P4x2 (C29) (C30) (C31) (C32) (C33) (C34) (C35) (C36) (C37) (C38) (C39) 184 Px6 = an5‘P5X2 (C40) ng = ngxz-ng2 (C -41) leo = anlO’Plox2 (C42) P1112 = -P12x2 (c.43) Y-Coordinate Squaring both sides of equation (C.16) we get d d2 d 2s2+c2c 2 +k2c2 2 d2 -1 y2::[31 151I ( 1 1 ) S ( 81 I] T (CA4) (c2c12+d2d1252312X1- 2312) Substituting equations (C.18)—> (C.21) in equation (C.44) gives (l-slz)(l-k’2812)312[(1-k252){(1-k’2512)52+(1-52)(l-512)}+ y2 = k2(1-s2)s2{(1-k2s2)s12-1)12 (c.45) [ (l-sz)(1-slz)+(l-k2s2)(1-k'2312)52312]2[1-(l-k252)512]2 2 4 6 8 10 12 = PnyO+Pny28 +Pny4s +Pny68 +Pny88 +Pnylos +Pny128 (C 46) P0+P282+P4s4+P636+P858+Plos104-P128” ' Where Pnyo, Pnyz, Pny4’ PHYO’ Pays, PHYIO’ Pny12, P0, P2, P4, P6, P8, P10, P12 are p01yn0mial functions of 31 viz: ny6 = ny8 = Pny10= Pny12 = 185 512-(4—k2)sl4+3(2-k2)s16-(4—3k2)s18+1-k2)sl10 —2[2s12-(9-2k2)sl4+(15-7k2)315-(11-8k2)518+3(1-k2)sllo]k2 [2(1+2k2)812-2(4+1 1k2-2k4)sl4+(12+43k2-20k4)516_ (8+36k2-29k4)s13+(2+13k2)(1-k2)s 110119 -4[812-(5+k2)s14+(3-k2)(3+2kz)s15-(7+3k2-5k4)s18+ (l-k2)(2+3k2)s11°]k4 [$12-(4+7k2)s14+(6+23k2-4k4)s16-(4+25k2-10k4-4k6)318+ (1+10k2+4k4)(1-k2)sllo]k4 -2[-sl4+(3+k2)sl6-{3+2k2(1-k2)]518+(1+2k2)(1-k2)5110]k5 [316-(2-k2)sls+(1-k2)sl10]k3 Equation (C.46) can be written as where 2 4 6 8 10 12 Py0+Py28 +Py4s +Py68 +Py85 +Py10$ "I'Py125 = O P .. p -P 2 , Yo - nyO 0y 2 (C47) (C48) (C49) (C30) (C31) (C32) (C33) (C34) (C35) (C36) 186 ‘34 = Pny4-P4y2 (C57) Pys = Pnyfi'P6yz (c158) Pys = Pny8'P8y2 (C39) Ple = Pnle'PlOy2 (C60) Pyl2 = Pny12'P12y2 (C61) Equations (C36) and (C34) are, respectively, of the form f(k,s,s 1X) II C (C.62) f(k9598 1),) II C (C.63) Equations (C62) and (C.63) are solved simultaneously for s and $1, for every x and y input, using a quasi N ewton-Raphson numerical method. APPENDIX D APPENDIX D EVALUATING THE STRENGTHS OF THE POINT VORTICES We separate N [u -iv]k = 2 -in[1/(zk- zj)'1/(Zk' a2/'z‘j)1 (13.1) '=l _ I... into velocity components to get : A02 A . u = fi[(@% - 'g—EIEI‘YO—gn'] = f(vortex strength, geometnc factors) (D.2) A02 . v = Kifi - (%J = g(vortex strength, geometnc factors) (D3) where Kj = Vortex strength A02 = sz + yj2 a = Radius of circle AY = yk(A02) - azyj AX = xk(A02) - azxj YMYO =yk-yj XMXO =xk-xj (m2 + (A102 DENOMI = (XMXO)2 + (YMYO)2 DENOM2 j,k = indices for vortex location and point of interest respectively Now, (V9)k = vCosek - uSinGk ' I , xmxo (A02)(_A_X) _ (A02)(A1_Q YMYO . I303312110111 ' DENOM2 Icosek [DENOMZ ' DENOMI IsmekI 187 188 m0 _ (AozxA_X) ] _ [(AozxAx) _ YMYO ] - . “’hm fkj = [DENOMI DENOM2 C0591: DENOM2 DENOMI smek (D 5) From Equation (D.4) we express the tangential velocity at k points of interest due to j vortices as (V9)1 = Klf“ + Iczfu + K3f13 + + xjflj (V9)2 = K1f21 + K2f22 + x3f23 + + 1(5ij (V0)3 = Klf31 + K2f32 + K3f33 + ... + fif3j (Ve)k = Klfkl ‘I' 'szu + K3fB ‘I' ... + Kjfkj Or, in matrix form we get "(V191 ' ' f11 f12 f13 f1)“ "‘1 ' (V192 f21 f22 f23 f2j K2 . K (V193 = f31 f32 f33 f3) 3. or (Ve)k=ij 5. (D6) .. (V9)k I .. fkl sz fB fkja - IS - Solving for Kj there results ' f11 f21 f31 fkl‘ "(V191 ‘ “‘1 ‘ f12 f22 f32 sz (V192 "2 f13 f23 f33 fk3 (V193 = "3 or Flgjl (V9)k=sj (D?) L f11 f21 f31 fkj~ -(Ve)k . “‘5 ~ APPENDIX E APPENDIX E LIST OF COMPUTER PROGRAMS Page Flow Chart ................................................................................................. 192 Program Nomenclature for some Arrays and Control Parameters ............. 193 Program Swirl .......................................................................................... 197 DVM Library ............................................................................................. 199 - Subroutine INITP ..................................................................................... 199 - Subroutine INITPO ................................................................................... 201 - Subroutine IPV ......................................................................................... 203 - Subroutine GFCT OR ................................................................................ 204 - Subroutine JETMIX ................................................................................. 206 - Subroutine DVM ...................................................................................... 208 - Subroutine JET ......................................................................................... 209 - Subroutine FULL ...................................................................................... 210 - Subroutine VIC ......................................................................................... 212 — Subroutine INITPV ................................................................................... 213 189 190 - Subroutine BORDER ............................................................................... 215 - Subroutine GRIDVZ ................................................................................. 215 - Subroutine GRDV21 ................................................................................ 217 - Subroutine EXACT .................................................................................. 217 - Subroutine VORT ..................................................................................... 219 - Subroutine LAYGRD ............................................................................... 221 - Subroutine PNTV 2 ................................................................................... 224 - Subroutine PRNT ..................................................................................... 226 - Subroutine ADVECT ............................................................................... 228 - Subroutine UPDATE ................................................................................ 231 - Function PHI ........................................................................................ 231 DMAP Library ........................................................................................... 233 - Subroutine SSNl ...................................................................................... 233 - Subroutine XYUV .................................................................................... 234 - Subroutine NASZ ...................................................................................... 234 - Subroutine NASNWR .............................................................................. 236 - Subroutine NASJME ................................................................................ 240 - Subroutine NASFNB ................................................................................ 241 - Function NASEPS ................................................................................ 242 - Subroutine NASABX ............................................................................... 243 - Subroutine DECOM ................................................................................. 243 - Subroutine SUBST ................................................................................... 244 191 - Subroutine MOVEDD .............................................................................. 245 - Subroutine NORMSD .............................................................................. 245 - Subroutine EXCHSD ................................................................................ 246 - Subroutine VIPASD ................................................................................. 246 - Subroutine F ............................................................................................. 247 - Subroutine INTPLU ................................................................................. 249 - Subroutine CNDNSN ............................................................................... 250 - Function CN .......................................................................................... 251 - Function DN ........................................................................................... 251 - Function SN ........................................................................................... 251 - Function FUN CTN ................................................................................ 252 - Function GAUSS ................................................................................... 252 POISSON Library ..................................................................................... 254 - Subroutine COSGEN ................................................................................ 254 - Subroutine GENBUN ............................................................................... 255 - Subroutine HWSCRT ............................................................................... 256 - Subroutine MERGE .................................................................................. 262 - Subroutine POISD2 .................................................................................. 263 - Subroutine POISN2 .................................................................................. 269 - Subroutine POISP2 ................................................................................... 279 — Subroutine TRIX ...................................................................................... 282 - Subroutine TRI3 ....................................................................................... 283 192 ._ 8.3.5 _ , LE... _ _ PUa :9:sz TE; AmEm 213—00mm m2<-UO~E mMPDmEOU mOh Pz<=0>>OHm 193 C*** C*********************************************************************** C*** PROGRAM NOMENCLATURE FOR SOME ARRAYS AND CONTROL PARAMETERS ** C*********************************************************************** C*** PROGRAM ** C*********************************************************************** C*** it C*** SWIRL : ** C*** ** C*** ** C*** Parameters : ** C*** IDIMAX - Maximum number of gridpoints (for PET) ** C*** - 4*NCOUNT + [13+INTGER(LN(NCOUNT))]*NCOUNT ** C*** JERT = 1/0: Injector on/off ** C*** MAXPV - Maximum.number of vortex elements ** C*** NCOUNT - 181: Subgrid resolution for interpolation purpose ** C*** NCASE - 1/2: Circle/Square domain ** C*** NCOUNT - Number of grid points in the x and y directions ** C*** NDATA - Total no. of grid points(NDATA=NCOUNT**2 for a square) ** C*** NPRNT - 1/0: Print data : yes/no ** C*** NWALL - Number of wall test points ** C*** SWRL - 1.0/0.0: Swirling motion : yes/no ** C*** TOL - Tolerance ** C*** ** C*********************************************************************** C*** SUBROUTINES ** C*********************************************************************** C*** ** C*** INITP : ** C*** ** C*** ** C*** Parameters : ** C*** LCODE - Flag for monitoring array size ** C*** MBDCND - 1 Prescribed BC's on the left and right boundary ** C*** NBDCND - 1 Prescribed BC's on the lower and upper boundary ** C*** SGMA - Vortex cut-off radius ** C*** ** C*** ** C*** INITPO : ** C*** ** C*** ** C*** Arrays : ** C*** P - Interpolation points ** C*** SN - Jacobian elliptic function (S or $1) ** C*** UVDATA - (u,v) coordinates in square domain ** C*** XYDATA - (x,y) coordinates in circular domain ** C*** ** C*** Parameters : ** C*** RMODZ - Modulus. One of the arguments of the elliptic fctns. ** C*** RK - Complete elliptic integral of the first kind ** C*** RW - Radius of chamber - 1.0 cm. ** C*** ** C*** **** C*** IPV: **** 194 C*** ** C*** ** C*** Arrays : ** C*** GPNT - Vortex strength ** C*** RC - Radial locations of concentric circles ** C*** RPNT - Radial-coordinates of point vortices ** C*** S - Geometric shape factors ** C*** SINV - Inverted geometric shape factors ** C*** SRD - Dyer's velocity test points ** C*** SVD - Dyer's velocity data at test points ** C*** TPNT - Tangential-coordinates of point vortices ** C*** XPNT - x-coordinates of point vortices ** C*** YPNT - y-coordinates of point vortices ** C*** ** C*** Parameters : ** C*** GPVMN - Minimum vortex strength ** C*** GPVMX - Maxmum vortex strength ** C*** MXPV - Maximum allowable size for arrays NPVCI), RC() ** C*** NPV - Number of point vortices ** C*** NCCl - Number of concentric circles inside of circular do ** C*** NCC2 - Number of reference points for Dyer’s data ** C*** NPVC - Number of point vortices per concentric circle ** C*** ** C*** ** C*** JETMIX : ** C*** ** C*** ** C*** Parameters : ** C*** DSPM - Dsplcmt from injector lip to inner edge of shear layer ** C*** JETON = Turn injector on ** C*** JETOFF - Turn injector off ** C*** NXDIR - Number of grid cells in the x-direction ** C*** NYDIR - Number of grid cells in the y-direction ** C*** VIN - Inlet jet velocity ** C*** XA - x-location of injector lip ** C*** XAB — Horizontal distance between injector lip and seat ** C*** XB - x-location of injector seat ** C*** XPA — x-distance between vortex particle and lip ** C*** XPB - x-distance between vortex particle and seat ** C*** YAB - Vertical distance between injector lip and seat ** C*** YB - y-location of injector seat ** C*** YPA - y-distance between vortex particle and lip ** C*** YPB - y-distance between vortex particle and seat ** C*** ** C*** ** C*** JET ** C*** ** C*** ** C*** Arrays : ** C*** FTAG - Tags for jet vortices ** C*** ** C*** Parameters : ** C*** GPA - Strength of vortex issuing from valve lip ** C*** GPB - Strength of vortex issuing from valve seat ** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** c*** C*** c*** C*** C*** C*** C*** C*** C*** 195 XLA - x-location of vortex particle issuing from valve lip XLB - x-location of vortex particle issuing from valve seat YLA - y-location of vortex particle issuing from valve lip YLB - y—location of vortex particle issuing from valve seat XYUV Parameters : 61,62 - Initial guesses for the iteration process NAS routine UPNT - x-coordinate of vortex in square domain VPNT - y-coordinate of vortex in square domain HWSCRT : Arrays - BDA,BDB,BDC,BDD - Boundary data F - Streamfunction W - Work space Parameters : IDIMF - Number of grid points in the x-direction IDIMW - 2*NCOUNT*NCOUNT IDIMWP - IDIMW+2*NCOUNT JDIMF - Number of grid points in the y-direction GRDV21 : Arrays : UGl - u-component velocity in the circular domain UGlS - u-comp. vel. in the circular domain (single precision) UGZ - u—component velocity in the square domain VGl - v-component velocity in the circular domain VGlS - v-comp. vel. in the circular domain (single precision) VGZ - v-component velocity in the square domain Parameters : E11,E12- Matrices for velocity transformation EXACT Arrays : UGlE - u-component velocity in the circular domain VGlE - v-component velocity in the circular domain Parameters : NAP - Non-analytical grid points ** ** 'k* ** ** *‘k ** *‘k *‘k *‘k ** *‘k ** ** *‘k ** ** *‘k ** ** *‘k *‘k ** ** ** ** ** ** ** *‘k *‘k ** ** *‘k *‘k ** ** *‘k *‘k ** *‘k *‘k ** ** *‘k ** *‘k ** ** ** ** *‘k C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 196 ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** VORT Arrays : Q - Vorticity LAYGRD : Arrays : IBOX - Domain tags NG - Singularity tags UGN - u-component velocity on uniform grid in circular domain** VGN — v-component velocity on uniform grid in circular domain** XY - uniform (x,y) coordinates in circular domain PNTV2 Arrays U V u-component velocity at vortex location v-component velocity at vortex location ** ** ** ** ** ** ** C******************************************t**************************** C*** FUNCTIONS ** C*********************************************************************** C*** C*** C*** C*** C*** c*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** C*** C*** PHI PHI DHU,DHV- Interpolation function (Filter) Weight values CN,DN,SN CN, DN, SN - Jacobian elliptic functions GAUSS Gauss-Legendre quadrature for integration purposes ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** C*********************************************************************** C*** ** C*** 197 C*********************************************************************** PROGRAM SWIRL c*********************************************************************** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** THIS PROGRAM EMPLOYS DISCRETE VORTICES TO SIMULATE THE MIXING PROCESS PRIOR TO CHEMICAL REACTION IN THE FUEL-INJECTED CYLINDER MIXING DESCRIPTIONS ARE RESTRICTED TO A TWO OF A PISTON ENGINE. DIMENSIONAL, INCOMPRESSIBLE, INVISCID FORMULATION OF THE BULK DIFFUSION PROCESS. THE NUMERICAL SIMULATION IS BASED ON A VORTEX-IN-CELL FORMULATION. POISSON INVERSION IS ACCOMPLISHED USING THE FAST FOURIER TRANSFORM PACKAGE "FISHPAK" FROM THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH (NCAR), BOULDER, COLORADO. PLOTIT AND NCAR GRAPHICS PACKAGES ARE USED. ** ** ** ** ** ** *a ** ** ** ** C*********************************************************************** C**** C*** c*** C*** DOUBLE PRECISION G1,GZ,P,SN,SSl,RK,RMOD2,XX,YY DOUBLE PRECISION UGl,UGlE,UVDATA,VGI,VG1E,XYDATA PARAMETER(JERT-1,MCOUNT=181,NPRNT-l,SWRL-l.0) PARAMETER (NCASE'1,NCOUNT=61,G1=0.3719D0,G2=0.3719D0) PARAMETER(IDIMAX=1300,MAXPV-12000,NDATA-NCOUNT*NCOUNT) COMMON/ENDARI/MBDCND,NBDCND COMMON/GRIDl/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/GRIDZ/ULNGTH,UMN,UMX,VLNGTH,VMN,VMX COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/GRIDB/DX,DXINVT,DY,DYINVT,XMN,XMX,YMN,YMX COMMON/JETPAR/DJ,PMSIGN,TJ,VIN,XA,XB,XPA,XPB,YAB,YB,YPB COMMON/MAP/RK,RMOD2,XX,YY COMMON/MIXPAR/JETOFF,JETON,NXDIR,NYDIR COMMON/PTPARI/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNST1,HZ,RW,SGMAZI,TOL COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME COMMON/WORKSP/RWKSP(2727472) VIRTUAL UVDATA(2,NDATA),XYDATA(2,NDATA) VIRTUAL P(MCOUNT),SN(MCOUNT),SSl(2,NDATA) VIRTUAL XPNT(MAXPV),YPNT(MAXPV) VIRTUAL GPNT(MAXPV),RPNT(MAXPV),TPNT(MAXPV) VIRTUAL UGlE(NCOUNT,NCOUNT),VG1E(NCOUNT,NCOUNTI VIRTUAL UPNT(MAXPV),VPNT(MAXPV) DIMENSION F(NCOUNT,NCOUNT),W(IDIMAX) VIRTUAL BDA(NCOUNT),BDB(NCOUNT),BDC(NCOUNT),BDD(NCOUNT) VIRTUAL UGZ(NCOUNT,NCOUNT),VG2(NCOUNT,NCOUNT) VIRTUAL Q(NCOUNT,NCOUNT) VIRTUAL UN(MAXPV),VN(MAXPV),XPN(MAXPV),YPN(MAXPV) VIRTUAL UG1(NCOUNT,NCOUNT),VG1(NCOUNT,NCOUNT) VIRTUAL UGlS(NCOUNT,NCOUNT),VG1$(NCOUNT,NCOUNT) VIRTUAL UPV1(MAXPV),VPV1(MAXPV) VIRTUAL IBOX(NCOUNT,NCOUNT),NG(NCOUNT,NCOUNT),XY(2,NDATA) VIRTUAL FTAG(2,MAXPV),UGN(NCOUNT,NCOUNT),VGN(NCOUNT,NCOUNT) 198 OPEN(UNIT‘1,FILE"SWIRLS.CLNS',TYPEa’NEW') OPEN(UNIT'2,FILE-'SWIRLS.CLNS’,TYPE-'NEW') OPEN(UNIT'3,FILE-'SWIRLS.CLNS',TYPEI'NEW’) OPEN(UNIT=4,FILE-'SWIRLS.CLNS',TYPE-'NEW’) OPEN(UNIT-5,STATUs-'OLD',FILE-’SYS$OUTPUT') OPEN(UNIT'6,FILE-’SWIRLS.CLNS’,TYPEI'NEW') OPEN(UNIT-7,FILE"SWIRLS.CLNS',TYPE='NEW') OPEN(UNIT-8,FILE-'SWIRLS.CLNS',TYPEa’NEW') C*** C*** C*** Initialize parameters : c*** C*** CALL IWKIN(2727472) CALL INITP(GPNT,JERT,MAXPV,MCOUNT,NCOUNT,P,RPNT,SN,SSl, *SWRL,TPNT,UVDATA,XPNT,XYDATA,YPNT) C*** C*** C*** Perform.discrete vortex vortex calculations : c*** C*** DO 1 NTIME - NSTART,NSTOP C*** CALL DVM(BDA,BDB,BDC,BDD,F,FTAG,G1,G2,GPNT,IBOX,JERT, *MAXPV,MCOUNT,NCASE,NCOUNT,NG,NPRNT,P,Q,RPNT,SN,SSl,SWRL, *TPNT,UGN,UGI,UGlE,UGIS,UG2,UN,UPNT,UPV1,UVDATA,VGN,VG1, *VGlE,VGIS,VGZ,VN,VPNT,VPV1,VR,W,XPN,XPNT,XY,XYDATA,YPN,YPNT) C*** IF(LCODE.EQ.0.0R.NTIME.EQ.NSTOP) GO TO 2 C*** 1 CONTINUE C*** C*** C*** Close the debugging file C*** C*** 2 CLOSE(UNIT-1) CLOSE(UNIT-2) CLOSE(UNIT-3) CLOSE(UNIT-4) CLOSE(UNIT-6) CLOSE(UNIT-7) CLOSE(UNIT-8) c*** STOP 'NORMAL TERMINATION OF PROGRAM SWIRLS' END C*** 199 *** g*********************************************************************** C*** DVM LIBRARY CONTAINS THE FOLLOWING SUBROUTINES/FUNTIONS : *** C*** ** C*** SUBROUTINES : ** C*** 1. INITP - Inititializes parameters ** C*** -INITPO - Initializes grid and elliptic parameters ** C*** -IPV - Initializes coords. & strengths of the vortices ** C*** -GFCTOR - Computes the geometric shape factors ** C*** -LSGRR - Inverts the geometric shape factors [IMSL] ** C*** -JETMIX - Initializes jet and mixing parameters ** C*** 2. DVM - Makes calls vortex routines : ** C*** -JET - Simulates fluid injection ** C*** -FULL - Checks if cylinder fluid is "completely mixed" ** C*** -VIC - Calls Vortex-In—Cell formulation routines ** C*** -XYUV - Maps coordinates from circle to square [Dmaplib] ** C*** -INITPV- Initializes vorticity on square grid ** C*** -BORDER- Borders active points according to bndary conds. ** C*** -HWSCRT- Does the Fast Fourier Transform.[Poissonlib] ** C*** -GRIDV2- Computes velocities on a square grid ** C*** -GRDV21- Maps grid velocities : square -> circle ** C*** -EXACT- Computes exact velocities for non-analytical pts. ** C*** -VORT - Computes vorticity at grid points in uv-plane ** C*** -LAYGRD- Computes velocities on imposed uniform grid ** C*** -PNTV2 - Interpolates velocities at off grid locations ** C*** -PRNT - Prints data ** C*** -ADVECT - Advects vortices via Heun's time integrtn.scheme ** C*** -UPDATE - Updates parameters/array size ** C****** ***** C*** FUNCTIONS : ** C*** ** C*** 1. PHI - Evaluates the interpolation functions ** C*** ** C*** SUBROUTINE(S)/FUNCTIONS CALLED BUT NOT IN THIS LIBRARY : ** C*** SUBROUTINES : ** C*** 1. CNDNSN - Evaluates Jacobian elliptic functions CN,DN and SN ** C*** 2. SSNl - Evaluates Jacobian elliptic function SN ** C*** 3. XYUV - (See above) ** C*** 4. HWSCRT - (See above) ** C*** 5. PRNT - (See above) ** C****** ** C*** FUCTIONS : ** C*** 2. FUNCTN - Evaluates integrand of elliptic integral K(or K') ** C*** [1/sqrt(l-k**23in**2(phi)]/[1/sqrt(l-k'**23in**2(phi)]** C*** 3. GAUSS - Performs Gaussian integration of elliptic integrals ** c*** ** C*********************i************************************************t SUBROUTINE INITP(GPNT,JERT,MAXPV,MCOUNT,NCOUNT,P,RPNT,SN,SSl, *SWRL,TPNT,UVDATA,XPNT,XYDATA,YPNT) C*********************************************************************** c*** ** C*** THIS SUBROUTINE INITIALIZES PARAMETERS ' ** C*** ** C***************************************************************t******* C*** 200 DOUBLE PRECISION P,RK,RMOD2,SN,SSl,UVDATA,XX,XYDATA,YY COMMON/BNDARY/MBDCND,NBDCND COMMON/GRIDl/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/GRIDZ/ULNGTH,UMN,UMX,VLNGTH,VMN,VMX COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/GRIDB/DX,DXINVT,DY,DYINVT,XMN,XMX,YMN,YMX COMMON/JETPAR/DJ,PMSIGN,TJ,VIN,XA,XB,XPA,XPB,YAB,YB,YPB COMMON/MAP/RK,RMOD2,XX,YY COMMON/MIXPAR/JETOFF,JETON,NXDIR,NYDIR COMMON/PTPARl/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNST1,H2,RW,SGMAZI,TOL COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME EXTERNAL GAUSS,FUNCTN C*********************************************************************** c*** C*** C*** C*** C*** C*** c*** C*** C*** C*** C*** c*** C*** C*** C*** C*** C*** C*** c*** Grid, elliptic and interpolation parameters/arrays : CALL INITPO(DU,DV,DX,DY,MCOUNT,NCOUNT,P,PI,RW,SN,SSl,TWOPI, *UVDATA,XYDATA) Time step parameters : NFULL - 0 NSTART - 1 NSTOP - 1001 DTIME - 0.001 Vortex parameters : LCODE - l NPV - 0 GPVMN - 1.0E+06 GPVMX - -GPVMN CNSTl - 2.0/3.0 H2 - 3.0 SGMAO - 1.0E-06 SGMA - 0.05 SGMAZI - 1.0/(SGMA*SGMA) TOL - 10.*AMACH(4) Square domain parameters : MBDCND - 1 NBDCND - 1 DUINVT - 1 O/DU C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 201 DXINVT = 1.0/DX DVINVT - 1.0/DV DYINVT - 1.0/DY HDUINV = 0.5*DUINVT HDVINV 8 0.5*DVINVT UMN = -0.5*SNGL(RK) UMX - 0.5*SNGL(RK) ULNGTH - SNGL(RK) VLNGTH - ULNGTH VMN - UMN VMX - UMX XMN - -RW XMX - RW YMN - XMN YMX = XMX Initialize coordinates and strengths of vortices IF(SWRL.EQ.1.0) CALL IPV(GPNT,RPNT,TPNT,XPNT,YPNT) Initialize jet and mixing parameters : CALL JETMIX(JERT,NCOUNT) RETURN END C*********************************************************************** SUBROUTINE INITPO(DU,DV,DX,DY,MCOUNT,NCOUNT,P,PI,RW,SN,SSl,TWOPI, *UVDATA,XYDATA) C*********************************************************************** C*** C*** C*** C*** C*** THIS SUBROUTINE INITIALIZES GRID AND ELLIPTIC PARAMETERS (i) S(U) : U in P() and S in SN() (ii) Grid point coordinates:(U,V) in UVDATA(); (X,Y) in XYDATA() ** ** ** ** ** C*********************************************************************** C*** C*** IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON/MAP/RK,RMOD2,XX,YY COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE EXTERNAL GAUSS,FUNCTN VIRTUAL P(*),SN(*),SSl(2,*),UVDATA(2,*),XYDATA(2,*) C*********************************************************************** C*** c*** c*** C*** Grid parameters : C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 202 IMIN - 1 IMAX - NCOUNT IMID - (IMAX+IMIN)/2 ISIZE - IMAX-l JMID - IMID JMIN - IMIN JMAX - IMAX JSIZE - ISIZE Elliptic parameters RMODZ ' 0.5D0 PI - DACOS(-1.0D0) TWOPI - PI+PI RK - GAUSS( 0.0D0,0.5D0*PI,15,FUNCTN,RMOD2 ) DU - RK/DFLOAT(ISIZE) DV - DU Pretabulate S elliptic function : DUU - RK/DFLOAT(MCOUNT-1) DO 1 I - 1,MCOUNT P(I) - DUU*DFLOAT(I-1) - 0.5D0*RK CALL SSN1( P(I), RMOD2, SN(I) ) CONTINUE Initialize grid point coordinates (U,V) and (X,Y) DO 2 J ' JMIN,JMAX V - DV*DFLOAT(J-1) - 0.5D0*RK CALL CNDNSN( V, RMOD2, C1, D1, 31 ) C12 - C1*C1 012 - D1*D1 $12 - 51*31 DO 2 I ' IMIN,IMAX K - I+(J-JMIN)*JMAX U - DU*DFLOAT(I-1) - 0.5D0*RK CALL CNDNSN( U, RMOD2, C, D, S ) C2 - C*C D2 - D*D $2 - S*S DENOM-(C2*C12+D2*D12*SZ*812)*(1.0DO-D2*SIZ) XNUMPC*D*S*(C12*D12*(1.0D0-D2*512)+RMOD2*512*(C2*C12+D12*52)) YNUMPC1*D1*S1*(D2*(D12*SZ+C2*C12)+RMOD2*C2*SZ*(D2*512-1.0D0)) UVDATA(1,K) - U UVDATA(2,K) - V XYDATA(1,K) - XNUM/DENOM XYDATA(2,K) - YNUM/DENOM C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 203 IF(I.EQ.IMAX.AND.J.EQ.JMAX) THEN RW - DSQRT(XYDATA(1,K)**2+XYDATA(2,K)**2) DX ' 2.0D0*RW/DFLOAT(ISIZE) DY - DX END IF Store non-analytical points : IF(DABS(U).EQ.DABS(V)) THEN SSl(1,K) = $1 SSl(2,K) = S END IF CONTINUE RETURN END C*********************************************************************** SUBROUTINE IPV(GPNT,RPNT,TPNT,XPNT,YPNT) C****************************************************t****************** C*** c*** C*** C*** C*** THIS SUBROUTINE INITIALIZES THE COORDINATES AND STRENGTHS OF THE POINT VORTICES USING DYER'S VELOCITY DATA [33] AND MATRIX INVERSION ROUTINE FROM IMSL ROUTINES. ** ** ** ** ** C*********************************************************************** C*** C*** c*** C*** C C C C C*** C*** c*** C*** c*** C*** c*** C*** C*** PARAMETER(MXPV-200,NCC1-8,NCC2=10) Number of vortices per concentric circle DATA NPVC/4,6,6,6,28,32,15,8/ !@100 ms. DATA NPVC/4,6,6,8,31,32,15,8/ !@150 ms. DATA NPVC/4,6,6,13,31,32,15,8/ !@200 ms. DATA NPVC/4,6,6,16,34,34,18,8/ !@300 ms. DATA NPVC/4,6,6,6,26,30,16,8/ !@500 ms. Radial location of concentric circles DATA RC/0.1,0.2,0.30,0.40,0.50,0.60,0.70,0.80/ Velocity test points : DATA SRD/0.00,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90/ Dyer’s velocity data at test points : DATA SVD/0.00,7.80,24.SO,33.10,37.20,39.00,37.60,33.70,!@100 ms. *30.00,28.20/ 204 DATA SVD/0.00,6.20,14.70,21.60,25.10,23.9o,22.2o,20.4o,!@150 ms. *18.80,17.40/ DATA SVD/0.00,S.30,12.60,16.60,17.00,16.30,15.10,14.00, !@2oo ms. *13.10,12.40/ C DATA SVD/0.00,3.20,10.oo,11.40,10.3o,9.6o,9.00,8.40, !@300 ms. c *7.80,7.60/ C DATA SVD/0.00,2.80,4.90,5.50,5.20,4.90,4.70,4.50, !@500 ms. C *4.30,4.30/ 00 COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/PTPARl/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNSTl,H2,RW,SGMAZI,TOL VIRTUAL S(NCC2,MXPV),SINV(MXPV,NCC2) VIRTUAL NPVC(NCC1),RC(NCC1),SVD(NCC2),SRD(NCC2) C*** C*********************************************************************** C*** C*** Determine the total number of point vortices C*** C*** DO 20 N - 1,NCC1 NPV - va + NPVC(N) 20 CONTINUE C*** C*** C*** Compute the geometric factors, perform matrix inversion C*** and then compute the strengths of the point vortices C*** C*** CALL GFCTOR(GPNT,NCC2,NPV,NCC1,NCC2,NPV,NPVC,NCC2, *RC,S,SINV,SRD,SVD,RPNT,TPNT,XPNT,YPNT) c*** RETURN END C*** c*** C*********************************************************************** SUBROUTINE GFCTOR(GPNT,LS,LSINV,NCC1,NCC2,NCS,NPVC,NRS,RC,S,SINV, *SRD,SVD,RPNT,TPNT,XPNT,YPNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES THE GEOMETRIC FACTORS, PERFORMS MATRIX ** C*** INVERSION AND THEN COMPUTES THE STRENGHTS OF THE POINT VORTICES ** C*** *t C*********************************************************************** C*** COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/PTPARllGPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNST1,H2,RW,SGMAZI,TOL VIRTUAL GPNT(*),RPNT(*),TPNT(*),XPNT(*),YPNT(*) VIRTUAL NPVC(*),RC(*),S(LS,*),SINV(LSINV,*),SRD(*),SVD(*) EXTERNAL LSGRR C*** Ct****************t*t****************************t********************** 205 C*** C*** Compute the geometric shape factor relative to the reference point C*** C*** DO 30 I - 1,NCC2 J - 0 T - 0.0 x - SRD(I) Y - 0.0 C*** DO 20 N - 1,NCC1 DELTO ' TWOPI/NPVC(N) DO 10 K - 1,NPVC(N) J = J + 1 T0 - DELTO*FLOAT(K) x0 - RC(N)*COS(T0) Y0 - RC(N)*SIN(T0) IF(I.EQ.1) THEN RPNT(J) - SQRT(X0*X0+Y0*Y0) TPNT(J) - T0 XPNT(J) ' X0 YPNT(J) = Y0 END IF C*** C*** C*** Initialize parameters for the shape factors C*** C*** XMXO - X-XO YMYO - Y-YO A02 - X0*X0+Y0*Y0 AX - A02*X-(RW*RW)*XO AY - A02*Y-(RW*RW)*Y0 AX2 - AX*AX AY2 = AY*AY XMXOZ - XMXO*XMXO YMY02 - YMYO*YMYO DENOMl - XMX02+YMY02 DENOMZ - AX2+AY2 IF(DENOM2.EQ.0.0) DENOMZ - SGMAO SQRTD - SQRT(DENOM1) H1 - (3.0*CNST1*SQRTD)/SGMA C*** IF(SQRTD.LE.SGMA) THEN SU ' (A02*AY/DENOM2)-SGMA2I*(H2-H1)*YMYO SV - -(A02*AX/DENOM2)+SGMAZI*(H2-H1)*XMXO ELSE SU - (A02*AY/DENOM2)-(YMYO/DENOM1) SV 8 -(A02*AX/DENOM2)+(XMXO/DENOMl) END IF C*** ' S(I,J) - SV*COS(T)-SU*SIN(T) C*** 10 CONTINUE 20 CONTINUE 30 CONTINUE 206 C*** C*** C*** Invert the 3 matrix : C*** C*** CALL LSGRR(NRS,NCS,S,LS,TOL,IRANK,SINV,LSINV) C*** C*** C*** (i) Compute the strengths of the point vortices. C*** (ii) Check for the minimum and maximum strengths. C*** C*** DO 50 J = 1,NPV GPNT(J) - 0.0 DC 40 I - 1,NCC2 GPNT(J) - GPNT(J) + SVD(I)*SINV(J,I) 40 CONTINUE GPVMN - AMIN1(GPNT(J),GPVMN) GPVMX - AMAX1(GPNT(J),GPVMX) 50 CONTINUE C*** RETURN END C*** C*********************************************************************** SUBROUTINE JETMIX(JERT,NCOUNT) C*********************************************************************** C*** C*** C*** C*** C*** C*** THIS SUBROUTINE INITIALIZES JET AND MIXING SIMULATION PARAMETERS JDIR - 0/1 - Same/opposite direction of flow SIV - 0.0/1.0 - Different/same initial velocity. ** ** ** ti ** ** C*********************************************************************** c*** c*** PARAMETER(JDIR-1,SIv-1.0) COMMON/GRID 1/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/JETPAR/DJ,PMSIGN,TJ,VIN,XA,XB,XPA,XPB,YAB,YB,YPB COMMON/MIXPAR/JETOFF,JETON,NXDIR,NYDIR COMMON/PTPARZ/CNST1,H2,RW,SGMAZI,TOL COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME C*********************************************************************** C*** C*** C*** C*** Initialize jet parameters IF(JERT.EQ. JETON - NS JETOFF - NS END IF 1) THEN TART TOP C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** c*** C*** C*** C*** C*** c*** c*** c*** C*** C*** 207 Specify injector orifice diameter : DJ = DU Specify injector angle. CCW from horizontal : TJ = PI/3.0 IF(JDIR.EQ.1) THEN XB = -0.5*DJ*SIN(TJ) ELSE XB 3 0.5*DJ*SIN(TJ) END IF IF(XB.GE.0.0) PMSIGN = -1.0 IF(XB.LE.0.0) PMSIGN - 1.0 YB = 0.5*DJ*COS(TJ) YAB - DJ*COS(TJ) XAB = DJ*SIN(TJ) XA = PMSIGN*XAB+XB Compute the initial minimum allowable displacement of jet particles from injector orifice IF(TJ.EQ.(PI/2.0)) DSPM - 3.633021*DJ IF(TJ.EQ.(PI/2.4)) DSPM - 4.067413*DJ IF(TJ.EQ.(PI/3.0)) DSPM - 4.938183*DJ IF(TJ.EQ.(PI/4.0)) DSPM - 6.794676*DJ IF(TJ.EQ.(PI/6.0)) DSPM - 16.39000*DJ XPB = DSPM*COS(TJ) YPB - DSPM*SIN(TJ) pr - xps YPA = YPB Initialize jet velocity : IF(SIV.EQ.1.0) THEN VIN - 6.794676*DJ/DTIME !If velocity at 90 deg. is desired ELSE VIN = DSPM/DTIME END IF Initialize mixing parameters : 208 C*** NXDIR I 8 NYDIR I 8 C*** RETURN END c*** C*** c*********************************************************************** SUBROUTINE DVM(BDA,BDB,BDC,BDD,F,FTAG,G1,G2,GPNT,IBOX,JERT, *MAXPV,MCOUNT,NCASE,NCOUNT,NG,NPRNT,P,Q,RPNT,SN,SSl,SWRL, *TPNT,UGN,UG1,UG1E,UGIS,UGZ,UN,UPNT,UPV1,UVDATA,VGN,VG1,VGlE,VG1$, *VG2,VN,VPNT,VPV1,VR,W,XPN,XPNT,XY,XYDATA,YPN,YPNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE MAKES CALLS TO DISCRETE VORTEX ALGORITHMS ** C*** ** C*********************************************************************** C*** DOUBLE PRECISION G1,GZ,P,SN,SSl,UG1,UGlE,UVDATA,VG1,VGlE,XYDATA COMMON/MIXPAR/JETOFF,JETON,NXDIR,NYDIR COMMON/PTPARI/GPVMN,GPVMX,LCODE,NO,NPV COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME C*** C*** C*** Turn injector on/Off C*** C*** IF(JERT.EQ.1) THEN IF(NTIME.GE.JETON.AND.NTIME.LE.JETOFF) THEN CALL JET(FTAG,GPNT,MAXPV,RPNT,TPNT,XPNT,YPNT) IF(LCODE.EQ.0) GO TO 3 END IF C*** C*** C*** Check if cylinder fluid is “completely mixed" C*** C*** CALL FULL(IBOX,NCOUNT,RPNT,XPNT,YPNT) c*** c*** C*** Compute advection velocities C*** C*** CALL VIC(BDA,BDB,BDC,BDD,F,GI,GZ,GPNT,IBOX,1,JERT,MCOUNT,NCASE, *NCOUNT,NG,P,Q,SN,SSl,UGN,UGl,UGlE,UGlS,UGZ,UPNT,UPV1,UVDATA,VGN, *VGl,VGlE,VG1$,VGZ,VPNT,VPV1,W,XPNT,XY,XYDATA,YPNT) C*** C*** C*** Check if time to print data : C*** C*** IF(NPRNT.EQ.1) THEN CALL PRNT(F,FTAG,NCASE,NCOUNT,Q,RPNT,UGlS,UGZ,UGN,UVDATA,VGlS, 209 *VG2,VGN,XPNT,XY,XYDATA,YPNT) END IF C*** C*** C*** Initialize vortex elements for the next time step : C*** C*** CALL ADVECT(BDA,BDB,BDC,BDD,F,FTAG,Gl,GZ,GPNT,IBOX,JERT,MCOUNT, *NCASE,NCOUNT,NG,P,Q,RPNT,SN,SSl,TPNT,UGl,UGlE,UGlS,UGZ,UGN,UN, *UPV1,UPNT,UVDATA,VGl,VGlE,VG1$,VGZ,VGN,VN,VPV1,VPNT,XPN,XPNT,YPN, *YPNT,W,XY,XYDATA) C*** C*** C*** Updata parameters and monitor array size : C*** C*** CALL UPDATE(MAXPV) C*** RETURN END C*** C*********************************************************************** SUBROUTINE JET(FTAG,GPNT,MAXPV,RPNT,TPNT,XPNT,YPNT) c*********************************************************************** C*** ** C*** THIS SUBROUTINE SIMULATES THE FLUID INJECTION PROCESS ** C*** Ref.: ** C*** "Piston-Cylinder Fluid Motion Via Vortex Dynamics", W.T. Ashurst ** C*** Fluid Mechanics of Combustion Systems, ed. by T. Morel et. a1. ** C*** ** C*********************************************************************** C*** COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/JETPAR/DJ,PMSIGN,TJ,VIN,XA,XB,XPA,XPB,YAB,YB,YPB COMMON/PTPARl/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNST1,HZ,RW,SGMAZI,TOL COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME VIRTUAL GPNT(*),FTAG(2,*),RPNT(*),TPNT(*),XPNT(*),YPNT(*) C*** C*********************************************************************** C*** SEAT CALCULATIONS (B) ** C*********************************************************************** C*** C*** Compute location of particle issuing from valve seat (region b) : C*** c*** XLB I XB+PMSIGN*XPB YLB I (YB+YPB)-RW GPB I PMSIGN*0.5*(VIN*VIN*DTIME) GPB I GPB/TWOPI C*** c*** C*** Update parameters and store values in point vortex arrays : C*** 210 C*** NPV I NPV + 1 IF(NPV.GT.MAXPV) GO TO 4 XPNT(NPV) I XLB YPNT(NPV) I YLB GPNT(NPV) I GPB RPNT(NPV) I SQRT(XLB*XLB+YLB*YLB) TPNT(NPV) I ATAN2(YLB,XLB) NJP I NPV c*** Ct********************************************************************** C*** LIP CALCULATIONS (A) ** C**********‘k‘k*********************************************************** C*** C*** Compute location of particle issuing from valve lip (region a) : C*** C*** XLA - XA+PMSIGN*XPA YLA - YLB-YAB C*** C*** C*** Update parameters and store values in point vortex arrays : C*** C*** NPV I NPV + 1 IF(NPV.GT.MAXPV) GO TO 4 XPNT(NPV) I XLA YPNT(NPV) I YLA GPNT(NPV) I IGPB RPNT(NPV) I SQRT(XLA*XLA+YLA*YLA) TPNT(NPV) I ATAN2(YLA,XLA) C*** C*** C*** Tag the jet vortices as fuel particles C*** C*** 2 DO 3 N I NJP,NPV FTAG(1,N) I 1.0 3 CONTINUE C*** 4 IF(NPV.GT.MAXPV) LCODE I 0 RETURN END c*** C*** C*********************************************************************** SUBROUTINE FULL(FTAG,IBOX,NCOUNT,RPNT,XPNT,YPNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE DETERMINES IF DOMAIN IS FILLED WITH VORTICES ** c*** ** Ctttt***§*************************************************************** c*** PARAMETER(NSKMI1) 211 COMMON/MIXPAR/JETOFF,JETON,NXDIR,NYDIR COMMON/PTPARl/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNST1,H2,RW,SGMAZI,TOL COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME VIRTUAL FTAG(2,*),IBOX(NCOUNT,*),RPNT(*),XPNT(*),YPNT(*) C*** C*********************************************************************** c*** C*** Initialize parameters : C*** C*** NCELLS I 0 R0 I 2.0*RW/FLOAT(NCIRS) DO 7 J I 1,NRAYS DO 7 I I 1,NCIRS IBOX(I,J) I 0 7 CONTINUE C*** C*** C*** Tag grid cells within Circle C*** C*** DO 1 J I 1,NRAYS RY I IRW+R0*FLOAT(J-l) RYP I RY+R0 HRYP I RY+(0.5*R0) DO 1 I I 1,NCIRS RX I -RW+R0*FLOAT(I-l) RXP I RX+R0 HRXP I RX+(0.5*R0) RCC I SQRT(HRXP*HRXP+HRYP*HRYP) C*** IF(RCC.LE.RW) THEN IBOX(I,J) I 1 NCELLS I NCELLS + 1 END IF C*** 1 CONTINUE C*** C*** C*** Identify cell C*** C*** DO 3 J I 1,NRAYS RY I -RW+R0*FLOAT(J-1) RYP I RY+R0 D0 3 I I 1,NCIRS NYES I 0 RX I IRW+RO*FLOAT(I-1) RXP I RX+R0 C*** C*** C*** Tag cell with at least NO vortex/vortices c*** 212 C*** NO I 2 IF(IBOX(I,J).EQ.1) THEN DO 2 N I 1,NPV IF(FTAG(1,N).EQ.1.0) THEN IF((XPNT(N).GE.RX.AND.XPNT(N).LE.RXP).AND. *(YPNT(N).GE.RY.AND.YPNT(N).LE.RYP)) THEN NYES I NYES + 1 IF(NYES.EQ.NO) THEN IBOX(I,J) I 2 GO TO 3 END IF END IF END IF 2 CONTINUE END IF 3 CONTINUE C*** C*** C*** Check if each cell contains at least NO vortex/vortices C*** C*** NYES I 0 DO 4 J I 1,NRAYS DO 4 I I 1,NCIRS IF(IBOX(I,J).EQ.1) GO TO 70 IF(IBOX(I,J).EQ.2) NYES I NYES + 1 4 CONTINUE C*** 5 IF(NYES.EQ.NCELLS) THEN NFULL I 1 END IF C*** 70 RETURN END C*** c*** C*********************************************************************** SUBROUTINE VIC(BDA,BDB,BDC,BDD,F,G1,G2,GPNT,IBOX,IFLG,JERT, *MCOUNT,NCASE,NCOUNT,NG,P,Q,SN,SSl,UGN,UG1,UG1E,UGlS,UGZ,UPNT,U, *UVDATA,VGN,VG1,VG1E,VG1$,VG2,VPNT,V,W,XPNT,XY,XYDATA,YPNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES GRID VELOCITIES IN SQUARE DOMAIN USING ** C*** VORTEX-IN-CELL (VIC) FORMULATION. OFF-GRID VELOCITIES CALCULATED ** C*** USING QUADRATIC SPLINE INTERPOLATION TECHNIQUE (COUET’S FILTER). ** c*** ** c*********************************************************************** C*** PARAMETER(NAPI1) DOUBLE PRECISION Gl,G2,P,SN,SSl,UGl,UGlE,UVDATA,VGl,VG1E,XYDATA COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME C*** c*********************************************************************** 213 C*** IF(NPV.EQ.0) RETURN C*** IF(NCASE.EQ.2.AND.NTIME.GT.NSTART) GO TO 1 CALL XYUV(G1,GZ,1,MCOUNT,0,P,SN,SSl,UPNT,UVDATA,VPNT,XPNT, *XYDATA,YPNT) IF(NAP.EQ.0) GO TO 2 1 CALL INITPV(F,GPNT,NCOUNT,UPNT,VPNT) CALL BORDER(F,NCOUNT) CALL HWSCRT(BDA,BDB,BDC,BDD,0.0,F,NCOUNT,PERTRB,IERROR,W) CALL GRIDV2(F,NCOUNT,UG2,VG2) C*** 2 IF(NCASE.EQ.1) THEN CALL GRDV21(GPNT,NAP,NCOUNT,UG1,UGlE,UGlS,UG2,VG1,VG1E,VG1$,VG2, *XPNT,XYDATA,YPNT) C*** IF(IFLG.EQ.1) THEN IF(NAP.EQ.0) GO TO 3 CALL VORT(NCOUNT,UG2,VG2,Q) 3 CALL LAYGRD (IBOX, NCOUNT, NG, UGlS , UGN, XY, VGl S , VGN, XYDATA) END IF C*** CALL PNTV2(1,NCOUNT,0,UG15,UPNT,U,VG1$,VPNT,V) ELSE CALL PNTV2(1,NCOUNT,0,UG2,UPNT,U,VGZ,VPNT,V) END IF C*** RETURN END c*** c*** C*********************************************************************** SUBROUTINE INITPV(F,GPNT,NCOUNT,UPNT,VPNT) C****t****************************************************************** C*** ** C*** THIS SUBROUTINE SPREADS VORTICITY TO GRID POINTS ** C*** it C*********************************************************************** c*** COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/GRIDl/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/GRIDZ/ULNGTH,UMN,UMX,VLNGTH,VMN,VMX COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/PTPARI/GPVMN,GPVMX,LCODE,NO,NPV DIMENSION F(NCOUNT,*) VIRTUAL GPNT(*),UPNT(*),VPNT(*) EXTERNAL PHI C*** C*********************************************************************** C*** C*** Initialize the charge array : C*** C*** DO 1 J = JMIN,JMAX 214 DO 1 I I IMIN,IMAX F(I,J) I 0.0 1 CONTINUE C*** C*** C*** Locate the point vortices : C*** C*** FACTR I ITWOPI*DUINVT*DVINVT DO 60 N I 1,NPV GPOINT I FACTR*GPNT(N) HUP I (UPNT(N)-UMN)*DUINVT+1. HVP I (VPNT(N)-VMN)*DVINVT+1. ICOL I IFIX(HUP) JROW I IFIX(HVP) DCOL I HUP-FLOAT(ICOL) DROW I HVP-FLOAT(JROW) C*** IPV I IFIX(HUP+.5)-2 JPV I IFIX(HVP+.5)-2 C*** DO 50 II I 1,3 IDELTA I IPV + II HUC I FLOAT(IDELTA) DHU I ABS(HUP-HUC) C*** C*** C*** Check if outside the computational domain in the u-direction : c*** C*** IF(IDELTA.LT.IMIN) IDELTA I IMIN IF(IDELTA.GT.IMAX) IDELTA I IMAX c*** DO 40 JJ I 1,3 JDELTA I JPV + JJ HVR I FLOAT(JDELTA) DHV I ABS(HVP-HVR) c*** c*** C*** Check if outside the computational domain v-direcction : C*** c*** IF(JDELTA.LT.JMIN) JDELTA I JMIN IF(JDELTA.GT.JMAX) JDELTA I JMAX C*** F(IDELTA,JDELTA) I F(IDELTA,JDELTA) + GPOINT*PHI(DHU,DHV) 40 CONTINUE 50 CONTINUE 60 CONTINUE C*** RETURN END C*** C*** 215 C*********************************************************************** SUBROUTINE BORDER(F,NCOUNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE SPECIFIES THE BOUNDARY CONDITION ** C*** ** C*********************************************************************** C*** COMMON/NSIZE /IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE DIMENSION F(NCOUNT,*) C*** C*********************************************************************** C*** C*** Initialize the border of the E-array : C*** C*** DO 1 I I IMIN,IMAX,ISIZE DO 1 J I JMIN,JMAX F(I,J) I 0.0 1 CONTINUE DO 2 J I JMIN,JMAX,JSIZE DO 2 I I IMIN+1,ISIZE F(I,J) I 0.0 2 CONTINUE C*** RETURN END C*** c*** C*********************************************************************** SUBROUTINE GRIDV2(F,NCOUNT,UGZ,VG2) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES VELOCITIES AT GRID POINTS ** C*** IN THE UV-PLANE USING A CENTERED DIFFERENCE TECHNIQUE ** C*** WHICH ASSURES LOCAL CONSERVATION OF MASS ** c*** ** C*********************************************************************** C*** COMMON/GRIDl/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE DIMENSION F(NCOUNT,*) VIRTUAL UG2 (NCOUNT , *) , VG2 (NCOUNT, *) C*** C*********************************************************************** C***— C*** Initialize arrays C*** C*** DO 50 J - JMIN,JMAX DO 50 I - IMIN,IMAX UG2(I,J) - 0.0 VG2(I,J) a 0.0 50 CONTINUE 216 C*** C*** C*** Interior grid points C*** C*** DO 1 J I JMIN+1,JSIZE DO 1 I I IMIN+1,ISIZE UG2(I,J) I UGZ(I,J)+(F(I,J+1)-F(I,J-1))*HDVINV VG2(I,J) I VGZ(I,J)+(F(I-1,J)-F(I+1,J))*HDUINV 1 CONTINUE C*** C*** C*** Lower boundary grid points C*** C*** DO 3 I -IMIN,IMAX . UGZ(I,JMIN) I UGZ(I,JMIN)+2.0*F(I,JMIN+1)*HDVINV IF(I.EQ.IMIN) VG2(I,JMIN)IVGZ(I,JMIN)-2.0*F(I+1,JMIN)*HDUINV IF(I.EQ.IMAX) VG2(I,JMIN)IVGZ(I,JMIN)+2.0*F(I-1,JMIN)*HDUINV IF(I.GT.IMIN.AND.I.LT.IMAX) VGZ(I,JMIN) I VGZ(I,JMIN) + * (F(I-1,JMIN)-F(I+1,JMIN))*HDUINV C*** 3 CONTINUE C*** C*** C*** Upper boundary grid points : C*** C*** DO 4 I I IMIN,IMAX UGZ(I,JMAX) I UG2(I,JMAX)-2.0*F(I,JMAXI1)*HDVINV IF(I.EQ.IMIN) VGZ(I,JMAX)IVGZ(I,JMAX)-2.0*F(I+1,JMAX)*HDUINV IF(I.EQ.IMAX) VG2(I,JMAX)IVGZ(I,JMAX)+2.0*F(I-1,JMAX)*HDUINV IF(I.GT.IMIN.AND.I.LT.IMAX) VGZ(I,JMAX) I VGZ(I,JMAX) + * (F(I-1,JMAX)-F(I+1,JMAX))*HDUINV 4 CONTINUE C*** C*** C*** Left boundary grid points : C*** C*** DO 5 J I JMIN+1,JSIZE UGZ(IMIN,J) I UGZ(IMIN,J)+(F(IMIN,J+1)-F(IMIN,J-1))*HDVINV VGZ(IMIN,J) I VG2(IMIN,J)-2.0*F(IMIN+1,J)*HDUINV 5 CONTINUE C*** C*** C*** Right boundary grid points C*** C*** DO 6 J I JMIN+1,JSIZE UG2 (IMAX, J) I UG2 (IMAX, J) + (F (IMAX, J+1) IF (IMAX, J-1))*HDVINV VG2(IMAX,J) I VG2(IMAX,J)+2.0*F(IMAX-1,J)*HDUINV 6 CONTINUE C*** 217 RETURN END C*** C*** C*********************************************************************** SUBROUTINE GRDV21(GPNT,NAP,NCOUNT,UG1,UG1E,UGlS,UG2,VG1,VG1E, *VG1$,VG2,XPNT,XYDATA,YPNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES GRID VELOCITIES IN THE CIRCLE DOMAIN ** C*** ** C*********************************************************************** C*** DOUBLE PRECISION E11,E12,UG1,UG1E,VG1,VG1E,XYDATA COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE VIRTUAL UG2 (NCOUNT, *) ,VGZ (NCOUNT, *) VIRTUAL UGl (NCOUNT, *) , UG1E(NCOUNT, *) , UGlS (NCOUNT, *) VIRTUAL VGl (NCOUNT, *) ,VG1E(NCOUNT, *) , VGIS (NCOUNT, *) C*** C*********************************************************************** C*** C*** Compute exact solution C*** C*** CALL EXACT(GPNT,NCOUNT,NAP,UGlE,UGlS,VG1E,VGlS,XPNT,XYDATA,YPNT) IF(NAP.EQ.0) RETURN C*** OPEN(UNITI3,FILEI'MATRCES.CLNS’,TYPEI'OLD’) DO 2 J I JMIN,JMAX DO 2 I I IMIN,IMAX READ(3,1) E11,E12 1 FORMAT(' ',T2,2(2X,1PD15.7)) UGl(I,J) I E11*DBLE(UG2(I,J))+E12*DBLE(VG2(I,J)) VG1(I,J) I IE12*DBLE(UG2(I,J))+E11*DBLE(VG2(I,J)) IF(NAP.EQ.1) THEN IF(I.EQ.IMIN.OR.I.EQ.IMAX) THEN IF(J.EQ.JMIN.OR.J.EQ.JMAX) THEN UG1(I,J) I UG1E(I,J) VG1(I,J) I VGlE(I,J) END IF END IF END IF UGIS(I,J) I SNGL(UG1(I,J)) VG1$(I,J) I SNGL(VG1(I,J)) 2 CONTINUE CLOSE(UNITI3) C*** RETURN END C*** C*** C*********************************************************************** SUBROUTINE EXACT(GPNT,NCOUNT,NAP,UG1E,UGIS,VGlE,VG1$,XPNT,XYDATA, *YPNT) 218 C*********************************************************************** C*** it C*** THIS SUBROUTINE COMPUTES EXACT VELOCITIES ** C*** USING A MODIFIED MILNE-THOMSON CIRCLE THEOREM ** C*** ** c*********************************************************************** c*** DOUBLE PRECISION UGlE,VG1E,XYDATA COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/PTPAR1/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNST1,H2,RW,SGMAZI,TOL VIRTUAL GPNT(*),XYDATA(2,*) VIRTUAL UG1$(NCOUNT,*),VGlS(NCOUNT,*) VIRTUAL UGlE (NCOUNT, *) ,VG1E(NCOUNT, *) , XPNT (*) , YPNT (*) C*** C*********************************************************************** C*** C*** Initialize the arrays : C*** C*** DO 1 J I JMIN,JMAX DO 1 I I IMIN,IMAX UG1E(I,J) I 0.0D0 VGlE(I,J) I 0.0D0 UGIS(I,J) I 0.0 VGlS(I,J) I 0.0 1 CONTINUE C*** IF(NAP.EQ.0) THEN JNC I 1 INC I 1 ELSE JNC I JSIZE INC I ISIZE END IF C*** DO 3 N I 1, NPV A02 I XPNT(N)*XPNT(N)+YPNT(N)*YPNT(N) C*** DO 2 J I JMIN,JMAX,JNC DO 2 I I IMIN,IMAX,INC K I I + (J-JMIN)*JMAX C*** C*** C*** Initialize parameters for u,v velocities : C*** C*** AX I A02*SNGL(XYDATA(1,K))-(RW*RW)*XPNT(N) AY I A02*SNGL(XYDATA(2,K))I(RW*RW)*YPNT(N) AX2 I AX*AX AY2 I AY*AY XMXO I SNGL(XYDATA(1,K))-XPNT(N) YMYO I SNGL(XYDATA(2,K))IYPNT(N) C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 219 XMXOZ I XMXO*XMXO YMYOZ I YMYO*YMYO DENOMl I XMX02+YMY02 DENOM2 I AX2+AY2 IF(DENOM2.EQ.0.0) DENOM2 I SGMAO SQRTD I SQRT(DENOM1) H1 I (3.0*CNST1*SQRTD)/SGMA Remove singularity according to Bald & Del Prete’s formulation IF(SQRTD.LE.SGMA) THEN UTMP I GPNT(N)*( (A02*AY/DENOM2)-SGMA2I*(H2-H1)*YMYO) VTMP I GPNT(N)*(-(A02*AX/DENOM2)+SGMA2I*(H2-H1)*XMXO) ELSE UTMP I GPNT(N)*( (A02*AY/DENOM2)-(YMYO/DENOM1)) VTMP I GPNT(N)*(-(A02*AX/DENOM2)+(XMXO/DENOMl)) END IF Compute the velocity components : N UG1E(I,J) I UGlE(I,J)+DBLE(UTMP) VG1E(I,J) I VG1E(I,J)+DBLE(VTMP) IF(NAP.EQ.0) THEN UG1$(I,J) I SNGL(UG1E(I,J)) VGlS(I,J) I SNGL(VG1E(I,J)) END IF CONTINUE CONTINUE RETURN END C*********************************************************************** SUBROUTINE VORT(NCOUNT,UGZ,VG2,Q) C*********************************************************************** C*** C*** C*** C*** C*** THIS SUBROUTINE COMPUTES VORTICITY AT GRID POINTS IN THE UV-PLANE USING A CENTERED DIFFERENCE TECHNIQUE WHICH ASSURES LOCAL CONSERVATION OF VORTICITY ** ** ** ** ** C*********************************************************************** C*** C*** COMMON/GRIDI/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE VIRTUAL Q(NCOUNT,*),UG2(NCOUNT,*),VGZ(NCOUNT,*) C***t******************************************************************* 220 C*** C*** Initialize array C*** C*** DO 50 J I JMIN,JMAX DO 50 I I IMIN,IMAX Q(I,J) - 0.0 50 CONTINUE C*** C*** C*** Interior grid points : C*** C*** DO 1 J I JMIN+1,JSIZE DO 1 I I IMIN+1,ISIZE DUDY I (UGZ(I,J+1)-UGZ(I,J-1))*HDVINV DVDX I (VGZ(I+1,J)-VGZ(I-1,J))*HDUINV Q(I,J) I DVDX-DUDY 1 CONTINUE C*** C*** C*** Lower boundary grid points C*** C*** DO 3 I I IMIN,IMAX DUDY I 2.0*UGZ(I,JMIN+1)*HDVINV IF(I.EQ.IMIN) DVDX I 2.0*VGZ(I+1,JMIN)*HDUINV IF(I.EQ.IMAX) DVDX I I2.0*VGZ(I-1,JMIN)*HDUINV IF(I.GT.IMIN.AND.I.LT.IMAX) DVDX I (VGZ(I+1,JMIN)-VGZ(I-1,JMIN))* *HDUINV Q(I,JMIN) I DVDX-DUDY C*** 3 CONTINUE C*** C*** C*** Upper boundary grid points : C*** C*** DO 4 I I IMIN,IMAX DUDY I I2.0*UGZ(I,JMAX-1)*HDVINV IF(I.EQ.IMIN) DVDX I 2.0*VGZ(I+1,JMAX)*HDUINV IF(I.EQ.IMAX) DVDX I I2.0*VGZ(I-1,JMAX)*HDUINV IF(I.GT.IMIN.AND.I.LT.IMAX) DVDX I (VGZ(I+1,JMAX)-VGZ(I-1,JMAX))* *HDUINV Q(I,JMAX) I DVDX-DUDY 4 CONTINUE C*** C*** C*** Left boundary grid points : C*** c*** ' DO 5 J - JMIN+1,JSIZE DUDY - (UG2(IMIN,J+1)IUG2(IMIN,J-1))*HDVINV DVDX - 2.0*VGZ(IMIN+1,J)*HDUINV 221 Q(IMIN,J) = DVDX-DUDY 5 CONTINUE C*** C*** C*** Right boundary grid points : C*** C*** DO 6 J I JMIN+1,JSIZE DUDY I (UG2(IMAX,J+1)-UG2(IMAX,J-1))*HDVINV DVDX I I2.0*VG2(IMAX-1,J)*HDUINV Q(IMAX,J) I DVDX-DUDY 6 CONTINUE C*** RETURN END C*** C*** C*********************************************************************** SUBROUTINE LAYGRD(IBOX,NCOUNT,NG,UGlS,UGN,XY,VG1$,VGN,XYDATA) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES VELOCITIES ON UNIFORM ** C*** GRID POINTS SUPERIMPOSED ON THE CIRCULAR DOMAIN ** C*** ** C*********************************************************************** C*** DOUBLE PRECISION XYDATA COMMON/GRID3/DX,DXINVT,DY,DYINVT,XMN,XMX,YMN,YMX COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME EXTERNAL PHI VIRTUAL UG1$(NCOUNT,*),UGN(NCOUNT,*),XY(2,*) VIRTUAL VGlS(NCOUNT,*),VGN(NCOUNT,*),XYDATA(2,*) VIRTUAL IBOX(NCOUNT,*),NG(NCOUNT,*) C*** C*********************************************************************** C*** C*** Initialize the arrays/parameters C*** C*** DO 1 J I JMIN,JMAX DO 1 I I IMIN,IMAX UGN(I,J) I 0.0 VGN(I,J) I 0.0 NG(IIJ) ' 0 1 CONTINUE C*** C*** C*** Assign radial location of singularity (RS) C*** and velocity reduction factor (RF) C*** C*** RS ll c3<> C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** c*** c*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 222 Tag grid points outside circle IF(NTIME.EQ.NSTART) THEN DO 2 J I JMIN,JMAX Y I DY*FLOAT(J-JMIN)+YMN DO 2 I I IMIN,IMAX X I DX*FLOAT(I-IMIN)+XMN K I I + (J-JMIN)*JMAX IBOX(I,J) I 0 XY(1,K) I X XY(21K) " Y R I SQRT(X*X+Y*Y) IF(R.GT.XMX) IBOX(I,J) I 1 CONTINUE END IF Interpolate velocities to grid points DO 70 NS I 1,2 DO 60 J I JMIN,JMAX DO 60 I I IMIN,IMAX K I I + (JIJMIN)*JMAX Locate vortex position HXP I (SNGL(XYDATA(1,K))-XMN)*DXINVT+1. HYP I (SNGL(XYDATA(2,K))IYMN)*DYINVT+1. ICOL I IFIX(HXP) JROW I IFIX(HYP) DCOL I HXP-FLOAT(ICOL) DROW I HYP-FLOAT(JROW) IPV I IFIX(HXP+.S)-2 JPV I IFIX(HYP+.5)-2 Determine argument, dx, for the weighting function : DO 50 II I 1,3 IDELTA I IPV + II IF(NS.EQ.2) GO TO 3 HXC I FLOAT(IDELTA) DHX I ABS(HXP-HXC) C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** c*** c*** C*** C*** c*** C*** c*** 223 Check if outside the computational domain IF(IDELTA.LT.IMIN) IDELTA I IMIN IF(IDELTA.GT.IMAX) IDELTA I IMAX Determine argument, dy, for the weighting function : DO 40 JJ I 1,3 JDELTA I JPV + JJ IF(NS.EQ.2) GO TO 4 HYR I FLOAT(JDELTA) DHY I ABS(HYP-HYR) Check if outside the computational domain r— IF(JDELTA.LT.JMIN) JDELTA I JMIN IF(JDELTA.GT.JMAX) JDELTA I JMAX Check if grid point is within circle IF(IBOX(IDELTA,JDELTA).EQ.1) GO TO 40 IF(NS.EQ.2) GO TO 5 Interpolate velocities to grid points UGN(IDELTA,JDELTA)IUGN(IDELTA,JDELTA)+UG18(I,J)*PHI(DHX,DHY) VGN(IDELTA,JDELTA)IVGN(IDELTA,JDELTA)+VGIS(I,J)*PHI(DHX,DHY) IF(NS.EQ.2) THEN Modify velocities within singularity zone Check (x,y) coordinates of grid points X I (IDELTA-1)*DX + XMN Y I (JDELTA-1)*DY + YMN Check radial location of grid points 224 C*** R I SQRT(X*X+Y*Y) C*** IF(R.GE.RS) THEN C*** C*** C*** Check angular location of grid points C*** C*** T I ATAN2D(Y,X) IF(Y.LT.0.0) T I T + 360.0 C*** IF((T.GE.35.0.AND.T.LE.55.0) .OR. *(T.GE.125.0.AND.T.LE.145.0) .OR. *(T.GE.215.0.AND.T.LE.235.0) .OR. *(T.GE.305.0.AND.T.LE.325.0)) THEN C*** C*** C*** Check if duplicating. If not tag grid point C*** C*** IF(NG(IDELTA,JDELTA).EQ.1) GO TO 40 NG(IDELTA,JDELTA) I 1 C*** C*** C*** Modify grid velocity : C*** C*** UGN(IDELTA,JDELTA) I RF*UGN(IDELTA,JDELTA) VGN(IDELTA,JDELTA) I RF*VGN(IDELTA,JDELTA) C*** END IF END IF END IF C*** 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C*** RETURN END C*** c*** C*********************************************************************** SUBROUTINE PNTV2(IFLG,NCOUNT,NTP,UGC,UPT,U,VGC,VPT,V) C*********************************************************************** C*** ** C*** THIS SUBROUTINE CALCULATES VELOCITIES AT OFF GRID LOCATIONS ** c*** ** C*********************************************************************** C*** COMMON/GRIDl/DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/GRIDZ/ULNGTH,UMN,UMX,VLNGTH,VMN,VMX 225 COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/PTPARl/GPVMN,GPVMX,LCODE,NO,NPV VIRTUAL UPT(*),VPT(*) VIRTUAL UGC(NCOUNT,*),U(*),VGC(NCOUNT,*),V(*) EXTERNAL PHI C*** C*********************************************************************** C*** IF(IFLG.EQ.1) NP I NPV IF(IFLG.EQ.2) NP I NTP DO 90 N I 1,NP U(N) I 0. V(N) I 0. UTEMP I 0. VTEMP I 0. HUP I (UPT(N)-UMN)*DUINVT+1. HVP I (VPT(N)IVMN)*DVINVT+1. C*** ICOL I IFIX(HUP) JROW I IFIX(HVP) DCOL I HUP-FLOAT(ICOL) DROW I HVP-FLOAT(JROW) C*** IPV I IFIX(HUP+.5)-2 JPV I IFIX(HVP+.5)-2 C*** DO 80 II I 1,3 IDELTA I IPV + II HUC I FLOAT(IDELTA) DHU I ABS(HUP-HUC) C*** C*** C*** Check if outside the computational domain : C*** C*** IF(IDELTA.LT.IMIN) IDELTA I IMIN IF(IDELTA.GT.IMAX) IDELTA I IMAX C*** DO 70 JJ I 1,3 JDELTA I JPV + JJ HVR I FLOAT(JDELTA) DHV I ABS(HVP-HVR) C*** C*** C*** Check if outside the computational domain : C*** C*** IF(JDELTA.LT.JMIN) JDELTA I JMIN IF(JDELTA.GT.JMAX) JDELTA I JMAX C*** UTEMP I UTEMP + UGC(IDELTA,JDELTA)*PHI(DHU,DHV) VTEMP I VTEMP + VGC(IDELTA,JDELTA)*PHI(DHU,DHV) 70 CONTINUE 80 CONTINUE 226 C*** C*** C*** Store the velocity C*** C*** U(N) - UTEMP V(N) - VTEMP 90 CONTINUE C*** RETURN END C*** C*** C*********************************************************************** SUBROUTINE PRNT(F,FTAG,NCASE,NCOUNT,Q,RPNT,UGlS,UG2,UGN, *UVDATA,VGlS,VG2,VGN,XPNT,XY,XYDATA,YPNT) C*********************************************************************** C*** ** C*** THIS SUBROUTINE PRINTS RESULTS OF VORTEX CALCULATION ** C*** ** C*********************************************************************** c*** DOUBLE PRECISION RK,RMOD2,UVDATA,XX,YY,XYDATA PARAMETER(NGRD-0,NPos-o,TANG-1.0,XDIR-1.0) COMMON/CONSNT/PI,SGMA,SGMAO,TWOPI COMMON/JETPAR/DJ,PMSIGN,STOPER,TJ,VIN,XA,XB,XPA,XPB,YAB,YB,YPB COMMON/MAP/RK,RMOD2,XX,YY COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE COMMON/PTPARI/GPVMN,GPVMX,LCODE,NO,NPV COMMON/PTPARZ/CNSTI,H2,RW,SGMAZI,TOL COMMON/TIMER/DTIME,NFULL,NSTART,NSTOP,NTIME,TIME DIMENSION F(NCOUNT,*) VIRTUAL UVDATA(2,*),XY(2,*),XYDATA(2,*) VIRTUAL FTAG(2,*),XPNT(*),YPNT(*),RPNT(*) VIRTUAL Q(NCOUNT,*),UG1$(NCOUNT,*),VGlS(NCOUNT,*) VIRTUAL UGN(NCOUNT,*),UG2(NCOUNT,*),VGN(NCOUNT,*),VG2(NCOUNT,*) C*** C*********************************************************************** C*** IF(NPV.EQ.0) GO TO 300 C*** C*** C*** NON-GRID POINT DATA : C*** C*** IF(NPOS.EQ.1) THEN IF(NCASE.EQ.1) THEN WRITE(1,108) NTIME WRITE(2,108) NTIME DO 1 N - 1,NPV D IF(FTAG(1,N).EQ.1.0) THEN D WRITE(1,110) XPNT(N),YPNT(N),NTIME*DTIME,RPNT(N) D ELSE D WRITE(2,110) XPNT(N),YPNT(N),NTIME*DTIME,RPNT(N) 227 D END IF D l CONTINUE D END IF D END IF 108 FORMAT(' ’,T2,’NTIME I',T10,I5,/) D 110 FORMAT(' ',T2,1PE15.7,T19,1PE15.7,T37,1PE15.7,T54,1PE15.7) C*** C*** GRID POINT VELOCITY DATA : C*** c*** IF(NGRD.EQ.1.0) THEN C*** DO 2 NLOOP I 1,2 IF(NLOOP.EQ.1) THEN J1 I JMID J2 I JMID J3 I 1 13 I 10 ELSE J1 I JMAX J2 I JMIN J3 I -2 I3 I 2 END IF C*** DO 131 J I J1,J2,J3 DO 131 I I IMIN,IMAX,I3 K I I + (J-JMIN)*JMAX IF(NLOOP.EQ.1) THEN IF(XDIR.EQ.0.0) K I IMID + (I-JMIN)*JMAX IF(TANG.EQ.1.0) THEN IF(SNGL(XYDATA(1,K)).EQ.0.0.AND.SNGL(XYDATA(2,K)).EQ.0.0) THEN T I ATAN2(SNGL(XYDATA(2,K))ySGMAO) ELSE T I ATAN2(SNGL(XYDATA(2,K))ISNGL(XYDATA(1,K))) END IF VTC I ABS(VGlS(I,J))*COS(T)-ABS(UG1$(I,J))*SIN(T) VRC I IABS(VG1$(I,J))*SIN(T)-ABS(UG1$(I,J))*COS(T) c*** ELSE C*** VTC I UGIS(I,J) VRC I VGlS(I,J) END IF C*** IF(I.EQ.IMIN.AND.J.EQ.J1) THEN IF(NCASE.EQ.1) WRITE(3,3) NTIME IF(NCASE.EQ.2) WRITE(7,3) NTIME END IF C*** IF(XDIR.EQ.1.0) THEN IF(NCASE.EQ.1) THEN WRITE(3,4) SNGL(XTDATA(1,K)),SNGL(XXDATA(2,K)),P(I,J), 228 *VTC,VRC,Q(I,J) ELSE WRITE(7,4) SNGL(UVDATA(1,K)/RK),SNGL(UVDATA(2,K)/RK), *F(I,J).VG2(I,J).UGZ(I,J),Q(I,J) END IF C*** ELSE C*** IF(NCASE.EQ.1) THEN WRITE(3,4) SNGL SQUARE (u,v) ** C*** ** C*********************************************************************** C*** DOUBLE PRECISION G1,GZ,P,SN,SSl,UVDATA,XYDATA DOUBLE PRECISION RK,RMOD2,X,XX,YY COMMON/MAP/RK,RMOD2,XX,YY COMMON/PTPARl/GPVMN,GPVMX,LCODE,NO,NPV DIMENSION X(2) VIRTUAL UPT(*),VPT(*),XPT(*),YPT(*) C*** C*********************************************************************** C*** IF(IFLG.EQ.1) NP I NPV IF(IFLG.EQ.2) NP I NL DO 1 I I 1,NP XX I DBLE(XPT(I)) YY I DBLE(YPT(I)) CALL NASZ(G1,GZ,SSl,UVDATA,X,XYDATA) CALL INTPLU(MCOUNT,P,UPT(I),SN,X(2)) CALL INTPLU(MCOUNT,P,VPT(I),SN,X(1)) 1 CONTINUE C*** RETURN END C*** C*** Ct********************************************************************** VSUBROUTINE NASZ(G1,G2,SSl,UVDATA,X,XYDATA) C*********************************************************************** C*** C*** THIS SUBROUTINE PREPARES INPUT PARAMETERS FOR SUBROUTINE NASNWR ** 235 C*** *** C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION (A-H,O-Z) EXTERNAL F INTEGER N, MAXN, MAXIT, IRSTAT PARAMETER (MAXNI42) LOGICAL AGIVEN DIMENSION X(2), R(2), A(MAXN,MAXN), WK(MAXN,MAXN+3) VIRTUAL SSl(2,*),UVDATA(2,*),XYDATA(2,*) COMMON/MAP/RK,RMOD2,XX,YY COMMON/NSIZE/IMAX,IMID,IMIN,ISIZE,JMAX,JMID,JMIN,JSIZE C*** C*********************************************************************** C*** C*** Determine the sign of S and $1 : C*** C*** AGIVEN I .FALSE. IF(XX.GE.0.0DO) THEN IF(YY.GE.0.0D0) THEN X(l) I G1 X(2) I 62 ELSE X(l) I -61 X(2) I GZ END IF ELSE IF(YY.GE.0.0D0) THEN X(l) I 61 X(2) I IGZ ELSE X(l) I -G1 X(2) I -G2 END IF END IF C*** TOL I 1.0D-6 MAXIT I 400 N I 2 C*** C*** C*** Call the subroutine that would solve C*** F(x)-0 using a quasi Newton method : c*** c*** CALL NASNWR(F,X,N,R,A,MAXN,TOL,WK,MAXIT,AGIVEN,IRSTAT) C*** C*** C*** Limiting cases : C*** C*** IF(XX.EQ.0.0DO) X(2) I 0.0D0 IF(YY.EQ. DO 22 J I DO 22 I I K - 236 0.0D0) X(1) I 0.0D0 JMIN,JMAX IMIN,IMAX I + (J-JMIN)*JMAX IF(XX.EQ.XYDATA(1,K).AND.YY.EQ.XYDATA(2,K)) THEN IF(DABS(UVDATA(1,K)).EQ.DABS(UVDATA(2,K))) THEN X(1) I X(2) I END IF END IF 22 CONTINUE C*** RETURN END C*** C*** SSl(1,K) SSl(2,K) C***************t******************************************************* SUBROUTINE NASNWR(F,X,N,R,A,MAXN,TOL,WK,MAXIT,AGIVEN,IRSTAT) C*********************************************************************** c*** C*** C*** c*** THIS SUBROUTINE SOLVES A SYSTEM OF NONLINEAR ALGEBRAIC EQUATIONS BASED ON THE NEWTON-RAPHSON METHOD ** ** ** ** C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION (A-H,O-Z) EXTERNAL INTEGER LOGICAL DIMENSION INTRINSIC EXTERNAL EXTERNAL INTEGER INTEGER INTEGER INTEGER LOGICAL PARAMETER PARAMETER C*** F N, MAXN, AGIVEN X(1), DABS, NASABX, MOVEDD, MAXJME, IT, NPLUSl, J77, J88, EVALJM, XFOUND (MAXJMEIIO) (MAXJMSIZ) R(l), DMINl, NASJME, VIPASD MAXJMS IT1, NPLUSZ, MAXIT, IRSTAT A(MAXN,1), WK(MAXN,1) DMAXl, DSQRT NASFNE JMS, J10 NPLUS3 J99, IDUMMY C*********************************************************************** C*** NPLUSI I N+1 NPLUSZ I NPLU51+1 NPLUS3 I NPLU82+1 EVALJM I .NOT.AGIVEN C*** IT I 0 ITl I 0 JMS I 0 DFRAC I .5000 C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** C*** c*** C*** C*** C*** C*** C*** 188 237 Investigate the initial guess : XNORMZ I 0.D00 CALL VIPASD(X,1,X,1,N,XNORM2) IF (XNORM2.EQ.0.D00) THEN IRSTAT I -101 RETURN ENDIF CALL NASFNE(F,X,N,R,RNORM) CALL MOVEDD(N,R,WK(1,NPLU51)) Reevaluate the function Jacobi matrix if necessary : CONTINUE IF (EVALJM) THEN CALL NASJME(F,X,N,R,A,MAXN) ENDIF Iterate with an estimate of the function Jacobi matrix CONTINUE DO 399 J99I1,MAXJME XFOUND I RNORM.LT.TOL IF (XFOUND) THEN The residual norm is within given tolerance. All done : IRSTAT I 0 RETURN ELSE Next iteration is necessary : TNORM I RNORM Compute the nominal iteration step : DO 188 J88I1,N WK(J88,NPLUS3) I -R(J88) WK(J88,NPLU82) I WK(J88,NPLUSI) CALL MOVEDD(N,A(1,J88),WK(1,J88)) CONTINUE CALL NASABX(WK,MAXN,N,WK(1,NPLU83),WK(1,NPLUSl),R,IDUMMY) C*** 238 IF (IDUMMY.EQ.0) THEN JMS - JMS+1 IF (JMS.GT.MAXJMS) THEN CALL NASFNE(F,X,N,R,RNORM) CALL NASJME(F,X,N,R,A,MAXN) IRSTAT - -10 RETURN ELSE CALL MOVEDD(N,WK(1,NPLUSZ),WK(1,NPLUSI)) ENDIF ELSE JMS I 0 ENDIF DNORM2 I 0.D00 CALL VIPASD(WK(1,NPLU81),1,WK(1,NPLUSI),1,N,DNORM2) IF (DNORM2.EQ.0D00) THEN IRSTAT - -102 RETURN ENDIF CORR - DMINl(DFRAC/DSQRT(DNORM2/XNORM2),1.D00) X8 I CORR J10 - o C*** C*** C*** Evaluate the function for new arguments : C*** 200 288 C*** C*** CONTINUE DO 288 J88I1,N WK(J88,NPLUSl) I CORR*WK(J88,NPLUSl) WK(J88,1) I X(J88)+WK(J88,NPLU31) CONTINUE CALL NASFNE (F, WK, N, R, RNORM) IF(RNORM.GE.TNORM) GOTO 4 IT1 I IT1+1 IF(IT1.EQ.1) GOTO 5 IT1 I 0 DFRAC I DMIN1(DFRAC+.2D00,2.D00) GOTO 5 CORR I .ZDOO X8 I X8*.2D00 IT1 I 0 DFRAC I DMAX1(DFRAC-.05D00,.1D00) J10 I J10+1 IF(J10.LE.2) GOTO 200 X8 I X8*5.D00 CONTINUE C*** C*** Extrapolate the function Jacobi matrix : C*** TEMPl - 1.D00/(DNORM2*X8**2) DNORM2 - 1.D00-X8 239 D0 388 J88I1,N X(J88) I WK(J88,1) TEMP2 I (R(J88)+WK(J88,NPLUS3)*DNORMZ)*TEMPl DO 377 J77I1,N A(J88,J77) I A(J88,J77)+TEMP2*WK(J77,NPLU31) 377 CONTINUE 388 CONTINUE ENDIF 399 CONTINUE XFOUND I RNORM.LT.TOL IF (XFOUND) THEN C*** C*** C*** The residual norm is within given tolerance, fortunatelly : C*** C*** IRSTAT I 0 RETURN ELSE C*** C*** C*** Next iteration is necessary : C*** C*** IT - IT+MAXJME IF (IT.GE.MAXIT) THEN C*** C*** C*** Iteration limit reached : C*** C*** IRSTAT I I1 RETURN ELSEIF (RNORM.LE..5D00*TNORM) THEN C*** C*** C*** The function Jacobi matrix need not be updated : C*** C*** GO TO 2 ELSE c*** C*** C*** The function Jacobi matrix must be updated : C*** C*** EVALJM I .TRUE. GO TO 1 ENDIF ENDIF c*** RETURN END c*** 240 C*** C*********************************************************************** SUBROUTINE NASJME(F,X,N,R,A,MAXN) C*********************************************************************** c*** ** C*** THIS SUBROUTINE EVALUATES THE FUNCTION JACOBI MATRIX ** C*** ** c***************************************************************************** C*** EXTERNAL F INTEGER N, MAXN DOUBLE PRECISION X, R, A DIMENSION X(1), R(1), A(MAXN,1) INTRINSIC DABS EXTERNAL NASEPS DOUBLE PRECISION DABS, NASEPS DOUBLE PRECISION DIFFR, XINCR, XPERT, XSAVE INTEGER ZEROS INTEGER J88, J99 LOGICAL FIRST, MODIFY SAVE DIFFR, FIRST DATA FIRST/.TRUE./ C*** c*********************************************************************** C*** C*** Deternime the nominal difference-to-argument ratio : C*** C*** IF (FIRST) THEN DIFFR I DSQRT(NASEPS(1)) FIRST I .FALSE. ENDIF C*** C*** C*** Compute the Jacobi matrix (by columns) : C*** C*** DO 199 J99-1,N MODIFY I .FALSE. 100 CONTINUE c*** C*** C*** Compute perturbed argument : C*** C*** XSAVE I X(J99) IF (XSAVE.EQ..0D00.0R.MODIFY) THEN XINCR I DIFFR ELSE XINCR I XSAVE*DIFFR ENDIF XPERT I XSAVE+XINCR X(J99) I XPERT c*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** 188 C*** c*** C*** C*** C*** c*** 241 Evaluate the function residuals CALL F(X,N,A(1,J99)) Restore the argument : X(J99) I XSAVE Compute an estimated partial derivative : ZEROS I 0 DO 188 J88I1,N A(J88,J99) I (A(J88,J99)-R(J88))/XINCR IF (A(J88,J99).EQ..0D00) THEN ZEROS I ZEROS+1 ENDIF CONTINUE Test for possibly misleading results. Deside whether to modify the finite difference and recompute the column : IF (ZEROS.EQ.N.AND..NOT.MODIFY) THEN MODIFY I .TRUE. GO TO 100 ENDIF 199 CONTINUE C*** C*** C*** RETURN END C*********************************************************************** SUBROUTINE NASFNE(F,X,N,R,NORM) Ci********************************************************************** C*** ** C*** THIS SUBROUTINE EVALUATES THE FUNCTION RESIDUALS AND THEIR NORM ** C*** ** C*********************************************************************** c*** C*** EXTERNAL F INTEGER N DOUBLE PRECISION X, DIMENSION X(1), INTRINSIC DSQRT DOUBLE PRECISION DSQRT R: R(l) NORM 242 C*********************************************************************** C*** C*** Evaluate residuals first C*** C*** CALL F(X,N,R) C*** C*** C*** Evaluate the norm of R : C*** C*** CALL NORMSD(R,1,N,NORM) C*** RETURN END C*** C*** C*********************************************************************** DOUBLE PRECISION FUNCTION NASEPS(IDUMMY) C*********************************************************************** C*** ** C*** THIS FUNCTION COMPUTES THE HOST ARITHMETIC UNIT ACCURACY ** C*** ** c*********************************************************************** C*** INTEGER IDUMMY DOUBLE PRECISION EPS, BUFFER INTEGER NR INTEGER J99 LOGICAL FIRST PARAMETER (NRI32) DIMENSION BUFFER(NR) SAVE EPS, FIRST DATA EPS/l.D00/,FIRST/.TRUE./ C*** C*********************************************************************** C*** IF (FIRST) THEN 10 CONTINUE EPS - EPS/2.D00 EUFFER(1) - 1.DOO+EPS Do 199 J99-1,NR—1 BUFFER(J99+1) - BUFFER(J99) 199 CONTINUE DO 299 J99-1,NR IF (BUFFER(J99).EQ.1.DOO) THEN FIRST - .FALSE. EPS - EPS*2.D00 NASEPS - EPS RETURN ENDIF 299 CONTINUE Go To 10 ELSE a—v 243 NASEPS - EPS ENDIF C*** RETURN END C*** C*** c*********************************************************************** SUBROUTINE NASABX(A,NA,N,B,X,MP,IDET) C*********************************************************************** C*** ** C*** THIS SUBROUTINE SOLVES A*XIB BY LU DECOMPOSITION ** C*** AND BACKWARDS SUBSTITUTION ** C*** ** C*********************************************************************** C*** INTEGER NA, N, IDET, MP DOUBLE PRECISION A, B, x DIMENSION A(NA,1),B(1), X(1), MP(1) C*** C*********************************************************************** C*** CALL DECOM(A,NA,N,X,MP,IDET) IF (IDET.EQ.O) RETURN CALL SUBST(A,NA,N,B,x,MP) RETURN END C*** C*** C*********************************************************************** SUBROUTINE DECOM(A,NA,N,VS,MP,IDET) C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(NA,1),VS(1),MP(1) DOUBLE PRECISION NASEPS C*** C*********************************************************************** C*** EPS - NASEPS(1)*2**MIN(4,N/2) DO 100 J00I1,N CALL NORMSD(A(J00,1),NA,N,TEMP) IF (TEMP.EQ..0D0) GO TO 4 VS(J00) I .lDl/TEMP MP(J00) I J00 100 CONTINUE IDET - 1 DO 200 J00I1,N IS I J00 X I .000 I00 I JOO-l Do 110 J10-J00,N TEMP I -A(J10,J00) IF (Ioo.GT.0) CALL VIPASD(A(J10,1),NA,A(1,J00),1,Ioo,TEMP) 244 A(J10,J00) I ITEMP TEMP I DABS(TEMP*VS (J10)) IF (X.GE.TEMP) GO TO 110 X I TEMP IS I J10 110 CONTINUE IF (IS.EQ.J00) GO TO 1 IDET I -IDET CALL EXCHSD(A(IS,1),NA,A(J00,1),NA,N) TEMP I VS(J00) VS(J00) I VS(IS) VS(IS) I TEMP JS I MP(J00) MP(J00) I MP(IS) MP(IS) I JS 1 CONTINUE IF (X.LE.EPS) GO TO 4 X I -.1D1/A(J00,J00) K00 I J00+1 IF (K00.GT.N) GO TO 200 D0 120 J20IK00,N TEMP I -A(J00,J20) IF (IOO.GT.0) CALL VIPASD(A(J00,1),NA,A(1,J20),1,I00,TEMP) A(J00,J20) I X*TEMP 120 CONTINUE 200 CONTINUE RETURN 4 IDET I 0 RETURN END C*** C*** C*********************************************************************** SUBROUTINE SUBST(A,NA,N,B,VS,MP) c*********************************************************************** C*** IMPLICIT DOUBLE PRECISION(A-H,OIZ) DIMENSION A(NA,1),B(1),VS(1),MP(1) C*** C*********************************************************************** c*** DO 100 J00I1,N VS(J00) I B(MP(J00)) 100 CONTINUE VS(l) I IVS(1)/A(1:1) IF (N.EQ.1) GO TO 1 DO 200 J00I2,N TEMP I VS(J00) CALL VIPASD(A(J00,1),NA,VS(1),1,J00-1,TEMP) VS(J00) I -TEMP/A(J00,J00) 200 CONTINUE 1 M00 I N DO 300 J00I1,N 100 I J00I1 245 TEMP I VS(M00) IF (I00.GT.0) CALL VIPASD(A(M00,M00+1),NA,VS(M00+1),1,I00,TEMP) VS(M00) I ITEMP M00 I MOO-1 300 CONTINUE RETURN END C*** C*** C*********************************************************************** SUBROUTINE MOVEDD(N,X,Y) C*********************************************************************** c*** ** C*** THIS SUBROUTINE MOVES ELEMENTS OF A VECTOR TO ANOTHER VECTOR ** C*** ** C*********************************************************************** C*** INTEGER N DOUBLE PRECISION X, Y DIMENSION X(1), Y(l) INTEGER J99 C*** C*********************************************************************** C*** DO 199 J99-1,N Y(J99) I X(J99) 199 CONTINUE C*** RETURN END c*** c*** C*********************************************************************** SUBROUTINE NORMSD(X,INCX,N,NORM) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES THE EUCLIDEAN NORM OF A VECTOR ** c*** ** C*********************************************************************** c*** INTEGER INCX, N DOUBLE PRECISION X, NORM DIMENSION X(1) INTEGER J99 cr** C*********************************************************************** C*** NORM I .0000 DO 199 J99I1,INCX*N,INCX NORM I NORM+X(J99)**2 199 CONTINUE NORM I DSQRT(NORM) C*** RETURN 246 END C*** C*** C*********************************************************************** SUBROUTINE EXCHSD(X,INCX,Y,INCY,N) c*********************************************************************** C*** ** C*** THIS SUBROUTINE SWAPS ELEMENTS OF TWO VECTORS ** C*** ** C*********************************************************************** C*** INTEGER INCX, INCY, N DOUBLE PRECISION X, Y DIMENSION X(1), Y(l) INTEGER I99, J99 DOUBLE PRECISION TEMP C*** C*********************************************************************** C*** I99 I 1 DO 199 J99-1,INCX*N,INCX TEMP I X(J99) X(J99) I Y(I99) X(199) I TEMP I99 I I99+INCY 199 CONTINUE C*** RETURN END C*** C*** C*********************************************************************** SUBROUTINE VIPASD(X,INCX,Y,INCY,N,VIP) c*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES A VECTOR INNER PRODUCT ** C*** AND ADDS IT TO A GIVEN VALUE ** C*** ** C*********************************************************************** C*** INTEGER INCX, INCY, N DOUBLE PRECISION X, Y, VIP DIMENSION X(1), Y(1) INTEGER I99, J99 c*** c*********************************************************************** C*** I99 I 1 DO 199 J99I1,INCX*N,INCX VIPI VIP+X(J99)*Y(I99) I99 I I99+INCY 199 CONTINUE C*** RETURN 247 END c*** C*** C*********************************************************************** SUBROUTINE F(X,N,R) Ci*********************************************************************t C*** ** C*** THIS SUBROUTINE PROVIDES THE SYSTEM OF NONLINEAR ** C*** ALGEBRAIC EQUATIONS FOR USE IN SUBROUTINE NASNWR ** C*** ** C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N DIMENSION X(1), R(l) COMMON/MAP/RK,RMOD2,XX,YY C*** C*********************************************************************** C*** C*** Initialize parameters : C*** C*** C*** C*** K parameter : C*** C*** RKZ I RMODZ RK4 I RK2*RK2 RK6 I RK4*RK2 RKB I RK6*RK2 C*** C*** C*** K prime parameter : C*** C*** RKPZ I RMODZ RKP4 I RKP2*RKP2 C*** C*** C*** S parameter : C*** C*** 32 I X(2)*X(2) $4 I 82*82 56 I 84*52 $8 I 86*82 $10 I 88*82 $12 I 510*82 C*** C*** C*** 51 parameter : C*** C*** 802 I X(1)*X(1) C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** C*** c*** C*** 248 814 I 502*802 $16 I 814*802 $18 I 516*802 $110 I 518*802 $112 I 3110*302 X,Y parameters : X2 I XX*XX Y2 I YY*YY Compute the 81 polynomials : Numerator : PNX2 I 1.-2.*(3.-2.*RK2)*SOZ+(15.-18.*RK2+4.*RK4)*814- & 4.*(S.-3.*RK2)*RKP2*Sl6+(15.-13.*RK2)*RKP2*818- & 6.*RKP4*5110+RKP4*8112 PNX4 I -((1.+RK2)-2.*(3.+2.*RK2*RKP2)*302+(15.+7.*RK2-22. & *RK4+4.*RK6)*Sl4-4.*(5.+7.*RK2-5.*RK4)*RKP2*516+ & (15.+22.*RK2-29.*RK4)*RKP2*818-2.*(3.+8.*RK2)* & RKP4*8110+(1.+3.*RK2)*RKP4*3112) PNXG I (l.-2.*(4.-RK2)*802+(25.-15.*RK2-4.*RK4)*Sl4-4.* & (lO.+RK2-2.*RK4)*RKP2*816+(35.-3.*RK2-20.*RK4)* & RKP2*818-2.*(8.+7.*RK2)*RKP4*5110+3.*(1.+RK2)* & RKP4*5112)*RK2 PNX8 I -(-2.*802+(11.-7.*RK2)*Sl4-4.*(6.-RK2)*RKP2*516+ & 2.*(13.-7.*RK2-2.*RK4)*RKP2*818-2.*(7.+2.*RK2)* & RKP4*8110+(3.+RK2)*RKP4*8112)*RK4 PNXlOI (Sl4-4.*RKP2*816+2.*(3.—2.*RK2)*RKP2*318-4.*RKP4 & *8110+RKP4*3112)*RK6 PNYO I $02-(4.-RK2)*814+3.*(2.-RK2)*816-(4.-3*RK2)*818+RKP2*8110 PNY2 I -2.*(2.*302-(9.-2.*RK2)*Sl4+(15.-7.*RK2)*816- & (11.-8.*RK2)*818+3.*RKP2*8110)*RK2 PNY4 I (2.*(1.+2.*RK2)*502-2.*(4.+11.*RK2-2.*RK4)*Sl4+ & (12.+43.*RK2-20.*RK4)*316-(8.+36.*RK2-29.*RK4)*318 & +(2;+13*RK2)*RKP2*5110)*RK2 PNY6 I -4.*($02-(S.+RK2)*Sl4+(3.-RK2)*(3.+2.*RK2)*816-(7.+ & 3.*RK2-5.*RK4)*818+(2.+3.*RK2)*RKP2*5110)*RK4 PNY8 I (SOZ-(4.+7.*RK2)*Sl4+(6.+23.*RK2-4.*RK4)*Sl6-(4.+25.*RK2 & I10.*RK4-4.*RK6)*818+(1.+10.*RK2+4.*RK4)*RKP2*3110)*RK4 PNYlOI I2.*(ISl4+(3.+RK2)*SlS-(3.+2.*RK2*RKPZ)*SlB+ & (1.+2.*RK2)*RKP2*8110)*RKG PNY12I ($16-(2.-RK2)*818+RKP2*5110)*RK8 Denominator : 249 C*** C*** PDO I 1.-4.*602+6.*Sl4-4.*Sl6+818 PD2 I I2.*(1.-(5.+RK2)*502+2.*(5.+RK2)*Sl4-10.*816+ & (5.-2.*RK2)*SlB-RKP2*3110) P04 I 1.-6.*(1.+RK2)*802+(15.+22.*RK2-RK4)*314- & 4.*(5.+7.*RK2-2.*RK4)*Sl6+3.*(5.+4.*RK2*RKP2)* & $18-2.*(3.+2.*RK2)*RKP2*8110+RKP4*8112 PD6 I -4.*(-SOZ+(5.+RK2)*Sl4-(10.+RK2-RK4)*Sl6+(lO.-3.* & RK2-2.*RK4)*818-5.*RKP2*5110+RKP4*5112)*RKZ PD8 I (6.*Sl4-4.*(6.-RK2)*Sl6+(36.-20.*RK2-RK4)*318- & 4.*(6.—RK2)*RKP2*5110+6.*RKP4*8112)*RK4 PDlO I -2.*(I2.*316+3.*(2.-RK2)*818-(6.-7.*RK2+RK4)*SllO & +2.*RKP4*3112)*RK6 P012 I (818-2.*RKP2*3110+RKP4*3112)*RKB C*** C*** C*** Compute the coefficients of the S polynomial C*** C*** PXO I -PDO*X2 PXZ I PNXZ-PD2*X2 PX4 = PNX4-PD4*X2 PX6 I PNX6-PD6*X2 PX8 I PNXB-PD8*X2 PX1O I PNXlO-PD10*X2 PX12 I -PD12*X2 C*** PYO I PNYO-PDO*Y2 PYZ I PNYZ-PD2*Y2 PY4 I PNY4-PD4*Y2 PY6 I PNYG-PD6*Y2 PY8 I PNYB-PD8*Y2 PYlO I PNYlO-PD10*Y2 PY12 I PNYlZ-PD12*Y2 C*** C*** C*** Compute the S polynomials C*** C*** R(l) I PYO+PY2*SZ+PY4*S4+PY6*SG+PY8*88+PY10*SlO+PY12*SlZ R(Z) I PXO+PX2*SZ+PX4*S4+PX6*$6+PX8*38+PX10*810+PX12*812 C*** RETURN END C*** C*** C*********************************************************************** SUBROUTINE INTPLU(MCOUNT,P,P0,SN,SO) C*********************************************************************** C*** ** C*** THIS SUBROUTINE PERFORMS LINEAR INTERPOLATION FOR U OR V ** C*** ** C*********************************************************************** 250 C*** DOUBLE PRECISION P,PD,SO,SN,SNQ VIRTUAL P(*),SN(*) C*** C*********************************************************************** C*** C*** Check the bounds on $0 : S/Sl(1) < SO < S/Sl(NPOINT) : c*** C*** IF(SO.LT.SN(1)) So - SN(I) IF(SO.GT.SN(MCOUNT)) so - SN(MCOUNT) C*** C*** C*** Determine the relative location of 80 C*** and interpolate for the value of U or V C*** C*** DO 1 I I 1,MCOUNT IF(SO.GT.SN(I)) GO TO 1 SNQ I (SO-SN(I))/(SN(I-1)-SN(I)) PD I P(I-1)IP(I) P0 I SNGL(SNQ*PD + P(I)) GO TO 2 1 CONTINUE 2 RETURN END C*** C*** C***********************************fi*********************************** SUBROUTINE CNDNSN( V0, CMO, CN1, DN1, SNl) C*********************************************************************** C*** ** C*** THIS SUBROUTINE EVALUATES JACOBIAN ELLIPTIC FUNCTIONS ** C*** USING A DESCENDING LANDEN TRANSFORMATION WITH THE REDUCED ** C*** FUNCTIONS APPROXIMATED BY CIRCULAR (TRIG) FUNCTIONS ** C***(See page 573, "Handbook of Mathematical Functions" by Abramowitz)** c*** ** C*********************************************************************** c*** IMPLICIT DOUBLE PRECISION ( A-H, O-z ) C*** c*********************************************************************** C*** C*** Perform three applications of reduction to reduce the parameter c*** C*** CMl I 1.000 - CMO RMU1 - ( ( 1.0Do-DSQRT(CM1) ) / ( 1.0D0+DSQRT(CM1) ) )**2 v1 - v0 / ( 1.0D0+DSQRT( RMU1 ) ) C*** CM2 - 1.0D0 - RMUl _ . RMU2 - ( ( 1.0D0-DSQRT(CM2) ) / ( 1.0D0+DSQRT(CM2) ) )**2 v2 - v1 / ( 1.0D0+DSQRT( RMU2 ) ) C*** 251 CM3 - 1.0D0 - RMU2 RMU3 - ( ( 1.0Do-DSQRT(CM3) ) / ( 1.0D0+DSQRT(CM3) ) )**2 v3 a v2 / ( 1.0D0+DSQRT( RMU3 ) ) C*** DN3 - DN(V3,RMU3)**2 + DSQRT(RMU3) - 1.0D0 DN3 - DN3/( 1.0Do + DSQRT(RMU3) - DN(V3,RMU3)**2 ) 8N3 - (1.0D0+DSQRT(RMU3))*SN(V3,RMU3) 8N3 - SN3/( 1.0D0 + DSQRT(RMUB)*SN(V3,RMU3)**2 ) CN3 - CN(V3,RMU3)*DN(V3,RMU3) CN3 - CN3/( 1.0Do + DSQRT(RMU3)*SN(V3,RMU3)**2 ) C*** DN2 - DN3**2 + DSQRT(RMUZ) - 1.0D0 DNz - DN2/( 1.0Do + DSQRT(RMU2) - DN3**2 ) SN2 - (1.0D0+DSQRT(RMU2))*SN3 SN2 - SN2/( 1.0D0 + DSQRT(RMU2)*SN3**2 ) CN2 - CN3*DN3 CN2 - CN2/( 1.0Do + DSQRT(RMU2)*SN3**2 ) C*** DN1 - DN2**2 + DSQRT(RMUl) - 1.0Do DN1 - DN1/( 1.0D0 + DSQRT(RMUl) - DN2**2 ) SN1 - (1.0D0+DSQRT(RMU1))*SN2 SN1 - SNl/( 1.0Do + DSQRT(RMU1)*SN2**2 ) CN1 - CN2*DN2 CN1 - CN1/( 1.0D0 + DSQRT(RMU1)*SN2**2 ) C*** RETURN END C*** C*** C*********************************************************************** FUNCTION CN( X, RM ) C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION ( A-H, OIZ ) CN I DCOS(X) + 0.25D0*RM*(X-DSIN(X)*DCOS(X))*DSIN(X) RETURN END C*** C*** C*********************************************************************** FUNCTION DN( X, RM ) C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DN I 1.0D0 - 0.SDO*RM*DSIN(X)**2 RETURN END C*** C*** ctttttttttt*********************************************t*************** FUNCTION SN( X, RM ) C***ttt*ttt*******************t***************************************** C*** IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) 252 SN I DSIN(X) - 0.25D0*RM*(X-DSIN(X)*DCOS(X))*DCOS(X) RETURN END C*** C*** C*********************************************************************** FUNCTION FUNCTN( DX, CMO ) C*********************************************************************** C*** ** C*** THIS FUNCTION EVALUATES THE VALUE OF THE INTEGRAND ** C*** [1/sqrt(1-k**23in**2(phi)] OR [l/sqrt(1—k'**23in**2(phi)] ** C*** IN THE COMPLETE ELLIPTIC INTEGRAL (K) OR ITS COMPLEMENTARY (K') ** C*** RESPECTIVELY ** C*** ** c*********************************************************************** C*** IMPLICIT DOUBLE PRECISION (A-H, O-z ) C*** C*********************************************************************** C*** FUNCTN - ( 1.0D0 - CMO*DSIN(DX)**2 ) FUNCTN - DSQRT( FUNCTN ) FUNCTN - 1.0D0/FUNCTN RETURN END C*** C*** C*********************************************************************** FUNCTION GAUSS( A, B, M, FUNCTN, CMO ) C*********************************************************************** C*** IMPLICIT DOUBLE PRECISION ( A-H, O-z ) DIMENSION NPOINT(7), KEY(8), 2(24), WEIGHT(24) DATA NPOINT / 2, 3, 4, S, 6, 10, 15 / DATA KEY / 1, 2, 4, 6, 9, 12, 17, 25 / DATA 2 / 0.577350269, 0.0, 0.774596669, 0.339981044, 0.861136312, .0, 0.538469310, 0.906179846, 0.238619186, 0.661209387, .932469514, 0.148874339, 0.433395394, 0.679409568, .865063367, 0.973906529, 0.0, 0.201194094, 0.394151347, .570972173, 0.724417731, 0.848206583, 0.937273392, 0.987992518 / DATA WEIGHT / 1.0, 0.888888889, 0.555555556, 0.652145155, 0.347854845, 0.568888889, 0.478628671, 0.236926885, 0.467913935, 0.360761573, .171324493, 0.295524225, 0.269266719, 0.219086363, .149451349, 0.066671344, 0.202578242, 0.198431485, .186161000, 0.166269206, 0.139570678, 0.107159221, .070366047, 0.030753242 / m m m m m CDCDCDCD m b m b m C)C>C)O C*** C*********************************************************************** C*** C*** Find the subscripts of the first 2 and weight values C*** C*** DO 1 II1,7 253 IF( M.EQ.NPOINT(I) ) GO TO 2 1 CONTINUE C*** C*** INVALID M USED C*** GAUSS I 0.0D0 RETURN C*** C*** C*** Setup initial parameters : C*** C*** 2 JFIRST I KEY(I) JLAST I KEY(I+1) - 1 C I 0.5D0*(B-A) D I 0.5D0*(B+A) C*** C*** C*** Accumulate the sum in the M-POINT formula C*** C*** SUM I 0.0D0 DO 3 JIJFIRST, JLAST IF( Z(J).EQ.0.0D0 ) SUM I SUM + WEIGHT(J)*FUNCTN(D,CMO) IF( Z(J).NE.0.0D0 ) SUM I SUM + WEIGHT(J)*(FUNCTN(Z(J)*C+D,CMO) * + FUNCTN(-Z(J)*C+D,CMO)) 3 CONTINUE C*** C*** C*** Make the interval correction and the return C*** C*** GAUSS I SUM*C RETURN END C*** 254 C*** C*********************************************************************** c********* POISSON LIBRARY CONTAINS THE FOLLOWING SUBROUTINES : ** C****** ** C*** 1. COSGEN - Computes requiredcosine values in ascending order ** C*** 2. GENBUN - Reorders unknowns, reverses columns and returns ** C*** 3. HWSCRT - Does the Fast Fourier Transform [In Poisson/Lib] ** C*** 4. MERGE - Merges two ascending strings of nos. in array TCOS ** C*** S. POISD2 - Solves Pooisson’s eqn. for Dirichlet bndry conds. ** C*** 6. POISN2 - Solves Pooisson's eqn. for Neumann bndry conds. ** C*** 7. POISP2 - Solves Pooisson's eqn. for Periodic bndry conds. ** C*** 8. TRIX - Solves a system of linear equations where the ** C*** coefficient matrix is a rational fctn. given by ** C*** Tridiagonal (..., A(I), B(I), C(I), ...) ** C*** 9. TRIB - Solves three linear systems whose common coeffs. ** C*** matrix is a rational fctn.in the matrix given by ** C*** Tridiagonal (..., A(I), 8(1), C(I), ...) ** C*** ** C*********************************************************************** SUBROUTINE COSGEN (N,IJUMP,FNUM,FDEN,A) C*********************************************************************** C*** ** C*** THIS SUBROUTINE COMPUTES REQUIRED COSINE VALUES IN ASCENDING ** C*** ORDER. WHEN IJUMP .GT. 1 THE ROUTINE COMPUTES VALUES ** C*** ** C*** 2*COS(J*PI/L) , JI1,2,...,L AND J .NE. 0(MOD N/IJUMP+1) ** C*** WHERE L - IJUMP*(N/IJUMP+1). ** C*** WHEN IJUMP - 1 IT COMPUTES ** C*** 2*COS((J-FNUM)*PI/(N+FDEN)) , J=1, 2, ... ,N ** C*** WHERE ** C*** FNUM - 0.5, EDEN - 0.0, FOR REGULAR REDUCTION VALUES ** C*** FNUM - 0.0, EDEN - 1.0, FOR B-R AND C-R WHEN ISTAG - 1 ** C*** FNUM I 0.0, EDEN I 0.5, FOR B-R AND C-R WHEN ISTAG - 2 ** C*** FNUM - 0.5, EDEN - 0.5, FOR B-R AND C-R WHEN ISTAG - 2 ** C*** IN POISN2 ONLY. ** C*** ** C*********************************************************************** C*** DIMENSION A(l) C*** C*********************************************************************** C*** PI I 3.14159265358979 IE (N .EQ. 0) GO To 105 IE (IJUMP .E0. 1) GO TO 103 K3 - N/IJUMP+1 K4 I K3-l PIBYN - PI/ELOAT(N+IJUMP) Do 102 KI1,IJUMP K1 I (K-1)*K3 K5 - (H-1)*K4 DO 101 II1,K4 x - K1+I K2 - KS+I 255 A(K2) I I2.*COS(X*PIBYN) 101 CONTINUE 102 CONTINUE GO TO 105 103 CONTINUE NPl I N+1 Y I PI/(FLOAT(N)+FDEN) DO 104 II1,N X I FLOAT(NPl-I)IFNUM A(I) I 2.*COS(X*Y) 104 CONTINUE 105 CONTINUE RETURN END C*** C*** C*********************************************************************** SUBROUTINE GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR,W) c*********************************************************************** C*** DIMENSION Y(IDIMY,1) DIMENSION W(1) ,B(1) ,A(l) ,C(1) C*** C*********************************************************************** C*** IERROR I 0 IF (M .LE. 2) IERROR I 1 IF (N .LE. 2) IERROR I 2 IF (IDIMY .LT. M) IERROR I 3 IF (NPEROD.LT.0 .OR. NPEROD.GT.4) IERROR I 4 IF (MPEROD.LT.0 .OR. MPEROD.GT.1) IERROR I 5 IF (MPEROD .EQ. 1) GO TO 102 DO 101 II2,M IF (A(I) .NE. C(1)) GO TO 103 IF (C(I) .NE. C(1)) GO TO 103 IF (B(I) .NE. 8(1)) GO TO 103 101 CONTINUE GO TO 104 102 IF (A(l).NE.0. .OR. C(M).NE.0.) IERROR I 7 GO TO 104 103 IERROR I 6 104 IF (IERROR .NE. 0) RETURN MP1 I M+1 IWBA I MP1 IWBB I IWBA+M IWBC I IWBB+M IWBZ I IWBC+M IWB3 I IWBZ+M IWW1 I IWB3+M IWWZ I IWW1+M IWW3 I IWW2+M IWD I IWW3+M IWTCOS I IWD+M IWP I IWTCOS+4*N 256 D0 106 II1,M K I IWBA+II1 W(K) I -A(I) K I IWBC+II1 W(K) I -C(I) K I IWBB+II1 W(K) I 2.IB(I) DO 105 JI1,N Y(I,J) - 'Y‘IIJ) 105 CONTINUE 106 CONTINUE MP I MPEROD+1 NP I NPEROD+1 GO TO (114,107),MP 107 GO TO (108,109,110,111,123),NP 108 CALL POISP2 (M,N,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2), 1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS), 2 W(IWP)) GO To 112 109 CALL POISD2 (M,N,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWW1), 1 W(IWD),W(IWTCOS),W(IWP)) GO To 112 110 CALL POISN2 (M,N,1,2,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWBZ), 1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS), 2 W(IWP)) GO To 112 111 CALL POISN2 (M,N,1,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWBZ), 1 W(IWBB),W(IWW1),W(IWW2),W(Iww3),W(IWD),W(IWTCOS), 2 W(IWP)) 112 IPSTOR - W(IWWl) IREv - 2 IF (NPEROD .EQ. 4) GO TO 124 113 GO TO (127,133),MP 114 CONTINUE C*** c*** C*** Reorder unknowns when MPIO C*** C*** MH I (M+1)/2 MHMl I MH-l MODD I 1 IF (MH*2 .EQ. M) MODD I 2 DO 119 JI1,N DO 115 II1,MHM1 MHPI I MH+I MHMI I MH-I W(I) I Y(MHMI,J)-Y(MHPI,J) W(MHPI) I Y(MHMI,J)+Y(MHPI,J) 11$ CONTINUE W(MH) I 2.*Y(MH,J) GO TO (117,116),MODD 116 W(M) I 2.*Y(M,J) 117 CONTINUE 118 119 120 121 122 C*** 257 D0 118 II1,M Y(I,J) I W(I) CONTINUE CONTINUE K I IWBC+MHM1~1 I I IWBA+MHM1 W(K) I 0. W(I) I 0. W(K+1) I 2.*W(K+1) GO TO (120,121),MODD CONTINUE K I IWBB+MHM1I1 W(K) I W(K) -W(I-1) W(IWBCIl) I W(IWBC-1)+W(IWBB-l) GO TO 122 W(IWBB-l) I W(K+1) CONTINUE GO TO 107 C*** C*** Reverse columns when NPERODI4 C*** C*** 123 124 125 126 127 128 129 130 131 132 133 C*** C*** IREV I 1 NBY2 I N/2 DO 126 JI1,NBY2 MSKIP I N+1-J DO 125 II1,M A1 I Y(I,J) Y(I,J) I Y(I,MSKIP) Y(I,MSKIP) I A1 CONTINUE CONTINUE GO TO (110,113),IREV CONTINUE DO 132 JI1,N DO 128 II1,MHM1 MHMI I MH-I MHPI I MH+I W(MHMI) I .5*(Y(MHPI,J)+Y(I,J)) W(MHPI) I .5*(Y(MHPI,J)-Y(I,J)) CONTINUE W(MH) I .S*Y(MH,J) GO TO (130,129),MODD W(M) I .5*Y(M,J) CONTINUE DO 131 II1,M Y(I,J) I W(I) CONTINUE CONTINUE CONTINUE C*** Return storage requirements for W array 258 C*** C*** W(l) - IPSTOR+IWP-1 RETURN END C*** C*** C*********************************************************************** SUBROUTINE HWSCRT(BDA,BDB,BDC,BDD,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C*************************************i*it****************************** C*** ** C*** W:A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR ** C*** WORK SPACE. W MAY REQUIRE UP TO 4*(N+1)+(13+INT(LOG2(N+1)))*(M+1)** C*** LOCATIONS. THE ACTUAL NUMBER OF LOCATIONS USED ID COMPUTED BY ** C*** HWSCRT AND IS RETURNED IN LOCATION W(l). ** C*** ** C*********************************t************************************* C*** COMMON/BNDARY/MBDCND,NBDCND COMMON/GRID1 /DU,DUINVT,DV,DVINVT,HDUINV,HDVINV COMMON/GRIDZ /ULNGTH,A,B,VLNGTH,C,D COMMON/NSIZE /MP1,IMID,IMIN,M,NP1,JMID,JMIN,N DIMENSION F(IDIMF,1), W(I) VIRTUAL BDA(*), BDB(*), BDC(*), BDD(*) C*** C*********************************************************************** C*** C*** Check for invalid parameters : C*** C*** IERROR - 0 IF (A .GE. B) IERROR - 1 IF (MBDCND.LT.0 .OR. MBDCND.GT.4) IERROR - 2 IE (C .GE. D) IERROR - 3 IF (N .LE. 3) IERROR - 4 IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR - 5 IE (IDIMF .LT. M+1) IERROR - 7 IE (M .LE. 3) IERROR - 8 IF (IERROR .NE. 0) RETURN NPEROD - NBDCND MPEROD - 0 IF (MBDCND .GT. 0) MPEROD - 1 TWDELx - 2./DU DELXSQ - 1./DU**2 TWDELY - 2./Dv DELYSQ I 1./DV**2 NP - NBDCND+1 MP - MBDCND+1 NSTART - 1 NSTOP - N NSKIP - 1 GO To (104,101,102,103,104),NP 101 NSTART - 2 . GO To 104 102 103 104 C*** C*** 259 NSTART I 2 NSTOP I NP1 NSKIP I 2 NUNK I NSTOP-NSTART+1 C*** C*** Enter boundary data for X-boundaries : C*** 105 106 107 108 109 110 111 112 113 114 115 116 117 C*** C*** MSTART I 1 MSTOP I M MSKIP I 1 GO TO (117,105,106,109,110),MP MSTART I 2 GO TO 107 MSTART I 2 MSTOP I MP1 MSKIP I 2 DO 108 JINSTART,NSTOP F(2,J) I F(2,J)-F(1,J)*DELXSQ CONTINUE GO TO 112 MSTOP I MP1 MSKIP I 2 DO 111 JINSTART,NSTOP F(1,J) I F(1,J)+BDA(J)*TWDELX CONTINUE GO TO (113,115),MSKIP DO 114 JINSTART,NSTOP F(M,J) I F(M,J)-F(MP1,J)*DELXSQ CONTINUE GO TO 117 DO 116 JINSTART,NSTOP F(MP1,J) I F(MP1,J)-BDB(J)*TWDELX CONTINUE MUNK I MSTOP-MSTART+1 C*** Enter boundary data for Y-boundaries : C*** C*** 118 119 120 121 122 123 124 GO TO (127,118,118,120,120),NP DO 119 IIMSTART,MSTOP F(I,2) I F(I,2)-F(I,1)*DELYSQ CONTINUE GO TO 122 D0 121 IIMSTART,MSTOP F(I,1) I F(I,1)+BDC(I)*TWDELY CONTINUE GO TO (123,125),NSKIP DO 124 IIMSTART,MSTOP F(I,N) I F(I,N)-F(I,NP1)*DELYSQ CONTINUE GO TO 127 260 125 D0 126 IIMSTART,MSTOP F(I,NPI) I F(I,NP1)-BDD(I)*TWDELY 126 CONTINUE C*** C*** C*** Multiply right side by DV**2 : C*** C*** 127 DELYSQ I DV*DV DO 129 IIMSTART,MSTOP DO 128 JINSTART,NSTOP F(I,J) I F(I,J)*DELYSQ 128 CONTINUE 129 CONTINUE C*** C*** C*** Define the A,B,C coefficients in W-array : C*** C*** ID2 I MUNK ID3 I ID2+MUNK ID4 I ID3+MUNK S I DELYSQ*DELXSQ 5T2 I 2.*S DO 130 II1,MUNK W(I) I S J I ID2+I W(J) I 'ST2+ELMBDA*DELYSQ J I ID3+I W(J) I S 130 CONTINUE IF (MP .EQ. 1) GO TO 131 W(I) I 0. W(ID4) I 0. 131 CONTINUE GO TO (135,135,132,133,134),MP 132 W(IDZ) I 8T2 GO TO 135 133 W(IDZ) I 5T2 134 W(ID3+1) I 3T2 135 CONTINUE PERTRB I 0. IF (ELMBDA) 144,137,136 136 IERROR I 6 GO TO 144 137 IF ((NBDCND.EQ.0 .OR. NBDCND.EQ.3) .AND. 1 (MBDCND.EQ.0 .OR. MBDCND.EQ.3)) GO TO 138 GO TO 144 C*** C*** C*** For singular problems must adjust data C*** to insure that a solution will exist C*** C*** 261 138 A1 I 1. A2 I 1. IF (NBDCND .EQ. 3) A2 = 2. IF (MBDCND .EQ. 3) A1 I 2. $1 I 0. MSPl I MSTART+1 MSTMl I MSTOP-1 NSPl I NSTART+1 NSTMl I NSTOP-1 DO 140 JINSP1,NSTM1 S I 0. DO 139 IIMSP1,MSTM1 S I S+F(I,J) 139 CONTINUE $1 I Sl+S*A1+F(MSTART,J)+F(MSTOP,J) 140 CONTINUE $1 I A2*Sl S I 0. DO 141 IIMSP1,MSTM1 S I S+F(I,NSTART)+F(I,NSTOP) 141 CONTINUE $1 I Sl+S*A1+F(MSTART,NSTART)+F(MSTART,NSTOP)+F(MSTOP,NSTART)+ 1 F(MSTOP,NSTOP) S I (2.+FLOAT(NUNK-2)*A2)*(2.+FLOAT(MUNK-2)*Al) PERTRB I 81/3 DO 143 JINSTART,NSTOP DO 142 IIMSTART,MSTOP F(I,J) I F(I,J)-PERTRB 142 CONTINUE 143 CONTINUE PERTRB I PERTRB/DELYSQ C*** C*** C*** Solve the equation : C*** c*** 144 CALL GENBUN (NPEROD,NUNK,MPEROD,MUNK,W(1),W(ID2+1),W(ID3+1), 1 IDIMF,F(MSTART,NSTART),IERR1,W(ID4+1)) W(l) I W(ID4+1)+3.*FLOAT(MUNK) C*** c*** C*** Fill in identical values when have periodic boundary conditions c*** C*** IF (NBDCND .NE. 0) GO TO 146 DO 145 IIMSTART,MSTOP F(IrNPl) I F(I,1) 145 CONTINUE 146 IF (MBDCND .NE. 0) GO TO 148 DO 147 JINSTART,NSTOP F(MplrJ) - F(I,J) 147 CONTINUE IF (NBDCND .EQ. 0) F(MP1,NP1) I F(I,NPI) 148 CONTINUE 262 RETURN END C*** C*** c*********************************************************************** SUBROUTINE MERGE (TCOS,I1,M1,I2,M2,I3) C*********************************************************************** C*** ** C*** THIS SUBROUTINE MERGES TWO ASCENDING STRINGS OF NUMBERS IN THE ** C*** ARRAY TCOS. THE FIRST STRING IS OF LENGTH M1 AND STARTS AT ** C*** TCOS(I1+1). THE SECOND STRING IS OF LENGTH M2 AND STARTS AT ** C*** TCOS(I2+1). THE MERGED STRING GOES INTO TCOS(I3+1). ** C*** ** C*********************************************************************** C*** DIMENSION TCOS(1) c*** C*********************************************************************** C*** J1 I 1 J2 I 1 J I I3 IF (M1 .EQ. 0) GO TO 107 IF (M2 .EQ. 0) GO TO 104 101 J I J+1 L I J1+Il X I TCOS(L) L I J2+I2 Y I TCOS(L) IF (X-Y) 102,102,103 102 TCOS(J) - x J1 - J1+1 IE (J1 .GT. M1) Go To 106 Go TO 101 103 TCOS(J) I Y J2 I J2+1 IF (J2 .LE. M2) GO TO 101 IF (J1 .GT. M1) GO TO 109 104 K - J-J1+1 DO 105 J-J1,M1 M - K+J L - J+Il TCOS(M) - TCOS(L) 105 CONTINUE GO TO 109 106 CONTINUE IF (J2 .GT. M2) GO TO 109 107 K - J-J2+1 DO 108 J-J2,M2 M - K+J L - J+I2 J TCOS(M) - TCOS(L) 108 CONTINUE 109 CONTINUE 263 RETURN END C*** C*** C*********************************************************************** SUBROUTINE POISD2 (MR,NR,ISTAG,BA,BB,BC,Q,IDIMQ,B,W,D,TCOS,P) C*********************************************************************** C*** ** C*** SUBROUTINE TO SOLVE POISSON'S EQUATION FOR DIRICHLET BOUNDARY ** C*** CONDITIONS. ** C*** ISTAG I 1 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A. ** C*** ISTAG I 2 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A+I. ** C*** ** C*********************************************************************** C*** DIMENSION Q(IDIMQ,1) ,BA(1) ,BB(1) ,BC(1) , 1 TCOS(1) ,B(1) ,D(1) ,W(1) , 2 P(l) C*** C*********************************************************************** C*** M I MR N I NR JSH I 0 PI I 1./FLOAT(ISTAG) IP I IN IPSTOR I 0 GO TO (101,102),ISTAG 101 KR I 0 IRREG I 1 IF (N .GT. 1) GO TO 106 TCOS(1) I 0. GO TO 103 102 KR I 1 JSTSAV I 1 IRREG I 2 IF (N .GT. 1) GO TO 106 TCOS(1) I I1. 103 DO 104 II1,M B(I) I Q(I,1) 104 CONTINUE CALL TRIX (1,0,M,EA,BB,EC,E,TCOS,D,W) DO 105 II1,M Q(I,1) I 3(1) 105 CONTINUE GO TO 183 106 LR I 0 DO 107 II1,M P(I) I 0. 107 CONTINUE NUN I N JST I 1 JSP I N C*** C*** 264 C*** C*** IRREG-1 when no irregularities have occurred, otherwise it is 2 C*** 108 C*** C*** L I 2*JST NODD I 2-2*((NUN+1)/2)+NUN C*** C*** NODD-1 when NUN is odd, otherwise it is 2 C*** 109 110 111 C*** C*** GO TO (110,109),NODD JSP I JSP-L GO TO 111 JSP I JSP-JST IF (IRREG .NE. 1) J5? I JSP-L CONTINUE C*** C*** Regular reduction : C*** 112 113 114 115 116 117 C*** CALL COSGEN (JST,1,0.5,0.0,TCOS) IF (L .GT. JSP) GO TO 118 DO 117 JIL,JSP,L JMl I J-JSH JPl I J+JSH JMZ I J-JST JP2 I J+JST JM3 I JMZ-JSH JP3 I JP2+JSH IF (JST .NE. 1) GO TO 113 DO 112 II1,M 3(1) I 20*Q(IIJ) Q(IIJ) - Q(IIJM2)+Q(IIJP2) CONTINUE GO TO 115 DO 114 II1,M T I Q(I,J)-Q(I,JMl)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2) B(I) I T+Q(I'J)-Q(IIJM3).Q(IIJP3) Q(IIJ) ' T CONTINUE CONTINUE CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W) ‘DO 116 II1,M Q(IIJ) - Q(IIJ)+B(I) CONTINUE CONTINUE c*** c*** C*** Reduction for last unknown : C*** 118 GO TO (119,135),NODD 119 C*** C*** 265 GO To (152,120),IRREC C*** c*** Odd number of unknowns : C*** 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 JSP - JSP+L J - JSP JMl - J-JSH JPl - J+JSH JMZ - J-JST JP2 - J+JST JM3 I JM2-JSH GO To (123,121),ISTAG CONTINUE IF (JST .NE. 1) GO To 123 D0 122 II1,M s(I) - Q(I,J) Q(I,J) - o. CONTINUE GO TO 130 GO To (124,126),NODDPR DO 125 II1,M IPl - IP+I 3(1) I .5*(Q(IpJM2)-Q(I:JM1)-Q(I'JM3))+P(IP1)+Q(I,J) CONTINUE GO To 128 DC 127 II1,M 8(1) I .5*(Q(I'JMZ)-Q(I:JM1)-Q(I,JM3))+Q(I,JP2)-Q(I,JP1)+Q(I,J) CONTINUE DO 129 II1,M Q(IIJ) I 05*(Q(IIJ)-Q(IIJM1).Q(IIJP1)) CONTINUE CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W) IP - IP+M IPSTOR - MAXO(IPSTOR,IP+M) DO 131 II1,M 1P1 - IP+I P(IPl) - Q(I,J)+B(I) 3(1) I Q(I,JP2)+P(IPI) CONTINUE IF (LR .NE. 0) GO TO 133 D0 132 I-1,JST KRPI - KR+I 'TCOS(KRPI) - TCOS(1) CONTINUE GO To 134 CONTINUE CALL COSGEN (LR,JSTSAV,0.,FI,TCOS(JST+1)) CALL MERGE (TCOS,0,JST,JST,LR,KR) CONTINUE CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS) CALL TRIX (KR,KR,M,BA,BB,BC,B,TCOS,D,W) DO 135 II1,M 135 C*** C*** 266 IPl - IP+I Q(I,J) - Q(I,JM2)+B(I)+P(IP1) CONTINUE LR - KR KR - KR+L GO To 152 C*** Even number of unknowns : C*** C*** 136 137 138 139 140 141 142 143 144 145 146 147 148 JSP I JSP+L J I JSP JMl J-JSH JPl J+JSH JMZ J-JST JPZ J+JST JM3 JM2-JSH GO TO (137,138),IRREG CONTINUE JSTSAV I JST IDEG I JST KR I L GO TO 139 CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS) CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1)) IDEG I KR KR I KR+JST IF (JST .NE. 1) GO TO 141 IRREG I 2 DO 140 II1,M 8(1) I Q(I,J) Q(I,J) I Q(I,JMZ) CONTINUE GO TO 150 D0 142 II1,M 8(1) I Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3)) CONTINUE GO TO (143,145),IRREG DO 144 II1,M Q(I,J) I Q(I,JM2)+.5*(Q(I,J)-Q(I,JMl)-Q(I,JP1)) CONTINUE IRREG I 2 GO TO 150 CONTINUE GO TO (146,148),NODDPR DO 147 II1,M IPl I IP+I Q(I,J) I Q(I'JM2)+P(IP1) CONTINUE IP I IP-M GO TO 150 DO 149 II1,M Q(I,J) I Q(I,JM2)+Q(I,J)-Q(I,JM1) 149 150 151 152 C*** C*** 267 CONTINUE CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W) DO 151 II1,M Q(I,J) I Q(IrJ)+B(I) CONTINUE NUN - NUN/2 NODDPR - NODD an - JST JST - 2*JST IF (NUN .GE. 2) GO TO 108 C*** C*** Start solution : C*** 153 154 155 156 157 158 159 160 161 162 163 164 c*** C*** J I JSP DO 153 II1,M 3(1) I Q(I,J) CONTINUE GO TO (154,155),IRREG CONTINUE CALL COSGEN (JST,1,0.5,0.0,TCOS) IDEG I JST GO TO 156 KR I LR+JST CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS) CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1)) IDEG I KR CONTINUE CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W) JMl I J-JSH JP1 I J+JSH GO TO (157,159),IRREG DO 158 II1,M Q(IrJ) I .S*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I) CONTINUE GO TO 164 GO TO (160,162),NODDPR DO 161 II1,M IPl I IP+I Q(I,J) I P(IP1)+B(I) CONTINUE IP I IP-M GO TO 164 D0 163 II1,M Q(I,J) I Q(I,J)-Q(I,JM1)+B(I) CONTINUE CONTINUE C*** Start back substitution : C*** C*** JST I JST/2 268 JSH I JST/2 NUN I 2*NUN IF (NUN .GT. N) GO TO 183 D0 182 JIJST,N,L JMl I J-JSH JP1 I J+JSH JM2 I J-JST JP2 I J+JST IF (J .GT. JST) GO TO 166 D0 165 II1,M 3(1) I Q(I,J)+Q(I'JP2) 165 CONTINUE GO TO 170 166 IF (JP2 .LE. N) GO TO 168 D0 167 II1,M 8(1) I Q(I,J)+Q(I'JM2) 167 CONTINUE IF (JST .LT. JSTSAV) IRREG I 1 GO TO (170,171),IRREG 168 D0 169 II1,M 8(1) I Q(LJ)+Q(I'JM2)+Q(I,JP2) 169 CONTINUE 170 CONTINUE CALL COSGEN (JST,1,0.5,0.0,TCOS) IDEG I JST JDEG I 0 GO TO 172 171 IF (J+L .GT. N) LR I LR-JST KR I JST+LR CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS) CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1)) IDEG I KR JDEG I LR 172 CONTINUE CALL TRIX (IDEG,JDEG,M,BA,BB,BC,B,TCOS,D,W) IF (JST .GT. 1) GO TO 174 DO 173 II1,M Q(I,J) I B(I) 173 CONTINUE GO TO 182 174 IF (JP2 .GT. N) GO TO 177 175 DO 176 II1,M Q(LJ) I .5*(Q(I,J)-Q(I,JM1)-Q(I,JPI))+B(I) 176 CONTINUE GO TO 182 177 GO TO (175,178),IRREG 178 IF (J+JSH .GT. N) GO TO 180 D0 179 II1,M IPl I IP+I Q(I,J) I B(I)+P(IP1) 179 CONTINUE IP I IP-M GO TO 182 180 DO 181 II1,M 269 Q(I,J) I B(I)+Q(IIJ)-Q(IIJM1) 181 CONTINUE 182 CONTINUE L - L/Z GO TO 164 183 CONTINUE C*** C*** C*** Return storage requirements for P vectors : C*** C*** W(l) I IPSTOR RETURN END C*** C*** C***************************************************************************** SUBROUTINE POISN2 (M,N,ISTAG,MIXBND,A,BB,C,Q,IDIMQ,B,B2,B3,W,W2, 1 W3,D,TCOS,P) C***************************************************************************** C*** ** C*** SUBROUTINE TO SOLVE POISSON’S EQUATION WITH NEUMANN BOUNDARY ** C*** CONDITIONS. ** C*** ** C*** ISTAG I 1 IF THE LAST DIAGONAL BLOCK IS A. ** C*** ISTAG I 2 IF THE LAST DIAGONAL BLOCK IS A-I. ** C*** MIXBNDI 1 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTH BOUNDARIES.** C*** MIXBNDI 2 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTTOM AND ** C*** DIRICHLET CONDITION AT TOP.(FOR THIS CASE, MUST HAVE ISTAG I 1.) ** C*** ** C*********************************************************************** C*** DIMENSION A(l) ,BB(1) ,C(1) ,Q(IDIMQ,1) , 1 3(1) 132(1) 133(1) :W(1) r 2 W2(1) ,W3(1) ,D(1) ,TCOS(1) , 3 K(4) :P(1) EQUIVALENCE (K(1),K1) :(K(2),K2) ,(K(3),K3) :(K(4),K4) C*** c*********************************************************************** c*** FISTAG I 3-ISTAG FNUM I 1./FLOAT(ISTAG) FDEN I 0.5*FLOAT(ISTAG-1) MR I M IP I IMR IPSTOR I 0 I2R I 1 JR I 2 NR I N NLAST I N KR I 1 LR I 0 GO TO (101,103),ISTAG 101 CONTINUE 102 103 104 105 106 107 108 C*** 270 D0 102 II1,MR Q(IIN) I .5*Q(I,N) CONTINUE GO TO (103,104),MIXBND IF (N .LE. 3) GO TO 155 CONTINUE JR I 2*IZR NROD I 1 IF ((NR/2)*2 .EQ. NR) NROD I 0 GO TO (105,106),MIXBND JSTART I 1 GO TO 107 JSTART I JR NROD I 1-NROD CONTINUE JSTOP I NLAST-JR IF (NROD .EQ. 0) JSTOP I JSTOP-IZR CALL COSGEN (I2R,1,0.5,0.0,TCOS) IZRBYZ I I2R/2 IF (JSTOP .GE. JSTART) GO TO 108 J I JR GO TO 116 CONTINUE C*** C*** Regular reduction : C*** C*** 109 110 111 112 113 DO 115 JIJSTART,JSTOP,JR JPI I J+I2RBY2 JP2 I J+I2R JP3 I JP2+IZRBY2 I JIIZRBYZ JM2 I J-I2R I JM2-I2RBY2 IF (J .NE. 1) GO TO 109 JMl I JP1 JM2 I JPZ JM3 I JP3 CONTINUE IF (IZR .NE. 1) GO TO 111 IF (J .80. 1) JM2 I JP2 DO 110 II1,MR 3(1) I 2-*Q(IoJ) Q(I'J) ' Q(Ilm)+Q(IIJP2) CONTINUE GO TO 113 CONTINUE DO 112 II1,MR FI I Q(I,J) Q(I,J) I Q(I,J)-Q(I,JMl)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2) 3(1) I FI+Q(I:J)“Q(I,JM3)-Q(I:JP3) CONTINUE CONTINUE 114 C*** c*** 271 CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W) DO 114 II1,MR Q(IuJ) I Q(I,J)+B(I) CONTINUE C*** C*** End of reduction for regular unknowns C*** 115 C*** C*** CONTINUE C*** Begin special reduction for last unknown : C*** C*** 116 C*** C*** J I JSTOP+JR NLAST I J JMl I J-IZRBYZ JM2 I J-IZR JM3 I JM2-IZRBY2 IF (NROD .EQ. 0) GO TO 128 C*** C*** Odd number of unknowns : C*** 117 118 119 120 121 122 123 124 125 126 IF (12R .NE. 1) GO To 118 D0 117 II1,MR 3(1) - FISTAG*Q(I,J) Q(I,J) I Q(I,JMZ) CONTINUE GO TO 126 D0 119 II1,MR 3(1) I Q(IIJ)+05*(Q(IIJM2)-Q(IIJM1)-Q(IIJM3)) CONTINUE IF (NRODPR .NE. 0) GO To 121 DC 120 II1,MR II - IP+I Q(I,J) I Q(I,JM2)+P(II) CONTINUE IP - IP-MR GO To 123 CONTINUE DO 122 II1,MR Q(IoJ) I Q(I,J)-Q(I'JM1)+Q(I,JM2) CONTINUE IF (LR .EQ. 0) GO TO 124 CALL COSGEN (LR,1,0.5,FDEN,TCOS(KR+1)) GO To 126 CONTINUE DO 125 II1,MR B(I) - FISTAG*B(I) CONTINUE CONTINUE 127 128 C*** C*** 272 CALL COSGEN (KR, 1 , 0 . 5, FDEN, TCOS) CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W) DO 127 II1,MR Q(IIJ) - Q(IIJ)+B(I) CONTINUE KR I KR+I2R GO TO 151 CONTINUE C*** C*** Even number of unknowns : C*** 129 130 131 132 133 134 135 136 137 138 139 JPl I J+IZRBY2 JP2 I J+I2R IF (I2R .NE. 1) GO TO 135 D0 129 IIl,MR B(I) I Q(I,J) CONTINUE CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) IP I 0 IPSTOR I MR GO TO (133,130),ISTAG DO 131 II1,MR P(I) I B(I) B(I) I B(I)+Q(I,N) CONTINUE TCOS(1) I 1. TCOS(Z) I 0. CALL TRIX (1,1,MR,A,BB,C,B,TCOS,D,W) DO 132 II1,MR Q(I,J) I Q(I,JM2)+P(I)+B(I) CONTINUE GO TO 150 CONTINUE DO 134 II1,MR P(I) I B(I) Q(I,J) I Q(I,JM2)+2.*Q(I,JP2)+3.*B(I) CONTINUE GO TO 150 CONTINUE DO 136 II1,MR B(I) I Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3)) CONTINUE IF (NRODPR .NE. 0) GO TO 138 D0 137 II1,MR II I IP+I B(I) I 3(I)+P(II) CONTINUE GO TO 140 CONTINUE DO 139 II1,MR B(I) I B(I)+Q(I,JP2)-Q(I,JP1) CONTINUE 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 C*** 273 CONTINUE CALL TRIX (12R,o,MR,A,EE,C,E,TCOS,D,W) IP - IP+MR IPSTOR - MAXO(IPSTOR,IP+MR) DO 141 II1,MR II - IP+I P(II) - B(I)+.5*(Q(I,J)-Q(I,JMI)-Q(I.JP1)) B(I) - P(II)+Q(I,JPZ) CONTINUE IE (LR .EQ. 0) GO To 142 CALL COSGEN (LR,1,0.5,FDEN,TCOS(IZR+1)) CALL MERGE (TCOS,0,12R,IZR,LR,KR) GO To 144 DO 143 I-1,I2R II - KR+I TCOS(II) - TCOS(I) CONTINUE CALL COSGEN (KR,1,0.5,FDEN,TCOS) IE (LR .NE. 0) GO To 145 GO To (146,145),I$TAG CONTINUE CALL TRIX (KR,KR,MR,A,33,C,3,TCOS,D,W) GO To 148 CONTINUE DO 147 II1,MR B(I) - FISTAG*3(I) CONTINUE CONTINUE DO 149 II1,MR II - IP+I Q(I,J) - Q(I,JM2)+P(II)+B(I) CONTINUE CONTINUE LR - KR KR - KR+JR CONTINUE GO To (152,153),MIXBND NR - (NLAST-1)/JR+1 E (NR .LE. 3) GO To 155 GO To 154 NR - NLAST/JR IE (NR .LE. 1) GO TO 192 I2R - JR NRODPR - NROD GO To 104 CONTINUE H c*** c*** Begin solution : C*** C*** J I 1+JR JMl I J-I2R JPl I J+I2R C*** C*** 274 JM2 I NLAST-12R IF (NR .EQ. 2) GO TO 184 IF (LR .NE. 0) GO TO 170 IF (N .NE. 3) GO TO 161 C*** C*** Case N-3 C*** 156 157 158 159 160 C*** GO TO (156,168),ISTAG CONTINUE DO 157 II1,MR B(I) I Q(IIZ) CONTINUE TCOS(1) I 0. CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 158 II1,MR Q(IIZ) I B(I) B(I) I 4.*B(I)+Q(I,1)+2.*Q(I,3) CONTINUE TCOS(1) I -2. TCOS(Z) I 2. I1 I 2 I2 I 0 CALL TRIX (I1,I2,MR,A,BB,C,B,TCOS,D,W) DO 159 II1,MR Q(IIZ) ' Q(II2)+B(I) B(I) I Q(I'1)+2.*Q(I:2) CONTINUE TCOS(1) I 0. CALL TRIX (I,0,MR,A,BB,C,B,TCOS,D,W) DO 160 II1,MR Q(Ipl) I B(I) CONTINUE JR I 1 IZR I 0 GO TO 194 C*** c*** c*** Case N - 2**P+1 C*** 161 162 163 164 CONTINUE GO TO (162,170),ISTAG CONTINUE DO 163 II1,MR 3(1) I Q(IIJ)+-5*Q(II1)_Q(IIJM1)+Q(IINLAST)-Q(IIJM2) CONTINUE CALL COSGEN (JR,1,0.5,0.0,TCOS) CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) DO 164 II1,MR Q(IIJ) ' 05*(Q(IIJ)-Q(IIJM1)-Q(IIJP1))+B(I) B(I) I Q(I,I)+2.*Q(I,NLAST)+4.*Q(I,J) CONTINUE 165 166 167 C*** C*** 275 JR2 I 2*JR CALL COSGEN (JR, 1, 0.0, 0.0,TCOS) DO 165 II1,JR 11 I JR+I 12 I JR+1II TCOS(II) I -TCOS(I2) CONTINUE CALL TRIX (JR2,0,MR,A,BB,C,B,TCOS,D,W) DO 166 II1,MR Q(I,J) I Q(I,J)+B(I) B(I) I Q(I,1)+2.*Q(I,J) CONTINUE CALL COSGEN (JR,1,0.5,0.0,TCOS) CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) DO 167 II1,MR Q(I,1) I .5*Q(I,1)-Q(I,JM1)+B(I) CONTINUE GO TO 194 C*** C*** Case of general N with NR I 3 C*** 168 169 170 171 172 173 174 175 176 177 DO 169 II1,MR B(I) I Q(IIZ) Q(IIZ) - 0. 32(I) I Q(I,3) 33(I) I Q(Irl) CONTINUE JR I 1 I2R I 0 J I 2 GO TO 177 CONTINUE DO 171 II1,MR B(I) I 05*Q(Ill)-Q(IIJM1)+Q(IIJ) CONTINUE IF (NROD .NE. 0) GO TO 173 DO 172 II1,MR II I IP+I B(I) I B(I)+P(II) CONTINUE GO TO 175 D0 174 II1,MR B(I) I B(I)+Q(I,NLAST)‘Q(I,JM2) CONTINUE CONTINUE DO 176 II1,MR T I .5*(Q(IIJ)-Q(IIJM1)“Q(I,JPl)) Q(I,J) - T 82(I) I Q(I,NLAST)+T B3(I) I Q(I,1)+2.*T CONTINUE CONTINUE 178 179 180 181 182 183 184 C*** K1 I K2 I 276 KR+2*JR-l KR+JR TCOS(K1+1) I I2. K4 I CALL K4 I CALL CALL K3 I CALL K4 I CALL CALL IF (LR .EQ. CALL CALL CALL K3 I K4 I CALL K1+3-ISTAG COSGEN (K2+ISTAG-2,1,0.0,ENUM,TCOS(K4)) K1+K2+1 COSGEN (JR-1,1,0.0,1.0,TCOS(K4)) MERGE (TCOS,K1,K2,K1+K2,JR-1,0) K1+K2+LR COSGEN (JR,1,0.5,0.0,TCOS(K3+1)) K3+JR+1 COSGEN (KR,1,0.5,FDEN,TCOS(K4)) MERGE (TCOS,K3,JR,K3+JR,KR,K1) 0) GO TO 178 COSGEN (LR,1,0.S,EDEN,TCOS(K4)) MERGE (TCOS,K3,JR,K3+JR,LR,K3-LR) COSGEN (KR,1,0.5,FDEN,TCOS(K4)) KR KR TRI3 (MR,A,BB,C,K,3,32,33,TCOS,D,W,W2,W3) DO 179 II1,MR B(I) I B(I)+32(I)+B3(I) CONTINUE TCOS(1) I 2. CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 180 II1,MR Q(I,J) - Q(I,J)+B(I) B(I) - Q(I,1)+2.*Q(I.J) CONTINUE CALL CALL IF (JR .NE. COSGEN (JR,1,0.5,0.0,TCOS) TRIX (JR, 0,MR,A,BB,C, B, TCOS,D,W) 1) GO TO 182 DO 181 II1,MR Q(I,1) I B(I) CONTINUE GO TO 194 CONTINUE DO 183 II1,MR Q(Ill) I .5*Q(I,1)-Q(I,JM1)+B(I) CONTINUE GO TO 194 CONTINUE IF (N .NE. 2) GO TO 188 C*** C*** C*** Case N I 2 C*** 185 D0 185 II1,MR 3(I) I Q(I,1) CONTINUE TCOS(1) I 0. CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 186 II1,MR 186 187 188 c*** C*** 277 Q(Icl) I B(I) B(I) I 2.*(Q(I,2)+B(I))*FISTAG CONTINUE TCOS(1) I IFISTAG TCOS(Z) I 2. CALL TRIX (2,0,MR,A,BB,C,B,TCOS,D,W) DO 187 II1,MR Q(Ill) " Q(I11)+B(I) CONTINUE JR I 1 IZR I 0 GO TO 194 CONTINUE C*** C*** Case of general N and NR I 2 C*** 189 190 191 192 193 194 C*** DO 189 II1,MR II I IP+I 33(1) I 0. B(I) I Q(I,1)+2.*P(II) Q(II 1) a .5*Q(II1)-Q(IIJM1) 32(I) I 2.*(Q(I,1)+Q(I,NLAST)) CONTINUE K1 I KR+JRI1 TCOS(K1+1) I I2. K4 I K1+3IISTAG CALL COSGEN (KR+ISTAG-2,1,0.0,FNUM,TCOS(K4)) K4 - K1+KR+1 CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K4)) CALL MERGE (TCOS,K1,KR,K1+KR,JR-1,0) CALL COSGEN (KR,1,0.5,FDEN,TCOS(K1+1)) K2 I KR K4 I K1+K2+1 CALL COSGEN (LR,1,0.5,FDEN,TCOS(K4)) K3 I LR K4 I 0 CALL TRI3 (MR,A,BB,C,K,3,32,B3,TCOS,D,W,W2,W3) DO 190 II1,MR B(I) I B(I)+32(I) CONTINUE TCOS(1) I 2. CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 191 II1,MR Q(I,1) I Q(I'1)+B(I) CONTINUE GO TO 194 DO 193 II1,MR B(I) I Q(I,NLAST) CONTINUE GO TO 196 CONTINUE C*** 278 C*** C*** Start back substitution : C*** 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 J I NLAST-JR DO 195 II1,MR B(I) I Q(I,NLAST)+Q(I,J) CONTINUE JM2 I NLAST-IZR IF (JR .NE. 1) GO TO 198 DO 197 II1,MR Q(I,NLAST) I 0. CONTINUE GO TO 202 CONTINUE IF (NROD .NE. 0) GO TO 200 D0 199 II1,MR II I IP+I Q(I,NLAST) I P(II) CONTINUE IP I IP-MR GO TO 202 DO 201 II1,MR Q(I,NLAST) I Q(I,NLAST)IQ(I,JM2) CONTINUE CONTINUE CALL COSGEN (KR,1,0.5,FDEN,TCOS) CALL COSGEN (LR,1,0.5,FDEN,TCOS(KR+1)) IF (LR .NE. 0) GO TO 204 DO 203 II1,MR B(I) I FISTAG*B(I) CONTINUE CONTINUE CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W) DO 205 II1,MR Q(I,NLAST) I Q(I,NLAST)+B(I) CONTINUE NLASTP I NLAST CONTINUE JSTEP I JR JR I I2R I2R I I2R/2 IF (JR .EQ. 0) GO TO 222 GO TO (207,208),MIXBND JSTART I 1+JR GO TO 209 JSTART I JR CONTINUE KR I KR-JR IF (NLAST+JR .GT. N) GO TO 210 KR I KR-JR NLAST I NLAST+JR JSTOP I NLAST-JSTEP GO TO 211 210 211 212 213 214 215 216 217 218 219 220 221 222 C*** C*** 279 CONTINUE JSTOP I NLAST-JR CONTINUE LR I KRIJR CALL COSGEN (JR,1,0.5,0.0,TCOS) DO 221 JIJSTART,JSTOP,JSTEP JM2 I J-JR JPZ I J+JR IF (J .NE. JR) GO TO 213 DO 212 III,MR B(I) I Q(IIJ)+Q(IIJP2) CONTINUE GO TO 215 CONTINUE DO 214 II1,MR 3(1) I Q(IIJ)+Q(IIJM2)+Q(IIJP2) CONTINUE CONTINUE IF (JR .NE. 1) GO TO 217 DO 216 II1,MR Q(I,J) I 0. CONTINUE GO TO 219 CONTINUE JMI I JIIZR JPI I J+I2R DO 218 II1,MR Q(IIJ) I .5*(Q(I,J)-Q(I,JMl)-Q(I,JP1)) CONTINUE CONTINUE CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) DO 220 II1,MR Q(IIJ) I Q(IIJ)+B(I) CONTINUE CONTINUE NROD I 1 IF (NLAST+I2R .LE. N) NROD I 0 IF (NLAST? .NE. NLAST) GO TO 194 GO TO 206 CONTINUE C*** C*** Return storage requirements for P vectors C*** C*** C*** W(I) I IPSTOR RETURN END C*********************************************************************** SUBROUTINE POISP2 (M,N,A,BB,C,Q,IDIMQ,B,B2,B3,W,W2,W3,D,TCOS,P) C*********************************************************************** C*** ** 280 C*** SUBROUTINE TO SOLVE POISSON EQUATION WITH PERIODIC BOUNDARY ** C*** CONDITIONS. ** C*** ** C*********************************************************************** C*** DIMENSION A(1) ,33(1) ,C(1) ,Q(IDIMQ,1) , 1 8(1) ,BZ(1) .83<1) .W(1) , 2 w2(1) ,w3(1) ,D(1) ,TCOS(1) , 3 9(1) C*** C*********************************************************************** C*** MR I M NR - (N+1)/2 NRMl I NR-l IF (2*NR .NE. N) GO To 107 C*** C*** C*** Even number of unknowns : C*** C*** DO 102 JI1,NRM1 NRMJ I NRIJ NRPJ I NR+J DO 101 II1,MR 5 I Q(IINRMJ)-Q(IINRPJ) T I Q(I,NRMJ)+Q(I,NRPJ) Q(I,NRMJ) I S Q(I,NRPJ) I T 101 CONTINUE 102 CONTINUE DO 103 II1,MR Q(IINR) I 2.*Q(I:NR) Q(I,N) I 2.*Q(I,N) 103 CONTINUE CALL POISD2 (MR,NRM1,1,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P) IPSTOR I W(l) CALL POISN2 (MR,NR+1,1,1,A,BB,C,Q(1,NR),IDIMQ,B,BZ,B3,W,W2,W3,D, l TCOS,P) IPSTOR I MAXO(IPSTOR,INT(W(1))) DO 105 JI1,NRM1 NRMJ I NRIJ NRPJ I NR+J DO 104 II1,MR S I .5*(Q(I,NRPJ)+Q(I,NRMJ)) T I .5* (Q(I,NRPJ)'Q(I,NRMJ)) Q(IINRMJ) - S Q(IINRPJ) I T 104 CONTINUE 105 CONTINUE DO 106 II1,MR Q(I,NR) I .5*Q(I,NR) Q(IpN) I 05*Q(IIN) 105 CONTINUE 281 GO TO 118 107 CONTINUE C*** C*** C*** Odd number Of unknowns : C*** C*** DO 109 J-1,NRM1 NRPJ - N+1-J DO 108 I-1,MR S - Q(I,J)-Q(I,NRPJ) T I Q(I,J)+Q(I,NRPJ) Q(IIJ) ' S 108 CONTINUE 109 CONTINUE DO 110 I-1,MR Q(I.NR) = 2.*Q(I,NR) 110 CONTINUE LH - NRMl/Z DO 112 J-1,LH NRMJ - NR-J DO 111 I-1,MR S ‘ Q(IIJ) Q(I,J) - Q(IpNRMJ) Q(IINRMJ) a S 111 CONTINUE 112 CONTINUE CALL POISDZ (MR,NRM1,2,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P) IPSTOR - W(l) CALL POISN2 (MR,NR,2,1,A,BB,C,Q(1,NR),IDIMQ,B,82,BB,W,W2,W3,D, 1 TCOS,P) IPSTOR - MAXO(IPSTOR,INT(W(1))) DO 114 J-1,NRM1 NRPJ - NR+J DO 113 I-1,MR S - .S*(Q(I,NRPJ)+Q(I,J)) T - .5*(Q(I,NRPJ)-Q(I.J)) Q(I,NRPJ) - T Q(IIJ) ' S 113 CONTINUE 114 CONTINUE DO 115 I-1,MR Q(I,NR) - .5*Q(I,NR) 115 CONTINUE DO 117 J-1,LH NRMJ - NR-J DO 116 I-1,MR S - Q(I,J) Q(I,J) - Q(IINRMJ) Q(IINRMJ) ' S 116 CONTINUE 117 CONTINUE 118 CONTINUE 282 C*** C*** C*** Return storage requirements for P vectors : C*** C*** W(l) I IPSTOR RETURN END C*** C*** C*********************************************************************** SUBROUTINE TRIX(IDEGBR,IDEGCR,M,A,B,C,Y,TCOS,D,W) C*********************************************************************** C*** ** C*** SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS WHERE THE ** C*** COEFFICIENT MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY ** C*** TRIDIAGONAL ( . . . , A(I), B(I), C(I), . . . ). ** C*** ** C*********************************************************************** C*** DIMENSION A(1) ,B(1) ,C(1) ,Y(l) , l TCOS(1) ,D(1) ,W(1) C*** C*********************************************************************** C*** MMl I M-l EB I IDEGBR+1 PC I IDEGCR+1 L I FB/FC LINT I 1 DO 108 KI1,IDEGBR X I TCOS(K) IF (K .NE. L) GO TO 102 I I IDEGBR+LINT XX I X-TCOS(I) DO 101 II1,M W(I) I Y(I) Y(I) I XX*Y(I) 101 CONTINUE 102 CONTINUE Z I 1./(B(1)-X) 0(1) I C(1)*Z Y(1) I Y(1)*Z DO 103 II2,MM1 Z I 1./(B(I)-X-A(I)*D(I-l)) 9(1) I C(I)*Z 'Y(I) I (Y(I)-A(I)*Y(I-l))*Z 103 CONTINUE z I B(M)—X-A(M)*D(MM1) IF (2 .NE. 0.) GO TO 104 Y(M) I 0. GO TO 105 104 Y(M) I (Y(M)-A(M)*Y(MM1))/Z 105 CONTINUE 283 DC 106 IP=1,MM1 I I M-IP Y(I) I Y(I)-D(I)*Y(I+1) 106 CONTINUE IF (K .NE. L) GO To 108 D0 107 II1,M Y(I) - Y(I)+W(I) 107 CONTINUE LINT - LINT+1 L - (FLOAT(LINT)*FB)/FC 108 CONTINUE RETURN END C*** C*** C*********************************************************************** SUBROUTINE TRI3 (M,A,B,C,K,Y1,Y2,Y3,TCOS,D,W1,W2,W3) C*********************************************************************** C*** ** C*** SUBROUTINE TO SOLVE THREE LINEAR SYSTEMS WHOSE COMMON COEFFICIENT** C*** MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY ** C*** TRIDIAGONAL (...,A(I),B(I),C(I),...) ** C*** ** C*********************************************************************** C*** DIMENSION A(1) ,B(l) ,C(1) ,K(4) , 1 TCOS(1) ,Y1(l) ,Y2(l) ,Y3(l) , 2 D(l) ,W1(1) ,W2(1) ,w3(1) C*** C*********************************************************************** C*** MMl I M-l K1 I X(1) K2 I X(2) K3 I K(3) K4 I K(4) F1 I Kl+1 F2 I K2+1 F3 I K3+1 F4 I K4+1 K2K3K4 I K2+K3+K4 IF (K2K3K4 .EQ. 0) GO TO 101 L1 I Fl/FZ L2 I F1/F3 L3 I F1/F4 LINTl I 1 LINTZ I 1 LINT3 I 1 KINTl I K1 KINTZ I KINT1+K2 KINT3 I KINT2+K3 101 CONTINUE DO 115 NI1,K1 X I TCOS(N) 102 103 104 105 106 107 108 109 110 111 112 113 284 IF (K2K3K4 .EQ. 0) GO TO 107 IF (N .NE. L1) GO To 103 D0 102 1=1,M w1(1) = 21(1) CONTINUE IF (N .NE. L2) GO TO 105 DO 104 1:1,M w2(1) a 22(1) CONTINUE 1E (N .NE. L3) GO TO 107 D0 106 II1,M W3(I) s 23(1) CONTINUE CONTINUE z - 1./(B(1)-X) 0(1) - C(1)*Z 21(1) - 21(1)*z 22(1) = 22(1)*z 23(1) = 23(1)*z DO 108 I=2,M z - 1./(B(I)-X-A(I)*D(I-l)) 0(1) = C(I)*Z 21(1) = (Y1(I)-A(I)*Y1(I-1))*Z 22(1) = (Y2(I)-A(I)*Y2(I-l))*z 23(1) - (Y3(I)-A(I)*Y3(I-1))*Z CONTINUE DO 109 IPI1,MM1 I - M-IP 21(1) = 21(1)-D(I)*21(1+1) 22(1) - 22(1)-D(1)*22(1+1) 23(1) - Y3(I)-D(I)*Y3(I+l) CONTINUE IF (K2K3K4 .EQ. 0) GO TO 115 IF (N .NE. L1) GO TO 111 1 - LINT1+KINT1 xx - x-TCOS(1) DO 110 II1,M 21(1) - xx*21(1)+w1(1) CONTINUE L1NT1 - L1NT1+1 L1 - (FLOAT(LINT1)*F1)/F2 IF (N .NE. L2) GO TO 113 1 - LINT2+KINT2 xx - X-TCOS(I) DO 112 1-1,M 22(1) - xx*22(1)+w2(1) CONTINUE LINTZ - LINT2+1 L2 - (FLOAT(LINT2)*F1)/F3 IF (N .NE. L3) GO TO 115 1 - LINT3+KINT3 xx - x-TCOS(1) DO 114 1=1,M 23(1) - xx*23(1)+w3(1) 285 114 CONTINUE LINT3 I LINT3+1 L3 I (FLOAT(LINT3)*F1)/F4 115 CONTINUE RETURN END BIBLIOGRAPHY 10. BIBLIOGRAPHY Abernathy, F. H. and Kronauer, R. E., "The formation of vortex streets", Journal of Fluid Mechanics, Vol. 13, p. l, (1962). Ashurst, W. T., "Numerical Simulation of Turbulent Mixing Layers Via Vortex Dynamics", Turbulent Shear Flows I, edited by F. Durst et. al., Springer- Verlag, Berlin, pp. 402-413, (1979). Ashurst, W. T., "Calculation of Plane Sudden Expansion Flow via Vortex Dynamics", 2nd Shear Flow, London, (1979). Ashurst , W. T., "Piston-Cylinder Fluid Motion via Vortex Dynamics", Fluid Mechanics of Combustion Systems, edited by T. Morel et. al., Fluid Engineering Division, ASME, Boulder Colorado, (1981). Ashurst, W. T., "Lagrangian-Eulerian Calculation of Turbulence Effects on a Premixed Flame Propagation", presented at the 1981 Spring Meeting of the Western States Section of the Combustion Institute, paper WSS/CI/81-23, (1981). Baker, G. R., "The Cloud-In-Cell Technique applied to the Roll-up of Vortex Sheet", Journal of Computational Physics, Vol.31, pp. 76-95, (1979) Bartholomew, R. W., Private communication. (1985) Beal, J. T. and Majda, A., "Vortex Methods II : Higher Order Accuracy in Two and Three Dimensions", Mathematics Of Computation, Vol. 39, pp. 29-52, (1982). Birdsall, C. K. and Fuss, D., "Clouds-In-Clouds, Clouds-In-Cells Physics for Many-Body Plasma Simulation", Journal of Computational Physics, Vol.3, pp. 494w511, (1969). Birkhoff, G. and Fisher, J., "DO Vortex Sheets Roll Up?", Rend. Circ. Mat. Palermo Ser. 2, 8, p. 77, (1959). 286 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 287 Bowman. F.. 'Immflmmwth—Amhm" Dover Publications, Inc., New York, (1961). Brown, G. L. and Roshko, A. J., "On Density Effects and Large Structures in Turbulent Mixing-Layer Growth at Moderate Reynolds Num ", Journal of Fluid Mechanics, Vol. 64, pp. 775-816, (1974). Burt, R. and Truth, K. A., "Penetration and Vaporization of Diesel Fuel Sprays", Symposium on "Diesel Engine Combustion", Institution of Mechanical Engineers, (1970). Carslaw, H. S. and Jaeger, J. C., Conduction in Solids, Oxford University Press, London, England, p. 256, (1959). Chorin, A. J., "Numerical Methods for use in Combustion Modeling", Proceedings of the International Conference on Numerical Methods in Science and Engineering, Versailles, France, pp. 229-236, (1979). Chorin, A. J ., "Numerical Studies of Slightly Viscous Flow", Journal of Fluid Mechanics, Vol. 57, pp. 785-796, (1973). Chorin, A. J ., "Vortex Sheet Approximation of Boundary Layers", Journal of Computational Physics, Vol. 27, pp. 428-442, (1978). Chorin, A. J., "Vortex Models and Boundary Layer Instabilities", SIAM Journal of Scientific and Statistical Computations, Vol. 1, pp. 1-24, (1980). Chorin, A. J. and Bernard, P. S., "Discretization of a Vortex Sheet, with an Example of Roll-Up", Journal of Computational Physics, Vol. 13, pp. 423-429, (197 3). Christiansen, J. P., "Numerical Simulation of Hydrodynamics by the Method of Point Vortices", J. Comp. Phys., Vol. 13, pp. 363-379, (1973). Clements, R. R. and Maull, D. J ., "The Representation of Sheets of Vorticity by Discrete Vortices", Prog. Aeronaut. Sci., Vol. 6, pp. 129-146, (1975). Corrsin, 8., "Simple Theory of an Idealized Turbulent Mixer", A.I.Ch.EJ. Vol. 3, p. 329., (1957) Couet, B., Buneman, O. and Leonard, A., "Simulation of Three-Dimensional Incompressible Flows with Vortex-In-Cell Methods", J. Comp. Phys., Vol. 39, pp. 305—328, (1981). 24. 25. 26. 27. 28. 29. 30. 31. 32 33. 34. 35. 288 Couet, B. and Leonard, A., "Exact Extension to the Three-Dimensional Incompressible Flows with Vortex Elements", Ann. Rev. Fluid Mech., Vol. 17. PP. 523-559, (?). Courant, R., Friedrichs, K., and Lewy, H., "Uber die Partiellen Differenzengleichungen der Matthematischen Physik", IBM J. p.215, (March 1967); Math. Ann., Vol. 100, p.32, (1928). Curl, R. L., "Dispersed Phase Mixing : 1. Theory and Effects of Simple Reactors", A.I.Ch.E.J., Vol. 9, pp. 175-181, (1963). Dimotakis, P. E. and Brown, G. L., "The Mixing Layer at High Reynolds Number : Large Structure Dynamics and Entrainment", Journal of Fluid Mechanics, Vol. 78, p.535, (1976). Dopazo, C. and O’Brien, E. E., Acta Astronaut, Vol. 1, p. 1239, (1974). Dopazo, C. and O’Brien, E. E., 4th. Int. Coll. on Gasdynamics of Explosions and Reactive Systems, La Jolla, (1973). Damn, C. and O’Brien, E. E., "Statistical Treatment of Non-isothermal Chemical Reactions in Turbulence", Combustion Science and Technology, Vol. 13, p. 99, (1976). Dopazo, C., "Probability Density Function Approach for a Turbulent Axissymmetric Heated Jet Centerline Evolution", Phys. Fluids, Vol. 22, pp. 20—30, (1979). Dorodnicyn, A. A., "Review of methods for solving Navier-Stokes equations", needing of 3rd htgmafional Conference on Numerical Methods in Fluid Mechanics, Lecture Notes in Physics, Springer-Verlag, pp. 1-11. (1973). Dyer, T. M. "Characterization of One- and Two-Dimensional Homogeneous Combustion Phenomena in a Constant Volume Bomb", Sandia National Labs, Report SAND78-8704. Einstein, A., Investigation on the Theog of the Brownian Movement, Translation, Methuen and CO., Ltd., London, (1926); Reprint, Dover Publications, Inc., New York, (1956). Evangelista, J. J. , Shinnar, R. and Katz, in Twelfth Symposium International on Combustion. (The Combustion Institute, Pittsburg, Pa.), p. 401, (1968). 36. 37. 38. 39. 40. 41. 42. 43. 45. 46. 47 289 Flagan, R. C. and Appleton, J. P., "A Stochastic Model of Turbulent Mixing with Chemical Reaction: Nitric Oxide Formation in a Plug-Flow Burner", Combustion and Flame, Vol. 23, p. 249, (1974). Fox, R. L., "Multipoint Distribution Function Hierarchy for Compressible Turbulent Flows", Physics of Fluids, Vol. 18, p. 1245, (1975). Gasman, A. D, Johns, R. J. R. and Watkins, A. P., "Development Of Prediction Methods for In-Cylinder Processes in Reciprocating Engines", WM edited by J N Martavi and C. A. Amann, Plenum Press, New York, pp. 69- 129, (1980). Hald, C. and Del Prete, V. M., "Convergence of Vortex Methods for Euler’s Equations", Mathematics of Computation, Vol. 32, pp. 791-809, (1978). Hammersley, J. M. and Handscomb, D. C., Monte Cglo Methods, Metheun, London, (1964). Hockney, R. W., Goel, S. P. and Eastwood, J. W., "Quiet High-Resolution Computer Models of a Plasma", J. Comp. Phys., Vol. 14, pp. 148-158, (1974). Janicka, J., Kolbe, W. and Kolmann, W., "The Solution of a PDF-Transport Equation For Turbulent Diffusion Flames", in Proceedings of the Heat Transfer in Fluid Mechanics Institute (Standford University Press, Standford), p. 296, (1978). Kraichnan, R. H. and Montgomery, D., "Two-dimensional turbulence", Rep. Prog. Phys., Vol. 43. PP. 547-619, ( 1980). Lamberti, J. Probability, Benjamin, New York, (1966) Leonard, A., "Vortex Methods for Flow Simulations", Journal of Computational Physics, Vol. 37, pp. 289-335, (1980). Levenspiel, O. and Spielman, L. A., "A Monte Carlo Treatment for Reacting and Coalescing Dispersed Phase Systems", Chem. Engng. Sci., Vol. 20, p. 247, (1965). Lundgren, T. 8., "Distribution Functions in the Statistical Theory of Turbulence", Physics Of Fluids, Vol. 10, p. 969, (1967). 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 61. 290 Lyn, W. T., "Study of Burning Rate and Nature of Combustion in Diesel Engines", 9th. Symposium of Combustion, pp. 1069-1082, (1960). Meng, J. C. S. and Thomson, J. A. L., "Numerical Studies of some Non-Linear Hydrodynamic Problems by Discrete Vortex Element Methods", Journal of Fluid Mechanics, Vol. 84, part 3, pp. 433-453, (1978). Milinazzo, F. and Saffman, P. G., "The Calculation of Large Reynolds Number Two-Dimensional Flow Using Discrete Vortices with Random W ", Journal of Computational Physics, Vol. 23, pp. 380-392, (1977). Milne-Thomson, L. M., "Hydrodynamical Images", Proceedings of the Cambridge Philosophical Society, Vol. 36, pp. 246-247, (1940). Moore, D. W., "The Discrete Vortex Approximation of a Finite Vortex Sheet", California Institute of Technology, Report AFOSR-1804-69, (1971). Newman, J. A. and Brzustowski, T. A., "Behaviour of a Liquid Jet Near the Thermodynamic Critical Region", AIAA Journal, Vol. 9, 1971, NO. 8. Pope, S. B., Philos. Trans. R. Soc. London, Ser. A291, 529 (1979). Pope, S. B., "The Probability Approach to the Modeling of Turbulent Reacting Flow", Combustion and Flame, Vol. 27, p. 299, (1976). Pope, S. B., "A Monte Carlo Method for the PDF Equations of Turbulent Reactive Flow", Combustion Science and Technology, Vol. 25, pp. 159-174, (1981). Pratt, D. T., Prog. "Mixing and Chemical Reaction in Continuous Combustion", Energy Combustion Sci., Vol. 1, p. 73, (1976). Rhines, P. B., "Geopstrophic turbulence", Ann. Rev. Fluid Mech., Vol. 11, pp. 401-41, (1979). Rife, J. and Heywood, J. B., "Photographic and Performance Studies of Diesel Combustion with a Rapid Compression Machine", SAE Transactions, Vol. 83, Section 4, pp. 2942-2961, (1974). Rosenhead, L., Proc. Royal Soc. Vol. A134, p. 170, (1931). Roshko, A. J., "Structures and Turbulent Shear Flows : a new 100 ", AIAA Journal, Vol. 14, pp. 1349-1357, (1976). 62. 63. 65. 66. 67. 68. 69. 70. 71. 72. 73. 291 Rossow, V. J., "Convective merging of vortex cores in lift-generated wakes", J. Aircr., Vol. 14. PP. 283-90, (1977). Shahed, S. M., Chin, W. S. and Lyn, W. T., "A Transient Spray Mixing Model for Diesel Combustion", A Mathematical model of diesel combustion. Conference on Combustion Engines, Cranfield, (July, 1975). Shestakov, A. 1., "Numerical Solution of the Navier-Stokes Equations at High Reynolds Numbers", Ph.D. Thesis, Department of Mathematics, University of California at Berkeley, (1975). Shreider, Yu A., (Ed.) The Mente Q10 Method, Pergamon, London, (1966). Spalding, D. B., "A General Theory of Turbulent Combustion, the Lagrangian Aspects", AIAA 15th. Aerospace Science Meeting, L. A., California, Report 77-141. Stratonovich, R. L., "Topics in the Theory of Random Noise", New York, Vol. 1, Chapters 1 & 4, (1963). Swartztrauber, P., Sweet, R. and Adams, J. Source listings for FISHPAK from the National Center for Atmospheric Research, Boulder, CO., (1986). Takami, H., "Numerical Experiment with Discrete Vortex Approximation, with Reference to the Rolling Up of a Vortex Sheet", Dept. of Aero. and Astro., Standford University, Report SUDAER 202, (1964). Wang, S. S., "Grid-Insensitive Computer Simulation of the Kelvin-Helmholtz Instability and Shear Flow Turbulence", Ph.D. Thesis, Standford University Institute for Plasma Research Report No. 710, (1977). Westwater, F. L., British A. R. C. Reports and Memoranda No. 1692, p.116, (1936). Winant, C. D. and Browand, F. K., "Vortex Pairing : the Mechanism of Turbulent Mixing Layer Growth at Moderate Reynolds Numbers", Journal of Fluid Mechanics, Vol. 63, part 2, pp. 237-255, (1974). Yanenko, N. N., The Methg ef befional Steps, Springer Verlag, New York, (1971). 292 74 Zabusky, N. J., "Coherent Structures in Fluid Dynamics". W Nonlinearity in the Natural Sciences. Ed. A. Perlmutter, L. F. Scott, pp. 45-205, New York : Plenum. 423 pp. (1977). 75. Zabusky, N. J ., "Recent Developments in Contour Dynamics for the Euler Equations", Ann. N. Y. Acad. Sci., Vol. 373. pp. 160-170, (1981).