Adams, Mark

DSpace@MIT

Short Characteristic Solution to Neutral and Photon TransportMark L. AdamsMassachusetts Institute of Technology, 167 Albany Street, Cambridge, MA 02139, USA(April 1999)Neutral particle transport is described by an integral form, which will lead to an inherentlystable simple numerical solution, of the Boltzmann equation with BGK collision operator.Our governing equation possesses one part depicting an initial (or boundary) conditionand another describing the causal nature of the. Using the short characteristic method,and retaining all terms to order of the mesh spacing squared, discretization is performedvia a very plane second-degree Taylor polynomial with second-order dierence formulas.Acclimating the numerical result to areas outside neutral particle transport the discretizedsolution is morphed into one describing the transport of photons and compared with resultstaken from the astrophysics community. What this paper brings to the literature is a new,accurate, unembellished nite dierence numerical solution to the transport of neutrals andphotons.I. INTRODUCTIONIf our system contains identical neutral particles, possesses no external forces and we may claim thatthe dominant eect of collisions is to restore local equilibrium, then a suitable governing equation is theBoltzmann equation with BGK collision operator:@f@t+ ~v  ~5f = 1(fM − f): (1) represents a collision interval and fM is the highly nonlinear Maxwellian velocity-distribution function. [1,2]By introducing the total derivative ( ddt =@@t + ~v  ~5) on the left-hand-side (LHS) of Eq. 1, a seemingly rst-order dierential equation is formed:dfdt+1f =1fM: (2)where the parameter t now describes a characteristic path in phase space and time.Integral equations allow for stable nite dierence numerical solutions. Solveing Eq. 2 by use of anexponential antiderivative [exp(R t01 dt0)] leads to the following integral equation:f(t) = e−R t01 dt0f(0) +Z t01e−R tt01 dt00fMdt0: (3)On the right-hand-side (RHS), the rst term is dependent on our initial (or boundary) condition and thesecond term appears to sum the contributions of the Maxwellian in such a manner as to include no sourceeect before its time. The essence of this solution is captured in Duhamel’s principle.In most of this paper we parry the physics and focus on establishing a discretization method. Section IIwill delineate a general short characteristic method that may be used for central and boundary regionsof a system; Section III acclimates the discretization method to the description of photon transport; andSection IV contains thoughts of interest for the reader concerned with applying this research.II. DISCRETIZATIONSeeking a balance between removing addle feelings and beginning to grate the reader with details, thissection is dense with information contained within an exterior shell. Our integral equation will be discretizedusing a second-degree Taylor polynomial with second-order dierence formulas (center for interior regionsand backward for boundaries). Much of the detail is directed toward the central region results, the boundaryregion results are presented for completeness.1A. Interior DiscretizationThe short characteristic solution takes a point of focus and considers the influence of local neighbors. Tobegin, we establish a local coordinate (indexed by ‘i’ in Fig. 1) and write our integral equation in the from:fi =Z titi−11(t)exp−R tit1(s) ds fM(t)dt + fi−1 exp−R titi−11(s) ds: (4)t tt i i+1i-1t ttFIG. 1. Mesh Spacing for Interior Discretization FormulasThen we approximate the Maxwellian as a second-degree Taylor polynomial:fM(t) = fM jt=ti +dfMdtjt=ti (t− ti) +12d2fMdt2jt=ti (t− ti)2 +    (5)and use the following centered second-order dierence formulas for Maxwellian derivatives:dfMdtjt=ti=fM(ti+1)− fM(ti−1)2(4t) −16fM(iv)()(4t)2 (6)d2fMdt2jt=ti=fM(ti+1)− 2fM(ti) + fM(ti−1)(4t)2 −112fM(iv)()(4t)2 (7)After placing Eqs. 6 and 7 into Eq. 5, inserting the latter into Eq. 4 and simplifying notation by lettingfM jt=ti= fMi we may write:fi =Z titi−11(t)exp−Rtit1(s) ds[fMi + ffMi+1 − fMi−124 t g(t− ti) +f12fMi+1 − 2fMi + fMi−1(4t)2 g(t− ti)2]dt + fi−1 exp−R titi−11(s) ds (8)The RHS of this equation consists of four separate integrals. Each term will be investigated apart from therest and collected at the end of this section, thus obtaining our desired solution.I) Term OneThis term may be easily solved after a change of variables. To perform the change of variables wedene z and a few properties:z  −Z tit1(s)ds (9)dzdt=ddtZ tti1(s)ds =1(t)(10)zi−1  −Z titi−11(s)ds (11)2Then the rst integral clearly becomes:=) fMiZ 0zi−1expz dz = fMi (1− expzi−1) (12)(While investigating individual terms we denote a step of evolution by \=)", equations that representa denition or some other type of simplication appear with no arrow.)II) Term TwoAfter the following change of variables:=) −ffMi+1 − fMi−124 t gZ titi−11(t)exp−R tit1(s) ds(ti − t)dt= −ffMi+1 − fMi−124 t gZ 0zi−1(ti − t) expz dz; (13)and then integrating by parts we may write:=) −ffMi+1 − fMi−12(4t) g[−4 t expzi−1 +Z 0zi−1(t) expz dz]: (14)The collision interval is treated in a fashion identical to the Maxwellian velocity-distribution function;it is expanded as a second-degree Taylor polynomial with centered second-order dierence formulas forthe derivatives. Expanding the collision interval in Eq. 14 allows us to write:=)−ffMi+1 − fMi−12(4t) g[−4 t expzi−1 +i(1− expzi−1)−fi+1 − i−12(4t) gZ 0zi−1(ti − t) expz dz + fi+1 − 2i + i−12(4t)2 gZ 0zi−1(ti − t)2 expz dz] (15)To nish the discretization, we need to address the last term on the RHS and ensure that our order ofaccuracy is retained. Consistent with the series expansions used thus far, we choose to drop terms oforder (ti − t)3. Working on this last term, we use integrating by parts and our collision interval seriesexpansion to arrive at the following equation:Z 0zi−1(ti − t)2 expz dz =−(4t)2 expzi−1 +Z 0zi−12(ti − t) expz dtdzdz =−(4t)2 expzi−1 +2iR 0zi−1(ti − t) expz dz1 + 2f i+1−i−124t g(16)And then we need to address the penultimate term on the RHS of Eq. 15 by useing integration byparts. Instead of dropping terms we substitute the result of Eq. 16 where necessary. So the fullydiscretized from of this integral is found to be:Z 0zi−1(ti − t) expz dz =− expzi−1(i +4t + [f i+1−2i+i−12(4t)2 g1+2f i+1−i−124t g](4t)2) + i1 + f i+1−i−124t g − 2i[f i+1−2i+i−12(4t)2 g1+2f i+1−i−124t g](17)3(In Term Three this equation will be substituted back into Eq. 16.) To nish we merely note that theabove integral is the solution of Term Two without the appropriate global factor in front. Adding thisfactor leads to the nal discretized second integral term:=) −ffMi+1 − fMi−12(4t) g[− expzi−1(i +4t + [f i+1−2i+i−12(4t)2 g1+2f i+1−i−124t g](4t)2) + i1 + f i+1−i−124t g − 2i[f i+1−2i+i−12(4t)2 g1+2f i+1−i−124t g]] (18)III) Term ThreeAll the ingredients necessary to solve this integral are contained within the previous terms. Here wedene the following terms to simplify future notation:M1  i+1 − 2i + i−1 (19)M2  i+1 − i−1 +4t; (20)then use these denitions to clean up the following messy fraction:[f i+1−2i+i−12(4t)2 g1 + 2f i+1−i−124t g] =M12(4t)M2 (21)Z 0zi−1(ti − t) expz dz =−2(4t) expzi−1(i +4t + 4tM12M2 ) + 2(4t)i4t + M2− 2iM1M2(22)and nally work the new denitions into Eq. 16:Z 0zi−1(ti − t)2 expz dz = −(4t)3 expzi−1M2+−4(4t)2i expzi−1(i +4t + (4t)M12M2 ) + (2(4t)i)2M22 + (4t)M2− 2iM1 (23)Eq. 23 is the essence of Term Three. Adding the necessary global factor to this integral yields thefollowing discretized term:=) 12(fMi+1 − 2fMi + fMi−1(4t)2 )f−(4t)3 expzi−1M2+−4(4t)2i expzi−1(i +4t + (4t)M12M2 ) + (2(4t)i)2M22 + (4t)M2− 2iM1 g (24)IV) Term FourThe fourth term has already been discretized via our denition of zi−1:=) expzi−1 fi−1 (25)Before collecting the terms above it is important to note that zi−1 itself contains an unevaluated integral.Using our collision interval approximation this term may be written:zi−1 = 4t( 112i+1 −23i− 512i−1) (26)Collecting the terms in Eq. 8 ( Eqs. 12, 18 (with slight modications), 24 and 25) leads to our discretizedintegral Boltzmann equation solution:4fi = fMi (1− expzi−1)−ffMi+1 − fMi−12(4t) g[− expzi−1(i +4t + M12(4t)M2 (4t)2) + i1 + f i+1−i−124t g − 2i M12(4t)M2]+12(fMi+1 − 2fMi + fMi−1(4t)2 )f−(4t)3 expzi−1M2+−4(4t)2i expzi−1(i +4t + (4t)M12M2 ) + (2(4t)i)2M22 + (4t)M2− 2iM1 g+ expzi−1 fi−1; (27)where M1, M2 and zi−1 are dened by Eqs. 19, 20 and 26, respectively.B. Boundary DiscretizationImagine a gas contained within a solid box. For this system our central discretization formula will beuseful throughout the interior of the box, however when considering a trajectory that is incident on a sideof the box our central formula becomes invalid. To see this, consider a trajectory whose location (ti) is onthe solid-gas interface (side of the box) and directed into the solid. The value of fi depends on knowledge oflocations inside the box (ti+1) which are assumed to be in local equilibrium and described by a Maxwellianvelocity-distribution. It is clear that the ordered lattice of the box is not represented by a Maxwellian. (A lessinteresting argument could be based on the improbable occurrence of a system with a solid box containinga gas of the same species.) Thus, for systems such as these, we have developed a backward nite-dierencediscretized solution which only considers locations within the gas | see Fig. 2 for a description of the newmesh.t i-3t ttt i-2 t i-1 t itFIG. 2. Mesh Spacing for Boundary Discretization FormulasTo arrive at the boundary discretization formula our technique of the previous subsection is still valid; newdenitions and those that have been changed are now presented before the boundary discretized solution.The central second-order dierence formulas are replaced by the following backward second-order dierenceformulas:dgdtjt=ti=gi−2 − 4gi−1 + 3gi2(4t) +13g(iii)()(4t)2 (28)d2gdt2jt=ti=−gi−3 + 4gi−2 − 5gi−1 + 2gi(4t)2 +1112g(iv)()(4t)2; (29)where g ! ffM; g. New M denitions are used:M3 = −i−3 + 4i−2 − 5i−1 + 2i (30)M4 = 4t + i−2 − 4i−1 + 3i (31)M12(4t)M2 −!M32(4t)M4 (32)and zi−1 takes the form:5zi−1 = −4t12 (7i+2i−1+5i−2− 2i−3) (33)Aware of these changes, one may derive the following boundary discretized integral Boltzmann equationsolution:fi = fMi (1− expzi−1)−(fMi−2 − 4fMi−1 + 3fMi )[expzi−1(i +4t + (4t)M32M4 ) + i4t + M4− 2iM3M4]+12(4t)2 (−fMi−3 + 4fMi−2 − 5fMi−1 + 2fMi )[−(4t)3 expzi−1M4+−4(4t)2i expzi−1(i +4t + (4t)M32M4 ) + (2(4t)i)2M42 + (4t)M4− 2iM3 ]+ expzi−1 fi−1 (34)III. DISCRETIZATION OF THE RADIATION TRANSFER EQUATIONTaking advantage of similarities between the Boltzmann equation with BGK collision operator and theradiation transfer equation, we morph our neutral transport numerical solution into a photon transportnumerical solution. Although this makes it easy to understand the procedure for discretization of photons,comparison with existing astrophysics literature [3{5] is very dicult. Confronting this quandary, an eortis made to allow future comparisons between the solution developed here and those exist in astrophysicsliterature.A. Radiation Transfer EquationPhoton transport may be characterized by an equation of the form:dfRdt= (DfRDt) jcoll; (35)where fR represents the photon distribution function. Most radiation transfer literature discusses the evo-lution of specic intensity, dened by the relation I(~r; n^; ; t)  chfR(~r; n^; ; t), rather than the photondistribution function; the specic intensity is more useful in descriptions of energy transport. In the casewhere the direction of observation and radiation propagation are along the same line, our governing radiationtransfer equation takes the form:dIdt+ I = S: (36)By comparison with Eq. 2 it is easy to see that the dominant inhomogeneity is now the source function (S),where to zero-order this term is behaves like the Planck distribution function.B. DiscretizationTo arrive at a numerical radiation transfer equation, and hence solve for the transport of photons, aquick variable re-identication is performed. The notation gets a little confusing here since in the previoussection we used  to represent the collision interval and here we dene the optical depth () as d = −dz.We change notation by letting: f −! I, 1 −! , fM −! S and zi−1 −! − 4 i−1. Placing the newnotation into Eq. 27 leads to the following central numerical radiation transfer equation solution (Findingthe boundary region discretized solution is trivial.):6Ii = Si(1− exp−4i−1)−fSi+1 − Si−12(4t) g[− exp−4i−1(( 1)i +4t + M12(4t)M2 (4t)2) + ( 1 )i1 + f (1 )i+1−( 1 )i−124t g − 2( 1)i M12(4t)M2]+12(Si+1 − 2Si + Si−1(4t)2 )f−(4t)3 exp−4i−1M2+−4(4t)2( 1)i exp−4i−1(( 1)i +4t + (4t)M12M2 ) + (2(4t)( 1 )i)2M22 + (4t)M2− 2( 1)iM1g+ exp−4i−1 Ii−1 (37)Although the equation above allows for one to understand how to acclimate our neutral particle transportsolution to photon transport it is not useful for comparison with the literature. In astrophysics the shortcharacteristic solution of radiation transfer will typically appear in the form:Ii = Ii−1 exp−4i−1 +Si−1 + Si + γSi+1; (38)After some algebra, Eq. 37 takes the above form with the following coecient:CD = 42i−12i+1 − 2ii−1i+1(i−1 + i+1) +2i (i−1 − i+1 + (4t)i−1i+1)(i−1 − i+1 + 2(4t)i−1i+1) (39) CD = i−1i+1f−2(exp−4i−1 −1)i−1i+1+(4t) exp−4i−1 2i [i+1 − 2i−1(1 + (4t)i+1)]+i[(exp−4i−1 −1)i+1 + i−1(1 − exp−4i−1 +(4t)i+1 − 2(4t) exp−4i−1 i+1)]g (40)  CD = ifi(i−1 − i+1 + 2(4t)i−1i+1)[(exp−4i−1 −1)i+1 + i−1(1 − exp−4i−1 +(4t)i+1)] + 2i−1i+1[(exp−4i−1 −1)i+1 + i−1(−1 + exp−4i−1 +2(4t) exp−4i−1 i+1)]g (41)γ  CD = i−1i+1fii−1(−1 + exp−4i−1 +(4t) exp−4i−1 i)−i+1[2(exp−4i−1 −1)i−1 + i(−1 + exp−4i−1 +(4t)i−1 + 2(4t) exp−4i−1 i−1)]g (42)where CD is just a common denominator. No further comparisons will be made at this time.IV. DISCUSSIONThis paper delineates a procedure for the discretization of general transport equations in a manner ductileenough to be applied to both neutral and photon transport. The procedure is based on the short characteristicmethod and uses a simple second-degree Taylor polynomial with second-order dierence formulas. It thissection we discuss some details of the neutral solution as well as the photon solution.A. Neutral TransportTo begin, it is instructive to look at the limits, both collisional (4t ! 1) and collisionless (4t ! 0), ofthe central discretization solution for neutrals. In the collisional limit we nd:fi = fMi [−2(4t)2 + 1] +     fMi ; (43)and in the collisionless we nd:fi = fi−1 + (fMi − fMi−1)−4t(fMi+1 − 2fMi + fMi−1): (44)These compare very well with the known rst-order limits of:7fi = fMi (45)fi = fi−1 + (fMi − fMi−1) (46)for the collisional and collisionless limits, respectively.Another comment, one that is probably obvious to most readers, is that our solution has the ability tohandle multi-species systems. Nothing prevents us from dening eective variables and updating their valuebetween iterations or applying a linearization technique as will be discussed in the next subsection.B. Radiation TransferIn the introduction we mentioned that the Boltzmann equation with BGK collision operator is valid whenthe dominant eect of collisions is to restore local equilibrium. Due to photon scattering, a process denedas a photon-matter interaction with the net eect being a change in photon direction, radiation transfer isfrequently non-local. In this case the short characteristic method breaks down in the sense that convergencemay take an innite amount of time. Many techniques have been developed to solve this problem and theinterested reader is referred to the literature. [6,7]ACKNOWLEDGMENTSThe author would like to thank S.I. Krasheninnikov (MIT) for his discussions on neutral particle transportand H.A. Scott (LLNL) for invaluable assistance in understanding radiation transfer.[1] S. Chapman and T.G. Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge MathematicalLibrary Series 1990. Cambridge University Press, Cambridge, 3rd. edition, 1995.[2] Carlo Cercignani. Mathematical Methods in Kinetic Theory. Plenum Press, New York, 2nd edition, 1990.[3] G.L. Olson, L.H. Auer, and J.R. Buchler. A rapidly convergent iterative solution of the non-LTE line radiationtransfer problem. Journal of Quantitative Spectroscopy and Radiation Transfer, 35(6):431, 1986.[4] G.L. Olson and P.B. Kunasz. Short characteristic solution of the non-LTE line transfer problem by operatorperturbation | I. the one-dimensional plane slab. Journal of Quantitative Spectroscopy and Radiation Transfer,38(5):325, 1987.[5] P.B. Kunasz and L.H. Auer. Short characteristic integration of radiative transfer problems: Formal solution intwo-dimensional slabs. Journal of Quantitative Spectroscopy and Radiation Transfer, 39(1):67, 1988.[6] Dimitri Mihalas. Stellar Atmospheres. W.H. Greeman and Company, San Francisco, 2nd edition, 1978.[7] D.C. Eder and H.A. Scott. The calculation of line transfer in expanding media. Journal of Quantitative Spec-troscopy and Radiation Transfer, 45(4):189, 1991.8

Short Characteristic Solution to Neutral and Photon Transport

http://dspace.mit.edu/bitstream/1721.1/93721/1/99rr007_full.pdf

Short Characteristic Solution to Neutral and Photon Transport

Abstract

Similar works

Full text

Available Versions

DSpace@MIT