Abstract Generated abstract
A numerical method is developed for calculating the propagation of blast waves from a point explosion in an inviscid perfect gas with counterpressure, for plane, cylindrical, and spherical symmetry and for various adiabatic exponents. The approach combines central asymptotic formulas with the method of integral relations, dividing the region between the central interval and the shock front into strips and reducing the gas-dynamic equations to a Cauchy problem for ordinary differential equations in a shock parameter. The resulting scheme is implemented for three strips on a BESM-2 computer and applied to explosions in air, with accuracy monitored by conservation of mass and energy and by an adiabaticity condition. The computed spherical results are compared with earlier finite-difference data and show practical adequacy of low-order strip approximations.
Full Text
HYDROMECHANICS
V. P. KOROBEINIKOV, P. I. CHUSHKIN
A METHOD FOR CALCULATING A POINT EXPLOSION IN GASES
(Presented by Academician L. I. Sedov on 6 VIII 1963)
Let us consider the problem of the propagation of blast waves in a quiescent unbounded gas in the presence of counterpressure, in the plane \((\nu = 1)\), cylindrical \((\nu = 2)\), and spherical \((\nu = 3)\) cases (for the formulation of the problem and the results of its investigation, see \((^{1-3})\)). In the case of a strong explosion, when the counterpressure is negligible, the problem under consideration has a self-similar analytical solution \((^1)\). To calculate the initial stage of an explosion with counterpressure, a linearization method is used. For solving the non-self-similar problem of an explosion with account taken of counterpressure, an approximate approach is possible, based on the introduction of interpolation formulas for the unknown functions, the unknown parameters in which are determined from special integral relations \((^{2,4})\). An exact solution of this problem can be obtained by numerical methods; for a spherically symmetric explosion in air at \(\gamma = 1.4\), such a solution was computed, for example, in \((^{5,6})\) by the method of finite differences.
In the present work, for a numerical solution of the problem of an explosion with counterpressure for \(\nu = 1, 2, 3\) and various \(\gamma\), a scheme of the method of integral relations \((^{7,8})\) will be applied.
The gas is assumed to be inviscid, non-heat-conducting, and to have a constant adiabatic exponent \(\gamma\). We take the equations of unsteady one-dimensional motion of a perfect gas in the form
\[ \frac{\partial \rho}{\partial t}+\nu \frac{\partial}{\partial \zeta}(\zeta \rho u)=0, \qquad \varepsilon=\frac{\rho v^{2}}{2}+\frac{p}{\gamma-1}, \]
\[ \frac{\partial \varepsilon}{\partial t} +\nu \frac{\partial}{\partial \zeta}[\zeta u(\varepsilon+p)]=0, \qquad \frac{\partial p^{1/\gamma}}{\partial t} +\nu \frac{\partial}{\partial \zeta}(\zeta u p^{1/\gamma})=0, \tag{1} \]
where \(\zeta = r^\nu\), \(u = v/r\); \(\rho, \varepsilon, p, v\) are, respectively, the density, energy, pressure, and velocity of the gas. The boundary conditions of the problem are the known conditions on the shock wave and the condition at the center \(v(0,t)=0\).
We introduce new dimensionless independent variables \(\xi=\zeta/\zeta_{\mathrm{в}}\), \(q=a_{\mathrm{н}}^{2}/c^{2}\) \((0\le \xi \le 1,\ 0\le q \le 1)\) and the unknown functions
\[ g=\rho/\rho_{\mathrm{н}},\qquad e=\varepsilon q/\rho_{\mathrm{н}},\qquad \psi=(pq/p_{\mathrm{н}})^{1/\gamma},\qquad \varphi=u r_{\mathrm{в}}/c. \tag{2} \]
The relation between the dimensionless time \(\tau=t/t^{0}\) and \(q\) is given by the differential equation
\[ \tau'=\sqrt{\frac{q}{\gamma}}\, \frac{r_{\mathrm{в}}}{r^{0}} \left( \frac{\delta'}{\delta}+\frac{1}{q} \right)\frac{1}{\nu}. \tag{3} \]
Here and everywhere below, a prime denotes differentiation with respect to \(q\); the subscript “н” refers to the undisturbed medium, and “в” to the shock wave. In addition, the following notation is used:
\[ r^{0}=\left(\frac{E_{0}}{\rho_{\mathrm{н}}}\right)^{1/\nu}, \qquad t^{0}=r^{0}\left(\frac{\rho_{\mathrm{н}}}{p_{\mathrm{н}}}\right)^{1/2}, \qquad \delta=\frac{\sigma_{\nu}\zeta_{\mathrm{в}}}{\nu q r^{0\nu}}, \]
\[ \sigma_{\nu}=2\pi(\nu-1)+(\nu-2)(\nu-3). \]
where \(a\) is the speed of sound; \(c=dr_{\mathrm{v}}/dt\) is the shock-wave velocity; \(E_0\) is the explosion energy, calculated for \(\nu=1\) per unit area of the charge, and for \(\nu=2\) per unit length of the charge.
For the sought functions \(g,e,\psi,\varphi\) the following asymptotic formulas hold:
\[ \begin{gathered} g=\frac{\gamma+1}{\gamma-1}\left(C\delta\psi_0^\gamma \xi\right)^{1/(\gamma-1)} +O\left(\xi^{(2d\nu-1)/\nu}\right), \qquad e=\frac{\psi_0^\gamma}{\gamma-1}+O\left(\xi^d\right), \\ \psi=\psi_0+O\left(\xi^d\right), \qquad \varphi=\frac{1}{\gamma}\left(1+q\delta^{-1}\delta'\right)^{-1} \left(1-\gamma q\psi_0^{-1}\psi_0'\right)+O\left(\xi^d\right), \\ d=\frac{\nu+2(\gamma-1)}{\nu(\gamma-1)}, \qquad C=\frac{\alpha\nu(\nu+2)^2(\gamma+1)^2}{8\gamma\sigma_\nu}, \end{gathered} \tag{4} \]
where \(\alpha\) is a constant from the self-similar solution.
In the integration domain \(0\leq \xi \leq 1,\ 0\leq q\leq 1\), we single out, bounded by the lines \(\xi=0\) and \(\xi=\xi_0(q)\), a central interval in which the asymptotic formulas (4) are applicable. To choose the value of \(\xi_0\), we fix the Lagrangian coordinate of a particle \(\eta=\eta_0=(r_{\mathrm{n}0}/r^0)^\nu\) (\(r_{\mathrm{n}}\) is the initial coordinate of the point) so that, in the initial stage of the explosion, the corresponding asymptotic formula (4) gives \(\psi_0\) with a specified accuracy. Then, using the equations and the conditions of the problem, we obtain for \(\xi_0\) the expression
\[ \xi_0= \frac{C^{-1}}{\delta\psi_0} \left[ \frac{\gamma\sigma_\nu\eta_0 C}{\nu q(\gamma+1)} \right]^{(\gamma-1)/\gamma}. \tag{5} \]
Let us note that other methods of choosing the central interval are also possible.
The solution in the domain lying between the boundary of the central interval \(\xi=\xi_0(q)\) and the shock wave \(\xi=1\) will be determined by the method of integral relations. In the \(n\)-th approximation we divide this domain into \(n\) strips, drawing \(n-1\) intermediate lines
\[ \xi_i(q)=\frac{n-i}{n}\,\xi_0(q)+\frac{i}{n} \qquad (i=1,2,\ldots,n-1). \]
We integrate system (1), transformed to the variables (2), with respect to \(\xi\) from each line \(\xi=\xi_i\) \((i=0,1,\ldots,n-1)\) to the shock wave \(\xi_n=1\). As a result we obtain \(3n\) independent integral relations.
If, in general form, the sought functions are denoted by \(\Omega_s\), with \(\Omega_1=g,\ \Omega_2=e\) and \(\Omega_3=\psi\), and, in addition, the values of \(\Omega_s\) on each line \(\xi=\xi_i\) are indicated by the corresponding index \(i\), then the system of integral relations, taking into account the boundary conditions at the shock wave and formulas (4), will have the form
\[ \begin{gathered} I'_{s0}+\Omega_{s0}\xi_0' +\varkappa_s\Omega'_{30} +\delta^{-1}\delta'X_s =q^{-1}Y_s, \\ I'_{si}+\Omega_{si}\xi_i' -\delta^{-1}\delta'W_{si} =\left(W_{si}+\beta_s I_{si}\right)q^{-1} \\ (i=1,2,\ldots,n-1;\ s=1,2,3), \end{gathered} \tag{6} \]
where
\[ I_{si}=\int_{\xi_i}^{1}\Omega_s\,d\xi, \qquad X_1=g_0\xi_0-1+I_{10}, \]
\[ X_2=e_0\xi_0-Y_2+I_{20}, \qquad X_3=\psi_0\xi_0-\psi_n g_n^{-1}+I_{30}, \]
\[ Y_1=\beta_3 g_0\xi_0-X_1, \qquad Y_2=q/(\gamma-1), \qquad Y_3=\beta_3(\psi_0\xi_0+I_{30})-X_3, \]
\[ W_{1i}=g_i(\varphi_i-1)\xi_i+1-I_{1i}, \qquad W_{2i}=e_i(\varphi_i-1)\xi_i+Y_2-I_{2i}+\psi_i^\gamma\varphi_i\xi_i, \]
\[ W_{3i}=\psi_i(\varphi_i-1)\xi_i+\psi_n g_n^{-1}-I_{3i}, \qquad e_0=\psi_0^\gamma/(\gamma-1), \]
\[ \varphi_i=\pm \left[ \left(e_i-\frac{\psi_i^\gamma}{\gamma-1}\right)\frac{2}{\gamma g_i} \right]^{1/2} \xi_i^{-1/\nu}, \qquad \varkappa_1=\frac{g_0\xi_0}{\psi_0}, \qquad \varkappa_2=\frac{\gamma\xi_0 e_0}{\psi_0}, \]
\[ \varkappa_3=\xi_0, \qquad \beta_1=0, \qquad \beta_2=1, \qquad \beta_3=1/\gamma, \]
where the quantities with subscript \(n\) are expressed in terms of \(q\) from the conditions at the shock wave.
We shall now approximate the integrand functions by Lagrange interpolation polynomials with interpolation nodes on the lines \(\xi=\xi_k\) \((k=0,1,\ldots,n)\). For the integral \(I_{sl}\) we obtain the expression
\[ I_{sl}=(1-\xi_0)\sum_{k=0}^{n} A_{lk}\Omega_{sk}, \qquad A_{lk}=\frac{(-1)^{\,n-k}}{k!(n-k)!\,n}\int_l^n \frac{x(x-1)\cdots(x-n)\,dx}{x-k}, \]
where \(A_{lk}\) are generalized Cotes coefficients.
Substituting these representations for the integrals into the system of integral relations (6), we arrive at an approximating system of \(3n\) ordinary differential equations in \(q\) for the unknown functions. We shall determine the quantity \(g_0\) not from this system, but from the corresponding asymptotic formula (4). Then, finally, to determine the \(3n-1\) unknown functions \(\psi_0\), \(\delta\), \(g_i\), \(e_i\), \(\psi_i\), we obtain a system of differential equations that can be written in the form
\[ (1-\xi_0)\left(A_{00}-\sum_{i=1}^{n-1}\omega_i A_{i0}\right)\Omega'_{m0} +K_m\xi'_0+x_m\Omega'_{30} +\left(X_m+\sum_{i=1}^{n-1}\omega_i W_{mi}\right)\frac{\delta'}{\delta} \]
\[ =\frac{1}{q}\left[Y_m-\sum_{i=1}^{n-1}\omega_i(W_{mi}+\beta_m I_{mi})\right] -(1-\xi_0)\left(A_{0n}-\sum_{i=1}^{n-1}\omega_i A_{in}\right)\Omega'_{mn}, \tag{7} \]
\[ \Omega'_{sj}=\frac{1}{1-\xi_0}\left[\Omega_{sj}\xi'_0+\frac{1}{\Delta_n}\sum_{i=1}^{n}\Delta_{ij}F_{si}\right] \]
\[ (m=2,3;\qquad s=1,2,3;\qquad j=1,2,\ldots,n-1). \]
Fig. 1
Here the following notation has been introduced:
\[ K_m=(1-A_{00})\Omega_{m0}-A_{0n}\Omega_{mn} +\sum_{i=1}^{n-1}\omega_i Z_{mi}, \]
\[ Z_{si}=A_{i0}\Omega_{s0}+A_{in}\Omega_{sn}-\frac{n-i}{n}\Omega_{si}, \]
\[ F_{si}=Z_{si}\xi'_0-(1-\xi_0)(A_{i0}\Omega'_{s0}+A_{in}\Omega'_{sn}) +\frac{\beta_s}{q}I_{si} +\left(\delta^{-1}\delta'+\frac{1}{q}\right)W_{si}; \]
\[ \omega_i=\frac{1}{\Delta_n}\sum_{j=1}^{n-1}A_{0j}\Delta_{ij},\qquad \Delta_n=|A_{ij}|, \]
\(\Delta_{ij}\) are the cofactors of the elements \(A_{ij}\) of the determinant \(\Delta_n\). The derivative \(\xi'_0\) entering (7) is eliminated by means of formula (5).
The system of differential equations (7), to which equation (3) is also adjoined, is integrated numerically on an electronic computer. The initial data for a certain \(q\) are prescribed from the linearized solution of the explosion problem \((^{2,3,9,10})\). After integration of system (7), the fields of the desired functions are readily computed from the adopted interpolation polynomials.
The solution of the problem is determined for different numbers of strips \(n\), which makes it possible to judge the accuracy of the computed solution. In practice, satisfactory accuracy is already obtained for \(n=3,4\). The numerical integration of the systems for such \(n\) is carried out comparatively simply, since here one has to solve a Cauchy problem with initial data. Overall control of the accuracy of the solution is effected by the fulfillment of the integral conservation laws of mass and energy, as well as of the adiabaticity condition in the particle with coordinate \(\eta_0\).
According to the scheme developed, for the case \(n=3\) a program was compiled for the BESM-2 electronic computer of the Computing Center of the USSR Academy of Sciences; it makes it possible to carry out calculations for \(\nu=1,2,3\) and various \(\gamma\). Using this program, an explosion in air \((\gamma=1.4)\) was computed in the plane, cylindrical, and spherical cases. The results of the calculations are shown in Fig. 1, where the dependence \(\tau(q)\) is plotted. There also, for \(\nu=3\), the data of the calculations \((^5)\) are plotted as points for comparison.
In conclusion, we express our sincere gratitude to V. P. Kardikov for valuable discussions and participation in the development of the initial versions of the computational scheme, and also to E. Bishimov for compiling the program and carrying out the calculations.
Mathematical Institute named after V. A. Steklov
Academy of Sciences of the USSR
Received
27 VI 1963
CITED LITERATURE
- L. I. Sedov, Dokl. Akad. Nauk SSSR, 52, No. 1 (1946).
- L. I. Sedov, Similarity and Dimensional Methods in Mechanics, 4th ed., Moscow, 1957.
- V. P. Korobeinikov, N. S. Melnikova, B. V. Ryazanov, Theory of Point Explosion, Moscow, 1961.
- V. P. Korobeinikov, Zh. Prikl. Mekh. Tekh. Fiz., No. 6, 45 (1962).
- D. E. Okhotsimskii, I. L. Kondrasheva, Z. P. Vlasova, R. K. Kazakova, Tr. Mat. Inst. im. V. A. Steklova AN SSSR, 50 (1957).
- H. Goldstine, J. Neumann, Comm. Pure and Appl. Math., 8, No. 2, 327 (1955).
- A. A. Dorodnitsyn, Tr. III Vsesoyuzn. matem. s”ezda, 1956, 4, 447 (1958).
- O. M. Belotserkovskii, P. I. Chushkin, Zh. Vychisl. Mat. Mat. Fiz., 2, No. 5, 731 (1962).
- N. S. Burnova (N. S. Melnikova), Investigation of the Problem of a Point Explosion, Dissertation, Moscow, 1953.
- V. P. Korobeinikov, P. I. Chushkin, Zh. Prikl. Mekh. Tekh. Fiz., No. 4, 48 (1963).