16.
VORONOV, N.A. and Reptint
KONYUKHOVA,
107
of the energy eqtdon
Linearization
N. B. Stability of particle-like solutions of some wave equations.
ITEF, No. 26,1978.
17. VORONOV, N. A. Stability of time-dependent ITEP-127, Moscow, 1978.
particle-like solutions of some wave equations. Preprint
1s. ABRAMOV, A.A.et al, Calculation of the eigenvatues and eigenfunctions of ordinary differenti with singularities. Zh. vj%hisi,Mut. mat. Fiz., 24 S,1155-1173.1980.
equations
19‘ VAKHITOV, N. G. and KOLOKOLOV, A. A. Stationary solutions of the wave equation in a medium w&h saturation of a non-linearity. ID. vuzov. Rndiofiz., 16,7,1020-1028,1973. 20. FLUGGE, S. PTactiml quantum me&nics (Zadachi po kvantovoi mekhanike), Vol. 1, “Mit”, Moscow, 1974. 21. GLASKO, V. 8. era/. investigation of the particle-like solutions of the non-linear equation of a scalar field. Zh. eks. rewsfiz., 35,2 (8),452-457,1958. 22. KONYUKHOVA, N. B. On numerical isolation of the solutions tending to 2~0 at infinity of certain twodimensional non-linear sets of ordinary differential equations, Zh. @c&l. Mat. mot. Fiz., 10, 1, W-07,191o. 23. BAKHVALOV, N. S. Numeric01 methods (Chislennye metody), “Nauka”, Moscow, 1973.
U.S.S.R. Comput. Maths. Muth. Phys. Vol. 21, No. 1, pp%iOf-112,198l. Printed in Great Britain
0041-5553/81/0101.07-06507.5010 01982. Pergamon Press Ltd.
LINEARIZATION OF THE ENERGY EQUATION WHEN IT IS SOLVED JOINTLY WITH THE RADIATIVE TRANSFER EQUATION*
A.V.ADEEV and A.S.MEL'NICHENKO Chelyabinsk (Received 20 February 1979; revised 14 JamaT
1980)
NEWTON’S method ifused for the joint solution of the unsteady spectral equation of radiative transfer and the energy equation. The energy equation is linearized over the temperatures at two adjacent points. In the solution of the transfer equation the derivatives are estimated by the Monte Carlo method. In the joint solution of the unsteady radiative transfer equation and the energy equation the question arises of the convergence of the intermediate iteratians between these equations and the rate of convergence for a given time-step. In this paper Newton’s method is used to accelerate the convergence of the iterative process. The integrals occurring in the energy equation are Iinearized over the temperatures at two adjacent points. It is known that the principal difficulties arising in the joint solution of the transfer equation and the energy equation are caused by the heating up of optically dense matter. In this case the approximation of the temperature - dependence of the absorption energy of adjacent points only is justtied; at the same time calculations show that this approximation does not worsen the convergence of the iterations at points with high temperatures. *zJr. v_%W. x?r. mu?. Fir.., 21,1,107-512,
19$1.
A. V. Adeev and A. S. Mel ‘nichenko
108
In the solution of the transfer equation the derivatives are estimated by the Monte Carlo method. The efficiencies of this scheme, of the iterationless scheme described in [ I] and of a scheme with simple iterations are compared by solving model problems. 1. The transfer of radiative energy in an immobile plane medium and when there is no scattering is described by the system of integrodifferential equations 1 az, --+p
c
at
$
=xu(z,~-z”) 2hva
I,, = 0
(1*I>
(
-1
c2
1 exp (hv/kT) - 1 ’
(1.a
Here Z, = 1, (x, t, cc)is the spectral intensity of radiation of frequency v, x,=xv (T) is the absorption coefficient allowing for stimulated emission, T is the temperature, ~1is the cosine of the angle between the x-axis and the photon flight direction, E is the internal energy of the substance, c is the velocity of light in a vacuum, h is Planck’s constant, k is Boltzmann’s constant, and Ivp is Planck’s function, The system (1 .l), (1.2) is closed by a system of boundary and initial conditions, whose specific form is unimportant for us. Various methods for the joint solution of Eqs. (1 .I), (1.2) are considered in [2-4]. Thus, in [2] the convergence of an iterative process for isolating the radiative heat-conduction term is investigated. In [3] Newton’s method is used to linearize the right side of the steady transfer equation. Here the derivatives are estimated by ordinary difference ratios. Some principles of the construction of implicit schemes, where Eq. (1.1) is considered in the diffusion approximation, are explained in [4]. In this paper the linearization of the integrals occurring in the energy equation is proposed. The temperature values at the preceding step are taken as the initial approximations for solving the unsteady transfer equation. 2. We consider Eq. (I .2) in the discrete approximation:
(2.1)
where
s is the index of the iteration, i is the spatial index, n is the time index, r is the time step, Ti are the temperatures of the spatial points, and m is the number of spatial points. We linearize the functionsJi, Bi, Ei over the temperatures:
(23
Linean~ation of the energy equation
E:+i--Erg + !$T,'+'-T,'). *
109
(2.2 Cont’d.)
Solving Eq. (1 .I) by some method for futed Tf, we cm compute.lis, 8 Jis/a 5 and from system (2.1) determine Tp+l with the aid of (2.2). System (2.1) in the approximations of (2.2) is a system of linear equations of order m in the Tis+l. The additional computation of au the derivatives 8 Jplaq, i, j = 1,2, . . . , m, can considerably reduce the efficiency of the scheme considered. The idea of the method consists of simplifying the approach. We will omit the relation Jfs+l = Ji (TlS+l , . . . , Tms+l) from all the cells except for i itself and the adjacent i- 1, i + 1. Then I;+’ -Ii’
+
g(T:,:‘-T:,,), j=-r,O,l I+1 c
(2.3)
Here we note that Eq. (1 .l) is solved in the complete formulation, where Jfl=J(( T,‘, . . . , T,‘) . We will use the approximation (2.3) only to solve Eq. (2.1) to solve the convergence of the iterative process. In particular, this approximation will be sufficient for small photon flights in the substance being heated. With condition (2.3) system (2.1) is solved by the pivotal condensation method. 3. We will consider the application of the iterative scheme of section 2 to solve Eq. (1 .I) by the Monte Carlo method. The estimation of the quantity 8Jp/i3 T when compared with the “standard” estimate OfJiS is changed as follows.
Let j = 1, if a photon enters the cell i from the left;j = 1 otherwise; j = 0, if a photon is radiated from the i-th cell; xl=xV (Ti+j) ; rj = ri+j is the distance between the boundaries of the (itj)-th cell in the direction of flight of the photon; I is the number of the source, L is the number of the trajectory; IVLris the weight of the I-th photon radiated from the source with number I; pL[ is the probability that the L-th photon from the source with the number I is not absorbed along the path to the cell i-l (or it 1);Nl is the number of photons radiated from the source with the number 1. HereI,i=1,2,.
. , ,m. Then
l-1
L-1
where ELI is the contribution of the L-th history of the I-th source to the quantity Ji.
A. V. Adeev and A. S. Mel ‘nichenko
110 IflZi-l,i,
it 1, then
ch=W~~p~
I,
[I-exp(-w-0)
eXp(--xjrj)
Similarly,
If I = i- 1, i + 1, we assume that the spectral cross section of absorption allowing for the stimulated emission %=xV (2’) is representable as
x,(T) = q[
I-exp(
-f)]
,
where f(r) is some function differentiable with respect to T. Then apart from a numerical coefficient xy (2’) I,, (T) -f (I’) exp (-V/T). From this for I = i-l, i + 1 &J=f(Ti+j)exp(-Y/Ti+j)exp(-Xjrj)
[1-exp(-Xoro)],
8x*,
+v-r,l-
. aTi,,
T2i*j
>
Similarly, for I = i au=f(Ti)exp(--y/T,)
+ro-r
8x0
exp
[I-exp(-~oro)
I,
(-xoro)
dTi’l-exp(-Xoro)
Therefore, when estimating the derivatives the results recorded must be taken with the weights given in curly brackets. The greatest efficiency of the method is attained in the simultaneous estimation of Ji and aJi/aT on the same trajectories, since these quantities are integrals of high multiplicity and their estimation by a Monte Carlo method requires considerable computational time. We note that it is possible to consider twodimensional problems, retaining in principle the given scheme for computing the derivatives and an approximation of the type (2.3). In this case the derivatives are estimated from all the adjacent cells surrounding the given i-th cell and having a common boundary with the i th cell. However in this system (2.1) will not degenerate into a tridiagonal system, therefore at each iteration instead of a pivotal condensation it is necessary to solve a system of linear equations of the m-th order in TiS+' .
Linearization of fhc energy equation
111
4. Example. Consider a plate of thickness 4 cm, heated to T = 1 keV by an absolutely black source at the point x = 0, cy = 0.81 (see [ 11). Equation (1) was solved by the Monte Carlo method by a scheme with “anticipation” [S] using the technique of dependent trials [6]. The spread of spectral radiation in the plate was calculated by the following iterative schemes for solving system(l.l),(1.2). 1. The iterative scheme of [5] with (Y= 0.5: a) we solve Eq. (1 .l) for TiS, b) for the replacement of the index s + 1 by s on the right side of Eq. (2 .l) we find Tic”+1)~~ c) T:+’ -0.5 (Tt(“+i)’+Tc”) , 2. The iterative scheme described in section 2. The iterations are regarded as terminated when IT,‘+’ -Tc’lc~T:+~, where c is the specified accuracy of convergence of the iterations It was assumed that in all the calculations AY= 0.4 cm, E = 0.01 ? the number of photons emitted in the whole system was N = 5000. The number of photons was distributed over the sources in proportion to the energy of radiation. On reducing the statistics by a factor of two the results changed negligibly. For e = 0.001 the number of iterations remained approximately the same as for E = 0.01. Roblem I
x.,T)=$[
l-exp(
-f)]
,
where P, Tare measured in keV. The calculation was carried out with r = 2.10-S psec. The results obtained by both schemes agree well with each other and with the results in [ 11. For the calculation by scheme 1 the number of iterations per step amounted to 6-7 for heating the cold matter, at high temperatures the iterations disappeared. In the calculation by scheme 2 the number of iterations reached 34 for heating the cold matter, and for ct > 15 cm the process converged after 2 iterations. The relative temperature increment per step for ct < 60 cm was greater than e. At high-temperature points the iterations converged practically immediately, and the entire iterative process was continued on the cold matter. For 7 = 10-4~~~ the number of iterations at the first step reached 20, and for 7 = 5.10-d bsec it reached 40. Subsequently the number of iterations fell sharply to 2. As for r = 2.10-S psec, the process converges slowly only at low-temperature points.
A. V. Adeev and A. S. Mel’nichenko
112
The iterations carried out by scheme 1 diverged for 7 > 10-s psec. Problem II. The effect of an optically thick layer on the distribution of spectral radiation was considered. The conditions are the same as in Problem I, except that for 2 cm
X”(T) = qq
i-exp
(-4)1.
Since the domain 2 cm Q x Q 2.4 cm is practically opaque in the basic spectrum of the radiation, the photons in this domain are not selected uniformly over its volume, but only close to its edges, where because of the small distances the matter becomes transparent. The variation of the sampling was taken into account by the statistical weight of a photon. The problem was computed by schemes 1 and 2,~ = 2.10-s psec,N = 2500, E = 0.01. Up to ct = 15 cm the calculations agree well with each other and with [ 11,the number of iterations then reached 8. For ct > 15 cm the iterations by scheme 1 disappeared, while the number of iterations by scheme 2 fell to 3-4. The computation of problems I and II showed that the approximation (2.5) does not worsen the convergence of the iterative process at high-temperature points. 5. The computations of problems I and II showed that scheme i converges slowly at high temperatures and also diverges for large values of 7. In the iterationless scheme of [I] it was shown that the number of fictitiously scattered photons depends linearly on the mean Planck cross section and on the step 7. Therefore, in this scheme an increase in 7 or of the optical density of the matter entails an increase in the number of scatterings, and consequently, of the computing time. a considerable extent scheme 2 avoids all the drawbacks of schemes 1 and [l ] mentioned: a) the iterations mainly occur at the temperature front, b) a calculation with large 7 is possible by increasing the number of iterations at only the fust steps. To
Translated by J.
Berry.
REFERENCES 1. FLECK, J. A. Jr. and CUMMINGS, J. D. An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport. J. Conrpur. P&s., 8,313-342, 1971. 2. CHETVERYSHKIN, B. N. A method for simultaneously solving the energy equation and the radiation transfer equation. Zh. v~chisl.Mat. mat. Fir., 10,5,1290-1292, 1970. 3. ZUEV, A. I. Application of Newton’s method for solving the problem of the propagation of non-equilibrium radiation. Zh. vychisl. Mat. mat. Fiz., 13,3,792-798, 1973. 4. CHARAKHCH’YAN, A. A. Calculation of the unsteady spherically symmetric flows of a radiating grey gas. In: Dy~mics of mdiotinggas (Dinamika izluchayushchego gaza). No. 1,54-74, VTs Akad. Nauk SSSR, Moscow, 1974. 5. MELWCHENKO, A. S. and OGIBXN, V. N. Application of the Monte Carlo method to the solution of spectral problems of radiative heatexchange. Zh. vj%hisLMr. mat. Fiz., 17,4,1068-1074,1977. 6. FROLOV, A. S. and CHENTSOV, N. N. The use of dependent trials in the Monte Carlo method to obtain smooth curves. In: Roceedings of zhe Sixth All-Union conference on the theory of probability and mathematicalstatistics(Tr VI Vses. soveshchaniya po teorii veroyatnostei i matem. statistike), 425-427, Gos. izd-vo LitSSR, Vil’nyus, 1962.