UDC 534.222.2

HYDROMECHANICS

V. P. KOROBEINIKOV

ON THE CALCULATION OF ONE-DIMENSIONAL FLOWS IN A CYLINDRICAL AND PLANE EXPLOSION IN AN IDEALLY CONDUCTING GAS, TAKING INTO ACCOUNT BACK PRESSURE AND A MAGNETIC FIELD

(Presented by Academician L. I. Sedov, 19 IV 1965)

1. Suppose that a cylindrical or plane point explosion has occurred in a stationary electrically conducting gas in the presence of a magnetic field \((^{1})\). We shall assume the conductivity of the gas to be infinite, and the gas to be perfect with a constant ratio of specific heats \(\gamma\).

Let \(E_0\) denote the energy released in the explosion, calculated per unit length in the cylindrical case \((\nu = 2)\) and per unit area in the plane case \((\nu = 1)\). Let the magnetic field be parallel to the axis of symmetry (the axis of the explosion) for a cylindrical explosion and parallel to the plane of symmetry (the plane of the explosion) for a plane explosion. In these cases the flow will be one-dimensional and will depend only on the time \(t\) and on one spatial coordinate \(r\)—the distance from the axis or plane of symmetry. As a result of the explosion, a magnetohydrodynamic shock wave will begin to propagate through the gas.

For the motion of the gas behind the front of the shock wave, we assume the equations of magnetohydrodynamics to be fulfilled:

\[ \frac{1}{\nu}\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \xi}(\xi \rho u)=0, \qquad \frac{1}{\nu}\frac{\partial \mathscr{E}}{\partial t} + \frac{\partial}{\partial \xi}\bigl[\xi u(\mathscr{E}+p^*)\bigr]=0, \]

\[ \frac{1}{\nu}\frac{\partial p^{1/\gamma}}{\partial t} + \frac{\partial}{\partial \xi}(\xi u p^{1/\gamma})=0, \qquad \frac{1}{\nu}\frac{\partial H}{\partial t} + \frac{\partial}{\partial \xi}(\xi uH)=0, \tag{1} \]

\[ \frac{1}{\nu}\frac{\partial}{\partial t}\frac{\rho u}{r^{\nu-2}} + \frac{\partial}{\partial \xi}(\rho v^2+p^*) + \frac{\nu-1}{\nu}\frac{\rho v^2}{\xi}=0, \]

where \(\xi=r^\nu,\ u=\dfrac{v}{r},\ \mathscr{E}=\dfrac{1}{2}\rho v^2+\dfrac{p}{(\gamma-1)}+\dfrac{1}{8\pi}H^2,\ v\) is the velocity, \(\rho\) the density, \(p^*=p+\dfrac{1}{8\pi}H^2,\ p\) the pressure.

From system (1) it follows (see, for example, \((^{2,3})\)) that it reduces to a system of equations of gas dynamics if \(p^*\) is regarded as the pressure and special thermodynamic relations are adopted. For \(\gamma=2\), equations (1) reduce to the usual equations of gas dynamics for a perfect gas.

The conditions at the shock waves have the form

\[ \rho_n(v_n-D)=-D\rho_\infty, \qquad \rho_n v_n(v_n-D)+p_n^*=p_\infty^*, \]

\[ (v_n-D)\mathscr{E}_n+p_n^*v_n=-D\mathscr{E}_\infty, \qquad H_n\rho_\infty=H_\infty\rho_n. \tag{2} \]

Here and below, quantities immediately behind the front of the shock wave are denoted by the subscript \(n\); the subscript \(\infty\) refers to the initial undisturbed parameters of the medium, with \(\rho_\infty=\mathrm{const},\ p_\infty=\mathrm{const},\ H_\infty=\mathrm{const}\).

In addition to conditions (2), as in the corresponding gas-dynamical problem \((^{1})\), we have the condition of symmetry of the flow

\[ v(0,t)=0. \tag{3} \]

To solve the problem of the propagation of the blast wave, it is necessary to integrate system (1) with conditions (2), (3).

We note that in system (1) not all equations are independent. One of the equations can be regarded as a consequence of the others, taking into account the relations between $\mathscr{E}$, $v$, $\rho$, $p$, $H$. The “extra” relations between the unknown functions can be used for an additional check on the accuracy of the solution.

2. Introduce new dimensionless variables

\[ q=\frac{\gamma p_{\infty}}{\rho_{\infty}D^2},\qquad \xi=\frac{\zeta}{\zeta_n},\qquad \psi=\left(\frac{q}{p_{\infty}}p\right)^{1/\gamma},\qquad \varphi=\frac{ur_n}{D}, \]

\[ g=\frac{\rho}{\rho_{\infty}},\qquad e=\frac{q}{p_{\infty}}\mathscr{E},\qquad h=\frac{q}{p_{\infty}}p^*,\qquad G=\frac{H}{H_{\infty}}. \tag{4} \]

In the variables (4), system (1) can be written as

\[ \frac{\partial g}{\partial q} +\mu\frac{\partial}{\partial \xi}(g\varphi\xi-g\xi) +\mu g=0, \]

\[ \frac{\partial\psi}{\partial q} +\mu\frac{\partial}{\partial \xi}(\psi\varphi\xi-\psi\xi) +\mu\psi-\frac{\psi}{\gamma q}=0, \]

\[ \frac{\partial e}{\partial q} +\mu\frac{\partial}{\partial \xi}(e\varphi\xi-e\xi+h\varphi\xi) +\mu e-\frac{e}{q}=0, \]

\[ \frac{\partial m}{\partial q} +\mu\left[ \frac{m}{\nu} +\frac{\partial}{\partial \xi} \left(g\varphi^2\xi^{2/\nu}+\frac{h}{\gamma}-m\xi\right) +\frac{\nu-1}{\nu}g\varphi^2\xi^{(2-\nu)/\nu} \right] -\frac{m}{2q}=0. \tag{5} \]

Since in the case under consideration system (1) has the integral

\[ G=g, \tag{6} \]

which satisfies the boundary conditions (2), then to solve the problem it is sufficient to find only the three functions $g$, $\psi$, $\varphi$.

From the kinematic condition $D=dr_n/dt$ there follows a relation between the dimensionless time $\tau=t/t^0$ and $q$,

\[ \tau'=R_n\frac{\mu}{\nu}\sqrt{\frac{q}{\gamma}}. \tag{7} \]

Here and below a prime will denote differentiation with respect to $q$.

In systems (5), (7) the following notation has been introduced:

\[ R=r/r^0,\qquad t^0=r^0(\rho_{\infty}/p_{\infty})^{1/2},\qquad r^0=(E_0/p_{\infty})^{1/\nu},\qquad \delta=\sigma_{\nu}R_n^{\nu}/\gamma q, \]

\[ m=\varphi g\xi^{(2-\nu)/\nu},\qquad \mu=\delta'/\delta+1/q,\qquad \sigma_{\nu}=2(\nu-1)\pi+(\nu-2)(\nu-3). \]

Conditions (2) can also be transformed to the new variables.

For the numerical solution of the problem we shall use a method analogous to that developed earlier [4, 5] for the problem of an explosion in gas dynamics.

Integrating system (5) with respect to $\xi$ from some $\xi_j=\xi_j(q)$ to $\xi_n=1$, taking into account the conditions at the shock wave, we obtain the system of integral relations:

\[ I'_{1j}+g_j\xi'_j-\mu\left[1+g_j\xi_j(\varphi_j-1)-I_{1j}\right]=0, \]

\[ I'_{2j}+e_j\xi'_j-\mu\left[ q\left(\frac{1}{\gamma-1}+\frac{\gamma}{2}\beta\right) +e_j\varphi_j\xi_j-e_j\xi_j+h_j\varphi_j\xi_j-I_{2j} \right] -\frac{1}{q}I_{2j}=0, \]

\[ I'_{3j}+\psi_j\xi'_j-\mu\left[ \psi_j\xi_j(\varphi_j-1)+\frac{\psi_n}{g_n}-I_{3j} \right] -\frac{1}{\gamma q}I_{3j}=0, \]

\[ I'_{4j}+m_j\xi'_j+\mu\left[ \frac{1}{\nu}I_{4j} +m_j\xi_j-g_j\varphi_j^2\xi_j^{2/\nu} -\frac{h_j}{\gamma} +\frac{\nu-1}{\nu}I_{5j} + \frac{q}{\gamma}\left(1+\frac{\gamma}{2}\beta\right) \right] -\frac{1}{2q}I_{4j}=0. \tag{8} \]

If we denote \(\Omega_1=g,\ \Omega_2=e,\ \Omega_3=\psi,\ \Omega_4=m\), then in system (8) we have

\[ \beta=\frac{H_\infty^2}{4\pi\gamma p_\infty},\qquad I_{sj}=\int_{\xi_j}^{1}\Omega_s\,d\xi,\qquad I_{5j}=\int_{\xi_j}^{1}m\varphi\,d\xi \quad (s=1,2,3,4). \]

Between the quantities entering system (8) there are the dependencies

\[ e_j=\frac{\gamma}{2}g_j\varphi_j^2\xi_j^{2/\nu} +\frac{\psi_j^\gamma}{\gamma-1} +\frac{\gamma\beta}{2}qg_j^2,\qquad m_j=g_j\varphi_j\xi_j^{(2-\nu)/\nu},\qquad h_j=\psi_j^\gamma+\frac{\gamma\beta}{2}qg_j^2. \tag{9} \]

Relations (7)—(9) are the starting point for constructing the method of computation.

1. As in papers \((^{4,5})\), in the integration region \(0\leq \xi \leq 1\) we introduce a central interval \(\xi_0(q)\), assuming that within the central interval the computation is to be carried out by special asymptotic formulas. The solution in the region \(\xi_0\leq \xi \leq 1\) will be determined using the scheme of the method of integral relations \((^6)\).

In the \(n\)-th approximation we divide this region into \(n\) strips, drawing \(n-1\) intermediate lines

\[ \xi_i(q)=\frac{n-i}{n}\xi_0+\frac{i}{n} \quad (i=1,2,\ldots,n-1). \]

If the integrands in (8) are approximated by Lagrange interpolation polynomials in \(\xi\), with interpolation nodes on the lines \(\xi=\xi_k\) \((k=0,1,\ldots,n)\), then one can obtain an approximating system of ordinary differential equations of first order with respect to the independent variable \(q\). For \(n>1\) the interpolation polynomials may have order \(n\) or lower; the chosen method of interpolation will affect only the coefficients of the quadrature formulas for the integrals entering (8). Assuming computations with low approximations, we shall take the degree of the polynomial equal to \(n\), and the interpolation “through,” with nodes at all points \(\xi_k\).

The approximating system can be solved for the derivatives of the unknown functions \(\Omega_{sj}\) and reduced to the form

\[ \Omega'_{sj}=F_{sj}(q,\Omega_{sj}) \quad (j=0,1,\ldots,n-1). \tag{10} \]

System (10), together with equation (7) and relations (6), (9), makes it possible to determine all the necessary quantities at the boundaries of the strips and the law of motion of the shock wave. Using the interpolation polynomials, one can find the fields of the unknown functions and other characteristics of the flow.

For the system of ordinary differential equations (7), (10), the Cauchy problem is posed if the solution is known in some small interval \(0\leq q\leq q_0\), which can be used to prescribe the initial data. The Cauchy problem is solved numerically using computers. General control of accuracy may be carried out by checking the satisfaction of the integral conservation laws of mass and energy, as well as the fulfillment of the additional relation (9) between \(e_j,\varphi_j,g_j,\psi_j\) when using in the computation the complete approximating system following from (8).

In what follows we shall assume that \(\beta<1,\ q_0\beta\ll 1\). This makes it possible, for prescribing the initial data at \(q=q_0\), to use the results of the computation of the linearized gas-dynamic problem with allowance for the back pressure (7). We choose \(\xi_0(q)\), fixing the Lagrangian coordinate of a particle at the instant \(q_0\). The computation of the unknown functions inside the zone of the central interval is carried out by the asymptotic formulas given in \((^4)\). We note that in this case the approximating system (10) will differ from the corre-

of the corresponding system for the gas-dynamics problem \((4,5)\) by the right-hand sides of the equations obtained from the second and last integral relations (8), and by the dependence of \(\Omega_{sn}\) on \(q\).

The computations were carried out in the fourth approximation \((n = 4)\), mainly for a cylindrical explosion \((\nu = 2)\) with \(\gamma = 1.4\), \(\gamma = 5/3\), and various \(\beta < 1\). For a plane explosion, the case \(\gamma = 1.4\), \(\beta = 0.2\) was computed. Some of the results of the computations are given in Figs. 1 and 2.

Fig. 1

Fig. 2

In Fig. 1 the dependence of \(\tau\) on \(R_n\) is shown, and in Fig. 2 the function \(p_0/p_n\) of \(\tau\) is given for the case \(\gamma = 1.4\), \(\nu = 2\) with \(\beta = 0\) and \(\beta = 0.2\), where \(p_0 = p(0,t)\), \(q < 0.77\).

A check of the accuracy by the integral law of conservation of mass gave, for \(q < 0.7\), a maximum error of no more than \(8\%\). In verifying the integral law of conservation of energy, the error did not exceed \(1.6\%\), and checking the calculation of \(\psi_j\) by the differential equations and by relation (9) gave a relative error not exceeding \(1.2\%\).

The computations showed that the accuracy of the calculation decreases as \(q\) increases. As \(q \to q_\infty = 1/(1+\beta)\), in order to improve the accuracy it is expedient to use the asymptotic laws of attenuation of shock waves. These laws can be obtained from the corresponding general formulas for the equations of gas dynamics (see, for example, \((8,9)\)), if one takes into account the corresponding dependence between the total pressure \(p^*\), the density \(\rho\), and the entropy \(S\).

For computing problems in the case \(\beta \geqslant 1\), it is desirable to solve accurately the linearized problem for a strong shock wave, taking into account the influence of the field \(H_\infty\). An approximate solution of this problem is considered in work \((10)\).

The author expresses sincere gratitude to E. Bishimov for assistance with programming, and to L. G. Straut and R. P. Agafonova for programming the problem for the BESM-2 computer and carrying out the computations.

Mathematical Institute named after V. A. Steklov
Academy of Sciences of the USSR

Received
7 IV 1965

CITED LITERATURE

1. L. I. Sedov, Methods of Similarity and Dimensionality in Mechanics, 4th ed., Moscow, 1957.
2. L. D. Landau, E. M. Lifshitz, Electrodynamics of Continuous Media, Moscow, 1957.
3. F. A. Baum, S. A. Kaplan, K. P. Stanyukovich, Introduction to Cosmic Gas Dynamics, Moscow, 1958.
4. V. P. Korobeinikov, P. I. Chushkin, DAN, 154, No. 3 (1964).
5. V. P. Karlickii, V. P. Korobeinikov, P. I. Chushkin, All-Union Congress on Theoretical and Applied Mechanics, Moscow, 1964, p. 103.
6. A. A. Dorodnitsyn, Proceedings of the Third All-Union Mathematical Congress, 1956, 4, 1958.
7. V. P. Korobeinikov, P. I. Chushkin, K. V. Sharovatova, Tables of Gas-Dynamic Functions for the Initial Stage of a Point Explosion, Publishing House of the Computing Center of the Academy of Sciences of the USSR, Moscow, 1963.
8. L. D. Landau, E. M. Lifshitz, Mechanics of Continuous Media, Moscow, 1953.
9. O. S. Ryzhov, Journal of Applied Mechanics and Technical Physics, No. 2 (1961).
10. J. D. Murrey, Mathematika, 8, pt. 2 (1961).