UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The discrete velocity method in plasma physics Irwin, Andrew J. 1993

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1993_fall_irwin_andrew.pdf [ 3.93MB ]
JSON: 831-1.0080528.json
JSON-LD: 831-1.0080528-ld.json
RDF/XML (Pretty): 831-1.0080528-rdf.xml
RDF/JSON: 831-1.0080528-rdf.json
Turtle: 831-1.0080528-turtle.txt
N-Triples: 831-1.0080528-rdf-ntriples.txt
Original Record: 831-1.0080528-source.json
Full Text

Full Text

THE DISCRETE VELOCITY METHOD IN PLASMA PHYSICSByAndrew J. IrwinB. Sc. (Chemical Physics), University of Toronto, 1991.A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE STUDIESINSTITUTE OF APPLIED MATHEMATICSANDDEPARTMENT OF MATHEMATICSWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAAugust 1993© Andrew J. Irwin, 1993In presenting this thesis in partial fulfilment of the requirements for an advanced degree atthe University of British Columbia, I agree that the Library shall make it freely availablefor reference and study. I further agree that permission for extensive copying of thisthesis for scholarly purposes may be granted by the head of my department or by hisor her representatives. It is understood that copying or publication of this thesis forfinancial gain shall not be allowed without my written permission.Institute of Applied MathematicsandDepartment of MathematicsThe University of British Columbia1984 Mathematics Road, Vancouver, CanadaV6T 1Z2Date: August 31, 1993.AbstractThe objective of this thesis is to apply a discrete velocity model to kinetic theoryproblems in plasma physics. Numerical approaches commonly used in kinetic theoryare described and compared with the discrete velocity models. The discrete Boltzmannequation (DBE) is a commonly used discrete velocity method for problems in rarefiedgas dynamics and is adapted for plasma physics in this thesis.The Boltzmann equation is used to model the relaxation to equilibrium of a pureelectron plasma. The first of two problems studied is the relaxation of test particles ina Maxwellian bath. This is a linear version of the second problem and serves as a testexperiment for the method and the computer code. The second problem is the nonlinearrelaxation of an anisotropic velocity distribution of electrons due to self-collisions. Thisphysical situation arises in many natural phenomena, such as atmospheric and spaceplasmas, as well as in many laboratory investigations.Pure electron plasmas have been the subject of many experiments, including studies oftime dependent transport properties and studies of the relaxation of anisotropic velocitydistributions. The Fokker-Planck and Boltzmann equations are commonly used as thetheoretical starting point for studies of plasmas.iiTable of ContentsAbstract^ iiTable of Contents^ iiiList of Figures vList of Tables^ viii1 Introduction 11.1 Background  ^11.2 Discrete Velocity Models  ^41.3 Kinetic Theory  ^51.3.1 Equations for the problems studied  ^71.3.2 Limitations of the Boltzmann equation ^  101.4 Plasma Physics ^  131.4.1 Confinement Theorem ^  141.4.2 Experiments with pure electron plasmas ^  161.5 Thermalization of a Test Particle in a Heat Bath  201.6 Relaxation of Anisotropic Distributions ^  211.7 Other Numerical Approaches^  221.7.1 Simulations  ^231.7.2 Related methods for plasma physics ^  252 Numerical Methods^ 291112.1 Introduction ^  292.2 The Discrete Boltzmann Equation ^  292.3 Numerical Integration ^  312.4 Distribution Function Considerations ^  342.5 Discrete Cross Sections ^  353 Calculations and Results 403.1 Introduction ^  403.2 Definitions of computational quantities ^  413.2.1 Cross-sections  ^443.2.2 Sample results ^  473.3 Comparison with expansion method ^  513.4 Computational Results ^  623.4.1 Isotropic Cross-sections  ^623.4.2 'Rutherford-type' cross-sections  ^703.5 Discussion ^  754 Summary and Conclusions^ 79Bibliography^ 81A Discrete Boltzmann Equation Program Listing^ 88ivList of Figures1.1 Comparison of a continuous and a discrete distribution function ^31.2 Geometry of a collision.  ^81.3 The Rutherford cross-section, truncated near where it diverges. ^ 121.4 Schematic of basic electron plasma experimental apparatus. ^ 151.5 Density profiles and rotation profiles of a column of electron plasma relax-ing to equilibrium. (From Malmberg et al. 41.) ^  171.6 Experimentally measured temperature relaxation from a bi-Maxwelliandistribution. Symbols are experimental data points and the curve is fit tothe theoretical predictions of Ichimaru and Rosenbluth. (From Malmberget al. 41.) 192.1 The discrete velocities lie on nodes in the grid. The picture in three di-mensions is similar^  322.2 Discrete velocity collision sphere.  ^352.3 Angular part of Rutherford cross-section, o-(x) = 1/ sin4(x/2). The solidline is the analytic form, symbols show typical sample points used in thediscrete velocity computation  373.1 Representation of a bi-Maxwellian distribution. ^  453.2 Temperature relaxation for the linear problem with the hard sphere cross-section at different velocity resolutions, r = 5, 6, and 7. At short timesthe curves do not overlap. The uppermost curve corresponds to r = 5. . . 48v3.3 The relaxation of an anisotropic distribution function, initially bi-Maxwellianwith a(0) — 5^3.4 The evolution of TH and T1 for different velocity grid resolutions, r = 6,7, and 8, with a(0) = 5, using the hard sphere cross-section ^3.5 Partial sums of SN(x0) = E0 a(0)S(x) for xo = 0.5 and ai(0) = ((7 —xo is a reduced speed ^3.5 Partial sums of SN(x0) = Es1-%.7=0 a(0)S(x) for xo = 0.5 and a2(0) = ((7 —1)17)2. xo is a reduced speed ^3.6 Convergence of temperature relaxation with increasing numbers of Soninepolynomials. The initial temperature ratio is T1(0)/T2 = 1.75. The curvescorrespond to N = 2, 3, 4, 5 ^3.7 Comparison of temperature relaxation calculated with Sonine polynomialexpansions and DBE for initial temperature ratios T1(0)/T2 = 0.5, 1.5,1.75. The DBE velocity resolution is r = 7, and N = 8 Sonine polynomialsare used in the expansion ^3.8 Log-linear plot of (T1(t)/T2 — 1) versus t from the Sonine polynomial ex-pansion, with T1(0)/T2 = 5. The dashed line is a linear extrapolation ofthe straight tail. ^3.9 Test particle temperature relaxation for the hard sphere cross-section. Avelocity resolution of r = 7 was used. The labels T indicate initial tem-perature ratios, T1(0)/T2 ^3.10 Log-linear plot of (T1(t)/T2— 1) versus t for hard sphere cross-section fromthe DBE calculation. The labels T indicate initial temperature ratios,T1 (0)/T2. ^vi3.11 Test particle temperature relaxation for Maxwell molecule cross-sectionfrom the DBE calculation. A velocity resolution of r = 6 was used. Thelabels T indicate initial temperature ratios, T1(0)/T2 .  3.12 Log-linear plot of (T1(t)/T2 — 1) versus t for the Maxwell molecule cross-section. The T labels indicate the initial temperature ratio, T1(0)/T2. • •3.13 Relaxation of temperature anisotropy, a, for the hard sphere cross-section A velocity resolution of r = 7 was used 3.14 Relaxation of temperature anisotropy, a, for the Maxwell molecule cross-section. A velocity resolution of r = 7 was used^3.15 Log-linear plot of (711(t) — T)IT and (711(t) — T)IT versus t with a = 5for Maxwell molecules^3.16 Test particle temperature relaxation for the 1/94 cross-section. A velocityresolution of r = 7 was used. The labels T indicate the initial temperatureratio, T1(0)/T2.  3.17 Relaxation of temperature anisotropy, a, for the 1/g4 cross-section. Avelocity resolution of r = 7 was used. ^3.18 Test particle temperature relaxation for the angular part of the Rutherfordcross-section. A velocity resolution of r = 10 was used ^3.19 Relaxation of temperature anisotropy, a, for the angular part of the Ruther-ford cross-section. A velocity resolution of r = 10 was used. ^3.20 Experimental relaxation rate for a pure electron plasma. The symbols arefrom experiments with pure electron plasmas and the line is from theory.(From Malmberg et al. 41.)  viiList of Tables2.1 Flowchart for the DBE numerical algorithm^  302.2 Cross-sections used in the DBE ^  393.1 Timings for one iteration at various grid sizes. The times are directlyrelated to the number of collision terms; the time to compute collisiontables and start-up information is not included.   41VI"Chapter 1Introduction1.1 BackgroundThis thesis presents a study of discrete velocity models applied to kinetic theoryand plasma dynamics. Discrete velocity models are currently used to solve problems inrarefied gas dynamics. The main objective of the work is to test the applicability ofdiscrete velocity models as a method of solution of the Boltzmann equation.In 1872 Ludwig Boltzmann derived a nonlinear integro-differential equation to de-scribe the evolution of a gas. The Boltzmann equation has been tremendously successfulfor a wide range of problems. The basis of the physical description is a probability distri-bution function for the gas particles. Physical observables, such as particle density, bulkvelocity and temperature are defined in terms of moments of the distribution function.Perhaps the most famous result in kinetic theory is the Boltzmann H-theorem. Boltz-mann defined an entropy function, H, and showed that for closed systems described byBoltzmann equation dynamics, H increased with time. The irreversible dynamics drivethe distribution towards an equilibrium 'Maxwell-Boltzmann' (or Maxwellian) distribu-tion. Section 1.3 describes the Boltzmann equation and the physical equations used inthis thesis.The central notion behind discrete velocity methods is that an adequate approxima-tion to the distribution function can be achieved by discretizing velocities and making the1Chapter 1. Introduction^ 2distribution function constant over regions of moderately large size in velocity space. Fig-ure 1.1 illustrates the difference between an analytic distribution function and a discretedistribution. In this investigation the physical property of interest is the temperature.The temperature of a gas is directly proportional to the second velocity moment of thedistribution function. The goal of the discrete velocity model is to represent the distri-bution function accurately enough to describe the lower moments but coarsely enoughto achieve a large computational simplification of the equations.Discrete velocity methods are of interest because they are conceptually simple andare adaptable to a wide range of problems. The approximations, given in detail inchapter 2, are motivated by the microscopic physics of kinetic theory. Discrete velocitymethods are used to solve a diverse range of problems such as classic hydrodynamicproblems involving fluid flow past a plate, chemically reacting fluid mixtures and theescape of a planetary atmosphere. The methods have been used for problems governedby hydrodynamic equations with complicated geometries and in the rarefied gas (lowpressure) transition regime where kinetic descriptions become necessary.There are a variety of approaches to solving the Boltzmann equation numericallywhich do not rely on a discrete velocity approximation. The distribution function canbe expanded over a set of basis functions. The basis set is chosen according to the ge-ometry and symmetries of the problem so that in some situations the operators simplifyconsiderably when applied to the various terms. Another technique discretizes variablesto facilitate the solution of the Fokker-Planck differential form of the Boltzmann opera-tor. These finite-difference methods are commonly used for problems with complicatedgeometries or for problems which have no convenient symmetries to make expansionmethods advantageous. Other methods, such as the discrete simulation Monte Carlo andlattice gases, simulate particles in a discrete velocity space. Although these methods donot directly solve the differential equations they can be considered equivalent in some(a) Continuous Maxwellian(b) Discrete MaxwellianChapter 1. Introduction^ 3Figure 1.1: Comparison of a continuous distribution function with a discretized versiontypical of one used in the DBE. The distribution function has been integrated over vy.The x and z labels on the axes refer to the vz. and vz components of velocity in reducedunits.Chapter 1. Introduction^ 4cases.1.2 Discrete Velocity ModelsBroadwell" initiated the widespread interest in discrete velocity models by using anapproximation of the Boltzmann equation with very few velocities (only 6 or 8) to solveproblems for rarefied gases. This model is known as the Discrete Boltzmann equation(DBE) or earlier as the Broadwell model. Broadwell, Gatignol" and Cabannes5 wereamong the first to use the DBE to investigate shock wave propagation and the Couetteand Rayleigh flow problems. These are fundamental problems in the kinetic theory ofgases which test our understanding of the macroscopic physical properties of a gas andhow they relate to the microscopic physics, especially the effects of collisions, convectiveflow and external forces. The Couette and Rayleigh flow problems concern a gas flowingeither above a moving plate or sandwiched between two moving plates with oppositevelocities. The problem is to discover the density and flow velocity of the gas near theplate-gas boundaries, the velocity profile above the plate, the shear stress, and the slip atthe plate. In the shock problem the goal is to find the density profile on each side of theshock. These investigations were concerned with the case where the density of the gasis so low that the standard hydrodynamic equations would no longer apply. The DBEhas also been applied to the heat conduction problem which models the temperature ofa gas between two plates of different temperatures. Much of this early work focused onobtaining analytic results and the computations done were quite simple and used onlya small number of velocities. A review of the mathematical development of the DBE isgiven by Platowski and Inner.' Much of the recent work in discrete kinetic theory hasfocussed on the mathematical aspects of the theory. Contributors to a recent symposium7discuss questions such as the fluid dynamic limits of the discrete kinetic theory, existenceChapter 1. Introduction^ 5of solutions under various conditions and extensions to the basic method for includingmany body collisions or chemical reactions."The numerical technique used in this thesis is the DBE. A large number of velocities(, 103) are used to define the discrete representation of the distribution function, butunlike the studies mentioned above the applications in this thesis are to problems withno spatial variation. The treatment follows the recent calculations by Inamuro andSturtevant,1° who were the first to attempt a calculation with a large number of velocities,and Merryfield,11 who applied the method to the problem of the escape of a planetaryatmosphere. The details of the method will be described in chapter 2.The goal of this thesis is to add to the theoretical investigations and recent computa-tions by applying the DBE to a new problem and to test the usefulness of discrete velocitymodels in kinetic theory and plasma physics. Good numerical methods which are easy toapply are needed to model plasmas to facilitate the comparison of experimental resultswith theory.1.3 Kinetic TheoryAs one of many possible applications, I choose to study the relaxation to equilibrium ofan electron plasma in the absence of external forces. As mentioned above the descriptionbeing used is a probabilistic description of the electron plasma. This section and thenext provide a brief introduction to the Boltzmann equation used to model the evolutionof plasmas and some basic facts about plasma physics. More details can be found in thestandard references. 12-17The electron plasma under consideration is modeled with a distribution function fChapter 1. Introduction^ 6defined byf(x, v, t) dx dv = the number of particles in a volume^(1.1)dx centred on x in position space anda volume dv centred on v in velocityspace.Macroscopic physical properties are defined as moments of the distribution function.For example, particle density, n, bulk velocity, c, and temperature, T, are given byn(x,t) = (1) =^f(x,v,t) (1.2)c(x,t) = jn-(v, ) =^f vf(x,v,t) dv, (1.3)—3 kT(x, t)2= —1( —1 mv2 ) =^I lmv2 f(x v t) dv ,n 2^n^2(1.4)where k is Boltzmann's constant and m is the mass of the particles. All of these integralsare taken over the entire velocity space. In the discrete velocity computations the integralsare evaluated as sums over a compact subset of velocity space.The Boltzmann equation is often used to describe the evolution of distribution func-tions f(x, v, t) in kinetic theory. It can be described structurally as having two parts— the convective derivative on the left hand side and the collision operator on the righthand side. The evolution of species 1 due to the influence of species 2 is described byDfi+ F • Vvi =Dt^at +vl• (1.5)where fi = fi(x, v, t) and f2 = f2(x, v, t). The collision operator, L, describes theinteraction of two distribution functions as an integral over all possible collisions. Theintegral over v2 accounts for the collision of fi particles with velocity v1 with all 12particles. The weighting of the relative importance of each collision is determined bythe differential cross-section cr(g, 0) and is a function of the primed velocities resultingChapter 1. Introduction^ 7from a collision. The integral over solid angles 0 includes all possible outcomes from acollision. The collision operator is defined byL[fi, hi = I [L(vi).1.2(v) — fi(vi)f2(v2)] g cr(g, (1) dfl dv2,^(1.6)where fi is the distribution function for species i, cr(g , 0) is the cross-section describingthe collision of the particles, g = v2 — v1 is the relative velocity, g = 41, and F isthe applied force. Primed velocities are velocities resulting from a collision, unprimedvelocities refer to the incident velocities. 0 = (x, 0) is the solid angle between g and g',and c10 = sin x dx do if x is the polar angle and 0 is the azimuthal angle in standardspherical coordinates. Collision dynamics are restricted to the plane of the two collidingparticles and thus the integral over 0 just produces a factor of 27r. Figure 1.2(a) showsthe collision between two particles with velocities v1 and v2. The conservation of mo-mentum requires that the relative velocity after a collision, g', be a diameter of the circle.Figure 1.2(b) shows the centre of momentum frame, including the relative velocities, thescattering angle and the impact parameter b.A complete description of the dynamics for a two species system couples the evolutionof each species in a pair of nonlinear equations:D fiDt — ii[fl' f21?LaDt -= LV2, fil.1.3.1 Equations for the problems studiedThe first problem is the relaxation of a test electron in an infinite bath of equilibriumelectrons. All the particles are electrons, but there are still two distribution functions:the test particles (species 1) with an initial Maxwellian distribution of temperature TI.fi (v, t = 0; TO = ni m^(27rkTi) e-mv2 /2kT1 (1.7)Chapter 1. Introduction^ 8(a) Lab frame,.••• .......0g'/2^,.,,821.1ide.1^Itog/2 ...., L...•.........7.•••......•,....0..•Particle 2,,...................•...(b) Centre of momentum frameFigure 1.2: Geometry of a collision.Chapter 1. Introduction^ 9and the bath particles (species 2) with a Maxwellian distribution of temperature T2(independent of time),f 2(V ; T2) = n2 ^(2-k71T2) -^.772^p-mv2/2kT2^ (1.8)These distributions look like the function plotted in figure 1.1. With f2 fixed the problemis linear in fi and Boltzmann equation isaii(v,t) _at^L[ii,12].The velocity and space derivatives have disappeared from the left hand side because noexternal forces or spatial variations are included in this problem. The cross-section inthe collision operator is the Rutherford cross-section,u(g, x) =e4 X77/2g4 sin-4 (-2-) (1.10)since there are only electron-electron interactions.18 The cross-section would be the sameif only one ionic species were present. The mass of the electron is given by m, e is theelectron charge, and the scattering angle in the centre of momentum frame, x, is givenbycos(x) = (g • g')/g2.Since the test particles are assumed to be so dilute that they have no influence on thebath distribution, there is no equation for the time evolution of f2(v, t) which is takento be a time independent Maxwellian, equation (1.8). As a result equation (1.9) is linearin Ii and this first problem is much simpler than the main problem. This test particleproblem is included to test the DBE implementation and an expansion method will beused to check the results.The second problem is the study of the relaxation of an anisotropic velocity distri-bution function. There is only one species, with an initial bi-Maxwellian distribution(1.9)Chapter 1. Introduction^ 10function parameterized by two temperatures,(m ^e -mq /2kT, -my?, pull . (1.11)21-kTII 1/3TJ-2/3The two temperatures give the function oblate or prolate contours, see for example fig-ure 3.1. The velocity is resolved into a component parallel to the z axis and the remainder,perpendicular to i, v = v11 + vi = v11.& Vi(i + .). The Boltzmann equation for thisproblem isa fw(vi,t) = L[f] = f [f(v)f(v) — f (vi)f (v2)] g o-(g , x) di? dv2 .^(1.12)1.3.2 Limitations of the Boltzmann equationBoltzmann's original derivation of equation (1.5) relied on a physically intuitive argu-ment based on simple collision dynamics. More precise derivations from the Liouvilleequation and the BBGKY hierarchy of statistical mechanics have been worked out.13These derivations highlight the assumptions necessary for the Boltzmann equation to bevalid. In particular the particles are required to form a low density gas with short rangeforces so that only binary collisions are important.14,15 This is equivalent to neglectingthree particle correlations, meaning that many-particle interactions have the same effectas a series of binary collisions. The Coulomb potential, proportional to 1/r, is a long-range force which means that every particle always interacts with every other particle.Therefore the Boltzmann equation appears not to apply to a plasma with only long-rangeCoulomb forces.A common assumption in plasma physics is that large angle collisions are unimpor-tant relative to the more common small angle collisions. Multi-body weak (small-angle)collisions can be considered to be equivalent to a series of weak binary collisions.12,17Shoub19'2° has recently questioned this approximation for high energy test particles in aMaxwellian bath.fBliff(VII, VI) =) 3/2Chapter 1. Introduction^ 11There is another problem with the Boltzmann equation description of a pure electronplasma. The Rutherford cross-section diverges in the limit of small relative velocity (g)and small scattering angle (x). A numerical computation must make some approximationto resolve this difficulty. The usual justification is that although the cross-section is largethe effect is small because the deflection angle is small and so the momentum changecaused by the collision should be small. The effects of large impact parameter collisionsare often discarded by cutting off the cross-section at a collision angle corresponding toan impact parameter of a Debye length.The total cross-section describes the total effect of scattering by a single scatteringcentre and is given by^atot(g)^olg, (1) di?.The Rutherford total cross-section diverges and the nature of this divergence can be seenin the following calculation of the change in relative velocity due to the scattering overall angles. This is the integral f a(g, 0)g Ag dli and the component along the relativevelocity g diverges because of the 1/ sin4(x/2) singularity.21f ir^e4^1Jo 4m2g4 sin4(x/2) ( 2g) sin2(x/2) • 27g sin x dx= ^re 1T sin x dxg2m2 io sin2(x/2)= 27re4 r  sin x dxg2m2 io cos x — 12re4= ,^ ln(cos x — 1)r,neThis diverges logarithmically at 0. It is customary12 to 'eliminate' this problem bychanging the lower limit of integration to cut off the integral at the Debye length, whichwill be defined in equation (1.16) in the next section. A small velocity makes the cross-section diverge for all angles (see figure 1.3) so a cut off at some xinin is not clearly theirChapter 1. Introduction^ 12Figure 1.3: The Rutherford cross-section, truncated near where it diverges. g is therelative velocity in reduced units.best choice. Chang has argued for a cut off at a minimum relative velocity to remove thedivergence instead of the usual small angle cut off.22 Solutions to this difficulty which aresuitable for the DBE are discussed in section 2.5.The scattering angle x is related to the impact parameter b and relative speed g bysin (-)2 = [1 +b8(g2)]b2  1 -1/2where b4(g), the impact parameter for 900 scattering is given byb4(g) . 2e2 (1.14)mg2 .(1.13)Chapter 1. Introduction^ 13(Since 900 scattering occurs when sin(x/2) = 1/12-.19) The angle )(min of the Debyelength truncation is given by1.4 Plasma PhysicsA2  I -1/2sin(xmin/2) = {1 + ^D^.bg(g2)(1.15)A common definition for a plasma is that it is an ionized or partially ionized quasineu-tral gas which exhibits 'collective behaviour'. This collective behaviour arises because ofthe long range electrostatic and magnetic forces, especially in a low density `collisionless'plasma in which short range forces are unimportant. This image can be given someprecision by listing several conditions. If a test charge is introduced into a plasma othercharges rearrange around it to shield the long range effect of the charge. The shieldingdistance, or Debye length, depends on the ratio of the thermal energy of the most mobilespecies (usually electrons) to the interaction energy and is given bykT, )1/2,617rne2)where 71, is the electron temperature.This definition suggests another criterion for a plasma: the dimensions of the plasmamust be much larger than AD or this shielding behaviour would be impossible. Further-more there must be many particles available to shield a test charge. The number ofparticles in a `Debye sphere' must be large: 77,17rA3D > 1. Finally, the density of neutralparticles must be low enough so that the oscillation time for the average charged particleis longer than the mean collision time with a neutral particle. If this were not true thegas would be governed by the usual hydrodynamic equations (unless the gas was verylow density and then a kinetic description would be appropriate).16A gas of electrons can act as a plasma. Since there are no positive ions (and noneutral particles either) the gas can't be quasineutral, but in all other respects the gasAD = (1.16)Chapter 1. Introduction^ 14can act as a plasma.23 The electrons will 'shield' out a test charge — the plasma willadjust over a sphere with a Debye length radius to smooth out the charge density andbeyond this length the only effect will be an increased charge of the whole sphere.24,25The most important characteristic, the collective behaviour, is very prominent in theform of oscillations, waves, modes and instabilities.Electron plasmas are composed of only one species so they are simpler and are quitedifferent from standard quasi-neutral plasmas, because, for example, there is only onekind of interaction between particles. The Debye length-scale in the plasma is very usefulbecause it can be used to cut off the divergence in the Rutherford cross-section.1.4.1 Confinement TheoremIn principle a pure electron plasma can be confined indefinitely.23 This contrasts stronglywith the case of a quasi-neutral plasma in which (a) confining the plasma is a hard prob-lem and (b) a confined plasma can not in principle be in thermal equilibrium. This isimportant because the problems studied are about the relaxation to equilibrium. Thisdoes not exclude the dynamics of quasi-neutral plasmas from DBE analysis, but thesystem of equations would need to be more complicated. A schematic of a typical con-finement apparatus is shown in figure 1.4. Axial confinement can be ensured simply bycharging the ends of the confinement tube and grounding the centre; the resulting elec-tric fields can be made strong enough to contain many electrons. An axial magnetic fieldprevents radial losses as can be seen from the following intuitive argument. The totalcanonical angular momentum for a system of N particles in cylindrical coordinates is thesum of its mechanical and vector potential parts:P0 = E (mve, + _c Ao(ri)lii)..7=1(1.17)Chapter 1. Introduction^ 15Magnetic field1.----electron^Gate^ confined plasma^ GategunFigure 1.4: Schematic of basic electron plasma experimental apparatus.For a uniform axial magnetic field, ile(r) = Br/2. The mechanical part is insignificantfor large B field. The N-particle Hamiltonian is invariant under rotations if the physicalsystem is perfectly rotationally symmetric, so the mean square radius is conserved:eB NPo ''-' — E r..^ (1.18)2c j.iIf there were oppositely charged particles this argument would fail because an electron andion could move together radially and still conserve the sum. There are details not coveredby this simple argument: the mechanical part increases as the particle moves radiallythrough the electric field, moving electrons radiate energy and this could contribute to anexpansion of the column, and collisions with neutral particles provide a torque and hencecause expansion of the plasma. O'Nei126 has worked out a detailed confinement theoremfor plasmas with a single sign of charge. In a real system losses occur, but the conservationof angular momentum and energy provide an upper bound on the fraction of electronsthat can escape. These losses are due to processes which transfer angular momentum tothe system, i.e., collisions with neutral particles, asymmetries in the cylindrical field andinteractions with the wal1.23 In practice the plasmas can be confined for about a day.Chapter 1. Introduction^ 161.4.2 Experiments with pure electron plasmasSingle species, nonneutral plasmas have attracted a great deal of attention in recentyears because of the fundamental physics involved and novel laboratory applicationssuch as the study of very cold electrons or positrons or even the possible synthesis ofanti-hydrogen.27-3° Malmberg, Davidson and many others have been working with elec-tron gases for more than two decades.25'31-37 Experimentalists have performed a widerange of experiments on a variety of nonneutral plasmas. These one-species plasmascomposed of, for example, electrons, positrons or antiprotons, can be confined for a longtime (,•-• 105s). Basic phenomena such as equilibrium properties, plasma waves and insta-bilities have been studied. Many of these results are summarized in the proceedings ofa recent conference' and in Davidson's monograph.' Current research is vigorous andnew studies are appearing regular1y.33-43 Good numerical methods are needed to studythese systems theoretically to complement experiments.A schematic of a basic containment apparatus is shown in figure 1.4.38 The majorcomponents are the confinement region, electric fields (the `gates'), magnetic fields andthe injection and release pathways. On a more complicated apparatus there may be along series of gates and containment regions of various sizes. The gates are used to expandand shrink the plasma as a way of cooling or heating it. There are also many probesused to perturb the plasma and to analyse the properties of the injected, confined andreleased plasma. Figure 1.5 shows typical properties of an electron plasma in the variousstages from injection to equilibrium.' At equilibrium, the radial density has become aconstant 107 cm-3, and the rotational frequency LOR is about 0.5 • 107 s-1 over the wholeof the plasma. The plasma radius is 1.5 cm, considerably larger than a typical Debyelength AD = 0.2 cm. More complicated phenomena such as oscillations and waves canarise in these plasmas but will not be discussed here.Chapter 1. Introduction^ 17Figure 1.5: Density profiles and rotation profiles of a column of electron plasma relaxingto equilibrium. (From Malmberg et al. 41.)Chapter 1. Introduction^ 18It is not possible to have detailed control of the radial density function of the injectedplasma since it comes from a source such as a hot wire, but, once a confined plasma hasequilibrated it can then be heated along one axis by compression using a series of electricfield gates. This will produce a bi-Maxwellian distribution which is the subject of thisstudy and much previous research.Experimenters have observed a number of relaxation processes:42• Electron-electron collisions cause the velocity distribution to equilibrate along in-dividual magnetic field lines with a relaxation time on the order of milliseconds.• On a longer time scale, of the order of seconds, electron-electron interactions causetransport of heat and particles across the magnetic field. The temperature of thewhole plasma becomes uniform. This temperature relaxation is shown in figure1.6.41 The particle transport causes the density profile to become constant in theinterior of the plasma. Once the plasma has achieved constant density it rotates asa bulk fluid under the influence of the magnetic field. Malmberg reported that thisdensity relaxation is 5000 times faster than the predictions of Boltzmann theory.43The problems in this thesis do not consider distribution functions with spatialvariation, so this result can't be observed in these calculations. More general dis-crete velocity models have a continuous spatial dependence and the DBE might beadaptable to study this phenomenon.• If the applied fields are sufficiently uniform the electrons can be contained for atleast many hours. However eventually the plasma expands radially, causing a lossof plasma to the walls due to imperfect azimuthal symmetry in the container.This description outlines the general behaviour of the electrons which make up theplasma. Charged particles gyrate around field lines. Thermal motions and collisionsW.--1 1^I11Tilt_.^.^I^.^.^' if^ '0 10^20 30^40^50^100Chapter 1. Introduction^ 191.5T(eV)1..st (msec)Figure 1.6: Experimentally measured temperature relaxation from a bi-Maxwellian dis-tribution. Symbols are experimental data points and the curve is fit to the theoreticalpredictions of Ichimaru and Rosenbluth. (From Malmberg et al. 41.)Chapter 1. Introduction^ 20with other electrons cause transport off magnetic field lines within the constraints ofthe confinement theorem. A nonuniform radial distribution causes shear forces as theelectrons move on the cylinders of constant magnetic field. Finally, asymmetries allowthe electron gas to expand and be lost due to collisions with the container.1.5 Thermalization of a Test Particle in a Heat BathThe first problem studied in this thesis is the relaxation of a test particle in aMaxwellian bath of identical particles. The problem was chosen for its similarities tothe major problem studied in this thesis. The physical situation is very similar: theparticles are electrons and the relaxation is due to Coulomb collisions. However, thisproblem is linear and thus much simpler. The DBE is an unnecessarily complicatedmethod for linear problems. The results from the calculation will be compared withresults from other methods as a means of testing the integrity of the computer code.Two Maxwellian distribution functions, each with its own temperature, are consid-ered. The test particles are so dilute that they cannot interact with each other. The bathelectrons are assumed to be in such a great extent that their Maxwellian is not affectedby the occasional collision with a faster (or slower) test electron. Thus only the effect ofthe bath electrons on the test electrons needs to be described.Shoub has studied the problem of a high energy test particle in a Maxwellian bathand concluded that for this case close encounters are very important. 19,20 Przybylski andLigou note that this problem is especially pronounced in the case of high energy chargedparticles incident on a plasma which has a highly 'forward peaked' cross-section.44 Ligou'smethod separates the cross-section into parts, one singular and one 'regular'. Landesmanand Morel use a Legendre basis expansion for the regular part and a finite difference codefor the troublesome singular part.'Chapter 1. Introduction^ 211.6 Relaxation of Anisotropic DistributionsThe bi-Maxwellian (equation (1.11)) is a common anisotropic distribution function.It has a form similar to the Maxwellian, except that there are two temperatures, oneassociated with a single direction in velocity space and the other associated with theperpendicular components. The study of the relaxation of anisotropic distributions withequation (1.12) is non-linear and much more difficult than the linear test particle problem.These distributions arise in nonneutral plasmas, the solar wind,46,47 terrestrial po-lar and ionospheric plasmas:18-53 In nonneutral plasmas created in the lab precise bi-Maxwellian distributions can easily be generated while in the solar wind they arise natu-rally. In the polar wind 0+ ions can be as much as ten times hotter along one directionthan in the perpendicular plane." The anisotropic distributions in the ionosphere areimportant for such practical considerations as the design of aerospace technology. In thenumerical scheme discussed in this thesis the significance of a bi-Maxwellian distributionis that it is anisotropic and not described well by a few terms in a spectral representationor as a small deviation from a Maxwellian. Far from equilibrium distributions whichcan be adequately modeled with the discrete Boltzmann equation are common in plasmaphysics, especially in fusion research and the study of high energy plasmas produced bylasers." This relaxation problem even arises in certain stellar systems because of thesimilarity between the gravitational and electrostatic potentials.55Many different approaches to the general relaxation problem are possible. Often atwo-species plasma is considered. Reimann and Toepffer56 investigated the dependenceof the temperature relaxation time on the mass ratio of the two species. Rouet andFelix57 studied the evolution of a distinguished subpopulation of electrons and ions in anidentical distribution of unlabeled electrons and ions. McCormack and Williams' studiedthe relaxation of a binary mixture of Maxwell molecules from an initially anisotropicChapter 1. Introduction^ 22velocity distribution. They modeled the temperature relaxation with the macroscopicequations derived from the appropriate moments of the Boltzmann equation. Schameland co-workers59 investigated the Coulomb relaxation due to like particle collisions ofan anisotropic distribution. They assumed that the distribution function was alwaysbi-Maxwellian and derived equations using the kinetic energy moment of the Boltzmannequation to get an expression for the relaxation of the temperature ratio T11/71. AFokker-Planck code using a Legendre polynomial expansion69,61 was used to test theapproximation that the distribution was bi-Maxwellian through the evolution. Theyconcluded that there was good agreement between the two methods for T1/711 > 1 butnot for 71/711 < 1. In this case the bi-Maxwellian model gave a faster relaxation rate.These results are supported by later work by McGowan and Sanderson.62The results of the calculations with anisotropic distribution functions are given inchapter 3.1.7 Other Numerical ApproachesThe description of other methods is divided into two parts. The first describes meth-ods used in kinetic theory or fluid dynamics which are related to the discrete velocitymethods. The second describes other standard methods used in solving plasma physicsproblems.The DBE used in this thesis is a particular example of a discrete velocity method.9The terminology in the literature is somewhat varied, but usually a discrete velocitymethod uses a finite set of velocity vectors but continuous time and space variables.There are two other classes of methods which are related to discrete velocity methods:the discrete simulation Monte Carlo (DSMC) and lattice gas automata and the latticeBoltzmann equation (LGA/LBE). They are similar because all the methods use discreteChapter 1. Introduction^ 23velocities and can be used to solve similar problems. The primary difference is thatthese other methods are phenomenological simulations of microscopic physics instead ofmethods of solving integro-differential equations. Only a tiny fraction of the particlesin a real distribution can be simulated. Any DSMC or LGA simulation will necessarilyhave noise in the computed moments due to small sampling of ensembles. The LBE andDBE avoid this difficulty by computing with probability densities instead of particles.1.7.1 SimulationsThe DSMC refers to a discrete velocity model which simulates particles subject to 'rea-sonable' microscopic rules, such as requiring conservation of mass, momentum and energyin a collision. The method divides naturally into two steps: convective flow and parti-cle collisions. At the beginning of an iteration particles are grouped together in cellsmuch smaller than the mean free path. Possible binary collisions are selected until thepredicted collision rate is reached over a given time interval. The particles which havenot collided are moved a distance v At and the collision step starts again. The methodis easily adaptable to a massively parallel computer architecture by dividing physicalspace according to particle density and assigning each region to a different processor.Goldstein63-65 uses a discrete velocity method which uses integer values for velocities toreduce demand on memory and computational time. Choosing specific discrete velocitiesand discrete time steps means that particles are always located on the nodes of a regularlattice which reduces the amount of memory required and improves the computationalefficiency further.This model does not solve the Navier-Stokes or Boltzmann equations but can 'reduce'to the Navier-Stokes or Boltzmann equations with the correct choices of velocities andcollisions for certain physical situations. The DSMC has also been used by Erwin et al.to study the shock wave problem for helium and neon gases.'Chapter 1. Introduction^ 24Another common discrete method is the Lattice Gas Automaton (LGA).87 In thismodel, time, as well as velocity and space, is discretized. Particles are located at thenodes of a lattice and move to a connected node at the next iteration. Usually onlyone speed is possible; however, some authors include particles with zero speed and morecomplicated lattices to allow for a third speed.88'89 Collisions obey simple microscopicrules designed only to conserve particles and momentum. This description is very dif-ferent from an integration of a partial differential equation, such as the Navier-Stokesequation, yet it can produce very similar and accurate results when compared with moretraditional hydrodynamic methods." This model is easily adaptable to parallel computerarchitectures, and simulations on specialized hardware are very fast. Cellular automataand LGA have been used to model several problems including shear flow between par-allel plates, shock waves, Rayleigh flow, chemically reacting systems," diffusion," flowsin porous media," phase transitions and magnetohydrodynamics.74 The proceedings oftwo recent conferences7" contain many more studies of LGA and the LBE.McNamara and Zanettin have developed a Boltzmann-Lattice Gas hybrid which sim-ulates average populations of particles on a lattice. This method eliminates the noisein macroscopic quantities calculated with LGA; simulating average populations insteadof single particles is equivalent to averaging over an ensemble of initial LGA configura-tions. Unfortunately, it also loses the computational advantage of all computations beingboolean operations. LGA simulations also suffer from a loss of Galilean invariance dueto the discrete lattice and discrete velocities. Chen, Chen and Matthaeus" show thatthe Lattice Boltzmann equation (LBE) eliminates these problems and it can be used tosimulate Navier-Stokes flows correctly if the Reynolds number is not too high." The LBEhas also been used for magnetohydrodynamics,8° and to describe the two phase flow ofmiscible fluids." Benzi, Succi and Vergassola82 have written a review of the lattice Boltz-mann equation methods and Monaco and Preziosi8 gave a general summary of methodsChapter 1. Introduction^ 25for many problems.1.7.2 Related methods for plasma physicsThere are many choices to be made when approaching a plasma physics problem describedby the Boltzmann equation and this is reflected in the array of solution methods. Finitedifference schemes and particle-in-cell (PIC) codes83-88 are common in plasma physics,especially in fusion research.86'87 Another common approach is to expand the distributionfunction in terms of some basis functions. Spectral methods are often used when thematrix representation of the problem's operator is diagonal for a special choice of basisfunctions. For example, Laguerre polynomials are ideal for the linear problem with aMaxwell molecule cross-section discussed in this thesis. In general, the solution of theBoltzmann equation and in particular the evaluation of the collision integral is verydifficult. A number of common approaches are outlined below.For problems involving only Coulomb forces the collision integral can be rewrittenin the Fokker-Planck form as a differential operator. This is done by making a Taylorexpansion in velocity increments due to collisions and truncating the expansion at secondorder. This discounts large changes in velocity resulting from a collision as unimportant.This was done by Chandrasekhar88 and written in the following form which is convenientfor plasma physicists by Rosenbluth, MacDonald and Judd:210^ 1 0L[f] = —^[Tv • A — ---FI, • ( Df)] .^ (1.19)The drift (dynamical friction vector), A, and diffusion (velocity diffusion tensor), D,coefficients can be written using the Rosenbluth potentials g and h:A = , c Oh4--M2 if9V'c 02gD = 2m2 thrthr'Chapter 1. Introduction^ 26g(v) =^dv' f (V)Iv — vi I,h(v) =^dv' f (vI) ^1 Iv — v' Iwhere c = 2re' in A, the Coulomb logarithm is in A = in (AD /bo), the Debye length isAD = (kT/4rnee2)112 and the classical distance of closest approach is bo = e2A.mg2 ,z-e2/3kT. (The approximation is justified since bo only appears in a logarithmic factorand usually in A>> 1.) This was derived using the Coulomb potential and ignoring smallangle scattering. This reformulation is advantageous because it cuts off the singularityin the integral operator due to the Rutherford cross-section by cutting off a above largeimpact parameter collisions. This is the standard starting place for many finite differenceschemes which are common methods in fusion plasma calculations.87Przybylski and Ligou44 used a hybrid Boltzmann-Fokker-Planck equation (BFPE) byseparating a divergent o- into a regular part and a singular part. Their goal was to includethe high energy collisions which are not included in the Boltzmann equation. The regularpart can be treated with the Boltzmann equation part and the singular part is computedusing the Fokker-Planck part of the equation. This means the distribution function ishighly anisotropic and that many more terms in a Legendre expansion will be needed.The BFPE circumvents this problem, providing better results with less effort.Another common simplification of problems studied is to consider distribution func-tions which have only small deviations from an equilibrium Maxwellian. This makes theproblem accessible to perturbation methods. The distribution function can be writtenas a Maxwellian plus an anisotropic term:f (v , t) =^(v) fi(v, 0).^ (1.20)fi(v) is Maxwellian and fi(v, 0) is small compared with f0M(v).89,9° Many importantproblems involve distribution functions which are not close to a Maxwellian and so theyChapter 1. Introduction^ 27cannot be attacked by perturbation methods or a linearization of the appropriate equa-tions.There is another common family of techniques which is very different from the onetaken here: assume that the distribution function has a certain functional form through-out the relaxation. Livi and Marsch' used special distribution functions such as 'kappa'and 'self-similar' distributions which they assumed asymptotically evolved to a Maxwellianwithout changing functional form just by varying one parameter in the function. Theyshowed that the collision operator could be rewritten in a form with an explicit relax-ation rate. This equation with a specified distribution function can be solved to computethe relaxation rate. No detailed information about the distribution function can be in-ferred this way, but often one is only interested in macroscopic properties of the wholesystem. McGowan and Sanderson have shown that far from equilibrium distributionsdo not retain their shape during relaxation.62 It is well known that the collision ratedepends on the detailed form of the distribution function. To estimate the deviations inthe computed collision rates the detailed evolution of the distribution function must becomputed. Schamel et al.59 used this approximation for a bi-Maxwellian distribution.Another standard method in plasma physics for solving the Boltzmann equation isto perform an expansion of the function and operators in appropriate basis functions.For example, a convenient way to represent anisotropic, but azimuthally symmetric,distribution functions is to write the deviations from a Maxwellian in a Legendre (h)function basis (using notation from equation (1.20)):Nfl(V) = fl(V 1 e) = E ii(v)Pi(A)1=0where it = cos 0. Legendre functions are suitable for the Boltzmann collision operatorbecause the operator L is a 'scalar' in the quantum mechanical sense of being invariantunder rotations. Legendre polynomials are the polar (0) factor of the spherical harmonics(1.21)Chapter 1. Introduction^ 28and are eigenvectors of L. For the method to be useful the basis functions should beeigenfunctions of the collision operator and the distribution should be well representedby a few terms. This is the standard method for linear problems, but it can workwell for non-linear operators too. In the latter case one must often resort to numericalcalculations. Using basis functions to solve the nonlinear Boltzmann equation requiresthe calculation of moments of the collision operator and these integrals can be difficultto evaluate accurately to high order.Chapter two describes the DBE in detail. The third chapter gives the details of thecomputations made and discusses the results for both problems. Chapter four gives abrief summary, with conclusions and suggestions for future work.Chapter 2Numerical Methods2.1 IntroductionThis chapter provides the details of the discrete Boltzmann equation discrete ve-locity model. The problem is to solve an adequate approximation to the Boltzmannequation (1.5). The most basic properties of a solution are that the distribution shouldbecome Maxwellian (since there are no external forces), and that particle density, mo-mentum and energy should be conserved. The first is guaranteed by an H theorem' andthe others should be apparent from the description of the model. A full solution will givethe temporal evolution of the distribution function, but the interesting result will be therelaxation of the temperature ratios.The discrete velocity method simplifies the Boltzmann equation by restricting thedomain of definition of the distribution functions to a finite set of velocity vectors. Theintegro-differential equation is reduced to a large set of coupled ordinary differentialequations. The problem is now to find values of Ni(t) ..-÷- f (vi,t) given the initial values.An outline of the algorithm is given in table 2.1, and this chapter describes the detailsof the various steps.2.2 The Discrete Boltzmann EquationThe difficult part in computing the solution to the Boltzmann equation for the prob-lems in this thesis is the evaluation of the collision integral. The time integration is29Set initial parametersr, AtiSet initial distribution functionMaxwellian or bi-MaxwellianiCompute time-saving look-up tablespossible collision outcomesand symmetry tablesITIME ITERATIONIIf the distribution function isnot close enough to equilibriumrepeat iteration^IFor eachr____,.particle zAdd gain term Eki AtiNkNIusing collision tablesICompute loss term AliciNiN.;for each jIAdvance in time using secondorder differencing schemeICompute macroscopic quantitiesCheck conservation of particles,momentum and temperatureIOutput results for this time stepChapter 2. Numerical Methods^ 30Table 2.1: Flowchart for the DBE numerical algorithm.Chapter 2. Numerical Methods^ 31straightforward. The essence of this method is that the domain of the distribution func-tion, the velocity space, is discrete. Specifically the velocities, v, used form a uniformlattice as shown in figure 2.1. The absolute scales are irrelevant and reduced units areused throughout. A discussion of the units accompanies the results in chapter 3. TheDBE approximation to the Boltzmann equation isE A IL!(Nk — NiNi ),^1, p = (203^(2.1)j,k,1=1= E RE AlicfNkNi) — (Eft) NJ/v.] ,^(2.2)j^k,1^ k,1where the Ni are values of the distribution function on the box centered on velocity i,andA151^ l- gij akii = Ivi — 14.Here^is the transition probability for the collision pair with initial velocities vi, vjresulting in vk, v1 velocities. The effect of the cross-section is given by alil, the probabilityof choosing the k,1 outcome from a i, j collision. The resolution of the velocity grid isspecified by, r, the number of speeds on one half-axis.2.3 Numerical IntegrationThe time evolution of the Boltzmann equation for the problems in this thesis should besmooth. A 'predictor-corrector' integration method should be able to follow the evolutionwithout falling into any pathological, e.g., oscillatory or discontinuous, behaviour of f(t).The computational bottleneck in the integration is the evaluation of the collision integral.Although the goal of the discrete velocity approximation is to make the collision term assimple as possible to evaluate the dimensionality suggests that a method which requiresvery few evaluations of the collision terms desirable. This rules out high-order Runge-Kutta or Richardson extrapolation methods.91aNiatChapter 2. Numerical Methods^ 32Figure 2.1: The discrete velocities lie on nodes in the grid. The picture in three dimen-sions is similar.Chapter 2. Numerical Methods^ 33The following is a standard derivation of differencing schemes, accurate to secondorder in time. The discrete Boltzmann equation can be writtend—dt N(t) = L(N) = C(N) — D(N), (2.3)where N(t) = (N1(t),N2(t), ... , Np(t)) is the discrete representation of the distributionfunction, and C(N) is the first set of sums in equation (2.2) and D(N) is the second partof the collision sum. A Taylor expansion of N(t) in t isdNiNi+' = Ni + At dt i12 d2Ni+ (At) dt2 + ORAt)1,i(2.4) where the abbreviation Ni = N(ti) has been used and Ni+1 = N(ti + At). The derivativeterms can be rewritten another way:dNi. Lii_ dLi(2.5)dtd2Nidt2 i^—^dtLi — Li-1(2.6)=Atwhere Li = L(Ni). This gives a second order (in At) explicit scheme for N:N1+1^•^AtArt •^•(3Lz — L'-') + 0[(At)3] (2.7)=^+ —2An implicit scheme can also be worked out. Forward and backward time steps giveAt dN1 ( At)2 d2NNN1"2 = Ni +^+^2 )^+ 0[(Atri,^(2.8)2 dt^2^dt2^i i. N1+1 _ At dN + 1 ( —At) 2 &N + 0[0,031.^(2.9)^2 dt^2^2 1 dt241 41Rearranging these equations yieldsN1+1 = Ni + —At(L1+1 + Li) + 0[(At)3],2(2.10)Chapter 2. Numerical Methods^ 34which is a second order implicit scheme. Computationally the C term appears to handledadequately by the explicit method, but the D term requires the implicit scheme to bestable. Approximate values for Di+1 are computed with an iteration. The estimateDi+l = Di is used to compute provisional values for N. Then better estimates forDi+1 are computed and Nil' is re-evaluated. In the linear test particle problem the samescheme is used but no iteration is needed for the D term.2.4 Distribution Function ConsiderationsThe initial distribution functions used — Maxwellian and bi-Maxwellian — have ro-tational symmetry about the z axis. The effects of collisions can't change this becausethe forces result from a radial potential; there is no azimuthal dependence in the cross-section. This means that a simplification in the collision sum, equation (2.2), can bemade. The discrete velocity grid breaks the rotational symmetry so instead of getting anexact reduction in the dimension of the problem, only a large fraction of the terms arerepeated and need to be calculated only once. The condition f(va) = f(vb) isIVazi = 1Vbziva23J .4_ va2y^2^2= Vbx + Vby(2.11)(2.12)Computationally this is achieved by checking the symmetry and using a recorded valueof the collision sum for each Ni.The representation of test particles with high (or low) energies relative to the bathparticles or the representation of a bi-Maxwellian with a not near 1 causes some problems.Figure 3.1 illustrates this problem — either the distribution function is not zero on theboundary of the grid or the distribution function is non-zero only on a small region of thevelocity grid. A possible solution to this problem is to increase the resolution in an effortto have a better representation of the 'narrow' distribution functions. This is undesirableChapter 2. Numerical Methods^ 35Figure 2.2: Collision sphere (represented by a circle) for particles with velocities v1, v2.Antipodal pairs of dots on the surface indicate the allowed resultant velocities.because the objective is to use as few velocities as possible. Another solution is to changethe velocity units while keeping the energy ratio constant.2.5 Discrete Cross SectionsSince the collision integral is evaluated on a discrete set of velocities, the cross-sectionis also discretized. This is quite straightforward for cross-sections which depend only on g,but difficulties arise for angle dependent interactions. A typical collision relative velocitymight only have 4 or 8 permissible outcomes, resulting in a very coarse sampling of o(x).Figure 2.3 shows a typical situation. The most obvious problem is the zero deflectionChapter 2. Numerical Methods^ 36collisions; the cross-section diverges and must be replaced with a finite approximation.Previous computations with the discrete velocity method have used hard sphere cross-sections. This is easy to implement: a normalization Ek1 CLIO = 1 is used and in thesource term the k,1 dependence of a1,51 is determined by weighting all possible collisionsequally. The relative velocity factor gii is simple to include in the sum on j. The otherisotropic models used (e.g. Maxwell molecules) are a simple change which is accomplishedby changing the velocity weighting factor in the j sum. The anisotropic models areimplemented by changing the weightings in the Eki alici term. A complication arises forthe divergent cross-sections. For every collision there is a head on (x = 0) collision.There is no obviously correct way to weight these terms which have an 'infinite' weight.It is necessary to verify that some properties of the Boltzmann equation are pre-served in the discrete Boltzmann equation approximation. The conservation of somemacroscopic quantities will be checked here, following the calculation done by Gatigno1.3One concern is that the discretization of the cross-section into itqj weights could resultin errors due to coarse sampling between the gain and loss terms, or worse, due to someasymmetry in the approximation. This possibility can be eliminated by checking somevelocity dependent averaged quantities. Consider a function 0(v1) = 0i which may de-pend on v but not on t. Then ( ) = Ei OiNi. If ( 0 ) is conserved then the collision sumshould be zero:o =^(ft i^jkla^= E^+ qANkNl —^—ijkl0i4NiNi)—0kANkNi)—E(ANkNl + (IsiAZNkN, — Ok4NkNIijkl1^/^i •— E vkiAk3iNkNi + ojAiliNkNi — okANkNI2 jkiChapter 2. Numerical Methods^ 37_500400300—Sb200100---1 2 3XFigure 2.3: Angular part of Rutherford cross-section, cr(x) = 1/ sin4(x/2). The solid lineis the analytic form, symbols show typical sample points used in the discrete velocitycomputation.Chapter 2. Numerical Methods^ 381.—2  ,, E. . A tt. ( 0 i + 0.; — 0 k — cb 1 ) Nk NIk,The first step requires microreversibility, illici = A. Subsequent manipulations involvethe interchange of dummy summation variables and the symmetry condition Aill = At!.If these constraints are satisfied then Oi = 1, Oi = vi and cbi = lvi 12 are conserved. Thismeans that particle density, momentum and energy are all conserved.3 Each conserved95i gives a transport equation for a macroscopic quantity.Previous applications of the discrete velocity model to the Boltzmann equation haveonly used the isotropic hard sphere cross-section. For problems involving elastic collisionsof neutral particles this is fine, but it is not satisfactory for collisions in a plasma. TheRutherford cross-section is very different; it is anisotropic and depends on the relativevelocity in a collision. Five cross-sections are used in this thesis to test the expressionof the different interactions within the coarse discretization. A list of the models usedis in Table 2.2. The hard sphere, Maxwell molecule and Rutherford cross-sections arecommonly used and implemented here to show the versatility of the DBE. The isotropicRutherford and simple anisotropic cross-sections are included because they are two factorsof the Rutherford cross-section which are implemented in different ways in the DBE.Modifying the program to use a different cross-section is much easier with the DBE thanwith other standard methods. For example, changing the cross-section in a polynomialexpansion method requires a difficult computation of matrix elements (see section 3.3).This is one of the main advantages of the DBE, and is a consequence of the simplicity ofthe model.Chapter 2. Numerical Methods^ 39Cross-section a kn...,ijHard Spheres 1 1Maxwell Molecules g i(v; — vi) • (vi — vi)Isotropic 'Rutherford' 1/g4 Vivi — v212f^.^\ 2g gAnisotropic 'Rutherford' 1/ sin4(x/2)• g — g • g')Rutherford lie sin4(x/2)21( g . g .. g . gi )Table 2.2: Cross-sections used in the DBE. The normalization constants have been omit-ted for clarity.Chapter 3Calculations and Results3.1 IntroductionIn this chapter the calculations performed and the results of the calculations are de-scribed. The details of the computations and the FORTRAN program used for the linearand the non-linear problems are very similar and so the description of both problemshas been integrated. The DBE is very versatile and a variety of calculations are possiblefor different cross-sections. However, the nature of the five cross-sections affects the use-fulness of the method as shall be seen from the calculations. The DBE is shown to beinefficient in the solution of the linear test particle problem, but this benchmark programserves as a test of the computer program's integrity.Results were generated for five different cross-sections with several different temper-ature ratios for each cross-section. For each case several runs with different time stepsizes and velocity grids were computed to verify that an accurate representation of theevolution of the distribution function was achieved. About 400 time steps were requiredto allow a typical distribution to relax. Each time iteration requires the computation ofmany terms in the collision integral. On our Silicon Graphics 4D/340S machine about180 000 collisions can be evaluated each second. Table 3.1 gives timings for a single timestep. These figures are from versions of the code that assume azimuthal symmetry ofthe distribution function about the z-axis. A listing of the program is provided in theAppendix.40Chapter 3. Calculations and Results^ 41r # of collision terms time for iteration4 362 740 < 2 s5 1 762 500 8s6 6 048 050 32s7 17 344 886 90s8 43 507 108 230s9 96 156 270 9.2 min10 197 351 718 17 min11 379 343 580 37 min12 835 730 228 77 minTable 3.1: Timings for one iteration at various grid sizes. The times are directly re-lated to the number of collision terms; the time to compute collision tables and start-upinformation is not included.3.2 Definitions of computational quantitiesReduced units for time and velocity were used for all the computations. This meansthat time, t, is scaled to a new time, r. This scaling depends on the units of the cross-section used and full details will be given with a description of each cross-section. Thissection describes the various units and computational simplifications made.The first problem is concerned with the relaxation of a test particle distribution, ini-tially Maxwellian at temperature T1(0), owing to elastic collisions with a second speciesof the same mass. This second species is in large extent and the distribution function isMaxwellian at temperature T2. The precise details of the representation of the distribu-tion function and the factors that will influence the computation are now established. Anormalized Maxwellian distribution function for a particle density no isf(v) = no ( m2irkT3/2) e-mv2/2kT (3.1)Chapter 3. Calculations and Results^ 42The dimensions in cgs units are[no] = particles/cm3= 9.1096 • 10-28 gk = 1.3806 • 10-18 erg/°K[T] = °K[v] = cm/s.A check shows that the exponential factor is dimensionless and that the normalizationconstant has dimensions ofg \ 3/2^g^\ 3/2 = (cm\ _3erg)^cm2/s2)^s )The choice of temperature then prescribes the velocity units (v' is the dimensionlesscomputational velocity),v = v'V2kT/m cm/s^ (3.3)v' • 5.5055 • 10/cm/s.^ (3.4)The particle density factor can be factored out of the Boltzmann equation and includedin the redefinition of time.For a normalized equilibrium (Maxwellian) distribution the temperature is defined by3^1^1-2 kT = (-2mv2) = -2 mv2f(v)dv,^ (3.5)and in the discrete velocity approximation is computed as3 =Ni(v! vy2 v!).^ (3.6)The expression for the bi-Maxwellian (BM) anisotropic initial distribution functionused in the second problem is given by2r kTo (0)a2/3(0) )^-3/2fBm^ eriz(i-1-0(0)v1)/2kTH(0) (3 .7)(3.2)Chapter 3. Calculations and Results^ 43where a(0) = To (0)/7'1(0) andT(0) = T1(0) (2a(0) + 13^) .^ (3.8)The temperature ratio at time t is denoted by a(t). It is dimensionless and needs noconversion.During the relaxation the distribution function does not necessarily remain Maxwellian.The evolution of a(t) is used as a measure of the relaxation of the distribution function.The parallel and perpendicular temperatures for the anisotropic distributions used inthe second problem are given by the following moments with the distribution function.Recall that the equipartition theorem indicates that each degree of freedom has -1kT ofenergy, thus the kinetic energy in the perpendicular plane is 2• i kTi.^1 ^rb1 ICTII = (2mvii )2= —1m E Niv!^2^ikTi = ( —1mv2 )2 11. -7,71E Ni(v! + v y2).^2 ^iThis gives an expression to relate the parallel and perpendicular temperatures to thepredicted equilibrium temperature,—3T = T1 + —1 TH2^2T =211 + To3^.(3.9)These definitions illustrate the problem with trying to represent distribution functions(or components of distribution functions) with significantly different temperatures: thesums are only computed over the range of discrete velocities used and not over all velocityspace. If the distribution function differs much from zero for velocities beyond the rangeChapter 3. Calculations and Results^ 44of the velocity grid the computed temperatures will be wrong. The broader distributionfunction can be used to determine the velocity scaling to make it 'fit' the grid. Thisplaces a lower limit on the thinness of the cooler distribution; if it is too steep thetemperature sum may be a poor approximation to the integral. The collision integral issimilarly affected. Figure 3.1 illustrates this difficulty for anisotropic distributions withtemperature ratios of a = 1/5 and 5. For r = 1, the distribution function is very poorlyrepresented — in fact, it is constant regardless of temperature. For r = 2, the 64 Ni takeon four different values. With r = 4 the temperature represented by a Maxwellian isonly 98% of the temperature of a Maxwellian described by a velocity grid with r = 5.For higher resolutions, the numerical temperature deviates by less that 0.1% from theexpected temperature, assuming that the distribution function is close enough to zero onthe boundary of the velocity grid.3.2.1 Cross-sectionsAs noted in the introduction, five different cross-sections (listed in table 2.2) are used withthe DBE. The program is very easily modified for different cross-sections. The changesrequired involve changing a simple weighting function. By contrast, the description of anexpansion method in section 3.3 makes it clear that changing the cross-section can be avery difficult task with the expansion method.The first two cross-sections are important models. The hard-sphere cross-section isoften used to describe non-reactive collisions in neutral fluids. It is the only cross-sectionwhich has been previously used with discrete velocity models. The Maxwell moleculecross-section arises from a potential function, V(r) oc 1/r5, and is important in kinetictheory because some calculations of transport properties simplify when it is used.' Thenext two cross-sections are the isotropic and anisotropic factors of the Rutherford cross-section and are included as stages towards the implementation of the fifth cross-section,(a) a = 1/5(b) a = 5Chapter 3. Calculations and Results^ 45Figure 3.1: Representation of a bi-Maxwellian distribution for two different temperatureratios showing computational resolution, r = 8. x and z are the vz and vz components ofthe reduced velocity. The distribution function has been integrated over vy; the functionplotted is f = f fBm dvy.Chapter 3. Calculations and Results^ 46the full Rutherford cross-section.The hard sphere cross-section isd217(9 9 X) = 7 (3.10)where d is the diameter of the colliding particle. The total cross-section is 47o- = rd2.An examination of the Boltzmann equation determines the scaling of time in thecomputation. Recall that0 f-t-(vi,t) = f [f(v1)f(v12) — f(v1)f(v2)]g cr(g,x) di? dv2.^(3.11)The dimensions on the right hand side are[f2 dv2] = [n(If]= (cm)-6(cm/s)-3[g] = cm/s[c(g, X) dill = cm2and the left hand side has units of [f] = cm-3(cm/s)-3.In the computation, no = 1, the velocity is scaled as described above and the totalcross-section is scaled to 1. A change of variables,^I = noVni/2kT7rd2r^ (3.12)= 5.70. 10-6n0d2T-112 T (3.13)transforms the reduced time in the computation, r, to real time, I. Different cross-sectionshave different normalizations, but they can all be written as a(x, g) = ao • crd(x, g) wheread is the dimensionless cross-section and cro has units of cm2, so the scaling of time forother cross-sections has the same form with Ird2 replaced with different constants.The isotropic Maxwell molecule cross-section is0.(x,g) . d2 ()V0 _..1(3.14)7r^it j gChapter 3. Calculations and Results^ 47where d and Vo are related to the polarizability and are taken from the potential function,^V(r) = Vo ( ) 4 .^ (3.15)rj mFor the DBE calculations I have set 2V0&= 1. The time scaling for this cross-it V 2kTm 4V0d2^2noV0d2 section is adjusted accordingly to t = no 2kT m^T = kT r3.2.2 Sample resultsIn this section some sample results are given. The purpose is to introduce the kinds ofplots drawn and to give a few examples of the calculations done to test the convergenceof the results.Figure 3.2 shows how the temperature relaxation for the linear problem convergeswith increasing velocity resolution. A large temperature ratio, T1(0)/T2 = 4, was usedbut the curves are still nearly indistinguishable except at very short times. Figure 3.3illustrates an example of an anisotropic distribution function at various times during therelaxation. Figure 3.4 shows the relaxation of 711 and TL for several different velocityresolutions for the second problem with the hard sphere cross-section and a(0) = 5.Three different types of plot are used to show the results of the computations. Thefirst is a plot of the temperature relaxation computed from the distribution functions. Forthe linear problem the relaxation of the test particle temperature to the bath temperatureis shown, e.g., figure 3.2. For the relaxation of an anisotropic distribution, two curves,TIM and TH(t), are shown relaxing from initial temperatures to the effective temperatureof the whole distribution function, equation 3.9 and figure 3.4. A second plot illustratesthe relaxation of the anisotropy by plotting a(t) = TH(t)ITL(t) versus time. The thirdplot, a log-linear plot of the temperature deviation from equilibrium versus time whichindicates the degree to which the decay to an equilibrium temperature is exponential, isused for both problems.Chapter 3. Calculations and Results^ 48T^I __--16^8Figure 3.2: Temperature relaxation for the linear problem with the hard spherecross-section at different velocity resolutions, r = 5, 6, and 7. At short times the curvesdo not overlap. The uppermost curve corresponds to r = 5.(a) Or = 00.0151.010.005(c) lir = 1.60.015,0.010.005(b) tIr = 0.80.015e0.010.005(d) tIr = 3.60.015f0.010.005Chapter 3. Calculations and Results^ 49Figure 3.3: The relaxation of an anisotropic distribution function, initially bi-Maxwellianwith a(0) = 5.Chapter 3. Calculations and Results^ 50Figure 3.4: The evolution of 711 and 711 for different velocity grid resolutions, r = 6, 7,and 8, with a(0) = 5, using the hard sphere cross-section.Chapter 3. Calculations and Results^ 513.3 Comparison with expansion methodAs noted in chapter one, the linear problem can be solved by expanding the distri-bution function in terms of a set of basis functions. In this section, I will describe anexpansion in Sonine polynomials (Laguerre polynomials with parameter 1/2) and com-pare the results with data from the previous section. This expansion and the use ofSonine polynomials is discussed in Chapman and Cowling.'First the test particle distribution is written as the product of a Maxwellian dis-tribution at the bath temperature and another term which will be expanded in Soninepolynomials,fi(v,t) = f(o)(v) [1+ 0(v, 0],where3/2 e-mv2/2kT,.4 (v) = 270)^m1-kT2)is the same as the bath distribution. The Boltzmann equation is(3.16)(3.17)a (0)—Oth (v)[1 + 0(v, i)] =Since L[e),f1°)] = 0 and fr)(vi)(fr)(vD[1+ 0(vi,t)iir("12)—.4(vi)[1+ 0(vi,i)].e(v2)) o g dO dv2. (3.18)= f(o) ^this simplifies to^11(0) ao^(fr) ovr) _ op)) g dfl dv2at^J= ANvie(v2)[0(vD - 0(vi)]cr g da dv2.(3.19)(3.20)The distribution functions are isotropic so it is convenient to write the velocity as adimensionless speed with the bath temperature.771y=2 )1/2X (-2kT2(3.21)Chapter 3. Calculations and Results^ 52so that4fi (X, t) =^x e_x22^[1 -I- 0(x,t)]-Nfir-Define an(t) by expanding 0 in Sonine polynomials,00(x, t) = E an(t)Sn(X2)n=1wheren^r(n + 1/2)^x2i S(x2) = E(—l)  .i=o^+ 1/2) (n — i)!z!and the polynomials satisfy the orthogonality relation(3.22)(3.23)(3.24)f: x2e-x2 Sn(x2)Sm(x2) dx = r(n + 3/2)UnM•2n! (3.25)With equation (3.24) the Boltzmann equation (3.20) becomesfr) ct° sn(4)dadnt(t)an(t) f f(0) f2[S,i(x?) — Sn(x?)]crg dt? dv2.^(3.26)n=1^ n=1Multiply this by Sm(xl) and integrate over v1 to reduce the Boltzmann equation to setof coupled equations for an(t),wheredam(t)^"Nm^ = E A,nnan(t)dt n=1(3.27)Nm =^S„,,2 (x)fr)(xi) dvi,^ (3.28)andArnn =^.110)(xi)f2(x2)Sm(4)[Sn(x2) — Sn(4)}cr g dO dv2^(3.29)The Nm can be evaluated using the orthogonality relation, giving2r(m + 3/2) Nm = Vim! (3.30)Chapter 3. Calculations and Results^ 53The matrix elements of the collision operator, A„,n, for a hard sphere cross-section aretaken from a computation by Lindenfeld and Shizgal,"'224kT min("1) min("1P- min("1)-P-a (8 p +I)! —^Ect2 — E^ E 22Pr(n n' — 2s — 2p — q — 1) p (1 )11."9:2+q-12 ).19!p=1^s=0^q=0(3.31)(n — q — s — p)!(ns — q — s — p)!q! 2 k2lwhere d is a constant from the interaction potential, see equation (3.15). The relaxationof the distribution function is studied with the solution of equation (3.27) subject tothe initial conditions, a(0). An initial Maxwellian distribution can be written in thefollowing waym  3/2\fl (V1I 0) = ni 2^7rkTi(0))^e-mv?/2kT1 (0)3/24 ( T2 )^2^2T2/71(0)= ni— — e -xTi (0)(3.32)(3.33)or, using the expansion of equation (3.16),m 3/2ii(vi, 0) = ni^-mvv2kT2E1 + 0(x, OA^(3.34)(27rkT2) e-, n1Nix2e-r2 [1 + 0(X, OM^ (3.35)A comparison of these two expressions shows that0(x, 0) = 73/2e(1 -")x2 — 1,^ (3.36)where 7 = T2/T1 (0). The generating function for the Sonine polynomials,00G(t, x2) = (1- — 03/2 = E0 sn(x2)tn.givesco73/2e^= E an(0)S(x2)n=0(3.37)(3.38)Ann, =Chapter 3. Calculations and Results^ 54(using the fact that So^1 and setting ao = 1.) Define t by (1 — 7) = t/(t — 1), then-y = 1/(1— t) and the left hand sides of equations (3.37, 3.38) are equal. Thus the valuesfor an(0) in the expansion of the initial test particle distribution, are given byan(0) = = (7 — 1)n7 )(3.39)It is convenient to define an = an\FAC and to rewrite equation (3.27) in a symmetricform,^d a m^Amn (3.40)dt^VN„Nman.Define the matrix A by (A)mn = AmniVNn/Vm, with n, m = 1, , N. Now considerhow to solve this N-term approximation to the Boltzmann equation,da— = • a.di (3.41)The matrix A can be diagonalized with eigenvalues Ai and eigenvectors ui. Define U tobe the matrix of eigenvectors, and A = U • a. Thenso thatdiTit = (U-1 • A • U)Aandan =dt(3.42)(3.43)an(t) = a(0)et ,^ (3.44)and the solution is1 Nan(t) = ^ E unkak(o)eAkt.k=0(3.45)The choice of cross-section affects the computation of matrix elements of the Boltz-mann collision operator in the chosen basis. Lindenfeld and Shizga194'92 computed theChapter 3. Calculations and Results^ 55matrix elements for the Boltzmann equation with a hard sphere cross-section and theirresults are used in my calculation. The eigenvalues, An, and eigenvectors, uni, are deter-mined from the numerical diagonalization of A.The time dependent temperature is given by3-2 nkTi(t)^= "rnAr?f(Nri,t) dvi (3.46)= -1-mv?f(°)(vi)(1 + OM) dvi (3.47)CC.4+^an(t)S,i(x2)) x4dx (3.48)eo_rnkT2 Jo3 Ti (t) x3e-22d(x2) - a(t) x2e-'2S„(x2)(-x2) dx] (3.49)2 T2 W7r 2 -1 03 Ti(t) 22 T211(5/2) Niai(t)— (3.50)(t)1 - ai(t) (3.51)T2Only ai(t) is used to compute the temperature of the evolving system, but all theeigenvalues contribute to the time evolution because of the transformation in equa-tion (3.42). The convergence of the temperature evolution was studied by increasingthe number of Sonine polynomials, N, used in the expansion.This method is severely limited by the representation of the initial distribution func-tion by the Sonine polynomial expansion. From equation (3.39), it is clear that if the testparticle is more than twice as hot as the bath particles (-y < 1/2) the expansion diverges.For moderate temperature ratios, e.g., 7 = 2/3, the series converges and only about 5terms are needed to approximate its sum. For smaller ratios, 7 = 1/1.75, 1/1.95 as manyas 40-50 terms are needed. The discrete velocity method can accommodate ratios of upto about -y = 1/5. Figure 3.5(a) shows how the expansion of fi(x, 0) for a specific valueof x converges for temperature ratios near, but greater than, 1/2. For a temperatureratio less than 1/2 the series diverges as shown in figure 3.5(b). Efficient methods forChapter 3. Calculations and Results^ 56(a) Convergent case, 7 = 100/190.Figure 3.5: Partial sums of SN(xo) = EO ai(0)S1(4) for xo = 0.5 andai(0) = ((7 — 1)/7)I. xo is a reduced speed.Chapter 3. Calculations and Results^ 5710^20^30^40^50N(b) Divergent case, 7 = 100/205.Figure 3.5: Partial sums of SN(xo) = EL a1(0)Si(x4) for xo = 0.5 anda(0) = ((V — 1)/7)1. xo is a reduced speed.Chapter 3. Calculations and Results^ 58computing the limit of the partial sums could be used (e.g., Richardson extrapolation95)but the usefulness of the method is still restricted because the ar,(0) coefficients divergefor 7 < 1/2.Figure 3.6 shows four different temperature relaxation curves. The uppermost curvecorresponds to an expansion with 2 Sonine polynomials. The curves for N = 4 and5 are visually indistinguishable and represent a converged solution to the expansionmethod for the linear test particle problem. Figure 3.7 compares the results from theexpansion method and the DBE calculation for several different initial temperature ratios.Figure 3.8 is a log-linear plot of (Ti(t)— T2)/T2 = ai(t) versus t for an initial temperatureratio of 5. The dashed line on the plot is a linear extrapolation from the last few timesteps. Notice that the temperature decay is not purely exponential, but has a changingdecay rate at early times. This is because al (t) evolves under the influence of severaldifferent eigenvalues. Similar behaviour is shown for the linear DBE calculation with thehard sphere cross-section in figure 3.10.This contrasts with the results for the Maxwell molecule cross-section which has apurely exponential decay for both the DBE calculation (figure 3.12) and the Sonineexpansion. (The lines in figure 3.12 are straight for short times. The curvature at longtime is a result of round-off error in the calculation.) The matrix A from equation (3.29)for a Maxwell molecule cross-section is diagonal. The relaxation is determined by thefirst eigenvalue and thus is purely exponential. Shizgal and Hubert92'96 have computedthese eigenvalues. Using 40 Legendre polynomials in the expansion gives Ai = —1.872,so the temperature evolution is given by1ai(t) = _ai(o)e-1.872tVA7 (3.52)Chapter 3. Calculations and Results^ 591^2^3^4t/-rFigure 3.6: Convergence of temperature relaxation with increasing numbers of Soninepolynomials. The initial temperature ratio is T1(0)/T2 = 1.75. The curves correspond toN = 2, 3, 4, 5.3 4ii1-Sonine expansion _-DBE _Chapter 3. Calculations and Results^ 601.50.5E24)4".......••••■4-1E:Figure 3.7: Comparison of temperature relaxation calculated with Sonine polynomialexpansions and DBE for initial temperature ratios T1(0)/T2 = 0.5, 1.5, 1.75. The DBEvelocity resolution is r = 7, and N = 8 Sonine polynomials are used in the expansion.0.010.002E=.1..."..E2°.1..,EI.'0.1Chapter 3. Calculations and Results^ 610^1^2^3^4tirFigure 3.8: Log-linear plot of (T1(t)/T2 — 1) versus t from the Sonine polynomial ex-pansion, with T1(0)/T2 = 5. The dashed line is a linear extrapolation of the straighttail.Chapter 3. Calculations and Results^ 623.4 Computational Results3.4.1 Isotropic Cross-sectionsFigure 3.9 shows the DBE solution to the linear hard sphere problem for various initialtest particle temperatures. The temperature plotted is the ratio T1(t)/T2, and the time isscaled as described in the previous section. Figure 3.10 is a log-linear plot of the differencein temperature between the test particles and the bath. The relaxation rate (the 'slope' ofthe log-linear plot) is clearly not constant, but depends on the instantaneous distributionfunction. The relaxation process is complicated, but a common explanation for thechanges in the relaxation rate is that the energetic tail of the distribution function relaxesfaster than the bulk of the distribution.46,62The relaxation for the linear problem with Maxwell molecule cross-sections, fig-ure 3.11, shows different results. Although the temperature relaxation curves have thesame general appearance as the hard sphere results, the relaxation rate for Maxwellmolecules is nearly constant at early times (figure 3.12). This agrees with the theoreticalprediction from the Sonine polynomial expansion.The results for the DBE calculation of the relaxation of anisotropic temperaturedistributions are in figures 3.13, 3.14, and 3.15. Figure 3.13 shows the relaxation ofa(t) for several different initial anisotropy ratios using the hard sphere cross-section.Figure 3.14 is a similar plot, but for the Maxwell molecule cross-section. A log-linearplot of the deviation of the parallel and perpendicular temperatures from equilibrium fora(0) = 5 for the Maxwell molecule cross-section is shown in figure 3.15. As in the plotfor the linear problem, figure 3.12, the lines are straight at early times.Chapter 3. Calculations and Results^ 63Figure 3.9: Test particle temperature relaxation for the hard sphere cross-section. Avelocity resolution of r = 7 was used. The labels T indicate initial temperature ratios,T1(0)/7' 2 4t/r6 81Chapter 3. Calculations and Results^ 64Figure 3.10: Log-linear plot of (T1(t)/T2 — 1) versus t for hard sphere cross-section fromthe DBE calculation. The labels T indicate initial temperature ratios, T1(0)/T2.Chapter 3. Calculations and Results^ 65Figure 3.11: Test particle temperature relaxation for Maxwell molecule cross-section fromthe DBE calculation. A velocity resolution of r = 6 was used. The labels T indicateinitial temperature ratios, T1(0)/T2.1EI'l.■,■••••••E_T 0. 11...,E-70.010.004Chapter 3. Calculations and Results^ 660^5^1 0t/rFigure 3.12: Log-linear plot of (T1(t)/T2 — 1) versus t for the Maxwell moleculecross-section. The T labels indicate the initial temperature ratio, 711(0)/712.Chapter 3. Calculations and Results^ 67Figure 3.13: Relaxation of temperature anisotropy, a, for the hard sphere cross-section.A velocity resolution of r = 7 was used.Chapter 3. Calculations and Results^ 685432100^2^4^6^8^10^12t/rFigure 3.14: Relaxation of temperature anisotropy, a, for the Maxwell moleculecross-section. A velocity resolution of r = 7 was used.Chapter 3. Calculations and Results^ 69Figure 3.15: Log-linear plot of (Ti(i) — T)/T and (Ti(t) — T)/T versus t with a = 5 forMaxwell molecules.Chapter 3. Calculations and Results^ 703.4.2 'Rutherford-type' cross-sectionsThe Rutherford cross-section ise47(g, x) — ,^, sin-4 (1)n ,v.g-.^2 ' (3.53)It consists of a velocity dependent part, 1/g4, times an angular part, 1/ sin4(x/2). Usefulinformation about the properties of the DBE method can be acquired by consideringdifferent versions of the cross-section. First, as a link to the earlier isotropic cross-sectioncalculations, I have ignored the angular part of the cross-section. Figures 3.16 and 3.17show the temperature relaxation for both problems. The time for the distribution torelax is much longer, in reduced time unit, with the 1/g4 cross-section than with either thehard sphere or Maxwell molecule cross-sections. The smaller relative velocity collisionsare given much more weight than the collisions between particles of larger relative velocity.This has the effect of reducing the energy exchange between the different distributionfunctions and the different velocity components in the anisotropic distribution, resultingin a slower relaxation.Next, the anisotropic factor was included, but the velocity part ignored. The rapidchange in this cross-section near x = 0 makes computation with the DBE more difficult.Special treatment (e.g. a cut off) is needed for head-on collisions (x = 0). For the plotsshown, head-on collisions were ignored. The rapidly varying cross-section results in animperfect conservation of particle density. This problem can be reduced by increasing thevelocity resolution from r = 7 used for the isotropic cross-sections to r = 10, the valueused for figures 3.18, and 3.19. Higher resolutions are possible but the time requiredfor the computation increases roughly as r6, so more accurate computations becomeincreasingly more difficult.Finally, the whole Rutherford cross-section was used in the DBE, using the same cut200 10 30 40543Fe:■......•••••■-I-).......El:210Chapter 3. Calculations and Results^ 71t/rFigure 3.16: Test particle temperature relaxation for the 1/g4 cross-section. A velocityresolution of r = 7 was used. The labels T indicate the initial temperature ratio, T1(0)/T2.Chapter 3. Calculations and Results^ 725432100^10^20^30^40t/TFigure 3.17: Relaxation of temperature anisotropy, a, for the 1/94 cross-section. Avelocity resolution of r = 7 was used.Chapter 3. Calculations and Results^ 735432100^2^4^6^8tirFigure 3.18: Test particle temperature relaxation for the angular part of the Rutherfordcross-section. A velocity resolution of r = 10 was used.Chapter 3. Calculations and Results^ 745432100^2^4^6^8t/rFigure 3.19: Relaxation of temperature anisotropy, a, for the angular part of the Ruther-ford cross-section. A velocity resolution of r = 10 was used.Chapter 3. Calculations and Results^ 75off as in the angular-Rutherford cross-section. At this point the DBE becomes unsatisfac-tory. The combination of the longer relaxation times and the sharply peaked cross-sectionact to make the computation very difficult and, in its present form, unusable.3.5 DiscussionThe DBE was applied to a variety of calculations. It works very well for isotropic cross-sections for both linear and non-linear Boltzmann equation problems. The qualitative(and quantitative) behaviour of the relaxation rates computed with the DBE for hardsphere and Maxwell molecule cross-sections agree well with the expansion method.The relaxation time, defined as the reciprocal of the relaxation rate, is an often usedbut somewhat vague measure of large changes in a distribution function. Two differ-ent measurements are of interest for the calculations in this thesis: the time for energyequipartition between two Maxwellian distributions and the time for an anisotropic dis-tribution to be significantly changed towards an isotropic distribution. Spitzer' definesthe latter time asT312t„ = 0.266 n ln A. (3.54)This is the "self-collision" time for electrons, with T the temperature in Kelvin, n thenumber density per cm3, and A = AD /bo, see section 1.7.2.For the test particle problem in which there are two groups of identical particles,both with Maxwellian distributions, but at different temperatures, Spitzer finds that therelaxation rate depends on the temperature difference from equilibrium:wheredTi (T2 — T1)dt^teg (3.55)377117n2k312 ^( Ti^T2'\3"28(2701/2n244e4ln A 7.7.4 + —rn2teq = (3.56)Chapter 3. Calculations and Results^ 764.84. 10-2= ^ (Ti + T2)3/2n2 In A (3.57)(3.58)In a conventional plasma of protons (H+) and electrons with a general non-equilibriumdistribution, the expected relaxation process is (i) first electron-proton collisions make theelectron distribution function isotropic without changing the energy distribution. Then(ii) electron-electron collisions force the electron distribution function to a Maxwellianand (iii) proton-proton collisions make the proton distribution function Maxwellian. Theelectron relaxation is much faster and the two Maxwellian distributions can have differenttemperatures. Finally, (iv) the distribution function of the whole system is driven to aMaxwellian with one temperature. In a test particle system, where both distributions arefor particles with the same mass, (iii) and (iv) don't occur, but it should be possible toobserve the isotropization of the distribution before the slower relaxation to a Maxwellian.Figure 3.20 shows the experimental relaxation rate for an anisotropic pure electronplasma versus nTe-q3/2 ln(r//b). This illustrates the n, T and magnetic field dependenceof the relaxation rate for the pure electron plasma relaxation experiments.41Inamurol° notes that the discrete velocity method is superior to Monte Carlo methodsin that there is no noise and there are no sampling errors. Unfortunately the computationtime increases roughly as the square of the number of velocities rather than just linearlywith the number of particles, and so if the collision integral is particularly difficult asin the case of varying cross-sections, the resolution required may be computationallyinaccessible.Shoub2° has shown that under highly non-equilibrium conditions, the effects of closeencounter collisions can dominate the effect of distant encounters for the test particleproblem. However, his calculations are for very large temperature ratios, e.g. T1(0)/T2 =10, 20, 30, and these are difficult to describe with a small velocity resolution. His initial100106^107^ 108nT312 In (11/b) (cm-3 — ev-3/2)eqlo'Chapter 3. Calculations and Results^ 77Figure 3.20: Experimental relaxation rate for a pure electron plasma. The symbolsare from experiments with pure electron plasmas and the line is from theory. (FromMalmberg et al. 41.)Chapter 3. Calculations and Results^ 78hot distribution function is a delta function. No tests have been done to show how wellthis problem could be approximated with the DBE. This would be an interesting areafor further exploration.Chapter 4Summary and ConclusionsThe discrete Boltzmann equation has been used in a new application to solve theBoltzmann equation for two kinetic theory problems in plasma physics. The method hasbeen shown to be more flexible than methods which expand the distribution function ina polynomial basis. The cross-section is easy to change because the DBE is a straightfor-ward discretization of the Boltzmann equation. Calculations with the hard sphere andisotropic Maxwell molecule cross-sections work well. Unfortunately the sharply peakedRutherford cross-section is not easy to represent when only a few velocities are used. Thisis because the regular lattice discretization of velocities results in only a few irregularlyspaced scattering angles being included in the collision sum. This difficulty can be partlyovercome by increasing the velocity resolution. Different choices for the velocity latticemight be able to overcome this difficulty, but some of the simplicity of the model wouldbe lost.The test particle problem was solved two different ways for the isotropic hard sphereand Maxwell molecule cross-sections and the agreement between the polynomial expan-sion method and the DBE was good. Further investigations could try to use delta func-tion initial conditions and examine the relaxation of large temperature ratios to try andreproduce the results of Shoub.2°Another interesting extension of the method would be to include spatial variationin the distribution function. Malmberg' observed particle relaxation in a pure electron79Chapter 4. Summary and Conclusions^ 80plasma which was much faster than predicted by standard Boltzmann equation argu-ments. It would be interesting to try to construct a theoretical model to reproduce hisexperimental results.Bibliography[1] J. E. Broadwell. Study of rarefied shear flow by the discrete velocity method. J.Fluid Mech., 19, 401-414 (1964).[2] J. E. Broadwell. Shock structure in a simple discrete veloctiy gas. Phys. Fluids, 7,1243-1247 (1964).[3] R. Gatignol. Kinetic theory for a discrete velocity gas and application to the shockstructure. Phys. Fluids, 18, 153-161 (1975).[4] R. Gatignol. Unsteady Couette flow for a discrete velocity gas. Rarefied Gas Dy-namics, 11, 195-204 (1978).[5] H. Cabannes. Couette flow for a gas with a discrete velocity distribution. J. FluidMech., 76, 273-287 (1976).[6] T. Platowski and R. Illner Discrete velocity models of the Boltzmann equation:a survey on the mathematical aspects of the theory. SIAM Review, 30, 213-255(1988).[7] R. Gatignol and Soubbaramayer, editors. Advances in Kinetic Theory and Contin-uum Mechanics. Springer-Verlag, 1991.[8] Roberto Monaco and Luigi Preziosi. Fluid Dynamic Applications of the DiscreteBoltzmann Equation. World Scientific, 1991.[9] Roberto Monaco, Mirian Pandolfi Bianchi, and Alberto Rossani. Thermodynamicsof a discrete velocity model of a gas mixture with bimolecular reactions. RarefiedGas Dynamics, 18, (1993).[10] T. Inamuro and B. Sturtevant. Numerical study of discrete velocity gases. Phys.Fluids A, 2, 2196-2203 (1990).[11] William J. Merryfield and Bernie D. Shizgal. Discrete velocity model for a single-component planetary exosphere. Rarefied Gas Dynamics, 18, (1993).[12] Fred L. Hinton. Collisional Relaxation in Plasmas, chapter 1.5. Volume 1 of Galeeveand Sudan,' 1983.[13] L. E. Reich!. A Modern Course in Statistical Physics. University of Texas Press,1980.81Bibliography^ 82[14] Stewart Harris. An Introduction to the theory of the Boltzmann Equation. Holt,Rhinehart and Winston, 1971.[15] M. B. Lewis. The Boltzmann and Fokker-Planck equations, chapter V, pages 115—139. In Hochstim,' 1969.[16] Francis F. Chen. Introduction to plasma physics and controlled fusion, volume 1.Plenum, second edition, 1974.[17] I. P. Shakarofsky, T. W. Johnston, and M. P. Bachynski. The particle kinetics ofplasmas. Addison-Wesley, 1966.[18] Herbert Goldstein. Classical Mechanics. Addison-Wesley, second edition, 1980.[19] Edward C. Shoub. Failure of the Fokker-Planck approximation to the Boltzmannintegral for (1/r) potentials. Phys. Fluids, 30, 1340-1352 (1987).[20] Edward C. Shoub. Close encounters in Coulomb and gravitational scattering. I.Relaxation of isotropic like-particle collisions. Astrophys. J., 389, 558-589 (1992).[21] Marshall N. Rosenbluth, William M. MacDonald, and David L. Judd. Fokker-Planckequation for an inverse-square force. Phys. Rev., 107, 1-7 (1957).[22] Yongbin Chang. Integration of the Boltzmann collision term over scattering angles.Phys. Fluids B, 4, 313-318 (1992).[23] Thomas M. O'Neil. Plasmas with a single sign of charge. In Roberson and Drisco11,29pages 1-25.[24] Ronald C. Davidson. An introduction to the physics of nonneutral plasmas. Addison-Wesley, 1990.[25] Ronald C. Davidson. Electrostatic shielding of a test charge in a nonneutral plasma.J. Plas. Phys., 6, 229 (1971).[26] Thomas M. O'Neil. A confinement theorem for nonneutral plasmas. Phys. Fluids,23, 2216-2218 (1980).[27] C. M. Surko, M. Leventhal, A. Passner, and F. J. Wysocki. A positron plasma inthe laboratory — How and why. In Roberson and Driscoll," pages 75-90.[28] William P. Kells. Trapping antiprotons, their cooling, and some experimental uses.In Roberson and Driscoll,' pages 118-127.[29] C. W. Roberson and C. F. Driscoll, editors. Nonneutral Plasma Physics, volume175. American Institute of Physics, 1988.Bibliography^ 83[30] Gerald Gabrielse. Extremely cold antiprotons. Scientific American, 267, 78-89(1992).[31] Ronald C. Davidson and Nicholas A. Kra11. Vlasov equilibria and stability of anelectron gas. Phys. Fluids, 13, 1543-1555 (1970).[32] S. F. Nee, A. W. Trivelpiece, and R. E. Pechacek. Synchrotron radiation from amagnetically confined nonneutral hot electron plasma. Phys. Fluids, 16, 502 (1973).[33] J. H. Malmberg and J. S. deGrassie. Properties of nonneutral plasma. Phys. Rev.Lett., 35, 577-580 (1975).[34] A. J. Theiss, R. A. Mahaffey, and A. W. Trivelpiece. Rigid-rotor equilibria of non-neutral plasmas. Phys. Rev. Lett., 35, 1436 (1975).[35] R. A. Mahaffey, S. A. Goldstein, R. C. Davidson, and A. W. Trivelpiece. Rigidrotation and surface envelopes of nonneutral-plasma columns. Phys. Rev. Lett., 35,1439-1441 (1975).[36] J. H. Malmberg and T. M. O'Neil. Pure electron plasma, liquid and crystal. Phys.Rev. Lett., 39, 1333 (1977).[37] M. H. Douglas and T. M. O'Neil. Transport of a nonneutral electron plasma due toelectron collisions with neutral atoms. Phys. Fluids, 21, 920-925 (1978).[38] B. R. Beck, J. Fajans, and J. H. Malmberg. Measurement of collisional anisotropictemperature relaxation in a strongly magnetized pure electron plasma. Phys. Rev.Lett., 68, 317-320 (1992).[39] Michael E. Glinsky, Thomas M. O'Neil, Marshall N. Rosenbluth, Kenju Tsuruta,and Setsuo Ichimaru. Collisional equipartition rate for a magnetized pure electronplasma. Phys. Fluids B, 4, 1156-1166 (1992).[40] S. N. Bhattacharyya and K. Avinash. Toroidal equilibrium of a non-neutral plasmawith toroidal current, intertia and pressure. J. Plasma Physics, 47, 349-359 (1992).[41] John H. Malmberg, C. F. Driscoll, B. Beck, D. L. Eggleston, J. Fajans, K. Fine, X.-P. Huang, and A. W. Hyatt. Experiments with pure electron plasma. In Robersonand Drisco11,29 pages 28-71.[42] C. F. Driscoll, J. H. Malmberg, and K. S. Fine. Observation of transport to thermalequilibrium in pure electron plasmas. Phys. Rev. Lett., 60, 1290-1293 (1988).[43] A. W. Hyatt, C. F. Driscoll, and J. H. Malmberg. Measurement of the anisotropictemperature relaxation rate in a pure electron plasma. Phys. Rev. Lett., 59, 2975-2978 (1987).Bibliography^ 84[44] K. Przybylski and J. Ligou. Numerical analysis of the Boltzmann equation includingFokker-Planck terms. Nuc. Sci. & Engr., 81, 92-109 (1982).[45] Mark Landesman and J. E. Morel. Angular Fokker-Planck decomposition and rep-resentation techniques. Nuc. Sci. & Engr., 103, 1-11 (1989).[46] S. Livi and E. Marsch. On the collisional relaxation of solar wind distributions.Annales Geophysicae, A4, 333-340 (1986).[47] E. Marsch and S. Livi. Coulomb collision rates for self-similar and kappa distribu-tions. Phys. Fluids, 28, 1379-1386 (1985).[48] A. R. Barakat and R. W. Shunk. Transport equation for multicomponent anisotropicspace plasmas. Plas. Phys., 24, 389-418 (1982).[49] A.R. Barakat and R.W. Schunk. Comparison of Maxwellian and bi-Maxwellianexpansions with Monte Carlo simulations for anisotropic plasmas. J. Phys. D, 15,2189-2203 (1982).[50] A.R. Barakat and R.W. Schunk Comparison of transport equations based onMaxwellian and bi-Maxwellian distributions for anisotropic plasmas. J. Phys. D,15, 1195-1216 (1982).[51] A. R. Barakat and It. W. Shunk. Effect of hot electrons on the polar wind. J.Geophys. Res., 89, 9771 (1984).[52] P. M. Banks and G. J. Lewak. Ion velocity distributions in a partially ionizedplasma. Phys. Fluids, 11, 804 (1968).[53] T. E. Holzer, J. A. Fedder, and P. M. Banks. A comparison of kinetic and hydro-dynamic models of an expanding ion-exosphere. J. Geophys. Res., 76, 2453 (1971).[54] J. C. Kieffer, J. P. Matte, H. Pepin, M. Chaker, Y. Beaudoin, T. W. Johnston,C. Y. Chien, S. Coe, G. Mourou, and J. Dubau. Electron distribution anisotropy inlaser-produced plasmas from x-ray line polarization mesurements. Phys. Rev. Lett.,68, 480 (1992).[55] John M. Retterer. Relaxation with close encounters in stellar systems. Astr. J., 84,370-382 (1979).[56] U. Reimann and C. Toepffer. Temperature relaxation in a strongly coupled plasma.Laser and Particle Beams, 8, 771-779 (1990).[57] J. L. Rouet and M. R. Felix. Relaxation for one-dimensional plasma: Test particlesversus global distribution behaviour. Phys. Fluids B, 3, 1830-1834 (1991).Bibliography^ 85[58] Francis J. McCormack and Anna M. Williams. Anisotropic relaxation in binarymixture of Maxwell molecules. Phys. Fluids, 15, 995-998 (1972).[59] H. Schamel, H. Hamnen, D. F. Diichs, T. E. Stringer, and M. R. O'Brien. Non-linear analysis of Coulomb relaxation of anisotropic distributions. Phys. Fluids B,1, 76-86 (1989).[60] M. R. O'Brien, M. Cox, and D. F. H. Start. Fokker-Planck calculations of RFheating in tokamaks. Comp. Phys. Comm., 40, 123-129 (1986).[61] S. Jorna and L. Wood. Fokker-Planck calculations on relaxation of anisotropicvelocity distributions in plasmas. Phys. Rev. A, 36, 397-399 (1987).[62] A. D. McGowan and J. J. Saunderson. On the relaxation of non-thermal plasmas.J. Plasma Physics, 47, 373-387 (1992).[63] D. Goldstein, B. Sturtevant, and J. E. Broadwell. Investigations of the motion ofdiscrete-velocity gases. In Rarefied Gas Dynamics: Theoretical and ComputationalTechniques, volume 118, pages 100-117,1989.[64] D. Goldstein. Near-continuum applications of a discrete velocity gas model. RarefiedGas Dynamics, 17, 846-853 (1990).[65] D. B. Goldstein. Discrete-velocity collision dynamics for polyatomic molecules.Phys. Fluids A, 4, 1831-1839 (1992).[66] Daniel A. Erwin, Gerald C. Pham-Van-Diep, and E. Philip Muntz. Nonequilibriumgas flows. I: A detailed validation of Monte-Carlo direct simulation for monatomicgases. Phys. Fluids A, 3, 697-705 (1991).[67] U. Frisch, B. Hasslacher, and Y. Pomeau. Lattice-gas automata for the Navier-Stokes equation. Phys. Rev. Lett., 56, 1505-1508 (1986).[68] B. T. Nadiga, J. E. Broadwell, and B. Sturtevant. A study of a multispeed cellularautomaton. Rarefied Gas Dynamics, 16, 155-170 (1989).[69] S. Chen, M. Lee, K. H. Zhao, and G. D. Doolen. A lattice gas model with temper-ature. Physica D, 37, 42 (1989).[70] U. Frisch, D. d'Humieres, B. Hasslacher, P. LaRemand, Y. Pomeau, and J.-P. Rivet.Lattice gas hydrodynamics in two and three dimensions. Complex Systems, 1, 648(1987).[71] Ray Kapral, Anna Lawniczak, and Paul Masiar. Lattice gas automaton approachto nonlinear reaction dynamics. Rarefied Gas Dynamics, 18, (1993).Bibliography^ 86[72] Y. H. Qian, D. d'Humieres, and P. Lallemand. Diffusion simulation with a deter-ministic one-dimensional Lattice-Gas model. J. Stat. Phys., 68, 563-572 (1992).[73] G. A. Kohring. The Cellular Automata approach to simulating fluid flows in porousmedia. Physica A, 186, 97-108 (1992).[74] Shiyi Chen, Daniel 0. Martinez, W. H. Matthaeus, and Hudong Chen. Magnetohy-drodynamics computations with lattice gas automata. J. Stat. Phys., 68, 533-556(1992).[75] 'Lattice gas methods for partial differential equations: Theory, applications andhardware', NATO Advanced Workshop, Sept 1989, editor Gary D. Doolen in PhysicaD, 47, 1-340 (1991).[76] 'Advanced Research Workshop on Lattice Gas Automata', Nice 1991, J. Stat. Phys.,68, 347-672 (1992).[77] Guy R. McNamara and Gianluigi Zanetti. Use of the Boltzmann equation to simu-late lattice-gas automata. Phys. Rev. Lett., 61, 2232-2235 (1988).[78] Hudong Chen, Shiyi Chen, and William H. Matthaeus. Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method. Phys. Rev. A, 45, R5339-5342 (1992).[79] Shiyi Chen, Zheng Wang, Xiaowen Shan, and Gary D. Doolen. Lattice Boltzmanncomputational fluid dynamics in three dimensions. J. Stat. Phys., 68, 379-400(1992).[80] Shiyi Chen, Hudong Chen, Daniel Martinez, and William Matthaeus. Lattice Boltz-mann model for simulation of magnetohydrodynamics. Phys. Rev. Lett., 67, 3776-3779 (1991).[81] Richard Holme and Daniel H. Rothman. Lattice-gas and lattice-Boltzmann modelsof miscible fluids. J. Stat. Phys., 68, 409-429 (1992).[82] R. Benzi, S. Succi, and M. Vergassola. The lattice Boltzmann equation: theory andapplications. Phys. Reports, 222, 145-197 (1992).[83] Alan Mankofsky. Numerical simulation of nonneutral plasmas. In Roberson andDriscoll,' pages 246-269.[84] Richard J. Procassini and Bruce I. Cohen. A comparison of Particle-in-Cell andFokker-Planck methods as applied to the modeling of auxiliary-heated mirror plas-mas. J. Comp. Phys., 102, 39-48 (1992).Bibliography^ 87[85] Richard J. Procassini and Charles K. Birdsall. Particle simulation model of transportin a bounded, Coulomb collisional plasma. Phys. Fluids B, 3, 1876-1891 (1991).[86] J. Killeen, G. D. Kerbel, M. G. McCoy, and A. A. Mirin. Computational Methodsfor Kinetic Models of Magnetically Confined Plasmas. Springer-Verlag, 1986.[87] Charles K. Birdsall and A. Bruce Langdon. Plasma physics via computer simulation.McGraw-Hill, second edition, 1990.[88] S. Chandrasekhar. Stochastic problems in physics and astronomy. Rev. Mod. Phys.,15, 1-89 (1943).[89] Gregory C. Statter. Collisional relaxation in plasmas. Master's thesis, University ofBritish Columbia, Department of Chemistry, 1988.[90] K. El Sayed, T. Ivicht, H. Haug, and L. Banyai. Study of the Coulomb Boltzmannkinetics in a quasi-twodimensional electtron gas by eigenfunction expansions andMonte-Carlo simulation. Z. Phys. B, 86, 345-357 (1992).[91] William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.Numerical Recipes. Cambridge University Press, second edition, 1992.[92] Michael J. Lindenfeld and Bernie D. Shizgal. Matrix elements of the Boltzmanncollision operator for gas mixtures. Chem. Phys., 41, 81-95 (1979).[93] Sydney Chapman and T. G. Cowling. The Mathematical Theory of Non-UniformGases. Cambridge University Press, third edition, 1970.[94] Bernie D. Shizgal and M. Karplus. Nonequilibrium contributions to the rate ofreactions. II. Isolated multicomponent systems. J. Chem. Phys., 54, 4345-4356(1971).[95] Carl M. Bender and Steven A. Orszag. Advanced mathematical methods for scientistsand engineers. McGraw-Hill, 1978.[96] Bernie D. Shizgal and D. Hubert. Nonequilibrium nature of ion distribution func-tions in the high latitude auroral ionosphere. Rarefied Gas Dynamics, 16, 3-22(1989).[97] L. Spitzer. Physics of Fully Ionized Gases. Wiley-Interscience, second edition, 1962.[98] A. A. Galeeve and R. N. Sudan, editors. Handbook of Plasma Physics, volume 1.North-Holland, 1983.[99] A. R. Hochstim, editor. Kinetic processes in gases and plasmas. Academic Press,1969.Appendix ADiscrete Boltzmann Equation Program Listing/* dbe.F -- Discrete Boltzmann Equation solver.copyright (c) 1993 Andrew J. Irwin^*//* This code solves the discrete Boltzmann equation for a spatiallyhomogeneous distribution of particles. Velocity is discretizedaccording to the scheme of Inamuro and Sturtevant.The temporal discretization is accomplished via a second-orderimplicit/explicit ("Adams-Bashforth-Crank-Nicholson") scheme.The particle velocities are measures in units of the mean speed of onevelocity component for a Maxwellian distribution -- i.e. the co-efficient\sqrt(2k/m) = 1 and not its real value.IMPLICIT REAL*8(A-H2O-Z)/* Set velocity resolution at compile time. This is the number of speedsone one half axis. For more than 9 increase NPTAB to, say, 80000 */#define NPU 7#define NPV NPU*define NPW NPU*define NPTAB 10000/* The cross-section to be used is also specified at compile timeby using a -Ddefine on the command line or #defining something in here.This allows unused code to be omitted for certain cross-sections,resulting in some (small?) computational advantage.The preprocessor variables used are the following:ifdef LINEAR then use linear code, else solve the non-linear problemHS, MAX, RUTH select the hard-sphere, Maxwell molecule orRutherford cross-sections. RUTH has further subdivisions,ISO and ANISO, indicating if the ISOtropic or ANISOtropic partsare to be used. If both ISO and ANISO are defined then the fullRutherford cross-section is used*/*define HSPARAMETER (NPZ=1,NPRU=2*NPU-1,NPRV=2*NPV-1,NPRW=2*NPW-1,NPL=900,NPL2=900)DIMENSION FN(NPZ,-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)DIMENSION FNm(NPZ,-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)DIMENSION FNP(-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)DIMENSION RHS(NPZ,-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)88Appendix A. Discrete Boltzmann Equation Program Listing^ 89DIMENSION RHSOLD(NPZ,-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)DIMENSION DSUM(-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)DIMENSION DSUMP(-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)DIMENSION dsumpl(0:NPRU,O:NPRV,0:NPRW)DIMENSION dsump2(0:NPRU,0:NPRV,0:NPRW)DIMENSION U(-NPU+1:NPU),V(-NPV+1:NPV),W(-NPW+1:NPW)REAL*4 TARR(2)DIMENSION USQ(0:NPRU),VSQ(0:NPRV),WSQ(0:NPRW)DIMENSION R(O:NPRU,O:NPRV,0:NPRW),AIJ(0:NPRU,0:NPRV,0:NPRW)DIMENSION FNR(0:NPL,0:8)DIMENSION VEL(0:NPRU,O:NPRV,0:NPRW),PROB(1000)INTEGER NR(0:NPL,0:8),NC(0:NPL,0:8),NCS(0:NPL,0:8),IS(0:8)INTEGER NA(0:NPRU,0:NPRV,0:NPRW),IPAR(0:NPRU,O:NPRV,0:NPRW)INTEGER ICOLL(NPTAB4O:7),IFILL(0:NPL,0:7)DIMENSION R2(0:NPRU,0:NPRV),FNR2(0:NPL2,0:4)INTEGER NR2(0:NPL2,0:4),NC2(0:NPL2,0:4),NCS2(0:NPL2,0:4)INTEGER NA2(0:NPRU,0:NPRV),IPAR2(0:NPRU,O:NPRV),IS2(0:4)INTEGER ISYMM(NPTAB4O:3),IFILL2(0:NPL2,0:3)INTEGER ID0(-NPU+1:NPU,-NPV+1:NPV,-NPW+1:NPW)NZ = NPZNU = NPUNV = NPVNW = NPWNRU= NPRUNRV= NPRVNRW= NPRWwrite(6,*) ' X, Y, Z resolution'write(6,1010)nru+1,nrv+1,nrw+11010 format (316)NL = NPLNL2= NPL2NTAB=NPTABNSQ= NPU*NPUIZ = 1IIMIN = -NU+1IIMAX = NUDO 1 IPR=1,1000PROB(IPR) = 1.0DO/DFLOAT(IPR)^1^CONTINUE/* Output (or optionally input) initial parameters */read(5,*)temp,temp2c1031 format (2e15)TEMP = 2.0D0temp2 = 1.0d0write(6,1030)temp1030^format('Alpha =',g17.9)CCON = 1.0D0UMAX = 3.0D0VMAX = 3.0D0WMAX = 3.0D0WRITE(6,*)' ENTER MAXIMUM U,V,W VELOCITIES'write(6,1011)UMAX,VMAX,WMAXAppendix A. Discrete Boltzmann Equation Program Listing^ 90^1011^FORMAT(3D15.6)NT = 800IWRIT =1WRITE(6,*)' ENTER NUMBER OF TIME STEPS, INTERVAL TO WRITEwrite(6,1001)NT,IWRIT^1001^FORMAT(2I6)1003^FORMAT(D15.6)DELT = 0.01D0DELT2 = 0.5DO*DELTWRITE(6,*)' ENTER TIME STEP DELTwrite(6,1003)DELT/* Compute the discrete velocities in UO, VO, WO */CALL SETV(U,V,W,UMAX,VMAX,WMAX,NU,NV,NW)/* write discrete velocities for plotting, axis labelingwrite(6,*)' xdo 75 I=-nu+1,nuwrite(6,1200)u(i)75^continue1200^format(d13.6,', ')write(6,*) ']'^*//* Compute 3d Collision table */CALL TABLE3D(USQ,VSQ,WSQ,R,AIJ,FNR,VEL,NR,NC,NCS,IS,NA,IPAR,ICOLL,IFILL,UMAX,VMAX,WMAX,NRU,NRV,NRW,NL,NTAB)/* Compute 2d symmetry table */CALL TABLE2D(USQ,VSQ,R2,FNR2,NR2,NC2,NCS2,IS2,NA2,IPAR2,ISYMM,IFILL2,NRU,NRV,NL2,NTAB)/* Set Initial Distribution Function, and compute properties at t=0 */CALL INIT(FN,fnm,U,V,W,TEMP,temp2,NZ,NU,NV,NW)DENS = 0.0D0tpar = 0.0d0tperp = 0.0d0WRITE(6,*)' INITIAL DIST:'DO 4 I=-NU+1,NUDO 4 J=-NV+1,NVDO 4 K=-NW+1,NWDENS = DENS + FN(IZ,I,J,K)tpar = tpar + fn(IZ,i,j,k)*(w(k)*w(k))tperp = tperp + fn(IZ,i,j,k)*(u(i)*u(i)+v(j)*v(j))4^CONTINUE#ifdef LINEARwrite(6,1096)0.0,dens,(tpar+tperp)/3*2#elseifwrite(6,1096)0.0,dens,tpar,tperp/2Iendif/* Write integrated version(s) of the density function for plots */Appendix A. Discrete Boltzmann Equation Program Listing^ 91CALL plotout(FN,O,NU,NV,NW)/* Start of Time iteration */CALL DTIME(TARR)DO 100 IT=1,NTDO 5 IZ=1,NZ^/* IZ can be usd for z spatial variation */ILOOP = 0CALL RESETIDO(IDO,NU,NV,NW)/* zero dsump1,2 used in implicit method */DO 8 I=0,NrUDO 8 3=0,NrVDO 8 K=0,NrWdsumpl(i,j,k)=0.0d0dsump2(i,j,k)=0.0d08^continue1* Compute collision terms */DO 10 I=-NU+1,NUDO 10 J=-NV+1,NVDO 10 K=-NW+1,NWIF(IDO(I,J,K) .EQ. 0)THENDSUM(I,J,K) = 0.0D0weight = 0.0d0iweight = 0DO 20 II=-NU+1,NUDO 20 JJ=-NV+1,NVDO 20 KK=-NW+1,NW/* Compute Relative Velocity (twice cm relative velocity) */IC = II-IJC = JJ-JKC = KK-KICA= IABS(IC)JCA= IABS(JC)KCA= IABS(KC)LC = IC*IC + JC*JC + KC*KCIP = IPAR(ICA,JCA,KCA)CSUM = 0.0D0DSUM1 = 0.0D0IPROB = 0weighta = 0.0d0INDEX = 3*NCS(LC,IP) - 2DO 30 NCOLL=1,NC(LC,IP)IICP = (IC+ICOLL(INDEX,IP))/2 + IIF(IICP .LT. IIMIN .0R. IICP .GT. IIMAX)GOTO 31JJCP = (JC+ICOLL(INDEX+1,IP))/2 + JIF(JJCP .LT. IIMIN .011. JJCP .GT. IIMAX)GOTO 31KKCP = (KC+ICOLL(INDEX+2,IP))/2 + KIF(KKCP .LT. IIMIN .0R. KKCP .GT. IIMAX)GOTO 31Appendix A. Discrete Boltzmann Equation Program Listing^ 92IICM = IICPIF(IICM .LT.JJCM = JJCPIF(JJCM .LT.KKCM = KKCPIF(KKCM .LT.- ICOLL(INDEX,IP)IIMIN .0R. IICM .GT.- ICOLL(INDEX+1,IP)IIMIN .0R. JJCM .GT.- ICOLL(INDEX+2,IP)IIMIN .0R. KKCM .GT.IIMAX)GOTO 31IIMAX)GOTO 31IIMAX)GOTO 31/* CROSS-SECTIONSNote that the factor of g in the A_{ij}-{k1} is included here *//* Hard Spheres */#ifdef HSweightl = 1* vel(ica,jca,kca)#endif/* Maxwell molecules (1/g) */#ifdef MAXweightl = 1/vel(ica,jca,kca) * vel(ica,jca,kca)weight1 = 1if (vel(ica,jca,kca).eq.0) weight1=0/* Rutherford */#ifdef RUTH#ifdef ANISO/* compute relative velocities in real v space not coordinate space */rlx = (u(I)-u(II))rly = (v(j)-v(jj))rlz = (w(k)-w(kk))r2x = (u(iicm)-u(iicp))r2y = (v(jjcm)-v(jjcp))r2z = (w(kkcm)-w(kkcp))weightl = rlx*rlx +rly*rly +rlz*rlz -rlx*r2x -rly*r2y -rlz*r2zif (vel(ica,jca,kca).ne.0) weightl = weightl/vel(ica,jca,kca)if(dabs(weight1).1t.le-3) thenweight 1=0elseweightl = 1/(weightl*weight1)endif#endif#endif/* diagnostic write statementswrite(6,*)i,j,k,weightl,u(ic),v(jc),w(kc)write(6,638)gg,gp,acos(gp/gg),weightl,vel(ica,jca,kca) *//* weight C term with (possibly) angular cross-section */#ifdef LINEARCSUM = CSUM + weightl * 0.5* (* FN(IZ,IICP,JJCP,KKCP)*FNm(IZ,IICM,JJCM,KKCM) +* FNm(IZ,IICP,JJCP,KKCP)*FN(IZ,IICM,JJCM,KKCM) )#elseCSUM = CSUM + FN(IZ,IICP,JJCP,KKCP)*FN(IZ,IICM,JJCM,KKCM)*weight1lendif#endifAppendix A. Discrete Boltzmann Equation Program Listing^ 93weighta = weighta + weightlif(weightl.ne .0)IPROB = IPROB+1^9999^FORMAT(7I6,2D12.4)31^CONTINUEINDEX = INDEX+3ILOOP = ILOOP+130^CONTINUEIF(IPROB .GT. 1000)THENWRITE(6,*)' IPROB TOO BIGSTOPENDIFwweight = vel(ica,jca,kca)#ifdef RUTH#ifdef ISOwweight = vel(ica,jca,kca)**(3)if (wweight.gt .le-5) wweight=1/wweight#endif#endifDSUMP1(ICA,JCA,KCA)=DSUMP1(ICA,JCA,KCA) + wweightif(wweight.ne .0)DSUMP2(ICA,JCA,KCA)=DSUMP2(ICA,JCA,KCA) + 1if (weighta.eqRHS(IZ,I,J,K)#ifdef LINEARDSUM(i,j,k) =.0) weighta=1= RHS(IZ,I,J,K) + CCON*CSUM*wweight/weightaDSUM(i,j,k)+FNm(IZ,II,JJ,KK)*ccon*wweight#else DSUM(i,j,k) = DSUM(i,j,k)+FN(IZ,II,JJ,KK)*ccon*wweight#endif20^CONTINUEWRITE(6,*)fn(iz,i,j,k)*dsum(i,j,k),rhs(iz,i,j,k),c^*^fn(iz,i,j,k)*dsum(i,j,k)-rhs(iz,i,j,k)write(6,1300)i,j,k,FN(IZ,i,j,k)*dsum(i,j,k),rhs(1,i,j,k)c^*^, FN(IZ,i,j,k)*dsum(i,j,k) - rhs(1,i,j,k)1300 format(3I4,3D14.6)/* If there is azimuthal symmetry, keep track of this result to save time */CALL DOSYMM(RHS,IDO,ISYMM,DSUM,NC2,NCS2,IZ,I,J,K,NZ,NU,NV,NW,NL,NL2,NTAB)ENDIF /* end of IDO if */10^CONTINUE/* Use Euler treatment of implicit terms on first time step */IF(IT .EQ. 1)THENDO 35 I=-NU+1,NUDO 35 J=-NV+1,NVDO 35 K=-NW+1,NWRHSOLD(IZ,I,J,K) = RHS(IZ,I,J,K)35^CONTINUEAppendix A. Discrete Boltzmann Equation Program Listing^ 94ENDIF/* Advance solution in time using ABCN with single iteration */*ifdef LINEARDO 40 I=-NU+1,NUDO 40 J=-NV+1,NVDO 40 K=-NW+1,NWFN(1,I,J,K) = ((1.0DO-DELT2*DSUM(I,J,K))*FN(IZ,I,J,K) +DELT2*(3.0DO*RHS(IZ,I,J,K)-RHSOLD(IZ,I,J,K)))/(1.0DO+DELT2*DSUM(I,J,K))40^CONTINUE*else DO 40 I=-NU+1,NUDO 40 J=-NV+1,NVDO 40 K=-NW+1,NWFNP(I,J,K) = ((1.0DO-DELT2*DSUM(I,J,K))*FN(IZ,I,J,K) +DELT2*(3.0DO*RHS(IZ,I,J,K)-RHSOLD(IZ,I,J,1C)))/(1.0DO+DELT2*DSUM(I,J,K))40^CONTINUEDO 50 I=-NU+1,NUDO 50 J=-NV+1,NVDO 50 K=-NW+1,NWDSUMP(I,J,K) = 0.0D0DO 55 II=-NU+1,NUDO 55 JJ=-NV+1,NVDO 55 KK=-NW+1,NWICA= IABS(II-I)JCA= IABS(JJ-J)KCA= IABS(KK-K)if (dsump2(ica,jca,kca).eq.0) dsump2(ica,jca,kca)=1DSUMP(I,J,K)=DSUMP(I,J,K)+CCON*FNP(II,JJ,KK)*DSUMP1(ica,jca,kca)/DSUMP2(ica,jca,kca)DSUMP(I,J,K)=DSUMP(I,J,K)+CCON*VEL(ICA,JCA,KCA)*FNP(II,JJ,KK)55^CONTINUEFN(IZ,I,J,K) = ((1.0DO-DELT2*DSUM(I,J,K))*FN(IZ,I,J,K) +DELT2*(3.0DO*RHS(IZ,I,J,K)-RHSOLD(IZ,I,J,K)))/(1.0DO+DELT2*DSUMP(I,J,K))50^CONTINUE*endif5^CONTINUE/* Compute results from this new distribution function ... */DO 95 IZ=1,NZDENS = 0.0D0tpar = 0.0d0tperp = 0.0d0vx = 0.0d0vy = 0.0d0vz = 0.0d0vv = 0.0d0DO 95 I=-NU+1,NUDO 95 J=-NV+1,NVDO 95 K=-NW+1,NWRHSOLD(IZ,I,J,K) = RHS(IZ,I,J,K)Appendix A. Discrete Boltzmann Equation Program Listing^ 95RHS(IZ,I,J,K) = 0.0D0DENS = DENS + FN(IZ,I,tpar = tpar + fn(IZ,i,tperp = tperp + fn(IZ,vx = vx + fn(IZ,i,j,k)vy = vy + fn(IZ,i,j,k)vz = vz + fn(IZ,i,j,k)vv = vv + fn(IZ,i,j,k)95^CONTINUEJ,K)j,k)*(w(k)*w(k))i,j,k)*(u(i)*u(i)+v(j)*v(j))*u(i)*v(j)*w(k)*(u(i)*u(i)+v(j)*v(j)+w(k)*w(k))TIME = DFLOAT(IT)*DELT/* Write out results */IF (JMOD(IT,IWRIT) .EQ. 0)THENCALL DTIME(TARR)WRITE(6,*)'Num collisions, Iterate, time, timingsWRITE(6,*)ILOOP,IT,TIME,TARR(1),TARR(2)DO 96 I=1,NUDO 96 J=1,NVDO 96 K=-NW+1,NWWRITE(6,1099)I,J,K,FN(1,I,J,K)WRITE(6,1097)I,J,K,DSUM(I,J,K),DSUMP(I,J,K),FNP(I,J,K),FN(1,I,J,K)write(6,*)u(i),v(j),w(k)CONTINUEWRITE(6,1198)DENSformat(' Total particle density=',e19.9)FORMAT(3I3,4e17.9)FORMAT(3I7, e19.9)format(g12.4,3e17.9)WRITE(6,1012)IT,TIMEINEARwrite(6,1096)time,dens,(tpar+tperp)/3*2write(6,1096)time,dens,tpar,tperp/2ENDIF100^CONTINUEEND/* End of the main code. Subroutines follow *//* set vectors of velocities */c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE SETV(U,V,W,UMAX,VMAX,WMAX,NU,NV,NW)c^THIS IMPLEMENTATION IS FOR EVEN MODELSc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++IMPLICIT REAL*8(A-H2O-Z)DIMENSION U(-NU+1:NU),V(-NV+1:NV),W(-NW+1:NW)961198109710991096*ifdef L*else#endif write(6,1095)vx,vy,vx,vv1095^format('<Vx>',g17.9,' <Vy>',g17.9,' <Vz>',g17.9,' <v.v>',g17.B)CALL plotout(FN,it,NU,NV,NW)Appendix A. Discrete Boltzmann Equation Program Listing^ 96UFAC = 1.0DO/DFLOAT(2*NU-1)DO 10 I=-NU+1,NUU(I) = -UMAX + 2.0DO*DFLOAT(I+NU-1)*UFAC*UMAX10^CONTINUEVFAC = 1.0DO/DFLOAT(2*NV-1)DO 20 3=-NV+1,NVV(J) = -VMAX + 2.0DO*DFLOAT(J+NV-1)*VFAC*VMAX20^CONTINUEWFAC = 1.0DO/DFLOAT(2*NW-1)DO 30 K=-NW+1,NWW(K) = -WMAX + 2.0DO*DFLOAT(K+NW-1)*WFAC*WMAX30^CONTINUEENDc++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE INIT(FN,FNM,U,V,W,TEMP,TEMP2,NZ,NU,NV,NW)c++++++++++++++++++++++++++++++++++++++++++++++++++++IMPLICIT REAL*8(A-H2O-Z)DIMENSION U(-NU+1:NU),V(-NV+1:NV),W(-NW+1:NW)DIMENSION FN(NZ,-NU+1:NU,-NV+1:NV,-NW+1:NW)DIMENSION FNm(NZ,-NU+1:NU,-NV+1:NV,-NW+1:NW)c TEMP/TEMP2 is alpha = T\par / T\perpdens=0.0d0densm=0.0d0c various distribution functions; velp, velz are (squares) parallel and vzc my fn is bi-Maxwellian.DO 10 IZ=1,NZDO 10 I=-NU+1,NUDO 10 J=-NV+1,NVDO 10 K=-NW+1,NWVEL = DSQRT(U(I)*U(I) + V(J)*V(J) + W(K)*W(K))VELP= DSQRT(U(I)*U(I) + V(J)*V(J) + (W(K)+1.0D0)**2)VELM= DSQRT(U(I)*U(I) + V(J)*V(J) + (W(K)-1.0D0)**2)velp= U(I)*U(I) + V(J)*V(J)velz= W(K)*W(K)FN(IZ,I,J,K)= DEXP(-TEMP*VEL*VEL)FN(IZ,I,J,K)= DEXP(-TEMP*VEL)FN(IZ,I,J,K)= TEMP/VEL**2FN(IZ,I,J,K)= 1.0D0FN(IZ,I,J,K)= DEXP(-TEMP*VELP*VELP)+DEXP(-TEMP*VELM*VELM)FN(IZ,i,j,k)= DEXP(-TEMP*velz-velp)^old styleFN(IZ,i,j,k)= DEXP(TEMP2*(-TEMP*velp-velz))FNm(IZ,i,j,k)= DEXP(-velp-velz)dens=dens+FN(IZ,i,j,k)densm=densm+FNm(IZ,i,j,k)10^CONTINUEDO 20 IZ=1,NZDO 20 I=-NU+1,NUDO 20 J=-NV+1,NVAppendix A. Discrete Boltzmann Equation Program Listing^ 97DO 20 K=-NW+1,NWFN(IZ,i,j,k)=FN(IZ,i,j,k)/densFNm(IZ,i,j,k)=FNm(IZ,i,j,k)/densmif ((i.eq.nu).and.fn(iz,i,j,k).gt.xmax) xmax=fn(iz,i,j,k)if ((k.eq.nu).and.fn(iz,i,j,k).gt.zmax) zmax=fn(iz,i,j,k)20 continuewrite(6,101)xmax,zmax101^format('max f over x border' ,d15.9,' z border' ,d15.9)ENDc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE DOSYMM(RHS,IDO,ISYMM,DSUM,NC2,NCS2,IZ,I,J,K,NZ,NU,NV,NW,NL,NL2,NTAB)c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++c^SETIDO SEES TO IT THAT THE SYMMETRY REQUIREMENTS ON THE DISTRIBUTIONc FUNCTION ARE TAKEN ADVANTAGE OF. FOR SYSTEMS VARYING IN ONE SPATIALc DIMENSION (HERE Z), THIS MEANS THAT COMPONENTS OF FN HAVING IDENTICALc VELOCITY Z-COMPONENTS AND PERPENDICULAR COMPONENTS ARE IDENTICAL, ANDc SO NEED NOT BE RECOMPUTED. SETTING IDO(I,J,K)=1 ENSURES THAT FN ISc NOT RECOMPUTED UNNECESSARILY FOR THE APPROPRIATE (I,J,K).c^THIS IMPLEMENTATION IS FOR EVEN MODELS; FOR ODD MODELS,c HAVE IS = 2*I, JS = 2*3, ETC.c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++IMPLICIT REAL*8(A-H2O-Z)DIMENSION RHS(NZ,-NU+1:NU,-NV+1:NV,-NW+1:NW)DIMENSION DSUM(-NU+1:NU,-NV+1:NV,-NW+1:NW)INTEGER NC2(0:NL2,0:4),NCS2(0:NL2,0:4)INTEGER ISYMM(NTAB4O:3)INTEGER ID0(-NU+1:NU,-NV+1:NV,-NW+1:NW)IS = 2*I - 1JS = 2*3 - 1LS = IS*IS + JS*JSIP = 3INDEX = 2*NCS2(LS,IP) - 1DO 30 NSYMM=1,NC2(LS,IP)IIS = (ISYMM(INDEX,IP)+1)/2JJS = (ISYMM(INDEX+1,IP)+1)/2RHS(IZ,IIS,JJS,K) = RHS(IZ,I,J,K)DSUM(IIS,JJS,K) = DSUM(I,J,K)IDO(IIS,JJS,K) = 1IIS = (-ISYMM(INDEX,IP)+1)/2RHS(IZ,IIS,JJS,K) = RHS(IZ,I,J,K)DSUM(IIS,JJS,K) = DSUM(I,J,K)WRITE(6,*)IIS,JJS,K,DSUM(I,J,K)IDO(IIS,JJS,K) = 1INDEX = INDEX + 230^CONTINUEENDc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE RESETIDO(IDO,NU,NV,NW)c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++IMPLICIT REAL*8(A-H2O-Z)Appendix A. Discrete Boltzmann Equation Program Listing^ 98INTEGER ID0(-NU+1:NU,-NV+1:NV,-NW+1:NW)DO 10 I=-NU+1,NUDO 10 J=-NV+1,NVDO 10 K=-NW+1,NWIDO(I,J,K) = 0^10^CONTINUEENDc+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE plotout(FN,It,NU,NV,NW)c++++++++++++++++++++++++++ output summed 2D table +++++++++IMPLICIT REAL*8(A-H2O-Z)DIMENSION FN(1,-NU+1:NU,-NV+1:NV,-NW+1:NW)write(6,20)it,nu*2, nw*220^format(' a',i2,' := matrix(',i4,',',i4,', [')DO 10 I=-NU+1,NUDO 10 k=-Nw+1,Nwsum = 0.0d0DO 11 j=-Nv+1,Nvsum = sum + FN(1,i,j,k)c x, y are indistinguishable, so sum on y.11 continuewrite(6,21)sum21^format (g15.7,',')c21^format(e24.16,',')100^format(e17.9)10^CONTINUEwrite(6,22)22^format ('0 ]):')ENDC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE TABLE2D(USQ,VSQ,R2,FNR2,NR2,NC2,NCS2,IS2,NA2,IPAR2,ISYMM,IFILL2,NU,NV,NL2,NTAB)c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++IMPLICIT REAL*8(A-H2O-Z)DIMENSION USQ(0:NU),VSQ(0:NV)DIMENSION R2(0:NU,0:NV)DIMENSION FNR2(0:NL2,0:4)INTEGER NR2(0:NL2,0:4),NC2(0:NL2,0:4),NC52(0:NL2,0:4)INTEGER NA2(0:NU,0:NV),IPAR2(0:NU,0:NV),I52(0:4)INTEGER ISYMM(NTAB4O:3),IFILL2(0:NL2,0:3)CHARACTER*60 FILEOUTPCON = 4.0D0/2.0D0C COMPUTE ALLOWED RELATIVE VELOCITIES R, NUMBER OF ALLOWED COLLISIONSAppendix A. Discrete Boltzmann Equation Program Listing^ 99C^NC FOR EACH (RELATIVE VEL)**2 L AND PARITY IPAR; SET UP INDEXINGC^TABLE NCS FOR COLLISION TABLE ICOLL.C DO 5 I=0,NUUSQ(I) = DFLOAT(I)**25^CONTINUEDO 6 J=0,NVVSQ(J) = DFLOAT(J)**26^CONTINUEC DO 10 I=0,NUDO 10 J=0,NVR2(I,J) = USQ(I) + VSQ(J)NA2(I,J) = 0IF(I .EQ. 0)NA2(I,J)=NA2(I,J)+1IF(J .EQ. 0)NA2(I,J)=NA2(I,J)+1IF(JMOD(I,2)^.EQ. 0 .AND. JMOD(J ,2) .EQ. 0)IPAR2(I,J)=0IF(JMOD(I,2)^.EQ. 0 .AND. JMOD(J ,2) .EQ. 1)IPAR2(I,J)=1IF(JMOD(I,2)^.EQ. 1 .AND. JMODO ,2) .EQ. 0)IPAR2(I,J)=2IF(JMOD(I,2)^.EQ. 1 .AND. JMOD(J ,2) .EQ. 1)IPAR2(I,J)=310 CONTINUEC DO 15 IP=0,4DO 15 L=0,NL2NR2(L,IP) = 0FNR2(L,IP) = 0.0D0IFILL2(L,IP) = 015 CONTINUEC DO 24 IP=0,41S2(IP)^= 024 CONTINUEC DO 25 L=0,NL2FL = DFLOAT(L)DO 20 I=0,NUDO 20 J=0,NVIP = IPAR2(I,J)IF(R2(I,J) .EQ. FL)THENIF(NA2(I,J) .EQ. 0)THENFNR2(L,IP)=FNR2(L,IP)+1.0D0FNR2(L,4) =FNR2(L,4) +1.0D0GOTO 30ELSE IF(NA2(I,J) .EQ. 1)THENFNR2(L,IP)=FNR2(L,IP)+0.5D0FNR2(L,4) =FNR2(L,4) +0.5D0GOTO 30ELSE IF(NA2(I,J) .EQ. 2)THENFNR2(L,IP)=FNR2(L,IP)+0.25D0FNR2(L,4) =FNR2(L,4) +0.25D0ENDIF30^CONTINUEENDIF20^CONTINUEDO 35 IP=0,4NR2(L,IP) = NINT(FNR2(L,IP))NC2(L,IP) = NINT(FNR2(L,IP)*PCON + 0.001D0)NCS2(L,IP) = 1S2(IP)+1IS2(IP) = 1S2(IP) + NC2(L,IP)Appendix A. Discrete Boltzmann Equation Program Listing^ 100C^WRITE(6,*)L,IP,NC2(L,IP),NCS2(L,IP)^35^CONTINUEC^WRITE(6,*)"25^CONTINUECC CONSTRUCT COLLISION TABLEC DO 50 I=0,NUDO 50 3=0,NVIP = IPAR2(I,J)L = NINT(R2(I,J))INDEX = 2*(NCS2(L,IP) + IFILL2(L,IP)) -1IF(NA2(I,J) .EQ. 0)THENISYMM(INDEX,IP) = IISYMM(INDEX+1,IP) = JINDEX = INDEX + 2ISYMM(INDEX,IP)^= IISYMM(INDEX+1,IP) = -JIFILL2(L,IP) = IFILL2(L,IP) + 2ELSE IF(NA2(I,J) .EQ. 1)THENISYMM(INDEX,IP)^= IISYMM(INDEX+1,IP) = JIFILL2(L,IP) = IFILL2(L,IP) + 1ELSE IF(NA2(I,J) .EQ. 2)THENISYMM(INDEX,IP)^= IISYMM(INDEX+1,IP) = JIFILL2(L,IP) = IFILL2(L,IP) + 1ELSEWRITE(6,*)' PROBLEM WITH NA ARRAY'STOPENDIF50^CONTINUECC^WRITE(6,*)' ISYMM:'C DO 55 IL=1,20C^WRITE(6,1009)IL,ISYMM(IL,0),ISYMM(IL,1),ISYMM(IL,2),ISYMM(IL,3)1009^FORMAT(9I6)C 55 CONTINUEC ENDC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SUBROUTINE TABLE3D(USQ,VSQ,WSQ,R,AIJ,FNR,VEL,NR,NC,NCS,IS,*^NA,IPAR,ICOLL,IFILL,UMAX,VMAX,WMAX,NU,NV,NW,NL,NTAB)C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++IMPLICIT REAL*8(A-H2O-Z)DIMENSION USQ(0:NU),VSQ(0:NV),WSQ(0:NW)DIMENSION R(0:NU,0:NV,0:NW),AIJ(0:NU,0:NV,0:NW)DIMENSION FNR(0:NL,0:8)DIMENSION VEL(0:NU,0:NV,0:NW)INTEGER NR(0:NL,0:8),NC(0:NL,0:8),NCS(0:NL,0:8),IS(0:8)INTEGER NA(0:NU,0:NV,0:NW),IPAR(0:NU,0:NV,O:NW)INTEGER ICOLL(NTAB4O:7),IFILL(0:NL,0:7)CHARACTER*60 FILEOUTC PCON = 8.0D0/2.0D0CC COMPUTE ALLOWED RELATIVE VELOCITIES R, NUMBER OF ALLOWED COLLISIONSAppendix A. Discrete Boltzmann Equation Program Listing^ 101C^NC FOR EACH (RELATIVE VEL)**2 L AND PARITY IPAR; SET UP INDEXINGC^TABLE NCS FOR COLLISION TABLE ICOLL.C NUO = (NU+1)/2NVO = (NV+1)/2NWO = (NW+1)/2DO 5 I=0,NUUSQ(I) = DFLOAT(I)**25^CONTINUEDO 6 J=0,NVVSQ(J) = DFLOAT(J)**26^CONTINUEDO 7 K=0,NWWSQ(K) = DFLOAT(K)**27^CONTINUEC DO 10 I=0,NUDO 10 J=0,NVDO 10 K=0,NW**R(I,J,K) = USQ(I) + VSQ(J) + WSQ(K)NA(I,J,K) = 0VEL(I,J,K) = DSQRT(R(I,J,K))*2.0DO*UMAX/DFLOAT(NU)IF(I^.EQ.^0)NA(I,J,K)=NA(I,J,K)+1IF(J^.EQ.^0)NA(I,J,K)=NA(I,J,K)+1IF(K .EQ. 0)NA(I,J,K)=NA(I,J,K)+1IF(JMOD(I,2)^.EQ. 0^.AND.^JMOD(.1,2)^.EQ.^0^.AND.JMOD(K,2)^.EQ. 0)IPAR(I,J,K)=0IF(JMOD(I,2)^.EQ. 0^.AND.^JMOD(.7,2)^.EQ. 0^.AND.JMOD(K,2)^.EQ.^1)IPAR(I,J,K)=1IF(JMOD(I,2)^.EQ. 0^.AND.^JMOD(J,2)^.EQ.^1^.AND.* JMOD(K,2) .EQ.^0)IPAR(I,J,K)=2IF(JMOD(I,2) .EQ. 0 .AND.^JMOD(J,2)^.EQ.^1^.AND.* JMOD(K,2) .EQ.^1)IPAR(I,J,K)=3IF(JMOD(I,2) .EQ. 1 .AND. JMOD(J,2)^.EQ. 0^.AND.* JMOD(K,2) .EQ. 0)IPAR(I,J,K)=4IF(JMOD(I,2) .EQ. 1 .AND.^JMOD(J,2)^.EQ. 0^.AND.* JMOD(K,2) .EQ.^1)IPAR(I,J,K)=5IF(JMOD(I,2) .EQ. 1 .AND.^JMOD(.1,2)^.EQ.^1^.AND.* JMOD(K,2) .EQ. 0)IPAR(I,J,K)=6IF(JMOD(I,2) .EQ. 1 .AND.^JMOD(J,2)^.EQ.^1^.AND.10 *^CONTINUE^ JMOD(K,2) .EQ. 1)IPAR(I,J,K)=7C DO 15 IP=0,8DO 15 L=0,NLNR(L,IP) = 0FNR(L,IP) = 0.0D0IFILL(L,IP) = 015^CONTINUEC DO 24 IP=0,8IS(IP)^= 024^CONTINUEC DO 25 L=0,NLFL = DFLOAT(L)DO 20 I=0,NUDO 20 J=0,NVDO 20 K=0,NWAppendix A. Discrete Boltzmann Equation Program Listing^ 102IP = IPAR(I,J,K)IF(R(I,J,K) .EQ. FL)THENIF(NA(I,J,K) .EQ. 0)THENFNR(L,IP)=FNR(L,IP)+1.0D0FNR(L,8) =FNR(L,8) +1.0D0GOTO 30ELSE IF(NA(I,J,K) .EQ. 1)THENFNR(L,IP)=FNR(L,IP)+0.5D0FNR(L,8) =FNR(L,8) +0.5D0GOTO 30ELSE IF(NA(I,J,K) .EQ. 2)THENFNR(L,IP)=FNR(L,IP)+0.25D0FNR(L,8) =FNR(L,8) +0.25D0GOTO 30ELSE IF(NA(I,J,K) .EQ. 3)THENFNR(L,IP)=FNR(L,IP)+0.125D0FNR(L,8) =FNR(L,8) +0.125D0ENDIF30^CONTINUEENDIF20^CONTINUEDO 35 IP=0,8NR(L,IP) = NINT(FNR(L,IP))NC(L,IP) = NINT(FNR(L,IP)*PCON + 0.001D0)NCS(L,IP) = IS(IP)+1IS(IP) = IS(IP) + NC(L,IP)C^WRITE(6,*)L,IP,NC(L,IP),NCS(L,IP)35^CONTINUEC^WRITE(6,*)"25^CONTINUECC CONSTRUCT COLLISION TABLEC DO 50 I=0,NUDO 50 J=0,NVDO 50 K=0,NWIP = IPAR(I,J,K)L = NINT(R(I,J,K))INDEX = 3*(NCS(L,IP) + IFILL(L,IP)) -2IF(NA(I,J,K) .EQ. 0)THENICOLL(INDEX,IP) = IICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = KINDEX = INDEX + 3ICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = -KINDEX = INDEX + 3ICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = -JICOLL(INDEX+2,IP) = KINDEX = INDEX + 3ICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = -JICOLL(INDEX+2,IP) = -KIFILL(L,IP) = IFILL(L,IP) + 4ELSE IF(NA(I,J,K) .EQ. 1)THENIF(K .EQ. 0)THENICOLL(INDEX,IP)^= IAppendix A. Discrete Boltzmann Equation Program Listing^ 103ICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = 0INDEX = INDEX + 3ICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = -JICOLL(INDEX+2,IP) = 0ENDIFIF(J .EQ. 0)THENICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = 0ICOLL(INDEX+2,IP) = KINDEX = INDEX + 3ICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = 0ICOLL(INDEX+2,IP) = -KENDIFIF(I .EQ. 0)THENICOLL(INDEX,IP)^= 0ICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = KINDEX = INDEX + 3ICOLL(INDEX,IP)^= 0ICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = -KENDIFIFILL(L,IP) = IFILL(L,IP) + 2ELSE IF(NA(I,J,K) .EQ. 2)THENICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = KIFILL(L,IP) = IFILL(L,IP) + 1ELSE IF(NA(I,J,K) .EQ. 3)THENICOLL(INDEX,IP)^= IICOLL(INDEX+1,IP) = JICOLL(INDEX+2,IP) = KIFILL(L,IP) = IFILL(L,IP) + 1ELSEWRITE(6,*)' PROBLEM WITH NA ARRAY'STOPENDIFIF(INDEX .GT. NTAB)THENWRITE(6,*)' TABLE ICOLL DIMENSIONED TOO SMALLSTOPENDIF^50^CONTINUEWRITE(6,*)' ICOLL:'DO 55 IL=1,120WRITE(6,1009)IL,ICOLL(IL,0),ICOLL(IL,1),ICOLL(IL,2),ICOLL(IL,3),ICOLL(IL,4),ICOLL(IL,5),ICOLL(IL,6),ICOLL(IL,7)1009^FORMAT(9I6)C 55 CONTINUEEND


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items