Stochastic separated flow models with applications in numerical computations of supersonic particle-laden turbulent flows

Abstract

In this article, three stochastic separated flow models were applied to predict the dispersion of inertial fuel particles in the supersonic turbulent flows. The flow field of continuous phase was simulated by means of Reynolds-averaged Navier–Stokes method with a two-equation turbulence model. Clift’s expression was used to modify the drag force on the particle considering the compressibility effects. The particle-phase statistics were obtained by a secondary-order time-weighed Eulerian method. The ability of those stochastic separated flow models was then compared for predicting the mean particle velocity and the particle dispersion. For obtaining a statistically stationary solution, the stochastic separated flow model required the largest number of computational particles, whereas the improved stochastic separated flow model was found to need the least. The time-series stochastic separation flow model lay in-between. Compared with the other two models, the particle dispersion was over-predicted by the stochastic separated flow model in the supersonic particle-laden boundary layer flow, while the improved stochastic separated flow model was less predictable for the particle spatial distribution in the particle-laden strut-injection flow. Three models could well predict the mean velocities of the particle phase. This study is valuable for selecting a validated model used for predicting the particle dispersion in supersonic turbulent flows.

Introduction

The low-speed two-phase turbulent flows are commonly encountered in abundant industrial applications, such as energy conversion devices and propulsion systems. These incompressible two-phase flows have been successfully studied by means of computation methods. As the scramjet engine achieves the supersonic combustion, in which the liquid fuel is atomized into spray droplets before the evaporation and the ignition occur, the supersonic particle-laden flows in the scramjet combustor have attracted increasing attention. ¹ Fuel droplets are quickly transported by high-speed air streams and reside in the combustor within several microseconds. Hence, the knowledge based on the dispersion of fuel droplets/particles in supersonic flows is fundamental for illustrating the combustion process in scramjet combustors. ² However, the research on the supersonic particle-laden flows is far from sufficient, compared with the low-speed two-phase flows.

Compared with the expensive experiments, the numerical simulation has been an effective tool to investigate the two-phase flows. When solving the two-phase flow equations, the continuous phase is best represented by an Eulerian description. The particle dispersion in the flow field can be modeled by either the Eulerian description or the Lagrangian description. ³ Since the particle pseudo-fluid model based on the closure of modified kinetic theory of the granular flow has not been adaptable for the supersonic flow, the Lagrangian particle trajectory models are usually utilized by many researchers. Ren et al. ⁴ studied the turbulent dispersion of non-evaporating droplets in the supersonic shear layer and found that the smaller the diameter of the droplets, the more rapidly momentum and heat exchanges between two phases are achieved. The turbulent mixing between fuel droplets and supersonic air streams dominates the gas-phase combustion process.^5,6 Wu et al. ⁷ investigated the characteristics of the droplet dispersion in supersonic shear vortexes based on the Eulerian–Lagrangian numerical simulations. However, the droplet dispersion in high-speed turbulent flows has not been investigated enough.^8,9

The Reynolds-averaged Navier–Stokes (RANS) model is an effective approach to predict the turbulent flows, and the stochastic separated flow (SSF) models are applied to simulate the trajectories of Lagrangian particles. Hennick and Lightstone ¹⁰ presented a comparison of the SSF models for two-phase incompressible flows. However, the applications of the SSF models in compressible/supersonic particle-laden flows have not been sufficiently discussed. The complicated flow structures, such as compression waves and expansion waves, can affect the particle dispersion. Furthermore, the model of the drag force acting on the Lagrangian particles should be modified, considering the compressibility effects. In this investigation, we attempt to show the advantages and disadvantages of three different SSF models through numerical predicting of the gas-particle supersonic flows. We take the supersonic boundary layer flow and the supersonic flow over a strut in a channel as the research objects.

Governing equations and numerical methods

Governing equations of continuous phase

The Lagrangian transport of particles through a continuous carrier gas flow is characterized by the following governing equations after applying the Favre filter

$$\frac{\partial {\overline{W}}_{\mathrm{g}}}{\partial t}+\frac{\partial {\overline{F}}_{\mathrm{g}i}^{(E)}}{\partial x_i}=\frac{\partial {\overline{F}}_{\mathrm{g}i}^{(V)}}{\partial x_i}$$

(1)

where ${\overline{W}}_{\mathrm{g}}=\left(\begin{array}{c}\hfill \tilde{\rho }\hfill \\ \multicolumn{1}{c}{\tilde{\rho }{\overline{U}}_i}\\ \multicolumn{1}{c}{\tilde{\rho }\overline{E}}\end{array}\right)$ , ${\overline{F}}_{\mathrm{gi}}^{(E)}=\left(\begin{array}{c}\hfill \tilde{\rho }{\overline{U}}_i\hfill \\ \multicolumn{1}{c}{\tilde{\rho }{\overline{U}}_i{\overline{U}}_j+p{\delta }_{\mathrm{ij}}}\\ \multicolumn{1}{c}{\tilde{\rho }{\overline{U}}_i(\overline{E}+(p/\tilde{\rho }))}\end{array}\right)$ , and ${\overline{F}}_{\mathrm{gi}}^{(V)}=\left(\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{{\overline{\tau }}_{\mathrm{ij}}}\\ \multicolumn{1}{c}{{\overline{U}}_i{\overline{\tau }}_{\mathrm{ij}}-{\overline{q}}_i}\end{array}\right)$ .

The total energy is given by

$$\overline{E}=\overline{e}+\frac{1}{2}({\overline{U}}^2+{\overline{V}}^2+{\overline{W}}^2)+k$$

(2)

where k is the turbulence kinetic energy and $\overline{e}$ is the internal energy. Hence, the perfect gas state equation can be expressed as

$$p=(\gamma -1)\tilde{\rho }(\overline{e}-k)$$

(3)

The viscous stress is ${\overline{\tau }}_{\mathrm{ij}}={\overline{\tau }}_{\mathrm{ij},\mathrm{L}}+{\overline{\tau }}_{\mathrm{ij},\mathrm{T}}$ , including laminar and turbulent components, and ${\overline{q}}_i=-\kappa \partial \overline{T}/\partial x_i$ is the heat conduction. The viscosity is $\mu ={\mu }_{\mathrm{L}}+{\mu }_{\mathrm{T}}$ and the thermal conductivity is $\kappa ={\kappa }_{\mathrm{L}}+{\kappa }_{\mathrm{T}}$ . Both the turbulence viscosity and the turbulent thermal conductivity are calculated by means of turbulence closures.

The turbulence model of k–ε proposed by Hwang and Lin ¹¹ is employed as follows

$$\frac{\partial \tilde{\rho }k}{\partial t}+\frac{\partial \tilde{\rho }{\overline{U}}_{j}k}{\partial x_j}=\frac{\partial }{\partial x_j}\left(({\mu }_{\mathrm{L}}+\frac{{\mu }_{\mathrm{T}}}{{\sigma }_k})\frac{\partial k}{\partial x_j}\right)-\frac{1}{2}\frac{\partial }{\partial x_j}\left(\frac{{\mu }_{\mathrm{L}}k}{\epsilon }\frac{\partial \hat{\epsilon }}{\partial x_j}\right)+{\overline{S}}_k$$

(4)

$$\frac{\partial \tilde{\rho }\tilde{\epsilon }}{\partial t}+\frac{\partial \tilde{\rho }{\overline{U}}_j\tilde{\epsilon }}{\partial x_j}=\frac{\partial }{\partial x_j}\left(({\mu }_{\mathrm{L}}+\frac{{\mu }_{\mathrm{T}}}{{\sigma }_{\epsilon }})\frac{\partial \tilde{\epsilon }}{\partial x_j}\right)-\frac{\partial }{\partial x_j}\left(\frac{{\mu }_{\mathrm{L}}\tilde{\epsilon }}{k}\frac{\partial k}{\partial x_j}\right)+{\overline{S}}_{\tilde{\epsilon }}$$

(5)

The reduction turbulence dissipations is $\tilde{\epsilon }=\epsilon -\hat{\epsilon }$ where $\hat{\epsilon }=2\tilde{\rho }{\mu }_{\mathrm{L}}{(\nabla \sqrt{k})}^2$ . The turbulence viscosity is then calculated as

$${\mu }_{\mathrm{T}}=\tilde{\rho }\left(\frac{k^2}{\tilde{\epsilon }}\right)C_{\mu }{f}_{\mu }(y_{\lambda })$$

(6)

where $y_{\lambda }=s/{\lambda }_{\mathrm{t}}$ . ${\lambda }_{\mathrm{t}}=\sqrt{\tilde{\rho }{\mu }_{\mathrm{L}}k/\tilde{\epsilon }}$ is the Taylor scale and $s$ is the minimum distance to the near wall. $f_{\mu }=1-\exp (-0.01y_{\lambda }-0.008y_{\lambda }^3)$ , ${\sigma }_k=1.4-1.1\exp [-(y_{\lambda }/10)]$ and ${\sigma }_{\epsilon }=1.3-1.0\exp [-(y_{\lambda }/10)]$ .

The source terms in equations (4) and (5) are ${\overline{S}}_k={\tau }_{\mathrm{ij},\mathrm{L}}\partial {\overline{U}}_i/\partial x_j-\tilde{\rho }(\hat{\epsilon }+\tilde{\epsilon })$ and ${\overline{S}}_{\tilde{\epsilon }}=C_{\epsilon 1}\tilde{\epsilon }/k{\tau }_{\mathrm{ij},\mathrm{T}}\partial {\overline{U}}_i/\partial x_j-C_{\epsilon 2}\tilde{\rho }{\tilde{\epsilon }}^2/k$ , respectively.

The empirical closure coefficients in the above models are taken as $C_{\mu }=0.09$ , $C_{\epsilon 1}=1.44$ , and $C_{\epsilon 2}=1.92$ , respectively.

The compressibility correction model postulated by Wilcox ¹² is applied for the compressibility effects, and the turbulent kinetic energy dissipation rate is given by

$$\rho \tilde{\epsilon }=\rho {\tilde{\epsilon }}_{\mathrm{s}}+\rho {\tilde{\epsilon }}_{\mathrm{d}}$$

(7)

where ${\tilde{\epsilon }}_{\mathrm{s}}$ is the solenoidal dissipation rate, and the dilatation dissipation rate ${\tilde{\epsilon }}_{\mathrm{d}}$ is defined as

$${\tilde{\epsilon }}_{\mathrm{d}}={\xi }^{*}F({\mathrm{M}}_{\mathrm{T}}){\tilde{\epsilon }}_{\mathrm{s}}$$

(8)

where ${\mathrm{M}}_{\mathrm{T}}=\sqrt{2k/a^2}$ is the turbulence Mach number and the empirical coefficient ${\xi }^{*}=1.5$ . $F({\mathrm{M}}_{\mathrm{T}})$ is expressed as follows

$$F({\mathrm{M}}_{\mathrm{T}})=[{\mathrm{M}}_{\mathrm{T}}^2-{\mathrm{M}}_{\mathrm{T}0}^2]H({\mathrm{M}}_{\mathrm{T}}-{\mathrm{M}}_{\mathrm{T}0})$$

(9)

where the function H is the Heaviside step function and ${\mathrm{M}}_{\mathrm{T}0}=0.25$ .

A finite difference methodology is used to solve the governing equations of gas phase. A two-step explicit Runge–Kutta time-integration methodology is applied, obtaining a second-order time-accurate computation. The inviscid flux is discretized by utilizing the Roe-type Riemann solver and the second-order spatial accuracy is obtained by the Monotonic Upstream Centered Scheme for Conservation Laws (MUSCL)-type scheme. A modified Harten–Hyman entropy condition is used to avoid unrealistic solutions caused by the Roe-type Riemann solver. The viscous flux is discretized by a second-order central difference scheme.

Governing equations of discrete phase

The particles are tracked individually in a Lagrangian manner. It is assumed that the density of the particles is much larger than that of the continuous phase such that only the drag force is significant. The particle collision and dense particulate effects are neglected. In addition, all the particles are assumed to collide with the walls elastically.

The Lagrangian particle equations for the position and the velocity are given by

$$\{\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{X}}_{\mathrm{p}}}{\mathrm{d}t}={\mathit{U}}_{\mathrm{p}}\hfill \\ \multicolumn{1}{c}{\frac{\mathrm{d}{\mathit{U}}_{\mathrm{p}}}{\mathrm{d}t}=\frac{C_{\mathrm{D}}{\mathrm{Re}}_{\mathrm{r}}}{24{\tau }_{\mathrm{rp}}}({\mathit{U}}_{\mathrm{g}@\mathrm{p}}-{\mathit{U}}_{\mathrm{p}})}\end{array}$$

(10)

Here, X _p is the particle position vector, U _p is the particle velocity vector and U _g@p is the gas velocity seen by the particle. ${\mathrm{Re}}_{\mathrm{r}}=\rho d_{\mathrm{p}}|U_{\mathrm{p}}-U_{\mathrm{g}@\mathrm{p}}|/{\mu }_{\mathrm{L}}$ is the particle Reynolds number and ${\tau }_{\mathrm{rp}}={\rho }_{\mathrm{p}}{d}_{\mathrm{p}}^2/(18{\mu }_{\mathrm{L}})$ is the aerodynamics relaxation time, where d _p is the particle diameter. The drag coefficient, C _D, corrected by Schiller and Naumann, ¹³ is defined as

$$C_{\mathrm{D}}=\frac{24}{{\mathrm{Re}}_{\mathrm{r}}}(1+0.15{\mathrm{Re}}_{\mathrm{r}}^{0.687})$$

(11)

A fourth-order Lagrange interpolation procedure is employed to obtain the gas-phase velocities at the particle locations as the particles do not locate at the grid points. When the particle is located at $(X_{\mathrm{p}},Y_{\mathrm{p}},Z_{\mathrm{p}})$ inside the computational cells where coordinates are $(x_0,y_0,z_0)\dots (x_l,y_l,z_l)$ , the Lagrangian interpolation for the particle velocity can be expressed as

$${\overline{\mathit{U}}}_{\mathrm{g}@\mathrm{p}}(X_{\mathrm{p}},Y_{\mathrm{p}},Z_{\mathrm{p}})={\displaystyle \sum _{i=1-n/2}^{n/2}}{\displaystyle \sum _{j=1-n/2}^{n/2}}{\displaystyle \sum _{k=1-n/2}^{n/2}}\left({\displaystyle \underset{\genfrac{}{}{0ex}{}{l=1-n/2}{(l\ne i)}}{\overset{n/2}{\Pi }}}\frac{X_{\mathrm{p}}-x_l}{x_i-x_l}{\displaystyle \underset{\genfrac{}{}{0ex}{}{l=1-n/2}{(l\ne j)}}{\overset{n/2}{\Pi }}}\frac{Y_{\mathrm{p}}-y_l}{y_j-y_l}{\displaystyle \underset{\genfrac{}{}{0ex}{}{l=1-n/2}{(l\ne k)}}{\overset{n/2}{\Pi }}}\frac{Z_{\mathrm{p}}-z_l}{z_k-z_l}\right){\overline{\mathit{U}}}_{i,j,k}$$

(12)

where the subscripts i, j, k, and l represent the cell indexes, different from those used in the previous sections, and n = 6.

The Lagrangian particle equations are integrated by a third-order Adams–Bashforth approach. The integration of the particle position is given by

$${\mathit{X}}_{\mathrm{p}}^{n+1}={\mathit{X}}_{\mathrm{p}}^n+\Delta t(\frac{23}{12}{\mathit{U}}_{\mathrm{p}}^n-\frac{16}{12}{\mathit{U}}_{\mathrm{p}}^n+\frac{5}{12}{\mathit{U}}_{\mathrm{p}}^{n-2})$$

(13)

and the particle velocity is integrated as

$${\mathit{U}}_{\mathrm{p}}^{n+1}={\mathit{U}}_{\mathrm{p}}^n+\Delta t(\frac{23}{12}{\mathbf{\varphi }}_{\mathrm{p}}^n-\frac{16}{12}{\mathbf{\varphi }}_{\mathrm{p}}^n+\frac{5}{12}{\mathbf{\varphi }}_{\mathrm{p}}^{n-2})$$

(14)

where $\mathbf{\varphi }$ denotes the right-hand side of the equation of motion for the particle and is divided by the particle mass, $m_{\mathrm{p}}$ .

Eulerian statistics of discrete phase

The discrete particles are tracked individually in a Lagrangian manner. Therefore, an Eulerian stationary statistics would be achieved when sufficient particles are calculated. The mean velocity of the discrete phase in each control volume is obtained by a second-order time-weighed Eulerian statistical approach

$${\overline{\mathit{U}}}_{\mathrm{p}}(i,j,k)=\frac{{\displaystyle \sum _m}{\mathit{u}}_{\mathrm{p},m}^L\Delta t_{(i,j,k),m}}{{\displaystyle \sum _m}\Delta t_{(i,j,k),m}}$$

(15)

where Δt is the characteristic time for the mth discrete particle residing in the computational volume, as shown in Figure 1.

Figure 1.

Schematic of the Eulerian statistical method.

Inter-phase coupling and compressibility effects

The one-way coupling is employed in this study, which refers to the dilute two-phase flows with relatively small particle mass loading ratios, and the presence of the particles does not significantly affect their surrounding continuous phase.

The effects of turbulence on the particle motion are concerned. In addition, the particle motion in the supersonic flows can be affected by the high compressibility. The effects of compressibility on the particle motion are considered via the correction of the drag coefficient

$$C_{\mathrm{D},\mathrm{Kn},{\mathrm{Re}}_{\mathrm{r}}}=\frac{24}{{\mathrm{Re}}_{\mathrm{r}}}(1+0.15{\mathrm{Re}}_{\mathrm{r}}^{0.687})f(\mathrm{Kn})$$

(16)

Clift et al.’s ¹⁴ expression, considering the compressibility effects, is given by

$$f(\mathrm{Kn})=\frac{1}{1+\mathrm{Kn}(2.514+0.8\exp (-0.55/\mathrm{Kn}))}$$

(17)

where $\mathrm{Kn}=\sqrt{\pi \gamma /2}({\mathrm{M}}_{\mathrm{r}}/{\mathrm{Re}}_{\mathrm{r}})$ is the Knudsen number and ${\mathrm{M}}_{\mathrm{r}}=|{\mathit{U}}_{\mathrm{p}}-{\mathit{U}}_{\mathrm{g}@\mathrm{p}}|/a_{\mathrm{g}}$ is the relative Mach number.

The particle spreading in a uniform supersonic air flow (U = 1000 m/s and T = 300 K) is utilized to illustrate the effects of compressibility on the drag coefficient, as shown in Figure 2. The particle density is 800 kg/m³ and the particle diameters are 2 and 20 µm, respectively. For the momentum relaxation of particles, the compressibility decreases the drag coefficient, and the compressibility effects become stronger with a decrease in the two-phase velocity slip. The two-phase velocity slip becomes smaller as the particle tends to be in the relaxation equilibrium.

Figure 2.

Temporal variation in the drag coefficient.

SSF models

Along the trajectory of Lagrangian particles, the continuous phase velocity seen by the discreet particle is decomposed into a time-averaged velocity and a fluctuating velocity

$${\mathit{U}}_{\mathrm{g}@\mathrm{p}}={\overline{\mathit{U}}}_{\mathrm{g}@\mathrm{p}}+\mathit{u}{\prime }_{\mathrm{g}}$$

(18)

Turbulent eddies create the velocity fluctuations, which is modeled via using the Monte Carlo method. The velocity fluctuations of the gas phase are generated from the random sampling of a zero-mean Gaussian probability density function with a standard deviation, ${\sigma }^2=2/3k=\overline{{{\mathit{u}\prime }_{\mathrm{g}}}^2}$ , that is

$$\mathit{u}{\prime }_{\mathrm{g}}=\xi (\overline{{{\mathit{u}\prime }_{\mathrm{g}}}^2})$$

(19)

The interval time during which the particle interacts with the local turbulence is identified by two time scales, that is, the eddy lifetime, ${\tau }_{\mathrm{fl}}$ , and the transit time required for the particle to cross the eddy, ${\tau }_{\mathrm{p}}$ , respectively. If the particle moves almost as fast as the local fluid, it could be captured by the turbulent eddy and remains in the eddy during the whole eddy lifetime. Whereas if there is a significant velocity slip between the two phases, the particle can cross the eddy before the eddy decays. Therefore, the interaction time, $\Delta \tau$ , is determined by the minimum of the two time scales, that is

$$\Delta \tau =min({\tau }_{\mathrm{fl}},{\tau }_{\mathrm{p}})$$

(20)

The eddy lifetime is then evaluated by

$${\tau }_{\mathrm{fl}}=\frac{l_{\mathrm{e}}}{\left|{\mathit{u}\prime }_{\mathrm{g}}\right|}$$

(21)

where the dissipation length scale of the eddy is given by

$$l_{\mathrm{e}}=C_{\mu }^{0.75}\frac{k^{1.5}}{\epsilon }$$

(22)

The transit time required for the particle crossing the eddy is expressed as

$${\tau }_{\mathrm{p}}=\{\begin{array}{c}\hfill {\tau }_{\mathrm{fl}},l_{\mathrm{e}}>{\tau }_{\mathrm{rp}}|{\mathit{U}}_{\mathrm{g}@\mathrm{p}}-{\mathit{U}}_{\mathrm{p}}|\hfill \\ \multicolumn{1}{c}{-{\tau }_{\mathrm{rp}}\ln (1-\frac{l_e}{{\tau }_{\mathrm{rp}}|{\mathit{U}}_{\mathrm{g}@\mathrm{p}}-{\mathit{U}}_{\mathrm{p}}|})}\end{array},l_{\mathrm{e}}\le {\tau }_{\mathrm{rp}}|{\mathit{U}}_{\mathrm{g}@\mathrm{p}}-{\mathit{U}}_{\mathrm{p}}|$$

(23)

Although several researchers^15–17 have proposed corrections for the SSF models, the difficulty of determining the action time of a given fluctuating continuous phase velocity seen by the particle still remains unsolved. The SSF model has long been used in predicting the two-phase turbulent flows because of its simplicity and robustness.

An improved stochastic separated flow (ISSF) model, considering the intermittent action of the turbulence, has been proposed to obtain a reasonable statistical characteristic of the dispersed phase. In this model proposed by Chan et al., ¹⁸ the velocity of the dispersed phase, including the mean velocity and the fluctuating velocity, is transported along its stochastic trajectory

$${\mathit{X}}_p=\int ({\mathit{U}}_p+\mathit{u}{\prime }_{\mathrm{p}})\mathrm{dt}$$

(24)

where the particle fluctuating velocity, $\mathit{u}{\prime }_{\mathrm{p}}$ , is obtained by a random number generator, ξ, obeying a zero-mean Gaussian probability density function with a standard deviation, ${\sigma }^2=\overline{{\mathit{u}{\prime }_{\mathrm{p}}}^2}$ , that is

$$\mathit{u}{\prime }_{\mathrm{p}}=\xi (\overline{{{\mathit{u}\prime }_{\mathrm{p}}}^2})$$

(25)

In order to involve the temporal variation in the particle fluctuating velocity, a Lagrangian transport equation is given as follows

$$\frac{\mathrm{d}\overline{{{\mathit{u}\prime }_{\mathrm{p}}}^2}}{\mathrm{d}t}=\frac{2\overline{{\mathit{u}\prime }_{\mathrm{p}}{u\prime }_{\mathrm{g}}}}{{\tau }_{\mathrm{rp}}}-\frac{2\overline{{{\mathit{u}\prime }_{\mathrm{p}}}^2}}{{\tau }_{\mathrm{rp}}}$$

(26)

where the two-phase fluctuation correlation term is closed by

$$\overline{{\mathit{u}\prime }_{\mathrm{p}}{\mathit{u}\prime }_{\mathrm{g}}}=\overline{{{\mathit{u}\prime }_{\mathrm{g}}}^2}\exp (-B_k\frac{{\tau }_{\mathrm{rp}}}{{\tau }_{\mathrm{T}}})$$

(27)

here, B_k  = 0.5 and ${\tau }_{\mathrm{T}}=C_{\mathrm{T}}k/\epsilon$ where C _T = 0.165. Wang et al. ¹⁹ have compared and evaluated different closure models of the two-phase fluctuating velocity correlation terms.

Due to the consideration of the intermittent action of the turbulence, the ISSF model predicts good results for the particle fluctuating velocity and uses less sample particles, compared with the SSF model. However, the two conventional SSF models, including the SSF and ISSF models, predict the concerned particle dispersion poorly. Therefore, the time-series stochastic separation flow (TSSSF) model has been developed. ²⁰

In the TSSSF model, the instantaneous fluctuating velocity of the gas phase, ${\mathit{u}{\prime }_{\mathrm{g}}}^{(n)}$ , is sampled using the Monte Carlo method at the calculation time step n. It is kept to act on the change in the particle instantaneous velocity until its action can be neglected according to the real-time calculation of the auto-correlation coefficient between the ${\mathit{u}{\prime }_{\mathrm{g}}}^{(n)}$ and ${\mathit{u}{\prime }_{\mathrm{g}}}^{(n-1)}$ . It means that the TSSSF model considers the temporal and spatial correlations of the continuous phase fluctuating velocity, achieved by the low-order model

$${\mathit{U}}_{\mathrm{g}@\mathrm{p}}^{(n)}={\overline{\mathit{U}}}_{\mathrm{g}@\mathrm{p}}^{(n)}+\mathit{u}{\prime }_g^{(n)}-{\mathbf{\theta }}_{\mathit{u}}^{\mathbf{(}\mathbf{n}\mathbf{)}}\mathit{u}{\mathbf{\text{'}}}_{\mathbf{g}}^{\mathbf{(}\mathbf{n}\mathbf{-}1\mathbf{)}}$$

(28)

$${\mathit{u}{\prime }_{\mathrm{g}}}^{(n)}=\xi \left[{\overline{{{\mathit{U}\prime }_{\mathrm{g}}}^2}}^{(n)}\right]$$

(29)

$${\theta }_u=\frac{\sqrt{{{\mathit{u}\prime }_{\mathrm{g}}}^{4(n-1)}-4{\mathbf{\rho }}_u^{2(n)}\mathit{u}{\text{'}}_g^{(n-1)}{{\mathit{u}\text{'}}_g}^{2(n)}-{{\mathit{u}\text{'}}_g}^{2(n-1)}}}{2{\mathbf{\rho }}_u^{(n)}{{\mathit{u}\text{'}}_{\mathrm{g}}}^{(n-1)}}$$

(30)

where ${\mathbf{\rho }}_u^{(n)}$ is the random number and distributes uniformly at the time step n

$${\mathbf{\rho }}_u^{(n)}=\mathbf{\zeta }[-\frac{{{\mathit{u}\text{'}}_{\mathrm{g}}}^{3/2(n-1)}}{2{{\mathit{u}\text{'}}_{\mathrm{g}}}^{(n)}},\frac{{{\mathit{u}\text{'}}_{\mathrm{g}}}^{3/2(n-1)}}{2{{\mathit{u}\text{'}}_{\mathrm{g}}}^{(n)}}]$$

(31)

Table 1 lists all the governing equations of the three SSF models.

Table 1.

SSF, ISSF, and TSSSF models.

SSF model	ISSF model	TSSSF model
$\frac{\mathrm{d}{U}_{i,\mathrm{p}}}{\mathrm{d}t}=\frac{{\overline{U}}_{i,\mathrm{g}@\mathrm{p}}+{u\prime }_{i,\mathrm{g}}-U_{i,\mathrm{p}}}{{\tau }_{\mathrm{r},\mathrm{p}}}$	$\frac{\mathrm{d}{U}_{i,\mathrm{p}}}{\mathrm{d}t}=\frac{{\overline{U}}_{i,\mathrm{g}@\mathrm{p}}-U_{i,\mathrm{p}}}{{\tau }_{\mathrm{r},\mathrm{p}}}$	$\frac{\mathrm{d}{U}_{i,\mathrm{p}}}{\mathrm{d}t}=\frac{U_{i,\mathrm{g}@\mathrm{p}}-U_{i,\mathrm{p}}}{{\tau }_{\mathrm{r},\mathrm{p}}}$
$\frac{\mathrm{d}{X}_{i,\mathrm{p}}}{\mathrm{d}t}=U_{i,\mathrm{p}}$	$\frac{\mathrm{d}{X}_{i,\mathrm{p}}}{\mathrm{d}t}=U_{i,\mathrm{p}}+u{\prime }_{i,\mathrm{p}}$	$\frac{\mathrm{d}{X}_{i,\mathrm{p}}}{\mathrm{d}t}=U_{i,\mathrm{p}}$
$u{\prime }_{i,\mathrm{g}}={\xi }_i(\overline{u{\prime }_{i,g}^2})$	$\frac{\mathrm{d}\overline{u{\prime }_{i,\mathrm{p}}^2}}{\mathrm{d}t}=\frac{2\overline{{u\prime }_{i,\mathrm{p}}{u\prime }_{i,\mathrm{g}}}}{{\tau }_{\mathrm{rp}}}-\frac{2\overline{u{\prime }_{i,\mathrm{p}}^2}}{{\tau }_{\mathrm{rp}}}$ $u{\prime }_{i,\mathrm{p}}={\xi }_i(\overline{u{\prime }_{i,p}^2})$ $\overline{{u\prime }_{i,\mathrm{p}}{u\prime }_{i,\mathrm{g}}}=\overline{u{\prime }_{i,\mathrm{g}}^2}\exp (-B_k\frac{{\tau }_{\mathrm{rp}}}{{\tau }_{\mathrm{T}}})$	$U_{i,\mathrm{g}@\mathrm{p}}^{(n)}={\overline{U}}_{i,\mathrm{g}@\mathrm{p}}^{(n)}+u{\prime }_{i,\mathrm{g}}^{(n)}-{\theta }_{i,u}^{(n)}u{\prime }_{i,\mathrm{g}}^{(n-1)}$ $u{\prime }_{i,\mathrm{g}}^{(n)}={\xi }_{i,u}({\overline{u{\prime }_{i,\mathrm{g}}^2}}^{(n)})$ ${\theta }_{i,u}=\frac{\sqrt{u{\prime }_{i,\mathrm{g}}^{4(n-1)}-4{\rho }_{i,u}^{2(n)}u{\prime }_{i,\mathrm{g}}^{(n-1)}u{\prime }_{i,\mathrm{g}}^{2n}-u{\prime }_{i,\mathrm{g}}^{2(n-1)}}}{2{\rho }_{i,u}^{(n)}u{\prime }_{i,\mathrm{g}}^{(n-1)}}$ ${\rho }_{i,u}^{(n)}={\zeta }_i[-\frac{u{\prime }_{i,\mathrm{g}}^{(n-1)}}{2u{\prime }_{i,\mathrm{g}}^{(n)}}\sqrt{u{\prime }_{i,\mathrm{g}}^{(n-1)}},\frac{u{\prime }_{i,\mathrm{g}}^{(n-1)}}{2u{\prime }_{i,\mathrm{g}}^{(n)}}\sqrt{u{\prime }_{i,\mathrm{g}}^{(n-1)}}]$

SSF: stochastic separated flow; ISSF: improved stochastic separated flow; TSSSF: time-series stochastic separation flow.

Results and discussions

Supersonic particle-laden boundary layer flow

The transverse spray-jet is one of the main fuel injection approaches in the scramjet combustor. As shown in Figure 3, the liquid fuel is injected into the supersonic turbulent boundary layer (TBL) flow and atomizes into the droplets because of the primary and secondary breakup mechanisms.^21,22

Figure 3.

Sketch of the injection of droplets into boundary layer flow: (a) transverse-jet flow and (b) simplified particle-laden boundary layer flow.

The atomization process is not concerned in this study. Therefore, the whole flow prototype, shown in Figure 3(a), is simplified as the supersonic particle-laden boundary layer flow, shown in Figure 3(b). Moreover, the flow is assumed to be fully developed due to the high speed of the inflow. The computation domain is specified as L_x  = 244 mm and L_y  = 5.57 mm in the streamwise and wall-normal directions, respectively.

The calculation parameters of two-phase flow are shown in Table 2. The droplet material density is taken as 800 kg/m³. The sizes of spray droplets are very small due to the high efficiency of the atomization in supersonic streams. Here, spray droplets are assumed as sphere particles with the same diameter of 20 µm. All particles are vertically injected into the supersonic flow with the same wall-normal velocity equal to 20 m/s. The air inflow Mach number is 2.25 and the flow Reynolds number is 25,000 per unit millimeter.

Table 2.

Two-phase flow parameters for simulation TBL.

Continuous phase		Discrete phase
Fluid medium	Air	Droplet material density ρ p (kg/m3)	800
Inflow Mach number	2.25	Diameter d p (µm)	20
Flow Reynolds number/mm	25,000	Injection velocity V j (m/s)	20

Figure 4 shows the distributions of velocity components, turbulence kinetic energy, and turbulence kinetic dissipation rate, respectively. It is found that the supersonic TBL is very thin. The turbulence fluctuations mainly exist close to the wall and show a very high level of turbulence kinetic energy.

Figure 4.

Continuous phase flow field of simulation TBL: (a) streamwise velocity $\overline{U}$ , (b) wall-normal velocity $\overline{V}$ , (c) turbulence kinetic energy k, and (d) turbulence kinetic energy dissipation rate ε.

Figure 5 shows the particle dispersion in the boundary layer. The contours represent the particle number density, which is obtained statistically by

$$C(x,y,z)=\frac{n_{\mathbf{cell}}(x,y,z)}{n_{\mathbf{total}}}$$

(32)

Figure 5.

Particle spatial dispersion predicted by different SSF models in TBL: (a) SSF model, (b) ISSF model, and (c) TSSSF model.

where $n_{\mathbf{cell}}(x,y,z)$ is the number of computational particles in one cell volume at the position (x, y, z) and $n_{\mathbf{total}}$ is the total number of computational particles.

It is shown that the dispersion is over-predicted by means of the SSF model. The results predicted by the ISSF and TSSSF models are consistent with each other, and the spatial distributions of particles are more concentrated, compared with the prediction by the SSF model.

At the profiles of y/D = 1 and 2, closing to the injection inlet, the velocities of discrete particles predicted by three models are the same in the wall-normal and streamwise directions, respectively. Furthermore, if the same number of particles is employed in the three models, the smoother statistics of particle velocities can be obtained by the ISSF model, compared with that of the other two models. The particle-phase velocities predicted by the three SSF models are different in the downstream region. Especially for the conventional SSF model, the statistical velocities are distributed extensively in the streamwise direction due to the large streamwise velocity of the particles (Figure 6).

Figure 6.

Predicted particle-phase velocities by different SSF models in TBL: (a) wall-normal component of particle-phase velocity and (b) streamwise component of particle-phase velocity.

For different particle sizes, the profiles of the predicted velocity at y/D = 4 are shown in Figure 7. With the increase in the particle diameter, particles accelerate more slowly in both wall-normal and streamwise directions and need more time to catch up with the surrounding gas. Hence, large-sized particles disperse in a relative narrow range and the three models predict nearly the same velocity distributions for the particles with the diameter of 50 µm. The penetration depth is a very important measurement parameter in the study of the high-speed jet flow and affects the mixing and atomization processes. The wall-normal velocity of the particle phase represents the penetration depth, and then, the large-sized particles with high wall-normal speed have a deep penetration depth, since the large-sized ones have long aerodynamic response time and motion almost with their initial speed. However, large-sized particles distribute in a small range along the streamwise direction since the initial streamwise velocity is 0, and the weak dispersion could hinder the mixing between the fuel droplets and the air.

Figure 7.

Predictions of the particle-phase velocity at y/D = 4 in TBL: (a) wall-normal component of particle-phase velocity and (b) streamwise component of particle-phase velocity.

The required particle numbers for obtaining the above steady statistical results are shown in Table 3. The ISSF model is found to need around 1000 computational particles to obtain a statistically stationary solution of the mean particle velocities, while the SSF model needs almost 6000 computational particles. The TSSSF model requires less than 5000 computational particles, which are much larger than that of the SSF model but still less than that of the ISSF model.

Table 3.

Particle number for different SSF models in TBL.

SSF model	Particle number
SSF	~6000
ISSF	~1000
TSSSF	~5000

SSF: stochastic separated flow; ISSF: improved stochastic separated flow; TSSSF: time-series stochastic separation flow.

Strut-injection particle-laden flow

The strut injection (SJ) is another important method for the fuel injection and flame holding in the scramjet combustor, as shown in Figure 8. In this study, we focus on the flow downstream the strut, and the computation domain is specified as L_x  = 300 mm and L_y  = 100 mm in the streamwise and wall-normal directions, respectively.

Figure 8.

Sketch of strut-injection particle-laden flow.

The calculation parameters of the two-phase flow are shown in Table 4. The inflow conditions are specified according to the experiment by Deutsches Zentrum für Luft- und Raumfahrt (DLR). ²³ The inflow velocity is 732 m/s, the static pressure 100 kPa, and the static temperature is 340 K. The strut has a wedge shape with an apex angle of 12° and a length of L = 32 mm.

Table 4.

Two-phase flow parameters for simulation SJ.

Continues phase		Discrete phase
Fluid medium	Air	Droplet material density ρ p (kg/m3)	800
Inflow Mach number	2.0	Diameter d p (µm)	20
Flow Reynolds number	270,000	Injection velocity U j (m/s)	100

The non-evaporating droplets (kerosene) are injected from the strut bottom plate with the streamwise velocity equal to 100 m/s. The particles are assumed to have the same diameter of 20 µm and the particle density is specified as 800 kg/m³. The inflow Mach number is 2.0. The flow Reynolds number is 270,000 based on the wedge thickness (H = 6.73 mm).

Figure 9 shows the flow field of the continuous phase downstream the strut. The supersonic inflow is deflected at the tip of the strut, and the two shocks are formed symmetrically on each side. After impacting the upper and down walls, the shock waves reflect back toward the centerline of the channel. The λ-shock waves are observed in the downstream region. On the other hand, the supersonic flow forms expansion waves at the strut bottom corners. The expansion waves interact with the shock waves, and they both interact with the shear layers generated by the strut. The recirculation zone forms behind the strut and could hold the flame stability in the supersonic combustion process.

Figure 9.

Continuous phase flow field of simulation SJ: (a) streamwise velocity $\overline{U}$ , (b) transverse velocity $\overline{V}$ , (c) turbulence kinetic energy k, and (d) turbulence kinetic energy dissipation rate ε.

It is found that the turbulence fluctuations mainly exist in the shear regions, which are illustrated by the high level of turbulence kinetic energy and the corresponding large turbulence kinetic dissipation rate, as shown in Figure 9(c) and (d).

Figure 10 shows the particle dispersion predicted by three different SSF models, and these figures focus on the particle dispersion in the recirculation region downstream the strut. The particles with the size of 20 µm have a short aerodynamic response time and follow the continuous phase quickly. Furthermore, the highest concentration of the particles is found in the recirculation region. The particles disperse widely downstream the strut due to the transportation of the supersonic streams. Three models predict similar spatial distributions of the particles, and therefore, the particle concentrations are similar. It is also found that the particles mainly accumulate in the shear layer. The concentrated distributions of the particles predicted by the SSF and TSSSF models are consistent with each other. Whereas the ISSF model predicts the particle dispersion more weakly, which is different from those predicted by the other two models.

Figure 10.

Particle spatial dispersion predicted by different SSF models in SJ: (a) SSF model, (b) ISSF model, and (c) TSSSF model.

The particle spatial distributions downstream the strut are depicted in Figure 11. The particles spread out of the recirculation regions and diffuse transversely. The SSF model predicts a similar particle spatial distribution, compared with that of the TSSSF model, whereas the ISSF model reports that the particles accumulate densely in the particular regions and not as uniformly as the results from other models. Figure 12 shows the statistical velocities of the particle phase at four different profiles downstream the strut. The predictions of the three models are almost the same. In the ISSF model, the turbulence kinetic energy transport equation of particle phase is solved along the particle path. Hence, the particle velocity can be obtained more smoothly with less sample particle amount in the statistics than the other two models. It is also noticed that the particle velocities in both streamwise and transverse directions yield self-similarity.

Figure 11.

Particle spatial dispersion downstream the recirculation region: (a) SSF model, (b) ISSF model, and (c) TSSSF model.

Figure 12.

Predicted particle-phase velocities by different SSF models in SJ: (a) streamwise component of particle-phase velocity and (b) transverse component of particle-phase velocity.

Compared with the supersonic particle-laden boundary layer flow, the SSF model needs more computational particles to obtain the steady statistical results in the SJ particle-laden flow, while the ISSF and TSSSF models remain almost unchanged, as shown in Table 5. Hence, the SSF model is more sensitive for the required tracking particles than the other two models.

Table 5.

Particle number for different SSF models in SJ.

SSF model	Particle number
SSF	~8000
ISSF	~1000
TSSSF	~5000

SSF: stochastic separated flow; ISSF: improved stochastic separated flow; TSSSF: time-series stochastic separation flow.

Conclusion

The accurate numerical simulations of supersonic particle-laden flows are the emerging topics, companying with the development of scramjet technologies. This study is aimed at the numerical solution of the discrete particles by successfully utilizing three different SSF models. Two simplification prototype flows, involved in fuel injection approaches for the scramjet combustor, that is, the supersonic particle-laden boundary layer flow and the SJ particle-laden flow, are numerically simulated by means of RANS coupled with the SSF, ISSF, and TSSSF models. The statistical particle-phase velocities and the particle spatial dispersions are obtained.

The particle dispersion in the supersonic boundary layer is over-predicted via the SSF model. Therefore, the statistical distribution of particle-phase velocities spreads widely in the streamwise direction downstream the injection. Large-sized particles concentrate narrowly along the streamwise direction but penetrate deeply. The ISSF model underestimates the particle dispersion in the SJ supersonic flow. In addition, three particle trajectory-tracking models predict similar statistical velocities of the particle phase downstream the strut. The ISSF model is capable of obtaining a statistically stationary solution with very few computational particles, whereas the SSF model requires the largest number of computation particles among the three models.