Numerical simulation of nozzle flow based on Euler equation

In this paper, the Euler-based numerical simulation method is employed to investigate the flow characteristics inside the Laval nozzle. In detail, the central difference method is used to discretize the flow field of the nozzle. And two kinds of flow conditions are simulated numerically: the isentropic flow from subsonic to supersonic, and the isentropic flow at full subsonic speed. The results show that the velocity increases and then decreases slowly, reaching the peak point near the nozzle throat. The pressure decreases non-linearly, and the density increases at first and then drops, following by a sharp rise later. The shock wave oscillation also appears in the flow from subsonic to supersonic, and this is because the artificial viscous term is not enough in discrete solution. Then we analyze the law of change of physical quantity by changing the shape of the nozzle. The results demonstrate that the larger the nozzle’s curvature, the faster the change rate of a physical quantity is. And the velocity at the nozzle inlet approaches zero when the area at the inlet and outlet approaches infinity.


Introduction
Computational fluid dynamics (CFD) has become a new branch integrating the disciplines of fluid mechanics and been one of the three basic approaches that can be used to solve problems related to fluid dynamics and heat transfer, etc., as shown in Figure 1. It is not only related to mathematics but also computer science. Besides, it plays an important role in industrial applications and academic research, especially in the field of high-technology engineering areas of aeronautics, as shown in Figure 2 [1]. A nozzle is often a pipe or tube of varying cross-sectional area, and it can be used to direct or modify the flow of a fluid. Nozzles are frequently used to control the rate of mass flow, speed, direction, or the pressure of the stream that emerges from them [2]. A nozzle is a relatively simple device, just a specially shaped tube. It is applied to different scenarios in a variety of shapes and sizes, like aircraft or rocket engines. Laval nozzle is mainly used to achieve supersonic speed, such as use in supersonic jet engines, and is widely used in several steam turbines and rocket engine nozzles. Laval's convergent-divergent nozzle was first applied in a rocket engine by Robert Goddard, and it could be simply shown in Figure 3 [2]. This nozzle in a rocket engine is used to accelerate a hot, pressurized gas to a supersonic speed and shape the exhaust flow through expansion. Then the heat energy can be maximally converted into directed kinetic energy [3]. Rockets typically use a fixed convergent section followed by a fixed divergent section for the nozzle design [4]. There exist several studies related to nozzle flows. For instance, Xu and Huang [5] analysed the characteristic design parameters of the rocket nozzles of the future launch vehicle system by using the Two-Dimensional Kinetic Program. They concluded that only when the nozzle area ratio and nozzle geometric properties meet a specific relationship, the efficiency of the rocket engine can be maximized. Liu et al. [6] optimized the nozzle profile of the rocket engine to improve the nozzle thrust with the help of CFD methodology by using three different optimization algorithms, including successive quadratic programming methods, genetic algorithm, and interdigitation strategy. For the axisymmetric Laval nozzle, Chen et al. [7] compared numerical simulations with the experimental data by using Non-oscillatory and non-free parameter Dissipation difference scheme and found a good agreement. Similarly, Wang et al. [8] used the methodology of CFD to simulate flow fields of nozzle combined with theoretical analysis, which demonstrated that it is necessary to use different wall functions for different operating regimes of Laval nozzle.
As a good agreement between the experimental data and numerical simulation, we try to apply the Euler equations to simulate the flow field of the Laval nozzle because it is time-consuming and expensive if we do physical experiments. In this paper, nozzle flow is obtained by using the Euler method. In detail, the influence of nozzle shapes on the flow is investigated schematically in both flow regimes of from subsonic at the inlet to supersonic at the outlet and full subsonic flow.

Euler equations
Euler equation is the inviscid form of the Navier-Stokes equation. For solving the less viscous fluid flow, we can omit the high-order partial derivative in the actual calculation processes to simplify the calculation process and increase the speed of calculation. This is the main advantage of employing the Euler equation to solve the aerodynamics problems and the main reason for adopting it to simulate the nozzle flow in this paper.
In the expansion section of the Laval nozzle, we select a small control body with the length of dx. in Figure 4. The left end of the control body is the initial physical quantity U, p, e,  and the right end isU dU + , p dp ＋ . Also, the area change is dS . We assume that the above physical quantities are uniformly distributed on the cross-section, so volume fraction in equation (2)   When dx  0, the product of two or three differentials in equation (2) quickly approaches zero, which has little influence on the original equation. To simplify the operation, we omit these terms and get equation (3): Divide dx at both ends of equation (3)and dx approaches zero. So ordinary differential terms can be written as partial differential form. Finally, we get one of the governing equations: The momentum equation and energy equation can be obtained similarly. Finally, we get the following equations, which are called the generalized Euler equation.
To facilitate our numerical discretization and iterative solution, we also need to write the equation in the form of a conservation equation.
The general formula of the conservation equation is: whereU , F and J are matrices of 3 1  , they are where J is the source term. And from the thermodynamic equation where γ=1.4 is the specific heat ratio of air. In this paper, the time advance method will be used, so the partial derivation of time is obtained:

Discretized equations
We use the following governing equations for the object to solute the discretization equations. With the same number of grid nodes, it is known that the central difference scheme is more accurate than the upwind differential scheme [9]. So, we use the central difference to discrete the governing equations. Firstly, we expand the governing equations into the following forms, we can derive the continuity equation, momentum equation, and energy equation.
Secondly, by employing the central difference, the partial derivative of U, ρ, and P can be discretized respectively as: where i represents the number of a node. By putting the discrete scheme into the governing equation, we can get a polynomial equation in several. The updated value of U is written as: The updated value of  is expressed as: The updated value of P is listed below:

Assumptions and shape of Laval nozzle 2.3.1 Assumptions
We assume that the flow field through the Laval nozzle is a steady isentropic flow. The fluid at the nozzle inlet comes from a container called the plenum. In real life, the velocity of gas inside the plenum tends to zero. In theory, accelerates through a pressure difference at both ends and reaches supersonic speeds at the nozzle outlet to propel the unit forward. We assume that the crosssectional area at any point of the nozzle is S, and S is a function of x, x is the axial distance of the nozzle. In this paper, to simplify the physical situation further and to increase the iterative speed of the computation, the fluid flow in the nozzle is assumed to be a one-dimensional flow, which means the variable density, velocity of the flow field, and the pressure are all functions of x. Also, the distribution of the flow field variables is uniform at any section, which is helpful for us to analyze the dynamic change of the flow field variables from the axial direction and increase the readability of the numerical experimental results. In this paper, we study the relationship between flow field variables and nozzle shape. With the same boundary conditions, we will set the nozzle shapes as follows: ( ) 2

Shape of Laval nozzle
Since S is the locus of motion of the boundary points of the cross-section, and the above five boundaries are quadratic functions of x, the size of the opening is determined by the coefficient of the quadratic term, so we change the shape of the nozzle by changing the coefficient of the quadratic term. Meanwhile, in the above analytic formula, (x-1.5) 2 and the constant term 1 is used to satisfy the need of grid generation in the system. We set the domain of x in the grid to [0,3], to make the physical distribution of the nozzle flow field in the center of the grid.

Subsonic -supersonic isentropic flow
In Figure 6-10, the variation trend of velocity is firstly increased and then gradually decreased, and the velocity becomes unstable and oscillates near the nozzle outlet. However, the whole trend is still a acceleration process because the velocity at the outlet is faster compared with the initial velocity obviously. And the source of acceleration is the pressure difference between the inlet and outlet. The flow field pressure is a process of gradual decrease. However, the density of the flow field increases first and then decreases, reaching the maximum near the throat, and then decreases gradually in the nozzle expansion section. Here, we can use the continuity equation of ideal fluid to assist in qualitative analysis of flow field changes: where 1 1 , S v are the cross-sectional area at one point of the nozzle and the local velocity, 2 2 , S v are the cross-sectional area at another point of the nozzle and the corresponding local velocity.
In the contraction section of the nozzle, the sectional area of the nozzle decreases gradually along the direction of the flow field, so the velocity of the flow field increases gradually. Theoretically, the maximum velocity of the flow field can be obtained at the minimum sectional area of the nozzle. However, in the expansion section of the nozzle, the velocity of the flow field decreases gradually because the cross-sectional area of the nozzle increases along the direction of the flow field. And for ideal fluid, the velocity of the flow field is equal at the same large cross-sectional area in theory. However, the fluid studied in this paper is not only controlled and affected by the continuity equation. Besides the continuity equation, the momentum equation and energy equation also have regulating effects on various physical quantities of fluid micro-clusters. These three equations together constitute the Euler equation of viceless isentropic flow. This is also the main reason why the results obtained through numerical experiments are different from the theories.

Relationship between flow field speed and nozzle shape.
By comparing the variation of the flow velocity of five Laval nozzles with different shapes shown in Figure 11, it can be found that the greater the curvature at the throat, the higher the maximum velocity of the flow field. The velocity reaches its maximum value at the expansion section of the nozzle, but not strictly at the throat. Meanwhile, with the increase of throat curvature, the nozzle inlet and outlet area gradually increases, and in this process, the initial velocity gradually decreases. It can be predicted that when the nozzle inlet area tends to infinity, the initial velocity at the nozzle inlet tends to zero indefinitely.  Figure 12, the larger the curvature of the nozzle throat, the larger the cross-sectional area of inlet and outlet, and the more nonlinear the pressure change of flow field is. In the nozzle, the pressure is almost always on a downward trend. When the curvature at the nozzle throat is large enough, for example, when the nozzle cross-section function satisfies S(x)=1+5.1(x-1.5) 2 , and the pressure ratio between inlet and outlet is 1/0.016, the pressure curve even shows an increase first in the contraction section. Then a rapid decrease near the throat, and finally goes through a short process of first increase and then decrease.  [9] detected that the mass flow distribution along the nozzle became a horizontal line after 700 iteration time steps, indicating that the mass flow at all positions in the nozzle had converged to the same constant, which proved that in onedimensional isentropic steady nozzle flow, the mass flow rate is constant.
Therefore, the first partial differential term of the equation (5) is equal to 0, and the equation becomes: Since Equation (22) is constant, so equation (23) can be further simplified into the following equation: We use this equation to analyze the change in momentum. Since the coefficient term on the left side of the equation is constant, if the velocity of the flow field increases, the partial derivative of the velocity concerning x is positive, and the partial derivative of the pressure concerning x is negative, so the pressure decreases. Therefore, the velocity continues to increase, and the pressure continues to decrease in the interval [1,2]. Figure 13, the larger the curvature of the throat, the larger the maximum density. Density decreases rapidly near the larynx. It can be seen from the figure that the greater the curvature of the throat, the faster the density in the range of [1,2] near the throat changes. At the same time, for nozzles of different shapes, the greater the curvature at the throat, the larger the cross-sectional area at the inlet and outlet, the smaller the density vibration near the nozzle outlet, and the more stable the density value.  Figures 14-18 show that, for the full subsonic isentropic flow, the flow field velocity increases and then decreases and reaches the maximum value at the throat. Similarly, the flow field density increases and then decreases, but it grows up again near the outlet. However, changes of flow field pressure when the shape curvature at the nozzle throat are small tend to be linear. As shown in the above Figures, with the increasing curvature at the nozzle throat, the area at the nozzle inlet and outlet also increases, the velocity at the nozzle inlet gradually decreases, and the velocity at the nozzle throat gradually increases. It indicates that if the sectional area at the nozzle entrance approaches infinity, the velocity at the nozzle entrance will also approach infinitesimal. In combination with the similar conclusion obtained in the subsonic-supersonic isentropic flow in the previous section, it can be concluded that when the cross-sectional area at both ends of the nozzle is large enough, the nozzle inlet velocity tends to zero.

Full subsonic isentropic flow
Meanwhile, for the Laval nozzle, the maximum velocity of the flow field is larger with the greater the curvature at the throat, which is the same as the conclusion obtained in the subsonic to supersonic isentropic flow in the previous section. For different types of flows, the larger the pressure difference between the inlet and outlet is, the larger the maximum velocity of the flow field will be.

Analysis of the oscillation
3.3.1 Shock wave analysis In Section 3.1, all the variables in the flow field have large or small numerical oscillations near the outlet of the nozzle expansion section from the flow field variable curve. This oscillation occurs because the fluid flow is accelerated to supersonic speed at a specific pressure ratio between the inlet and outlet. This may cause a positive shock wave to form somewhere in the nozzle expansion section. The velocity, density, and pressure of the flow field will change before and after the shock wave. In general, in explosions, supersonic flows, there will be shock waves. The shock wave of an ideal gas has no thickness and is a discontinuity in the mathematical sense. However, the actual gas has viscosity and heat transfer property, making the shock wave of the actual gas, not a discontinuous "truncation", but a continuous type. Therefore, the actual shock wave has a thickness, though the value is only a certain multiple of the free path of gas molecules, but the process is still rapid. 15 higher partial derivative. In Navier-Stokes equations, the high order partial derivative usually represents the physical viscous dissipation effects of flow. Still, the high order partial derivative is obtained by Taylor expansion, no similar physical meaning, so often referred to as "numerical dissipation". These higher-order partial derivative coefficient due to the physical properties much like N-S equation viscous effect, so collectively referred to as "artificial viscosity" [9].

Numerical Dissipation and Artificial Viscosity
The flow studied in this paper is an inviscid flow. Still, the artificial viscosity due to numerical behavior cannot be avoided, which makes the artificial viscosity appear implicitly in our numerical solution. Artificial viscosity, however, although can reduce the precision of the solution, but can improve the stability of the solution, as the numerical oscillation that appears in section 3.1. It is just because of the implicit artificial viscosity in numerical solution is not enough, so the solution will still be unstable. To eliminate this unstable shock, we need to add more artificial viscosity to stabilize the solution, eventually getting a stable and smooth solution.

Conclusions and future work
By solving the Euler equation with the finite difference method, we obtained two different flows in Laval nozzles. In detail, we investigated the influence of the nozzle shapes on the flow field inside the nozzle. By comparing the numerical results of flows from the subsonic to supersonic isentropic and full subsonic isentropic flows, the following conclusions can be drawn: (1) in the flows from the subsonic to supersonic isentropic, the velocity increases firstly and then decreases. In the flow field area [0,3], the velocity reaches the maximum near the abscissa coordinate of 2. The pressure shows a trend of gradual decline, and the density also shows a trend of increase and then decrease before increasing again; (2) in the full subsonic isentropic flow, the velocity also shows a trend of increase at first and then decrease. In the flow field area [0,3], the velocity reaches its maximum value at the abscissa coordinate of 1.5, that is, at the throat. The flow field density at the throat is equal to the flow field density at the inlet and outlet, and the pressure shows a trend of gradual decline. The speed decreases the fastest near the throat where it also has the largest pressure gradient; (3) for Laval nozzles of different shapes, as the curvature at the nozzle throat increases, the cross-sectional area at the nozzle inlet and outlet also increases, the velocity at the nozzle inlet decreases and the maximum velocity of the nozzle increases. As the nozzle inlet cross-sectional area approaches infinity, the initial velocity of the nozzle also approaches zero. Meanwhile, the larger the curvature at the nozzle throat, the faster the pressure and density change around the throat.
In the flows from the subsonic to supersonic isentropic, we think that the oscillation may be caused by the inadequacy of the hidden artificial viscosity term caused by the central difference method. In the future, we will further improve the finite difference method to solve the oscillation problem.