Theoretical and numerical analysis on choked multiphase flows of gas and solid particle through a convergent–divergent nozzle

Abstract

Particle–gas two-phase flows show significantly different behaviours compared to single-gas flow through a convergent–divergent nozzle. Non-equilibrium effects, thermal and velocity lags lead to the inefficiency of nozzle performance. In the present studies, theoretical analysis and numerical simulations were carried out to study particle–gas flows through a convergent–divergent nozzle. Homogeneous and equilibrium model that no slip in velocity and temperature occurs between particle phase and gas phase was considered to derive mass flow rate and sound speed of multiphase flows. Two-phase flows are regarded as isentropic flows that isentropic relations can be used for homogeneous equilibrium model. Particle volume fraction, specific heat ratio and gas constant parameter were also rewritten. Discrete phase model that includes Lagrangian-Eulerian tracking method was used to calculate particle–gas flows in computational fluid dynamics (CFD) studies. Particle mass loading was varied to investigate its effect on choking phenomena for particle–gas flows at different nozzle pressure ratios. Mass flow rate and sound speed of two-phase flows were theoretically calculated and compared with numerical results. Particle number density and particle velocity vectors were also obtained to explain particle motion through convergent–divergent nozzle.

Introduction

When aircrafts or rockets start taking off, the residuals from reaction of fuel mixtures developing in combustion chamber will be moved outside as the exhaust gas. Residuals are always liquid drops or solid particles, which lead to the erosion damage sustained by walls of rocket nozzle. Particle–gas flows through supersonic nozzle show significantly different behaviours from single-gas flow. Solid particles in the exhaust result to inefficiency in expansion process in propulsive nozzle due to the non-equilibrium effects and thermal lags between gas phase and solid particles as well as velocity slip. Interaction between gas phase and particle phase, particle collision and thermal transfer result to the loss in momentum and energy of multiphase flows.

Needle drug delivery devices have been widely used for drug injection in medical fields. Solid drug powders can be directly injected to human bodies without any sharp metal needles and skin injury. Supersonic nozzles are used to accelerate momentum of drug powders. However, researchers are facing the same problems that momentum and mass of drug powders are significantly difficult to be controlled due to non-equilibrium powder–gas flows. Therefore, investigations of particle–gas two-phase flows induced by supersonic nozzles have significantly practical meaning on improvement in performance of aircrafts and needle-free drug delivery devices in aerospace and medical engineering fields.

Numerical studies were performed on multiphase flows through supersonic micro nozzles at different diameters of liquid drop. Liquid drops with smaller diameter and Stokes number were observed to track supersonic gas flows well, and mass loading of liquid drops was obtained to decrease the nozzle preference in thrust force. ¹ A theoretical model was derived to predict the maximum mass flow rate of single-gas flow and equilibrium two-phase mixture flows through a convergent nozzle. ² The maximum mass flow rate was obtained to depend on local slip ratio and stagnation pressure at initial state. Fauske theoretically investigated mass flow rate and heat transfer rate of multiphase flows through a convergent nozzle and compared obtained results with previous experimental ones. The new theoretical model contained the most available data to predict reliable heat transfer process instead of using frozen or complete equilibrium mass and heat transfer process. Critical mass flow rate was calculated to agree experimental results well.^3,4

Richer derived a separated model to calculate critical mass flow rate of steam–water mixture flows. Hydrodynamic and non-equilibrium thermal effects caused by rapid depressurization were investigated. Interaction terms of interphase for mass, momentum and energy were also considered. The relative velocity between two phases was observed to be affected by interfacial forces. ⁵ Theoretical and numerical studies on particle–gas flows were carried out by Ling et al. Frozen and equilibrium models were used for theoretical investigation and compared with results from numerical simulations. Different time scales with which the solution is closed to the equilibrium state and particle mass loading (PML) were investigated. The interaction between particles and expansion fan showed that the longest equilibrium scale time was associated with expansion fan. ⁶ Different models for calculating multiphase flows in mass, momentum and heat transfer were summarized by Dauria and Vigni. Equilibrium, non-equilibrium and frozen theories on calculation of particle–gas flows were discussed in detail. ⁷

Akmandor and Nagashima investigated vapor–liquid flows through a convergent–divergent nozzle using equilibrium and non-equilibrium models. Sound speed of mixture was homogenously calculated. Mass flow rate and sound speed of mixture were obtained under choked and unchoked conditions at throat, respectively. Theoretical calculations agreed with experimental results well. ⁸ Choked mass flow rate of particle–gas flows at nozzle throat was also analytically calculated using homogeneous and equilibrium model by Wang et al. Particle–gas flows were regarded as isentropic flows which all isentropic equations can be used. Experimental studies were also done to compare with theoretical results. ⁹

Rudinger carried out theoretical analysis on sound speed and mass flow rate of particle–gas flows. PML and particle volume fraction were investigated to have obvious influence on sound speed of mixture and choked mass flow rate through a convergent–divergent nozzle. Particle velocity and thermodynamics were also considered.^10,11 Hwang and Chang performed numerical studies on gas–particle flows through a rocket nozzle. Particle Mach number and temperature were obtained, and different particle diameters were investigated. ¹² Chang et al. applied flux vector-splitting scheme to calculate single-gas flow and dilute particle–gas flows through a jet propulsion laboratory (JPL) nozzle. The accuracy of the scheme was tested to be reliable for simulating dilute two-phase flows. ¹³

Neilson and Gilchrist carried out analytical and experimental studies on effects of the ratio of particle mass flow rate to total mass flow rate on gas flows through ducts with constant cross-sectional area and convergent–divergent nozzle. For subsonic particle–gas flows, the ratio of particle mass flow rate to total mass flow rate had no influence on gas velocity if it was less than 0.8. For supersonic particle–gas flows through a convergent–divergent nozzle, gas velocity was observed to be not affected by particle loading if this ratio was less than 0.3. ¹⁴ Quasi-one-dimensional model was proposed to calculate supersonic particle–gas flows through JPL nozzle by Forde. Shock wave location and structure were observed to be influenced by particle loading ratios and particle diameter. In addition, shock wave strength was strongly determined by particle loading ratios. ¹⁵ Chang numerically simulated single-gas flow and particle–gas flows through a JPL nozzle at different particle diameter and particle mass fractions. Larger-sized particles and higher particle loading ratios decreased gas velocity faster. The assumption of constant lag in velocity between gas phase and particle phase was shown to be not reasonable for simulating transonic particle–gas flows. ¹⁶ Particle behaviours were shown to be greatly different in supersonic and subsonic flows based on Henderson's studies. He proposed different drag force models for calculating subsonic and supersonic particle–gas flows, respectively. ¹⁷

In this article, theoretical and numerical studies were carried out to investigate choked phenomena and sound speed of particle–gas flows. Mass flow rate and sound speed of particle–gas flows were theoretically calculated by homogeneous and equilibrium model and compared with numerical results at different nozzle pressure ratios (NPRs) and PMLs. Particle number density (PND) and particle concentration were also considered.

Theoretical analysis

Mass flow rate

There is no slip in velocity and temperature between particle phase and gas phase, and particle–gas flows can be regarded as isentropic flows, which are described by equilibrium model.^8,9 Particle–gas flows are choked at nozzle throat as operating pressure ratio is high enough. Mass flow rate is an important parameter to indicate choked characteristics of particle–gas flows. PML x is the main factor that affects mass flow rate of particle–gas flows. It is defined as

$$x=\frac{{m^{•}}_p}{{m^{•}}_p+{m^{•}}_g}$$

(1)

As PML increases, particle volume fraction ɛ also increases. Based on different particle volume fractions, particle–gas flows can be described as dilute and dense flows. Particle volume fraction ɛ that is related to PML, particle and gas density is calculated by equation (2).

Due to large density of solid particle, particle volume fraction is calculated to be low at high PML. At PML of 5%, particle volume fraction is obtained to be 0.05%, and gas–particle flows are regarded as dilute flows. At PML of 51%, particle volume fraction is calculated to be approximate to 1%.

$$\varepsilon =\frac{\frac{x}{1-x}\frac{{\rho }_g}{{\rho }_p}}{1+\frac{x}{1-x}\frac{{\rho }_g}{{\rho }_p}}$$

(2)

Specific heat ratio of mixture γ _e and a new two-phase constant parameter R _e are different from properties of single-gas phase, which are calculated based on equations (3) and (4)

$${\gamma }_e={\gamma }_g\frac{1+\frac{x}{1-x}\frac{c_{\mathit{pp}}}{c_{\mathit{pg}}}}{1+{\gamma }_g\frac{x}{1-x}\frac{c_{\mathit{pp}}}{c_{\mathit{pg}}}}$$

(3)

$$R_e=\frac{1-\varepsilon }{1-x}{R}_g$$

(4)

As particle–gas flows are choked, the isentropic relation between static and total pressure can be used as shown in equation (5)

$$\frac{p^{*}}{p_0}=(\frac{{\gamma }_e+1}{2}{)}^{-{\gamma }_e/{\gamma }_e{\gamma }_e-1{\gamma }_e-1}$$

(5)

The area-Mach number relation is also considered as

$$\frac{A}{A^{*}}=\frac{1}{M}\sqrt{[\frac{2}{{\gamma }_e+1}(1+\frac{{\gamma }_e-1}{2}{M}^2){]}^{{\gamma }_e+1/{\gamma }_e+1{\gamma }_e-1{\gamma }_e-1}}$$

(6)

As p _t  ≤ p*, particle–gas flows are not choked at nozzle throat and mass flow rate of two-phase flows is calculated as shown in equation (7)

$$m^{•}=\frac{P_{0}A}{\sqrt{R_e{T}_0}}\sqrt{\frac{2{\gamma }_e}{{\gamma }_e-1}(\frac{P_0}{P_t}{)}^{-2/2{\gamma }_e{\gamma }_e}-(\frac{P_0}{P_t}{)}^{-{\gamma }_e+1/{\gamma }_e+1{\gamma }_e{\gamma }_e}]}$$

(7)

As pt > p*, particle–gas flows are choked at nozzle throat and mass flow rate of two-phase flows keeps constant if operating conditions at upstream are fixed. The change of downstream flow conditions has on effect on upstream flow of nozzle throat. The mass flow rate of particle–gas flows is calculated as a constant based on equation (8)

$$m^{•}=\frac{P_{0}A}{\sqrt{R_e{T}_0}}\sqrt{{\gamma }_e}(\frac{2}{{\gamma }_e+1}{)}^{{\gamma }_e+1/{\gamma }_e+12({\gamma }_e-1)2({\gamma }_e-1)}$$

(8)

Mass flow rate of particle–gas flows was calculated at different PMLs and compared with mass flow rate of single-gas flow as shown in Figure 1. Mass flow rate of two-phase flows was higher than that for single-gas flow due to the presence of particles. As PML gradually increases, mass flow rate of two-phase flows also increases. Based on choked NPRs at different PMLs, particle–gas flows were much easier to be choked than single-gas flow and particle–gas flows at higher PML were observed to be more easily choked.

Figure 1.

Comparison on theoretical mass flow rate at different PMLs.

Sound speed of mixture

The equilibrium sound speed of two-phase flows is another important parameter to indicate choking phenomena of gas–particle flows. Based on foundations of gas dynamics, for single-gas flow, sound speed of gas phase is calculated by following equation

$$a_g^2=R{\gamma }_{g}T$$

(9)

For gas–particle isentropic flows, sound speed of mixture is obtained by equation (10) that is derived by Rudinger ¹⁰

$$(\frac{a_m}{a_g}{)}^2=\frac{{\gamma }_e}{{\gamma }_g(1+\frac{x}{1-x})(1-\varepsilon {)}^2}$$

(10)

Specific heat ratio of mixture in the above equation is calculated by equation (3). Equation (10) indicates sound speed of mixture increases with the increase of particle volume fraction. At low PML, ɛ is quite low and can be neglected. Sound speed of mixture was plotted at ρ_p/ρ_g = 0.01 and c_pp/c_pg = 1 for neglecting and considering particle volume fraction, respectively, as shown in Figure 2. At low PML, the difference of sound speed of mixture was quite small. As x gradually increased, ɛ also increased and the difference became quite large. As x became 1, solid particles replaced gas–particle flows as working fluid, and sound speed of mixture was infinite when ɛ is considered. However, if ɛ was neglected, sound speed of mixture was approximate to 0.

Figure 2.

Equilibrium sound speed of two-phase flows at ρ_p/ρ_g = 0.01 and c_pp/c_pg = 1.

Numerical methods

Computational domain

Two-dimensional axisymmetric nozzle was used to calculate particle–gas flows for all simulations. A convergent–divergent nozzle was simulated, and the throat diameter of nozzle is 25.4 mm with the design Mach number of 1.5. The detailed size is shown in Figure 3.

Figure 3.

Two-dimensional axisymmetric computational domain.

Numerical schemes

Gambit 2.4.6 was used to generate computational grids. Structure grids were used for full computational domain, and boundary layer meshes were also created near nozzle walls. ANSYS Fluent 15.0 was used as the solver to calculate particle–gas flows. The flow properties were mathematically analyzed by unsteady Reynolds-averaged Navier-Stokes equations. The turbulence model was solved by k-ω shear stress transport and Sutherland viscosity model showing variable viscosity with temperature change was used as viscosity model. Advection Upstream Splitting Method (AUSM) was used as flux type, and second-order upwind schemes were used for spatial discretization. Discrete phase model in which particle phase is regarded as discrete phase and gas phase is continuous phase was used to simulate gas–particle flows, and discrete random walk was considered.

Boundary conditions

The working fluid was assumed to be ideal gas and initialized at constant total temperature of 293.5 K. All nozzle walls were considered as adiabatic walls with constant temperature of 293.5 K. By considering particle–wall interaction, ‘reflect’ boundary condition was used. As particles met the divergent wall of nozzle, particles rebounded without any loss in particle total momentum. y + value was close to 0.5 near nozzle walls. Inlet stagnation pressure was fixed as 3.5 atm. The detailed boundary conditions for computational domain are shown in Figure 3. Particles were seeded after the calculation was convergent using single-gas flow. Anthracite was chosen as seeding particles, and the diameter of particles is 1 µm. Velocity magnitude of seeding particles was kept same as gas velocity at nozzle inlet, and PML was varied from 5 to 30%. The group injection of seeding particles was set at the position of nozzle inlet.

Mesh independent studies

In order to choose the most suitable grids for numerical simulations, mesh independence studies were carried out. Three groups of computational grids were created and simulated at same operating conditions. Mach number distributions were obtained at nozzle axis for different computational grids as shown in Figure 4. Mach number distributions were almost similar, and mesh size of 57,000 was used for all simulations on particle–gas flows through present convergent–divergent nozzle.

Figure 4.

Mach number distributions along nozzle axis for different computational grids.

Results and discussion

Gas and particle velocity contours were obtained at different PMLs as shown in Figure 5(a) to (e). Gas velocity contours were got by simulating single-gas flow at NPR = 0.5714 as shown in Figure 5(a). Shock wave location and structure were clearly seen, and a slip line occurred under the interaction of oblique shock wave and normal shock wave. Supersonic flows induced behind oblique shock wave interacted with subsonic flows developed behind normal shock wave. Slight difference between gas and particle velocity, so particles were regarded to track gas flows properly based on particle distributions. Gas and particle velocity gradually decreased with the increase of PMLs. In addition, shock wave structure and location were also affected by PMLs as shown in Figure 6. Normal shock wave and oblique shock wave showed different structures at different PMLs.

Figure 5.

Particle and gas velocity contours at different PMLs (NPR = 0.5714). (a) Single-gas flow, (b) PML = 5%, (c) PML = 10%, (d) PML = 20% and (e) PML = 30%.

Figure 6.

Particle and gas velocity contours at different PMLs (NPR = 0.5714).

Mass flow rate of particle–gas flows at nozzle throat was theoretically calculated and compared with CFD results at PML of 10% as shown in Figure 7. Mass flow rate was observed to be slightly higher in theoretical calculation compared to that in CFD study. This mainly resulted from that interaction with the continuous phase was considered in numerical simulations. Particle–gas flows were regarded as equilibrium and isentropic flows, and no slip occurred in gas and particle velocity from theoretical calculations. However, gas flow velocity was observed to be slightly different from particle velocity under this interaction effect in CFD study. Particle–gas flows were choked at almost same NPR with back pressure of 3.0 atm in both theoretical and CFD studies.

Figure 7.

Comparison on theoretical and numerical mass flow rate at nozzle throat at PML = 10%.

Figure 8.

Comparison on theoretical and numerical mass flow rate at nozzle throat at NPR = 0.8857.

Mass flow rate of particle–gas flows at nozzle throat was theoretically and numerically calculated at different NPRs and PMLs as shown in Figures 8 to 11. Mass flow rate was observed to be lower in CFD study compared to theoretical calculation as shown in Figures 8 and 10. This is mainly due to that homogeneous equilibrium model was used for theoretical calculation and the slip of velocity and temperature was not considered. In numerical simulation, interaction with continuous phase was considered, which means discrete phase is affected by continuous phase in momentum and energy transfer. Theoretical and numerical mass flow rate of two-phase flows was obtained to be higher that single-gas flow due to the presence of particles. At lower NPR, choked mass flow rate of two-phase flows was observed to be much higher than that for single-gas flow. It indicates that compared to subsonic flows, particle behaviours were more greatly influenced by supersonic flows even though different drag models were used for calculating supersonic and subsonic particle–gas flows, respectively, based on Henderson's conclusion ¹⁷ in present studies.

Figure 9.

Comparison on mass flow rate at nozzle throat at different PMLs (NPR = 0.8857).

Figure 10.

Comparison on theoretical and numerical mass flow rate at nozzle throat at NPR = 0.5714.

Figure 11.

Comparison on mass flow rate at nozzle throat at different PMLs (NPR = 0.5714).

Mass flow rate calculated from theoretical and numerical studies at different PMLs is shown in Figures 9 and 11. Mass flow rate was also observed to be higher in theoretical calculation compared to that in CFD study. It was shown to be higher at higher PMLs as well. Higher PML was initialized by higher particle mass flow rate as particle injection conditions were set. As PMLs became higher, mass flow rate of two-phase flows in CFD study was more approximate to theoretical results.

Sound speed of particle–gas flows was calculated based on equation (10) at different PMLs and NPRs as shown in Figures 12 and 13. Particles were initialized in same velocity and temperature as gas phase at nozzle inlet, so sound speed of gas phase in single-phase flow and two-phase flows was almost same. Due to the presence of particles, sound speed of two-phase flows was observed to be lower than that of single-phase flow. As PML increased, much lower sound speed of two-phase flows was obtained. Based on equation (10), PML and particle volume fraction reduce sound speed of two-phase flows.

Figure 12.

Sound speed of single-phase flow and two-phase flows at NPR = 0.5714. (a) PML = 5%, (b) PML = 10%, (c) PML = 20% and (d) PML = 30%.

Figure 13.

Sound speed of single-phase flow and two-phase flows at NPR = 0.8857. (a) PML = 5%, (b) PML = 10%, (c) PML = 20% and (d) PML = 30%.

As particle–gas flows moved through nozzle convergent section, sound speed of gas phase in single-phase flow and two-phase flows became different. Particle–gas flows were accelerated in nozzle convergent section, and temperature of two-phase flows was higher due to lower gas and particle velocity. The difference between sound speed of gas phase and two-phase flows became larger in supersonic flow regions compared to that in subsonic flow regions, which indicated that two-phase flows were more greatly influenced by particles in supersonic regions. As PMLs increased, sound speed of two-phase flows decreased, and the turbulence of particle–gas flows increased which induced the fluctuation of sound speed of two-phase flows.

Shock wave location and strength were observed to be affected by PML at NPR = 0.5714. Shock wave location lied much further downstream in two-phase flows compared to that in single-phase flow at higher PMLs. As particles were injected into single-gas flow at nozzle inlet, inlet total pressure was influenced by particle momentum and it was increased. Higher NPR induced shock wave location at much further downstream. In addition, shock wave strength was decreased as PMLs were increased.

Particle distributions were obtained at NPR =0.5714 as shown in Figure 14. Shock wave structure and slip line induced by the interaction of oblique shock wave and normal shock wave occurring in nozzle divergent section were clearly seen. Nine lines along flow streamline were created to plot the distribution of PND.

Figure 14.

Different lines for plotting particle number density at NPR = 0.5714. (a) Particle velcoity vectors adjacent to nozzle throat. (b) Particle velocity vectors behind nozzle throat.

Particle velocity vectors were obtained adjacent to and behind nozzle throat by tracking particles with skip ratio of 3 as shown in Figure 15(a) and (b). Particles concentrated along the wall in nozzle convergent section and particle velocity reached to be transonic at nozzle throat due to that the nozzle was choked at NPR = 0.5714. Particles were also observed to concentrate along the slip line. Particles lied upside of slip line were supersonic behind oblique shock wave, and subsonic particles occurred downside of slip line. Interaction between high momentum particles and low momentum particles resulted to high PND along slip line.

Figure 15.

High particle concentration regions at NPR = 0.5714. (a) Particle velcoity vectors adjacent to nozzle throat, (b) Particle velocity vectors behind nozzle throat.

PND was plotted along nine lines at NPR = 0.5714 and NPR = 0.8857, respectively, as shown in Figures 16 and 17. PML was fixed to be constant of 10% for both cases. PND gradually decreased from nozzle axis to nozzle wall or flow free boundary, especially behind nozzle exit. As particles moved from nozzle inlet to nozzle throat (X = 0), PND at each interval gradually increased due to the effect of convergent structure. PND adjacent to nozzle axis decreased behind nozzle throat. PND was observed to be higher in subsonic flows compared to that in supersonic flows at same PML.

Figure 16.

Particle number density at different lines at NPR = 0.5714.

Figure 17.

Particle number density at different lines at NPR = 0.8857.

Conclusions

Particle–gas two-phase flows were theoretically and numerically investigated at different PMLs and NPRs. Mass flow rate of particle–gas flows was observed to be higher compared to that of single-gas flows due to the presence of particles. Mass flow rate at nozzle throat was observed to be higher in theoretical calculation compared to that in CFD study, which mainly due to that equilibrium and isentropic particle–gas flows were considered in theoretical calculation. At higher PML, mass flow rate was also higher than that at lower PML. Sound speed of two-phase flows was obtained to be lower at higher PML. The difference between sound speed of single-gas flow and two-phase flows became larger in supersonic flow regions compared to that in subsonic flow regions, which indicated that two-phase flows were more greatly influenced by particles in supersonic regions. Shock wave location and strength were observed to be strongly determined by PML. Particles were observed to concentrate along the wall in nozzle convergent section and the slip line. PND gradually decreased from nozzle axis to nozzle wall or flow free boundary.