Dispersion of Heavy Particles in a Turbulent Channel Flow

These parameters yield a value of Re = vah/ν = 2100 for the Reynolds number of the flow. The authors predict a shear velocity of uτ = 0.11775 m/s for a shear Reynolds number of approximately Reτ = uτh/ν = 150.

Turbulent Flow

This example solves the Reynolds-averaged Navier-Stokes (RANS) equations and uses the standard k-ε turbulence model (Ref. 2), one of the most frequently used turbulence models in computational fluid dynamics. The k-ε model introduces dependent variables and transport equations for two new quantities:

•

The turbulent kinetic energy k (SI unit: m2/s2) represents the energy per unit mass associated with eddies in the flow.

•

The turbulent dissipation rate ε (SI unit: m2/s3) indicates the rate at which turbulent kinetic energy in the eddies is converted to thermal energy.

These new dependent variables provide insight into the size and lifetime of eddies in the flow:

•

The ratio k/ε has units of time. The average eddy lifetime is often estimated by multiplying this ratio by a dimensionless constant of order unity.

•

The ratio k3/2/ε has units of length and indicates the length scale of the largest eddies in the flow.

In addition to being a very common turbulence model in its own right, the k-ε model is relatively easy to combine with Lagrangian particle tracking methods because the turbulence variables immediately give estimates of the amplitude of the velocity perturbations due to turbulent eddies (proportional to

) and the average eddy lifetime (proportional to k/ε). For a stationary, incompressible flow, the transport equations solved when using the k-ε turbulence model are

where the dependent variables are the fluid velocity u (SI unit: m/s), pressure p (SI unit: Pa), and the aforementioned transport variables k and ε.

The dimensionless constants in these equations have the following default values:

For more details on the different turbulence models that are available, see the Single Phase Flow Interfaces chapter in the CFD Module User’s Guide. The k-ε model is usually a convenient choice when coupling the turbulent flow to a particle tracing simulation because of its relatively low memory requirements and fast convergence, and because it directly allocates degrees of freedom for k and ε, which are necessary to accurately model turbulent dispersion of the particles.

Lagrangian Particle Tracking

Marchioli et al. (Ref. 1) provide a wide range of DNS results for simulations with and without the lift and gravity forces. In this example, only the drag force is considered.

Because the particle density is several orders of magnitude greater than the air density, the buoyant force and added-mass effect can safely be neglected. The particles are also large enough that the Brownian Force can be ignored.

The equation of motion for each particle is therefore

where

•

mp (SI unit: kg) is the particle mass,

•

q (SI unit: m) is the particle position, and

•

FD (SI unit: N) is the drag force.

In general, the drag force is defined as

(1)

where u (SI unit: m/s) is the fluid velocity at the particle’s position and v (SI unit: m/s) is the particle velocity. The particle relaxation time τp (SI unit: s) is defined as

(2)

where

•

ρp (SI unit: kg/m3) is the particle density,

•

dp (SI unit: m) is the particle diameter, and

•

CD (dimensionless) is the drag coefficient.

The choice of drag law, which determines the appropriate definition of CD, depends on the relative Reynolds number Rer (dimensionless) of the particle in the fluid. For a spherical particle the relative Reynolds number is

For

the Stokes drag law is applicable,

so that Equation 2 reduces to

Particle Tracing with Turbulent Dispersion

In the Reynolds-Averaged Navier-Stokes (RANS) approach, the turbulent eddies are only solved for in a statistical sense, by estimating their energy and dissipation rate. A Direct Numerical Simulation (DNS) could resolve the individual eddies but is often too memory-intensive and time-consuming for many practical applications. Therefore, in Equation 1 the fluid velocity u is not given deterministically. Instead, the fluid velocity is treated as a linear combination of a deterministic part based on the mean flow (which is the field u solved for by the RANS equations) and a turbulent perturbation term Δu,

The amplitude and direction of the velocity perturbation Δu is derived from the turbulence variables k and ε that are solved for by the k-ε turbulence model.

The Particle Tracing for Fluid Flow interface supports two different formulations for the turbulent dispersion term:

•

The Discrete Random Walk (DRW) model is similar to the modified Eddy Interaction Model of Gosman and Ioannides (Ref. 3). In the DRW model, unique velocity perturbations are sampled and added to the mean velocity field at discrete times based on the estimated eddy lifetime.

•

The Continuous Random Walk (CRW) model, or Continuous Filter White Noise (CFWN) model (Ref. 4). In the CRW model, the velocity perturbations are integrated over time. The eddy velocity and lifetime are built into the time derivatives of the evolving velocity perturbation components.

In this example, the CRW model will be used. A brief overview of this model is given below. For more comprehensive details on both turbulent dispersion models, see the Particle Tracing for Fluid Flow chapter of the Particle Tracing Module User’s Guide.

The classical Langevin-equation model for homogeneous isotropic stationary turbulence (HIST) is (Ref. 5)

(3)

where the subscript

indicates a component of the fluid velocity field. In isotropic turbulence, the RMS fluid velocity fluctuation in any direction is equal,

Because the turbulence is isotropic, the velocity perturbations can be aligned with any orthonormal coordinate system. In the Particle Tracing for Fluid Flow interface they are simply aligned with the Cartesian coordinates.

Inhomogeneous and Anisotropic Turbulence

In wall-bounded flows the assumption of homogeneous isotropic turbulence no longer applies. The turbulence becomes inhomogeneous because the turbulent kinetic energy is heavily damped close to the walls, giving it a nonzero gradient in this region. The turbulence is anisotropic because the eddies in the region close to the wall are not equally likely to point in any direction; the instantaneous velocity component normal to the wall is more heavily damped than those in the streamwise and spanwise directions.

In the Particle Tracing for Fluid Flow interface, anisotropic and inhomogeneous turbulent dispersion are modeled by selecting the check box in the settings for the node. This check box is only available when the Continuous Random Walk model for turbulent dispersion is used. When this check box is selected, corrections for inhomogeneous and anisotropic turbulence are applied in the region y+ < 100, where y+ (dimensionless) is the wall distance in viscous units,

where x2 (SI unit: m) is the normal distance to the nearest wall, ν (SI unit: m2/s) is the kinematic viscosity of the fluid, and uτ (SI unit: m/s) is the friction velocity,

where τw (SI unit: N/m2) is the wall shear-stress. The friction velocity is usually defined on the boundaries by one of the turbulent flow interfaces. In the region of inhomogeneous, anisotropic turbulence, the components of the turbulent velocity perturbation are defined in the coordinate system given by

The normalized Langevin equations in these directions are

To account for the anisotropy of the flow, the following definitions of the σi terms are given. These expressions are taken from Dehbi (Ref. 5), wherein they are attributed to DNS fits of channel flow as described by Dreeben and Pope (Ref. 6).

In the boundary layer, the Lagrangian time scale τi from Equation 3 is approximately the same in all directions:

Kallio and Reeks Ref. 7 give the following polynomial fit for the time scale:

where ν (SI unit: m2/s) is the fluid kinematic viscosity. Away from the wall, the Lagrangian time scale is simply

where CL is a dimensionless constant.

Effect of Particle Inertia

To investigate the effect of particle inertia on anisotropic turbulent dispersion in the channel, a is performed over the Stokes number St (dimensionless) of the particle:

where τp (SI unit: s) is the particle relaxation time,

Marchioli et al. (Ref. 1) use six values of the Stokes number: 0.2, 1, 5, 15, 25, and 125. For a given value of St, the corresponding particle diameter is

Results and Discussion

The results of the turbulent flow simulation are shown in Figure 1 and Figure 2. The channel has an extremely high aspect ratio, so to better visualize the geometry, has been selected from the list in the settings. This allows the coordinate axes to be scaled independently of each other so that the plot fits the Graphics window.

Figure 1 shows the fluid velocity magnitude and velocity streamlines. Figure 2 shows the ratio of turbulent kinetic energy to turbulent dissipation rate in the modeling domain. This is greatest near the left boundary, which is the axis of symmetry. It decreases near the wall, suggesting that eddy lifetimes are much shorter there. It is useful to plot this ratio before proceeding with the particle tracing simulation because it shows the minimum resolution in time needed to accurately capture the particle-eddy interactions.

The particle trajectories for the greatest simulated value of the Stokes number, St = 125, are shown in Figure 3. The particles were released uniformly along the cross section halfway down the channel, to prevent them from hitting any or boundaries as a result of the turbulent diffusion. For intermediate values of the Stokes number, the particles cluster near the wall, as shown by the histograms in Figure 4. At very low Stokes number, St = 0.2, the effect of anisotropic turbulence on the number density of particles is less pronounced because such particles don’t have enough inertia to frequently cross the eddies. At very high Stokes number, St = 125, the effect is less pronounced compared to intermediate values because the inertia of such particles is so high that they often reflect off the wall and back into the bulk medium.

Figure 1: Fluid velocity in the channel. Velocity streamlines are shown as white lines.

Figure 2: Turbulence time scale, or ratio of turbulent kinetic energy to turbulent dissipation rate, in the channel.

Figure 3: Particle trajectories in the channel.

Figure 4: Comparison of histograms of particle position for different values of the Stokes number.

References

1. C. Marchioli, M. Picciotto, and A. Soldati, “Influence of gravity and lift on particle velocity statistics and transfer rates in turbulent vertical channel flow”, International Journal of Multiphase Flow, vol. 33, no. 3, pp. 227–251, 2007.

2. D.C. Wilcox, Turbulence Modeling for CFD, 2nd ed., DCW Industries, 1998.

3. A. D. Gosman and E. Ioannides, “Aspects of computer simulation of liquid-fueled combustors”, Journal of Energy, vol. 7, no. 6, pp. 482–490, 1983.

4. L. Tian and G. Ahmadi, “Particle deposition in turbulent duct flows–comparisons of different model predictions”, Aerosol Science, vol. 38, 2007, pp. 377–397.

5. A. Dehbi, “Turbulent particle dispersion in arbitrary wall-bounded geometries: A coupled CFD-Langevin-equation based approach”, International Journal of Multiphase Flow, vol. 34, no. 9, pp. 819–828, 2008.

6. T. D. Dreeben and S. B. Pope, “Probability density function and Reynolds-stress modeling of near-wall turbulent flows”, Physics of Fluids, vol. 9, no. 1, pp. 154–163,1997.

7. G. A. Kallio and M. W. Reeks, “A numerical simulation of particle deposition in turbulent boundary layers”, International Journal of Multiphase Flow, vol. 15, no. 3, pp. 433–446, 1989.

Particle_Tracing_Module/Fluid_Flow/flow_channel_turbulent_dispersion

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select ε .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

Load the model parameters from a file.

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file flow_channel_turbulent_dispersion_parameters.txt.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type halfWidth.

In the text field, type height.

Click

The geometry has a very high aspect ratio. Adjust the View settings to make it easier to see.

Definitions

In the window, expand the node.

Axis

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Click

Click the

button in the toolbar.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rhof	kg/m³	Basic
Dynamic viscosity	mu	muf	Pa·s	Basic

Turbulent Flow, k-ε (spf)

Inlet 1

In the window, under right-click ε and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 2 only.

Locate the section. In the Uav text field, type va.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 3 only.

Symmetry 1

In the toolbar, click

and choose .

Select Boundary 1 only.

Mesh 1

Mapped 1

In the toolbar, click

Distribution 1

Right-click and choose .

Select Boundaries 1 and 4 only.

In the window for , locate the section.

In the text field, type nelemHeight.

Distribution 2

In the window, right-click and choose .

Select Boundaries 2 and 3 only.

In the window for , locate the section.

From the list, choose .

In the text field, type nelemWidth.

In the text field, type 15.

Click

Study 1

In the toolbar, click

Results

Velocity (spf)

Add some streamlines to the default plot of the fluid velocity.

Streamline 1

Right-click and choose .

Select Boundary 2 only.

In the window for , locate the section.

Find the subsection. From the list, choose .

In the toolbar, click

Compare the resulting plot to Figure 1.

Turbulence Time Scale

In the toolbar, click

and choose .

In the window for , type Turbulence Time Scale in the text field.

Surface 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose ε .

In the toolbar, click

Compare the resulting plot to Figure 2.

Next, solve for the particle trajectories in the turbulent flow.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

At the bottom of the section, clear the check box next to . The particle trajectories are not solved for in the study step.

Click in the window toolbar.

In the toolbar, click

to close the window.

Add Study

In the toolbar, click

to open the window.

Go to the window.

At the bottom of the section, clear the check box next to the Turbulent Flow, k-ε interface, which will not be solved for in the study step.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Particle Tracing for Fluid Flow (fpt)

Particle Properties 1

In the window, under click .

In the window for , locate the section.

From the ρp list, choose . In the associated text field, type rhop.

In the dp text field, type dp.

Assign the condition to both the top and bottom boundaries. Assign the condition at the symmetry axis. The distance from the remaining boundary will be used to compute the anisotropic turbulent velocity perturbations.

Outlet 1

In the toolbar, click

and choose .

Select Boundaries 2 and 3 only.

Symmetry 1

In the toolbar, click

and choose .

Select Boundary 1 only.

Wall 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Add the drag force, using the fluid velocity and turbulence variables computed in the previous study.

Drag Force 1

In the toolbar, click

and choose .

Select Domain 1 only.

In the window for , locate the section.

From the u list, choose .

Locate the section. From the list, choose .

From the k list, choose .

From the ε list, choose .

Select the check box.

In the u∗ text field, type ustar_exp.

Release particles from the middle of the channel. Initially, the number density of particles is uniform over the width of the channel.

Release from Grid 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog box, choose from the list.

In the text field, type halfWidth/(2*Np).

In the text field, type halfWidth*(1-1/(2*Np)).

In the text field, type Np.

Click .

In the window for , locate the section.

In the qy,0 text field, type height/2.

Locate the section. Specify the v0 vector as


u	x
v	y

Study 2

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
St (Particle Stokes number)	0.2 1 5 15 25 125

Step 1: Time Dependent

In the window, click .

In the window for , locate the section.

In the text field, type 0 t1 t2.

Click to expand the section. Find the subsection. From the list, choose .

From the list, choose .

In the toolbar, click

Results

Particle Trajectories (fpt)

Compare the default trajectory plot to Figure 3.

Number Density, St = 0.2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Number Density, St = 0.2.

Locate the section. From the list, choose .

In the list, select .

From the list, choose .

In the list, choose and .

Click to expand the section. From the list, choose .

In the text area, type St = 0.2.

Locate the section. Select the check box.

In the associated text field, type y<sup>+</sup>.

Histogram 1

Right-click and choose .

In the window for , locate the section.

In the text field, type fpt.df1.yplus.

Locate the section. In the text field, type 50.

Click to expand the section. Select the check box.

From the list, choose .

In the table, enter the following settings:


Legends
t<sub>1</sub>
t<sub>2</sub>

In the toolbar, click

Number Density, St = 0.2

The histogram shows that the number density throughout the cross section is nearly uniform. Use a logarithmic scale and manual axis limits to more easily compare such histograms for each value of the Stokes number.

Click the

button in the toolbar.

Click the

button in the toolbar.

In the window, click .

In the window for , locate the section.

Select the check box.

In the text field, type 1.4.

In the text field, type 150.

In the text field, type 50.

In the text field, type 1e4.

In the toolbar, click

Duplicate the 1D Plot Group containing the plot and select other values of the Stokes number to observe how particle inertia affects the particle number density in the channel cross section. As the Stokes number increases, the particles begin to accumulate in the boundary layer close to the wall. All six sets of histograms are shown in Figure 4.