Planar shear flow effects on particle dispersion over a normal flat plate

Abstract

We investigate the effect of planar shear inflow on the dispersion of inertial particles in the unsteady laminar wake behind a normal flat plate and the corresponding particle deposition on the plate. Eulerian–Lagrangian simulations are performed to solve the dilute particulate suspension flow over the plate. To explore the effect of shear, we perform a comparative assessment between two different cases, distinguished based on their inflow profiles, namely uniform and linear shear flow (with shear parameter K = 0.1). Flow Reynolds number based on the mean inflow velocity is 60 in each case, where the unsteady wake is characterized by von Kármán vortex street. Three different classes of particles are tracked in each case, where each particle class is characterized by its Stokes number ($\mathrm{St}$), namely $\mathrm{St}=0.1$ (low), $\mathrm{St}=1$ (medium), $\mathrm{St}=10$ (high). Flow visualizations show that planar shear alters the vortex shedding pattern, especially at far downstream positions. Such changes in the wake flow lead to significant streamwise stretching and rotation of particle-depleted void regions and adjacent particle clusters. With the help of Voronoï diagrams, we characterize fluid shear driven changes in inhomogeneous particle distribution patterns and lateral shifts in positions of particle clusters observed for each $\mathrm{St}$-class of particles. The quantitative effect of shear on the dispersed phase statistics of each $\mathrm{St}$-class of particles is further analyzed. Such effects are not appreciable near the plate, but become evident further downstream away from the plate. At such positions, the impact of shear on the concentration of medium $\mathrm{St}$-class particles is the most significant among all classes, where the mean particle concentration near the wake centerline is significantly increased. Furthermore, local peaks in particle concentration profiles and minima in particle velocity profiles shift laterally in their positions toward the low-stream region due to the shear. The shear also causes a marginal reduction in particle deposition efficiency.

Appendices

Appendix A: LaParT details

See Fig. 11.

Fig. 11

Algorithm for the proposed Eulerian–Lagrangian solver MGLET-LaParT

1.1 1. Identifying the nearest stencil (neighbor search)

The first step in the algorithm is to identify the appropriate stencil (containing the velocity component data and their corresponding positions) encapsulating a particle for each velocity component. We again emphasize that the positions used to compute and store the velocity component field (i.e., $u_\mathrm{f}$, $v_\mathrm{f}$, and $w_\mathrm{f}$) are staggered in relation to those for the pressure field (p) in the respective directions (i.e., x, y, and z). The sequential search method is used to identify the host stencil.

1.2 2. Three-dimensional (3D) interpolation

We obtain the fluid velocity (${\varvec{u}}_\mathrm{f}({\varvec{x}}_\mathrm{p},t)$) at the particle position (${\varvec{x}}_\mathrm{p}$) by interpolating the velocity field data available at nodes of the identified stencil. The choice of the interpolation scheme is important since it affects the accuracy of numerically computed particle variables (e.g., particle trajectory) and the computation cost. We expect that the interpolation scheme should be at least second-order accurate.⁷ The scheme should be easy to implement and affordable in terms of computational cost. Thus, we have implemented the trilinear interpolation based on $2^3$ points [42] in the LaParT.

Fig. 12

Stencil for trilinear interpolation

The interpolated value $g^{TL}(x_\mathrm{p},y_\mathrm{p},z_\mathrm{p})$ is constructed from scalar values (i.e., $u_\mathrm{f}$, $v_\mathrm{f}$, and $w_\mathrm{f}$) located at $2^3$ grid nodes at any interior position $(x_\mathrm{p},y_\mathrm{p},z_\mathrm{p})$ within a rectangular 3D stencil. This stencil is formed by a bounding pair of constant x-, y-, and z-planes in each direction as shown in Fig. 12. The interpolated value $g^{TL}(x_\mathrm{p},y_\mathrm{p},z_\mathrm{p})$ by making use of the trilinear scheme is given as

$$\begin{aligned} g^{TL}(x_\mathrm{p},y_\mathrm{p},z_\mathrm{p}) = \sum _{i = 1}^8 W_if_i \end{aligned}$$

(8)

where the index of the corner point i varies from 1 to 8, $f_i$ is the scalar value at the corner point, and $W_i$ is its corresponding weight used in the interpolation (shown in Fig. 12). This scheme is computationally inexpensive and second-order accurate.

Here, we bring to the reader’s attention that the velocity component is interpolated independently of the other components for any interpolation scheme. In other words, the incompressibility of the velocity field is not enforced in the interpolation scheme. This non-conservative nature of the interpolation can induce the non-physical clustering of particles due to numerical artifacts (see [16, 43]). Nevertheless, we use a fine grid resolution to minimize such effects.

1.3 3. Numerical integration

An explicit second-order Adams–Bashforth scheme is used to integrate Eq. (5) numerically. Particle velocity at the new, i.e., ${(n+1)}{^{th}}$ time step by making use of Adams–Bashforth scheme is given using Eq. (9). Here, the solution at the first time step is computed using the first-order Euler’s scheme,

$$\begin{aligned} {\varvec{u}}_{\mathrm {p}}^{{n+1}} = {\varvec{u}}_\mathrm {p}^{n} + \frac{\gamma \Delta t_\mathrm {p}}{\tau _\mathrm {p}}\left[ \frac{3}{2}({\varvec{u}}_\mathrm {f}^{n}-{\varvec{u}}_\mathrm {p}^{n})-\frac{1}{2} ({\varvec{u}}_\mathrm {f}^{n-1}-{\varvec{u}}_\mathrm {p}^{n-1})\right] . \end{aligned}$$

(9)

The particle position is updated according to Eq. (10) as

$$\begin{aligned} {\varvec{x}}_\mathrm{p}^{n+1}={\varvec{x}}_\mathrm{p}^{n}+{\Delta t_\mathrm{p}}\left[ \frac{3}{2}{\varvec{u}}_\mathrm{p}^{n+1}-\frac{1}{2}{\varvec{u}}_\mathrm{p}^{n} \right] . \end{aligned}$$

(10)

The scheme is conditionally stable. We can not exactly estimate the $\Delta t_\mathrm{p}$ value to ensure that the numerical integration of Eqs. (5) and (6) is stable since these equations contain a term ${\varvec{u}}_\mathrm{f}({\varvec{x}}_\mathrm{p},t)$, which is not known apriori (see [44]). Jacobs et al. [45] mention that an approximate value of $\Delta t_\mathrm{p}$ (for stable numerical integration) can be given by

$$\begin{aligned} \Delta t_\mathrm{p} < \sigma _\mathrm{p} \tau _\mathrm{p} \end{aligned}$$

(11)

where the constant $\sigma _\mathrm{p}$ is equal to unity. Using Eq. (11), the value of $\Delta t_\mathrm{p}$ can be selected according to the $\tau _\mathrm{p}$ value. However, as noted by Kontomaris et al. [46], $\Delta t_\mathrm{f}$ required to satisfy the CFL criterion to solve the carrier phase equations is lower than $\Delta t_\mathrm{p}$ required for the stable computation of the dispersed phase equations. We use the same time step size value for both the carrier and dispersed phase (i.e., $\Delta t_\mathrm{f} = \Delta t_\mathrm{p}$) in turbulent flow simulations. This conservative approach leads to a $\Delta t_\mathrm{p}$ value which is at least one order of magnitude smaller than the particle relaxation time $\tau _\mathrm{p}$ and, thus, sufficiently small for the stable computation of the dispersed phase equations.

Regarding the boundary condition, if any particle approaches the body or a domain boundary, then an appropriate boundary condition is applied. For instance, if any particle crosses the boundary of the computation domain, then it would not be tracked further for the computation.

Appendix B: Verification and validation

1.1 1. Verification of LaParT

The verification of the numerical accuracy of schemes implemented in the Lagrangian solver using canonical numerical tests is discussed in test cases 1.1 and 1.2. The Lagrangian solver is isolated from the Eulerian solver, and numerical tests are performed on the former. The carrier phase velocity field (described using analytical expressions) is directly assigned to the Eulerian grid nodes without solving the Navier–Stokes equations. The motion of a particle in such tests can be determined analytically. Therefore, a direct comparison between numerical and analytical solutions can be made to quantify the numerical errors associated with any scheme. In the verification test 1.3, the carrier and dispersed phases are solved using MGLET-LaParT.

1.1.1 1.1. Particle in a steady sinusoidal flow

We track an inertia-free particle in a steady 2D sinusoidal flow (in xz-plane). This test problem is the same as that used by Marchioli et al. [27] to verify the order of accuracy of the interpolation based on the Taylor-series expansion. Prior to them, Kontomaris et al. [46] also had used the present problem for testing various interpolation schemes. The purpose of this test is to quantify the errors associated with an interpolation scheme and confirm its order of accuracy. The sinusoidal velocity field in this problem is described by $u_\mathrm{f} = A_1$ and $w_\mathrm{f} = A_2 sin(k_x x)$ (in m/s) where $A_1 = A_2 = 0.025$ are constants (in m/s) and $k_x = \pi$ is the wavenumber. The particle is located at $x_\mathrm{p}(0) = -0.8$ and $z_\mathrm{p}(0) = 0.75$ (both in m) in the beginning. The particle motion is tracked in the time interval of $0{-}100$ s. An equidistant grid is used. The time step size used for the computation is $\Delta t_\mathrm{p} = 0.01$ s to ensure that the numerical integration is stable and time stepping errors are small. A comparison of particle trajectory⁸ during the period of $0-100$  s for the interpolation scheme shows a good agreement (in Fig. 13a) with that determined analytically using expressions given in Marchioli et al. [27]. Let $w_{fn,i}$ and $w_{fa,i}$ be the numerically computed and analytically determined value of the undisturbed fluid velocity at particle position at $i{^{th}}$ time step, then the root-mean-square (rms) error in the interpolation of the fluid velocity⁹$\epsilon _{\mathrm{rms}, wf}$ over N number of time steps is given as

$$\begin{aligned} \epsilon _{\mathrm{rms}, wf} = \sqrt{\frac{\sum _{i = 1}^N (|w_{fn,i} - w_{fa,i}|)^2}{N}}. \end{aligned}$$

(12)

Furthermore, we define a grid resolution parameter R to describe the characteristic length scale of the flow $\lambda$ in comparison to the grid cell thickness $\Delta x$ as $R = {\lambda }/{\Delta x}$, where $\lambda = 2\pi /k_x$. We consider four different grids: coarse $R = 8$, coarse-medium $R = 16$, medium $R = 32$, and fine $R = 64$. As shown in Fig. 13b, the rms error in the fluid velocity interpolation $\epsilon _{\mathrm{rms}, wf}$ decreases as the grid resolution improves. We confirm that the trilinear interpolation scheme implemented is second-order accurate (i.e., $\epsilon _{\mathrm{rms}, wf}$ decreases at the rate $R^{-2}$).

Fig. 13

A single inertia-free particle in a steady, spatially varying sinusoidal flow. a Comparison of particle paths computed from the simulation against those using analytical expressions, b verification of the order of the fluid velocity interpolation scheme

1.1.2 1.2. Particle in a spatially uniform but temporally varying flow

This test comprises tracking the motion of an inertial particle in a spatially uniform but temporally varying velocity field in a sinusoidal fashion, as done in Jacobs et al. [45]. The flow velocity field is given by $u_\mathrm{f} = sin(k_t t)$ where $k_t = 16\pi$. The particle is displaced only in the x-direction. We set the particle initial position $x_\mathrm{p}(0)$ and velocity $u_\mathrm{p}(0)$ to zero. The particle is tracked during the time window between 0 and 1 s. Five different values of the particle Stokes number $\mathrm{St}$ = 0.1, 0.5, 1, 5, and 10 are considered. Here, the particle relaxation time $\tau _\mathrm{p}$ is normalized using the characteristic flow time scale $\tau _\mathrm{f}$ to obtain $\mathrm{St}$, where $\tau _\mathrm{f}$ = 1 s in this problem. The numerically computed instantaneous particle positions $x_\mathrm{p}(t)$ match those computed using analytical expressions, as shown in Fig. 14a. (We refer to Jacobs et al. [45] for the analytical expressions of $x_\mathrm{p}(t)$.) This match explains that the effect of $\tau _\mathrm{p}$ on the particle trajectory is captured accurately in the numerical integration. We compute the rms error ($\epsilon _{\mathrm{rms}}$) in the calculation of the particle velocity ($u_\mathrm{p}$) and position ($x_\mathrm{p}$) for different time step sizes $\Delta t_\mathrm{p}=10^{-2}$, $10^{-3}$, $10^{-4}$, and $10^{-5}$ for a particle with $\mathrm{St} = 1$, as shown in Fig. 14b. It is noticed that $\epsilon _{\mathrm{rms}}$ decreases with a decrease in $\Delta t_\mathrm{p}$, implying the implemented numerical scheme is convergent. The numerical solutions of particle position and momentum are first- and second-order accurate, respectively.

Fig. 14

A single inertial particle in a uniform, transient flow. a Comparison of particle trajectory computed from the simulation against those using analytical solution, b global error ($\epsilon _{\mathrm{rms}}$) in particle velocity ($u_\mathrm{p}$) and displacement ($x_\mathrm{p}$) versus time step size ($\Delta t_\mathrm{p}$)

1.1.3 1.3 Low inertia particles in the unsteady non-uniform flow

The hydrodynamic response of very low inertia particles to the non-uniform flow is expected to be similar to fluid particles (tracers). Therefore, such tracer-like particles should be distributed uniformly in the unsteady non-uniform flow. The present test verifies that the spatial distribution of tracer-like particles obtained from the numerical solution is consistent with the expected physical behavior. The confirmation would assure that the schemes implemented in LaParT respect the incompressibility constraint. To this end, the unsteady laminar wake behind a normal flat plate is considered a representative non-uniform test flow case. The Reynolds number for the flow based on the uniform inflow velocity ($U_\mathrm{c}$) and plate height (h) is 60, where the unsteady wake behind the plate is characterized with a von Kármán vortex street. Particles with $\mathrm{St}$ = 0.005 (representing significantly low particle inertia) are injected homogeneously near the inlet. A snapshot of the instantaneous distribution of particles in the computation domain taken at $t = 40$ $h/U_\mathrm{c}$ (from the beginning of the particle injection) is shown in Fig. 15. We confirm that the instantaneous spatial distribution of such tracer-like particles is homogeneous throughout the domain, which agrees with the expected physical behavior.

Fig. 15

Distribution of particles with very low inertia (i.e., $\mathrm{St}$ = 0.005) in the unsteady laminar flow over a normal flat plate

1.2 2. Validation of MGLET-LaParT

We discuss the validation of MGLET-LaParT using the benchmark particle-laden turbulent flow case of Marchioli et al. [47], where the results of five research groups using different numerical procedures have been documented. We consider the results of TUD and UUD groups for the comparison purpose. The flow configuration and test setup used in our simulation are the same as those of the benchmark validation test and summed up briefly again for the reader’s convenience. The validation test investigates the pressure-driven turbulent channel flow of a dilute suspension that comprises many monodispersed, small-sized inertial particles. The friction Reynolds number of the flow, based on the half-channel height h, and friction velocity $u_*$, is $\mathrm{Re}_* = \rho _\mathrm{f} u_* h/\mu _\mathrm{f}$ = 150. Here, the friction velocity is obtained using $\rho _\mathrm{f}$ and the mean wall shear stress $\tau _\mathrm{w}$ as $u_* = \sqrt{\tau _\mathrm{w}/\rho _\mathrm{f}}$. The channel is positioned such that the Cartesian co-ordinate¹⁰ axes x, y, and z are aligned in the streamwise, spanwise, and wall-normal directions, respectively. Hereafter, simulation parameters given in the wall-units will be denoted using symbols accompanied by a superscript $+$. Such simulation parameters are scaled using the viscous length $l_\mathrm{v} = \nu _\mathrm{f}/u_*$, time $\tau _\mathrm{v} = \nu _\mathrm{f}/u_*^2$, and velocity scales $u_*$. Periodic boundary conditions are imposed in the homogeneous x- and y-directions. The smooth flat walls positioned at $z = 0$ and $z = 2h$ in the inhomogeneous z-direction are treated by imposing the no-slip and non-penetration condition on the carrier fluid. The channel dimensions are $L_x=12.56h$ ($L_x^+=1885$), $L_y=6.28h$ ($L_y^+= 942.5$), and $L_z=2h$ ($L_z^+=300$) in the $x-$, $y-$, and $z-$ directions, respectively. The constant grid-cell dimensions in the $x-$ and $y-$ directions are $\Delta x^+=9.82$ and $\Delta y^+=4.91$, respectively. The grid is stretched in the $z-$ direction, where $\Delta z^+_{\min }=0.36$, and $\Delta z^+_{\max }=4.06$ are, respectively, the minimum (near the wall) and maximum grid-cell dimensions (channel center).

The dispersed phase comprises small-sized, denser-than-fluid, spherical particles. We consider three different cases, each with a unique particle Stokes number $\mathrm{St}_\mathrm{v}$ (i.e., 1, 5, and 25), where $\mathrm{St}_\mathrm{v}$ is given as $\mathrm{St}_\mathrm{v} = \tau _\mathrm{p}/\tau _\mathrm{v}$. Particle properties in each case are documented in Table 5. We have homogeneously injected a total of $10^5$ particles (in each $\mathrm{St}_\mathrm{v}$ case) into the statistically stationary turbulent carrier flow. Particles are randomly positioned in x- and y-directions with a uniform spacing in the z-direction. The periodic domain is considered for a particle in the streamwise (x) and spanwise (y) direction. If a particle approaches the smooth flat wall such that the distance between its geometric center and the wall is less than its radius, then this condition is considered as the particle-wall collision event. Such a particle is rebounded back elastically in the domain after the collision. In Fig. 16, a comparison of root-mean-square (rms) fluctuations in the streamwise ($u_{i, \mathrm{rms}}^+$), spanwise ($v_{i, \mathrm{rms}}^+$), and wall-normal ($w_{i, \mathrm{rms}}^+$) velocity components of our simulation with benchmark groups is shown for the fluid ($i = f$; see Fig. 16a) and particle phases ($i = p$; see Fig. 16b–d). With an increase in particle inertia from $\mathrm{St}_\mathrm{v}$ = 1 to 25, $u_{p, \mathrm{rms}}^+$ values of benchmark results in the near-wall region become larger compared to $u_{f, \mathrm{rms}}^+$. On the other hand, $v_{p, \mathrm{rms}}^+$ and $w_{p, \mathrm{rms}}^+$ values get smaller in magnitude than $v_{f, \mathrm{rms}}^+$ and $w_{f, \mathrm{rms}}^+$, respectively, with an increase in $\mathrm{St}_\mathrm{v}$. This trend is captured in our simulations. Furthermore, a quantitative comparison of profiles for $\mathrm{St}_\mathrm{v}$ = 1 (Fig. 16b), $\mathrm{St}_\mathrm{v}$ = 5 (Fig. 16c), and $\mathrm{St}_\mathrm{v}$ = 25 (Fig. 16d) suggests that our simulation results are in a good agreement with the benchmark data.

Table 5 Particle properties

Fig. 16

Root-mean-square of velocity fluctuations. a fluid, and particles with b $\mathrm{St}_\mathrm{v}$ = 1, c $\mathrm{St}_\mathrm{v}$ = 5, d $\mathrm{St}_\mathrm{v}$ = 25