UDC 533.697.4

HYDROMECHANICS

Yu. N. D’YAKONOV, L. V. PCHELKINA

ON THE DIRECT PROBLEM FOR A LAVAL NOZZLE

(Presented by Academician G. I. Petrov on 18 VII 1969)

The system of equations of gas dynamics in the case of a nonviscous and non-heat-conducting gas can be represented in the form

\[ d\vec W/dt+\chi^{-1}a^2\operatorname{grad}\pi=0,\qquad d\pi/dt+\chi\operatorname{div}\vec W=0, \]

\[ d\pi/dt-\chi\,d\varepsilon/dt=0,\qquad \chi=\chi(p,\rho), \tag{1} \]

where \(\pi=\ln p,\ \varepsilon=\ln\rho,\ \chi=(d\pi/d\varepsilon)_s=a^2 e^{\varepsilon-\pi}\).

For a perfect gas \(\chi=\gamma=\mathrm{const}\); in calculations of flows with allowance for equilibrium physicochemical reactions, the effective isentropic exponent \(\chi\) and the other necessary thermodynamic functions will be determined with the aid of the approximation proposed in \((^1)\).

Consider the problem of unsteady and steady motions of gas in Laval nozzles of prescribed geometrical shape. We shall seek the solution in a region including the sub-, trans-, and supersonic portions of the flow.

Let us write system (1) in the form

\[ T\,\partial X/\partial t+A\,\partial X/\partial \xi +B\,\partial X/\partial \eta+C=0, \tag{2} \]

where \(X=[u,\ v,\ \pi,\ \varepsilon]'\) is the column vector of the unknown functions; \(u,\ v\) are the projections of the velocity on the \(x,\ y\) axes of the cylindrical coordinate system; \(T,\ A,\ B,\ C\) are coefficient matrices.

For simplicity let the variables be \(\xi=x/F,\ \eta=y/G(x)\) (\(F\) is the length, and \(y=G(x)\) is the equation of the nozzle contour).

We shall arrange the exit section \((\xi=1)\) in such a way that the inequality \(u>a\) is satisfied everywhere on it. Then, as is known, no additional conditions need be imposed on this boundary. The entrance to the nozzle \((\xi=0)\) will be discussed below. Suppose that at the instant \(t=0\) the initial state of the gas is prescribed, \(t_\infty\) denotes the final instant of time and, consequently, the region of interest to us is as follows: \(0<t\le t_\infty,\ 0\le \xi\le 1,\ 0\le \eta\le 1\).

To solve the problem in the indicated region we apply the grid method in the form in which it is set forth in \((^2)\). Denote

\[ t^n=t^{\,n-1}+\Delta t,\qquad \xi_m=\xi_{m-1}+\Delta\xi,\qquad \eta_l=\eta_{l-1}+\Delta\eta,\qquad t^0=\xi_0=\eta_0=0, \]

\[ t^N=t_\infty,\qquad \xi_M=\eta_L=1,\qquad f(t^n,\xi_m,\eta_l)=f^n_{m,l}, \]

\[ n=1,2,\ldots,N\qquad m=1,2,\ldots,M,\qquad l=1,2,\ldots,L. \]

Replacing the derivatives at the points \((n+\tfrac12,\ m+\tfrac12,\ l)\) by ratios analogous to those given in \((^3)\), and passing from system (2) to a system of difference equations, we write it in the form

\[ a^{\,n+\frac12}_{m+\frac12,l}X^{n+1}_{m+1,l} + b^{\,n+\frac12}_{m+\frac12,l}X^{n+1}_{m,l} = c^{\,n+\frac12}_{m+\frac12,l}. \tag{3} \]

We give a brief description of the numerical algorithm for computing the vector \(X^{n+}\) on the ray \(\eta_l=\mathrm{const}\). For simplicity we shall assume the existence of only one point \(m^*\) such that for all \(m<m^*\) the inequa-

inequality \(u_{m+1/2,l}^{\,n+1/2} < a_{m+1/2,l}^{\,n+1/2}\), which ensures the stability of the sweep (3), for all \(m \geqslant m^*\) one has \(u_{m+1/2,l}^{\,n+1/2} \geqslant a_{m+1/2,l}^{\,n+1/2}\) (here \(a\) is the speed of sound), and everywhere \(u \geqslant 0\). Applying the auxiliary difference scheme (2), for \(m=m^*\) we compute the initial sweep relation

\[ \mu_{m^*,l}^{\,n+1} X_{m^*,l}^{\,n+1}=g_{m^*,l}^{\,n+1}, \]

where \(g\) is a scalar, after which, with the aid of (3), we carry out a direct sweep for all \(m=m^*-1, m^*-2, m^*-3,\ldots,0\):

\[ \mu_{m,l}^{\,n+1}= \frac{\mu_{m+1,l}^{\,n+1}(a^{-1}b)_{m+1/2,l}^{\,n+1/2}} {\left|\mu_{m+1,l}^{\,n+1}(a^{-1}b)_{m+1/2,l}^{\,n+1/2}\right|}, \qquad g_{m,l}^{\,n+1}= \frac{\mu_{m+1,l}^{\,n+1}(a^{-1}c)_{m+1/2,l}^{\,n+1/2}-g_{m+1,l}^{\,n+1}} {\left|\mu_{m+1,l}^{\,n+1}(a^{-1}b)_{m+1/2,l}^{\,n+1/2}\right|}. \]

Next the gas-dynamic parameters are determined at \(m=0\) (see below). Knowing the vector \(X_{0,l}^{\,n+1}\), one can successively find \(X_{m,l}^{\,n+1}\) for \(m=1,2,\ldots,m^*\) from equations (3) and the relations

\[ \mu_{m,l}^{\,n+1}X_{m,l}^{\,n+1}=g_{m,l}^{\,n+1}, \]

forming from them a system of the form

\[ X_{m+1}=DX_m+E. \]

For \(m=m^*+1, m^*+2, m^*+3,\ldots,M\), the calculation is carried out by an analogous scheme constructed with the aid of (3).

The computations are performed for all \(l=0,1,2,\ldots,L\). Since the original system (1) is nonlinear, the method of iterations is used in solving the difference equations on a layer. In computing the vector \(X_{m,l}^{\,n+1}\) on the nozzle contour \((l=L)\), the no-penetration boundary condition is, naturally, used. Passing from one time layer to another, we find the flow field in the entire region \(0<t\leqslant t_\infty,\ 0\leqslant \xi\leqslant 1,\ 0\leqslant \eta\leqslant 1\).

Fig. 1. \(a\)—geometrical shape of the nozzles, sonic lines and parameter values on the wall of nozzle I according to the authors’ results; \(b\)—parameter values on the axis of nozzle I according to the authors’ results; \(v\)—geometrical shape of the sonic line and parameter values on the axis and wall for nozzle II according to the solution of the inverse stationary problem [4]. \(1—t=0;\ 2—t\approx 0.05t_\infty;\ 3—t=t_\infty\)

In the present work some computational results are presented, relating mainly to the established flow regimes. Nozzle I, shown in Fig. 1, was obtained as a result of solving the inverse stationary problem by the method of [4]. The equation of the contour of nozzle II is as follows: for \(0\leqslant x\leqslant 71.825\), \(y=G(x)=2.32\); for \(71.825<x\leqslant 75.425\),

\[ y=G(x)=1.66+0.66\cos\frac{\pi}{3.2}(x-71.825). \]

All linear dimensions are referred to the radius of the critical section. The figure shows the portions of the nozzles near the exit, and for convenience their origins of reference along the \(x\)-axis are shifted in comparison with the equations used in the calculation. The critical sections of both nozzles coincide. For contour I the distance from the inlet section to the throat is \(\lambda=9.897\).

Figure 1 shows the positions of the sonic lines and the distribution of the parameters along the axis and the contour, obtained by the authors as a result of solving the direct nonstationary problem. The data practically correspond to the established flow regime. Separately, in the left part of the figure, for nozzle II the sonic lines are plotted at three instants of time \(t\).

Figure 2 shows, in several sections \(x=\mathrm{const}\), the transverse fields of the gas-dynamic parameters. Pressures everywhere are referred to the value \(p^*\), velocities to the value \(\sqrt{p^*/\rho^*}\), where \(p^*, \rho^*\) are the total pressure and density; \(\gamma=1.4=\mathrm{const}\). For comparison, the results obtained for nozzle I by solving the inverse problem \((^4)\) and by one-dimensional theory are also presented.

Fig. 2. Solid lines—according to the authors’ results for nozzle I; dashed lines—according to the authors’ results for nozzle II; black dots—according to the results of solving the inverse stationary problem \((^4)\) for nozzle I; light circles—pressure values in the calculation according to one-dimensional theory

The initial state of the gas for nozzle I was determined from the calculation by method \((^4)\), and for nozzle II from one-dimensional theory, with the additional assumption of a linear law of variation of the inclination angle of the velocity vector between the contour and the axis in the sections \(x=\mathrm{const}\). In both cases, in the course of the solution, stationary values of the parameters were maintained at the inlet \((\xi=0)\), i.e.
\[ X_{0,l}^{n}=X_{0,l}^{0} \]
for all \(n=1,2,\ldots,N\). Since the main purpose of the flow calculation in contour I was comparison with the data of \((^4)\), the indicated procedure is natural.

For nozzle II, the conditions at \(\xi=0\) are essentially equivalent to conditions at infinity, since the length of the subsonic part was chosen in such a way that disturbances during the time \(t=t_\infty\) did not have time to reach, from the inlet, the section of the contour that interests us and is shown in Fig. 1. In a certain sense, the solution obtained imitates (with accuracy up to the initial data at \(t=0\)) the operation of a certain peculiar aerodynamic tube.

Let us note one more formulation of the problem under consideration, based on the model of a tank-type installation. Attach, in one way or another, a sufficiently large reservoir filled with gas to the contour under study. If \(G\gg g\), where \(G\) is the mass of gas contained in the vessel, and \(g\) is the mass of gas that has passed through the nozzle during the time \(t_\infty\) necessary for reaching the established regime, then as a result of the calculation at \(t=t_\infty\) we obtain the sought stationary solution. In this case the problem proves to be completely closed, and no additional conditions, apart from the boundary conditions of impermeability, are required.

Among the features of the flow pattern, we note a slight, approximately \(0.1\text{–}0.4\%\) of \(p^*\), increase in pressure on the wall (compared with its value at \(x=0\)), which takes place before the beginning of intensive acceleration of the gas. In Fig. 1 it is almost imperceptible because of the smallness of the scale.

In the regime of established flow, the surface of transition through the speed of sound for both nozzles differs substantially from a straight line. The change in its position with time for contour II is shown in Fig. 1.

The sonic line very quickly (in comparison with the flow as a whole) assumes a form close to stationary.

From the data presented it is clear that already for \(x \geqslant 0.6\) the flow differs from one-dimensional flow (especially in the vicinity of the critical section), and the character of the dependence of the parameters on the variable \(\eta\) deviates noticeably from a linear law. The insignificant difference between the results obtained in solving the direct nonstationary problem and the data of (4), amounting in the vicinity of the throat to a quantity of the order of \(1\%\), is explained mainly by certain inaccuracies in specifying the nozzle contour.

Moscow State University
named after M. V. Lomonosov

Received
9 VII 1969

REFERENCES

1. Yu. N. Dlyakonov, DAN, 157, No. 4 (1964).
2. Yu. N. Dlyakonov, Z. S. Pakhomova et al., Reports at the III All-Union Congress on Theoretical and Applied Mechanics, Abstracts of Reports, Moscow, 1968, p. 122.
3. K. I. Babenko, G. P. Voskresenskii et al., Spatial Flow Past Smooth Bodies by an Ideal Gas, Moscow, 1964.
4. U. G. Pirumov, DAN, 176, No. 2 (1967).