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.
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.
The 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.
are each phase’s volume fraction, mass density, velocity, internal energy and total energy per unit mass, and the temperature, respectively; is the gas phase pressure, and here is also the mixture pressure.
The volume fraction meets the condition:
The gas phase adopts the ideal gas equation of state:
Let and be the specific heat coefficients of gas and solid particles; then .
are the interphase forces and the heat transfer in a unit volume, respectively.
The interphase forces: , where is the particle concentration, is the mass of a single particle, and is the drag resistance acting on a single particle [13].
Within 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: , where is convective heat transfer of a single particle [13]:where and are viscosity and thermal conductivity of the gas.
3. DG Method
If we setthen the above two-dimensional axisymmetric equations can be abbreviated as the following vector form.
Assuming that the solution area of the above equations is , the area is divided into nonoverlapped quadrilateral elements . Define a piecewise polynomial space:where represents the set of polynomials not exceeding degree on element .
Solution procedure of the above equations with DG method is as follows: find a numerical solution , such that for any experimental function and all , the following integral equation holds.
Applying Green formula to (11), we obtainwhere is the outward unit normal vector of edge ; is the numerical flux on the cell boundary.
The gas phase equations here adopt the numerical HLLC flux [24]:where the subscripts L and R indicate the left and right element of the edge. The intermediate state variables are as follows.
Wave velocity is defined aswhere and 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.
Here we chose first order polynomials as the basis functions for , and then an arbitrary component of vector can be expressed aswhereare the basis functions; is the center of element . , , and , and are the maximum and minimum coordinates of element in the X and Y directions, respectively.
Because formula (17) holds for all experimental functions , substituting it into (12) and making , formula (11) is equivalent towhere is the solution vector, is the mass matrix, and 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 is the throat radius and is the throat diameter, .
4.1. Computation Parameters
The combustion chamber stagnation condition is selected as , and ; the particle phase material density is , and its mass fraction is ; specific heat is . Two different-sized particles with radius and 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 Figure 2, 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 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.
The convergence history of the residuals in two-phase flow is presented in Figure 3. 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.
The numerical solution of two-phase flow on a fine mesh (380×100) was also computed. Figure 4 shows that there is only a slight difference between the flow fields with the coarse and fine meshes after convergence.
Figures 5 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 . 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 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.
Due to the presence of particles, the gas temperature also changes significantly from Figures 7 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.
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].
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]:where is the jet mass flow rate, is the jet velocity at the exit, is the outlet cross-sectional area, and 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 is , the engine thrust is 1635N, and thrust loss is 22%.(ii)When the particle radius is , 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 is , the total thrust is 2200N.(ii)When the particle radius is , 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 Figure 10, 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.
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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this manuscript.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grants Nos. 11571002, 11772067, and 11702028).