Real gas flow about a round cone

We propose a numerical method for a supersonic van der Waals gas flow about a round cone. This method enables one to calculate shock waves attached to the cone and flow parameters between the shock wave and the surface area of the cone. The advantages of the algorithm are demonstrated on some examples.


Introduction
We consider the supersonic van der Waals gas flow about a round cone. We assume that the upstream Mach number M ∞ > 1 is large enough and the shock wave is attached to the cone (see Figure 1).
The state equation of van der Waals gas has the form where p is the pressure, V is the molar volume of the gas, R is the universal gas constant, and T is the temperature. The constants V 0 and A provide the correct description of intermolecular forces. It is well known [1], [2] that the van der Waals gas can exist not only in the gaseous or liquid state but also in the two-phase state ( Figure 2, see also [3]).
We take e = c * V T 1 − as the equation for the internal energy. Here e is the internal energy, c * V is the heat capacity, A 1 is a constant. Usually, the equation of the internal energy is used with A 1 = 0, but the method proposed below enables one to use this equation in the more complicate form (1).

Equations and their reduction to a dimensionless form
The problem of a supersonic van der Waals gas flow about a round cone can be described by the system of PDE's of gas dynamics. But since the problem is axisymmetric this system is transformed into a system of ODE's. Following [4], we introduce some new constants for its reduction to a dimensionless form: Then the system in a dimensionless form reads (s = θ − a, 0 ≤ s ≤ b − a = δ, see Figure 1): Here we use the following notations: u(s), v(s), p(s), ρ(s) are the components of the dimensionless velocity, the dimensionless pressure and the dimensionless density respectively; It is convenient to perform the additional transformation In these terms system (2)-(5) becomes This system is supplemented with the condition of impermeability on the surface area of cone and the four conditions on the shock front which follow from the Rankine-Hugoniot conditions. Here k is the root of the function Q(k, b), 0 < k < 1 − β, 0 < β < 1:

Numeric method
Let M ∞ , a, α, β, γ ∞ and γ * be given constants. Because of the nonlinearity of the problem the proposed algorithm has an iterative structure.

Iterations
If the parameters k and b are fixed, then the boundary value problem (6)-(13) has a usual form and can be solved in any convenient way. For this propose we use a special form of equations (6)-(7). Assuming that the function M θ is known (in practice, we get M θ from the previous iteration), the boundary value problem (6)-(7), (10)-(11) can be solved explicitly: Note that in these terms the boundary condition (14) is reduced to the equation We then use the calculated functions U (s) and W (s) by solving numerically the Caushy plroblem (8)-(9), (12)-(13).
For calculating the values k and b for the next iteration we find the point of the intersection of the graphs of the functionk(b) defined by (16) and the implicit functionk(b) defined by (15).

First approximation
To built the first approximation for the solution we use the same fact connected with the intersection of curves (15) and (16) The typical examples of such functions (15) and (16) are given in Figure 3. A lot of numerical

Numerical examples
Let us put M ∞ = 5, γ ∞ = 1.8, γ * = 1.9, α = 4, β = 1/3, a = 10 • . We perform the algorithm described in sections 3.1, 3.2 with the discretization step 0.01 • . As the result we get the two Using the algorithm modification from section 3.3 we find pairs (b, M ∞ ) for several angles of the cone (see Figure 4). The minimal upstream Mach number for which the shock wave is attached to the cone can be defined by calculating of the minimums of the function M ∞ (b) (see Figure 5).
The parameter A 1 is not defined in the literature, and we thus put A 1 = 0. The influence of a nonzero value A 1 is demonstrated below as well.
The pressure and the density of the gas on the height 0 m arê p ∞ = 1.01325 10 5 Pa,ρ ∞ = 1.225 kg/m 3 .
We consider the case a = 3 • , M ∞ = 2. The discretization step is 0.01 • . As the result of modeling we get two shock waves with the angles b 1 = 30.0134 • (weak shock wave) and b 2 = 89.9240 • (strong shock wave). The solutions are showed in Figure 6. Other simulations show the following. The numerical solution corresponding to the height of 10000 m (p ∞ = 2.64999 10 4 Pa,ρ ∞ = 0.41351 kg/m 3 ) is slightly different from the described above. The increase of the Mach number decreases the angle of the weak shock wave (to b 1 = 7.6019 by M ∞ = 8) and increases the jump of the pressure and the density corresponding to the strong shock wave. A positive nonzero value A 1 decreases the angle of the strong shock wave (to b 2 = 75.9763 • by A 1 = 4 · 10 −6 ) and significantly changes the values of the pressure and the density (see Figure 7).