Aerodynamic Design of a Laval Nozzle for Real Gas Using Hodograph Method

We propose an approach for the design of the subsonic part of plane and axisymmetric Laval nozzles for real gases. The proposed approach is based on the hodograph method and allows one to solve the inverse design problem directly. Real gas effects are taken into consideration using the chemical equilibrium model. We present nozzle contours computed with the proposed method for a stoichiometric methane-air mixture. Results confirm that real gas effects have a strong influence on the nozzle shape. The described method can be used in the design of nozzles for rocket engines and for high-enthalpy wind tunnels.


Introduction
Supersonic nozzles are used in jet engines and in high-enthalpy wind tunnels for acceleration of a hot gas mixture in order to generate uniform supersonic flow.In the case of a rocket engine, the overall efficiency of the engine largely depends on the quality of the flow produced by the nozzle.In the case of a wind tunnel, flow uniformity directly influences the accuracy of experimental measurements.That is why in both cases there is a need for an accurate and robust method for the design of the nozzle.
Typically, it is desired to obtain a uniform supersonic flow with prescribed velocity and pressure values at the output cross-section of the nozzle.Sometimes, additional requirements should be satisfied, such as the absence of flow separations and shocks, limitations on the length of the nozzle, etc.Therefore, the design of a nozzle naturally leads to an inverse problem: to find the unknown nozzle shape knowing the flow parameters in some regions.It is worth mentioning that since transonic flow is very subtle and is prone to instability, such an inverse problem should be solved with very high accuracy.
If there is no separation in the boundary layer on the nozzle wall, then viscous effects have a negligible effect and the inviscid compressible gas model (Euler equations) can be used.In the supersonic region, the governing equations are of the hyperbolic type, and consequently, the inverse problem for plane and axisymmetric nozzles can be solved relatively easily by means of characteristic-based methods [1,2].The main difficulty is the solution of the inverse problem in the subsonic region where the governing equations are elliptic.
Apparently, the first approach for this problem was based on the solution of the Cauchy problem for an elliptic second-order equation for the stream function [3].Initial data (velocity distribution) were set on the axis of symmetry (or plane of symmetry), and then, a marching method in the orthogonal direction was applied.Such a Cauchy problem for an elliptic equation is the well-known Hadamard's example of an ill-posed problem.Thereby, in order to obtain a generalized solution, some regularizations were applied, which led to an unavoidable error in the solution.Besides, it has been noted that nozzles computed using this approach are very long.
Modern methods for the design of a Laval nozzle are mostly based on the solution of the direct problem, like the majority of methods for aerodynamic optimization.The direct problem is to find a flow in a nozzle with a given shape and given inflow conditions.In contrast with the inverse problem, the direct problem is well posed and in most cases can be solved numerically.Usually in such optimization methods, the nozzle shape is parameterized using any convenient form (splines, Bezier curves), and the problem is formulated as a functional minimization problem for some objective function, for instance the norm of the deviation of output flow parameters from the desired ones [1,4,5].
One of the main advantages of such methods is their versatility: one can use any appropriate optimization procedure equipped with any black-box solver for the direct problem.Nevertheless, in many cases, the direct problem should be solved many times (from hundreds to thousands), which leads to large computational effort.Apart from that, at intermediate steps of the optimization procedure, the solution of the direct problem may produce complicated flow structures containing, for instance, flow separations and instabilities.In such cases, numerical methods for the direct problem can still produce large errors.That is why there is a need for methods that directly solve the inverse problem.
In the present study, we extend the approach for aerodynamic optimization, described in [6], which is applicable to inviscid plane and axisymmetric potential compressible flows of a perfect gas with constant adiabatic index.This approach utilizes a transformation from the spatial (x, y) plane to the hodograph plane.After such a transformation, velocity components become independent variables, while space coordinates become unknowns.This technique was successfully applied to a variety of problems including the design of spike nozzle [2], a nozzle guide vane [7] and Laval nozzles for experimental facilities [6] in "TsNIIMash" (Central Research Institute of Machine Building).Experimental measurements in these supersonic wind tunnels show that deviations from required the Mach number at the outlet of nozzles do not exceed 1% and confirm the robustness of the method.The problem statement in the hodograph plane provides the monotonicity of the velocity magnitude along the nozzle wall.In accordance with boundary layer theory, this property guarantees the separation of the free wall [8,9], which in turn justifies the applicability of the inviscid model.Another advantage of this method is that it provides the opportunity to design nozzles with a plane sonic surface.This allows one to design subsonic and supersonic parts of the nozzle separately.
The jet in rocket nozzles consists of a hot multicomponent mixture, and the typical temperature in the subsonic part of a nozzle is about 3000-4000 K.The same is appliable for high-enthalpy wind tunnels.Under such conditions, the perfect gas model is not valid, and real gas properties should be taken into account.In this paper, the described approach is extended to real gases.The extension is based on the fact that assuming chemical equilibrium or frozen reactions of the multicomponent mixture has a two-parameter equation of state [10].Using this property, we obtain the general form of the equation for the stream function in the hodograph plane.To demonstrate the applicability of the extended approach, we use the two-parameter equation of state arising from the chemical equilibrium of the stoichiometric methane-air mixture.
The outline of the paper is the following: In Section 2, the main assumptions and governing equations are briefly described.Section 3 is devoted to the formulation or the inverse design problem in the hodograph plane.Section 4 contains short description of the chemical equilibrium model.Computational results and discussion are presented in Section 5.

Potential 2D Flows of Real Gases
We consider plane or axisymmetric steady potential flow of an inviscid compressible gas obeying an arbitrary two-parameter equation of state.We will use the following notation throughout this paper: β, the angle between the velocity vector and x-axis • M = V/a, the Mach number, a being the speed of sound, • e, internal energy, It was proven in [11] that flow potentiality implies that the entropy and stagnation enthalpy are constant, just like in the case of ideal gas with a constant adiabatic index: Using these first integrals, the system of Euler equations can be reduced to the system of two equations: The first equation is the continuity equation, and the second equation reflects the fact that the flow is irrotational.
Due to the existence of the first integrals, the density, Mach number and other parameters of potential flow are functions of the velocity magnitude only: I 0 and S 0 are uniquely determined from the inlet conditions.Equations ( 3) and ( 4) can be rewritten using natural coordinates φ and ψ for unknowns V and β [6]: The hodograph transformation is a mapping from (φ, ψ) or (x, y) to (V, β).Properties of this transformation had been established in [6].The change of variables can be performed easily with the help of expressions like: Similarly, all derivatives in the equations can be expressed via derivatives with respect to β, V and the Jacobian of hodograph mapping J = ∂(β, V)/∂(φ, ψ).This leads to the following equation: Here: Coordinates x, y are connected with ψ through the relations:

Inverse Problem Formulation in the Hodograph Plane
Since isolines of the ψ function correspond to stream lines, we can formulate the inverse nozzle design problem as a well-posed boundary problem for second-order Equation ( 8) in the hodograph plane [6].
In this paper, we utilize one of many possible formulations of the inverse problem (which is depicted schematically in Figure 1 due to its simplicity.We consider an infinitely long nozzle; uniform flow with velocity V 0 , which is parallel to the x axis, comes from infinity and approaches sonic velocity on line D 1 D 2 .As mentioned, such a scheme allows one to design subsonic and supersonic parts of the nozzle separately, which is very convenient.
We require that along the curvilinear part O 1 B of the nozzle wall, the velocity magnitude is constant and equal to V 0 .Consequently, O 1 and O 2 are mapped to the point O in the hodograph plane, and O 1 B is mapped to the vertical line segment.
Furthermore, we require that along the rectilinear part BC, the velocity magnitude is monotonically increased from V 0 to critical value V cr such that M(V cr ) = 1.Since β = const along BC, it is mapped to the horizontal line segment.Finally, we assume that along the CD 2 part, the velocity magnitude equals V cr , so it is mapped to the vertical line segment.Both ends of the straight sonic line are mapped to the point D.
As mentioned, the OBCD part corresponds to one stream line and the OD part to another, which leads to a straightforward statement of the boundary conditions.
The described flow scheme results in a very simple boundary problem for elliptic Equation ( 8) with degeneration on the CD boundary, which can be solved easily using, for instance, the five-point finite difference scheme on a uniform Cartesian grid.It should be noted that discontinuity in the boundary conditions does not complicate the numerical procedure, since corner points are not used in the difference scheme.

Chemical Equilibrium Model
At high temperatures, chemical reactions significantly change the chemical composition of the gas mixture and, consequently, all the thermodynamic properties.That is why reactions should be taken into account when designing a high-enthalpy nozzle.
In general case of chemical nonequilibrium additional equations for the concentrations of each component should be solved, which are coupled with the equations of motion: Fixing p and T locally and linearizing the right-hand side at the equilibrium concentrations C eq , for which F( C eq , T, p) = 0, we come to the following system: which has the solution: From this expression, it is clear that the relaxation time, e.g., the characteristic time for a mixture to reach the equilibrium, is determined by the smallest in absolute value eigenvalue of the Jacobian matrix: τ chem ≈ 1/|λ| min .In the particular case of a simple reaction, this estimate is reduced to τ chem ≈ 1/k r (T, p), where k r is the reaction rate constant.In the general case, τ chem depends both on the reaction rate constants of all individual backward and forward reactions and the concentrations of each component.
If this time is very small compared to characteristic time τ V associated with gas motion, we can assume that at each point, the mixture is in local chemical equilibrium.Another extreme case of frozen reaction corresponds to large values of the ratio τ chem /τ V .The advantage of both models is that the mixture can be described by the two-parameter equation of state, which is easily incorporated in the framework of Section 2.
Here, we give a simple estimate for τ chem based on the relaxation times of individual reactions.To demonstrate the applicability of the described hodograph-based method, we consider a stoichiometric methane-air mixture with 16 components: Equilibrium concentrations for T = 2000 K and 4000 K and p = 1 atmosphere are depicted in Figure 2. The relaxation times of dissociation-recombination reactions, which play the most important role in our case, lying below 10 −5 s for the considered temperature range [12].There are a few slow reactions, for instance the nitrogen oxidation reaction.The relaxation time of this reaction varies from one second at T = 2000 K to about 7 × 10 −7 s at 4000 K [12].However, the molar fractions of both NO and NO 2 are very small at low temperatures, and the effect of such reactions on the thermodynamic properties of the mixture is negligible.
The characteristic time of macroscopic motion in the subsonic part of a nozzle τ V > L/V cr , where L is the nozzle length and V cr is the speed of sound at the sonic surface.For L = 1 m, V cr = 1000 m/s, we obtain the estimate τ V = 10 −3 .Thus, in our case, τ V /τ chem > 10 2  1, and the chemical equilibrium model can be considered.We should note that for more complex mixtures, one should use more careful estimates based on the detailed kinetic scheme and reaction rate constants.
The requirement τ V τ chem is often not met in the supersonic part of the nozzle, if the density is small.Nevertheless, the standard approach is to use the equilibrium approximation in the subsonic part and the frozen flow approximation in the supersonic part of the nozzle [10].Since the frozen flow is a two-parameter gas, the proposed hodograph-based method can be applied to the supersonic part as well.The calculation of the equilibrium composition in a gas is performed by means of a numerical solution of the law of mass action.The chemical equilibrium conditions can be written in the form: where x i is the molar concentration of the i-th component, K pi the equilibrium constant of the i-th component and ν ij the stoichiometric coefficient.
With equations for the molar and mass concentration of components, the law of mass action can be formulated in the form of the material balance system of equations with respect to the molar concentrations of the components consisting of L linearly independent equations.When the system is solved, it is possible to determine N unknown molar concentrations x i .The system is solved by Newton's method.In this case, the accuracy of the solution depends on the accuracy of the equilibrium constants, which can be expressed through the free Gibbs energy.For the Gibbs energy, we use approximations from [13,14].The approximations have the following form: where T is the temperature.Therefore, the equilibrium constants can be expressed through the Gibbs energy from the following equations: where R A is the gas constant and m i q i (0) the tabular values of the molar heat of the reaction.From the obtained molar concentrations of the components, the molecular weight of the mixture can be calculated.The density of the mixture is calculated from the equation of state of an ideal gas, taking into account the molecular weight of the mixture.
To calculate the molar enthalpy, the following formula is used: where m k h k (0) is the tabular values of the molar enthalpy of formation of the k-th component.
The specific enthalpy of the mixture is found from the specific enthalpies of the components: The dependence of molar mass on T for the equilibrium stoichiometric methane-air mixture at p = 1 atm is shown in Figure 3.This dependence demonstrates that at high temperatures (>1800 K), real gas effects should be taken into account.The described procedure yields all thermodynamic parameters of the chemically-equilibrated mixture as functions of p and T.However, the coefficients in Equation ( 8) depend on enthalpy and entropy.That is why an intermediate procedure is necessary to connect these two parameterizations.
Let us assume that we want to compute some quantity at (I * , S * ).At first, we find the equations of entropy isolines S = const in the (p, T) plane: Initial condition T(p 0 ) = T 0 for this ODE, which follows from the inlet condition for the nozzle, defines one curve T(p), corresponding to some implicit value of entropy S * .After that, solving the equation: we obtain pair (p, T).It should be noted that Equation ( 12) is integrated only once for the given problem, yet ( 13) is solved for each value of V in each mesh node.

Results and Discussion
The boundary problem for Equation (8) depicted in Figure 1 is solved using the standard second order five-point finite-difference scheme on a Cartesian mesh.In the axisymmetric case, the nonlinear right-hand side term in Equation ( 8) is updated iteratively.
The color map of the solution ψ(V, β) for T 0 = 300 K, p 0 = 1 atm, M 0 = 0.05 is illustrated in Figure 4.The x axis corresponds to velocity magnitude V and the y axis to velocity angle β.Nozzle contours are computed from the solution by integrating Equation ( 10) along the right, upper and left boundaries (where ψ = 1) with the second order explicit Runge-Kutta method.Serial computations show that 200 nodes along each direction are sufficient to obtain the mesh independent solution.Computed contours for different values of inlet temperature and β max for p = 10 atm are plotted in Figure 5.We present contours only for the plane case as it is more convenient for comparison.
As expected, the nozzle shape strongly depends on T since the chemical composition significantly changes along the nozzle.In the case of T = 500 K, there are no reactions, and the computed contour coincides with that computed using the perfect gas model.It is worth noting that the proposed method allows one to compute contours with any β max including β max > π/2.
The presented contours demonstrate that the longest part of the contour is that on which the velocity magnitude is constant.
Two velocity profiles along the nozzle wall for T 0 = 3500, p 0 = 100 atm are depicted in Figure 6.As has been already mentioned, on of the advantages of the proposed method is that it provides monotone acceleration of the flow along the nozzle wall, which in turn guarantees the separation of the free flow.It can be seen from Figure 6 that acceleration can be very fast if the straight part of the nozzle contour is short.This drawback can be eliminated by using non-rectangular domains in the hodograph plane.In this case, we can provide a strictly increasing velocity profile.

Conclusions
In this paper, a novel approach for the design of the subsonic part of the axisymmetric and the plane Laval nozzle is presented.The approach is based on the hodograph method, developed earlier for the potential flow of perfect gases, and the chemical equilibrium model.The extension to real gases is straightforward since real gas properties, usually available in tabular form, are encapsulated in the coefficients of equations on the stream function.
The applicability of the proposed approach is demonstrated by the computation of plane nozzle contours for different values of inlet temperature and pressure.
In our future work, we will investigate the applicability of the proposed method using a comparison of the flow obtained from the inverse problem with the solution of the direct problem with computed contours.Such a comparison will allow us to estimate the influence of uncertainty in the input data on the flow in the nozzle.Furthermore, we will try to use polygonal domains in the hodograph plane in order to provide a strict increase of the velocity along the nozzle wall.

Figure 1 .
Figure 1.Hodograph transformation for the flow in the nozzle.

Figure 2 .
Figure 2. Molar fractions of the components for temperatures of 2000 K and 4000 K.

Figure 3 .
Figure 3. Molar mass of the methane-air mixture.

Figure 4 .
Figure 4. Color map of the solution of the boundary value problem in the hodograph plane.

Figure 6 .
Figure 6.Velocity magnitude along the nozzle wall.