# Numerical Simulation of Gas-Particle Two-Phase Flow in a Nozzle with DG Method

Authors: Duan Maochang; Yu Xijun; Chen Dawei; Qing Fang; Zou Shijun
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

In this paper, the discontinuous Galerkin (DG) method is applied to solve the governing equations of the dispersed two-phase flow with the two-fluid Euler/Euler approach. The resulting governing equations are simple in form and the solution process is very natural. The characteristics of the gas-particle two-phase flow in an engine nozzle are mainly analyzed, and the impacts of the particle mass fraction and particle size on the flow field and engine performance are evaluated. Because of the addition of particles, the gas flow field undergoes significant modifications. Increase in the mass fraction leads to a significant thrust loss in the gas phase, and the impact of the particles on the gas phase could be substantial. Therefore, a quantitative study of thrust loss in the nozzle due to the particle impact is made. It is found that the gas thrust in the two-phase flow is reduced, but the total thrust of the two-phase flow increases to a certain extent.


---

## Body

## 1. Introduction

Gas-particle nozzle flows are typical examples of two-phase flows, widely used in aviation, aerospace, and other modern industry fields. Understanding the underlying physics of these flows is very important for normal operation and high performance of these devices.Metal-containing composite propellant, especially with aluminum particles, is widely used in the modern solid rocket engine in order to achieve higher specific impulse. Adding metal powder to the propellant enhances the energy of the composite solid propellant and inhibits unstable combustion. The gas-particle two-phase flow formed by the mixture of incomplete combustion particle and gas can cause thrust loss. The optional design of these engine nozzles is dependent on accurate knowledge of the flow behavior and is important to the attainment of high thrust efficiencies for launch vehicles [1]. Because the measurement of two-phase flow field is not an easy task, it is difficult to carry out experimental research on two-phase flow thrust loss. Early research mainly relied on empirical relations and simplified analytical models, which were too simple to consider, so the accuracy was relatively low. For example, in the early two-phase flow model, the interaction between phases was just one-way coupling, and this assumption is only valid for two-phase flow with mass fraction less than 10% [2]; Houglund [3] assumes that the two-phase flow in the nozzle is a quasi-one-dimensional flow, whereas the cross-section of the nozzle used varies dramatically at the throat.In recent years, numerical simulation has become an effective research tool and has achieved relatively abundant results. At present, there are two approaches for numerical simulation of gas-particle two-phase flow: two-fluid model and particle-trajectory model [2]. CJ Hwang and GC Chang [4] computed the gas-particle two-phase flow in a nozzle by using time-dependent explicit MacCormark scheme to solve the gas equations and tracking the particle trajectories with Lagrangian approach. Golafshani and Loh [5] also used particle-trajectory model to compute the mixed flow of viscous gas and particles in a JPL nozzle and estimated the deposition of particulate debris in the shuttle engine based on the computation results. Madabhushi et al. [6] studied the Aft-Dome two-phase flow field in solid rocket motor with submerged nozzle and a turbulence model was included into the computation. IS Chang [7] used two-fluid model and time marching method to solve the two-phase flow field in the inclined nozzle of a Titan engine. HT Chang et al. [8] solved the TLNS equations of two-fluid model by adopting a MUSCL-type flux splitting scheme and calculated the thrust loss caused by the two-phase flow in the nozzle. Many other researchers have used Eulerian-Lagrangian approach or Eulerian-Eulerian approach to compute complex subsonic-supersonic two-phase flow fields such as plumes in nozzles [9–11].So far, the numerical methods for solving two-phase flow governing equations are mainly finite difference methods, which are relatively mature and rich in results [12, 13]. In 1973, Reed and Hill [14] first introduced the discontinuous Galerkin method to solve the neutron transport equation. Until the 1990s, Cockburn and Shu [15–17] combined it with TVD time discretization to solve the hyperbolic conservation laws. Then the RKDG method slowly attracted researchers' attention and began to be applied to the computational fluid dynamics. The DG method not only has the same characteristics as the finite element method to construct high-order schemes easily, but also has high-resolution shock-capturing capability and can well adapt to increasingly popular unstructured grids; its adaptive and parallel processing is also flexible and convenient. As a result, it is more adaptable to the future trend of large-scale scientific computing. DG method is also natural to be applied as LES on these complex flows. This will help the industry deviate from RANS to LES or ILES hopefully.In recent years, many researchers worked on DG method to solve multiphase flow. Franquet and Perrier [18] successfully solved the governing equations of Baer-Nunziato model, and it was also extended to interface capturing and front tracking algorithm of multiphase flow [19, 20]. The latest progress is made by O. Ejtehadi et al. [21], who introduced a novel treatment of two-phase flow source terms, making the results very robust and capable of calculating complex wave systems in two-phase flows.In this paper, Runge–Kutta DG method is first applied to solve the governing equations of two-dimensional axisymmetric dispersed two-phase flow in a JPL nozzle. No operator splitting is introduced to handle the source terms, and time advancing is not limited by the small time step. In particular, the influence of particle phase on the gas flow in the nozzle and the effect on performance of the nozzle are investigated.


## 2. Physical Model and Governing Equations

The two-fluid model is based on continuum assumption under Eulerian-Eulerian framework, and the discrete particles are averaged as a new fluid that fills the same space with the real fluid at the same time [22]. The interaction between the two phases in the model is two-way coupled. The obtained two-phase flow governing equations are unified in form and similar to the single-phase flow equations, so many numerical methods for single-phase flow computational fluid dynamics can be used for reference.The interaction between two phases is very complex, including many forms of force [23]. In order to make the model as simple as possible, suppose the two-phase flow we want to study has the following characteristics:(1)
Gas phase is ideally compressible.(2)
Viscosity and heat conduction only exist in the interphase interaction.(3)
Only the drag resistance of particles is considered.(4)
All particles are ideal spherical and solid and have a uniform diameter.(5)
There is no collision between particles.(6)
Particle phase is incompressible in materiel density and the pressure of the mixture is equal to the gas pressure.(7)
Interphase interaction is the same as that of a single particle.The two-dimensional axisymmetric governing equations approximately describe the gas-particle two-phase flow with the above characteristics in the nozzle, and its form is as follows.(1)
∂
∂
t
α
g
ρ
g
α
g
ρ
g
u
g
α
g
ρ
g
v
g
α
g
ρ
g
E
g
α
p
ρ
p
α
p
ρ
p
u
p
α
p
ρ
p
v
p
α
p
ρ
p
E
p
+
∂
∂
x
α
g
ρ
g
u
g
α
g
ρ
g
u
g
2
+
P
α
g
ρ
g
u
g
v
g
α
g
ρ
g
E
g
+
P
u
g
α
p
ρ
p
u
p
α
p
ρ
p
u
p
2
α
p
ρ
p
u
p
v
p
α
p
ρ
p
u
p
E
p
+
∂
∂
y
α
g
ρ
g
v
g
α
g
ρ
g
u
g
v
g
α
g
ρ
g
v
g
2
+
P
α
g
ρ
g
E
g
+
P
v
g
α
p
ρ
p
v
p
α
p
ρ
p
u
p
v
p
α
p
ρ
p
v
p
2
α
p
ρ
p
v
p
E
p
+
1
y
α
g
ρ
g
v
g
α
g
ρ
g
u
g
v
g
α
g
ρ
g
v
g
2
α
g
ρ
g
E
g
+
P
v
g
α
p
ρ
p
v
p
α
p
ρ
p
u
p
v
p
α
p
ρ
p
v
p
2
α
p
ρ
p
v
p
E
p
=
0
F
g
x
F
g
y
Q
g
+
V
p
·
F
g
0
F
p
x
F
p
y
Q
p
+
V
p
·
F
pThe first four equations are the gas phase governing equations, the last four are the particle phase governing equations, the fourth left term of the equations is the geometric source term, and the right term is source term taking into account the phase interaction; the subscript g represents the gas phase, and the subscript p represents the particle phase.α g , ρ g , u g , v g , e g , E g , T g , α p , ρ p , u p , v p , e p , E p , T p are each phase’s volume fraction, mass density, velocity, internal energy and total energy per unit mass, and the temperature, respectively; P is the gas phase pressure, and here is also the mixture pressure.The volume fraction meets the condition:(2)αg+αp=1The gas phase adopts the ideal gas equation of state:(3)P=γ-1ρgegLet(Cpg,Cvg) and Cp be the specific heat coefficients of gas and solid particles; then Eg=eg+(1/2)ug2=CvgTg+(1/2)ug2,Ep=ep+(1/2)up2=CpTp+(1/2)up2.F g = ( F g x , F g y ) , F p = ( F p x , F p y ) , Q g , Q p are the interphase forces and the heat transfer in a unit volume, respectively.The interphase forces:Fp=-Fg=(σ/m)D, where σ=αpρp is the particle concentration, m=(4/3)πrp3ρp is the mass of a single particle, and D is the drag resistance acting on a single particle [13].(4)D=12πrp2ρgVg-VpVg-VpCD=12πrp2ρgVg-VpVg-Vp0.48+28Re-0.85(5)Re=ρgdpug-up2+vg-vp2μgWithin the range of Reynolds number 0-1000, the drag coefficient in (4) gives a good fit to the standard drag coefficient curve. If the Reynolds number exceeds this range, different forces such as the history force would come into question.Interphase convective heat transfer:Qp=-Qg=(σ/m)Q, where Q is convective heat transfer of a single particle [13]:(6)Q=πdpμCpgPr-1Tg-TpNu=πdpμCpgPr-1Tg-Tp2.0+0.6Pr1/3Re1/2(7)Pr=Cpgμgkgwhere μg and kg are viscosity and thermal conductivity of the gas.


## 3. DG Method

If we set(8)U=αgρgαgρgugαgρgvgαgρgEgαpρpαpρpupαpρpvpαpρpEp,FU=αgρgugαgρgug2+PαgρgugvgαgρgEg+Pugαpρpupαpρpup2αpρpupvpαpρpupEp,GU=αgρgvgαgρgugvgαgρgvg2+PαgρgEg+Pvgαpρpvpαpρpupvpαpρpvp2αpρpvpEpRU=1yαgρgvgαgρgugvgαgρgvg2αgρgEg+Pvgαpρpvpαpρpupvpαpρpvp2αpρpvpEp,SU=0FgxFgyQg+Vp·Fg0FpxFpyQp+Vp·Fp,then the above two-dimensional axisymmetric equations can be abbreviated as the following vector form.(9)∂U∂t+∂FU∂x+∂GU∂y+RU=SUAssuming that the solution area of the above equations isΩ, the area Ω is divided into nonoverlapped quadrilateral elements Ω=⋃Ωij∣i=1,N,j=1,M. Define a piecewise polynomial space:(10)Vhk=vh∈L2Ω:vhΩij∈PkΩij;1≤i≤N,1≤j≤Mwhere Pk(Ωij) represents the set of polynomials not exceeding degree k on element Ωij.Solution procedure of the above equations with DG method is as follows: find a numerical solutionUh∈Vhk, such that for any experimental function vh∈Vhk and all 1≤i≤N,1≤j≤M, the following integral equation holds.(11)∫Ωijvh∂Uh∂tdΩ+∫Ωijvh∇·FUh,GUhdΩ+∫ΩijvhRUhdΩ=∫ΩijvhSUhdΩApplying Green formula to (11), we obtain(12)∫Ωijvh∂Uh∂tdΩ=-∫∂ΩijF^Uh,G^Uh·nvhdΩ+∫ΩijFUh,GUh·∇vhdΩ+∫ΩijvhSUh-RUhdΩwhere n is the outward unit normal vector of edge ∂Ωij; F^(Uh),G^(Uh)·n=H^Uh-,Uh+, is the numerical flux on the cell boundary.The gas phase equations here adopt the numerical HLLC flux [24]:(13)H^=F^L,ifSL>0F^UL∗,ifSL≤0<SMF^UR∗,ifSM≤0<SRF^R,ifSR≤0where the subscripts L and R indicate the left and right element of the edge. The intermediate state variables are as follows.(14a)UL∗=αρL∗αρVL∗αρEL∗=1SL-SMSL-ugnLαρLSL-ugnLαρVL+P∗-PL·nSL-ugnLαρEL-PLugnL+P∗SM(14b)UR∗=αρR∗αρVR∗αρER∗=1SR-SMSR-ugnRαρRSR-ugnRαρVR+P∗-PR·nSR-ugnRαρER-PRugnR+P∗SM(14c)F^UL∗=SMαρL∗SMαρVL∗+P∗·nSMαρEL∗+P∗SM(14d)F^UR∗=SMαρR∗SMαρVR∗+P∗·nSMαρER∗+P∗SM(14e)P∗=ρLugnL-SLugnL-SM+PL=ρRugnR-SRugnR-SM+PR(14f)SM=ρRugnRSR-ugnR-ρLugnLSL-ugnL+PL-PRρRSR-ugnR-ρLSL-ugnLWave velocity is defined as(15a)SL=min⁡ugnL-cL,V·n-c(15b)SR=minugnR-cR,V·n+cwhere V and c are the Roe averages of fluid velocities and sound velocities in the left and right element, respectively.Since the particle phase equations are weakly hyperbolic [25], many classical algorithms are not suitable. The Lax-Friedrich numerical flux is used for the particle phase equations.(16)H^=12FUL,GUL·n+FUR,GUR·n-λUR-ULHere we chose first order polynomials as the basis functions forVhk, and then an arbitrary component uh of vector Uh can be expressed as(17)uhx,t=∑l=1,3uijltBlijx,y,x,y∈Ωijwhere(18)B1ij=1,B2ij=x-xcΔx,B3ij=y-ycΔyare the basis functions; (xc,yc) is the center of element Ωij. Δx=0.5(xmax-xmin), Δy=0.5(ymax-ymin), and xmax,xmin,ymax, and ymin are the maximum and minimum coordinates of element Ωij in the X and Y directions, respectively.Because formula (17) holds for all experimental functions vh∈Vhk, substituting it into (12) and making vh=Bmij,m=1,2,3, formula (11) is equivalent to(19)ddtMij·uhij=LUhwhere uhij=(uij1,uij2,uij3) is the solution vector, Mij=∫ΩijBlijBmijdΩij(l=1,2,3;m=1,2,3) is the mass matrix, and L is the space discrete operator on the right.In this paper, the second-order accurate TVD Runge–Kutta time discretization is adopted, the integration in the element is approximated by a five-point Gaussian integration. In order to guarantee the stability of the scheme and suppress the potential spurious oscillations in the numerical solution, the minmod limiter with TVB correction [26] is used to limit the slope of the numerical solution.


## 4. Results Analysis

The JPL nozzle is a 45°-15° convergent-divergent conical nozzle [27]. The ratio of throat radius of curvature to throat radius is 0.625. The dimensions are shown in Figure 1, where R is the throat radius and D is the throat diameter, D=0.0406m.Figure 1
The nozzle configuration and grid system.


### 4.1. Computation Parameters

The combustion chamber stagnation condition is selected asT0=555K, and P0=10.34×105N/m2; the particle phase material density is ρp=4004.62kg/m3, and its mass fraction is ϕ=30%; specific heat is Cpp=1.38×103J/(kg·K). Two different-sized particles with radius 1.0μm and 10.0μm are considered.The subsonic inlet condition is used for the gas phase at the entrance, and the exit boundary condition is simple linear extrapolation from interior points since the gas is assumed to be supersonic at the exit plane. The nozzle wall is impermeable and adiabatic, so slip boundary condition is used. On the centreline, the symmetric boundary conditions are required.The particle phase has the same velocity and temperature as those of the gas phase at the inlet, with a mass fraction of 30%. Other boundary conditions are similar to the gas phase.The computation was carried out on a 190×50 structured grid with a single processor considering the low grid number. The time step is determined by CFL=0.2.


### 4.2. Comparison of Single-Phase to Two-Phase Flow

Since the JPL nozzle experiment only measured the single-phase flow, there is no experimental data for the two-phase flow [28]. Therefore, the single-phase flow in the nozzle is first simulated, to validate the numerical modeling.As shown in Figure2, the gas is accelerating along the centerline and the wall of the nozzle, the acceleration near the wall is faster, and the speed of sound is reached at the throat first. Due to the small throat radius of curvature, a weak shock wave is formed near the wall as the gas flows through the throat to divergent section, propagates downstream of the centerline, and reflects approximately at X=2 on the centerline. As the gas passes through the shock wave here, the Mach number drops significantly, then the gas accelerates again, and the Mach number continues to increase. The numerical results are basically consistent with the experimental data [28], except that there is a slight difference in resolution near the shock wave.Figure 2
One-phase flow field and Mach number profiles.The convergence history of the residuals in two-phase flow is presented in Figure3. It can be seen that after about 6000 time steps, the continuity equation residuals of gas and particle phase reach 10−6 and 10−4, respectively, and the serial execution of the code takes less than 5 minutes.Figure 3
Convergence history in two-phase flow.The numerical solution of two-phase flow on a fine mesh (380×100) was also computed. Figure4 shows that there is only a slight difference between the flow fields with the coarse and fine meshes after convergence.Figure 4
Comparison of Mach number contour lines between two grid systems.Figures5 and 6 show that some characteristics of the single-phase flow field are also reflected in the two-phase flow field. The gas still achieves sound velocity at the throat, but because of the addition of particle phase, the strength of the downstream shock wave is significantly weakened. The shock wave is flattened in the case of particle radius being 1.0μm. The Mach number of two-phase flow field is less than that of pure gas at the same position. This is because the solid particles cannot be accelerated by a convergent-divergent process, and their inertia is large, so their movement state cannot be changed rapidly. Their velocity lags behind the gas hinders acceleration and expansion of the gas. In the two-phase flow field containing 10.0μm particles, the shock wave downstream of the throat is dissipated, and the gas expansion is somewhat low, compared with HT Chang et al.’s results [8]. This may be related to the numerical viscosity of TVB modified limiter.Figure 5
Mach number distribution of two-phase flow field.Figure 6
Mach number profiles of two-phase flow.Due to the presence of particles, the gas temperature also changes significantly from Figures7 and 8; it can be seen that the gas rapidly expands after passing through the throat and the temperature decreases substantially, while the particles maintain their initial temperature due to the thermal inertia, and thus the gas phase is heated downstream. In particular, on the centerline, the gas temperature in two-phase flow field is approximately 200K higher than that in single-phase flow. Near the wall of the divergent section, the particle concentration is very low and has little influence on the gas temperature.Figure 7
Temperature profiles along centerline.Figure 8
Gas temperature profiles along wall.After the above analysis, it can be said that, under the same mass fraction, the smaller the particle size, the more the number of particles contained in a unit mass of propellant, the larger the contact area between the particles and the gas, and therefore the stronger the interaction between the particle phase and the gas phase. The stronger the interaction, the greater the friction and heat exchange between the two phases, and the greater the influence of particles on the gas flow.


### 4.3. Particle Distribution

After passing through the narrow throat, the gas can quickly change direction and accelerate along the wall of the divergent section. However, due to the inertia, particles cannot change their direction and velocity easily. Therefore, the particles near the wall of the convergent section detach from the wall after passing through the throat, slowly change their directions during the gas-induced acceleration, and form a limiting particle streamline. Outside this streamline, there is a region where the particle trajectory cannot reach, namely, a particle-free zone.Smaller particles have good tracking ability, which can be easily affected by gas movement, and are more likely to fill the entire nozzle with smaller particle-free zone. The larger particles have larger inertia, and their initial momentum cannot be modified without difficulty. There are large particle-free zones, and the particles are relatively concentrated near the centerline. The particle volume fraction at the entrance is about 7.0×10−4 and is lower than 10-3 in most areas of the flow field. The results in Figure 9 are consistent with those obtained using the Lagrangian method [4].Figure 9
Particle trajectory in the nozzle.


### 4.4. Propulsion Performance

The above study clearly shows that due to the addition of particles, the outlet jet velocity is greatly reduced. The effect of this change on engine performance is undoubtedly enormous, and the extent of this impact must be quantified.The engine thrust formula is as follows [29]:(20)F=m˙Ve+AePe=ρeVe2+PeAewhere m˙ is the jet mass flow rate, Ve is the jet velocity at the exit, Ae is the outlet cross-sectional area, and Pe is the outlet pressure.In this paper, the engine thrust obtained by the numerical computation of pure gas flow is 2097N. In order to examine the effect of particles on the gas flow, the thrust generated by the gas mass flow in the two-phase flow is first calculated and the results are as follows:(i)
When the particle radius is1.0μm, the engine thrust is 1635N, and thrust loss is 22%.(ii)
When the particle radius is10.0μm, the engine thrust is 1818N, and thrust loss is 13%.HT Chang et al. [8] also got similar results in their numerical study. It can be found that the interaction between the two phases caused by the small-sized particles is stronger, leading to a thrust loss of up to 22% in the gas phase, but a loss of 13% for the large-sized particles.For the space launch mission, every Newton of thrust is valuable. If the addition of particles would cause such a huge thrust loss, then even if there are many other benefits, adding metal powder to the propellant is not worth the cost. In fact, when the particle phase is accelerated by the gas phase, it also gains momentum. This part of momentum also generates additional thrust.After considering this additional particulate mass flow, the total thrust generated by the two-phase flow mixture is as follows.(i)
When the particle radius is1.0μm, the total thrust is 2200N.(ii)
When the particle radius is10.0μm, the total thrust is 2168N.The above computation results show that the gas outlet velocity is reduced significantly and the gas phase thrust is greatly lost in the case of two-phase flow; however, due to the existence of the particles, the outlet mass flow is increased, and additional nozzle thrust is generated. The total thrust of the engine does not drop; it increases. This is part of the reason why metal powder is added to the propellant even if there is a serious two-phase loss. It can also be seen that although the influence of different-sized particles on the gas phase varies greatly, the total thrust is almost the same under the same particle mass fraction.As shown in Figure10, the thrust loss of the gas phase in two-phase flow and the total thrust considering the additional thrust of the particle phase are finally calculated under different particle mass fractions and different particle sizes. The dotted line represents the thrust of pure gas flow, the thrust of gas phase in two-phase flow is below the dotted line, and the y-coordinate on the right shows the scale of the loss. The variation of thrust has a strong regularity: the larger the particle mass fraction is and the smaller the particle size is, the more the thrust loss of the gas phase is, and the total thrust variation does not have a clear trend. It should be noted that the particle mass fraction is not always the larger the better, and the particle size is not the smaller the better. This point requires more work to be clarified.Figure 10
Effect of particle mass fraction and size on thrust.


## 4.1. Computation Parameters

The combustion chamber stagnation condition is selected asT0=555K, and P0=10.34×105N/m2; the particle phase material density is ρp=4004.62kg/m3, and its mass fraction is ϕ=30%; specific heat is Cpp=1.38×103J/(kg·K). Two different-sized particles with radius 1.0μm and 10.0μm are considered.The subsonic inlet condition is used for the gas phase at the entrance, and the exit boundary condition is simple linear extrapolation from interior points since the gas is assumed to be supersonic at the exit plane. The nozzle wall is impermeable and adiabatic, so slip boundary condition is used. On the centreline, the symmetric boundary conditions are required.The particle phase has the same velocity and temperature as those of the gas phase at the inlet, with a mass fraction of 30%. Other boundary conditions are similar to the gas phase.The computation was carried out on a 190×50 structured grid with a single processor considering the low grid number. The time step is determined by CFL=0.2.


## 4.2. Comparison of Single-Phase to Two-Phase Flow

Since the JPL nozzle experiment only measured the single-phase flow, there is no experimental data for the two-phase flow [28]. Therefore, the single-phase flow in the nozzle is first simulated, to validate the numerical modeling.As shown in Figure2, the gas is accelerating along the centerline and the wall of the nozzle, the acceleration near the wall is faster, and the speed of sound is reached at the throat first. Due to the small throat radius of curvature, a weak shock wave is formed near the wall as the gas flows through the throat to divergent section, propagates downstream of the centerline, and reflects approximately at X=2 on the centerline. As the gas passes through the shock wave here, the Mach number drops significantly, then the gas accelerates again, and the Mach number continues to increase. The numerical results are basically consistent with the experimental data [28], except that there is a slight difference in resolution near the shock wave.Figure 2
One-phase flow field and Mach number profiles.The convergence history of the residuals in two-phase flow is presented in Figure3. It can be seen that after about 6000 time steps, the continuity equation residuals of gas and particle phase reach 10−6 and 10−4, respectively, and the serial execution of the code takes less than 5 minutes.Figure 3
Convergence history in two-phase flow.The numerical solution of two-phase flow on a fine mesh (380×100) was also computed. Figure4 shows that there is only a slight difference between the flow fields with the coarse and fine meshes after convergence.Figure 4
Comparison of Mach number contour lines between two grid systems.Figures5 and 6 show that some characteristics of the single-phase flow field are also reflected in the two-phase flow field. The gas still achieves sound velocity at the throat, but because of the addition of particle phase, the strength of the downstream shock wave is significantly weakened. The shock wave is flattened in the case of particle radius being 1.0μm. The Mach number of two-phase flow field is less than that of pure gas at the same position. This is because the solid particles cannot be accelerated by a convergent-divergent process, and their inertia is large, so their movement state cannot be changed rapidly. Their velocity lags behind the gas hinders acceleration and expansion of the gas. In the two-phase flow field containing 10.0μm particles, the shock wave downstream of the throat is dissipated, and the gas expansion is somewhat low, compared with HT Chang et al.’s results [8]. This may be related to the numerical viscosity of TVB modified limiter.Figure 5
Mach number distribution of two-phase flow field.Figure 6
Mach number profiles of two-phase flow.Due to the presence of particles, the gas temperature also changes significantly from Figures7 and 8; it can be seen that the gas rapidly expands after passing through the throat and the temperature decreases substantially, while the particles maintain their initial temperature due to the thermal inertia, and thus the gas phase is heated downstream. In particular, on the centerline, the gas temperature in two-phase flow field is approximately 200K higher than that in single-phase flow. Near the wall of the divergent section, the particle concentration is very low and has little influence on the gas temperature.Figure 7
Temperature profiles along centerline.Figure 8
Gas temperature profiles along wall.After the above analysis, it can be said that, under the same mass fraction, the smaller the particle size, the more the number of particles contained in a unit mass of propellant, the larger the contact area between the particles and the gas, and therefore the stronger the interaction between the particle phase and the gas phase. The stronger the interaction, the greater the friction and heat exchange between the two phases, and the greater the influence of particles on the gas flow.


## 4.3. Particle Distribution

After passing through the narrow throat, the gas can quickly change direction and accelerate along the wall of the divergent section. However, due to the inertia, particles cannot change their direction and velocity easily. Therefore, the particles near the wall of the convergent section detach from the wall after passing through the throat, slowly change their directions during the gas-induced acceleration, and form a limiting particle streamline. Outside this streamline, there is a region where the particle trajectory cannot reach, namely, a particle-free zone.Smaller particles have good tracking ability, which can be easily affected by gas movement, and are more likely to fill the entire nozzle with smaller particle-free zone. The larger particles have larger inertia, and their initial momentum cannot be modified without difficulty. There are large particle-free zones, and the particles are relatively concentrated near the centerline. The particle volume fraction at the entrance is about 7.0×10−4 and is lower than 10-3 in most areas of the flow field. The results in Figure 9 are consistent with those obtained using the Lagrangian method [4].Figure 9
Particle trajectory in the nozzle.


## 4.4. Propulsion Performance

The above study clearly shows that due to the addition of particles, the outlet jet velocity is greatly reduced. The effect of this change on engine performance is undoubtedly enormous, and the extent of this impact must be quantified.The engine thrust formula is as follows [29]:(20)F=m˙Ve+AePe=ρeVe2+PeAewhere m˙ is the jet mass flow rate, Ve is the jet velocity at the exit, Ae is the outlet cross-sectional area, and Pe is the outlet pressure.In this paper, the engine thrust obtained by the numerical computation of pure gas flow is 2097N. In order to examine the effect of particles on the gas flow, the thrust generated by the gas mass flow in the two-phase flow is first calculated and the results are as follows:(i)
When the particle radius is1.0μm, the engine thrust is 1635N, and thrust loss is 22%.(ii)
When the particle radius is10.0μm, the engine thrust is 1818N, and thrust loss is 13%.HT Chang et al. [8] also got similar results in their numerical study. It can be found that the interaction between the two phases caused by the small-sized particles is stronger, leading to a thrust loss of up to 22% in the gas phase, but a loss of 13% for the large-sized particles.For the space launch mission, every Newton of thrust is valuable. If the addition of particles would cause such a huge thrust loss, then even if there are many other benefits, adding metal powder to the propellant is not worth the cost. In fact, when the particle phase is accelerated by the gas phase, it also gains momentum. This part of momentum also generates additional thrust.After considering this additional particulate mass flow, the total thrust generated by the two-phase flow mixture is as follows.(i)
When the particle radius is1.0μm, the total thrust is 2200N.(ii)
When the particle radius is10.0μm, the total thrust is 2168N.The above computation results show that the gas outlet velocity is reduced significantly and the gas phase thrust is greatly lost in the case of two-phase flow; however, due to the existence of the particles, the outlet mass flow is increased, and additional nozzle thrust is generated. The total thrust of the engine does not drop; it increases. This is part of the reason why metal powder is added to the propellant even if there is a serious two-phase loss. It can also be seen that although the influence of different-sized particles on the gas phase varies greatly, the total thrust is almost the same under the same particle mass fraction.As shown in Figure10, the thrust loss of the gas phase in two-phase flow and the total thrust considering the additional thrust of the particle phase are finally calculated under different particle mass fractions and different particle sizes. The dotted line represents the thrust of pure gas flow, the thrust of gas phase in two-phase flow is below the dotted line, and the y-coordinate on the right shows the scale of the loss. The variation of thrust has a strong regularity: the larger the particle mass fraction is and the smaller the particle size is, the more the thrust loss of the gas phase is, and the total thrust variation does not have a clear trend. It should be noted that the particle mass fraction is not always the larger the better, and the particle size is not the smaller the better. This point requires more work to be clarified.Figure 10
Effect of particle mass fraction and size on thrust.


## 5. Conclusion

In this paper, the two-dimensional axisymmetric governing equations of gas-particle two-phase flow are solved by the DG method with two-fluid approach. The interaction between the two phases in the model is two-way coupled, and the numerical flux is designed according to the properties of the gas phase and particle phase equations. Gas-particle two-phase flow field in a JPL nozzle is studied. The differences of flow field Mach number and temperature distribution between single-phase flow and two-phase flow are analyzed; the influence of different-sized particles on the flow field and the distribution and trajectory of different-sized particles in the nozzle are also investigated. At last, the two-phase flow engine thrust loss caused by the presence of solid particles is studied, and the influence of different particle mass fraction and particle sizes on the rocket engine propulsion performance is evaluated. Results show that the presence of the particle phase will hinder the acceleration and expansion of gas phase. The smaller the particles are, the stronger the interaction is, and the greater the gas thrust loss will be; specific impulse of the engine will be much decreased. However the particles also generate an additional thrust, which has a certain effect on increasing the total engine thrust.This study makes a preliminary attempt to apply DG method to the numerical simulation of gas-particle two-phase flow, which lays a foundation for further research on more complex two-phase flows. In the future, we would like to choose a less dissipative limiter to further improve the resolution of the algorithm, so that we can capture finer flow structures. In addition, we are also interested in unsteady two-phase flow outside the nozzle. This is the two-phase flow problem in real flight environment. The influence of particles on the shape of the external jet and the shock structure is very significant. We hope that the DG method can play a leading role in the solution of these problems.



---

*Source: 7060481-2019-07-14.xml*