Investigation of the Horizontal Motion of Particle-Laden Jets

Abstract

Particle-laden jet flows can be observed in many industrial applications. In this investigation, the horizontal motion of particle laden jets is simulated using the Eulerian–Lagrangian framework. The two-way coupling is applied to the model to simulate the interaction between discrete and continuum phase. In order to track the continuum phase, a passive scalar equation is added to the solver. Eddy Life Time (ELT) is employed as a dispersion model. The influences of different non-dimensional parameters, such as Stokes number, Jet Reynolds number and mass loading ratio on the flow characteristics, are studied. The results of the simulations are verified with the available experimental data. It is revealed that more gravitational force is exerted on the jet as a result of the increase in mass loading, which deflects it more. Moreover, with an increase in the Reynolds number, the speed of the jet rises, and consequently, the gravitational force becomes less capable of deviating the jet. In addition, it is observed that by increasing the Stokes number, the particles leave the jet at higher speed, which causes a lower deviation of the jet.

1. Introduction

Particle-laden jet flows have various engineering applications, such as fuel injection in internal combustion engines. In this regard, the jets laden with particles have been investigated by many researchers experimentally and computationally [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. Many approaches have been used to experimentally analyze the jets laden with particles like Phase-Doppler Anemometry (PDA), Laser-Doppler Anemometry (LDA) and high-speed camera image processing. Fundamental experiments on dilute particle jets by the use of LDA can be found in [1,2,3,4]. Fan et al. investigated the effect of dispersed phase on the flow characteristics of a coaxial jet with particles laden. They showed when the particle ratio increases, the spreading ratio decreases [2]. Puttinger et al. analyzed the dilute particle-laden flow by using high-speed camera and image processing [5]. Due to the high spatial and temporal resolution of the latter method, it can be employed to study the turbulent flow with a high concentration ratio. This can be used for a wide variety of data and does not depend on the preset parameters. Generally, there are two main methods to simulate particle-laden flow, namely Eulerian–Eulerian and Eulerian–Lagrangian [6]. The former is used when the continuum assumption can be applied for both phases, while the latter is applied when particles constitute dispersed phases and fluid constitutes continuum phases. Both techniques have pros and cons [6]. A major problem of the first method is that some information such as the motion of the individual particles, and the collision between particles are not considered. In addition, applying boundary conditions for the second phase is challenging. However, the advantage of the Eulerian–Eulerian method is that the same equation and numerical methods are applied to both phases. In the Eulerian–Lagrangian method, more detailed information about the motion of particles and their interaction with the surrounding can be obtained. However, some computational difficulties arise when the number of particles increases or the size of particles decreases. For instance, when the size of particles decreases, the effect of turbulent eddies on particles should be calculated even in smaller sizes; therefore, the cost of the computational method will be increased. Wang et al. [7] simulated the particle-laden jet using the Eulerian–Lagrangian method. In their study, the gas-phase and fluid-phase velocities show good agreement with the experimental results. They investigated the effects of different parameters on the flow behavior, including particle size, Stokes number and mass loading ratio. The Stokes number ($Stk={\tau }_{p,st}/{\tau }_F$) is defined as the ratio of the characteristic time of a Lagrangian phase (${\tau }_{p,st}={\rho }_p{d}_p^2/18\mu$) to a characteristic time of the continuum phase (${\tau }_F=D_{jet}/U_{jet}$). The Eulerian–Eulerian method was applied to simulate particle-laden jet turbulent flow in Ref [8]. In their research, the major roles of the particle size and loading ratio on the flow characteristic of particle-laden jets are shown. The Eulerian–Lagrangian method was employed by Jebakumar and Abraham [9] in which the changes in centerline velocities and spreading rate were investigated by varying parameters such as fluctuating particle velocity at the jet inlet, gas-phase turbulence intensity. Lue et al. showed the major effect of Stokes numbers on the clustering of particles by means of one-way coupling direct numerical simulation of a turbulent particle-laden jet [10]. The effect of the Stokes number $\left(stk\right)$ on the particle motion and turbulent characteristics of particle-laden flow in the channel was investigated by Lee [11]. They found, for Stokes number of 5 that the particles increase the flow turbulence, while increasing the Stokes number of particles leads to decreasing the turbulence. The coupling behavior of solid particles and fluid makes modeling of this process complex. The coupling between solid and fluid phases has a substantial impact on multi-phase flows. If fluid exerts effects on particle motion but no reverse influence exists, the flow is called one-way coupling regime. If both phases have impacts on each other, the flow is called two-way coupling regime. In addition, if the collision between particles becomes important, the flow is called 4-way coupling regime. The Stokes number $\left(Stk\right)$ and concentration or loading ratios are important parameters in particle-laden flows to determine the coupling between solid and fluid phases. For more details, the reader is referred to reference [12]. Numerical simulations were carried out by Wang et al. [7] using one-way coupling in order to estimate the evolution of particles in a plane turbulent wall jet. Kazemi et al. [13] performed a CFD simulation with the assumption of two-way coupling between two phases. They applied three different dispersion models. Their results suggested that the particle loading ratio should be considered in choosing the dispersion model. The Eddy Life Time (ELT) model is considered as one of the conventional dispersion models which can be used to investigate particle dispersion resulted from instantaneous turbulent velocity fluctuations [14]. In the present study, two-way coupling simulation of horizontal turbulent jet laden with particles in Eulerian–Lagrangian framework is conducted employing Reynolds averaged Navier–Stokes (RANS) model and ELT model is employed for modeling particle dispersion. The particle trajectory resulting from the simulation is compared with the experimental results reported by Fleckhaus et al. [15].

The heavy jets, caused by the high difference in density between the jet and surrounding fluid, were investigated by some authors to study the mixing characteristics of wastewater discharge [16,17]. However, there are many industrial applications that contain the discharging process of heavy particles and fluid, which form heavy jets. For example, in industrial wastewater and an exhaust gas flow, which in the effect of solid particles on the jet profile of disposal, wastes must be taken into consideration. In this paper, the effects of three non-dimensional parameters are investigated on the characteristics of heavy jets, which result from heavy particles. The influence of Stokes number, mass loading ratio and Reynolds number (Re = ${\mathrm{U}}_{\mathrm{jet}}{\mathrm{D}}_{\mathrm{jet}}/\mathsf{\upsilon }$), where $\mathsf{\upsilon }$ is the Kinematic viscosity of the continuum phase, and mass loading ratio, which is defined mass fraction of the dispersed phase to continuum phase of the jet on the trajectory distribution of jet with particle-laden are investigated. The definition of the Stokes number is based on the Stokes’ relaxation time scale, and this number represents the diameter of the particles with a density of 2590kg/m³. In addition, turbulence modulation is not considered in this work. It shows that the effect of gravity on the motion of the jet can be seen clearly by the trajectory of particles.

2. Governing Equation

In this study, the carrier flow is considered as a continuous phase and particles as dispersed phases, which are simulated using the Eulerian–Lagrangian framework. The two-way coupling is applied to study the interaction between two phases.

2.1. The Eulerian Equations

The gas phase is described by the Reynolds averaged Navier–Stokes equations (RANS) for incompressible flow. The continuity and momentum equations can be written as follows [18]:

$$\frac{\partial u_j}{\partial x_j}=0$$

(1)

$$\frac{\partial u_i}{\partial t}+\frac{\partial u_i{u}_j}{\partial x_j}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\frac{\partial {\tau }_{ij}}{x_j}+\frac{1}{\rho }{S}_{u_i}$$

(2)

with $\rho$ being the density of the gas, $P$ the pressure, $u_i$ the velocity in the $i^{th}$ direction; $S_{u_i}$ defines momentum source term because of the interaction between Eulerian and Lagrangian parts; ${\tau }_{ij}$ denotes the effective stress tensor which can be defined as below:

$${\tau }_{ij}=\left(\nu +{\nu }_t\right)\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right]$$

(3)

where $\nu$ and ${\nu }_t$ are kinematic viscosity and turbulent kinematic viscosity, respectively; $\delta$ is the Kronecker delta. Here the turbulent viscosity ${\nu }_t$ is computed using the k-ω SST (Shear Stress Transport) turbulence model [19]. In the present study, a passive scalar equation is added to the solver as a tracer to track the gas phase, to investigate different fluid dynamic behavior between the tracer (particles with zero inertia) and inertial particles, which can be described as follows:

$$\frac{\partial c}{\partial t}+\frac{\partial u_{j}c}{\partial x_j}=\left(\frac{\nu }{\mathrm{Sc}}+\frac{{\nu }_t}{{\mathrm{Sc}}_t}\right)\frac{{\partial }^{2}c}{\partial x_j\partial x_j}$$

(4)

In the above equation, Schmidt number $\left(\mathrm{Sc}\right)$ and turbulent Schmidt number $\left({\mathrm{Sc}}_{\mathrm{t}}\right)$ are 0.4 and 0.7, respectively.

2.2. The Lagrangian Equations

The dispersed phase is composed of parcels. Each parcel is defined by the group of properties including velocity $\left(u\right)$, diameter $\left(d\right)$, density $\left(\rho \right)$, position $\left(x\right)$, computational weight (ω). The diameter, density and computational weight of parcels are constant after injection. It is assumed that the drag and gravity forces have a dominant effect compared to other forces. Therefore, other forces are neglected. In the case of the heavy particle, i.e., ρ_{2 >>} ρ₁, the volume force applied by surrounding fluid flow can be assumed to be equal to the drag force [20].

The governing equations for the disperse phase can be written as [12,21]:

$$\frac{d{x}_{p,i}}{dt}=u_{p,i}$$

(5)

$$m_p\frac{d{u}_{p,i}}{dt}=F_{D,i}+m_p{g}_i$$

(6)

The drag force $(F_{D,i})$ acts on a particle $i$ which is equal to:

$$F_{D,i}=m_p\frac{u_{r,i}}{{\tau }_p}$$

(7)

where relative velocity $\left(u_{r,i}\right)$ and particle relaxation time $({\tau }_p)$ can be described as:

$$u_{r,i}=u_{s,i}-u_{p,i}$$

(8)

$${\tau }_p=\frac{d_p^2{\rho }_p}{18\mu f_d}$$

(9)

Relaxation time is applied to show the resistance of a dispersed phase to adjust to continuum phase motion and $f_d$ is the correction factor, which is defined as below:

$$f_d=\{\begin{array}{l}1+0.15R{e}_p^{0.687};R{e}_p<1000\\ 0.44R{e}_p/24;otherwise\end{array}$$

(10)

$R{e}_p$ is the particle Reynolds number base on the relative velocity $(u_r)$ which is defined as:

$$R{e}_p=\frac{\rho u_r{d}_p}{\mu }$$

(11)

where $u_{s,i}$ denotes the instantaneous gas velocity. It is the summation of mean gas flow velocity and fluctuating gas flow velocity. The mean velocity is interpolated at the parcel location from RANS equations; the eddy lifetime $\left(\mathrm{ELT}\right)$ dispersion model is applied to calculate the fluctuating velocity. In the ELT model, the fluctuating velocity $\left(u_{s,i}^{″}\right)$ is defined as [13]:

$$u_{s,i}^{″}={\xi }_j\sqrt{\frac{2k}{3}}$$

(12)

where ${\xi }_j$ is the normally distributed random number of the standardized Gaussian distribution and $k$ denotes the turbulent kinetic energy. The updated value fluctuating velocity should be calculated after the characteristic lifetime $({\tau }_e)$:

$${\tau }_e=-{\tau }_L\ln \left(r\right)$$

(13)

where $r$ is a random number between 0 and 1 and Lagrangian integral time scale $({\tau }_L)$ calculated as:

$${\tau }_L=C_1\frac{k}{\epsilon }$$

(14)

where $C_1$ is an empirical factor whose value is adjusted based on the experimental results. In many studies, the value for this parameter is chosen as 0.3, which is also applied in the work of Refs [22,23].

The momentum transfer term $S_{u_i}$ in Equation (2) calculated for $l^{th}$ cell in each time step as:

$${\left[S_{u_i}\right]}^l=\frac{1}{\Delta t{\forall }_C^{\left[l\right]}}{\displaystyle \sum }_{i=1}^{N_p}\left(G_l\left(x_{p,i}\right){\omega }_i{F}_{D,i}\right)$$

(15)

where $\Delta t$ is the time step, ${\forall }_{\mathrm{C}}^{\left[l\right]}$ is the volume of the $l^{th}$ cell, $G_l\left(x\right)$ is the Kernel function of the $l^{th}$ cell, ${\omega }_i$ is the statistic weight of the parcel and $F_{D,i}$ is the drag force applied on the ith parcel. The momentum transfer term is obtained by the summation of product term for all parcels. The Kernel function is considered as a linear weighted interpolation operator in OpenFoam [24].

3. Analysis and Modeling

To solve the governing equations, open-source OpenFOAM package (www.openfoam.org) is employed. In OpenFOAM, the equations of the continuum phase are solved with the discrete cell-centered finite volume method on the discrete computational grid. The boundary condition at the inlet is a fixed velocity profile, zero gradient at the outlet and no-slip conditions are applied for the walls. The pressure is zero gradient for the inlet, and fixed value (0) applied to the outlet. The Concentration field is of a fixed value for the inlet and zero gradients for all others. The turbulent intensity is 10% for the inlet, and standard wall function k-$\omega$ are used for walls. The parcels are injected at a rate of 1,000,000 per second. Each parcel will be considered omitted after reaching the outlet, and in approaching the wall, they are considered rebounded.

To consider grid independency for our domain, three different grid sizes (1,100,000, 700,000, 400,000) were compared based on the velocity profile of the continuum phase along three sections $Z/D=10,20$ and $30$. The maximum error was 4% for the 700,000-cell grid. Therefore, the simulations were done using 700,000 cells in the domain. The computational domain and its grid quality are shown in Figure 1. A block of the structured grid is applied between the jet and cylindrical enclosure instead of using the polar grid in order to keep the number of the grid elements as low as possible. The PISO-SIMPLE (PIMPLE) algorithm with exterior loops is used to solve the pressure-velocity coupling of the Navier–Stokes equations and increase the time intervals [25]. The time interval is determined dynamically based on the Courant number. This number is kept less than 1 for the entire simulation. Diffusion terms are discretized following the standard second-order central differencing scheme, and the convection terms are discretized by the second-order “Gauss limited Linear” scheme, which based on Total Variation Diminishing (TVD) method. Lagrangian equations are integrated by analytical methods [26]. Particle tracking is performed using the Face-to-face particle tracking procedure approach [27].

In the present study, different cases are considered (summarized in

Table 1. The parameters of four different cases.

Mass Loading	$\mathbf{Stokes}\mathbf{Number}$$\left(\mathbf{Stk}\right)$	$\mathbf{Reynolds}\mathbf{Number}\left(\mathbf{Re}\right)$	Computational Weight of Parcels
0.53	5	5000	40
1.06	5	5000	80
2.13	5	5000	160
2.13	5	10,000	160
2.13	5	20,000	160
2.13	1	5000	17, 20
2.13	10	5000	55

) to investigate the effects of two-way coupling on the motion of a particle-laden jet. Since in all three mass loadings, and three investigated sizes of the particles, the same number of the parcels was injected, different statistical weights were obtained for the particles. The following table presents the statistical weight of the particles. Stokes number ($Stk$) is defined as the characteristic time of Lagrangian (${\tau }_{p,st}$) to the characteristic time of the flow (${\tau }_F$), where each is as the following.

$$Stk=\frac{{\tau }_{p,st}}{{\tau }_F}$$

(16)

$${\tau }_{p,st}=\frac{{\rho }_p{d}_p^2}{18\mu }$$

(17)

$${\tau }_F=\frac{D_{jet}}{U_{jet}}$$

(18)

The computational weight of parcels in each simulation is calculated as below:

$$\omega =\frac{{\dot{N}}_p}{{\dot{N}}_c}$$

(19)

where ${\dot{N}}_p$ and ${\dot{N}}_c$ are the number of particles and computational parcels respectively that inject per second during the simulation. Since ${\dot{N}}_c$ is one million per second, we have two million parcels after two seconds in our simulations. ${\dot{N}}_p$ is calculated from the injected mass flow rate (${\dot{m}}_{inj}$) as below:

$${\dot{N}}_p=\frac{{\dot{m}}_{inj}}{m_p}$$

(20)

where ${\mathrm{m}}_{\mathrm{p}}$ is the particle mass.

4. Validation

The validation of the simulation results was conducted by comparing the experimental normalized axial velocity observations reported by Fleckhaus et al. [15]. They used glass beads with properties of ${\rho }_p=2590kg/m^3,$ $d_p=64\mu m$, $Re=$ $20,000$ and with a mass loading ratio of 0.3. Figure 2 presents the predicted velocity of the disperse phase at three locations ($Z/D=10,20$ and $30$) along the centerline. The resulting values compare well with the published experimental data in [15]. The results were normalized by the maximum velocity of continuum phases along section $Z/D=10$. In Figure 2, the results of simulations are compared with the experimental data from [15], which are in a very good agreement. In Figure 2, dispersed and continuum phase are indicated using dashed and solid lines, respectively. Two million parcels are injected to reduce the statics errors to an acceptable value to obtain the graphs of the disperse phase for each case. Nevertheless, perturbations can be seen in the diagrams. These high-frequency oscillations are omitted using a high pass filter after applying Fourier transform in the frequency domain. There are different approaches such as filtering, parcel-number-density control algorithms which are used to achieve numerically convergent solutions in Eulerian–Lagrangian (EL) simulations [21,28,29]. As it is shown in Figure 2, the decay rate obtained from the simulation results for the continuum phase is higher than those obtained from the experimental data. It can be seen that at larger distances from the nozzle, the normalized velocity of the dispersed phase falls below that which is presented by the experimental data. In the case of the Lagrangian phase, the results show overestimated value for section Z⁄D = 10. However, the simulation results at Z⁄D = 20, 30 shows lower value comparing experimental result, due to a higher decay rate. The numerical dissipation can be a reason for the high decay rate of velocity profiles.

5. Results and Discussion

In the present study, to investigate the influence of the presence of particles on the jet motion direction, a passive scalar $\left(c\right)$ as a tracer of concentration of the continuum phase is injected at the nozzle (jet) entrance. This scalar shows the concentration of a scalar tracer which comes out from the jet. Figure 3, from left to right, demonstrates the influence of mass loading, Reynolds number of the jet and Stokes number of the particles on the concentration distribution of the continuous phase and distribution of the axial velocity of the continuum phase, respectively. The jet velocity profile at the top of the jet is quite similar to that of the spread of the jet tracer. In comparison, at the bottom of the jet, where the particles leave this region of the jet, the rate of spread of the jet tracer is higher. The amplitude of the velocity shows a general form of the jet. However, it cannot present the distribution of the flow rather the extension of the momentum in the calculation domain. In this study, one aim is to determine the influence of the particles on the mass motion of the continuous phase. The contour of the velocity is shown as line profiles in Figure 3. Velocity, unlike concentration, does not accumulate at one point. This is due to the presence of convection inherent in locations where there is flow velocity. Moreover, velocity does not accumulate over time at any given point in the domain as a result of the dissipation. Using the contour of concentration, it has been shown that the released flow is concentrated at the bottom of the main trajectory line of the jet. This is of critical importance in waste disposition, where it may contain particles that should not be accumulated at one point.

This lies in the fact that the downward motion of the particles changes the momentum direction downward and deviates from the original path. Moreover, it counteracts the influence of the ambient entrainment toward the jet in the lower region of the jet, which increases the rate of spread of the jet tracer as the particles leave the jet. In Figure 4, the path of the motion of the continuum phase is drawn by connecting the maximum points of the tracer. Furthermore, the rate of spread of the jet tracer at the bottom of the jet increases with growth in mass loading. Initially, the particles and jet have the same velocity. However, after a while, the particles tend to push jet in the downward direction as a result of gravitational force. It can be said that, gradually, the particles have more effects on the jet streamline. At higher Stokes number, the accumulation of continuum phase occurs at larger distance from the nozzle. This is due to the presence of larger particles, having a lower area/volume ratio, resulting in smaller drag forces per volume. It can be said that particles with a larger diameter and higher stokes number will decay their momentum at larger distances from the nozzle. On the other hand, particles will exit the jet and transfer their downward velocity, as a result of gravity, to the surrounding fluid. Due to this, smaller particles can transfer more downward momentum to the fluid. As a result, jets laden with smaller particles will have a downward trend. Initially (Z/D = 0), the particles and jet have the same velocity. However, after a while, the particles tend to push jet in the downward direction as a result of gravitational force. It can be said that, gradually, the particles have more effects on the jet streamline. The decay rate of the axial velocity and the concentration of the continuum phase (tracer) in the direction of the jet (Figure 4) are depicted in Figure 5.

The results reveal that the impact of mass loading on the gas-phase concentration decay rate is considerably low. However, the decay rate of the axial velocity decreases with a rise in mass loading. The total momentum of the particle laden-jet becomes larger, with an increase in mass loading. Meanwhile, the velocity of the gas-phase, which is deteriorated and scattered, is compensated by the transfer of the momentum from particles to the gas.

Figure 3b and Figure 6 indicate that the deflection of the jet path is diminished by increasing the Reynolds number. This is due to an increase in jet momentum, while the gravitational force is not strong enough to deflect the jet trajectory. Moreover, the increase in jet velocity enhances the ambient entrainment, and the initiated radial velocity toward the jet hinders the release of the particles from the jet, which leads to a reduction in the spread rate. However, by decreasing the Reynolds number, the gravitational force can deviate the particles from the jet and shift the momentum to the ambient, giving rise to the decay rate of the axial velocity.

Figure 3c shows that the spread of the jet tracer suddenly increases at the bottom of the jet. Considering Figure 7, this increase is due to the particles being released from the jet bottom. These particles drag down part of the continuous phase, which increases the depth of the continuous phase penetration beneath the jet. Moreover, it can be seen that particles move out of the jet faster and farther away due to the Stokes number rise. This leads to an increase in the depth of the continuous phase beneath the jet. Nevertheless, since the coupling amongst particles is weaker for higher Stokes number, they cannot deflect the jet trajectory and gradually depart the jet. In comparison, as Figure 8 demonstrates, particles with stronger coupling, due to lower Stokes number, remain in the jet resulting in jet deflection.

6. Conclusions

The motion of a particle-laden jet was simulated with the particle size of $d_p=64\mu m$ and density of ${\rho }_p=2590kg/m^3$. The direction of the motion of the jet was traced by adding a tracer to the solver. The simulation results were validated with available experimental data. Moreover, the impacts of Stokes numbers, mass loading ratio and Reynolds number were studied. The following statements can be inferred from the results:

The total momentum of the particle-laden jet becomes larger with an increase in mass loading and the velocity of the gas-phase, which is deteriorated and scattered, is compensated by the transfer of the momentum from particles to the gas.

With a decrease in Reynolds number, the gravity force results in deflection of the particles from the jet and shifts the momentum to the ambient, giving rise to the decay rate of the axial velocity.

As the Stokes number increases, the particles leave the jet faster and are therefore less capable of bending the continuous phase trajectory of the jet. However, particles with lower Stokes numbers have less slip velocity than the continuous phase and their longer presence in the jet results in greater gravitational force effects and consequently more jet deflection.