Spherical Particle Migration Evaluation in Low Reynolds Number Couette Flow Using Smooth Profile Method

Solid–liquid multiphase flow analysis is an interdisciplinary research area with various technological applications. From the sediment transport in rivers, fillers motion within polymers, to high-performance coolants, having a correct understanding of particle–fluid flow interactions plays a key role in research and development. Fluid flows containing small-size particles occur in biological and engineering systems, including aerosol transport, air pollution, paper-making, and targeted drug delivery processes. Analyzing fluid flow with suspended particles is important not only for predicting the motion of particles within the primary fluid phase but also for evaluating their impact on the flow of the primary phase. Computational methods have been developed to analyze and evaluate the multiphase flow in different scales. Multiphase flow modeling using computational fluid dynamics goes back to the 1950s and the 1960s when Particle-in-Cell (PIC) and Marker-and-Cell (MAC) methods were developed at Los Alamos National Laboratories [1–4]. A recent review of the MAC method can be found in Mckee et al. [5]. Although the MAC method provided significant progress in numerical modeling of the multiphase flows, it was relatively inaccurate. The volume of fluid (VOF), front tracking, as well as level set methods, were more robust and accurate techniques proposed in the late 1980s to investigate multiphase flows numerically [6–8]. The VOF approach uses a color/marker function representing the phase fraction of each fluid within a computational cell. The main problem with the VOF method is capturing the sharp boundary between different phases and interfacial properties such as surface curvature and surface tension [9]. Front Tracking Method introduced in the late eighties uses the Eulerian approach to solve a single set of governing equations, that is, Navier–Stokes and continuity equations for incompressible fluid flow, in the entire computational domain. Sharp changes in material properties such as viscosity and density can be captured, while the effect of surface tension is included by adding an appropriate source term at the interface between different phases in the front tracking technique [10]. The interface is directly captured by tracking the marker points on the boundaries between different phases. The direct representation of the boundaries between different fluids makes this method suitable and accurate to investigate and capture interfacial phenomena such as surface tension, force balance, phase change, and surface curvature with a high order of accuracy. Recently, other analytical and computational investigations were performed on particle–fluid interaction. Ingber et al. [11] evaluated different particle inter-actions, including particle/particle, wall/particle, and particle migration, in different types of nonlinear shear fields. A semi-analytical solution was developed for the motion of two spherical particles suspended in an unbounded arbitrary shear flow by the same group [12]. The results of the discrete phase element method (DPM) against the analytical solution were evaluated, and the verification and validation (V&V) for the single and double particle tra-jectories in a rectangular and cylindrical domain were conducted in [13–16]. Previously, our progress on applying the DPM using the ANSYS/Fluent code to model the particle migration was reported in [17]. The current investigation is dedicated to studying the particle behavior in a low Reynolds number Couette shear fluid flow under different boundary conditions using our recently developed Eulerian–Lagrangian computational code using a similar scheme as the smooth profile (SP) method [18, 19] to simulate the coupled solid and fluid phases simultaneously.

A finite volume CFD code has been developed to study the migration of a solid particle in low Reynolds number Couette flows between two parallel plates. Figure 1 shows the compu-tational domain used in the numerical simulations. The left and right boundaries are assigned as a periodic boundary condition.

The incompressible Navier-Stokes (N-S) equation is given as follows:

where $u_f$ is the fluid velocity, $P$ is the relative pressure with the hydrostatic pressure as the reference, and $\rho_f$ and $\mu_f$ are the density and viscosity of the fluid phase, respectively. Since the fluid is assumed to be incompressible, the continuity equation simplifies to the following:

The translational and rotational motions of the solid particle within the flow field are considered using the Lagrangian approach. These are as follows:

where $m_s$ is the particle mass, $u_s$ is the particle velocity, $F_f$ is the hydrodynamic force imposed by the fluid on the particle, and $F_b$ is the buoyancy force acting on the spherical particle. Also, $I_s$ is the rotational inertia of the particle, $\omega$ is the particle angular velocity, and $T_f$ is the torque inserted by the fluid to the particle at the interface. The hydrodynamic force and torque exerted by the fluid flow on the solid particle are given in Eq. (4):

To overcome discontinuity at the interface in this Eulerian-Lagrangian approach, the interface between solid particle and primary fluid phase is replaced by a kernel function, $\psi$, creating a SP from solid into the liquid with a predefined thickness [1,2]. The only input parameter for this smoothed kernel function is the thickness of the interface, $\delta$. The corresponding kernel function for the smoothed profile that is used in the present work is given by [2] as follows:

where $r$ is the radius of the particle, $x_c$ is the particle's center of mass, $d x$ is the spatial discretization size, and $\delta$ is the interface thickness, which is the only input parameter to define the smoothing kernel function. As it can be inferred from Eq. (5), the kernel function $\psi$ is 1 within the solid particle region, $\left|x-x_{\mathrm{c}}\right|<(\mathrm{r}-\delta / 2)$, zero in the fluid region, $\left|x-x_{\mathrm{c}}\right|>(\mathrm{r}+\delta / 2)$, while varying smoothly from 1 to 0 within the interface thickness. Using the kernel function, the velocity field in the entire solution domain can be expressed as follows:

where $u_s$ designates the velocity of the solid particle. With a well-defined time step and initial velocities within solid and fluid domains, Eqs. (1-6) are solved numerically using the smoothed profile method proposed by Nakayama and Yamamoto [16, 17]. The finite volume technique is applied to the N-S equations along with the continuity equation, given by (1) and (2). The derived system of equations is solved simultaneously to obtain fluid velocity distribution as well as pressure field with the use of an explicit second-order projection method [19,20] using spatially central-difference approximation on a fixed, staggered grid configuration. Obtained results are considered as intermediate velocity and pressure fields $\left(u f^*, P^*\right)$. Total velocity is equal to the velocity of the particle within the particle domain, $\left|x-x_{\mathrm{c}}\right|<(\mathrm{r}-\delta / 2)$.

At the interface, due to the no-slip boundary condition, the velocity of both phases is equal. These two constraints are applied by correcting the intermediate pressure and velocity fields. The corrected pressure field, $P_c$ and fluid velocity are used to study the disturbing effect of the solid particle on the fluid velocity field.

The total velocity field defined in Eq. (6) should be divergence-free. By taking the divergence of Eq. (6) and considering Eq. (7), the Poisson equation for the correction pressure field is derived as follows:

Consequently, Eq. (4) converts to the following equation:

After obtaining the hydrodynamic force and torque, particle translational and angular velocities are updated using Eq. (3). This procedure repeats for the next time step.

The initial motivation was to develop a new particle model originated from our experimental analyses, which have been conducted for several years. Figure 2 shows the experimental setup. As can be seen, the experimental setup consists of a tank and a computer-controlled stepper motor (compumotor). The compumotor is used to displace the belts shown in Figure 2. In the setup to determine the behavior of a spherical particle suspended in linear shear flow, one spherical polymethyl-metacrylate particle with a diameter of 6.35 mm is suspended between two layers of fluids with approximately equal viscosity and different density. The particle was placed into a tank between two belts, as shown in Figure 3.

To assess the accuracy of the developed CFD code, in the first scenario, a single-phase Couette flow between two parallel plates with a distance of 0.1 m is considered while the top plate moves with a velocity of 1 m/s toward the right and the bottom plate moves with the same velocity in the opposite direction (toward left). The density of the fluid has been considered as 1 kg/m³, and the fluid viscosity is 0.001 Pa. S. The simulations have been conducted for mesh grid sizes of 3 × 3, 6 × 6, 12 × 12, 24 × 24, and 48 × 48 cells, respectively. Figure 4 and Figure 5 show the numerical results as well as exact solutions for u velocity in the middle of the domain for a mesh grid of 6 × 6 cells and 48 × 48 cells.

As it can be observed from Figure 4 and Figure 5, the developed CFD solver gives accurate results even with the coarse mesh of 6 × 6 mesh elements. This is mainly because of the simplicity of the problem since there are no advection and pressure terms, and the flow field is mainly controlled by viscous diffusion.

The L₁, L₂, and L_∞ norms for single-phase Couette flow are presented in Table 1 in order to investigate the observed order of accuracy of the implemented CFD solver. Figure 6 also shows the obtained observed order of accuracy.

As shown in Table 1 and Figure 6, the observed order of accuracy is 4.5 for a mesh grid of 6 × 6 elements, and then decreases to 1.5 for a mesh grid of 12 × 12 elements, while it reaches 0.1 for a mesh grid of 48 × 48 elements. The formal order of accuracy is 2, both for spatial and temporal variables. This observed order of accuracy seems strange, while it reduces continuously by mesh refinement.

The exact solution of this flow is a straight line, and so the discretization error would be zero if there were no round-off errors and iterative errors. Round-off errors are unavoidable, and the results of Table 1 suggest that single precision was used. It is observed that the summation of truncation and round-off error for this simple flow becomes larger than the discretization error after the mesh grid size of 6 × 6 elements. As a result, the observed order of accuracy reduces continuously with further mesh refinement. Although the observed order of accuracy is less than 2 for Couette flow, Figure 5 shows that the obtained result for velocity matches very well the exact (analytical) solution for a mesh grid size of 48 × 48 elements.

The second verification scenario for a pressure-driven flow between two parallel plates is considered to further verify the code accuracy. The channel is 0.2 m long and 0.01 m wide with a periodic boundary condition at left and right ends, and an applied pressure gradient of 240 Pa. Fluid density is 1 kg/m³, and its viscosity is 0.001 Pa.S. Numerical simulations have been conducted for grid sizes of 8 × 8, 16 × 16, 32 × 32, 64 × 64, and 128 × 128 mesh elements.

Figure 7 compares velocity profiles predicted by the present CFD code for the primary phase with mesh grids of 16 × 16 and 128 × 128 elements with the exact solution. The CFD results match the exact solution fairly well, confirming the accuracy of the single-phase flow solver. Since the primary phase is a fully developed flow (solution does not change in the x-direction), the error was calculated along the y-direction. Note that the errors were calculated for cross-sections of the domain and they differ only at the level of round-off error.

The observed order of accuracy from the CFD results of a pressure-driven flow between two parallel plates for different norms is presented in Figure 8. As can be seen, the observed order of accuracy for the developed Fortran CFD code approaches the formal order of accuracy, which is 2. In contrast with the single-phase Couette flow, pressure and viscous terms affect the flow field in this simulation scenario. The truncation and round of errors are much smaller than discretization error, and therefore, the observed order of accuracy is 2, which is consistent with the formal order of accuracy of the applied numerical algorithm in the developed code.

The migration of a solid particle in a Couette flow is studied. The above-mentioned periodic and the no-slip wall boundary conditions in the x- and y-direction, respectively, are used (Figure 3). The following nondimensional numbers are also used in the study:

where $D$ is the particle diameter, $H$ is the distance between two parallel walls in the Couette flow, and $V$ is the velocity magnitude of the upper and lower walls. Table 2 summarizes the corresponding nondimensional numbers used in this study.

The experimental data on the migration of particles in a low Reynolds number Couette flow from our experimental setup was used to perform validation. The top belt moves with a velocity of $0.002 \mathrm{~m} / \mathrm{s}$ to the left, while the bottom belt moves with the same velocity to the right. The particle is initially located near the top belt, while its center of mass is 1.25 $D$ below the belt. To study the effect of the slip boundary condition, a different set of simulations, including the slip boundary condition instead of no-slip boundary conditions were performed. Slip boundary condition can be modeled numerically using the concept of slip length, $\beta$ [16]. Slip length measures a distance from the interface where the fluid velocity reaches to zero. Eq. (11) is implemented in the performed numerical simulation in order to find the fluid slip velocity at the walls of the computational domain (see Figure 1).

Figure 9 compares the experimentally measured and numerical predicted migration rates. It is observed that the particle migrates slowly toward the channel centerline. The migration rate is quite small, and the uncertainty in the experimental data is considerable. However, the CFD prediction and the experimental data follow the same trend. The ratio between the particle’s center of mass relative velocity and belt velocity as the applied CFD boundary condition is presented in Figure 10. It is seen that there is a considerable difference between the

experimental data and the CFD prediction for the case with no-slip boundary condition. The CFD result with a no-slip boundary condition at the belts suggests that the relative velocity ratio is around 0.21, while experimental data gives a value of 0.7.

In the performed numerical simulations with the no-slip boundary condition for the walls, the velocity profile within the fluid domain is high near the top and bottom belts such that the particle centerline velocity is around 79% of belt velocity, as shown in Figure 10. However, the velocity ratio obtained from the conducted experiment in Figure 10 is around 0.3. The deviation between the obtained numerical results and experimental data for the particle centerline velocity illustrated in Figure 10 is due to fluid slippage at the belts in the experiments.

It is observed in Figure 9 and Figure 10 that by considering a slip length of 4.5 μm, the velocity ratio and particle migration rate approach to the experimental data. Therefore, the current numerical simulations predict a slip length of 4.5 μm for the fluid near the belts in the performed experiments. Further numerical simulation with an apparent slip length of $\beta$ = 4.5 mm was conducted, and the result is presented in Figure 9 and Figure 10. It is observed that by considering apparent slip, the velocity ratio and particle migration rate get closer to the exper-imental data. Further numerical simulations at different Reynolds numbers were performed to investigate the effect of shear rate on particle migration in a Couette flow. Figure 11 presents the trajectory of the center of the mass of the particle as a function of time for different Reynolds numbers of 0.04–0.32. As it is seen, the particle migrates toward the center of the

flow domain $\left(y_c / H=0.5\right)$, where the shear stress is the lowest. The initial position of the solid particle is the same for all Reynolds numbers, and simulations are terminated once the particle leaves the solution domain from left or right boundaries (see Figure 2). At high Reynolds numbers, the particle travels faster, and therefore, as presented in Figure 11, the corresponding simulation time is shorter.

The inertial lift force on a migrating particle is exerted within a shear flow. This force is first observed and formulated by Saffman in the 1960s [21]. In this work, the combination of imposed inertia and lift forces from the fluid flow resembles the Saffman lift force and pulls the particle toward the channel centerline, as can be inferred from migration rates in Figure 12. The streamlines around the particle at Reynolds number of 0.32 are shown in Figure 13. A counterclockwise rotation of the particle can be observed, which tends to push the particle upward toward the top plate.

Figure 14 illustrates the angular velocity of the particle as a function of time for the Reynolds number in the range of 0.04–0.32. Angular velocity increases as the Reynolds number gets higher, and so does the fluid velocity near the particle surface. However, the observed trend for particle migration rate in Figure 11 suggests that the inertia force from fluid flow dom-inates in the case of low Reynolds number Couette flow (Re < 0.32) and, therefore, particle migrates toward the centreline. The average angular velocity correlates with the Reynolds number and can be fitted to Eq. (12):

As it can be inferred from Figure 14 and Eq. (12), in the low Reynolds number Couette flow ($\operatorname{Re}<0.32$), the induced angular velocity of the particle varies almost linearly with the Reynolds number.

The results of the newly developed Eulerian–Lagrangian model for a spherical single-particle migration in a low Reynolds number shear flow were presented. Code verification and solution verification were conducted for the primary phase. The solid particle (secondary phase) migration was evaluated by performing validation against the experimental data obtained from our low Reynolds apparatus. The obtained computational results of the solid particle migration rate and relative transport velocity showed good agreement with the experimental data confirming the accuracy and reliability of the Eulerian–Lagrangian approach. It was observed that the solid particle tends to migrate toward the channel centerline at all Reynolds numbers (0.04–0.32). The obtained induced average angular velocity of the solid particle was observed to vary almost linearly with the Reynolds number.