Numerical Simulation of Nozzle Flow Based On Euler Equation: Journal of Physics: Conference Series

Download as pdf or txt
Download as pdf or txt
You are on page 1of 17

Journal of Physics: Conference Series

PAPER • OPEN ACCESS You may also like


- Thermal stress analysis and structural
Numerical simulation of nozzle flow based on optimization of ladle nozzle based on finite
element simulation
Euler equation Zichao Rong, Jianhong Yi, Fengxian Li et
al.

- Effect of chevron nozzle penetration on


To cite this article: Qiaoxin Li et al 2021 J. Phys.: Conf. Ser. 2012 012078 aero-acoustic characteristics of jet at
M = 0.8
S R Nikam and S D Sharma

- 2D numerical modeling for plasma-


assisted CO2 pooling in supersonic
View the article online for updates and enhancements. nozzles: importance of a proper nozzle
contour design
Maryam Khaji, Kim Peerenboom, Joost
van der Mullen et al.

This content was downloaded from IP address 153.33.36.177 on 07/08/2022 at 01:25


ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

Numerical simulation of nozzle flow based on Euler equation

Qiaoxin Li1, a, *, †,Xu Qin2, b, *, †, Ruixin Qu3, c, *, †,Jiaqi Wang4, d, *, †


1
Department of Environment, Hohai University, Nanjing, 211100, China;
2
Naval Architecture and Ocean Engineering College, Jiangsu University of Science
and Technology, Jiangsu, Zhenjiang, 212100, China;
3
Mechanical engineering department, College of Beijing Institute of Petrochemical
Technology, Beijing, 102627, China;
4
Oceanography and Atmospheric Sciences College, Ocean University of China,
Qingdao, 266000, China.
*
Corresponding email: a1914050139@hhu.edu.cn, b182210101322@stu.just.edu.cn,
c
2018310603@bipt.edu.cn, dwangjiaqi4462@stu.ouc.edu.cn

These authors contributed equally.

Abstract. 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.

1. 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].

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution
of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Published under licence by IOP Publishing Ltd 1
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

Figure 1. The different disciplines contained within computational fluid dynamics

Figure 2. The three basic approaches to solving problems in a fluid dynamic and heat transfer

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].

Figure 3. Diagram of a De-Laval nozzle, showing approximate flow velocity (v), together with the
effect on temperature (t) and pressure (p)

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

2
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

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.

2. Methodology

2.1. 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.
 (1)

 dV +  U  dS  0
t

V S

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
is U + dU , p + dp , e + de , +d  . 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) is obtained as follows:

 ρU  dS = -ρUS + ρUS + USdρ + ρUdS + (2)


S

ρSdV + ρdUdS + UdSdρ + SdUdρ + dρdUdS


Figure 4. The control body used to derive the Euler equation

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):
(3)

ρUdS = ρUdS + ρSdU + USdρ = d (ρSU )
S

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:
 (  S )  ( US ) (4)
+ =0
t x
The momentum equation and energy equation can be obtained similarly. Finally, we get the
following equations, which are called the generalized Euler equation.
 ( US )  ( U 2 + p ) S dS (5)
+ = p
t x dx

3
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

 (  ES )  ( UHS ) (6)
+ =0
t x
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:
U F
+ J
(7)
t x
where U , F and J are matrices of 3 1 , they are
 S    uS 
    (8)
U =  uS ,F =  (  u 2 + p) S 
 
  ES    uHS 
   
where J is the source term. And from the thermodynamic equation
m2 (9)
p  ( -1)( - )
2p
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:
U
J-
F (10)
t x

2.2. 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.
(  S ) U S  (11)
+ S + U + SU 0
t x x x
U U  T (12)
 + U  -R (  +T )
t x x x
e U U  (ln S )
+  -  RT ( +U ) (13)
t x x x
Secondly, by employing the central difference, the partial derivative of U, ρ, and P can be
discretized respectively as:
U (i +1) - U (i - 1)
U x (i )  2 (14)
dx
 (i +1) -  (i - 1)
 x (i )  2 (15)
dx
p (i +1) - p (i - 1)
p x (i )  2 (16)
dx
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:
1 (17)
U next (i)  U - dt  (U (i) U x (i) +  Px (i)
 (i)
The updated value of  is expressed as:
U (i )
 next (i )   (i ) - dt  (U (i )   x (i ) +  (i )  U x (i ) +  (i )   S x (i  dx) (18)
S (i  dx)
The updated value of P is listed below:

4
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

P (i )
Pnext (i )  P (i ) - dt  (U (i )  Px (i ) +  (i )     U x (i )
 (i ) (19)
P (i )

 (i )
+  (i )  U (i )   S x (i  dx )
S (i  dx )

2.3. 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 cross-
sectional 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.

2.3.2 Shape of Laval nozzle In one-dimensional scenarios, because the space variable is x only, the
shape of the nozzle can be simulated by setting a function about x. The shape of the nozzle is made up
of the motion of the points on the cross-sectional area boundary, so the S(x) we set is essentially the
trajectory of this boundary point. For example, we set S(x)=1+0.3(x-1.5)2, we get the Laval nozzle and
its distribution as shown Fig. 5:

Figure 5. Nozzle shape


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:
S ( x )1  1 + 0.3  x  1.5 
2

S ( x ) 2  1 + 0.6  x  1.5 
2

(20)
S ( x )3  1 +1.4  x  1.5 
2

S ( x ) 4  1 + 2.2  x  1.5 
2

S ( x )5  1 + 5.1 x  1.5 
2

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.

5
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

3. Results and discussion

3.1. 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.

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 6. The velocity, density, and pressure in the nozzle with the shape equation of the nozzle being
S(x)=1+0.3(x-1.5)2

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 7. When the shape equation of the nozzle is S(x)=1+0.6(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure is 0.048

6
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 8. When the shape equation of the nozzle is S(x)=1+1.4(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure is 0.048

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 9. When the shape equation of the nozzle is S(x)=1+2.2(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure is 0.048

7
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 10. When the shape equation of the nozzle is S(x)=1+5.1(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure is 0.048

Here, we can use the continuity equation of ideal fluid to assist in qualitative analysis of flow field
changes:
S1v1  S 2 v2 (21)

where S1 , v1 are the cross-sectional area at one point of the nozzle and the local velocity, S 2 , v2 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.

3.1.1 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.

8
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(a) S(x)=1+0.3(x-1.5)2 (b) S(x)=1+0.6(x-1.5)2

(c) S(x)=1+1.4(x-1.5)2 (d) S(x)=1+2.2(x-1.5)2

(e) S(x)=1+5.1(x-1.5)2
Figure 11. Flow field rate comparison diagram in the nozzle of different shapes

3.1.2. Relationship between flow field pressure and nozzle shape According to the comparison of
pressure change images of different nozzles shown in 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.

(a) S(x)=1+0.3(x-1.5)2 (b) S(x)=1+0.6(x-1.5)2

9
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(c) S(x)=1+1.4(x-1.5)2 (d) S(x)=1+2.2(x-1.5)2

(e)S(x)=1+5.1(x-1.5)2
Figure 12. Pressure comparison diagram of the flow field in the nozzle of different shapes

We can also qualitatively explain the pressure change by analyzing the momentum equation. We
introduce a new physical quantity, called mass flow, which is defined as follows:
m  US (22)
Through the time advance solution, John D. Anderson [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 one-
dimensional 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:
 (  u S )  ( pS )
2
dS
+ = p (23)
x x dx
Since Equation (22) is constant, so equation (23) can be further simplified into the following
equation:
m U p
-  = (24)
S x x
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].

3.1.3 Relationship between flow field density and nozzle shape By comparing the curves of density in
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.

10
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(a) S(x)=1+0.3(x-1.5)2 (b) S(x)=1+0.6(x-1.5)2

(c) S(x)=1+1.4(x-1.5)2 (d) S(x)=1+2.2(x-1.5)2

(e) S(x)=1+5.1(x-1.5)2
Figure 13. Flow field density comparison diagram in the nozzle of different shapes

3.2. Full subsonic isentropic flow


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.

(a) Velocity (b) Density

11
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(c) Pressure (d) Shape of nozzle


Figure 14. When the shape equation of the nozzle is S(x)=1+0.3(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure 2.790

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 15. When the shape equation of the nozzle is S(x)=1+0.6(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure 2.790

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 16. When the shape equation of the nozzle is S(x)=1+1.4(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure 2.790

12
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 17. When the shape equation of the nozzle is S(x)=1+2.2(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure 2.790

(a) Velocity (b) Density

(c) Pressure (d) Shape of nozzle


Figure 18. When the shape equation of the nozzle is S(x)=1+5.1(x-1.5)2, the velocity, density, and
pressure of the flow field change. The nozzle boundary conditions are inlet pressure 3.000 and outlet
pressure 2.790

(a) S(x)=1+0.3(x-1.5)2 (b) S(x)=1+0.6(x-1.5)2

13
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

(c) S(x)=1+1.4(x-1.5)2 (d) S(x)=1+2.2(x-1.5)2

(e) S(x)=1+5.1(x-1.5)2
Figure 19. Comparison of full subsonic isentropic flow rate distribution when the nozzle inlet area is
arranged from small to large

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.
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.

3.3. 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.

3.3.2 Numerical Dissipation and Artificial Viscosity Numerically, the solution we get is not the exact
solution of the Euler equation, but the differential equation solution, and there is a discrete error
between the difference equation and the original solution of the Euler equation. So, we solve another
equation called the "correction equation". This correction equation is obtained by Taylor expansion of
the difference term of time and space in the discrete equation. So, this correction equation must have a

14
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

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].
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.

4. 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.

References
[1] [Tu J, Yeoh G H, Liu C. Computational fluid dynamics: a practical approach[M]. Butterworth-
Heinemann, 2018.[p1-p3]
[2] Wikipedia contributors. (2021, January 11). Nozzle. In Wikipedia, The Free Encyclopedia.
Retrieved 14:50, March 23, 2021, from
https://en.wikipedia.org/w/index.php?title=Nozzle&oldid=999696444
[3] Boyanapalli R, Vanukuri R SR, Gogineni P, et al. Analysis of composite De-Laval nozzle
suitable for rocket applications[J]. International Journal of Innovative Technology and
Exploring Engineering, 2013, 2: 336-344.
[4] "Nozzle design (converging/diverging - CD nozzle)". NASA.gov. Archived from the original on
20 March 2009. Retrieved 19 January 2009.
[5] Junmin Xu, Chongxi Huang.The influence of rocket design parameters on engine nozzle
efficiency.Rocket propulsion.1995,13(2):12
[6] Liu Minghao, Xu Xu, Fang Jie, Ma Chaokai. Optimization of nozzle profile of liquid rocket
engine based on CFD [A]. Beijing University of Aeronautics and Astronautics. Proceedings

15
ICMMAP 2021 IOP Publishing
Journal of Physics: Conference Series 2012 (2021) 012078 doi:10.1088/1742-6596/2012/1/012078

of the First National Symposium on Mechanics in the Aerospace Field (Volume 1) [C].
Beijing University of Aeronautics and Astronautics: Chinese Society of Mechanics, 2004:4.
[7] Cheng Yongsheng, Liu Hua, Xue Leiping. The NPD method is used to calculate the three-
dimensional nozzle airflow field[J]. Quarterly Journal of Mechanics, 2005 (04): 529-533.
[8] Wang Ping, Liu Xueshan, Qiao Limin. Flow field analysis of axisymmetric Laval nozzle[J].
Aircraft Design, 2013, 33 (02): 23-26.
[9] Anderson J D, Wendt J. Computational fluid dynamics[M]. New York: McGraw-Hill, 1995.

16

You might also like