Discrete Dynamics in Nature and Society

Volume 2019, Article ID 7060481, 12 pages

https://doi.org/10.1155/2019/7060481

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

^{1}Graduate School of China Academy of Engineering Physics, Beijing 100088, China^{2}Institute of Applied Physics and Computational Mathematics, Beijing 100088, China

Correspondence should be addressed to Yu Xijun; nc.ca.mcpai@jxuy

Received 8 April 2019; Revised 31 May 2019; Accepted 13 June 2019; Published 14 July 2019

Academic Editor: Francisco R. Villatoro

Copyright © 2019 Duan Maochang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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