CC BY
4Russian Journal of Nonlinear Dynamics, 2019, vol. 15, no. 4, pp. 543-550. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd190413
MSC 2010: 93B18, 93B52
Mathematical Modeling of Plasma Dynamics for Processes in Capillary Discharges
A statement of the problem is presented and numerical modeling of plasma-gas-dynamic processes in the capillary discharge plume is performed. In the developed model, plasma dynamic processes in a capillary discharge are determined by the intensity, duration of plasma formation processes in the capillary discharge channel, and thermodynamic parameters in the surrounding gaseous medium. The spatial distribution of temperature, density and pressure, radial and longitudinal velocities of pulsed jets of several capillary discharge channels is presented.
Keywords: capillary discharge, numerical method, plasma dynamics
1. Mathematical model
Pulse capillary discharge is one of the relatively simple methods for producing plasma [1, 2]. It is well known that this type of discharge is characterized by a long, fairly stable plasma structure of the pulsed jet in the atmosphere. The pulsed capillary discharge — a plasma flow generator in the experiments — was an interelectrode insert from a textolite cylinder with a diameter of 150 mm and a height of 50 mm, with an axial hole that is a working channel of a capillary discharge, electrodes, and housing. The electrodes were made in the form of flat steel plates, one of which covered the capillary discharge channel on one side. The initial breakdown of the plasma was carried out using electric explosion inside the capillary of metal conductors. Aluminum, copper or lead were used as plasma-forming substances.
In this paper, we consider individual results of a numerical study of the processes of formation and development, as well as broadband radiation from large-scale toroidal plasma vortex
Received May 28, 2019 Accepted October 10, 2019
This research was financially supported by the Russian Ministry of Science and Higher Education (Project No. 13.5240.2017/8.9) and the Bauman Moscow State Technical University Target Program for 2018-2020.
Victor V. Kuzenov vik.kuzenov@gmail.com Sergey V. Ryzhkov svryzhkov@bmstu.ru
Bauman Moscow State Technical University ul. 2-ya Baumanskaya 5, Moscow, 105005 Russia
V. V. Kuzenov, S. V. Ryzhkov
structures (which form near the boundary of a pulsed jet) and a capillary discharge torch in the air atmosphere.
When carrying out these calculations, the capillary discharge channel was not considered, and the gas-dynamic parameters flowing into the flooded space of the pulsed plasma jet were estimated as follows. In a first approximation, the mathematical model of plasmodynamic processes inside a channel of a capillary discharge can be constructed from the condition that the energy deposited from a capacitive storage is emitted by an optically dense plasma, and all electrical energy is transferred to the thermal energy of the plasma, which expires with sound velocity through a section of a capillary discharge.
Plasma-dynamic processes occurring in a plasma torch of a capillary discharge can be determined using the system of equations of viscous two-temperature radiation plasma dynamics. The numerical method for the equations of plasma dynamics is based on the method of fractional steps. The numerical solution is also described in [3, 4].
2. Numerical algorithm
The mathematical formulation and solution of the "hyperbolic" part of the system of equations is based on the divergent form and is formulated in the following two ways:
dU_ dF(U) ~dt +
= Fo
or
dU_ ~dt
= L( U), L = -
dF{U)
+ Fo
where £ is the parameter that may take one of the values of the set of values (r, z), the solution vector is F = (p, pu^, pE)T, and the vector of the right side is F2 = (Fp, Fpu, FE)T.
We present the vector version of the system of Euler equations to normal form with the time derivative selected on the left dJJi/dt: dFi/dt = L( Fi), where L is the right side of the system of Euler equations that does not contain time derivatives. The solution obtained at the previous time step is used as the initial approximation. Then the four-step version of the Runge-Kutta method is implemented as the following sequence of steps:
F(D =
Ff) =
[7f) + A
F<2) =
Fi4) =
'u^ + fL^
' F(0) + AtL (F(3)
It is known that such a method of finding a solution JJi with respect to time t solves one of the problems of numerically solving the Euler equations — the need to ensure that the desired functions are positive (if at tn the solution is positive, then it remains positive at tn+1). In the first fractional step, the divergent form of the Euler equations of the following form is used:
dp dpu.
= FP
d (Pu?)
+
d (pu2 + P)
d(pE) dt
+
c)t c)£
d(pEu^ + Pu^)
= F,
pai
d£
=F
E ■
This system of equations can be written in the vector form: dF/dt + dF(F)/d£ = F2, where u^ = (u, v), the solution vector has the form F = (p, pu^, pE)T, and the vector of the right side
is represented as F2 = (Fp, Fpu, FE)T. For a fractional time step t G [t, t + At/2], a nonlinear quasi-monotonic compact difference scheme of a higher order of accuracy is used:
| m+i/2)-m-i/2) -p
gt + ^ - P 2i ^ - iSi-1/2 ~ Si+1/2 ~ SiJ-
The gas-dynamic parameters Uif+1, U™ relate to the centers of the calculated cells, while the flows Ff±1/2 must be determined on the surface of these cells. At the same time, in order to
increase the order of approximation of the difference scheme, the gas-dynamic parameters Y^^ should be "restored" to the right (index — R) and "left" (index — L) from the boundaries of
the calculated cells. Then any reconstructed function Y(x), [x = {{}], £ G represented by piecewise polynomial distributions of the form:
no = F?(0 = Yi + vw ('IT) e - + ^ (S-) K -
2 ' 2
is
dU ^ ^ 2! VdC2 + at\i - &]3+bt\i - e/ + ct[e - + dt\i - e/+et [e - e/ = y + sr.
These piecewise-polynomial distributions should be limited (to make them monotonous) to some function p(Y) — the limiter, which has the following form:
<P(Yù = min (1, \Yi - max(Yk)\/\Yi - max( Y—1/2 ' ^+1/2)1'
Y - min(Yfc)\/\Yi - min(Yfc-1/2, Yfc+V2)\) . The function Y(x) satisfies the conditions of smooth conjugation:
i?&-i) = Y-i, Fn(Ci+i) = Y+i,
1) _ ^ 1) _ y»
and the condition of conservatism of the reconstructed function Y(x):
^ 2
2
The above conditions for smooth conjugation can be formulated as a system of linear algebraic equations. The spatial derivatives (dY/d£)i j included in the piecewise-polynomial distributions of Y(£), are calculated as follows: first, for the discrete function Yi we determine the approximate value Fi of the first partial derivative in the spatial variable £ with the eighth order of accuracy. To do this, in each cell with the number i for each restored value Yi j the nonmonotonicity index is calculated Ind(Y):
r>i - yi+2,j + m+u - + -
5I - n+-2j + - + - + y._2i.I + 9
where d is the small parameter.
Next, we find the first derivative fi with respect to the variable £ by the usual approximation formula of the second order of accuracy and make its "monotonous restriction" on the grid: f = vY)fi, fi = (Yi+1 — Yi-1)/(2Ag) + O(A2), where A is the step of the spatial grid in the direction £. Then the approximate "monotonized" value Fi of the first partial derivative with
d A6
respect to spatial variables with an approximation error at the level Fi = — + tjjqq + 0(A|)
can be found using the formula of the form (i. e., by solving a system of equations with a tridi-agonal matrix): Qt = (E + A2/30)fi, Fi = {(E + A2/66)-1Qt}t, where Aft = fi+i - f-i, A2fi = fi+1 — 2fi + 1, M is the unit operator. This formula is a symmetric finite difference of the sixth order of accuracy:
5 • A2 \ T f 2 • A2 A2x -1
( b2y
In piecewise-polynomial distributions there are second-order spatial derivatives I "^pr which are calculated with the eighth order of accuracy [4]:
— (si+i + Si_i) + s, = - 342^2 + <{Yi+i + Yi-1) +
51 23
+ 380A2 (1*+2 + Yl~2) ~ 6840A2 (1*+3 + Yl~3)'
Then an "antidiffusion" correction of the "recoverable" parameters Y (£) is performed at the edges of the cell Y^^ [5, 6].
It is known from the theory of approximation of the functions Y(£) by a truncated Taylor series (some polynomial) that oscillations of the interpolated function occur in the neighborhood of discontinuities of the original function Y(£) (and in areas of large solution gradients). However, the function Y(£) can be expanded into a series of a more general form (into a Lagrange -Burman series) in powers of a certain function f (£ — £i). In this case, it is possible to choose the function f (£) so as to reduce the amplitude of the parasitic oscillations of the numerical solution in the region li & 1. Let the main part of the reconstructed function Y(x), [x = {£}], £ G [—Ag/2, Ag/2] be expanded in a series in powers of the function. f (x) by means of the Lagrange-Burman decomposition [7], then the reconstructed function Y(£) can be written as
Y(t) = = + m){m-qpi +f\ ;-:V> +
+ aM - Ci]3 + bM - Ci]4 + cJC - Ci]5 + dite - + eje - Ci] 'dY\ f82y\ f8y\ [d2f
Pi = "TTTx ) P-2 =
where p1, p2 are the first coefficients of the expansion of the function Y(C) into a truncated La-grange-Burman series [8, 9]. Note [7] that, under the condition f (C) = C, the coefficients p1, p2
pass into the coefficients of the ordinary Taylor series for the corresponding powers of the monomial (£ — In the discontinuity region the expansion li œ 1 in the Lagrange-Burman series is carried out in powers of the monotonic function:
HO = Aç th ^(e - j , /? = 4- Ind{Y) + 6 • (1 - Ind(Y)), + 0.
For the reconstructed function Y(£) in the region li œ 1 (near strong discontinuities), the conditions of smooth conjugation and conservatism can be formulated as a system of linear algebraic equations.
3. Results of mathematical modeling
Figures 1-5 show two-dimensional spatial distributions of temperature, density, velocity and pressure at time t = 58.2 /is and t = 94.6 /s for W0 = 2.7 kJ, dk = 10 mm, P^ = 1 atm. As we can see from the spatial distributions, the primary toroidal vortex structures are developed at the initial stage.
A particular case of vortex structures is a plasma toroidal vortex or a ring vortex [10-13]. One of the main and important properties of the toroidal vortex is that it travels in an unbounded medium before its decay long distances compared to a cloud (plasma, gas, liquid) of the same size as the vortex. Thus, the distance traveled by the toroidal vortex to its decay can reach a value zmax & (60 ^ 150)Rin depending on their initial parameters (Rin is the initial radius of the toroidal vortex).
The form of plasma retractable Al erosion vapors (with a longitudinal velocity ~ 7 km/s) into the air can conditionally be described as follows: the "head" part is a spherical surface, and the "side" part has a conical surface.
This form of Al erosion vapor plasma is close to the flow form in the "stagnant" (bottom) region, where as a result of dissipation of the kinetic energy of the flow, the volume of gas adjacent to the surface of erosive vapor Al loses its kinetic energy and the remaining kinetic
Z, cm
15
R, cm
Fig. 1. Temperature distribution in a pulsed jet at t = 58.2 /j,s: 1
— 11000
10500
10000
— 9500
— 9000
— 8500
— 8000
— 7500
_ 7000
— 6500
— 6000
— 5500
— 5000
— 4500
— 4000
— 3500
— 3000
— 2500
— 2000
— 1500
— 1000
— 5000
— 390.093
acceleration vortex area.
R, cm
Fig. 2. Spatial distribution of the longitudinal (projection on the Z axis) velocity of a pulsed jet of a capillary discharge at the time moment t = 94.6 ps (W0 = 2.7 kJ, dk = 10 mm, = 1 atm).
8000
6000
4000
2000
T, K
P, kg/m3
Fig. 3. Radial distribution of density and temperature in a pulsed jet of a capillary discharge passing through the center of the accelerating vortex at the time moment t = 94.6 ps.
energy is not enough to achieve its back surface and the flow stops. In the external potential flow, the pressure is restored (increases downstream) in accordance with the Bernoulli law, since there is no dissipation here. This pressure applies to the entire depth of the mixing layer.
As a result, a pressure minimum occurs in the flow zone from the separation point to the attachment, and a positive pressure gradient exists downstream from it (that is, a positive pressure gradient VP acts against the flow), which leads to the appearance of a return flow.
1500
1000
500
u, m/s
v, m/s
-500
-2000
Fig. 4. Distribution of radial and longitudinal velocities in a pulsed jet of a capillary discharge, passing through the center of the accelerating vortex at the moment of time t = 94.6 /s.
10000
8000
6000
4000
2000
0
T, K
P, atm
20
15
10
0
0 2 4 6 8 10 12 14 16 18 20
Z, cm
Fig. 5. The longitudinal distribution of temperature and pressure in a pulsed jet of a capillary discharge (passing through the axis of symmetry of the plasma formation) at the time t = 94.6 js.
4. Conclusion
The statement of the problem is presented and numerical modeling of plasma-gas-dynamic processes in the capillary discharge plume and other applications is performed as the development of Ref. [14-17]. In the developed model, plasma dynamic processes in a capillary discharge are determined by the intensity, duration of plasma formation processes in the capillary discharge channel, and thermodynamic parameters in the surrounding gaseous medium. Comparison of
the results of calculations in a single plume of a capillary discharge with known and available
experimental data was made. Their satisfactory compliance is noted. Calculations of pulsed jets
arising from adjacent channels of a high-current capillary discharge were carried out.
References
[1] Pashchina, A. S., Efimov, A. V., Chinnov, V. F. and Ageev, A. G., Specific features of the radial distributions of plasma parameters in the initial segment of a supersonic jet generated by a pulsed capillary discharge, Plasma Phys. Rep., 2017, vol.43, pp. 796-800.
[2] Poniaev, S. A., Reznikov, B.I., Kurakin, R. O., Popov, P. A., Sedov, A. I., Shustrov, Y. A. and Zhukov, B.G., Prospects of use of electromagnetic railgun as plasma thruster for spacecrafts, Acta Astronautica, 2018, vol. 150, pp. 92-96.
[3] Kuzenov, V. V., The usage of regular development of mathematical model of processes in the torch of the capillary category, Physical and chemical kinetics in gas dynamics, 2011, vol. 11 (Russian).
[4] 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, 092704.
[5] 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.
[6] Xu, Zh. and Shu, Ch.-W., Anti-Diffusive Finite Difference WENO Methods for Shallow Water with Transport of Pollutant, J. Comput. Math., 2006, vol.24, no. 3, pp. 239-251.
[7] Vorozhtsov, E. V., Application of Lagrange-Burman Expansions for Numerical Integration of Invis-cid Gas Equations, Vychisl. Metody Programm., 2011, vol. 12, pp. 348-361 (Russian).
[8] Ovsyannikov, L. V., Lectures on the Basics of Gas Dynamics, Moscow: Nauka, 1981 (Russian).
[9] Kuzenov, V. V. and Ryzhkov, S. V., Approximate Method for Calculating Convective Heat Flux on the Surface of Bodies of Simple Geometric Shapes, J. Phys. Conf. Ser, 2017, vol.815, no. 1, 012024, 8 pp.
[10] Zarubin, V. S., Kuvyrkin, G. N., and Savel'eva, I. Y., Radiative-Conductive Heat Transfer in a Spherical Cavity, High Temp., 2015, vol. 53, no. 2, pp. 234-239; see also: Teplofiz. Vys. Temp., 2015, vol. 53, no. 2, pp. 243-249.
[11] 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.
[12] 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.
[13] Varaksin, A. Yu., Air Tornado-Like Vortices: Mathematical Modeling, High Temp., 2017, vol.55, no. 2, pp. 286-309; see also: Teplofiz. Vys. Temp., 2017, vol. 55, no. 2, pp. 291-316.
[14] Kuzenov, V. V. and Ryzhkov, S.V., Radiation-Hydrodynamic Modeling of the Contact Boundary of the Plasma Target Placed in an External Magnetic Field, Prikl. Fiz., 2014, no. 3, pp. 26-30 (Russian).
[15] 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.
[16] Romadanov, I., Smolyakov, A., Raitses, Y., and et al., Structure of nonlocal gradient-drift instabilities in Hall ExB discharges, Physics of Plasmas, 2016, vol.23, pp. 122111.
[17] Ryzhkov, S. V. and Kuzenov, V. V., New Realization Method for Calculating Convective Heat Transfer near the Hypersonic Aircraft Surface, Z. Angew. Math. Phys., 2019, vol.70, no. 2, Art. 46, 9 pp.
