Particle behavior in a turbulent flow within an axially corrugated geometry

Abstract

In this numerical study particle behavior inside a sinusoidal pipe geometry is analyzed. The 3D geometry consists of three identical modules, with a periodic boundary condition applied to the flow in the stream wise direction. The incompressible, turbulent gas flow is modeled using a Large Eddy Simulation (LES) approach. Furthermore, the particle dynamics are simulated using a Lagrangian point force approach incorporating the Stokes drag and slip correction factor. Four different sizes of particles, corresponding to a Stokes number less than unity, are considered along with two different inflow conditions: continuous and pulsatile. The pulsatile inflow has an associated flow frequency of 80 Hz. The fluid flow through the sinusoidal pipe is characterized by weak flow separation in the expansion zones of the sinusoidal pipe geometry, where induced shear layers and weak recirculation zones are identified. Particle behavior under the two inflow conditions is quantified using particle dispersion, particle residence time, and average radial position of the particle. No discernible difference in the particle behavior is observed between the two inflow conditions. As the observed recirculation zones are weak, the particles are not retained within the cavities for a long duration of time, thereby reducing their likelihood of agglomerating.

Introduction

Corrugated geometries have been studied extensively with regards to their effect on fluid mixing,^1–3 mass transfer,^4,5 and heat transfer.^6–10 By corrugated geometries we are referring to geometries with successive expansions and contractions. These geometric configurations have sometimes been referred in the literature as wavy-walled, furrowed, and sinusoidal channels as well. More broadly, it has been reported that both the mass and heat transfer characteristics of such geometrical configurations were better than their straight channel counterparts based on the Sherwood and Nusselt number calculated for different Reynolds number. This has been credited to the intrinsic flow structures found within such corrugated geometries which contribute toward better fluid mixing. However, a larger pressure drop has also been observed for corrugated geometries as compared to straight channel geometries. ⁷

Studies primarily focusing on understanding the flow behavior in corrugated geometries can be found in the past literature as well. Guzmàn and Amon^11,12 have reported that as the Reynolds number was increased, the laminar flow underwent successive bifurcations to transition and eventually arrive at a chaotic regime. Various studies based on linear stability analysis have been conducted to identify modes leading to flow instability inside corrugated geometries. The presence of a centrifugal instability has been identified both experimentally^13,14 and numerically. ¹⁵ Furthermore, a Kelvin-Helmholtz like traveling wave instability has also been found within corrugated pipes, which develops at the margin of the recirculation zones and the induced shear layer.^16–18

Studies focusing on the distribution of solid particles in such corrugated pipe like geometries have been reported in the literature as well. A particular case would be the so called wavy-walled geometry, consisting of wavy bottom wall and a flat top wall. The flow over wavy-walled geometries has been studied both experimentally^19,20 and numerically.^21–23 Studying particle distribution with such a geometric configuration can be used to approximate sediment transport. In the Direct Numerical Simulation (DNS) study by Boersma, ²⁴ particles with relaxation time on the order of $O({10}^{-2})$ were tracked in a wavy-wall configuration using a Lagrangian approach. It was reported that the maximum particle concentration was found slightly downstream of the peak of the wavy wall. This was credited to the high wall shear stress fluctuations observed upstream of the wave peak which were redirecting the particles back into the flow. Chang and Scotti ²⁵ further conducted a study analyzing the influence of coherent flow structures on the particle distribution in a wavy wall geometry. A Large Eddy Simulation (LES) was used to model the flow and the particles were tracked using a Lagrangian approach. The particle Stokes number was in the range of $O({10}^{-7}-{10}^{-4})$ . It was observed that the particle behavior depended on their initial location. Particles initially located in the upslope region followed the wall until they reached the wall crest. Subsequently, they dispersed into the outer flow downstream of the wall crest. Particles located close to the wall crest however, moved more slowly and were eventually pushed into the furrows. This difference in particle behavior was further linked with the difference in the coherent flow structures observed in the up-slope and down-slope region of the wavy wall.

Marchioli et al. ²⁶ further analyzed particle behavior inside a wavy wall configuration. Their setup (while similar in terms of methodology) differed from with Chang and Scotti ²⁵ in the type of fluid, the total number of particles in the domain, the size of the particles, and the ratio of the particle to fluid density utilized. The study by Marchioli et al. ²⁶ highlighted that there were two distinct flow regions affecting the particle behavior. The first region extended from the point of flow separation until the reattachment of the flow. Similarly, the second region consisted of the region downstream of the point of flow reattachment until the subsequent point of flow separation. It was reported that in the first region the particles were entrained by the recirculation zone and eventually deposited. In the second region the dominant coherent flow structures expelled the deposited particles toward the outer flow. Furthermore, it was observed that particles with more inertia tended to cluster around regions of high strain and low vorticity near the wall. A DNS study was performed by Lee et al. ²⁷ further analyzing particle behavior within various types of wavy wall configurations and a flat wall configuration. Particles with larger inertia had the tendency to cluster around the up-slope region of the wavy wall. This was again due to the dominant flow structures found in the up-slope region. The particle clustering close to the wall shifted from a strip pattern to a mixture of a strip and streak pattern as the ratio between the corrugation amplitude and the wall wavelength was successively reduced. A LES study performed by Chang and Park ²⁸ looked into the behavior of sediments in a sinusoidally shaped geometry under unsteady flow conditions. It was reported that particle dispersion was influenced by the presence of the coherent flow structures. In particular, during the deceleration phase of the flow, vortices were generated downstream of the crest of the wavy geometry that contributed to the convection of the sediments. Fonias and Grigoriadis ²⁹ have also considered a particle laden flow over a dune-like geometry using a LES approach for the flow and a four-way coupled Lagrangian approach for the particles. The main conclusions were that the presence of particles attenuated the mean flow velocity and the turbulent kinetic energy. Particles also tended to cluster around the upslope region of the dune as opposed to where the recirculation zones were found.

Scope of this work

Most of the previous work based on corrugated geometries have focused primarily on heat and mass transfer. The behavior of particles within such geometries has mostly been considered in geometries consisting of a lower wavy wall and a flat upper wall. Past studies^30,31 have explored the possibility of particle agglomeration using such corrugated geometries. These geometries consist of successive expansions and contractions meant to manipulate the flow field. If the particles being carried by the flow do not have significant inertia then a suitable flow manipulation may lead to enhanced interaction between the particles. Particle agglomeration may potentially aid with improving the efficiency of exhaust after treatment systems, specifically particulate filters, in heavy duty vehicles.

The current study considers a $3D$ idealized sinusoidal pipe like geometry consisting of expansion and contraction regions (see Figure 1). The relevant geometrical parameters associated with the geometry, such as geometrical wavelength and cross sectional diameter, are based on past studies.^31,32 Two different inflow conditions are considered. The Reynolds number for both inflow scenarios is $O({10}^5)$ . The primary aim of this study is to gain insight into the behavior of particles in such a geometry under varying inflow conditions. Firstly, the flow field found under different inflow conditions is presented. This is followed by quantifying the observed residence time and radial position for particles of different sizes under different inflow conditions.

Figure 1.

Illustration of the wavy pipe geometry with annotations highlighting the important geometrical dimensions. Periodic flow boundary conditions are applied in the streamwise direction.

Geometry and methodology

Geometry and mesh

In Figure 1 the wavy pipe geometry considered for this numerical study is presented. The minimum and maximum diameter of this geometry is $D_{\min }=0.05m$ and $D_{\max }=0.085m$ , respectively. Furthermore, the geometric wavelength, that is, the length between two successive minimum cross sectional planes is $\lambda =0.12m$ . It is indicated in Figure 1, that the geometry in its entirety consists of three wavelengths. The wavy wall is given by the expression:

$$y(x)=r_{\mathrm{mean}}-a\cos \left(\frac{2\pi x}{\lambda }\right),$$

(1)

where $r_{\mathrm{mean}}=0.5(r_{\max }+r_{\min })$ is the mean hydraulic radius with $r_{\min }$ and $r_{\max }$ denoting the minimum and maximum radius, respectively. The final term in equation (1) contains the corrugation amplitude $a=0.5(r_{\max }-r_{\min })$ . The main properties of the different grid resolutions considered for this numerical study are provided in Table 1. Three different grids are considered for the grid resolution study. For all three grids, the surface averaged non dimensional wall normal distance $y^{+}<1$ , where $y^{+}=y{u}_{\tau }/\nu$ with $u_{\tau }$ and $\nu$ denoting the friction velocity and kinematic viscosity, respectively. For all three grids the average cell spacing is calculated based on the formulation ³³ :

$$\mathrm{\Delta }x={\left(\frac{1}{N}\sum _{i=1}^N\mathrm{\Delta }{V}_i\right)}^{\frac{1}{3}},$$

(2)

where $N$ is the total number of cells in the computational domain and $V_i$ represents the volume of the $i$ ’th cell. This is provided in the second column of Table 1. The third and fourth columns of Table 1 denote the non-dimensional cell spacing in the axial and circumferential directions, respectively. They were estimated using the surface averaged friction velocity. In the final column of Table 1 the total number of cells in the computational domain are stated for each individual grid. The estimated Taylor micro length scale for this case is ${\lambda }_g\approx 1.3·{10}^{-3}m$ . Furthermore, based on Piomelli and Chasnov, ³⁴ the non dimensional cell spacing highlighted for the computational meshes $G2$ and $G3$ is sufficiently resolved only for a wall modeled LES.

Table 1.

Summary of the different grid resolutions considered.

Grid	$\mathrm{\Delta }x(\times {10}^{-3}\mathrm{m})$	$\mathrm{\Delta }{x}^{+}$	$\mathrm{\Delta }{\Theta }^{+}$	$N(\times {10}^6)$
$G1$	$2.2$	$644$	$408$	$0.1$
$G2$	$1.4$	$410$	$250$	$0.4$
$G3$	$0.9$	$270$	$170$	$1.4$

Numerical methods

The numerical simulations are performed using the finite volume based solver STAR-CCM+. ³⁵ The 3D incompressible continuity and Navier–Stokes equation are solved using a Large Eddy Simulation (LES) approach given by the following equations:

$$\nabla ·\mathbf{u}=0,$$

(3)

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}+\nabla ·{\tau }_{\mathrm{SGS}},$$

(4)

where the gas density $\rho$ is constant and the variables $\mathbf{u}$ and $p$ are the filtered velocity and pressure, respectively. A subgrid scale model is needed in order to evaluate the subgrid scale stress tensor ${\tau }_{\mathrm{SGS}}$ . A wall adapting local eddy viscosity (WALE) ³⁶ subgrid scale model has been applied in the current numerical setup. A bounded central differencing scheme is used as a convection scheme. This scheme has been utilized previously in past studies^37,38 utilizing the LES approach in STAR-CCM+. An implicit temporal scheme is used with a fixed time step $\mathrm{\Delta }t={10}^{-6}s$ which ensures a Courant number less than unity. In this study two different flow field cases are considered. For both cases periodic boundary conditions are applied in the streamwise direction by specifying a mass flow rate. In the first case a constant mass flow rate is imposed at the inlet plane. The corresponding bulk velocity at the inlet plane is denoted as $U_b$ Hereinafter, this case would be referred to as the continuous mass flow case.

For the second case a pulsatile mass flow rate based on a sinusoidal pulse form is applied at the inlet plane. The pulse form is characterized by a mean velocity $U_0$ , amplitude $A=0.125U_0$ and frequency $f$ , respectively. At the solid walls of the geometry no-slip boundary conditions are imposed. The wall treatment in the numerical approach utilizes a blending function based on Reichardt’s law. ³⁹

The Reynolds number for both flow cases based on the minimum diameter and the bulk velocity at the inlet is approximately $2.5\times {10}^5$ . Furthermore, the Womersley number associated with the pulsatile flow case is $283$ . Thus, the flow is highly turbulent and viscous effects carry significance only in the boundary layer.

Data reduction of the flow field data is achieved using time averaging and phase locked averaging. For a given flow variable $\psi (x,y,z,t)$ , the associated time averaged variable is calculated as follows:

$$\overline{\psi }(x,y,z,t)=\frac{1}{T}{\int }_{t_0}^{t_0+T}\psi (x,y,z,t)\mathrm{dt},$$

(5)

where the time $t_0$ denotes the instance at which the averaging procedure starts and $T$ denotes the total time interval for which the averaging procedure is performed. The fluctuating part of the flow variable $\psi \prime$ is the difference between the instantaneous value and the time averaged value associated with $\psi$ . Based on the fluctuations of the velocity field the turbulent kinetic energy (TKE) is given by $k=0.5\overline{u_i^{\prime }{u\prime }_i}$ , where the index notation is based on the Einstein convention. In an analogous manner the phase locked averaging can be calculated for a flow variable $\psi (x,y,z,t)$ as:

$$\begin{array}{c}<\psi (x,y,z,t_{\varphi })>=\frac{1}{N}\sum _{n=0}^{N-1}\psi (x,y,z,t_{\varphi }+\mathrm{nP}),\hfill \\ \multicolumn{1}{c}{\varphi =0,\dots ,P-1,}\end{array}$$

(6)

where $N$ is the total number of samples of the selected phase, $P$ is the total number of samples taken per period and the phase $\varphi$ is the specific instance of time such that $t=t_{\varphi }$ . The difference between the instantaneous value and the phase locked averaged quantity results in the fluctuating part $\hat{\psi }$ associated with the phase under consideration. For a given phase the turbulent kinetic energy can be calculated using the fluctuations of the velocity field as $k_p=0.5<{\hat{u}}_i{\hat{u}}_i>$ . The normalized pulse form applied as an inflow condition in the current setup is displayed in Figure 2.

Figure 2.

Normalized sinusoidal pulse form imposed for the pulsatile inflow scenario.

In the forthcoming results section the phase locked averaged data is presented for the phases indicated in Figure 2. The collection of statistics started after 11 flow through time, where the flow through time is given by $3\lambda /U_b$ . In the continuous mass flow case these statistics are collected for at least $173$ flow through time. For the pulsatile mass flow the statistics are collected for $40$ pulse cycles.

A Lagrangian point particle treatment is used to track the solid particles as described in Crowe. ⁴⁰ Four different sizes of the particles are considered: $2.5\mu m$ , $800\mathrm{nm}$ , $500\mathrm{nm}$ , and $180\mathrm{nm}$ . The corresponding Stokes number varies from $O({10}^{-5}-{10}^{-3})$ . Particle density is $1000\mathrm{kg}/m^3$ which gives a fluid density to particle density ratio significantly smaller than unity. For such a fluid density to particle density ratio an order of magnitude estimate (see e.g. Rudinger, ⁴¹ Hjelmfelt and Mockros ⁴² ) indicates that the most significant force governing the dynamics of the particles is the drag force. In the current setup the drag coefficient is based on the Schiller–Naumann correlation ⁴³ and the Cunningham slip correction factor^44,45 is also incorporated in the formulation. Once the flow field data has statistically converged a total of approximately $75000$ particles are released into the domain from the inlet plane. The particles are tracked until they exit the computational domain. For each particle size, a mono-disperse setup is considered for the two different inflow conditions. One way coupling is assumed and particle-particle interactions are not considered. The initial velocity of the particles is the same as that of the flow velocity. Three different parameters associated with the particles are reported: residence time, radial position, and dispersion. The dispersion is calculated based on the formulation presented by Yan et al. ⁴⁶ :

$$\mathrm{Dis}{p}_x=\frac{1}{N}\sqrt{\sum _{j=1}^N{(p_{x,j}-p_{x0,j})}^2},$$

(7)

where $N$ is the total number of particles in the computational domain. The variable $p_{x,j}$ represents the position of the $j$ ’th particle in the streamwise direction and $p_{x0,j}$ represents the initial position of the $j$ ’th particle in the streamwise direction. The dispersion in other coordinate directions can be defined in a similar fashion.

Results

Solver validation

In order to validate the solver the backward facing step case presented in the experimental study by Fessler and Eaton ⁴⁷ is considered. It should be noted that this particular validation case does not consider particles of the same size as our present setup. As the Stokes number in our setup is small the particles will follow the flow streamlines. Thus, a benchmark case which exhibits similar flow characteristics, such as flow separation, to our geometry is selected for the purpose of validation. The step height of the geometry in the aforementioned experimental study ⁴⁷ is $h=0.0267$ . Furthermore, the Reynolds number based on the step height (h) and maximum center line velocity (uc) is 18400. The extent of the computational domain in the numerical study is $30h$ , $4h$ , and $2.49h$ in the streamwise, spanwise, and crossstream directions, respectively. The distance leading up to the step in the streamwise direction is $10h$ and after the step is $20h$ . On the solid walls no slip boundary conditions are imposed. In the spanwise direction periodic boundary conditions are applied. This is because in the experimental setup ⁴⁷ the aspect ratio of the geometry is large enough to assume the flow in the spanwise direction is homogeneous. Lastly, in the streamwise direction a turbulent velocity profile is imposed which has the same bulk velocity as the maximum center line velocity in the experimental setup ⁴⁷ . An unstructured polyhedral mesh is used for this numerical study. The average cell spacing, based on equation (2), for the mesh is $0.7·\times {10}^{-3}m$ and the total number of cells is $10.9\times {10}^6$ . Furthermore in the aforementioned mesh, the surface averaged non dimensional distance $y^{+}<1$ . The data is averaged for 18 flow through time, where a single flow through time is given by 30h/uc. In Figure 3, the flow profiles at various streamwise stations downstream of the geometry step are compared between the experimental study ⁴⁷ and LES based results. Qualitatively the time averaged flow profiles appear to be similar. For the root mean square quantity deviations are observed in regions of high shear with a maximum deviation of about $7\%$ . However, the flow profile comparison outside of these regions are qualitatively similar. Next, the particle tracking capabilities of the solver are validated by tracking copper particles of size $D_p=70\mu m$ . The corresponding Stokes number is $6$ . Furthermore, the mass loading ratio is $3\%$ . In the numerical setup only the Stokes drag force and pressure gradient force are incorporated. The particle velocity u_p is considered across various streamwise stations. Fifty equidistant bins in the wall normal direction are considered for each station in order to obtain particle statistics. In Figure 4 the comparison of the particle velocity profile obtained in the experimental ⁴⁷ and numerical results is presented. The data comparison shows slight differences between the two sets of results in the region associated with the developed shear-layer downstream of the step edge. Overall a reasonable agreement between the flow and particle profiles results obtained by numerical simulation and the ones presented by Fessler and Eaton. ⁴⁷

Figure 3.

Normalized flow profile comparisons (Experiments ⁴⁷ vs LES) at different stream-wise flow locations: (a) time averaged streamwise component of the velocity field and (b) root mean square of the streamwise component of the velocity field.

Figure 4.

Normalized averaged stream wise particle velocity comparisons between experimental ⁴⁷ and numerical (LES based) results at different downstream locations from the step.

Grid resolution study

In order to determine the optimal grid the absolute normal deviation ${\Pi }_{\mathrm{ij}}$ , as defined in Bodin et al., ⁴⁸ of selected flow quantities is considered. The flow case considered is the continuous inflow scenario. Time averaged streamwise component of the velocity field and the root mean square (RMS) of the streamwise velocity fluctuations are the selected flow variables used to assess the computational grid resolution and its impact on the fluid flow solution. In Table 2 the summary of the absolute normal deviation results are presented. The data are extracted across the orthogonal line corresponding to the maximum cross sectional diameter in the first module. A decrease in the absolute normal deviation is observed as the grid is refined. This trend is more obvious for the RMS of axial velocity fluctuations. Furthermore, the energy spectra calculated in monitoring points located along the symmetry axis of the geometry is presented in Figure 5. The lower plot in Figure 5 indicates the precise position of the point probes which are used to extract the velocity component based time signal data. These time series data are needed to calculate the spectra. The post processing is performed using the pwelch command in MATLAB, ⁴⁹ which computes the power spectral density using the Welch ⁵⁰ power spectral density estimate. From Figure 5 it is evident that as the grid resolution is increased a larger portion of the inertial sub-range is resolved. Based on the results observed for the absolute deviation and the energy spectra, grid $G3$ is selected as the most optimal mesh. All subsequent results that are presented are calculated using this grid.

Table 2.

Summary of the absolute normal deviation for the results obtained using different grids.

Line location	Flow quantity	${\Pi }_{12}$ (%)	${\Pi }_{23}$ (%)
$x/\mathrm{\lambda }=0.5$	$\overline{u_x}$	5.2	5
$x/\mathrm{\lambda }=0.5$	$\sqrt{\overline{{u\prime }_x{u\prime }_x}}$	16	5

Figure 5.

Comparison of the energy spectra taken for the three grids under consideration.

The flow profiles

Firstly, the flow field data of the continuous inlet mass flow case is considered. In Figure 6(a) and (b) the contour plots of the time averaged streamwise component of the velocity field and the turbulent kinetic energy are displayed. The planes being displayed are the mid longitudinal plane and selected circumferential planes. In particular, the precise position of the displayed circumferential planes are indicated by the outlines shown in Figure 1 and the planes are placed above or below the expanded or contracted section they correspond to, respectively. The velocity contours in Figure 6(a) highlight that the flow successively decelerates and accelerates moving from the expanded to contracted section. Furthermore, as the flow decelerates the boundary layer separates and thus giving rise to shear layers. Weak recirculation zones are also observed in the expansion region of the geometry. As highlighted in the TKE contour plots (see Figure 6(b)) the largest TKE content is found within the induced shear layers. We consider a single line profile for the flow quantities of interest since the flow field does not appear to vary significantly moving from one segment to the next. The precise position of the line segment is indicated by the dash dotted line in the contour plots within Figure 6(a) and (b). In Figure 6(c) the left plot presents the normalized time averaged streamwise component of the velocity field. An extremely small negative velocity is obtained in the line profile extracted at an orthogonal line located in the middle of the second module. The right plot in Figure 6(c) presents the normalized TKE profile, where the two peaks in the profile correspond to the shear layer formed downstream of the contracted section of the module due to flow separation.

Figure 6.

(a) Contour plots of time averaged streamwise component of the velocity field, (b) contour plots of turbulent kinetic energy, and (c) flow profiles of the time averaged streamwise component of the velocity field (left) and turbulent kinetic energy (right) taken across the maximum cross sectional plane in the second segment. The bulk velocity $U_b$ and the square of the bulk velocity are used to normalize the flow quantity of interest respectively.

In the pulsatile inlet mass flow scenario phase locked line profile data are considered at various phases. These phases are highlighted in Figure 2. The precise position of the line segment under consideration is the same as the one considered previously for the continuous inlet mass flow scenario. In Figure 7 the normalized phase locked averaged streamwise velocity and turbulent kinetic energy line profiles are presented in the top and bottom panel, respectively. In each plot the results for two different phases has been grouped together. Phase $1$ and $3$ correspond to the same mass flow rate but are placed in the decelerating and accelerating part of the pulse, respectively. The same holds true for phase $4$ and $6$ . Finally, phase $2$ and $5$ correspond to minimum and maximum mass flow rates associated with the pulsation. Due to the difference in the mass flow rate it is evident why the phase locked velocity in phase $5$ is larger than phase $2$ . In the other cases, it appears that the phase associated with the decelerating part of the pulse has a slightly larger phase locked velocity. In the case of $P1$ and $P3$ the pulse has to pass through $P2$ , which corresponds to the minimum mass flow rate of the pulse, before reaching $P3$ . Thus, the inertial effects of the pulse form passing through the minimum mass flow rate are evident in phase $3$ which is why the velocity at $P3$ is slightly smaller than $P1$ . A similar situation is observed between phase $4$ and $6$ . This time before reaching the decelerating phase ( $P6$ ) the pulse form passes through the phase corresponding to the maximum mass flow rate. Thus, in this instance inertial effects contribute toward the decelerating phase having a velocity almost equal to the corresponding accelerating phase $P4$ . In the line profiles associated with the turbulent kinetic energy it is evident that the phases which have the larger phase locked velocity have a larger turbulent kinetic energy. A larger flow velocity contributes toward a higher shear which in turn leads to a larger turbulent kinetic energy. For each turbulent kinetic energy profile two peak values can be observed across the normalized coordinates. These peaks are associated with the induced shear layers, in a similar fashion to the continuous mass flow case. In the phases associated with the upper half of the pulse ( $\stackrel{·}{m}/{\stackrel{·}{m}}_{\max }>0.88$ ) the asymmetry of the flow profile becomes more pronounced. Previous studies^51,52 have reported the development of flow asymmetry under axi-symmetric geometric configurations. The increase of velocity above the mean velocity of the pulsations may perturb the flow field in such a way that leads to more flow asymmetry.

Figure 7.

Line profiles of the phase locked streamwise component of the velocity field (top panel) and pulsatile turbulent kinetic energy (bottom panel). The mean pulsation velocity $U_0$ and the square of the mean pulsation velocity are used to normalize the flow quantity of interest respectively. The orthogonal line across which the data are extracted from is the same as the one considered in the continuous flow scenario.

Particle behavior

The first thing considered in the particle behavior within the sinusoidal geometry is the particle dispersion in the $x$ , $y$ , and $z$ -direction. In Figure 8(a) and (b) the particle dispersion observed for continuous and pulsating inflow conditions are displayed, respectively. The $x$ -axis in each plot denotes the non-dimensional time. In the continuous inflow scenario, the flow through time (3λ/Ub) is used to normalize time. However, in the pulsating flow scenario time is normalized using the flow pulsation frequency f. It is evident that the particle with the largest size has the largest amount of dispersion. Furthermore, the streamwise dispersion under the pulsatile inflow scenario for the largest size particle is larger than the streamwise dispersion obtained under continuous inflow conditions for the very same particle size. Despite all particles having Stokes number significantly smaller than unity, it appears that inertia along with flow pulsations play a role in increasing flow dispersion, particularly in the streamwise direction. In the next step the residence time of an individual particle within each module is quantified. Figures 9 and 10 display histogram plots associated with the residence time of particles within an individual module under continuous and pulsatile inflow conditions. The abscissa of the histogram plots denotes the normalized residence time of the particles, where the time scale based on the geometric wavelength and bulk inlet velocity is used to normalize the residence time. Clearly, the majority of particles correspond to the three smallest bin associated with the residence time in the histogram plots. Furthermore, this is observed for all particle sizes under both inflow conditions. This means that particle retention in this geometry is very weak. While the fraction of particles associated with the second bin are approximately the same between the second and third module, there is however a difference in the histogram plots associated with the first module. In particular, the residence time of the particles appears to be slightly more spread out in terms of fraction of particles when compared to the last two modules. This is because after injection some particles get trapped in the boundary layer within the first module. There is also a slight difference in the residence time histogram plots associated with the first module under continuous and pulsatile inflow conditions. The fraction of particles associated with the second bin in the first module is larger in the pulsatile case when compared with the continuous flow case. The imposed pulse form (see Figure 2) is such that it decelerates first and then accelerates. As the convective velocity is reduced initially, a larger fraction of particles is retained, for a slightly longer time, in the first module as compared to the continuous flow scenario.

Figure 8.

Particle dispersion in the $x$ , $y$ , and $z$ -direction under: (a) continuous and (b) pulsatile inflow conditions.

Figure 9.

Histogram plots of particle residence time for each module under continuous inflow conditions for particle size: (a) $2.5\mathrm{\mu }\mathrm{m}$ , (b) $800\mathrm{nm}$ , (c) $500\mathrm{nm}$ , and (d) $180\mathrm{nm}$ .

Figure 10.

Histogram plots of particle residence time for each module under pulsatile inflow conditions for particle size: (a) $2.5\mathrm{\mu }\mathrm{m}$ , (b) $800\mathrm{nm}$ , (c) $500\mathrm{nm}$ , and (d) $180\mathrm{nm}$ .

By considering the particle radial positions across each individual module the tendency of the particles to be trapped in the boundary layer becomes more evident. In Figures 11 and 12 the histogram plots associated with the average radial position of a particle within an individual module are presented for continuous and pulsatile inflow conditions, respectively. The $x$ -axis of the histogram plots represents the normalized average radial position of the particle, where the maximum cross sectional radius is used to normalize the average particle radial positions. The average radial position of an individual particle in the last two modules is far less skewed as compared to the average radial position of particles in the first module. Most of the particles in the first module tend to occupy the radial position slightly larger than the contracted section. Some of these particles are amongst the earlier mentioned particles which get caught in the boundary layer.

Figure 11.

Histogram plots of the average radial position of a particle within each module under continuous inflow conditions for particle size: (a) $2.5\mathrm{\mu }\mathrm{m}$ , (b) $800\mathrm{nm}$ , (c) $500\mathrm{nm}$ , and (d) $180\mathrm{nm}$ .

Figure 12.

Histogram plots of the average radial position of a particle within each module under pulsatile inflow conditions for particle size: (a) $2.5\mathrm{\mu }\mathrm{m}$ , (b) $800\mathrm{nm}$ , (c) $500\mathrm{nm}$ , and (d) $180\mathrm{nm}$ .

Discussion and conclusions

The present numerical study considers mono-disperse simulations of low Stokes number particles under continuous and pulsating inflow conditions inside a sinusoidal geometry. Periodic boundary conditions are imposed in the streamwise direction. The flow field under different inflow conditions is firstly characterized. Particles are then introduced inside the computational domain. A Lagrangian point force approach in a one way coupling manner is utilized. Only the Stokes drag force is incorporated in the numerical setup. In order to get an insight into the particle behavior within the sinusoidal pipe geometry the particle dispersion, residence time within each module and the average radial position of a particle within a module are considered. In all dispersion plots an initial increase in the dispersion value is observed up to approximately one non-dimensional unit of time. At later times fluctuations are observed in the dispersion plots. Due to the weak recirculation regions the vast majority of the particles leave the computational domain within the first five non dimensional units of time. The fluctuations observed at later non-dimensional units of time in the plots are based on the remaining particles. Histogram plots associated with the particle residence time highlight that the vast majority of the particles are not retained in the domain for a long duration of time. This can be linked to the weak recirculation zones within the expansion region of the cavities. As most of the particles do not reside within any individual module for a long duration of time, the likelihood of particles interacting with other particles is small. However, in the pulsatile flow scenario the flow decelerates first and as a result of the reduced convective velocity particles stayed in the first module for a slightly longer duration as compared to the continuous flow case. Since the Stokes number associated with all the particles considered in this study is less than unity, once the particles move out of the first module they follow the flow and move out of the computational domain. This is further highlighted in the histogram plots associated with the average radial position of the particles, where the average radial position of the particles in the last two modules is nearly identical. There are some particles that get trapped in the boundary layer of the first module upon their release from the inlet plane of the computational domain. This behavior is exhibited in the slightly skewed distribution of radial position of particles in the histogram plots corresponding to the first module of the geometry.

There are two main limitations associated with the present analysis. The grid resolution utilized for this study is only appropriate for a wall modeled LES turbulence modeling approach. Secondly, the current setup has a mean free path of approximately $65\mathrm{nm}$ . This means that some of the particles, particularly the smallest sized particles, fall out of even the slip regime. In order to further quantify the effects of particle inertia it may be useful in future studies to consider significantly larger sized particles and also poly-disperse simulations. Furthermore, in future studies additional flow pulsation frequencies should also be considered such that the time scale associated with the pulsations is smaller than the flow through time (convection time scale).