Lift on Janus and stick spheres in laminar channel flow: a computational study

Abstract

The instantaneous motion of a spherical particle in a channel flow is governed by the forces experienced by the particle. The magnitude and direction of the forces depend on the particle to channel size ratio, particle position, nature of the sphere surface (sticky/slippery), fluid properties and relative velocity between the fluid and the particle. In this work, we report the lift, the force component directed normal to the streamwise direction, on two classes of spheres, sticky and Janus, in a channel of square cross-section. The Janus spheres considered have both sticky and slippery hemispheres with the boundary between the two hemispheres parallel to the channel midplane. The effect of particle to channel size ratio, dimensionless particle position and particle Reynolds number on the lift are studied. The Janus sphere placed at the channel centerline is observed to experience the lift directed from the sticky to the slippery hemisphere. A correlation is proposed to predict the lift on the Janus sphere placed at the centerline of the channel. A sticky sphere positioned close to the channel wall experiences a significant lift directed away from it. For the Janus sphere placed at an off-center position two possibilities arise—slippery hemisphere facing the channel centerline (case A) or sticky hemisphere facing the channel centerline (case B). For case A, the lift is always directed away from the wall. For case B, the direction of lift depends on the particle position as well as particle Reynolds number. The moment coefficients for the sticky and Janus sphere are also presented.

Appendices

Appendices

Appendix 1. Grid Independence

The three-dimensional computational domain is discretized with a structured, hexahedral mesh. The mesh is refined near the sphere, as can be observed from Fig. 19a. The first cell height near the sphere is a/100. Three different meshes having the number of elements ranging from 0.05, 0.20 and 0.40 million for $a/H = 0.20$ and $Y/H = 0$ are used for simulations. The thickness of the boundary layer on the wall decreases with an increase in the particle Reynolds number; therefore, the highest particle Reynolds number ($\text {Re}_{\mathrm {P}} = 80$) used in this work is chosen for the grid independence study. It can be observed that there is a good agreement between the results obtained using grids having 0.20 and 0.40 million elements, and therefore, further investigations are done using 0.2 million mesh size. Figure 19(b) and (c) show the plot of x-component of velocity on a line at the center of the midplane at a distance of 21H from the inlet and at the center of the sphere. The distance, 21H, is close to the rear sphere and that is why there is a dip in the x-velocity just after the sphere due to the vortex formation in that region. Here, $V_{\mathrm {x}}$ is the velocity in the x direction and V is the channel inlet velocity.

Fig. 19

(a) Mesh distribution on a cross section of the channel passing through the middle of the sphere. Distribution of x component of velocity on a line passing through the middle of the cross-sectional plane, (b) at a distance of 21H from the inlet (c) on the line passing through sphere for three different mesh sizes

Appendix 2. Validation

1.1 Velocity profile in the channel

In the simulations, the flow at the inlet boundary is defined to have a uniform velocity over the entire inlet area. The velocity profile becomes fully-developed downstream in the channel for the particle Reynolds numbers investigated. The analytical expression for the velocity in the laminar, fully-developed flow in a rectangular channel [45] is given by Eq. A1.

$$\begin{aligned} \frac{\mathrm {v}_{\mathrm {x}}}{\mathrm {V}}\mathrm {=} \frac{\mathrm {48}}{{\uppi }^{\mathrm {3}}}\mathrm {}\left\{ \frac{\sum \nolimits _{\mathrm {n}=1,3,\ldots }^{\infty } {\mathrm {}\frac{\mathrm {1}}{\mathrm {n}^{\mathrm {3}}}\left( \mathrm {-1} \right) ^{\frac{\mathrm {n-1}}{\mathrm {2}}}\left[ \mathrm {1-}\frac{\cosh \left( \frac{\mathrm {n}\uppi \mathrm {y}}{\mathrm {W}} \right) }{\cosh \left( \frac{\mathrm {n}\uppi \mathrm {H}}{\mathrm {2W}} \right) } \right] \mathrm {cos}\left( \frac{\mathrm {n}\uppi \mathrm {z}}{\mathrm {W}} \right) } }{\left[ \mathrm {1-}\frac{\mathrm {192}}{{\uppi }^{\mathrm {5}}}\left( \frac{\mathrm {W}}{\mathrm {H}} \right) \, \sum \nolimits _{\mathrm {n}=1,3,\ldots }^{\infty } {\mathrm {}\frac{\mathrm {1}}{\mathrm {n}^{\mathrm {5}}}\mathrm {tanh}\left( \frac{\mathrm {n}\uppi \mathrm {H}}{\mathrm {2W}} \right) } \right] } \right\} \end{aligned}$$

(A1)

Here, V is the mean velocity, W is the channel width, and H is the channel height. The velocity profile for four different $\text {Re}_{\mathrm {P}}$ is compared with the analytical expression. It can be seen in Fig. 20 that the profile is in good agreement with the analytical expression.

Fig. 20

Comparison of velocity profile obtained from CFD simulation with analytical solution (Eq. A1) on a line passing through the middle of the cross-sectional upstream of the sphere for different particle Reynolds numbers for $a/H = 0.20$

1.2 Pressure distribution on the sphere

Analytical solution for flow around a sphere can be obtained for Stokes flow around a sphere [2, 46]. To establish veracity of the simulations, the pressure distribution for creeping flow around the sphere is compared with that obtained from analytical solution given by Eq. A2.

$$\begin{aligned} P-P_{{\infty }}= \frac{\mathrm {3}}{\mathrm {4}} \frac{\mu aVCos \uptheta }{r^{2}} \end{aligned}$$

(A2)

where $P_{{\infty }}$ is the pressure far away from the sphere, and r is the distance from the center of the sphere. The pressure, denoted in terms of non-dimensional pressure coefficient, as defined in Eqs (A3, A4), is plotted in Fig. 21.

$$\begin{aligned} C_{\mathrm {p}}= & {} \frac{P\mathrm {-}P_{{\infty }}}{\frac{\mathrm {1}}{\mathrm {2}}\rho V^{2}} \end{aligned}$$

(A3)

$$\begin{aligned} C_{\mathrm {p}_{\mathrm {0}}}= & {} \frac{P_{\mathrm {0}}\mathrm {-}P_{{\infty }}}{\frac{\mathrm {1}}{\mathrm {2}}{\uprho }\mathrm {V}^{\mathrm {2}}} \end{aligned}$$

(A4)

The values obtained from the CFD simulations are in excellent agreement with those obtained from analytical solution.

Fig. 21

The pressure distribution on the sphere in creeping flow conditions ($\text {Re}_{\mathrm {P}}= 0.01$). The values are in excellent agreement with the analytical expression. The inset figure shows the angular direction of the flow, where $\theta = 0$ means front stagnation point of the sphere and $\theta = \pm \pi$ means rear stagnation point

1.3 Uniform flow over a Janus Sphere

We have also investigated the uniform flow over the Janus particle where the front hemisphere of the Janus sphere has zero shear (slip) condition, and no slip (stick) on the back hemisphere. The drag coefficient obtained through simulations is compared with the analytical solution, given by Eq. (A5) proposed by Swan and Khair [47] for stokes conditions where $\lambda$, slip length $= 2a/3$. It can be seen from Fig. S4 that the drag coefficient values obtained from CFD, by and large, are in agreement with the expression [3] for Janus sphere Fig. 22.

$$\begin{aligned} F_{\mathrm {D}}=3\pi \upmu aV\left( 1-\frac{\lambda }{2a} \right) \end{aligned}$$

(A5)

Fig. 22

The comparison of obtained drag coefficient with the analytical expression for the case of uniform flow around the Janus sphere

Fig. 23

Lift coefficient variation with particle Reynolds number when the Janus sphere is placed in a uniform flow

We have also obtained the lift coefficients in a uniform flow over the Janus sphere, as shown in Fig. 23. It can be seen that with the increase in particle Reynolds number, lift coefficient decreases. Unsurprisingly, the magnitude of lift coefficient remains same when free slip and no slip are applied interchangeably. The qualitative nature of the plot is linear suggesting a proportionality in the obtained values with particle Reynolds number.

Appendix 3. Correlation description

The values of the constants; A, B and C mentioned in Eq. (A6), are determined for all the particle Reynolds numbers.

$$\begin{aligned} C_{L}=A+B\,exp\left( C\frac{a}{H} \right) \end{aligned}$$

(A6)

The constant A is then fitted in a form given by Eq. (A7):

$$\begin{aligned} A=c_{1}{Re}_{P}^{2}+c_{2}{Re}_{P}+c_{3} \end{aligned}$$

(A7)

The $\hbox {R}^{\mathrm {2}}$ value of 0.993 suggests that the polynomial fit works fine for the term A. Similarly, the values obtained for B are plotted against particle Reynolds number which can be observed from Fig.  24a, B term. The $\hbox {R}^{\mathrm {2}}$ value of 0.994 suggests that the power fit is appropriate.

      The equation for constant, B, is of the form:

$$\begin{aligned} B=d_{1}exp\left( d_{2}{Re}_{P} \right) \end{aligned}$$

(A8)

In the same manner the values for constant C are varied against local particle Reynolds number as can be seen from Fig. 24b. The $\hbox {R}^{\mathrm {2}}$ value of 0.998 also reveals that the polynomial fit is good.

      The equation for constant, C, is of the form:

$$\begin{aligned} C=f_{1}+f_{2}{Re}_{P} \end{aligned}$$

(A9)

Thus the correlation is established between A, B and C and is then substituted in Eq. 5. The variation of $C_{L}$, with $\text {Re}_{\mathrm {P}}$ appears to vary in the form as mentioned below; to obtain the correlation in Eq. (5):

$$\begin{aligned} C{}_{L}\mathbf {=}\left( -0.0002{Re}_{P}^{2}+0.026{Re}_{P}+0.385 \right) +\left( 0.136\exp \left( -0.143{Re}_{P} \right) \mathrm {exp}\left( \left( 0.156{Re}_{P}+4.52 \right) \frac{a}{H}\right) \right) \nonumber \\ \end{aligned}$$

(A10)

Fig. 24

(a) The constant, A, gets a trend line through polynomial fit against particle Reynolds number giving a high $\hbox {R}^{\mathrm {2}}$ value, (b) The constant, B, gets a trend line through exponential fit curve against particle Reynolds number giving a high $\hbox {R}^{\mathrm {2}}$ value, (c) The constant C gets a linear trend line for the range of particle Reynolds number investigated