Numerical solution of an inverse problem of the dynamics of a radiating gas with axial symmetry

U.S.S.R. Comput.Maths.Math.Phys.,Vol.27,No.4,pp.78-82,1987

OO41-5553/87 $10.OO+O.OO 0 1 9 8 8 Pergamon Press plc

Printed in Great Britain

NUMERICALSOLUTION OF AN INVERSE PROBLEMOF THE DYNAMICSOF A RADIATING GAS WITH AXIAL SYMMETRY" V.I. GRYN,

V.I. ZUBOV and V.M. KRIVTSOV

A numerical method is proposed for solving an inverse, non-stationary, axially symmetric problem of the dynamics of a radiating gas which is associated with the direct problem of the action of laser radiation on an aluminium vessel and its vapour. Results of numerical calculations are presented. A formulation and a numerical method of solving an inverse, diagnostic problem in the dynamics of a radiating gas with spherical symmetry were put forward in /i/. In the same paper it was also suggested that it may be possible to solve analogous problems when there is axial symmetry. The algorithm described here also enables one to demonstrate the numerical solution of an inverse problem of the dynamics of a radiating gas when there is axial symmetry, and is a synthesis of an algorithm for solving the direct axially symmetric problem of the dynamics of a radiating gas /2-5/ and an algorithm for solving the inverse problem for the characteristic radiation transport equation. The algorithm for solving the inverse problem for the characteristic radiation transport equation and, also, the construction of the synthesis of the algorithms were proposed by the first of the authors. The inverse problem of the dynamics of a radiating gas which is investigated here is associated with the following direct problem. At t=0, there is solid aluminium in the halfspace x < 0 and a vacuum in the half-space x > 0 (t is the time, and x is the axial coordinate). When t>0, laser radiation acts on the aluminium and its vapour. The direction of the radiation is opposite to that of a unit vector ix directed along the x-axis. The frequency VL of the laser radiation is such that hvL~1.16 eV (where h isPlanck's constant). This corresponds to a neodymium laser. The intensity ]L of the laser radiation is specified by the formula

f

3J o (r) tt-L

0 -~
JL (t /3, r),

t./3.<~ t ~ t /2,

SL(t,r)~--- ~ .TL(t --t,r), t /2~<t.5t, r/>O, J., 0 ~ r ~ r J2, ]o(r)= 2 ( l - - r / r . ) J . , r./2r.; where

t.=10 ns, r.=5 mm, J.=72 (7n)-tx 109 J/cm2.s and r is the radial coordinate. The axially symmetric motion of the equilibrium vapour is described by a complete system of equations for the single temperature dynamics of a radiating gas without any allowance for viscosity and thermal conduction. Transport of the characteristic radiation is allowed for under the assumption that there is local thermodynamic equilibrium and no scattering. Vaporization is modelled by an isothermal discontinuity with the condition that the vapours are saturated. The mathematical description of the phenomenon corresponds to that given in /2-5/. The system of equations for the direct problem has the form (see /2/)

0(r I 0(ruI 0(rw Ot u c)x w ~ r + vo~ = O, Ow Ow Ow Op --+u--+w-g~/+ =% c)t Ox r V~r 7/'rt--v [rw(uZ+w~+E+pv)]2

(2) (3)

,+~-x,--v0 [ r u ( -a2+w~ ~ - - - + E+pv)]+

~[ ~( u2+w: + E)l=r[ k = J - d i v 2 OJ±/c3x=-4-k~]~,

J=J++]-,

*Zh.vychisl.Mat.mat. Fiz.,27,7,1078-1084,1987

78

(4)

I~d~dv)] o4n

O
(5)

79

Yo

Ot

Ox K--~-x i ,

x<0,

p=p(E,v),

(QV)I+× (v, T, v)l(t, x, r, ~, v)=×B(v, T),

T=T(E,v),

(6)

0
(7)

where ~ is a unit vector of solid angle, I is the intensity of the characteristic radiation, × is the absorption coefficient of the characteristic radiation, v is the frequency, T(t,x,r) is the temperature of the vapour, u(t,x,r) is the specific volume of the vapour, B is Planck's function, C is the specific heat capacity, u0 is the specific volume of the condensate, T.(t, x, r) is the temperature of the condensate, K is the thermal conductivity, E(t,x,r) is the internal energy of the vapour, p is the pressure of the vapour, J+ is the intensity of the incident laser radiation, ] - i s the intensity of the laser radiation which is reflected as by a mirror, J(t, x, r) is the intensity of the laser radiation, ks(u , T) is the absorption coefficient of the laser radiation, and u(t, x, r) and w(t, x, r) are the components of the velocity vector of the vapour along the x and r-axes respectively. Here and subsequently (where it does not give rise to any misunderstanding) the symbols g, g(x). g(z),g(x,z) etc. are employed instead of the symbol g(x, y, z). The boundary and initial conditions and the conditions of matching with (1)-(7) have been described in detail in /2/. The method for solving the direct problem and its implementation on a digital computer is the same as that presented in /2-5/ while the conservative scheme from /6/ was employed to calculate the laser radiation. The inverse problem differs from the direct problem in that, in the direct problem, u, w, v, p, T, E, To, ]±. I are the required functions while, in the inverse problem, the function × and the required functions of the direct problem are to be found. Moreover, the boundary conditions for the inverse problem contain one additional condition

I(t, x, r÷(t, x), Q, v)=l+(t, x, Q, v),

(Qix)=O,

(Qir)~O,

(8)

and, when {r>r+(t,x),x>O,t>O}, there is a vacuum. Here l+ is the intensity of the emitted characteristic radiation and ir is the unit vector along the r-axis.

Remark. The boundary conditions on (7) in the direct problem are as follows. When {x=X0(D, rO}U{r=R,OO}, zero incident radiation from without is specified and, when {x=O.rO, r>O, t>O}\{OO}. It is assumed that instead of the point function I+ the functions 1+6, 6 are known such that

II+(t, x, ~, v ) - l ÷ q t , x, ~, v)l~6(t, x, Q, v)l÷ a,

(~&)=O,

(ei~)~O.

The functions 6, I+~ can be found from experimental measurements /7/. In this paper the functions l+~, 6 are determined from the solution of the direct problem (as in /i/). In /i/ and /8/ a formulation of the inverse problem, corresponding to (7) and (8), for the characteristic radiation transport equation (for the problem of determining x and I from specified u, T, B, ]+ and the intensity 1_ of the characteristic radiation which is incident from without) is presented. Problems of the existence, uniqueness and stability of its solution are also considered in the same papers. Let us now describe the proposed numerical method for solving the inverse problem in the dynamics of a radiating gas. At each instant of time t=%j~l ({~,]=:0, I,...; % = 0 } is a time mesh), the absorption coefficient x(u, T, v) is determined from the solution of the inverse problem for the characteristic radiation transport (when doing this, the functions u and T in (7) are taken from the preceding time layer t=ij-~). Eq. (7) is then solved for I with the calculated function x(u, T, v) using a certain modification of the method from /9/ and /2/. The functions u, W , u , p , T , E , To, J ± at the instant of time t=Tj are then calculated on the basis of the method from /2, 4t 6/. The above-mentioned modification of the method from /8/ and /2/ involves the following. The matrix (g(×))(v,T) and the vector (~(~))(v,T) for the momentary treatment of Eq. (7) with respect to ~ which are defined in /2/, were calculated in /2/ prior to solving Eqs. (i)-(7). In solving the inverse problem, corresponding to (1)-(8), the absorption coefficient for the characteristic radiation x~(u, T, ~), which is recovered using 6, I+8, depends on t since the set of points (u, T) at w h i c h x 6 is recovered depends on t. (When 6 > 0 , the quantity z6 also depends on t so that, generally speaking, ×6~×8(u, T, v; l)). In accordance with what has been said, the matrix M and the vector i are recalculated at the instant of time t=Ti using the algorithm from /9/ and /2/ on the set of points (y, T), on which the absorption coefficient ~(u, T, v) changes upon passing from the instant of time Tj-~ to the instant of time Tj even if it changes for just some subset of frequencies v (with a non-zero Lebesgue measure). The logarithm of the mass absorption coefficient ~=In(×u/×.u.) is sought in the form of a table

~(v,T, vh)=~p, m-t<~-where

P~(m,n,k),

In (v/v.)-l.54 1.15

< m,

m , n , k = l , 2. . . . . n-t<~--

In (T/T.) +0.45 0.t5

(9) < n,

x.=l cm -1, v.=1 cm3/g and T.=leV. We note that a table of the logarithm of the mass absorption coefficient is known for

80

aluminium in the form of (9) and this was used in /2-5/ and is employed in the present paper to solve direct problems. A comparison is made between the table of ~ values recovered from the solution of the inverse problem corresponding to (1)-(8) and the known table of values. We will now describe the algorithm for calculating the elements ~,~ of the table of values. The solution of the inverse problem for the characteristic radiation transport equation depends on t,x and ~' are parameters (see /i/). For each (~,~a), ],k=l, 2,..., an approximate characteristic radiation absorption coefficient x~(r; ~,x~,va) is determined from the solution of the given inverse problem by a regularization method /i/ at certain x=xt(r¢)> 0, /=I, 2 .... , L(%, v~) . Using ~(r; ~, X. ~), v(r; %_~, x3, T(r; %_t, x~) for certain integers m =m(x~, xl, ~), n=n(T~, x,, ~ ) , ~" (the P-element of the ~-table) , is determined using the method described in /10/ and /ii/. The previously calculated element ~ of the ~ -table is modified on the basis of ~* using two techniques. The first technique corresponds to that used in /iO/ and /Ii/: where

~-,-~ (~vN.+[/) (N~+ t) -', ~.

N.

N~N~+I,

(10)

is the number of calculations of The second technique has the form

(lla)

c ~ [ (t-c¢/)(x.:+ (1-c~p) (~/)~1 (ct.+cz/-2a.ct/)-'. The quantity ap" is calculated on the basis of the function by the equalities a=exp[--56(1+4] ~--0.51~)],

a=0.1 exp[-106s(l+4[~-0.5[~)] t<'~<2, 6s=6, cz=10 -~, where

0<~8~
(llb)

~(r; t, x, %')

(t2a)

6s=5,

0 <'~<2,

which is defined

6s>6

or

-r°>2,

(12b)

(12c)

5(r; t, x, v)=6(t, x, ~, v), (Q&)=0, (Qi~)=[l--(r/r+)~] '~',

• ~(r;t,x,v)= i ×~(~;t,x,v)dg,

~=(r~+g~) '~', s=(r+~--r~)'~:;

o

and the quantity 6s, which is the overall error in the input data for the inverse problem in the case of the transport equation, is defined in /i/. The magntiude of ~p" is calculated from a comparison of the functions ~(r; t, ~ ~), v(r; t, x), T(r; t, z) using exactly the same technique as that employed to calculate ~p* using v, T, ~°(r; t, x, v), where ~6(r; t, x, v)=In (×sv/M.v.), M°=M6(r; t, Z, ~). When t=0, it is necessary to specify in some way the values of ~p, Np, ~p for any permissible P (see (IO) and (Ii)). It is assumed in the numerical calculations that, when t = 0 , the values are ~p----]0, Np~0, ~p--0. The magnitude ~p~--10 corresponds to x6=0. When ~p~--i0, the characteristic radiation thereby has no effect on the motion of the vapour. The quantity ~ corresponds to a lower estimate for min (×6, z)/max(~, ~) which holds in numerical calculations on analogous inverse problems of the dynamics of a radiating gas when there is spherical symmetry /I/. The case (12a) corresponds to an inverse problem of "good quality", (12b) to a problem of "satisfactory quality" and (12c) to a problem of "poor quality" In the numerical calculations presented below, 6--0.1. If ~p~exp(-1), ~p is thereby determined from the solution of inverse problems of good quality. We note that the quantity ~p, strictly speaking, does not characterize the accuracy of the solution x 6 of the inverse problem. The numbers m and n, at which ~n~ was calculated for precisely one k ~ i prior to the instant of time ti, i=I,...,7 are shown in Fig.l. The minimum of these numbers i occurs in the cell m - - 1 < ~ < m , n--1<7
81 to the function ~m, which is recovered by the technique (ii), (12) with an accuracy ~P~ exp (--I) up to the instant of time t2. The accuracy with which ~ is recovered using techniques (i0) and (ll), (12) is practically identical. At the same time the technique (ll), (12) indicates the points P at which the quantity ~r is determined mainly from a good-quality solution of inverse problems concerning the characteristic radiation transport equation, and consequently, with better possible accuracy. The set of points P, at which an element ~/ is calculated precisely once in the case of technique (iO), is practically identical with the analogous sets of points P in the case of the technique represented by (ii), (12) with ~p~10-' (for any t>0). Out of approximately 3000 elements of the ~ -table recovered using the technique (ii), (12), only 323 were of good quality (up~exp(--1)). ' I 0 1 2 3 2 q

I ii 2 2i2 5 7

~'

Z!! q 2 2 5 ! 1122256 !

q 1 I

[ ~

Z 21 q q

22 q 33225 5

~

II

~"~.

4'

"~ .......... .~" -

B

fO

0.3

~

Fig.1

I

08

}...

. L3

h~

Fig.4

"'%"

"

I

0

o

Z5

I

r} z

Fig.2 Fig.3 The results of the numerical calculations which have been presented ahow that, on the basis of the non-stationary, axially symmetric, inverse problems of the dynamics of a radiating gas of the type proposed in the present paper, it is possible to find the gas-dynamic functions in the case when the characteristic radiation has a substantial effect on the motion of the gas, the absorption coefficient of the characteristic radiation is unknown and the characteristic radiation emerging in planes which are orthogonal to the axis of symmetry is known with a certain error (not greater than 20%) from experimental measurements. Let ~ be a set of points (u, T, ~) such that variations in the absorption coefficient of the characteristic radiation on this set on substantially change the flow of gas. Moreover, let the absorption coefficient be known on a large part ~' of the set ~ (in the sense of the Lebesgue measure). Then, on the basis of the inverse problems being considered, it is possible to find the gas dynamic functions and to recover partially the absorption coefficient on /~\~' from the good quality of the inverse problems in the case of the characteristic radiation transport equation. At the same time, the error 8=Ix6--xl/min (x6, x) in establishing the absorption coefficient does not exceed 1 for the most part. such an accuracy in the specificaticn of the characteristic radiation absorption coefficient is acceptable in a number of problems on the dynamics of a radiating gas. On the basis of the inverse problems of the dynamics of a radiating gas being considered, it is also possible to compare fragments of the absorption coefficient table which has been recovered with tables calculated on the basis of quantum mechanical theories. Such a comparison may yield additional information on the limits of applicability of these quantum mechanical theories for determining the characteristic radiation absorption coefficient. The authors thank I.N. Naumova for help in debugging the program and Yu.D. Shmyglevskii for useful advice and interest. REFERENCES 1. GRYN V.I., On inverse diagnostic problems of the dynamics of a radiating gas, Zh. vychisl. Mat. mat. Fiz., 25, iO, 1506-1525, 1985.

USSR 27:4-F

82 2. ZUBOV V.I., KRIVTSOV V.M., NAUMOVA I.N. and SHMYGLEVSKII YU.D., Calculation of the motion of the vapour of a solid under the action of laser radiation, in: Dynamics of a Radiating Gas, Vychisl. Tsentr Akad. Nauk SSSR, Moscow, Issue 3, 76-104, 1980. 3. ZUBOV V.I., KRIVTSOV V.M., NAUMOVA I.N. and SHMYGLEVSKII YU.D., Calculation of the interaction of laser radiation with an aluminium vessel and its vapour, Zh. vychisl. Mat. mat. Fiz., 20, 6, 1513-1524, 1980; Correction, Ibid., 21, 2, 521, 1981. 4. ZUBOV V.I., KRIVTSOV V.M., NAUMOVA N.I. and SHMYGLEVSKII YU.D., Interaction of the radiation from a neodymium laser with an aluminium vessel and its vapour, in: Dynamics of a Radiating Gas (Dinamiki izluchayushchego gaza), Vychisl. Tsentr Akad. Nauk SSSR, Moscow, 56-73, 1981. 5. ZUBOV V.I., KRIVTSOV V.M., NAUMOVA N.I. and SHMYGLEVSKII YU.D., Calculation of the action of a laser on a planar obstruction and its vapour, Zh. vychisl. Mat. mat. Fiz., 23, 6, 1520-1522, 1983. 6. GRYN V.I., Schemes for calculating radiation transport, in: Dynamics of a Radiating Gas, Vychisl. Tsentr Akad. Nauk SSSR, Moscow, 17-55, 1981. 7. KEY M.H., LEWIS C.L.S. and LUNNEY J.G. et al., Time-resolved X-ray spectroscopy of laserproduced plasmas, Phys. Rev. Letts., 44, 25, 1669-1675, 1980. 8. GRYN V.I., On inverse problems in the theory of radiation transport in two-dimensional geometries, Zh. vychisl. Mat. mat. Fiz., 27, 3, 441-455,. 1987. 9. SHMYGLEVSKII YU.D., A version of the method of moments for calculating the transport of selective radiation, Zh. vychisl. Mat. mat. Fiz., 17, 3, 785-790, 1977. i0. GRYN V.I., On the regularization of the problem of determining the absorption coefficient from the radiation emitted from a moving gas, Chisl. Metody Mekhan. sploshnoi sredy, Nauka, Novosibirsk, 15, 3, 41-59, 1984. ii. GRYN V.I., Determination of the absorption coefficient from the radiation emitted from a moving gas, VNTITsentr, Moscow, GosFAP. POO5502 from 04.06.82.

Translated by E.L.S.

U.S.S.R. Comput.Maths.Math.Phys.,Vol.27,No.4,pp.82-89,1987 Printed in Great Britain

OO41-5553/87 $iO.OO+O.00 0 1 9 8 8 Pergamon Press plc

THE CHOICE OF THE NUMBER OF PARTICLES FOR THE "PARTICLES IN CELLS" METHOD*

L.B. ZHUK

The step operator error for a finite number of particles in a cell N is estimated for the "particles in cells" method in steady flows. Conclusions on the convergence of the method are drawn from this. The value of necessary to obtain a required accuracy in the solution is indicated. i. I n t r o d u c t i o n . The "particles in cells" method (p.c.m.) proposed in /1-5/ for calculating hydrodynamic problems, in addition to a fixed Euler cell grid, uses a Lagrangian particle grid whose mass is concentrated at points. The motion of a material is viewed as the displacement of mass points through the cell grid. This method was proposed for calculating large media deformations. It is usual to substantiate p.c.m. /2-5/ in terms of the splitting of the equations of hydrodynamics, written in a special form for continuous media, i.e. where the particle number 'N ~ . When N is sufficiently large, by virtue of a passage to the limit, the p.c.m. must be close in its properties of convergence and stability to an asymptotic difference scheme derived from the splitting. However, when N is finite, the "quantized" flow of material, being a source of additional error in the step operator compared with a continuous medium could strongly distort the convergence picture, compared with an asymptotic difference scheme. In /2/ it was shown that a small particle number (of the order of four or five) in a cell in the p.c.m, could lead to a significant non-physical fluctuations in calculated values of the flow field. These fluctuations are due to the motion of particles across the cell boundaries and are inevitable due to the quantization of material mass. In /2/ it was noted that fluctuations could entail a subsequent deviation from the expected value of the mean pressure, and to improve accuracy it was ~ro~osed to correct the pressure at certain time

*Zh.vychisl.Mat.mat.Fiz.,27,7,1085-1094,1987