CC BY
7Russian Journal of Nonlinear Dynamics, 2020, vol. 16, no. 2, pp. 325-341. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd200207
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 65Yxx, 65Zxx, 65Mxx, 78-04
Development of a Mathematical Model and the Numerical Solution Method in a Combined Impact Scheme for MIF Target
V. V. Kuzenov, S. V. Ryzhkov, A. V. Starostin
This work is devoted to the theoretical calculation of the processes of compression and energy release in the target by a combined action of a system of pulsed jets and intense laser radiation using a magnetic inertial plasma confinement method. A mathematical model, a numerical method, and a computational algorithm are developed to describe plasma-physical processes occurring in various types of high-temperature installations with high density. The results of the calculation of the hybrid effect of intensive energy flows on a cylindrical target are presented. The main gas-dynamic and radiative parameters of the compressed target plasma are found.
Keywords: computer simulation, magnetized target, mathematical modeling, numerical algorithm
1. Introduction
Magneto-inertial fusion (MIF) [1-5] is based on the heating of a plasma formation (the central part of the MIF target) consisting of deuterium (D) and tritium (T) by target shell implosion with high-energy laser beams (laser driver) and/or high-speed plasma jets (plasma liner). The heating process can be described qualitatively as follows: compression of the plasma target (to the stage of nuclear reactions) occurs after the passage of a strong shock wave through
Received July 17, 2019 Accepted November 15, 2019
This research is partially supported by the Ministry of Science and Higher Education of the Russian Federation (Project No. 0705-2020-0044).
Viktor V. Kuzenov
vik.kuzenov@gmail.com
Sergei V. Ryzhkov
svryzhkov@bmstu.ru
Andrei V. Starostin
starostinav@gmail.com
Bauman Moscow State Technical University
ul. 2-ya Baumanskaya 5, Moscow, 105005 Russia
the target, accompanied by the irreversible heating of the ablator material (provided that the total external pressure exceeds the thermal pressure of the MIF target shell material, and the speed of thermal or ionization waves is below the speed of the shock wave (SW)), acceleration of the target shell, the generation of this shock wave in a target fusion plasma and their cumulation on a geometric axis of the system.
This paper uses an element of a multilevel computational model (described in [6]), which allows analysis of the evolution of compression and energy release in the MIF target under a combined sequential action of intense laser radiation (in the endwise direction at the final stage of compression and in the direction perpendicular to the generatrix) and an external (relative to the target) system of pulsed jets.
One of the goals of this work is to assess the possibility of achieving "thermonuclear" parameters in the MIF plasma by means of relatively "weak" intense effects (combined scheme), as well as to estimate the possibility of creating compact neutron generators based on the combined effect scheme. In principle, such a scheme can form a neutron yield at the level of 1012-1016 neutrons per pulse.
2. Mathematical model and computer simulation
2.1. Mathematical model
Assuming that the supply of energy by laser radiation and by a system of pulsed jets is cylindrically symmetric and the shell of a single-layer cylindrical MIF target is thin, we develop a one-dimensional mathematical model of the compression and energy release processes in the MIF target. This model is based on the radiation plasma dynamics equations written in the central symmetric coordinate system: the Euler equations, the equation of transport of broadband self-radiation, the equation of magnetic induction, the equation of transport of laser radiation, the equations of state and laser radiation absorption coefficients, and determines the conditions of the thermonuclear reaction:
dp d(pu) _ ^ ^ u(v - 1)
dt dr
^ + = Fpi Fp = (2.1)
d(pu) d(pu2 + P) = dt dr r(n^Jr>
Fpu = ~(pu2 - fr = -clJx H]r,
d(pE) d(pEu + Pu + qs)
~~di~ -dr'-~ E 9r
qr = jr Er + jz Ez, = qe + qi + qlaz,
Fe = -(pEu + Pu) P = Pe + Pi,
2qt
Qbutti7') = — exp(-xw(r)5),
where t is time, r is the radial coordinate, p is density, u is the velocity along the coordinate r, P = P(p, e) is the static pressure, e is the specific internal energy, E = (e + u2/2) is the total energy of gas flow, F = (Fp, Fpu, FE) is the vector of the sources in the orthogonal coordinate system, Fp is the mass flow density, Fpu is the pulse flux density, FE is the energy flow density, q, qv are the total and spectral radiation fluxes, Te, Ti are the temperatures of plasma electrons
and ions (T = Te = Ti), xv is the spectral absorption coefficient, fr is the electromagnetic force, qr is the energy flow from the electromagnetic field, qe = —Ae grad Te, qi = —\i grad Ti, qlaz is the laser beam, Ae, Xi are the thermal conductivity coefficients of electrons and ions, jr is the current density, H(r) is the magnetic induction vector, Pe is the pressure of the electrons, Pi is the pressure of the ions, v = (1, 2) are the indices which are related to the cases of planar and axial symmetry, q^utt (r), qlaz are the densities of the laser radiation flux, respectively, in the end direction (in this case, it is assumed that the plasma parameters of the MIF target are homogeneous along the longitudinal direction, i.e., along the Z axis) and perpendicular to the single-layer cylindrical target.
The symmetry condition (on the axis of symmetry) is imposed as the boundary conditions of the Euler equations (2.1). On the outer boundary, nondisturbing conditions are imposed on
d2f
the flow leaving the computational domain: —-
dr2
= 0, where f = {p, u,v, e}.
r=i
It should be noted that the noncollinearity (Vne xVT = 0) of temperature VT and gradients of electron density Vne can be a source of a spontaneous magnetic field with components Bz or B^ under the influence of intense energy flows on the MIF target shells. It should also be noted that, if the specified inhomogeneity is present when the azimuthal angle p changes, then the spontaneous magnetic field for the target shape under consideration has a component Bz.
Needless to say that the disturbance (due to the uneven irradiation) of the intensity of the energy flow along the surface (along the z-axis) of the MIF target leads to the formation of nonzero gradients of electron concentration, Vne = (dne/dr, dne/dz) = 0, and of temperature, VT = (dT/dr, dT/dz) = 0, and thus may give rise to a spontaneous magnetic field with a component B^.
The assessment of the spontaneous magnetic field is based on the formula dB/dt & ck[Vne x VT]/ene & ckTx/(er2Laz). For the characteristic values of the plasma motion scale rLaz && 10 cm, plasma te^Qperature rTx && 1 x 106 K, flux density q0az & 1014 W/cm2 and the time of laser radiation exposure t = 10 ns the characteristic value of the magnetic field will be B ^ 1 MGs. Magnetic fields with a strength of the order of megagauss and above can significantly affect the physical processes associated with the transfer of energy by the process of thermal conductivity.
The magnetic induction equation for the field component B^, which takes into account the
continuity equation for density p, the law of conservation of form div(B) = 0 and the source of the spontaneous magnetic field, can be written as follows [7, 8]:
d{Bjp) + 1 dJrjuBjp) = c2 d {JdrB\ _ ck x dt ¡xJr dr Air^pJrdr \a dr J pe 6
drBv [Vlnf?? ) x VTl (dMnJdT d\n(ne) dT\
rdr ~J" lVln(ne)xVT]„-^ Qr dz- Qz Qry (2.2)
The magnetic induction equation for the field component Bz (external magnetic field) is written as
d{Bz/p) | 1 dJjvBjp) = c2 d fJrrJBz\ (hB± = dt ¡jlJ dr Air^ipJr dr \ a dr J' rdr ip'
The boundary conditions for solving the equation of magnetic induction (2.2) and (2.3) can be described as follows: on the axis of symmetry by the symmetry condition and on the outer
boundary by dH(r)/dr = 0. The electrical conductivity is determined by the Spitzer formula [8] taking into account the possible magnetization of the plasma. The transfer of broadband radiation can be considered using a multigroup diffusion approximation with boundary conditions, which are described in [9].
It should be noted that it is quite difficult to provide an energy-efficient mode of isentropic compression of the central part of MIF target substance: the time form of the laser pulse qlaz(t)\r=£ must be spelled out in a special way, the pulse duration must be consistent with the radius and mass of the MIF target shell. However, in this case, it is in principle possible to obtain the noticeable degrees of compression c = pf /pi, where pi and pf are the initial and final density values of the central part of MIF target substance, but the final compression temperature Tf is relatively small. It is important to note that the practical implementation of
MIF requires a high level of temperature Tf = 108 K. The fact that the "effective" adiabatic index y changes significantly (compared to the beginning of the process) during the compression process is important (the ratio T/p1-1 = const, y = const ceases to be observed), and in this case it is necessary to use the actual equation of state. Thus, to achieve this level of final temperatures Tf for a time that is much less than the time of the sound wave path through the uncompressed gas, the compression of the MIF target should be started with the help of an SW. Therefore, the hot-spot concept was applied in the work. In this case the substance of the MIF target is compressed in the first series of shock waves, and thermonuclear heating (ignition) is carried out in the last shock wave with large amplitude. For this purpose, it is possible to use a profiled (according to a special law) time pulse of the laser radiation flux density qlaz(t)\r=r In this case, the parameters of the laser radiation incident (along the r axis) on the MIF target can be found by solving the equation of laser radiation transfer [9, 10]
[q°iaz (t/t0 )n, t < t0, taking into account the corresponding boundary conditions qlaz \r=£ = <
I 0, tg < t < t(0 + At,
t , ¿0
where f qlaz(t)|r= dt f qlaz(t)\r= dt = 0.5 is the functional condition for finding the time o ' o
of pulse exacerbation t£ = = 7.17 ns (usually, an additional condition is also used
(t0 - t£)/t0 < 1), q0az = 1 x 1014 W/cm2, n = 2.8, t0 = 8.6 ns, At = 1.48 ns.
The value of the energy flux density q at the boundary of the computational domain
[0, t < t0, 0 in the pulse jet is determined by the equation q\r=i = < 0 where
{q0, t0 <t<t0 + At,
[1011, 1012] W/cm2. The possibility of creating pulsed plasma jets with these parameters follows from the results of [10]. The laser radiation flux coming to the end part of the target is
0, t < t0, Qlaz, t0 <t<t0
determined as follows: qt = < 1 q[az = 1 x 10 W/cm2.
The coefficients of electron and ion thermal conductivity Xe i in the case of magnetized plasma can be calculated using the formulas [8]. The calculation of the thermodynamic e(T, p), P(T, p) and optical xi(T, p) environment parameters included in this system of equations was conducted within the framework of local thermodynamic equilibrium approximation using the computer system ASTEROID, developed by RAS academician S. T. Surzhikov [11], the Thomas-Fermi model with quantum and exchange corrections and models of the average charge [12, 13].
2.2. Initial conditions
The initial intensity value H(r) of the "seed" magnetic field in a rarefied environment is a fraction of Tesla. The spectral flux and volume density of broadband radiation qv, Uv, laser radiation flux qlaz for r G [0, i] at the initial time t = 0 are equal to zero.
2.3. Computational algorithm
The numerical solutions of the one-dimensional equations of the plasma dynamics of MIF targets are based on the method of fractional steps and are described in [9]. Despite the one-dimensional nature of the problem, it places high demands on the numerical method used in its solution. First of all, the calculation scheme should possess improved dispersion and dissipative properties, be economical and algorithmically simple, have the property of monotonicity and approximate smooth solutions (preferably with the highest order of accuracy). These requirements are satisfied by the method of numerical solution of quasi-one-dimensional two-temperature one-liquid equations of gas dynamics, which is based on the method of fractional steps, consisting in this case of two steps [14].
The first fractional step takes into account gas-dynamic processes (these processes correspond to the "hyperbolic" part of the system of equations involved). In this case, the processes of radiation transfer and electromagnetic processes occurring in the devices of the MIF system are considered at the second fractional step. The mathematical formulation of the first fractional step and the solution of the "hyperbolic" part of the system of equations are based on the divergent form and formulated in the following two ways:
dU_ dF(U) ~dt +
= Fo
or
dU_ ~dt
= L(U ), L = —
dF{U)
+ Fo
Here the parameter £ can take one of the values of a set of values (r, z), the solution vector has the form U = (p, pu^, pE)T, the stream variable vector is written as mm, and the right part vector is F2 = (Fp, Fpu, Fe)T.
We reduce the vector form of the system of Euler equations to normal form with the extracted time derivative on the left-hand side dUj,/dt = L(Ui), where L is the right-hand side of the system of Euler equations which contains no time derivatives. The solution obtained at the previous time step is used as an initial approximation. Then the four-step version of the Runge-Kutta method is applied in the form of the following sequence of steps:
f(D =
F3 =
Û}0) + —rL (û\0)
F(2) =
f(4) =
Û}0) + ^L f ûU
' ff0) + AtL (F(3)
It is well known that this method of finding a solution JJi with respect to time makes it possible to handle one of the problems of numerically solving the Euler equations — the need to ensure the positivity of the required functions (if the solution is positive at time tn, then it remains positive at time tn+1). The first fractional step uses the divergent form of the Euler
equations:
dp dpui = dt d( p' d(pu^) djpul + P) dt + ~ ^ d(pE) djpEu g + ~df~+ ~ E'
This system of equations can be written in the vector form dU/dt + dF(U)/d£ = F2, where u^ = (u, v), the parameter £ can take one of the values of a set of values (r, 2), the solution
vector has the form U = (p, pu^, pE)T, the stream variable vector is written as mm, and the
right part vector is F2 = (Fp, Fpu, FE)T. Here (for the time fractional step t G [t, t + At/2]) we use a nonlinear quasi-monotone compact difference scheme of high order of accuracy:
9Uj | F (^+1/2)-^(^-1/2, .
~gf ^--A-= t2, = L4i—1/2 - 4i> 4i+i/2 -
Gas-dynamic parameters Un+1, Un refer to the centers of the computational cells, while the flows Fnt1/2 must be determined on the surface of these cells. To increase the order of approximation
R L
of the difference scheme, it is necessary to "restore" the gas-dynamic parameters Xi±'1/2 "to the right" (index R) and "to the left" (index L) from the boundaries of the calculated cells. Then any
reconstructed function Y(x), [x = {£}], £ G distributions of the form
h ^
2 ' 2
, is represented by piecewise polynomial
no = = * + vw je - & + ^ (^p) K - e,]2 +
+ a[C - &]3 + b№ - £i]4 + q[C - Ci]5 + di[( - C/ + ei[C - C/ = Fi + ¿F.
These piecewise polynomial distributions should be limited (to give them a monotonic form) by some function p(Y) — a limiter, which has the following form [15]:
<ß{Yi) = min (1, ^
\Yi - max(Yk-i/2, Ffc+i/2)r \Yi - min(Yk-1/2, Yfc+1/2)l ) ' The function Y(x) satisfies the conditions of smooth conjugation:
FT&-i) = Y-i, Fn ) = Y+i,
x) _ _ dF^t+1)
^ 2
and the condition of conservatism of the reconstructed function Y(x) is f = F(^).
2
The conditions of smooth conjugation mentioned above can be formulated as a system of
linear algebraic equations.
The spatial derivatives (dY/d£)ij included in the piecewise polynomial distributions Y(£) are calculated as follows: first, for the discrete function Yi the approximate value Fi of the first partial derivative of the spatial variable £ with the eighth order of accuracy should be determined. To this end, in each cell with number i for each recoverable value Yi ■ the nonmonotonicity index Ind(Y) is calculated:
Ind(Y)i U
- Yi+2J + - 3^-1 + i|3FM. - 4Yi_1j + Yt_2jl + 9
where the value d is a small parameter.
Next, we find the first derivative fi with respect to the variable £ using the usual approximate formula of the second order of accuracy and perform its "monotone restriction" on the grid fi = p(Yi)fi, fi = (Yi+1 — Yi-1)/(2A^)+O(A2), where A is the spatial grid step in the direction
of £. Then the approximate "monotonized" value Fi of the first partial derivative with respect
d A|
to the spatial variables £ with an approximate error at the level Fi = — + tjj^ + 0(A|) can be found using the formula (i.e., by solving a system of equations with a tridiagonal matrix)
Qi = (E + A2/30)fi, F = {(E + A2/6)-1Qi}i, where Af = fi+1 — f-v A2fi = fm — 2f + + 1, M is the unit operator. This formula is a symmetric finite difference of the sixth order of accuracy [16]:
In piecewise polynomial distributions atm, there are second-order spatial derivatives
d 2yA
—- I = si, which are denoted by the symbol s^ and calculated with the eighth order of d£ ' i accuracy [17]:
— (si+i + st_x) + s, = - 342^2 + <{Yi+i + Yi-1) +
+ 380A2(Ft+2 + ~ 6840A2 +
R L
Next, using the reconstructed function Y(£), gas-dynamic parameters Y^^ are "restored" "to the right" (index R) and "to the left" (index L) from the boundaries of the calculated cells. Then the "antidiffusion" correction of "recoverable" parameters Y(£) at the edges of the cell Y'R±ltL/2 is carried out [14].
From the theory of approximation of functions Y(£) by a truncated Taylor series (some polynomial) it is known [17] that oscillations of the interpolated function occur in the vicinity of discontinuities of the original function Y(£). However, the function Y(£) can be decomposed into a series of more general form (in the Lagrange-Burman series [15, 16]) by degrees of some function f (£ — £i). It is possible to choose a function f (£) in such way that the amplitude of the parasitic oscillations of the numerical solution in the area £i & 1 reduces. Let the main part of the reconstructed function Y(x), [x = {£}], £ G [—A^/2, A^/2], be decomposed into a series of
degrees of the function f (x) by the Lagrange -Burman expansion [18]. Then the reconstructed function Y(£) can be written as
no = mo = yt+ip(Xi) {/(* - opi + f2{i~il)P2 +
+ - &]3+bj\ - et]4 + C( - +( - +e\ -
'dY\ (d2Y\ (dY\ fd2f
dU,. \d£2 J, \dUi \d£
Pi = » P-2
I), " •(»: *(C
where p1, p2 are the first coefficients of function expansion Y(£) into a truncated Lagrange-Burman series [19]. It should be noted [18] that under condition f (£) = £ the coefficients p1, p2 turn into the coefficients of the usual Taylor series at the corresponding powers of the monomial (£-£i). In the area of the break (where the break indicator is li & 1) the decomposition in a Lagrange-Burman number is carried out in powers of monotonic functions:
/(£) = A^e ^ _ , /5 = 4. Ind(Y) + 6 • (1 - Ind(Y)), + 0.
For the reconstructed function Y(£) in the area li & 1 (near the strong discontinuities) the conditions of smooth conjugation and conservatism can be formulated as a system of linear algebraic equations.
The processes of radiation transfer and electromagnetic processes occurring in the devices of the MIF system are considered at the second fractional step. In this case, the one-dimensional equation is used to calculate the thermal processes of the area 0 ^ x ^ 1,0 ^ t ^ tlim:
dT d (~ dT \
To determine the magnetic field strength one can in general consider the equation
dH | 1 9J(vH) _ c2 d (JdH\ (t,_1}
M pj dC 47T^JdCKadU1
In this case, for the "hyperbolic" part of the equation, the above nonlinear quasi-monotone difference scheme of increased order of accuracy is used. Then, in accordance with the method of splitting by physical processes, a numerical method for solving the "parabolic" part of the equation should be described:
To construct a numerical solution of the magnetic field diffusion equation and the heat conduction equation using a monotonized difference scheme of high-order accuracy, we replace them with an equivalent system of two first-order equations:
_ ,dH dw dH 1 .
H(£,t = 0) = (£), H(£ = a, t) = (t), H(£ = b, t) = Vb(t).
2
For the numerical construction of the solution, let us introduce a grid wT x uh, uh = {xi = = iA£, i = 0, N}, ujt = {tn = nAt, i = 0, 1, ...}, where At and A£ are the steps in the time and spatial variables. Thus, a system of equations on a wT x uh grid can be represented in the following finite-difference form:
Hn+1 - Hn dwn+1 dHn+1 wn+1
C&+1/2.U Ai - dt ' +
dt ' dt k(ti+1/2» in)
0.
H + /2, )
Then this system of equations with respect to functions w, hn+1 = — ' (assuming
1/2' *n)At
i+1/2' in) and C (Cj+1/2, in) constant) can be written as
dWn+1 hn+1 hn
+
dt
dhn+1 wn+1
+ / = = o.
dt y a(ti+i/2» in)Ai
_ ^(^t+l/2' i'n) r ,
C (ti+1 /2 » in)
A solution within the calculation cell ti+1] can be bound in the following form:
wn+1 (0 = C, exp ((/ ^ai+l/2 nAt) + C2 exp (-(/ ^+1/2,,Ai), = " {pi exP (C/^+1/2 ,Ai) - C2 exp (-C/^+1/2 ,Ai)) + h'\0,
where C = t - ti, C G [0» AC = ti+1 - tj-The
constants C1, C2 can be found (here the functions w and H are assumed to be continuous) from the boundary conditions for the flow w\^=0 = wi, w\Az=ç. ^ = wi+1 and the
functions H\c=o = Hi, H\AC=ç+ ^ = Hi+1. Equating pairwise the resulting formulas for
the coefficients C1, C2, we obtain a system of algebraic relations between the functions w(t) and h(t):
" <+1 + 0 ( ] + = \^ai+1/2Ai )
= 2Q J (^¿+1 ~ Cj)/2 j ^n
v ^+1/2Ai j
(Kti - + 20 I ^ll^2 ] «+1 + = 0, V Vai+1/2Ai )
where the values h and w are given at the nodal points ti_RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2020, 16(2), 325 341
Let us write on the segments [£i-1/2, £i] and [£i, £i+1/2] the last equation of the system of algebraic relations between the functions W(£), h(£) in the following form:
- hr1 + 20 (<1/2 + <+1) = 0. ^ Ei. W•
Assuming that the relation w™+1 & wn+l1/2 & w™+11/2 holds, the equation for the flow
ira+1) can be written in a relatively simple form: = 4^/fljAi, Pi_1/2 = 4^/ai_1Ai,
a+1/2 ft . , , . ,ai-1/2
H(£i+l/2> t)k{tii+1/2' "g H(£i-l/2> t)k{tii-1/2' ^¿T
^ / £¿+1/2 ^ | ( |
If we do not take into account the influence of boundary conditions and consider the relation
\/ a At ^ (Ci+i — to be valid, then the solution t) can be approximated by the formula
^exp(-(/7,/2)2/KAi))
i VaiAt
Hence, in the first limiting situation (h/2)2 fa^t ^ x> follows the formula H(£i, tn +
+ At) & Hi (h/'2)/iTsi/aiAt, the flow limiter ai+1/2 = 1, ai_1/2 = 1. i
For another limiting situation (h/2)2 faiAt ^ 0, the limiter a of the flow w can be described by the formula
<W = * + [exp H^/h M))W+
oM/2 = l+Sign(ft/2-CT) [exp Hh/2f/iai_iAt))/M + l-sign(/,./2-CT)
To obtain a finite difference representation (in the case of the condition
wn+1 & w—1/2 &
wn+11/2) of the equation (C(£, t)dH/dt = -dw/d£) of the divergent energy conservation law, one can integrate it on the segment [£i, £i+1]:
H n+1 — Hn
tn+1) - tn+1) + 1+1/2At 1+1/2 / cxe, w = o.
Substituting into the balance ratio the formula obtained for the flux w(£i, tn+1), one can get the following three-point difference scheme:
A Hn+1 _ C Hn+1 I B Hn+1 = F Ai+1/2 Hi-1/2 Ci+1/2 Hi+1/2 + Bi+1/2 Hi+3/2 = Fi+1/2,
where the three-point difference scheme coefficients Ai+1/2 ; Ci+1/2; Bi+1/2 ; Fi+1/2 are determined by the formulas:
k(£i+3/2' trn)
ai+3/2
Bi+1/2
i+3/2
q I ¿¿+3/2 ¿»+1 ] g) / ^+1/2 ^
t J \ 2y^Âi
v a
a
in)
Ai+1/2
A
i—1/2
0 / ¿»+1/2 ¿A / - - : |
^ 2^/a~At J \2yfi-JSI.)
lie f \a'i+1/2 ite + \a'i+1/2 ?
^i+1/2' ln) o ^¿+1/2' ln) o il1 „
^ ^ _ _ P»+V2__!__P»+l/2__^ f C_
1+1/2 e ( ¿i+3/2 ~ 1 \ Q( ¿¿+1/2 0 / ¿¿+1/2 - ¿A Q( Zj ~ ¿¿-1/2 A ^ Ai
^ 2y^At J ^ 2v^At J ^ 2y^At J py^At)
Fi+l/2 = Hi+1/2 J ti
The magnetic field diffusion equation or the equation of thermal conductivity represents the system of ordinary differential equations (ODE) of the 1st order concerning a time derivative. These equations can be integrated with the second order of accuracy using a two-stage monotone version of the Runge-Kutta method:
H{n+1)* = H{n) + AtR(H{n)\ R(wn+1) = -
V J C (^i+1/2,tn)
H{n+1) _ Hin) H(H+1)* +AtR ('"+1)*) ~ 2 + 2 '
where R is the right-hand side of the system of equations.
There are many well-known schemes and equations [20, 21], but the authors have developed their own compact polynomial difference scheme and, for example, the broadband radiation transfer calculating method is considered on the basis of a multigroup diffusion approximation [9]. The time step At required for the integration of the above compact polynomial difference scheme is chosen using the Courant-Friedrichs-Lewy stability condition. The developed numerical method was verified by solving a number of sample (model) problems.
The "hyperbolic" (convective) part of the MIF targets computer model was tested using a one-dimensional version of the Riemann problem (Sod's problem) of the decay of an unstable gap of a given configuration. The results of the numerical solution of the one-dimensional version of the Riemann problem (for certain initial data, this problem is sometimes called Sod's problem) of the decay of an unstable discontinuity of a given configuration are given below. In this case, the integration region is a single segment x G [0, 1], which is divided into two parts (z1/2 = 1/2) and in each part at the initial time t = 0 its values of density p, pressure P and velocity of the gas inside the computational domain u = 0 are set. The calculation results are given for time t^ = 0.4. Thus, relative to the initial position of the gap z1/2 = 1/2 the left pL, PL and the
0.4 0.6
(a)
Fig. 1. Spatial distributions of pressure (a) and density (b) in the Soda problem (the number of computational cells is 350). The solid line is the exact solution.
right values pR, PR of thermodynamic parameters are given:
(p, u, P)L = (p = 1, u = 0, P = 1)L, (p, u, P)r = (p = 0.125, 0, P = 0.1)r,
f(P,u,P)L, z<z1/2,
(p, u, P) =
(p,u,P)r, z ^ zl/2,
z G [0, 1].
Comparison of the exact and approximate solutions (Fig. 1) has shown that the difference is no more than one per cent [9]. Verification calculations (with the aim of assessing the damping system of the reflected shock waves) was carried out for the experimental conditions implemented in one-diaphragm aerodynamic shock tube at IPMech RAS. The calculations have shown that the calculation error is within the experimental accuracy of the results and can reach the level of 10%. As additional test verification calculations, we considered the air flow around the wedge conjugate to the plate and also around the cone conjugate to the cylinder of the following incident flow parameters: pressure P = 2060 Pa, velocity V = 1860 m/s, temperature T = 223 K, and the Mach number M = 6. These results are also in good agreement with the above calculations (the relative error is 0.4%).
The "thermal" ("parabolic") part of the model was tested using some problems that allow accurate analytical solutions: the heating of a continuous environment filling a flat semibounded space r > 0, by a heat flow through the left fixed boundary r = 0 (the relative error is less then 1%).
3. Calculation results
The results of numerical modeling are discussed below. To describe the numerical results
t i
obtained, we use the notation Nfus = // QeFus(r, t) x 3.567 x 104 x 2nr dr dt [n/cm] to denote
0 0
the number of neutrons (per unit length of the MIF target) that have left the computational domain by time t.
The calculated area and the MIF target consist of a central part and one coaxial layer:
• the central part of the target (core radius Rc = 0.05 cm) is filled with a D-T mixture with the density p = 5 x 10_2 g/cm3 and the initial temperature equal to T = 297 K. It is
surrounded by a coaxial layer (outer radius Rc = 0.1 cm), consisting of the metal (Al) with the density p = 2.7 g/cm3 and the initial temperature equal to T = 297 K.
• the computational domain has the outer radius l = 0.2 cm.
The initial intensity value of the "seed" magnetic field in an environment is a fraction of T. The spectral flux and volume density of broadband radiation for r £ [0, I] at the initial time t = 0 are equal to zero. The spatial derivatives dT/dz, d ln(ne)/dz (calculated along the one-layer cylindrical target MIF) necessary to determine the source of the spontaneous magnetic field were calculated approximately as follows: dT/dz & (T(r = L/2, t) — T(r = 0, t))/L/2, dln(ne)/dz & (ln(ne)(r = L/2, t) — ln(ne)(r = 0, t))/L/2, where L is the length of the MIF target along its axis of symmetry (L = 0.2 cm), T(r = L/2, t) = T(r = 0, t)/2, lnne(r = = L/2, t) = lnne(r = 0, t)/2.
Let us introduce additional notation necessary for a correct discussion of the numerical results obtained. It is well known [22] that the intensity of the shock wave can be characterized by a dimensionless value called the amplitude of SW Z = (P2 — Pi)/pic2, where the indices 1, 2 correspond to the gas-dynamic parameters of the plasma before and after the front of SW. In addition, it is known that the most significant feature of the SW interaction with the geometric axis of the MIF target is the enhancement of SW after its reflection. It is convenient to characterize the process of SW reflection from the axis of the MIF target by the corresponding gain K = AP3/AP2, where AP3 = P3 — P1 is the excess pressure in the reflected SW and AP2 = P2 — P1 is the excess pressure in the falling SW. The velocity of the shock wave VD going through the target material can be estimated in the first approximation through the speed of sound VD ~ ^/7R0T/p, and the velocity of the evaporation wave Vs follows from the approximate formula Vs(t) & Lvp/qiaz or from the relation [23], where Lv is the heat of sublimation of the ablator material. The velocity of the "heat" wave front Vf is determined by the equation [24] V^ = rj/(r(3m + 2)), where r^ = y7xj; r is the duration of the laser pulse, % = 16AaTm/(3pCv) = aTm is the ratio of the radiant thermal diffusivity, a = 5.67 x 10"8 W/m2K4 is the Stefan-Boltzmann constant, A is the inverse of the absorption coefficient of the substance, and Cv is the heat capacity of a continuous medium.
From the results of [9, 25-27] and given calculations it follows that the "compression" of the MIF target is possible in several modes.
At relatively low values of the laser pulse flux density qiaz(t)\r=e a shock wave spreads through the target shell (assuming that the external pressure exceeds the thermal pressure of the material of the MIF ablator). As the intensity of the laser pulse flux density increases qiaz(t)|r= an evaporation wave is formed on the outer surface of the target, pushing the target to the axis of symmetry. In this case, in front of the ablator a SW which spreads along the central part of the target in the direction of its axis of symmetry is formed. If one continues to increase the intensity of the laser pulse flux density qiaz (t)|r=, then a "heat" wave occurs in the target shell material, which moves only the energy (but not the mass, the MIF ablator) and the "compression" of the target is noticeably weakened.
At the same time, the process of combined compression, relative to time t, can be divided into three characteristic stages, within which the specified modes of "compression" of the MIF target are realized:
• the "initial" stage of MIF target compression in the first series of shock waves,
• the stage of MIF target compression during the "aggravation" of the time impulse;
• the stage of limitation of plasma formation extension by the system of external pulsed jets.
The results of compression of the MIF target by a system of pulsed jets can be described in the form of several structural phases (the main phases are the phase of "collapse", "re-reflection", and "expansion" of the compressed and heated MIF target plasma). All these phases are explained by the hydrodynamic nature of the target compression and the movement of "thermal" waves, when the flow incident on the target deforms and accelerates the target plasma, generates in the plasma of the MIF target central part a system of shock and "thermal" waves, interacting with each other and with the axis of symmetry, leading to the compression of plasma and the "frozen" magnetic field in it.
To illustrate the "expansion" phase, the graphical relationships are shown in Figs. 2-5. At the "expansion" phase of the central part of the target, a laser beam with an intensity of qfaz = 5 x 1015 W/cm2 enters through the end part of the cylindrical target. For the figures below q0az = 5 x 1014 W/cm2.
The energy supplied to the plasma, which is formed during the three stages of "compression" of the MIF target, to the central part of the target leads to an increase in the plasma-dynamic parameters in the region of their maximum values (Fig. 2).
T (kK)
22000 20000 18000 16000 14000 12000
0.05 0.1 0.15 Target radius (cm) (a)
0.2
120 100 80 60 40 20 0
p (kg/m3)
0.05 0.1 0.15 Target radius (cm) (b)
0.2
Fig. 2. Distribution of temperature (a) and density (b) (the number of neutrons per unit length is Nfus = = 0.19 x 1015 n/cm) at time t = 10.2 ns.
u (km/s)
1.4E+08 1.2E+08 1E+08 8E+07
P (atm)
0.05 0.1 0.15 Target radius (cm) (a)
0.05 0.1 0.15 Target radius (cm) (b)
Fig. 3. Spatial distribution of velocity (a) and static pressure (b) for the number of neutrons per unit length Nfus =0.19 x 1015 n/cm at time t = 10.2 ns.
T (kK)
0.05 0.1 0.15 Target radius (cm) (a)
160 140 120 100 80 60 40 20 0
P (kg/m3)
0.05 0.1 0.15 Target radius (cm) (b)
0.2
Fig. 4. Dependence of temperature (a) and density (b) for the number of neutrons per unit length Nfus = = 0.208 x 1015 n/cm at time t = 10.3 ns.
7E+06 6E+06 5E+06 4E+06 3E+06 2E+06 1E+06 0
mag
(atm)
0.05 0.1 0.15 Target radius (cm) (a)
0.2
4E+08 3.5E+08 3E+08 2.5E+08 2E+08 1.5E+08 1E+08 5E+07 0
P (atm)
0.05 0.1 0.15 Target radius (cm) (b)
0.2
Fig. 5. Dependence of the magnetic (a) and static pressure (b) for the number of neutrons per unit length Nfus = 0.208 x 1015 n/cm at time t = 10.3 ns.
As a result of the heating process and intense expansion, a shock wave is formed whose velocity is directed toward the axis of symmetry of the MIF target (Figs. 2, 3). Before the shock wave, a "thermal wave" arises, which moves in the same direction as the shock wave. At time t = 10.3 ns (Figs. 4 and 5), these waves accumulate on the axis of symmetry. This process of interaction of the shock wave with the axis of symmetry is accompanied by a significant amplification of the reflected waves. At the next instant of time, the "thermal wave" propagates to the front SW.
4. Conclusions
In this paper we have presented an assessment of the plasma-dynamic parameters of a single-layer cylindrical MIF target for a combined scheme of the impact of intense energy flows. In this scheme, we carried out a sequential action of intense laser radiation (with a maximum radiation flux density q0az = 1014 W/cm2) along the forming single-layer cylindrical target (up to time t ^ t0), then, from time t0 < t, the external (relative to the target) system of
pulsed jets (with energy flux density 1012 W/cm2) and laser radiation (in the end direction) with intensity qfaz = 1 x 1015 W/cm2. The paper presents a new version of the method of numerical solution of one-dimensional equations of plasma dynamics. This method is based (for the "hyperbolic" part of the system of equations (2.1)) on a nonlinear quasi-monotonic compact polynomial difference scheme of high order of accuracy, and for the "parabolic" part of the system of equations on a new numerical method (which allows calculations in cases of intense discontinuities in transport coefficients) for solving the diffusion equation of the magnetic field and the thermal conductivity equation using a monotonized difference scheme of high order of accuracy. The paper analyzes all the stages of compression of the MIF target for the combined scheme of action and shows that the temperature of the central part of the target can reach at certain times the level of T = 250 million K. With the help of calculations, the magnitude of the magnetic field induced by the compression of the target is estimated. The possibility of creating neutron generators on the basis of a combined scheme of action (the number of neutrons per unit length Nfus = 1.03 x 1015 n/cm by the time of termination of the impact) is shown. We note that, regardless of the chosen type of MIF target, hydrodynamic instabilities can develop during the compression process, which prevent the achievement of optimal parameters of thermonuclear fuel. This paper does not take into account the influence of asymmetry in the calculation of the compression of the MIF target, which can lead to a difference in the observed neutron yield (in two orders of magnitude) from the predictions of one-dimensional calculations. It also follows from the above that analysis taking into account the inhomogeneity of physical processes and problems related to hydrodynamic instabilities (arising during the compression of the MIF target), of which the main role is played by Rayleigh - Taylor and Richtmyer - Meshkov instabilities, requires a further separate detailed study.
References
[1] Hsu, S. C. and Thio, Y. C.F., Physics Criteria for a Subscale Plasma Liner Experiment, J. Fusion Energ, 2018, vol.37, pp. 103-110.
[2] Aleksandrov, V. V., Branitski, A.V., Gasilov, V. A., Grabovskiy, E.V., Gritsuk, A.N., Mitro-fanov, K.N., Olkhovskaya, O.G., Sasorov, P. V., and Frolov, I.N., Study of Interaction between Plasma Flows and the Magnetic Field at the Implosion of Nested Wire Arrays, Plasma Phys. Control. Fusion, 2019, vol.61, no. 3, 035009, 16pp.
[3] Gomez, M.R., Slutz, S. A., Knapp, P.F., Hahn, K.D., Weis, M.R., Harding, E.C., Geissel, M., Fein, J.R., Glinsky, M.E., Hansen, S.B., Harvey-Thompson, A. J., Jennings, C. A., Smith, I. C., Woodbury, D., Ampleford, D.J., Awe, T.J., Chandler, G.A., Hess, M.H., Lamppa, D. C., Myers, C.E., Ruiz, C.L., Sefkow, A.B., Schwarz, J., Yager-Elorriaga, D.A., Jones, B., Porter, J.L., Peterson, K. J., Mcbride, R. D., Rochau, G. A., and Sinars, D. B., Assessing Stagnation Conditions and Identifying Trends in Magnetized Liner Inertial Fusion, IEEE Trans. Plasma Sci., 2019, vol.47, no. 5, pp. 2081-2101.
[4] Chirkov, A. Yu., Ryzhkov, S.V., Bagryansky, P. A., and Anikeev, A. V., Plasma Kinetics Models for Fusion Systems Based on the Axially-Symmetric Mirror Devices, Fusion Sci. Technol., 2011, vol.59 (1T), pp. 39-42.
[5] Kuzenov, V. V. and Ryzhkov, S. V., Calculation of Plasma Dynamic Parameters of the Magneto-Inertial Fusion Target with Combined Exposure, Phys. Plasmas, 2019, vol. 26, no. 9, 092704.
[6] Kuzenov, V. V., Ryzhkov, S. V., and Frolko, P. A., Numerical Simulation of the Coaxial MagnetoPlasma Accelerator and Non-Axisymmetric Radio Frequency Discharge, J. Phys. Conf. Ser., 2017, vol.830, 012049, 6pp.
[7] Kuzenov, V. V. and Ryzhkov, S. V., Numerical Simulation of the Effect of Laser Radiation on Matter in an External Magnetic Field, J. Phys. Conf. Ser., 2017, vol.830, 012124, 7pp.
[8] Braginskii, S. I., Transport Processes in a Plasma, in Reviews in Plasma Physics: Vol. 1, M. A. Leon-tovich (Ed.), New York: Consultants Bureau, 1965, p. 205.
[9] Kuzenov, V. V. and Ryzhkov, S. V., Numerical Modeling of Laser Target Compression in an External Magnetic Field, Math.. Models Comput. Simul, 2018, vol. 10, no. 2, pp. 255-264.
[10] Volosevich, P. P., Degtiarev, L. M., Levanov, E. I., Kurdiumov, S. P., Popov, Iu. P., Samarskii, A. A., and Favorskii, A. P., Ultrahigh Compression of Matter and Ignition of a Thermonuclear Reaction by an Intense Laser Pulse, Sov. J. Plasma Phys., 1976, vol.2, pp.491-498; see also: Fiz. Plazmy, 1976, vol. 2, pp. 883-897.
[11] Surzhikov, S. T., Computing System for Solving Radiative Gasdynamic Problems of Entry and ReEntry Space Vehicles, in Proc. of the 1st Internat. Workshop on Radiation of High Temperature Gases in Atmospheric Entry (Lisbon, Portugal, Oct 2003), B. Warmbein (Ed.), European Space Agency (Special Publication), vol.533, Noordwijk: ESA Publ. Div., 2003, pp. 111-118.
[12] Kuzenov, V. V., Ryzhkov, S. V., and Shumaev, V. V., Numerical Thermodynamic Analysis of Alloys for Plasma Electronics and Advanced Technologies, Probl. Atom. Sci. Tech., 2015, no. 4(98), pp. 53-56.
[13] Kuzenov, V. V., Ryzhkov, S. V., and Shumaev, V. V., Application of Thomas-Fermi Model to Evaluation of Thermodynamic Properties of Magnetized Plasma, Probl. Atom. Sci. Tech., 2015, no. 1(95), pp. 97-99.
[14] Xu, Z. and Shu, C.-W., Anti-Diffusive High Order WENO Schemes for Hamilton-Jacobi Equations, Methods Appl. Anal, 2005, vol. 12, no. 2, pp. 169-190.
[15] Special Functions Handbook, M. Abramovits, I.Stigan (Eds.), Moscow: Nauka, 1979 (Russian).
[16] Lavrentiev, M. A. and Shabat, B.V., Methods of the Theory of Function of Complex Variable, Moscow: Nauka, 1987 (Russian).
[17] Pinchukov, V. I., Numerical Simulation of Unsteady Flows with Transient Regimes, Comput. Math. Math. Phys., 2009, vol.49, no. 10, pp. 1765-1773; see also: Zh. Vychisl. Mat. Mat. Fiz., 2009, vol.49, no. 10, pp. 1844-1852.
[18] Vorozhtsov, E.V., Application of Lagrange-Burmann Expansions for the Numerical Integration of the Inviscid Gas Equations, Vychisl. Metody Programm., 2011, vol.12, no. 3, pp. 348-361 (Russian).
[19] Fabrikant, N.Ya., Aerodynamics: The General Course, Moscow: Nauka, 1964 (Russian).
[20] Cardona, R. and Miranda, E., On the Volume Elements of a Manifold with Transverse Zeroes, Regul. Chaotic Dyn, 2019, vol. 24, no. 2, pp. 187-197.
[21] Kudryashov, N. A. and Sinelshchikov, D.I., Special Solutions of a High-Order Equation for Waves in a Liquid with Gas Bubbles, Regul. Chaotic Dyn., 2014, vol. 19, no. 5, pp. 576-585.
[22] Ovsyannikov, L. V., Lectures on the Basics of Gas Dynamics, Moscow: Nauka, 1981 (Russian).
[23] Zhurovich, M.A., Zhitkova, O.A., Lebo, I.G., Mikhailov, Yu. A., Sklizkov, G.V., Starodub, A.N., and Tishkin, V. F., Smoothing of Ablation Pressure Nonuniformities in the Laser-Plasma Corona during Heating of Laser Fusion Targets, Quantum Electron., 2009, vol. 39, no. 6, pp. 531-536; see also: Kvantovaya Elektronika, 2009, vol. 39, no. 6, pp. 531-536.
[24] Zeldovich, Ya.B. and Raizer, Yu. P., Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena, Dover, 2002.
[25] Ryzhkov, S. V. and Kuzenov, V. V., Analysis of the Ideal Gas Flow over Body of Basic Geometrical Shape, Int. J. Heat Mass Transf, 2019, vol. 132, pp. 587-592.
[26] Romadanov, I., Smolyakov, A., Raitses, Ye., Kaganovich, I., Tian, T., and Ryzhkov, S., Structure of Nonlocal Gradient-Drift Instabilities in Hall E x B Discharges, Phys. Plasmas, 2016, vol. 23, no. 12, 122111, 69pp.
[27] Kuzenov, V. V. and Ryzhkov, S. V., Radiation-Hydrodynamic Modeling of the Contact Boundary of the Plasma Target Placed in an External Magnetic Field, Applied Physics, 2014, no. 3, pp. 26-30 (Russian).
