Numerical solution of isotropic relaxation problems by the method of Maxwellian expansion

U.S.S.R. Comput.Maths.Math.Phys Printed in Great Britain

.,Vo1.28,No.3,pp.187-198,988

0041-5553/88 $1o.00+0.oo 011989 Pergamon Press plc

NUMERICAL SOLUTION OF ISOTROPIC RELAXATION PROBLEMS BY THE METHOD OF MAXWELLIAN EXPANSION* I.N. KOLYSHKIN,

A.YA. ENDER and I.A. ENDER

A numerical analytic method of solving the Boltzmann equation based on Maxwell expanding the distribution function and the collision integral is perfected. Two test problems are solved by this method: the temperature relaxation of pseudo-Maxwellian molecules on a Maxwellian background, and the non-linear relaxation of the sum of two Maxwellian distributions. The effectiveness of the numerical scheme is demonstrated. Different varieties of non-monotonic relaxation are investigated. In /l, 2/ an integral transformation of the non-linear Boltsmann equation (BE) was in proposed based on expanding the distribution function (DF) and the collision integral Maxwell distributions ((a-u) and a-representation of the BE). A numerical scheme for solving the BE in the a-representation was developed there. In the example of the solution of the temperature relaxation problem for the hard sphere model, the high accuracy of the method with small expenditures of computer time was demonstrated, Interest in the isotropic relaxation problem is increasing. A study of the physico-chemical processes in gases and in a plasma demands a knowledge of the distribution function for particles whose velocity considerably exceeds the thermal velocity. Classical methods of solving the BE do not allow us to solve such a problem, so information on the behaviour of the DF can be obtained only from a numerical solution. The problem is SO complicated in the general case that it is best to carry out an investigation of the accuracy of the numerical scheme for distribution functions that are isotropicinterms of the velocities. The possibilities of a numerical solution of the BE in the isotropic case for Maxwellian molecules by the statistical modelling method were investigated in detail in /3/. This method does not guarantee sufficient accuracy of the distribution function in the high-velocity region, since even with a large number of operative molecules, only a small portion of them will fall in this region. The behaviour of the distribution function tails is closely connected with their high-order moments. A detailed analysis of the moments in /3/ has shown that the error in calculating the high-order moments increases rapidly with time. A solution of the relaxation problem is constructed in /4/ with the help of the Fouriertransform of the BE. Such a transformation leads to simplifications only for Maxwell-molecule a-representation, type models. In the the structure of the equation does not depend on the interaction potential, and the kernel of the collision operator can be calculated for any stepwise potentials /5, 6/. In the present paper, we perfect a numerical method of solving isotropic problems in the An analytic solution of the non-linear problem with arbitrary deviations a-representation. from equilibrium was obtained in /7/ for pseudo-Maxwellian molecules. Here this solution is Moreover, as also in used as a test problem for investigating the accuracy of the method. /3/, we make extensive use of the comparison with the exact moments of the distribution in We devote particular attention to the creation of an both linear and non-linear settings. internally-conservative numerical scheme, that is, schemes that do not violate the conservation laws without some kind of corrections to the collision integral or the distribution function.

a-representation of the Boltzmann equation. 1. (f(v,r,t)=f(lvi,t), the Boltzmann equation In the isotropic case space into a-space: dffdt=f (f, 1)-++=A with the use of Maxwell distributions f and

cp are connected

with arbitrary by the relation (see /l/j

is transformed

from

(cp, cp), temperatures

v-

(1.1) as basis

functions,

that is,

0

f(u,t)=SM(cr,v)cp(a,t)da. 0 Here

In (l.l), the M(a, u)-(a/n)Ytexp(-auZ) is the Maxwell distribution and a=m/(2kT). kernel of the operator A represents the mapping into a-space of the collision integral of two Maxwell distributions with parameters a1 and CL*,that is,

*Zh.vychisl.Mat.mat.Fiz.,28,6,901-916,1988

(1.2)

188

u), M(cQ, v) )-

J(M(al,

s 0

m

Mb, Wb,

al,

aAda,

(1.3)

~(cp,~)-SS~(a,a,,a,)cp(a,,t)cp(ar,t)dalda*.

(1.4)

0 Irrespective of the form of the section of the collision integral, the kernel has the following four properties /l/: 1)

A(a,

al, al)+O,

af+a2

(since J(M(al,v),Jf(a,,

A (a, al, at)

v))=O);

OD 2)

s

of the law of particle-number conservation);

(a consekence

0

3)

OD A (a7aCL1,01a) + A (a,a2,a ) 1 a s 0

dcz=O

(a consequence

conservation of energy);

(in some cases, for instance, for 4) A(a,arra2)=0 if acmin (a,,al) or aBalSa Maxwellian molecules, the region in which the kernel is specified is narrowed to the section [a,, 4 ). here, the functions 'In future,itwillbe convenient to transfer from a to T'-l/a; q'(T')=a*cp(a) and A'(T,',T,',T,')=a*A(a, a,, ax). Then, omitting the (0 and A are replaced by bars, from (1.1) and (1.4) we obtain d~(T,t)ldt=jSA(T,T,,T,)cp(T,,t)cp(T,,t)dT,dT,.

(1.5)

0 For isotropic distribution functions, the moments of the function f(u) are connected with the moments of the function q(T) as follows (in the general case, the connection between

the distribution function moments

f(v) and

m(a,a)

is shown in /5/j:

cc

v"(4nu'f(u))du,

M&-j

(1.6)

0 that is, to calculate any moment ofdegreekof f(v)wedonotrequireareturntothe v-space and it, in its turn, isthecorresponding moment of degree k of the function q(T). When developing numerical methods of solving Eq.tl.5) it is necessary to take account of the fact that q(T) can be a generalized function (so, in equilibrium, cp(T)=6(T-T,)). In this connection, it is best to change from the function 'p itself to some integral of it. We can, for instance, construct a numerical scheme by expanding the distribution function cp(T, t) in orthonormal functions oftherectangular momentum /8/:

W-T,-,) W’rT) (TrT,-,I”

a(Tl_

e(z)={;

that is, we seek

7

;=;

7

’

k=1,2,...,N,

cp in the form

cp(T,t)=

g,(t)%(T).

Substitutingthis representation of m(T, t) into the kinetic Eq.Cl.5) and integrating it with for the q,(T), we can obtain simple equations for the expansion coefficients g,(t) or populations of the temperature cells:

apThese equations take the form

I cpV)dT - G”,-T,-,)%n.

189

This scheme has, however, a drawback: in conserving the number of particles in the system, (~a,-l), it violates the law of conservation of energy. The scheme proposed in /l/, where the temperature relaxation problem is solved for the hard-sphere model, turns out to be simpler and more accurate. In /l/, at each instant of time the distribution function cp(T,t) is approximated by a collection of a-functions:

=v cp(T+t)=-

a,W6V-T,‘(~)).

a, =

F--f

5cp(WT.

(1.7)

Tp-’

T,’ is some temperature in the range Here (T,-,,T,). is easily made. For this purpose, it is sufficient which gives s

In fact, the transition to the u-space to integrate (1.7) with the Maxwellian,

that is, the distribution function in u-space is given by a finite collection Analogously, for moments Al,, using (1.71, we have distributions.

M. = x

T,‘“a,. P

of Maxwellian

(1.9)

in the form of (1.7) into (1.5) and integrating both sides of the Substituting q(T) equation in the interval (T,-,,T,), we obtain an expression for the time derivative of the population of the k-th cell at a time t:

N (1.10) where

A,,,= From

(1.11) and the second property

For sufficiently

small

At,

j A(T,T;,T,‘)dT. G.,

of the kernel,

using Euler's method,

(IAl)

it is clear that

we have

a,,(t+At)=a,(t)+(da,ldl)I,At.

(1.12)

T,‘(t+At). The presence of an In order to transfer to the new time step, we have to define additional free parameter enabled us /l/ to satisfy the law of the conservation of energy In /l/, we sought the temperature that is, to realise a fully-conservative numerical method. in the form T,‘(t+At) T,‘(t+At)=T,‘(t)+AT, The quantity AT is found from the conwhere AT was taken to be the same for all cells. dition that the energy of the system does not change. It should be noted that an analogous external correction when solving more complicated In /lo/ a special correcting polynomial kinetic problems was used widely in /9/ and /lO/. was introduced for correcting the distribution function, and the parameters of this polynomial The greater effectiveness of were chosen to be such that the conservation laws were obeyed. the conservative scheme was demonstrated. The question of the choice of correction method in problemswhere energy is not conserved The linear relaxation of a small addition on a background gas or the calculation remains open. In of the distribution functions of the components in gas mixtures can serve as examples. connection with these issues, we have to work out schemes that possess internal conservativeness.

Internally-conservative numerical method. A numerical method of integrating the Boltzmann equation is called internallyconservative

2.

if it guarantees an exact calculation of the rate of increase of the number of particles, of (We should the momentum and of the energy in each elementary cell at each instant of time. (a--u)-space /ll/, an internally-conservative remark that in solving the quasishock problem in (rr-T)method was used that guaranteed the satisfaction of the conservation laws in each cell). In such a scheme, when we are solving problems concerning the isotropic relaxation Of gas mixtures, the error in calculating the increase in energy can only be connected with a

190

discrete time step, and the satisfaction of the conservation laws will be guaranteed in isotropic relaxation in a one-component gas. scheme is constructed for problems of isoIn this section, an internally-conservative a-representation. tropic relaxation in the The method is verified in the problem of linear temperature relaxation - that is, the relaxation of a small admixture on an equilibrium backThis problem is a convenient test for trying-out the accuracy of the numerical ground gas. scheme, since an analytic solution has been obtained for it in the case of pseudo-Maxwellian molecules /7/. To guarantee internal conservativeness, we propose to introduce "energetic" cell populations Tr

b,=

Multiplying

both sides of (1.5) by T

j

TT(T)dT=a,T,*.

and integrating

over the interval

(T,-,,Th) we obtain

(2.1) Here h

Bw== .f Using Euler's

scheme by analogy

with

TA (T, T;, T,‘) dT.

(2.2)

(1.12), we have

b,(t+At)=b,(t)+(db,/dt)

ItAt.

(2.3)

After we have determined the number of particles ah and the energy bk in each temperature cell at the time t+At, according to formulae (1.12) and (2.3), we can determine the new coordinates of the /g-function T,'(tSAt):

b,(t+At) T,'(t+At)=-. n,(t+At)

(2.4

a,(t+At) and T,'(t+At) that are constructed in such a manner are theinitial The quantities parameters for the following time step. From the definition of Ba,, in (2.21, and from the third property of the kernel, it is clear that

This equality leads to the automatic non-linear problem.

3.

satisfaction

of the law of conservation

of energy

in the

Linear temperature relaxation.

In the linear temperature relaxation problem, a small admixture of particles relaxes on an equilibrium gas background (F,,=~lf(T~,u), Tb is the temperature ofthebackground, the molecular masses of the admixture and the background are taken to be equal), and at the initial instant the distribution function of the admixture is chosen to be Maxwellian with a temCorrespondingly, perature T..

cp(T, 0)=6(T-T.).

(pb=g(T-G),

Eq.(1.5) is simplified for such a linear problem, tends to the constant background temperature Tb:

dTU-9t) -=

dt

For pseudo-Maxwellian in the form /5, 6/

molecules

since, on the right, the temperature

s

AtT,T,,T,)cp(T,,t)dT,.

(3.1)

0

(a(g,O)=C/(4ng)), the kernel

T1

A(T,

T,,

T,) can be written

(3.2)

191

The analytical

solution

of problem

(3.1),

(3.2)

in

a-space

G(T-T,)+t’A(T,T.,Tb) y=ln-

T.-Tb

takes the form /?/

~,(2(fY)‘“) (Uy)”

(3.3)

1’

t’=tlT,

T--Tb ’

where

Z,(z) is a first-order modified Bessel function. In /7/, with reference to the case of Maxwellian molecules, expressions are also introduced for the distribution function expansion coefficients in terms of Sonine polynomials (Sonine moments) at any instant of time:

C.-S(1- T/Td"cp(T,t)dt-((1_T,/T,)"exp(h,t), where h, are the eigenvalues of the linear collision integral. For pseudo-Maxwellianmolecules, For any interaction sections, the moments M, are connected with the Sonine A,==-n/(n+l). moments C, by the simple relationship

The calculation scheme is so simple that it was implemented on a DZ 28 computer using BASIC. We solved the linear problem with initial conditions T.=5,T,==-l. Mote that the solution using standard Sonine polynomials diverges for these parameters /7/. The calculation was carried out by Euler's method using formulae (1.121, (2.3) and (2.4). We call such a scheme Scheme 1. The interval (T., Tb) was split into N equal parts. Using formulae (1.7)(1.9), we calculated the distribution function in a- and v-spaces, and the approximate The dependence of the moments moments. M,(t)/M,(O), n=l ,...,5 on time is shown in Fig.la. For N=lO and a time step ht=O.l, the approximate moment (1.9) and exact moment (1.6) on the graph are practically indistinguishable. Fig.lb shows the relative error 6M,, in percentage terms. The specific characteristic of 6M, is the presence of a maximum with At this intermediate time, the values of the high moments are small, and so the 4GtG8. absolute error turns out to be insignificant. We can assume that the presence of a maximum in 6M, is, basically connected with the discreteness of the time step. In fact, if ht is reduced by a factor of five, the error is reduced by a factor of about five. In connection with this, we implemented a transition to an internally-conservative scheme with second-order accuracy in time.

a

Fig.1 The easiest is to take account, analytically speaking, Taylor time series in the case of linear relaxation:

A’

Here the kernel of the operator -

of the second derivative

in the

takes the form _ _

AIP’(T,T,,

T,)== j A(T,T’,T,)A(T’,T,,T,)dT’.

(3.4)

. Substituting the expression for A from for pseudo-Maxwellian molecules:

(3.2) into

(3.4) and integrating

over

T’

we

have,

192

a'"(T.T,,T.)-A(T,T,,T.)(In~-Z)+I(T-T,). b As a result,

scheme

1 is replaced

by Scheme

2:

P-l

ck=ak, b,,

k=l,

P-L

2,. ..,N,

T,(t+At)=b,(t+At)/u,(t+At).

Here

‘ Tk

s

A(‘) (T, T,, Tb) dT = Af;,

cw=

(3.5)

( Tk-l =k s =k-1

A”‘(T,

TA(‘) (T, T,, Tb) dT = B$j, ,

T,, T,)=A(T,

i=l,2,

Tp, Tb).

in a calculation using scheme 2 (h,=0.1,N=10), the error GM..fallsbyl-2orderscomparedwith a calculation using scheme 1. Table 1 shows the influence on the error 6Ms and (IM, of changing the temperature step (t-1, k-0.1; the upper figure refers to a schemelcalculation, and the lower one refers to a scheme 2 calculation). We can see that, starting from N=1020, the size of the error has already become insensitive to further subdivision in T. The inclusion of the second time derivative plays a significantly greater role. Table OM,. %

N

T -

-2.40

4

-0.550

6

-1.92

-0.089 -1.07

10

-0.643 16

-1.82

20

-1.80

40

-1.79

1

0.008 0.019

0.017

-6.25 -2.78 -3.01 -6.43 -3.57 -0.22 -3.36 -0.027 -3.31 0.027 -3.26 -0.007

The second Sonine moment C,-

j (I-T/T,,)‘rp(T)dT,

is an important characteristic of the process - in the case of a monotonic distribution tion fq(T)>O, it corresponds to the square of the distribution function dispersion in space or to the excess in v-space: Ex=C,(T)‘-D’~(T-(T>)~=M,M,-M,‘. Here, for any quantity

q(T),

WV))-

funca-

(3.6)

j~U’MW’.

Figure 2,a shows the function q(T) ,for t-5. The solid curve corresponds to the exact solution (3.3). On the graph, the exact solution can hardly be distinguished from scheme 1 calculations with N-10 (dashed curve). (The approximate distribution function was calculated according to the formula q(T,‘)=a,lhr, where hr is the length of the temperature cell). Figure 2,b shows the relative errors &p for calculations with h*=O.l. The solid curve corresponds to a scheme 1 calculation with N=10, the dashed curve corresponds to a scheme 2 calculation with N-10, and the dot-dash curve corresponds to a scheme 2 calculation with N-20. The transition to scheme 2 significantly reduces the error, even with N=lO. Analogous curves for the function are shown in Fig.3. The exact distribution function f(r, 5), r-v' in v-space is calculated using the analytic solution (3.3). A careful analysis showed that the error in distribution function calculations (scheme 2, N=20) does not exceed 0.5% even in the far tails, (UCO,...,lO) 'although the distribution function varies by 12 orders of magnitude in this velocity interval.

193

For calculating more complex processes, in particular non-linear relaxation, an analytic account of the higher time derivatives becomes too cumbersome. So it is better to take them into account numerically. In the linear problem, a calculation was carried out using an improved Euler method that -by % 12 -

IO -

B-

6-

Y-

Z-

Fig.3

Fig.2

is one of the second-order runge-Kutta methods. The method of calculation (Scheme 3) in this case reduces to two consecutive Euler schemes with computation at each stage of the new points of application T,':at the first stage

cIb=afh,

b ,&,

k=l,

T,, . =b, ,+/a$,,;

2,. . . , N,

at the second stage

(3.7) czh--aah, b Ikr

k-i,

2,. . . , N.

The complete scheme per step is c,(t+L)-_(c,,+c*,)/2, T,‘(t+At)-bJa,.

cl-ah,

br,

k--l,

2,. . . , N,

We can easily see that such a method is internally conservative at each stage. A comparison of the calculation according to this scheme with the calculation according to scheme 2 shows that the results are practically identical. 4.

Non-linear temperature relaxation. We now consider non-linear temperature relaxation, that is, the process of temperature equalization in a one-component isotropic gas. In computing the non-linear process, the br do not differ in principle from those for linear formulae for the increase in a,,and relaxation. But in the non-linear variety, another summation appears in (3.6) and (3.7), C&, which are identical in form with Cc,l,except and the matrices CL? are replaced by Since a non-linear process is, in essence, a superthat Tb in (3.5) is replaced by T,. position of linear processes, an identical partition of the temperature interval guarantees the same accuracy of calculation as in the linear case. The accuracy of the numerical calculation for the non-linear temperature relaxation process can be checked using the moments. As is well-known, for Maxwellian molecules there is an analytical solution of the system of moment equations. For pseudo-Maxwellianmolecules, the system of Sonine moment equations takes the form

This system is recursively solvable. For the first six SOnine moments: G--l,

c,=o,

C,(t)-G(O)e-‘I’,

C,(t) = [C,(O) +3C,*(O) ]e--S”S-3Cp’(0)e-*“S, c,(t)-[C,(0)+2C,(0)Cs(0)]e-*“‘-2C*(0)C~(0)e-J~”.

C,(t) -C, (0) ewr’*,

194 Subsequently,thenumberof exponents with different indices grows rapidly - for instance, for it reaches 51, and for C,, it is 280. Hence it is best to use a computer to obtain cu the high-order moments. The solution can be represented in the form @("+1)-l c.=r,

k-8,“,

n-i, 2,. . . .

A "ksxP[o",(t) I,

To find On* we choose the following pairwise sums, all of which differ from each other: Op,+On_-p,,, S(p)(i
(4.la) (4.lb)

p(T,0)=+0.25)+$6(T-2.5).

Table 2 c,, CT,0)-0.5 16(r- 0.2) + a(r- nl

T

,

M,

Mm

rb(T- 2.5)

p(T,O)=+b(T-0.25)+

1

‘Ml0

1.302~10' 3.179.1033.104.10' 5.cW.10-'5.ooo-10-l 1.018.i03 6.422 1.940.10-' 2.169.10-' 1.018~iO~ 2.0X-40- 1.93610- 6.420 1.627.10*i.Oi2.10‘ 2.700 3.871.iO-2 5.937.10-2 1.625.W 1.010.i0’ 2.697 3.853.1O-2 5.919.10-2 2.426102 9.034 3.437.10-31.273 1.215.10~* 2.415.10* 3.393.10-s1.269 I 8.969 1.205.10-2 Table 3 ,&I. O)-0.5 [M(0.2. 0)+ Y (i. Oil

N,

2 1; exact dist. fn : 17 scheme 3 dist. fn N,

f (0.0)

7.084.10-' 9.957.W' 1.091 1.094

2.i57.1O-3 1.592.10-s ;M;I:

f(O.3) 5.049.10-' 5.737.10-1 ;9W&-

f(2.3) 1.089.10-3 9.519.10-‘ 9.625.10-' 9.621.10-'

.

t

.

f toI 0) -+

f (0.0)

2 178 scheme 3 dist. fn.

M (0.25, 9) + f

f(5.0) 1.058.10-" 7.039.10-" 1.249.10-'* 1.247.10-'2 f(5.3) ,.~ 3.913.10~1' 1.498.10~" I

2.i86.10-'3 2.172.10-'3

M (2.5. 0)

FlI.3:“o-~

! 17 exact dist. fn.

f (2.0)

-1.692.10*

9.729.10-’

f(0.3) 3.i89.1O-L 5.800.10-' -1.77310' 4.362.10~'

f (2.0)

f (5.0)

2.827.1O-3 7.092.10-"' -i.i18.1O-2 5274.10-7 -1.079.10-' 1.306 3.057.10-3 6.875.10-' f(5.3) f (2.3) 3.119.10-3 ;%:625;-:" 1.065.10-2 1.4t8.1O-i -4:186.10-8 1.232.10-' 3.296.1O-9

The calculations were carried out according to scheme 3 with N-40, h,-0.1. In addition, the same problems were solved using the standard expansion in Sonine polynomials.

(4.2)

here c,(t) was calculated analytically, as described above. Table 2 shows the exact moments (the upper figure) and the approximate moments M. (the lower figure). In Table 3, we compare the Sonine polynomial expansion with .~a=-&7 and 18 with the numerical calculation according to scheme 3 (N-40, h,=O.i).As we might have expected, the first problem, in which Grad's criterion is satisfied, gives good convergence, while in the second problem series

195

TJTeq~,=2.5>2. Figure 4a shows a distribution (4.2) diverges_, since at the initial instant a-space with the initial conditions (4.lb). The relative distribution function function in u-space with the same initial condtions is shown in Fig.4b. in

II

a

0

Fig.4

5. Relaxation of non-monotonic distribution functions. In the temperature relaxation problems described above, the distribution function cp(T) turned out to be positive. It is obvious from (1.2) that monotonic corresponds to f(v) rp(T) by the positive q(T). Moreover, since the excess Ex (see (3.6)) is connected with relationship

Ex=

i (T-T.)‘q(T,t)dT, 0

corresponds to positive cp(T,t). For Maxwellian molecules Ex(t)-Ex(O)e+', and then Ex>O is determined by the sign of Ex(0). The sign of Explays so the sign of Ex(t) an important role in the behaviour of the distribution function in v-space. It was pointed out in /12/ that for sufficiently large t and finite v the relative distribution function F(z, t)=f(v. t)/j(tT, m), is

a quadratic

trinomial,

whose branches

s=mvz/(2kT~),

are directed

upwards

when

Ex>O

and downwards

when

Ex+I. If we do not confine ourselves to functions m(T) of one sign, then not all the a, from (1.7) are positive and, moreover, in the interval (T,-,,T,) two &functions can appear with If we replace such a pair with weights that are different in sign and large in magnitude. h-function and conserve the sum-density and the energy, the relaxation process is one artificially accelerated in the given cell, and the point of application of the total-function can move far outside the boundaries of the cell. So, in the numerical scheme, one should make b-functions of different signs in one provision for the possibility of the appearance of We will write the two b-functions in the form cell. a,,6(T-T,‘)

-a&

(T-T,“)

-a,b(T-T,‘)

6(T-T,‘+e)-6(T-T,‘)

-

(5.1)

I*

d4 e In the case when e is small, the expression in square where ar=alr-ark, eIT,“-T,‘
(1.7) to

m(T)= In such well as analogy T,). In

V(T)

of alternating

sign as follows:

i [a,6(T-T,‘(t))-d,b’(T-T,‘(t))]. *-1

we have to determine the dipole moments in the cells, 4(t), as a representation, a*(t) and the points of application the populations T,(t). To this end, and by P and integrate over the interval (T,_,, with (1.10) and (2.1) we multiply (1.5) by the general form, all these three equations can be written as follows:

(5.2) Here

196

are the momenta

in the k-th cell, (5.3a)

M::;,- A::'CT,',T,'),

(5.3b)

(5.3c)

A,,, = j T’A(T, T,, T,)dT. I.4 (5.3) and the first two properties interaction sections

of the kernel

From

i=O,1,2,

Md&.=O,

it follows

p=1,2 ,...,

that for

any

particle-

N,

(5.4a)

(5.4b)

Eq.(5.2) It follows from (5.4b) that the new scheme also possesses internal conservativeness. is solved in the same way as (l.lO), (2.1). To start with (as also occurredinthetemperature relaxation problems we considered above) increases of the moments in the cells occur, and then, assuming the %"'(t+At)to be known, from the solution of the system a~=rnkco’, we determine

the parameters

ak,& and

T;

d,‘=m;“=

a,=m:“, (The temperature cell is possible

a,T,‘+d,=m,“),

a,Tk’2+2d,T,‘=m

at the instant

-m;‘m~‘=-_E4,

L‘*’

(5.5)

t+At: Tk’=-.

(1) mk -d,

(5.6)

4

T, may prove to be outside the k-th cell - that is, motion Here to non-monotonic relaxation problems).

from cell to

Exh= 7 (T-T,)$(T)dT Q-1 is to of we

The condition Ex,CO has to be satisfied for system (5.51 the excess in the k-th cell. be solvable. But in just this case it is also meaningful to take account of the derivative the b-function. Ex20, it is natural to set d,=O, that is, For those cells for which can use (1.7) inthesecells. For pseudo-Maxwellian molecules T. A'" 1~ -

and system

s T'A(T,T,, TJdT-e(T,-T,)9(T,-T,_,)T,’ %_I

(5.2) takes the form

(5.7)

Here

hi,-T,‘-T,‘,&

is the Kronecker

delta and

M:;; =

f A (T, T,‘, T,‘) T’ dT. TL,

The last term in (5.7) takes account of the "selfaction" in the cell - that is, an effect which cannot occur given the presence of one 6-function in the cell. Scheme (5.7), (5.6) was tested in the solution of non-monotonic isotropic relaxation problems. A numerical analysis showed that such problems can be solved with the same effectgives good possibilities' of iveness as monotonic ones, and the scheme so-constructed analysing the distribution function tails where Sonine-polynomial expansions arenotsuitable. the We consider a relaxation process that is described by the WKB-solution /13/. In a-representation, such a solution has the form /7/ ~e(T,t)=GIT-T,ft)l--d(T)G'(T-Tpftf). In the numerical

calculation

for the initial distribution

functions

in the form

(5.81 vB(T,O),

197

the moment

equations

for each cell

(5.7) take the form

which, as can easily be proved, corresponds to motion of the mapping (5.8) with a possible transition from cell to cell. In other words, in this problem the error of the numerical calculation is connected only with the size of the time step. We now consider the relaxation process where the initial distribution function is a sum of WEB and Maxwellian distributions (5.9)

cP(T,O)-(cpn(T,O)+cpaa)/2, where

(p&i(T-0.6)-0.46'(2'-0.6), ~JM==-~(T--T~).

It is interesting to trace how the relaxation of the distribution function tails changes, depending on the sign of Ex and the magnitude of the equilibrium temperature l'n for a change in

t and sufficiently large v, the behaviour of the distribution function is determined by the sign of the excess. Thus, a transition of the relative distribution function through unity up to the achievement of equilibrium is possible in the relaxation process - that is, we can have non-monotonic relaxation of the tails. ES(Ol j f/f,

e

a

f/f_0

‘O

‘E 5

’

5

I

I

I

02”

a

2

Y

6

8

Ia

5

laYI! lF

IO

0

2

4

6

8

IO

Fig.5 Figure 5 shows the corresponding relaxation processes at t= 0, 2, 5, 10, next to the The points a and e correspond to Es>O, and the points b, c, characteristic points a-e. tail, which d correspond to ExCO. Non-monotonic relaxation with superheating of the , was first noted in /lx/. The other cases of noncorresponds to case a (TM=0.36, T~~0.68) (T,=l.l,Tnmonotonic relaxation (with preuminary cooling of the tail) are c (T,=Tn--1)d the tail and e (T,=1.66, Tn=1.33)the direction of 1.05). In case b (T,=0.56,T,=0.78) relaxation, estimated from the sign of the excess, agrees with the initial behaviour and nonmonotonic relaxation does not occur. We may remark that in all the five cases we have considered, Sonine polynomial series the of must converge (TdTd2). However, practically speaking, the convergence region abbreviated series may prove to be insufficient for analysing the behaviour of the tails. v==2,3 In the problem described here, the Sonine polynomial series (4.2) converged up to No=7, and up to v-4,5 for No-20. with a number of terms REFERENCES 1. ENDER I.A. and ENDER A.YA., On a representation Nauk SSSR, 193, 1, 61-64, 1970.

of the Boltzmann

equation.

Dokl. Akad.

198

2. ENDER A.YA. and ENDER I.A., The Boltsmann equation in the e-u-representation. Izv. Akad. Nauk SSSR, Mekhan. Zhidkosti i Gasa, 4, 117-123, 1972. 3. NANBY K., Direct simulation scheme from the Boltsmann equation. Rept. Inst. High Speed Mech. 45, 348, 19-41, 1982. 4. GRIGOR'YEV YU.N. and MIKHALITSYN A.N., Numerical investigation of isotropic relaxation in a gas with Maxwellian interactions. Zh. vychisl. Mat. mat. Fiz., 25, 5, 742-756, 1985. 5. ENDER A.YA. and ENDER I.A., Integral transformation of the Boltzmann equation for different particle interaction laws. Preprint 605, FTI Akad. Nauk SSSR, 1979. 6. ENDER A.YA. and ENDER I.A., Boltzmann equation integral transform method. Analytical researches, In: Molecular Gas Dynamics, Nauka, Moscow, 1982. 7. ENDER A.YA. and ENDER I.A., Formation of the distribution function in temperature relaxation problems. Zh. Tekh. Fis. 54, 9, 1671-1680, 1984. 8. HARMUT KH., Theory of Sequential Analysis, Mir, Moscow, 1980. 9. YEN S., HICKS B. and AUSTIN R., Further development of the Monte-Carlo method for calculating Boltsmann collision integrals, in Mechanics and Dynamics and Rarified Gases, Mir, Moscow, 1976. 10. ARISTOV V.V. and CHEREMISIN F.G., A conservative disintegration method for solving the Boltsmann equation. Zh. vychisl. Mat. mat. Fiz., 20, 1, 191-207, 1980. 11. KOLYSHKIN I.N., ENDER A.YA. and ENDER I.A., Calculation of the Boltsmann equation collision integral and exact solution of the relaxation problem. Izv. Akad. Nauk SSSR, Mekhan. Zhidkosti i Gasa, 5, 132-141, 1977. 12. ALEKANIAN M., Integral representation for non-Maxwellian distribution and the approach to equilibrium. Phys. Lett. 4A, 1, 2, 1-5, 1979. 13. BOBYLEV A.B., The Fourier transform method in the theory of the Boltrmann equation for Maxwellian molecules. Dokl. Akad. Nauk SSSR, 225, 5, 1041-1044, 1975. 14. TJON J.A., On the approach to Maxwellian distribution. Phys. Lett. 7OA, 5, 6, 369-371, 1979.

Translated

U.S.S.R. Comput.Maths.Math.Phys., vo1.28,No.3,pp.198-205,l988 Printed

by H.T.

0041-5553/88 $lO.OO+O.CC 0'1989 Pergamon Press plc

in Great Britain

SEIDEL'S METHOD FOR THE TRANSPORT EQUATION WITH ANISOTROPIC SCATTERING IN R=* YU.V. KNYAZIKHIN

Seidel's method for solvinq the transport equation is considered. This method belongs to the class of iterative processes which is completely described Banach's theorem on the solution of operator equations of the second kind: the convergence and rateofconvergence are determined by the magnitude of the spectral radius of a certain operator. Thedependences of the magnitude of the spectral radius of this operator on the scattering indicatrix are studied.

Introduction. Consider

the stationary

single velocity

transport

equation

Here,

the unknown function cp(s,p) is the density of the particles which are travelling in the direction S-(st.s2,s1), ~s~‘--si’+s~*-i-ss’-=~ from the given point P==(s,,sa, z,);P is the surface of a unit sphere in R' with its centre at the origin of coordinates (we shall subsequently interpret the unit vector as a point on the unit sphere Q), the function which a(P), characterizes the absorption of the medium, is measurable, positive and bounded almost everywhere in the domain G where particle transport occurs: O-k(P)< dP, I4 is the scattering indicatrix on GX(-1,1) and the integral

sup vraia(P)==c&
which is a measurable,

lZh.vychisl.Mat.mat.Fiz.,28,6,917-925,1988

non-negative

function,

defined