# Stochastic Lagrangian wall deposition model for RANS prediction of deposition in turbulent gas-solid flows

Cheikhna Talebmoustaph, Pascal Fede, Olivier Simonin, Maxime Pallud, Priyank Maheshwari

► To cite this version:

Cheikhna Talebmoustaph, Pascal Fede, Olivier Simonin, Maxime Pallud, Priyank Maheshwari. Stochastic Lagrangian wall deposition model for RANS prediction of deposition in turbulent gas-solid flows. International Journal of Multiphase Flow, 2024, 178, pp.104900. 10.1016/j.ijmultiphaseflow.2024.104900. hal-04660037

HAL Id: hal-04660037

https://ut3-toulouseinp.hal.science/hal-04660037v1

Submitted on 23 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stochastic Lagrangian wall deposition model for RANS prediction of deposition in turbulent gas-solid flows

Cheikhna Talebmoustaph<sup>a,b</sup>, Pascal Fede<sup>b,*</sup>, Olivier Simonin<sup>b</sup>, Maxime Pallud<sup>a</sup>, Priyank Maheshwari<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Wall-bounded flows

Particle deposition

Particle-fluid interactions

RANS

## ABSTRACT

In industrial-scale Reynolds Averaged Navier-Stokes (RANS) simulations, a wall-modeled approach is often employed, specifically choosing a dimensionless wall-cell height much greater than one. However, in gas-solid flows, a coarse description of the boundary layer may lead to inaccurate results in terms of particle dispersion and deposition. In this study, we conducted an investigation on stochastic RANS dispersion models with deposition on smooth walls using a wall-modeled approach. These models were applied to turbulent channel flows in both vertical and horizontal orientations. The modeling approach involved the adoption of the Euler-Lagrange formalism and the RANS $k-\epsilon$ turbulence model in OpenFOAM®. To assess particle behavior first, a simple homogeneous isotropic stationary turbulence (HIST) case was examined. Fluid elements were also tracked in an inhomogeneous turbulence for dispersion analysis. These tests gave the good formulation for the dispersion model to be used for particle tracking coupled with turbulent channel flow simulation. However, the deposition results in turbulent channel using such stochastic dispersion model in the wall-cell were found to be consistently overestimated across various scenarios for low-inertia particles. To address these limitations arising from the coarse description of the boundary layer, we examined, implemented, and evaluated a Lagrangian stochastic wall deposition model. The activated deposition model demonstrated good agreement with experimental data.

## 1. Introduction

Transport of solid/liquid particles in turbulent flows is present in a wide range of environmental, industrial, and medical applications. Understanding the mechanisms behind particle-laden flows is crucial for understanding various phenomena, such as the atmospheric dispersion of pollutants, filtration, fooling/erosion of turbine blades, deposition on solar panels, etc. Thus, a thorough understanding of particle transport and deposition in turbulent flow fields is important for the successful development and application of reliable computational models in addressing these complex problems.

In wall-bounded flows, turbulence is anisotropic and heterogeneous in the near wall region (Dreeben and Pope, 1997), complex particle-fluid interactions are at play (Arcen et al., 2005). How accurate the boundary layer is captured by RANS simulations is crucial for particle dispersion and deposition. In practical CFD simulations for industrial applications, a wall-modeled approach (wall-cell dimensionless height $y_h^+ \gg 1$) is often adopted to reduce the computational cost. As a result, accurate particle dispersion and deposition becomes challenging

as particles see a coarse description of the boundary layer. In the context of particle-laden flows using RANS turbulence models, particles only see average quantities of the fluid phase. In order to account for particle dispersion, models must be employed. There is a variety of ways to model turbulence effects on particles (Huilier, 2021). How the fluid velocity “seen” or “experienced” $\mathbf{u}_f$ along the trajectory can be simulated is the crucial question in all dispersion models.

Stochastic models are often used to consider the missing fluctuating turbulent fluid velocity components. These quantities are obtained at the particle position from the Eulerian field given by the RANS results such as the turbulent kinetic energy, $k$, and dissipation rate, $\epsilon$. They also rely on the different time, velocity and length scales used. Two types of stochastic models are widely used to model particle dispersion: Eddy Interaction Models (EIM) and Langevin models.

Eddy Interaction Models have been popular due to their relative ease of implementation and reasonable computational expense. It was first introduced by Hutchinson et al. (1971) and extended by Gosman and Ioannides (1981). Kallio and Reeks (1989) employed such model

<footnote><sup>a</sup> TotalEnergies OneTech, 2 place Jean Milier, 92400, Courbevoie, France</footnote>

<footnote><sup>b</sup> Institut de Mécanique des Fluides de Toulouse (IMFT), Allée du Professeur Camille Soula, 31400, Toulouse, France</footnote>

<footnote>* Corresponding author.
E-mail address: pascal.fede@imft.fr (P. Fede).</footnote>
along with a Root Mean Square (RMS) wall-normal velocity and a Lagrangian timescale fitted to simulate particle deposition in a pipe flow. Matida et al. (2000) studied the dependence of different fits in the near wall region and additional forces. Parker et al. (2008) used an EIM coupled with RANS two-equation models ($k-\epsilon$) and Reynolds-Stress Model. The authors concluded that even when the boundary layer is fully resolved, namely for $y_h^+ \sim 1$ with $y_h^+ = y_h u_\tau / v_f$ and $y_h$ the height of the first cell, $u_\tau$ the wall friction and $v_f$ the kinematic fluid viscosity, deposition is overpredicted for low-inertia particles using two-equation models. Aguinaga et al. (2009) proposed a Lagrangian stochastic wall deposition model coupled with an EIM in the bulk flow to account for the coarse description of the boundary layer ($y_h^+ \gg 1$).

Langevin models provide a more physical representation of turbulence as particles see continuous fluctuating fluid velocity components. Several types of Langevin models are found in the literature (Simonin et al., 1993; Sommerfeld et al., 1993; Arcen and Tanière, 2009; Minier, 2015). Dehbi (2008) used a normalized Langevin model coupled with near wall Direct Numerical Simulation (DNS) fits to study particle deposition in vertical turbulent channel flows. Guingo and Minier (2008) employed a standard Langevin model and their results show an overprediction of deposition for low-inertia particles, as a results, a stochastic coherent structures model was employed for correction. Their model improved the predictions on deposition when compared with experiments. In addition to a Langevin model, Chibbaro and Minier (2008) also proposed an ad-hoc boundary condition for particles hitting the wall, to introduce features of coherent structures based on particle residence time in the near-wall region ($y^+ < 30$). A particle deposits when its particle residence time is greater than a characteristic time scale, whose parameters are obtained from DNS. Their findings demonstrated a good agreement with experimental data when the boundary condition is used.

Vertical channel flows have been widely utilized to evaluate dispersion models, whereas studies on horizontal channel flows are less common. Furthermore, there is a limited availability of research on deposition using wall-modeled approaches. This study aims to evaluate dispersion models and their predictions in terms of particles deposition when a coarse description of the boundary layer is adopted. These models are evaluated using the OpenFOAM®v2006 open-source CFD for turbulent channel flows (ESI Group, 2020), both in the vertical and horizontal orientations. Wall-modeled approaches are employed in the whole study and a stochastic wall deposition model (Aguinaga et al., 2009) is revisited, implemented and tested. In contrast with the work of Chibbaro and Minier (2008), this stochastic wall deposition model is applied above the wall in fully turbulent zones ($y^+ \sim 60$) where the Langevin models are applicable.

The article is structured as follows. First, we introduce the dispersion model currently available in OpenFOAM®, along with the single-step Langevin model proposed by Sommerfeld et al. (1993). These models are tested in both a simplified Homogeneous Isotropic Stationary Turbulence (Section 3) and turbulent channel flow scenarios (Section 4) to assess their capabilities in particle dispersion and deposition analysis. In Section 5 a stochastic wall deposition model is introduced. The model is then applied in the Section 6 for analyzing the deposition in horizontal and vertical channel. Conclusions are drawn in the last section.

## 2. Dispersion models

### 2.1. Eddy-Interaction Model (EIM)

The EIM also called “Eddy-Life-Time” model consists of modeling particle dispersion in turbulence. It is widely used in commercial CFD softwares when using RANS simulations. The interaction between a particle and turbulence is modeled by a series of eddies based on the statistical turbulent properties of the flow at its position. It is a stochastic discrete random walk treatment where each particle experiences a

fluctuating fluid velocity based on a random number sampled from a Gaussian distribution and scaled by the RMS of the fluctuating velocity component.

The interaction time between a particle and an eddy is limited by two criteria, namely, the eddy life time $\tau_e$ and the time required for a particle to cross the eddy (also called transit time $\tau_R$). The later takes into account the so-called “crossing trajectory effect”. It is assumed that at time $t$, a particle with velocity $\mathbf{u}_p$ is captured by an eddy which moves with a velocity composed of the mean fluid velocity, augmented by a random “instantaneous” component which is piecewise constant in time. When the lifetime of the eddy is over or the particle crosses the eddy, another interaction is generated with a different eddy, and so forth.

An EIM is implemented in OpenFOAM®v2006 as StochasticDispersionRAS. The interaction time $\tau_{int}$ is obtained from the following

$$ \tau_{int} = \min(\tau_e, \tau_R) \qquad \tau_e = \frac{k}{\epsilon} \qquad \tau_R = \frac{\lambda_e}{\|\mathbf{u}_r\|} \qquad (1) $$

$$ \lambda_e = \frac{C_\mu^{3/4} k^{3/2}}{\epsilon} $$

with $\mathbf{u}_r = \mathbf{u}_f@p - \mathbf{u}_p$, the relative fluid-particle slip velocity and $\lambda_e$ the eddy length scale. The isotropic turbulence modeling is adopted and during the interaction between the eddy and the particle, the sampled velocities are left unchanged and are obtained using

$$ \mathbf{u}'_{f@p} = \sigma_f |N_w| \mathbf{d}_r \qquad (2) $$

where $\mathbf{d}_r$ is a random isotropic direction vector and $N_w$ a random variable which follows a standard normal distribution $N(0, 1)$,

$$ \mathbf{d}_r = \begin{pmatrix} a \cos(\theta) \\ a \sin(\theta) \\ u \end{pmatrix} \qquad \theta = \omega_{r1} 2\pi \qquad u = 2\omega_{r2} - 1 \qquad a = \sqrt{1 - u^2} \qquad (3) $$

where $\omega_{r1}$ and $\omega_{r2}$ are two independent random numbers uniformly distributed between 0 and 1. For an isotropic turbulence, adopted here, the three RMS components of the fluid fluctuating velocity vector $\mathbf{u}'_{f@p}$ are equal to $\sigma_f = \sqrt{2k/3}$. This way of sampling the fluctuating velocity differs from standard EIMs gathered in Huilier (2021) review on dispersion models.

### 2.2. Langevin model

According to Pope (1994), in high Reynolds number, Langevin equation can be used to model fluid velocity fluctuations in homogeneous stationary turbulence. It represents the change in the velocity of a fluid element over time and can be broken down into two components: a damping (or ‘drift’) term, that is proportional to the velocity of the fluid and a random forcing term that has an average value of zero. The Langevin equation is a stochastic differential equation which uses Markov chains (Obukhov, 1959) for the fluctuating fluid velocity $\mathbf{u}'_f$

$$ \delta \mathbf{u}'_f = \mathbf{u}'_f(t + \delta t) - \mathbf{u}'_f(t) = -\frac{\mathbf{u}'_f}{\tau_f^t} \delta t + \sigma_f \sqrt{\frac{2}{\tau_f^t}} \delta \mathbf{W}. \qquad (4) $$

In (4), $\tau_f^t$ is the fluid Lagrangian timescale and $\delta t$ is the time increment. $\delta \mathbf{W}$ is a Wiener process having the following properties $\overline{\delta \mathbf{W}} = 0$ and $\overline{\delta \mathbf{W}^2} = \delta t$ where the overline bar is the time average operator. The discrete form of (4) correlates the fluid velocity fluctuation $\mathbf{u}'_f(t + \Delta t)$ to its value $\mathbf{u}'_f(t)$ at the previous time step (Sommerfeld et al., 1993)

$$ \mathbf{u}'_f(t + \Delta t) = \mathbf{u}'_f(t) R_L(\Delta t) + \frac{\sigma_f}{\sqrt{\Delta t}} \sqrt{1 - R_L^2(\Delta t)} \delta \mathbf{W}. \qquad (5) $$

Here $R_L(\Delta t)$ is the Lagrangian fluid velocity auto-correlation function given by

$$ R_L(\Delta t) = \exp\left(-\frac{\Delta t}{\tau_f^t}\right). $$
In the case of modeling the fluctuating fluid velocity at the particle position $\mathbf{u}'_{f@p}$, inertial and crossing-trajectory effects must be taken into account for inertial particles under gravity. In fact, they no longer follow the fluid elements and the Lagrangian autocorrelation function $R_f(\Delta t)$ in (5) must be modified to include space and time correlations through the following proposed model by Sommerfeld et al. (1993).

$$R_{EL,ij}(\Delta t, \Delta \mathbf{s}) = R_L(\Delta t) R_{E,ij}(\Delta \mathbf{s}).$$

With the Eulerian space correlation given by

$$R_{E,ij}(\Delta \mathbf{s}) = [f(\Delta \mathbf{s}) - g(\Delta \mathbf{s})] \frac{\Delta s_i \Delta s_j}{\|\Delta \mathbf{s}\|^2} + g(\Delta \mathbf{s}) \delta_{ij} \quad (6)$$

where $\Delta \mathbf{s}$ is the distance between a fluid element and the particle over a timestep $\Delta t$ (i.e. $\Delta \mathbf{s} = (\mathbf{u}_p - \mathbf{u}_{f@p})\Delta t$). Here $f(\Delta \mathbf{s})$ and $g(\Delta \mathbf{s})$ are respectively the longitudinal and transversal correlation coefficients given by

$$f(\Delta \mathbf{s}) = \exp\left(-\frac{\|\Delta \mathbf{s}\|}{L_E}\right) \quad (7)$$

$$g(\Delta \mathbf{s}) = \left(1 - \frac{\|\Delta \mathbf{s}\|}{2L_E}\right) \exp\left(-\frac{\|\Delta \mathbf{s}\|}{L_E}\right) \quad (8)$$

with the Eulerian length scale $L_E = C_L \tau_f' \sigma_f$ and the model constant $C_L$ set to 3. Sommerfeld et al. (1993) proposed $\tau_f' = C_T \sigma_f^2 / \epsilon$ where $C_T = 0.24$ is a fitted parameter. As described in Huilier (2021), the single step or the M.S. (Martin Sommerfeld) Langevin model (Sommerfeld et al., 1993) is written in this form

$$u'_{f@p}(t + \Delta t) = u'_{f@p}(t) R_{EL,li}(\Delta t, \Delta \mathbf{s}) + \frac{\sigma_f}{\sqrt{\Delta t}} \sqrt{1 - R_{EL,li}^2(\Delta t, \Delta \mathbf{s}) \Delta W_i} \quad (9)$$

where the Einstein convention of summing up over repeated indices is not adopted. To use Eq. (9), it is necessary for the matrix $R_{EL}$ to be diagonal. According to (6), this happens when $\Delta s_i \Delta s_j = 0$ for $i \neq j$ which correspond to $\Delta \mathbf{s}$ aligned along any axis of the reference system. In particular, this is verified when the drift induced by gravity is dominant, described by $\langle \|\mathbf{u}_p - \mathbf{u}_{f@p}\| \rangle_p = \tau_p g$, and when one axis of the reference system is aligned with gravity. To be used in any reference frame, we propose to use a tensorial form to generalize the formulation, so that

$$\mathbf{u}'_{f@p}(t + \Delta t) = \mathbf{u}'_{f@p}(t) \mathbf{R}_{EL}(\Delta t, \Delta \mathbf{s}) + \mathbf{A}_{f@p} \Delta \mathbf{W} \quad (10)$$

where $\mathbf{A}_{f@p}$ is a tensor, equal to a positive symmetric matrix given by

$$\mathbf{A}_{f@p} = \frac{\sigma_f}{\sqrt{\Delta t}} \sqrt{\mathbf{I} - \mathbf{R}_{EL}^2(\Delta t, \Delta \mathbf{s})} \quad (11)$$

and here $\mathbf{I}$ is the 3 × 3 identity matrix. It must be noticed that according to (6), (7) and (8) the eigenvalues of $\mathbf{I} - \mathbf{R}_{EL}^2$ are always positives allowing the computation of the square root.

The M.S. Langevin model, chosen for simplicity, presents some limitations, in particular in inhomogeneous turbulence and is limited to statistically stationary turbulence. In this paper, a specific work is conducted to enhance its predictions in turbulent channel flows. Other sophisticated Langevin models (Arcen and Tanière, 2009; Minier et al., 2014) can also be used but are not covered here. In addition, these dispersion models can also be coupled with the stochastic wall deposition model presented in this paper.

### 3. Dispersion in Homogeneous Isotropic Stationary Turbulence (HIST)

To assess particle dispersion models, it is common to use HIST and compare particle statistics with theoretical or numerical simulations obtained from DNS. Other test cases depending on the application can be used also such as round jets, mixing shear layers and vertical/horizontal channel flows. In this section, the different dispersion models will be investigated using HIST with and without gravity.

<figure><img src="image_3.png" /><figcaption>Fig. 1. View of a collection of 150 µm particles within the studied domain. The particles are colored by their velocities. Their sizes are not on scale. Periodic boundary conditions are applied on each face for the particles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Several monodispersed spherical particles are initially released with Stokes number $St$ roughly varying from $10^{-4}$ to 2. This correspond to diameters $d_p$ ranging from 2 to 150 µm, a fixed density $\rho_p = 2800$ kg/m$^3$, a turbulent kinetic energy $k = 0.07$ m$^2$/s$^2$ and a timescale $\tau_f' = 0.06$ s. Here the fluid is air under standard condition (i.e. density $\rho_f = 1.2$ kg/m$^3$ and kinematic viscosity $\nu_f = 1.568 \times 10^{-5}$ m$^2$/s). These parameters have been selected to realistically represent the scenario of soiling, focusing on dust particles depositing on solar panels. The numerical simulation is done in a three-dimensional domain where periodic boundary conditions are applied for particles. A total of 10,000 particles are initially released uniformly in the domain ($L_x = L_y = L_z = 1$ m) illustrated in Fig. 1.

#### 3.1. Particle tracking

From a numerical point of view, particles are considered, under the point-particle approximation, as dimensionless points of concentrated mass $m_p = \rho_p \pi d_p^3 / 6$ that corresponds to the mass of a sphere of diameter $d_p$ and density $\rho_p$. We consider a dilute regime (one-way coupling) where particles do not influence the fluid motion. We also neglect particle-particle collisions. Each particle is tracked individually using the following equation:

$$\frac{d\mathbf{u}_p}{dt} = -\frac{\mathbf{u}_p - \mathbf{u}_{f@p}}{\tau_p} + \left(\frac{\rho_p - \rho_f}{\rho_p}\right) \mathbf{g} \quad (12)$$

where the fluid velocity “seen” by the particle is expressed as

$$\mathbf{u}_{f@p} = \mathbf{u}_{f}(\mathbf{x}_{p}) + \mathbf{u}'_{f@p} \quad (13)$$

In (13), $\mathbf{u}_f(\mathbf{x}_p)$ is mean fluid velocity interpolated linearly at the particle position $\mathbf{x}_p$ and $\mathbf{u}'_{f@p}$ is the turbulent fluctuation modeled by a Lagrangian stochastic contribution given by a dispersion model. The non-Stokesian particle relaxation time (non-linear drag force) reads

$$\tau_p = \frac{\rho_p d_p^2}{18 \rho_f \nu_f} \frac{1}{1 + 1/6 Re_p^{2/3}} \quad (14)$$

where the particle Reynolds number is defined as $Re_p = d_p \|\mathbf{u}_{f@p} - \mathbf{u}_p\| / \nu_f$. The particle equation of motion in (12) is solved using an Euler implicit scheme implemented in OpenFOAM®. Each particle trajectory is obtained by advancing particle position:

$$d\mathbf{x}_p = \mathbf{u}_p dt. \quad (15)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. PDF of the fluctuating fluid velocity components normalized by $\sigma_f = \sqrt{2k/3}$ for 2 µm particles (left) and the normalized fluid kinetic energy “seen” $q_f^2_{@p}$ as a function of particle inertia (right) where the expected value should be equal to one. The solid black line on the left panel represents a Gaussian function with zero mean and a variance of $2k/3$. Statistics obtained using the OpenFOAM® EIM model without gravity.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. PDF of the fluctuating fluid velocity components normalized by $\sigma_f = \sqrt{2k/3}$ for 2 µm particles (left) and the normalized fluid kinetic energy “seen” $q_f^2_{@p}$ as a function of particle inertia (right). The solid black line on the left panel represents a Gaussian function with zero mean and a variance of $2k/3$. Statistics obtained using the M.S. Langevin model given by (9) without gravity.</figcaption></figure>

The OpenFOAM® solver used for the particle tracking is the IcoUnCoupledKinematicParcelFoam solution already implemented. In OpenFOAM®, a Face-to-Face algorithm is used to move the particles from initial position to final position in a time step $\Delta t$ based on a Lagrangian maximum Courant number $Co_{max}$ and a tracking time $\tau_{track}$ given by the user. These two parameters will set a maximum integration time step $\Delta t_{max} = Co_{max}\tau_{track}$. This way of obtaining the integration time step $\Delta t$ differs from what is usually done in the literature where it is updated in an efficient manner locally based on characteristic timescales. In OpenFOAM®, the maximum time step $\Delta t_{max}$ is chosen before running the simulations and in this paper we choose a value smaller than one tenth of the particle relaxation time in the Stokes regime $\tau_p^{St}$ (namely given by (14) with $Re_p = 0$) and the minimum value of $\tau_f$ in the domain.

## 3.2. Particle statistics

The average particle relaxation time, $\tau_{fp}^F$, is defined as

$$ \frac{1}{\tau_{fp}^F} = \left\langle \frac{1}{\tau_p} \right\rangle_p . \qquad (16) $$

In (16), $\langle \cdot \rangle_p$ is the arithmetic average over an ensemble of $N_p$ particles

$$ \langle \cdot \rangle_p = \frac{1}{N_p} \sum_{n=1}^{N_p} (\cdot)^{(n)}. \qquad (17) $$

Here $N_p$ is defined as the cumulative number of particles whose center is in a control volume $V$ at different times instants. In the present analysis the fluid–particle Stokes number $St = \tau_{fp}^F/\tau_f'$ is used to evaluate particle inertia. The fluid kinetic energy “seen” by the particles, i.e. $q_f^2_{@p} = \langle u_f'^2 \rangle_p/2$, by assumption should be equal to the prescribed turbulent kinetic energy $k$. To evaluate particle dispersion, we verify the Tchen–Hinze theory (Hinze, 1975) where the particles are expected to be in-equilibrium with stationary turbulence

$$ \langle u_p'^2 \rangle_p = f(St) \langle u_f'^2 \rangle_p $$

where $f(St)$ is a function of the Stokes number. Its derivation is based on the form of the Lagrangian autocorrelation function $R_f$. For an exponentially decaying function corresponding to a Langevin model it leads to

$$ f(St) = \frac{1}{1 + St} . \qquad (18) $$

This is also valid for an EIM where the eddy lifetime is sampled from an exponential distribution with a mean value $\tau_f'$ (Kallio and Reeks, 1989).
<figure><img src="image_6.png" /><figcaption>Fig. 4. Normalized fluid kinetic energy “seen” by the particle $q_{f@p}^2$ as a function of particle inertia using (9) left and (10) right.</figcaption></figure>

### 3.3. Results

Fig. 2 shows the PDF of the fluctuating velocity components for the OpenFOAM® dispersion model and a Gaussian function with zero mean and a variance of $2k/3$. The model fails to reproduce Gaussian turbulence. This discrepancy comes from (2) where the modulus of a Gaussian variable is multiplied by a vector uniformly distributed in space leading to a non Gaussian distribution. The fluid kinetic energy “seen” is also computed and compared with $k$ for all particle classes. The OpenFOAM® model is not able to simulate the correct fluid agitation and underestimate it by 70%. Taking the absolute value of a random number sorted from a centered Gaussian distribution changes the variance to $1 - 2/\pi \approx 0.36$.

For the M.S. Langevin model, the Gaussian shape is found (see Fig. 3). The fluid kinetic energy “seen” $q_{f@p}^2$ is computed using both formulations of the model, namely (9) and (10), with and without gravity. The results are shown on Fig. 4. We can clearly see the inability of the model (9) to retrieve the correct fluid kinetic energy “seen” when the Stokes number increase without gravity. In fact, in this scenario, the instantaneous fluid-particle relative velocity $u_r = u_f@p - u_p$ and hence $\Delta s$ can be in any direction: $R_{EL}$ is not diagonal. When gravity is added in the y axis, the error is reduced as $\Delta s$ tends to be aligned with the y-direction. On the other hand, the generalized formulation (10) behaves well in both cases as it uses the complete $R_{EL}$ tensor and will be adopted in this paper.

Fig. 5 illustrates the particle-fluid variance ratio and its comparison with the Tchen-Hinze theory extension with gravity. Overall, the model matches the theory, particles are in equilibrium with turbulence. The decrease in particle agitation is well captured. It is important that the Stokes number $\tau_{fp}^F/\tau_f^f$ here is modified to $\tau_{fp}^F/\tau_{f@p}^f$. Following Csanady (1963), in case where a mean drift induced by gravity is present and aligned with the y-axis

$$ \tau_{f@p,x}^t = \tau_{f@p,z}^t = \frac{\tau_f^t}{\sqrt{1 + 4\beta_c^2 \left(\frac{\tau_{fg}^g}{\sigma_f}\right)^2}} \quad (19) $$

in the transversal directions and

$$ \tau_{f@p,y}^t = \frac{\tau_f^t}{\sqrt{1 + \beta_c^2 \left(\frac{\tau_{fg}^g}{\sigma_f}\right)^2}} \quad (20) $$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Comparison of the particle-particle kinetic stress diagonal components with Tchen-Hinze theory given by (18). Results are obtained using the M.S. Langevin model with gravity acting in the y-direction. The stokes number is computed with the timescale given by (19) and (20).</figcaption></figure>

in the longitudinal direction where $\beta_c = 1/C_L$ is a constant. These quantities are computed in each direction in the case where gravity is acting in the y-direction and used to normalize $\tau_{fp}^F$. The discrepancy in the intermediate region ($0.1 < \tau_{fp}^F/\tau_f^f < 1$) is due to the crossing trajectory effect modeled through $R_{EL}$ which is based on the instantaneous relative velocity rather than using a Csanady (1963) form in terms of the mean relative velocity (Simonin et al., 1993; Huilier, 2021).

### 4. Fluid elements dispersion in turbulent channel flow

In this section, we investigate fluid elements statistics within a two-dimensional turbulent channel flow. It is a common test case to evaluate the capability of a dispersion model to retrieve the correct fluid agitation when turbulence is inhomogeneous in the tracer limit
<table><caption>**Table 1**
Properties of the fluid and of the examined channel flow.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Simulation track time</td><td>τ<sub>track</sub></td><td>10<sup>-4</sup></td><td>s</td></tr><tr><td>Lagrangian maximum Courant number</td><td>C<sub>max</sub></td><td>0.1</td><td>-</td></tr><tr><td>Fluid kinematic viscosity</td><td>ν<sub>f</sub></td><td>1.568 × 10<sup>-5</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Fluid density</td><td>ρ<sub>f</sub></td><td>1.2</td><td>kg m<sup>-3</sup></td></tr><tr><td>Fluid bulk velocity</td><td>U<sub>∞</sub></td><td>6</td><td>m s<sup>-1</sup></td></tr><tr><td>Fluid shear velocity</td><td>u<sub>τ</sub></td><td>0.27</td><td>m s<sup>-1</sup></td></tr><tr><td>Friction Reynolds number</td><td>Re<sub>τ</sub></td><td>1530</td><td>-</td></tr><tr><td>Channel Reynolds number</td><td>Re</td><td>30 612</td><td>-</td></tr><tr><td>Channel length</td><td>L</td><td>5</td><td>m</td></tr><tr><td>Channel width</td><td>w</td><td>0.08</td><td>m</td></tr><tr><td>Cell size</td><td>Δ<sub>x</sub> = Δ<sub>y</sub> = y<sub>h</sub></td><td>0.005</td><td>m</td></tr><tr><td>Wall cell height in dimensionless units</td><td>y<sub>h</sub><sup>+</sup></td><td>100</td><td>-</td></tr></tbody></table>

($τ_p → 0$) (MacInnes and Bracco, 1992). The well-mixed criterion is also investigated: very low-inertia particles which are uniformly mixed in the domain should remain well mixed as time evolves. Minier et al. (2014) set the well-mixed condition as a fundamental criterion for dispersion models. In OpenFOAM®, fluid elements tracking can be performed by prescribing directly the fluid velocity “seen” to the particles (i.e. $\mathbf{u}_p = \mathbf{u}_{f@p}$), (12) is not solved in this case.

## 4.1. Fluid flow field

The fluid flow is obtained from converged RANS $k-ε$ turbulence model where periodic boundary conditions are applied to obtain a fully developed flow. The OpenFOAM® incompressible steady state solver simpleFOAM was used. Uniform hexahedral meshes are employed. The desired wall cell dimensionless height $y_h^+$ is calculated using the Fanning coefficient:

$$y_h^+ = \frac{y_h u_τ}{\nu_f} \qquad u_τ = U_{f,x} \sqrt{\frac{\Lambda}{8}} \qquad \Lambda = 0.3164 Re^{-1/4} \\ 
Re = \frac{w U_{f,x}}{\nu_f} \qquad (21)$$

The wall-cell height can be obtained from a given value of $y_h^+$. Table 1 shows the fluid properties and the examined channel flow for $y_h^+ = 100$. The cell size is uniform in both directions $Δ_x = Δ_y$ and has been set to $5 × 10^{-3}$ m. It leads to 15 cells in the y-direction (wall-normal). In x-direction the domain length is 5 m. As explained in Section 6.1 the accurate measurement of the particle deposition requires a very long channel for the low Stokes number. Consequently, in x-direction the mesh has 1000 cells. Dimensionless fluid streamwise velocity, $U^+ = U_{f,x}/u_τ$, is shown by Fig. 6 and dimensionless turbulent kinetic energy, $k^+ = k/u_τ^2$, and dissipation rate, $ε^+ = εν_f/u_τ^4$, by Fig. 7.

## 4.2. Particle statistics

Fluid particles ($N_p = 10,000$) are initially released uniformly within the channel where periodic boundary condition is applied. Elastic rebound of particles is also applied at the walls. For averaging, the channel is binned in the y-direction where each bin height is equal to $y_h$. The particle number density $n_p$ (number of particles present in the bin divided by its volume) is computed in each bin and normalized by the channel number density $n_{p,0}$. The fluid element wall-normal agitation $\langle u_{p,y}^2 \rangle_p$ is also computed in each bin and compared with the prescribed one (i.e. $\sigma_f^2 = 2k/3$).

## 4.3. Results

As shown by Fig. 8, the fluid particle concentration is found to be not uniform across the channel. Fluid particles tend to accumulate in the channel center where turbulence levels are the lowest. This

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Streamwise fluid velocity $U^+$ profile as a function of the channel width $w$. Each data point is taken from the cell-centroids.</figcaption></figure>

“spurious drift” is an unwanted effect that needs to be corrected by the dispersion models in the tracer limit. In addition, the particle agitation is not fully retrieved. An increase is seen at the center of the channel where turbulence is almost homogeneous (see Fig. 7). This problem was in fact identified by MacInnes and Bracco (1992): when turbulence is inhomogeneous, dispersion models may produce an artificial transfer of fluid turbulence from the region of high intensity to low intensity. This is the case here, turbulence levels are at the lowest in the center of the channel.

It is possible to account for the anisotropy of the fluid phase by using LRR closure (Innocenti et al., 2019) or by using a normalized dispersion model as proposed by Dehbi (2008). The normalized version of the M.S. Langevin can be written following Dehbi (2008) and Minier et al. (2014) as

$$\delta \left( \frac{\mathbf{u}'_{f@p}}{\sigma_f} \right) = - \left( \frac{\mathbf{u}'_{f@p}}{\sigma_f} \right) \frac{\delta t}{\tau_f'} + \sqrt{\frac{2}{\tau_f'}} \delta \mathbf{W} \qquad (22)$$

and the discretized form following (10)

$$\mathbf{u}'_{f@p}(t + \Delta t) = \frac{\sigma_f(\mathbf{x}_p(t + \Delta t))}{\sigma_f(\mathbf{x}_p(t))} \mathbf{u}'_{f@p}(t) \mathbf{R}_{EL}(\Delta t, \mathbf{\Delta s}) + \mathbf{A}_{f@p} \mathbf{\Delta W} \qquad (23)$$

Fig. 8 shows the comparison between the two versions of the dispersion model. The agitation is fully retrieved and no increase in the center of the channel is seen. In addition, when adopting the normalized version, a slight improvement is seen for the particle concentration. The non-uniformity of the particle concentration is in fact attributed to the absence of a drift term in (22). Dehbi (2008) added a drift term to overcome the “spurious drift” effect

$$\delta \left( \frac{\mathbf{u}'_{f@p}}{\sigma_f} \right) = - \left( \frac{\mathbf{u}'_{f@p}}{\sigma_f} \right) \frac{\delta t}{\tau_f'} + \sqrt{\frac{2}{\tau_f'}} \delta \mathbf{W} + \frac{1}{3\sigma_f} \cdot \frac{\partial k}{\partial \mathbf{x}} \cdot \frac{\delta t}{1 + St} \qquad (24)$$

and the discretized form

$$\mathbf{u}'_{f@p}(t + \Delta t) = \frac{\sigma_f(\mathbf{x}_p(t + \Delta t))}{\sigma_f(\mathbf{x}_p(t))} \mathbf{u}'_{f@p}(t) \mathbf{R}_{EL}(\Delta t, \mathbf{\Delta s}) + \mathbf{A}_{f@p} \mathbf{\Delta W} \\ + \frac{\sigma_f(\mathbf{x}_p(t + \Delta t))}{3\sigma_f(\mathbf{x}_p(t))} \cdot \frac{\partial k}{\partial \mathbf{x}} \cdot \frac{\Delta t}{1 + St} \qquad (25)$$

where $\mathbf{A}_{f@p}$ is given by (11) with the fluid velocity RMS velocity $\sigma_f(\mathbf{x}_p(t + \Delta t))$. Fig. 9 shows the results when adding the drift term. It is clearly seen that fluid elements exhibit now a uniform concentration profile. In addition, the wall-normal agitation is unchanged. This normalized form of the M.S. Langevin with a drift term will be used in paper.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Profiles of turbulent kinetic energy $k^+$ profile (left) and dissipation rate $\epsilon^+$ (right) as a function of the channel width $w$. Each data point is taken from the cell-centroids.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Profiles of fluid particle number density (left) and the wall-normal agitation (right) for the two versions of the M.S. Langevin model. The black empty symbols are given by (10) and the filled symbols by (23).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Profiles of fluid particle number density (left) and the wall-normal agitation (right) for the M.S. Langevin model. The black filled symbols are given by (23) and the half filled symbols by (25).</figcaption></figure>
## 5. Particle deposition in horizontal and vertical channel

In this section the prediction of particle deposition in horizontal and vertical turbulent channel flows is analyzed. The deposition velocity $v_d$ can be measured for each class of particles with a relaxation time $\tau_p$. These quantities are made dimensionless using the shear velocity $u_\tau$ and the kinematic viscosity of the fluid $\nu_f$.

$$v_d^+ = \frac{v_d}{u_\tau} = \frac{J}{C_0 u_\tau} \qquad \tau_p^+ = \frac{\tau_p^{St} u_\tau^2}{\nu_f} \qquad \tau_p^{St} = \frac{\rho_p d_p^2}{18 \rho_f \nu_f} \\ u_\tau = \sqrt{\frac{\tau_\omega}{\rho_f}}.$$

Here, $J$ is the deposition rate, $C_0$ is the bulk concentration of particles, $\tau_\omega$ the wall shear stress and $\tau_p^{St}$ is the particle Stokes response time. Deposition of particles on vertical channels (gravity acting in the direction of the flow) has been studied experimentally and numerically (Sippola and Nazaroff, 2002). Three regimes depending on $\tau_p^+$ have been identified (Papavergos and Hedley, 1984):

* Diffusion regime ($\tau_p^+ \ll 1$): Particles with low-inertia respond to the turbulent boundary layer. Their deposition is controlled by Brownian motion.

* Inertial regime ($\tau_p^+ \gg 40$): Particles are very inertial and do not interact with the turbulent boundary layer.

* Diffusion-impaction regime ($1 < \tau_p^+ < 40$): Transitional regime which corresponds to a partial interaction of particles with the turbulent boundary layer.

Nerisson et al. (2011) derived the following analytical model for $v_d^+$ at a given dimensionless distance $y_h^+$ of the wall

$$v_d^+ = \frac{J}{u_\tau C(y_h^+)} = \left[ \frac{Sc_l^f \ln(y_h^+) + \frac{\omega}{\tau_p^{+2}}}{\kappa} \right]^{-1} \quad (26)$$

where $Sc_l^f$ is the Schmidt turbulent number, $\kappa$ the Von Karman constant, $y_h^+$ wall computational cell height and $\omega$ a constant. Fig. 10 illustrates the deposition velocity $v_d^+$ from the experiments of Liu and Agarwal (1974) and model in (26) for $Sc_l^f = 1$ and $\omega = 2.2 \times 10^3$ with $y_h^+ = 100$ and $y_h^+ = 30$. The model shows an overall good agreement with the experiments. Here, only the diffusion-impaction and inertial regimes are considered.

In horizontal channels, the direction of fluid flow is perpendicular to the direction of gravitational forces. As a result, gravity assists in deposition within the channel. An analysis of the deposition velocity $v_d^+$ in a horizontal channel indicates the existence of three distinct regimes (Gao et al., 2012). The relative widths of these regimes differ from their counterparts in a vertical channel, with the diffusion regime being narrowed ($\tau_p^+ < 10^{-3}$), and the diffusion-impaction regime being broadened due to gravity ($10^{-3} < \tau_p^+ < 20$). In the inertial regime ($\tau_p^+ > 20$), the deposition velocity of particles increases as a function of particle size. In addition, in horizontal channels deposition becomes a function of an additional number $g^+ = \tau_p^{St} g / u_\tau$ that is a Froude number. With this parameter, Wood (1981) proposed the following empirical equation

$$v_d^+ = 4.5 \times 10^{-4} \tau_p^{+2} + \|g^+\|. \quad (27)$$

Later, Nerisson et al. (2011) derived a generalized formula for any wall orientation:

$$v_d^+ = -g^+ \cdot n \left[ 1 - \exp \left( \frac{g^+ \cdot n Sc_l^f}{\kappa} \ln(y_h^+) + g^+ \cdot n \frac{\omega}{\tau_p^{+2}} \right) \right]^{-1} \quad (28)$$

where $n$ is the wall normal unit vector. Fig. 11 compares these two models with experiments (Montgomery and Corn, 1970; Sehmel, 1973; Kvasnak et al., 1993) where the friction velocity is imposed ($u_\tau \approx 0.3$ m/s) and $\tau_p^{St}$ varies. Both models show an overall good agreement

<figure><img src="image_12.png" /><figcaption>Fig. 10. Dimensionless deposition velocity $v_d^+$ as a function of particle inertia $\tau_p^+$ for the vertical channel. The symbols are the experimental data (▲) from Liu and Agarwal (1974) and the solid line (—) the prediction given by (26) from Nerisson et al. (2011).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Dimensionless deposition velocity $v_d^+$ as a function of particle inertia $\tau_p^+$ for the horizontal channel. Experiments (▲), Nerisson et al. (2011) model (—) and Wood (1981) model (— - -). The dotted line (---) here represents the particle settling velocity. In this case with gravity, the predictions of Nerisson et al. (2011) model and the settling velocity are superimposed.</figcaption></figure>

with the experiments. For larger values of $\tau_p^+$, a deviation is seen, Wood (1981) model gives higher values of $v_d^+$. Furthermore, if we draw the gravity settling velocity for these particles represented by the dotted line, we see a clear match with data and both models. This indicates that for this experiment, gravity is the dominant deposition mechanism. Other experiments with larger $u_\tau$ or lower $\tau_p^{St}$ values may reveal the contribution of turbulence in deposition.

Lagrangian stochastic simulations of particle deposition in vertical turbulent channel flows with $y_h^+ \gg 1$ have been studied previously. Parker et al. (2008) used the Reynolds-Stress Model and $y_h^+ = 120$ while Aguinaga et al. (2009) used the $k-\epsilon$ model and $y_h^+ = 30$ and 60. Both simulations were conducted using the dispersion model of Fluent (2006). Fig. 12 shows the numerical predictions of deposition for different particle inertia. In both simulations, the tendency of an
<figure><img src="image_14.png" /><figcaption>Fig. 12. Dimensionless deposition velocity $v_d^+$ as a function of particle inertia $\tau_p^+$ for the vertical channel. Experiments (▲) from Liu and Agarwal (1974) and numerical simulations from Aguinaga et al. (2009) and Parker et al. (2008).</figcaption></figure>

horizontal straight line is visible. Higher values of $\tau_p^+$ seem to give good agreement with experiments. However, for low-inertia particles a deviation is seen. In fact, this is due to the coarse description of the boundary layer as particles do not see the correct fluid agitation in the near-wall region $y^+ < 5$. There is a need to overcome this issue in applications where the mesh size cannot be refined to fully resolve the boundary layer.

A Lagrangian stochastic wall deposition model was proposed by Aguinaga et al. (2009) and consists of treating particles entering the computational wall cell. The model has already been implemented and assessed for vertical channels in Fluent (2006) coupled with an EIM (Aguinaga et al., 2009). The model assigns a deposition probability for any incoming particle based on its intrinsic properties and flow characteristics. This probability $r = J/J^-$ is calculated theoretically from the particle flux ratio between the incoming ($J^-$) and the deposited flux ($J$) as shown by Fig. 13. The latest can be either obtained from DNS or experiments. The model uses the particle deposited flux obtained from Nerisson et al. (2011) for any wall configuration:

$$J = C(y_h)u_{\tau}v_{d}^{+} \qquad (29)$$

In (29), $C(y_h)$ is the particle concentration at the height of wall computational cell. The incoming particle flux $J^-$ is obtained by assuming a Gaussian PDF of wall-normal particle velocities where its statistical moments are defined below:

$$J^{-} = -\frac{C(y_h)}{2} \left[ \sqrt{\frac{2}{\pi} \langle u_{p,y}^2 \rangle_p} \exp\left(-\frac{U_{p,y}^2}{2\langle u_{p,y}^2 \rangle_p}\right) - U_{p,y} \left(1 - \operatorname{erf}\left(\frac{U_{p,y}}{\sqrt{2\langle u_{p,y}^2 \rangle_p}}\right)\right) \right].$$

The mean wall-normal particle velocity for any wall orientation is given by

$$U_{p,y} = \frac{J}{C(y_h)} = u_{\tau}v_{d}^{+}.$$

The wall-normal particle velocity variance is obtained from the Tchen-Hinze theory, particles agitation is in equilibrium with the fluid turbulence. In the case of particle dispersion modeled by a Langevin type

<figure><img src="image_15.png" /><figcaption>Fig. 13. Sketch of the incident ($J^-$), reflected ($J^+$) and deposited ($J$) particle fluxes at the wall computational cell.</figcaption></figure>

approach

$$\langle u_{p,y}^2 \rangle_p = \frac{1}{1 + St} \langle u_{f@p,y}^2 \rangle_p.$$

Assuming a homogeneous isotropic turbulence at $y_h$, the wall-normal fluid velocity “seen” variance is given to be $2k/3$. The model requires $k$ and $\epsilon$ values given by the RANS simulation at $y_h$ and the dimensionless intrinsic parameters of the particle, namely $\tau_p^+$ and $g^+ \cdot n$. A Lagrangian stochastic deposition model is employed where a uniform random variable $0 < \omega_r < 1$ is sorted and its comparison with $r$ will determine the fate of each particles entering the wall cell and have a velocity vector pointing towards the wall, i.e. $u_p \cdot n < 0$:

* $r < \omega_r$: Particle is ejected back to the bulk flow modeled by an elastic rebound;

* $r > \omega_r$: Particle is deposited (see Fig. 14).

In this study, the deposition model is applied at an arbitrary critical height $y_c^+$ instead of the wall cell height $y_h^+$. This approach reduces the dependency of particle treatments on the mesh resolution and type. The potential outcomes for particles are visually represented in Fig. 14. The algorithmic process is depicted in the flowchart presented by Fig. 15. During the integration of particle trajectories at any given time $t$, the particle’s distance to the nearest wall and its corresponding $y_p^+$ value are calculated. If the particle is located within the near wall region, namely $y_p^+ < y_c^+$ and $u_p \cdot n < 0$, the deposition probability $r$ is determined using Eulerian fields at the particle’s location. The fate of the particle is determined by the random variable $\omega_r$. In this study, we added an elastic rebound for the fluid velocity “seen”. By doing this, the particle fluid covariance $\langle u_p u_f @ p \rangle_p$ is left unchanged.

# 6. RANS simulation of particle deposition in horizontal and vertical channel

## 6.1. Particle injection and deposition measurement

Monodispersed particles are injected uniformly at the inlet (see Fig. 16). The deposition velocity $v_d$ can be measured for each class of particles with a relaxation time $\tau_p$. These quantities are made dimensionless using the friction velocity $u^*$ estimated from the value the turbulent kinetic energy at the first cell point $k_p$ and the kinematic viscosity of the fluid $\nu_f$. In fact, the shear/friction velocity $u_\tau$ can be approximated by $u^*$ obtained from the log-region assumption. In a channel flow both quantities are close, but in more complex flows $u^*$ is
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Illustration of the potential scenarios when a particle crosses $y_c^+$ during its trajectory. $\bar{u}_p^-$ represents the particle velocity prior to the application of the deposition model, while $\bar{u}_p^+$ denotes the particle velocity after the treatment.</figcaption></figure>

more relevant of the turbulence in the near wall region. For example at a reattachment point $u_τ \to 0$, not $u^*$

$$v_d^+ = \frac{v_d}{u^*} = \frac{J}{C_0 u^*} \qquad \tau_p^+ = \frac{\tau_p^{St} u^{*2}}{\nu_f} \qquad u^* = \sqrt{k_p \sqrt{C_\mu}} \quad (30)$$

Here $J$ is the deposition rate and $C_0$ is a reference concentration of particles. Two standard channel cases will be used in this study, one vertical and the other horizontal where gravity is assisting deposition. The deposition velocity $v_d^+$ can also be calculated in a fully developed turbulent channel flow for each particle class (Kallio and Reeks, 1989; Matida et al., 2000) with

$$v_d^+ = \frac{1}{N_s} \sum_{i=1}^{N_s} \frac{1}{u^*} \frac{S_a U_\infty}{S_d} \ln \left( \frac{M_{i+1}}{M_i} \right) \quad (31)$$

where $S_a$ the channel section area, $S_d$ the deposition surface area, $U_\infty$ the bulk fluid velocity, $M_i$ the particle mass flux at the $i$th section and $N_s$ the number of sections (see Fig. 16). Particles are injected with the fluid mean velocity, the settling velocity $\tau_p g$, a stochastic process that is based on the dispersion model $u'_{f@p}$ and a new stochastic process based on the fluid-particle covariance. The initial particle velocity is given by

$$u_{p,i} = U_{f,i} + \alpha u'_{f@p,i} + \beta \gamma_i + \tau_p g_i \quad (32)$$

where $\gamma$ is a vector of three independent Gaussian random numbers with a mean value of zero and a standard deviation of unity. The two model's constant, $\alpha$ and $\beta$, reads

$$\alpha = f(St) \qquad \beta = \sqrt{f(St) \left[ \langle u'_{f@p,i} \rangle_p - \langle u'_{f@p,i} \rangle_p f(St) \right]} \quad (33)$$

The first mass flux measurement probe is located at $L/2$ to ensure that particles reach an equilibrium (i.e. $\partial n_p / \partial x = 0$ and $\partial(u'_p)^2 / \partial x = 0$).

## 6.2. Deposition results

### 6.2.1. Vertical channel

The experiment of Liu and Agarwal (1974) are done with gravity acting in the flow direction. Some authors (Kallio and Reeks, 1989; Parker et al., 2008; Dehbi, 2008; Guingo and Minier, 2008) do not consider gravity in their simulations while others (Matida et al., 2000) include gravity but neglect crossing-trajectory effects (not including the transit time $\tau_R$) in their dispersion models. For the sake of clarity, the two cases will be investigated. This will reveal the true behavior of such dispersion models in predicting deposition. The particle properties and dimensionless quantities are shown in Table 2. The deposition results of the different dispersion models are shown by Fig. 17.

In the case where gravity is neglected, the M.S. Langevin is able to match the experiments in the inertial regime. On the other hand, the OpenFOAM® EIM model underestimate deposition. In fact, the

<caption>**Table 2**</caption>

Material properties of the particles in the vertical channel. The dimensionless particle size is defined as $d_p^+ = d_p u_τ / \nu_f$ and is very small compared to the cell size $y_h^+$.

<table><thead><tr><th>Class</th><th>d<sub>p</sub> [µm]</th><th>τ<sub>p</sub><sup>+</sup> [-]</th><th>d<sub>p</sub><sup>+</sup> [-]</th></tr></thead><tbody><tr><td>1</td><td>2</td><td>0.158</td><td>0.034</td></tr><tr><td>2</td><td>3</td><td>0.355</td><td>0.052</td></tr><tr><td>3</td><td>4</td><td>0.631</td><td>0.069</td></tr><tr><td>4</td><td>5</td><td>0.987</td><td>0.086</td></tr><tr><td>5</td><td>8</td><td>2.526</td><td>0.138</td></tr><tr><td>6</td><td>10</td><td>3.947</td><td>0.172</td></tr><tr><td>7</td><td>12</td><td>5.683</td><td>0.207</td></tr><tr><td>8</td><td>20</td><td>15.79</td><td>0.344</td></tr><tr><td>9</td><td>30</td><td>35.52</td><td>0.517</td></tr><tr><td>10</td><td>40</td><td>63.14</td><td>0.689</td></tr><tr><td>11</td><td>50</td><td>98.66</td><td>0.861</td></tr><tr><td>12</td><td>60</td><td>142.08</td><td>1.033</td></tr><tr><td>13</td><td>90</td><td>319.67</td><td>1.55</td></tr><tr><td>14</td><td>150</td><td>887.98</td><td>2.583</td></tr></tbody></table>

underestimation of deposition is attributed to its inability to reproduce the correct fluid agitation (underestimated by 70% in HIST simple case).

In the case where gravity is acting in the flow direction which is the actual experiment of Liu and Agarwal (1974), deposition is still underpredicted by the OpenFOAM® EIM in the inertial regime whereas the M.S. Langevin model gives better agreement. In the presence of gravity or other body forces, the relative particle-fluid slip velocity $u_τ$ becomes significant for large $\tau_p^+$ particles. This effect decreases the transit time $\tau_R$ for the EIMs, eddy-particle interactions are more frequent which leads to less deposition.

Both models fail to reproduce the diffusion-impaction (low-inertia) regime where $v_d^+$ scales with $\tau_p^{+2}$. In fact, when using industrial scale mesh size namely $y_h^+ \gg 1$, particles see an inaccurate description of the turbulent boundary layer. Therefore the complex interactions between the particle and the fluid in the near wall region are not captured. This is given to the values of $k$ and $\epsilon$ along the particle trajectory between the first cell point and the wall which are not accurate enough. OpenFOAM® does not interpolate the values of $k$ and $\epsilon$ at the particles positions which make them see uniform values within a computational cell (even the wall-cell). This will overpredict particle dispersion and hence deposition.

### 6.2.2. Horizontal channel

The experiments of Montgomery and Corn (1970), Sehmel (1973) and Kvasnak et al. (1993) ($u_τ \approx 0.3$ m/s) are grouped and compared by Fig. 18 with numerical results. The particle properties and dimensionless quantities are shown in Table 3. As expected, deposition is consistently overestimated in all dispersion models when dealing with low-inertia particles ($\tau_p^+ < 10$), primarily due to the coarse representation of the boundary layer. However, in this experimental study,
<figure><img src="image_17.png" /><figcaption>Fig. 15. Algorithm for the deposition model handling arbitrary 3D geometries.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Schematic of the channel used in the numerical simulation. The vertical red lines are the mass flux probes for deposition velocity calculations.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Dimensionless deposition velocity $v_d^+$ with respect to the particle inertia $\tau_p^+$. The black-filled symbols are the experiments from Liu and Agarwal (1974). The left panel correspond to the case without gravity and the right-panel with gravity added in the flow direction. Results are obtained in the vertical channel configuration.</figcaption></figure>

<caption>**Table 3**</caption>

Material properties of the particles in the horizontal channel. The dimensionless particle size is defined as $d_p^+ = d_{ph}/\nu_f$ and is very small compared to the cell size $y_h^+$.

<table><thead><tr><th>Class</th><th>d<sub>p</sub> [µm]</th><th>τ<sub>p</sub><sup>+</sup> [-]</th><th>d<sub>p</sub><sup>+</sup> [-]</th></tr></thead><tbody><tr><td>1</td><td>0.75</td><td>0.022</td><td>0.013</td></tr><tr><td>2</td><td>1</td><td>0.039</td><td>0.017</td></tr><tr><td>3</td><td>1.5</td><td>0.089</td><td>0.026</td></tr><tr><td>4</td><td>2</td><td>0.16</td><td>0.034</td></tr><tr><td>5</td><td>3</td><td>0.36</td><td>0.052</td></tr><tr><td>6</td><td>4</td><td>0.63</td><td>0.07</td></tr><tr><td>7</td><td>5</td><td>0.98</td><td>0.087</td></tr><tr><td>8</td><td>6</td><td>1.42</td><td>0.1</td></tr><tr><td>9</td><td>10</td><td>3.95</td><td>0.17</td></tr><tr><td>10</td><td>12</td><td>5.68</td><td>0.21</td></tr><tr><td>11</td><td>15</td><td>8.88</td><td>0.26</td></tr><tr><td>12</td><td>20</td><td>15.79</td><td>0.35</td></tr><tr><td>13</td><td>25</td><td>24.67</td><td>0.44</td></tr><tr><td>14</td><td>30</td><td>35.52</td><td>0.52</td></tr></tbody></table>

gravity dominates the deposition process, where particles settle with a velocity determined by $\tau_p g$. This phenomenon becomes evident when the dispersion models are deactivated (i.e. no wall-normal particle agitation), resulting in even more accurate outcomes.

## 6.3. Results with the stochastic wall deposition model

The stochastic wall deposition model is now applied at a critical height $y_c^+ = 60$. Two meshes are used to test the sensitivity of the model: $y_h^+ = 50$ and $y_h^+ = 100$, respectively 30 and 15 cells across the width of the channel. The M.S. Langevin dispersion model is employed at the bulk. Fig. 19 shows the results for both the vertical ($g_x = 9.81 \text{ m/s}^2$) and horizontal channels ($g_y = 9.81 \text{ m/s}^2$) with and without the deposition model. When the stochastic wall deposition model is used, we have a better agreement with the experiments for both channel configurations for the two mesh cases ($y_h^+ = 50$ and $y_h^+ = 100$). The diffusion-impaction regime (i.e. $1 < \tau_p^+ < 40$) is now well reproduced. These particles have lower probability of deposition $r$ given by the stochastic wall deposition model than inertial ones ($\tau_p^+ > 40$). When they cross the critical height $y_c^+$, they are more likely to be ejected back into the bulk flow. In addition, the results show also low dependencies on the mesh size. Hence, the coupling between the M.S. Langevin model and the deposition model gives accurate deposition results for wall-modeled approaches. This stochastic wall deposition model can also be coupled with other complete Langevin dispersion models (Minier et al., 2014).

Fig. 20 shows the impact on the deposition results when not employing the normalized version of the M.S. Langevin from (24). An

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Dimensionless deposition velocity $v_d^+$ with respect to the particle inertia $\tau_p^+$. The black-filled symbols are the experiments from Montgomery and Corn (1970), Sehmel (1973) and Kvasnak et al. (1993). Results are obtained in the horizontal channel configuration.</figcaption></figure>

overall decrease in deposition is seen for all particle classes. In fact this is attributed to the deficiency of the model for maintaining the correct fluid agitation in inhomogeneous turbulence seen previously.

## 7. Conclusions

RANS dispersion models have been assessed in Homogeneous Isotropic Stationary Turbulence (HIST) and for deposition in turbulent channel flows using a wall-modeled approach ($y_h^+ \gg 1$). The investigation focused on dilute particulate flows and employed a combination of steady-state RANS $k-\epsilon$ turbulence modeling and Lagrangian particle tracking. The numerical simulations were conducted using the open-source solver, OpenFOAM®.

The analysis of particle statistics in HIST reveals that the OpenFOAM® EIM dispersion model, namely *StochasticDispersionRAS*, fails to accurately predict the expected fluid kinetic energy “seen” by the particles, accounting for only 30% of the total kinetic energy. This inadequacy can be attributed to the expression of $u'_f@p$, which
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Dimensionless deposition velocity $v_d^+$ with respect to the particle inertia $\tau_p^+$. The black-filled symbols are the experiments and blue symbols are the results with the stochastic wall deposition model. The simulations results using the normalized M.S. Langevin with a drift term (25) are shown for vertical (left) and horizontal (right) channels.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Dimensionless deposition velocity $v_d^+$ with respect to the particle inertia $\tau_p^+$ for the vertical channel. The black-filled symbols are the experiments. The simulations results using the M.S. Langevin in (9) are shown for vertical channel.</figcaption></figure>

employs the absolute value of a random variable following a Gaussian distribution. Furthermore, the observed PDF of sampled velocities deviates from the expected Gaussian distribution, exhibiting more of an exponential behavior. A generalized formulation of the M.S. Langevin model was implemented and gave good agreement.

Investigations into fluid element trajectories within a turbulent channel flow revealed significant limitations in the M.S. Langevin model, particularly its inability to accurately retrieve the fluid kinetic energy throughout the channel. Additionally, the model exhibited a notable 'spurious' drift phenomenon, deviating from expected physical behaviors. To address these deficiencies, a series of enhancements to the model have been implemented.

In the context of a wall-modeled approach, the deposition of low-inertia particles is consistently overestimated by the M.S. Langevin model.

Consequently, in order to address the limitations arising from the coarse description of the boundary layer, a Lagrangian stochastic wall deposition model was revisited, developed, and implemented to handle 3D complex geometries. The evaluation of this model demonstrated

promising results, particularly when it was activated and coupled with the M.S. Langevin model, showing favorable agreement between the computed deposition velocity and experimental observations.

In contrast to the models proposed by Guingo and Minier (2008) or Chibbaro and Minier (2008) whose simulations simulate the particle trajectories in the boundary layer by applying the Langevin modelling approach (originally developed for high Reynolds number flows) in the transition and viscous sublayers, the proposed method leads to apply Langevin model only in the fully developed turbulent region where they are legitimately applicable ($y_h^+ \gg 1$). The stochastic wall deposition model presented here is suitable for boundary layers in equilibrium. However, its validity may be questionable for impaction zones and regions with reverse pressure gradients such as detached boundary layers.

These findings strongly advocate for the inclusion of such deposition model in CFD codes, as they enable more accurate characterization of deposition in various geometries. By incorporating the coupling between the stochastic wall deposition model and other Langevin models (Arcen and Tanière, 2009; Minier et al., 2014), CFD simulations can provide improved predictions and a better understanding of particle behavior in turbulent flows, offering valuable insights across a wide range of practical applications.

### CRediT authorship contribution statement

**Cheikhna Talebmoustaph:** Writing – review & editing, Writing – original draft, Validation, Methodology, Conceptualization.
**Pascal Fede:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization.
**Olivier Simonin:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization.
**Maxime Pal-lud:** Funding acquisition.
**Priyank Maheshwari:** Writing – review & editing, Supervision, Funding acquisition.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.
Acknowledgments

Part of the financial support of this research project was provided by ANRT (Association Nationale de la Recherche et de la Technologie, Ministère de la Recherche) through the TotalEnergies/ CIFRE Contract No. 2020/1277.

References

Aguinaga, S., Simonin, O., Borée, J., Herbert, V., 2009. A Lagrangian stochastic model for droplet deposition simulations in connection with wall function approaches. In: Fluids Engineering Division Summer Meeting. Vol. 43727, pp. 795–805.

Arcen, B., Tanière, A., 2009. Simulation of a particle-laden turbulent channel flow using an improved stochastic Lagrangian model. Phys. Fluids 21 (4), 043303.

Arcen, B., Tanière, A., Oesterlé, B., 2005. Influence of the gravity field on the turbulence seen by heavy discrete particles in an inhomogeneous flow. In: Engineering Turbulence Modelling and Experiments 6. Elsevier, pp. 949–958.

Chibbaro, S., Minier, J.-P., 2008. Langevin PDF simulation of particle deposition in a turbulent pipe flow. J. Aerosol Sci. 39 (7), 555–571.

Csanady, G., 1963. Turbulent diffusion of heavy particles in the atmosphere. J. Atmos. Sci. 20 (3), 201–208.

Dehbi, A., 2008. A CFD model for particle dispersion in turbulent boundary layer flows. J. Nucl. Eng. Des. 238, 707–715. http://dx.doi.org/10.1016/j.nucengdes.2007.02.055.

Dreeben, T.D., Pope, S., 1997. Probability density function and Reynolds-stress modeling of near-wall turbulent flows. J. Phys. Fluids 9, 154–163.

ESI Group, 2020. OpenFOAM version 2006. https://www.openfoam.com/releases/openfoam-v2006/.

Fluent, 2006. FLUENT 6.3 User's Guide. Fluent Inc., Lebanon, NH, USA.

Gao, N., Niu, J., He, Q., Zhu, T., Wu, J., 2012. Using RANS turbulence models and Lagrangian approach to predict particle deposition in turbulent channel flows. J. Build. Environ. 48, 206–214. http://dx.doi.org/10.1016/j.buildenv.2011.09.003.

Gosman, A., Ioannides, E., 1981. Aspects of computer simulation of liquid fuelled combustors. In: AIAA 19th Aerospace Sci. Meeting Paper AIAA-81-0323. St. Louis, MO.

Guingo, M., Minier, J.-P., 2008. A stochastic model of coherent structures for particle deposition in turbulent flows. Phys. Fluids 20 (5), 053303.

Hinze, J., 1975. Turbulence, second ed. McGraw-Hill, New York.

Huilier, D.G., 2021. An overview of the Lagrangian dispersion modeling of heavy particles in homogeneous isotropic turbulence and considerations on related LES simulations. J. Fluids 6, 145. http://dx.doi.org/10.3390/fluids6040145.

Hutchinson, P., Hewitt, G., Dukler, A., 1971. Deposition of liquid or solid dispersions from turbulent gas streams: a stochastic model. Chem. Eng. Sci. 26 (3), 419–439.

Innocenti, A., Fox, R.O., Salvetti, M.V., Chibbaro, S., 2019. A Lagrangian probability-density-function model for collisional turbulent fluid-particle flows. J. Fluid Mech. 862, 449–489.

Kallio, G., Reeks, M., 1989. A numerical simulation of particle deposition in turbulent boundary layers. Int. J. Multiph. Flow 15, 433–446.

Kvasnak, W., Ahmadi, G., Bayer, R., Gaynes, M., 1993. Experimental investigation of dust particle deposition in a turbulent channel flow. J. Aerosol Sci. 24, 795–815.

Liu, B.Y.H., Agarwal, J.K., 1974. Experimental observation of aerosol deposition in turbulent flow. J. Aerosol Sci. 5, 145–155.

MacInnes, J., Bracco, F., 1992. Stochastic particle dispersion modeling and the tracer-particle limit. J. Phys. Fluids 4 (12), 2809–2824.

Matida, E.A., Nishino, K., Torii, K., 2000. Statistical simulation of particle deposition on the wall from turbulent dispersed pipe flow. Int. J. Heat Fluid Flow 21, 389–402, URL: www.elsevier.com/locate/ijh.

Minier, J.-P., 2015. On Lagrangian stochastic methods for turbulent polydisperse two-phase reactive flows. Prog. Energy Combust. Sci. 50, 1–62.

Minier, J.-P., Chibbaro, S., Pope, S.B., 2014. Guidelines for the formulation of Lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent flows. Phys. Fluids 26 (11).

Montgomery, T.L., Corn, M., 1970. Aerosol deposition in a pipe with turbulent airflow. J. Aerosol Sci. I, 185–213.

Nerisson, P., Simonin, O., Ricciardi, L., Douce, A., Fazileabasse, J., 2011. Improved CFD transport and boundary conditions models for low-inertia particles. Comput. & Fluids 40 (1), 79–91.

Obukhov, A., 1959. Description of turbulence in terms of Lagrangian variables. Adv. Geophys. 6, 113–116.

Papavergos, P., Hedley, A., 1984. Particle deposition behaviour from turbulent flows. Chem. Eng. Res. Des. 62, 275–295.

Parker, S., Foat, T., Preston, S., 2008. Towards quantitative prediction of aerosol deposition from turbulent flows. J. Aerosol Sci. 39, 99–112. http://dx.doi.org/10.1016/j.jaerosci.2007.10.002.

Pope, S.B., 1994. Lagrangian PDF methods for turbulent flows. Annu. Rev. Fluid Mech. 26 (1), 23–63.

Sehmel, G.A., 1973. Eddy diffusivities and deposition for isothermal flow and smooth surfaces. J. Aerosol Sci. 4, 138.

Simonin, O., Deutsch, E., Minier, J., 1993. Eulerian prediction of the fluid/particle correlated motion in turbulent two-phase flows. Appl. Sci. Res. 51, 275–283.

Sippola, M.R., Nazaroff, W.W., 2002. Particle deposition from turbulent flow: Review of published research and its applicability to ventilation ducts in commercial buildings.

Sommerfeld, M., Kohnen, G., Rueger, M., 1993. Some open questions and inconsistencies of Lagrangian particle dispersion models. In: Proceedings of the 8th Symposium on Turbulent Shear Flow. Kyoto, Japan.

Wood, N., 1981. A simple method for the calculation of turbulent deposition to smooth and rough surfaces. J. Aerosol Sci. 12 (3), 275–290.