2.1 Multiple-relaxation-time (MRT) lattice Boltzmann method

The MRT-LBM (Lallemand & Luo Reference Lallemand and Luo2003) is used to solve the fluid flow governed by the incompressible Navier–Stokes equations. The lattice Boltzmann equations (LBE) (d’Humiéres et al. Reference d’Humiéres, Ginzburg, Krafczyk, Lallemand and Luo2002) can be written as

(2.1) $$\begin{eqnarray}|\!\,f(\boldsymbol{x}+\boldsymbol{e}_{i}\unicode[STIX]{x1D6FF}t,t+\unicode[STIX]{x1D6FF}t)\!\rangle -|\!\,f(\boldsymbol{x},t)\!\rangle =-\unicode[STIX]{x1D648}^{-1}\hat{\unicode[STIX]{x1D64E}}[|\!\,m(\boldsymbol{x},t)\!\rangle -|\!\,m^{eq}(\boldsymbol{x},t)\!\rangle ],\end{eqnarray}$$

where the Dirac notation of ket $|\!\,\cdot \!\rangle$ vectors symbolize the column vectors. $|\!\,f(\boldsymbol{x},t)\!\rangle$ represents the particle distribution function, which has 19 components $f_{i}$ with $i=0,1,\ldots ,18$ because of the D3Q19 model used in our 3D simulations. The collision matrix $\hat{\unicode[STIX]{x1D64E}}=\unicode[STIX]{x1D648}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\unicode[STIX]{x1D648}^{-1}$ is diagonal with

(2.2) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D64E}}\equiv (0,s_{1},s_{2},0,s_{4},0,s_{4},0,s_{4},s_{9},s_{10},s_{9},s_{10},s_{13},s_{13},s_{13},s_{16},s_{16},s_{16}),\end{eqnarray}$$

where the parameters of $\hat{\unicode[STIX]{x1D64E}}$ are chosen as (d’Humiéres et al. Reference d’Humiéres, Ginzburg, Krafczyk, Lallemand and Luo2002): $s_{1}=1.19$ , $s_{2}=s_{10}=1.4$ , $s_{4}=1.2$ , $s_{9}=1/\unicode[STIX]{x1D70F}$ , $s_{13}=s_{9}$ , $s_{16}=1.98$ . $|\!\,m^{eq}\!\rangle$ is the equilibrium value of the moment $|\!\,m\!\rangle$ , where the moment $|\!\,m\!\rangle =\unicode[STIX]{x1D648}\boldsymbol{\cdot }|\!\,f\!\rangle$ , i.e.  $|\!\,f\!\rangle =\unicode[STIX]{x1D648}^{-1}\boldsymbol{\cdot } |\!\,m\!\rangle$ . $\unicode[STIX]{x1D648}$ is a $19\times 19$ linear transformation matrix which is used to map the column vectors $|\!\,f\!\rangle$ in discrete velocity space to the column vectors $|\!\,m\!\rangle$ in moment space. The matrix $\unicode[STIX]{x1D648}$ and $|\!\,m^{eq}\!\rangle$ are the same as those used by d’Humiéres et al. (Reference d’Humiéres, Ginzburg, Krafczyk, Lallemand and Luo2002) and Huang et al. (Reference Huang, Yang, Krafczyk and Lu2012b ). In (2.1), $\boldsymbol{e}_{i}$ are the discrete velocities. For the D3Q19 velocity model,

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{i}=c\left[\begin{array}{@{}ccccccccccccccccccc@{}}\,0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0\\ \,0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\ \,0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & -1 & 1 & -1\end{array}\!\!\right]\!, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $c$ is the lattice speed, defined as $c=\unicode[STIX]{x0394}x/\unicode[STIX]{x0394}t$ . In our study $\unicode[STIX]{x0394}x=1lu$ and $\unicode[STIX]{x0394}t=1ts$ , where $lu$ and $ts$ represent the lattice unit and time step, respectively. $mu$ is used to denote the mass unit. The macro-variables of fluid flow can be obtained from

(2.4a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D70C}=\mathop{\sum }_{i}f_{i},\quad \unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D701}}=\mathop{\sum }_{i}f_{i}e_{i\unicode[STIX]{x1D701}},\quad p=c_{s}^{2}\unicode[STIX]{x1D70C},\end{eqnarray}$$

where subscript $\unicode[STIX]{x1D701}$ denotes three coordinates. The parameter $\unicode[STIX]{x1D70F}$ is related to the kinematic viscosity of the fluid: $\unicode[STIX]{x1D708}=c_{s}^{2}(\unicode[STIX]{x1D70F}-0.5)\unicode[STIX]{x0394}t$ , where $c_{s}=c/\sqrt{3}$ is the sound speed.

The numerical method used in our study is based on the MRT-LBM (d’Humiéres et al. Reference d’Humiéres, Ginzburg, Krafczyk, Lallemand and Luo2002) and the dynamic multi-block strategy (Huang et al. Reference Huang, Yang and Lu2014).

Figure 1. Schematic diagram of the combination of coordinate transformation from $(x,y,z)$ to ( $x^{\prime },y^{\prime },z^{\prime }$ ) with three Euler angles ( $\unicode[STIX]{x1D711},\unicode[STIX]{x1D703},\unicode[STIX]{x1D713}$ ). Line ‘ON’ represents the pitch line of the $(x,y)$ and ( $x^{\prime },y^{\prime }$ ) coordinate planes. Two coordinate systems are overlapping initially. First the particle rotates around the $z^{\prime }$ axis with a recession angle $\unicode[STIX]{x1D711}$ and then the particle rotates around the new $x^{\prime }$ axis (i.e. line ‘ON’) with a nutation angle $\unicode[STIX]{x1D703}$ . Finally the particle rotates around the new $z^{\prime }$ axis with an angle of rotation $\unicode[STIX]{x1D713}$ .

2.2 Solid particle dynamics and fluid–solid boundary interaction

In our simulation, the ellipsoidal particle is described by

(2.5) $$\begin{eqnarray}\frac{x^{\prime 2}}{a^{2}}+\frac{y^{\prime 2}}{b^{2}}+\frac{z^{\prime ^{2}}}{c^{2}}=1,\end{eqnarray}$$

where $a$ , $b$ and $c$ are the lengths of the three semi-principal axes of the particle in the $x^{\prime }$ , $y^{\prime }$ and $z^{\prime }$ axis of a body-fixed coordinate system, respectively (see figure 1). The aspect ratio of the ellipsoidal particle is defined as $a/b$ . The body-fixed coordinate system can be obtained by a combination of coordinate transformation around the $z^{\prime }-x^{\prime }-z^{\prime }$ axis with Euler angles ( $\unicode[STIX]{x1D711},\unicode[STIX]{x1D703},\unicode[STIX]{x1D713}$ ) from the space-fixed coordinate system $(x,y,z)$ which initially overlaps the body-fixed coordinate system. The combination of coordinate transformation is illustrated in figure 1. The evolution axis always overlaps the $x^{\prime }$ direction. The migration and rotation of the particle are determined by the Newton equation and Euler equation, respectively,

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle m\frac{\text{d}\boldsymbol{U}(t)}{\text{d}t}=\boldsymbol{F}(t), & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D644}\boldsymbol{\cdot }\frac{\text{d}\unicode[STIX]{x1D6FA}(t)}{\text{d}t}+\unicode[STIX]{x1D6FA}(t)\times [\unicode[STIX]{x1D644}\boldsymbol{\cdot }\unicode[STIX]{x1D6FA}(t)]=\boldsymbol{T}(t), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the inertial tensor, and $\unicode[STIX]{x1D6FA}(t)$ and $\boldsymbol{T}(t)$ represent the angular velocity and the torque exerted on the particle in the body-fixed coordinate system, respectively. In the frame, $\unicode[STIX]{x1D644}$ is diagonal and the principal moments of inertia can be written as

(2.8a-c ) $$\begin{eqnarray}I_{x^{\prime }x^{\prime }}=m\frac{b^{2}+c^{2}}{5},\quad I_{y^{\prime }y^{\prime }}=m\frac{c^{2}+a^{2}}{5},\quad I_{z^{\prime }z^{\prime }}=m\frac{a^{2}+b^{2}}{5},\end{eqnarray}$$

where $m=4/3\unicode[STIX]{x1D70C}_{p}\unicode[STIX]{x03C0}abc$ is the mass of the particle and $\unicode[STIX]{x1D70C}_{p}$ is the density of the particle. It is not appropriate to solve (2.7) directly due to an inherent singularity (Qi Reference Qi1999). Thus four quaternion parameters are used as generalized coordinates to solve the corresponding system of equations (Qi & Luo Reference Qi and Luo2003). A coordinate transformation matrix with four quaternion parameters (Qi & Luo Reference Qi and Luo2003) is applied to transform corresponding items from the space-fixed coordinate system to the body-fixed coordinate system. With four quaternion parameters, equation (2.7) can be solved using a fourth-order-accurate Runge–Kutta integration procedure (Huang et al. Reference Huang, Yang, Krafczyk and Lu2012b ).

In the simulations, the fluid–solid boundary interaction is based on the schemes of Aidun et al. (Reference Aidun, Lu and Ding1998) and Lallemand & Luo (Reference Lallemand and Luo2003). The accurate moving-boundary treatment proposed by Lallemand & Luo (Reference Lallemand and Luo2003) is applied to solve the problem caused by the moving curved-wall boundary condition of the ellipsoid.

The general schemes for calculation of interactive force between fluids and particles in the LBM include stress integration, momentum exchange and volume fraction models, which has been summarized and analysed by Chen et al. (Reference Chen, Cai, Xia, Wang and Chen2013). The stress integration scheme may be good but it is not so efficient (Chen et al. Reference Chen, Cai, Xia, Wang and Chen2013). Here the momentum exchange scheme is used to calculate the force exerted on the solid boundary, which is accurate and efficient (Chen et al. Reference Chen, Cai, Xia, Wang and Chen2013). The forces due to the fluid nodes covered by the solid nodes and the solid nodes covered by the fluid nodes (Aidun et al. Reference Aidun, Lu and Ding1998) are also considered in the study.

To prevent overlap of the particle and the wall, usually the repulsive force between the wall and particle should be applied (Huang et al. Reference Huang, Yang and Lu2014). Here the lubrication force model is identical to that we used in Huang et al. (Reference Huang, Yang and Lu2014) and the validation of the force model has been tested extensively by Huang et al. (Reference Huang, Yang and Lu2014).

Our previous study on the tumbling mode of an ellipsoidal particle suspended in shear flow partially validated our three-dimensional LBM code (Huang, Wu & Lu Reference Huang, Wu and Lu2012a ) because the rotational period or orbit is very consistent with Jeffery’s analytical solution (Jeffery Reference Jeffery1922). The LBM code has also been validated by Yang, Huang & Lu (Reference Yang, Huang and Lu2015) for the case of migration of a neutrally buoyant sphere in a tube Poiseuille flow. The lift force, angular and migration velocities agree well with those in Yang et al. (Reference Yang, Wang, Joseph, Hu, Pan and Glowinski2005). Here three more cases are simulated to validate our LBM code. The first one concerns migrations of a neutrally buoyant sphere in tube Poiseuille flows (Karnis et al. Reference Karnis, Goldsmith and Mason1966). The second deals with an ellipsoid particle sedimentation in a vertical circular tube under gravity. The third one describes the rotation of a prolate spheroid in shear flow. The result is shown in appendix A. The accuracy of the LBM is also investigated in § 3.1.

3.1 Accuracy of the LBM

The accuracy of the LBM is investigated here. Grid independence and time-step independence studies are performed. As an example, cases of a prolate ellipsoid with $Re=162$ and $D/A=2$ were simulated. The key parameters for the cases are shown in table 1, where $\unicode[STIX]{x0394}t^{\ast }$ is a normalized time and $\unicode[STIX]{x0394}p$ is the pressure difference between the two ends of the tube. The corresponding results for grid independence and time-step independence are shown in figure 3. In figure 3(a), the grid size is labelled, e.g. the legend ‘ $A=40lu$ ’ means 40 grids are used to discretize the particle’s major axis. From figure 3(a), it is seen that the result of $A=60lu$ is very close to that of $A=72lu$ . However, for the coarse mesh, i.e.  $A=40lu$ , the period of the rotation has a significant discrepancy with respect to that of $A=72lu$ . The grid size with $A=60lu$ seems small enough to obtain an accurate result. Hence, for cases with $Re\approx O(100)$ , the grid size with $A=60lu$ is used.

For the time-step independence study, three cases with $\unicode[STIX]{x1D70F}=0.575,0.6,0.64$ were simulated (see table 1). The corresponding $\unicode[STIX]{x0394}t^{\ast }=\unicode[STIX]{x0394}t\unicode[STIX]{x1D708}/A^{2}=6.94\times 10^{-6},9.26\times 10^{-6},1.30\times 10^{-5}$ , respectively. It is seen from figure 3(b) that the curve for $\unicode[STIX]{x1D70F}=0.6$ is very close to that for $\unicode[STIX]{x1D70F}=0.575$ . Hence, the time step with $\unicode[STIX]{x1D70F}=0.6$ is small enough to get accurate results.

Figure 3. Effect of grid size (a) and time step (b) on the orientation of the prolate ellipsoid ( $\cos \unicode[STIX]{x1D6FE}$ ). The parameters for the cases are shown in table 1. Log–log plots of the error described in (3.2) as a function of grid size (c) and time step (d), respectively.

To test the accuracy of the present numerical method, the case with the finest mesh $A=90lu$ and the smallest time step $\unicode[STIX]{x0394}t^{\ast }=3.70\times 10^{-6}$ (Case C) is also simulated. The result is referred to as the accurate result. The $\cos \unicode[STIX]{x1D6FE}$ of Case A3 is very close to that of Case C. Since the error can vary with time, a time average over one period is required. After the tumbling state of the ellipsoid reaches a stable periodic state at $t_{1}$ with a period of $T$ , the relative error is defined as

(3.2) $$\begin{eqnarray}E_{2}=\frac{\displaystyle \sqrt{\int _{t_{1}}^{t_{1}+T}[\cos \unicode[STIX]{x1D6FE}(t)-\cos \unicode[STIX]{x1D6FE}_{0}(t)]^{2}\,\text{d}t}}{\displaystyle \sqrt{\int _{t_{1}}^{t_{1}+T}[\cos \unicode[STIX]{x1D6FE}_{0}(t)]^{2}\,\text{d}t}},\end{eqnarray}$$

where $\cos \unicode[STIX]{x1D6FE}_{0}(t)$ is the accurate result.

The least-squares data fit in figures 3(c) and 3(d) shows that the fitted slopes are 2.06 and 1.93, respectively. Hence, the present LBM solver has an approximately second-order accuracy in both space and time. That is consistent with the conclusion on accuracy of the LBM in Mei et al. (Reference Mei, Luo and Shyy1999).

5.1 Particle’s rotational energy and inertial effect of the fluid

In this section we mainly discuss both inertial effects of the fluid and the particle. First, the effect of fluid inertia on the particle behaviour is discussed. Figure 6 shows $|\text{cos}\,\unicode[STIX]{x1D6FE}|_{max}$ as a function of  $Re$ . $\unicode[STIX]{x1D6FE}$ is the angle between the $x^{\prime }$ -axis and the $z$ -direction. $\cos (\unicode[STIX]{x1D6FE})_{max}$ represents the maximum $\cos \unicode[STIX]{x1D6FE}$ in one period in the kayaking mode (see figure 5 b). It quantifies the extent of deviation from the log-rolling state. When the particle is in the log-rolling state, $|\text{cos}\,\unicode[STIX]{x1D6FE}|_{max}$ should be equal to zero because the major axis ( $x^{\prime }$ -axis) is perpendicular to the $z$ -direction.

It is seen from figure 6 that for cases $D/A=1.7$ , $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ , $|\text{cos}\,\unicode[STIX]{x1D6FE}|_{max}$ is not continuous as a result of the mode transition. When $Re$  is less than 80, $|\text{cos}\,\unicode[STIX]{x1D6FE}|_{max}\approx 0.9$ , the ellipsoid is in the kayaking state. As $Re$  increases above 90, the ellipsoid adopts the log-rolling mode. Hence, as the inertia of the fluid increases, the kayaking mode may transit to the log-rolling mode. This character can be understood as follows. When $Re$  increases, both the shear stress acting on the particle and the inertia of the fluid increase. The larger shear stress on the particle increases the angular velocity of the particle ( $\unicode[STIX]{x1D714}_{x^{\prime }}$ ) around the $x^{\prime }$ -axis (see figure 6). The energy of the rotational body is $(I_{x^{\prime }x^{\prime }}\unicode[STIX]{x1D714}_{x^{\prime }}^{2})/2$ . Hence the rotational energy increases. On the other hand, the larger fluid inertia may prevent the nutation or the flapping-like movement. Both the increasing rotational energy and fluid inertia are believed to assist the particle reach the log-rolling state without nutation. Although $\unicode[STIX]{x1D714}_{x^{\prime }}^{\ast }$ increases with $Re$ , the equilibrium radial position $r^{\ast }=r/R$ is approximately 0.46, and increases only slightly with  $Re$  (not shown).

Figure 6. Transition from the kayaking mode (the shaded area) to the log-rolling mode. $|\text{cos}\,\unicode[STIX]{x1D6FE}|_{max}$ and normalized angular velocity $\unicode[STIX]{x1D714}_{x^{\prime }}^{\ast }$ as functions of $Re$  for $D/A=1.7$ , $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ .

Figure 7. Phase diagram for a prolate ellipsoid ( $a/b=4$ ) inside Poiseuille flow with $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ .

The ellipsoidal particle with aspect ratio $a/b=4$ is also considered in our study. The phase diagram for a prolate ellipsoid with $a/b=4$ inside Poiseuille flow is shown in figure 7. For $D/A>1.9$ , the tumbling mode still occurs. For $1<D/A<1.9$ , there is only the kayaking mode and the log-rolling mode is not observed. The significant difference between figures 4 and 7 lies in the log-rolling mode. The disappearance of the log-rolling mode in figure 7 may be attributed to the significantly smaller rotational energy in cases with $a/b=4$ . A typical example is presented below. Suppose two cases $D1$ and $D2$ have identical key parameters $D/A=1.5$ , $Re=78.4$ , $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ , $a=30lu$ , and $\unicode[STIX]{x1D70F}=0.6$ , but in Case $D1$ , $b=c=a/2$ and in Case $D2$ , $b=c=a/4$ . That means in Cases $D1$ and $D2$ , the particle aspect ratios are 2 and 4, respectively. It is seen from figures 4 and 7 that the particles in Cases $D1$ and $D2$ adopt the log-rolling and kayaking modes, respectively. The behaviour may be understood as follows. In Case $D2$ the particle is more rod-like and the principal moment of inertia $I_{x^{\prime }x^{\prime }}$ in Case $D2$ is only $1/16$ of that in Case $D1$ because $I_{x^{\prime }x^{\prime }}=m((b^{2}+c^{2})/5)=4/3\unicode[STIX]{x1D70C}_{p}\unicode[STIX]{x03C0}abc((b^{2}+c^{2})/5)$ . In the equilibrium state, the normalized average angular velocities $\unicode[STIX]{x1D714}_{x^{\prime }}^{\ast }=(\unicode[STIX]{x1D714}_{x^{\prime }}A^{2})/\unicode[STIX]{x1D708}$ for Cases $D1$ and $D2$ are 36.5 and 22.3, respectively. It is seen that the energy of the rotational body $(I_{x^{\prime }x^{\prime }}\unicode[STIX]{x1D714}_{x^{\prime }}^{2})/2$ in Case $D2$ would be much smaller than that in Case $D1$ . A smaller energy of the rotational body may induce nutation motion when it rotates around the $x^{\prime }$ -axis and interacts with fluid. Hence, the particle in Case $D2$ more easily adopts the kayaking state instead of the log-rolling state. Again, we see that the increased rotational energy indeed assists the particle reaching the log-rolling state without nutation.

5.2 Inertial effect of the particle

To investigate the inertial effect of the particle, twelve cases listed in table 3 were simulated. They are classified into four groups. Here $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ is not limited to be unity, i.e. the cases with different $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ were performed. Group I represents a wide tube case with $D/A=3$ . Groups III and IV denote narrow tube cases ( $D/A=1.5$ ). Group II denotes cases with $D/A=1.9$ . The results of these cases with different density ratios are presented in table 3 and figure 8(b,c).

Table 3. The effect of the inertia of the particle, $r^{\ast }=r/R$ is the equilibrium radial position. $|\text{cos}\,\unicode[STIX]{x1D6FE}|_{max}$ is only applicable for the kayaking mode.

In Group I, only the tumbling state is observed. The normalized tumbling rotational angular velocity increases with the density ratio $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ . The particle’s lateral migration is slightly closer to the wall as $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ increases. Specifically, for the case with $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=3.0$ , the lateral migration is periodically oscillating and $r^{\ast }\in (0.566,0.570)$ , i.e. the centre of mass of the ellipsoid takes a wavy trajectory in the radial direction inside an axi-symmetric plane.

Group II in table 3 shows that for the tube $D/A=1.9$ , again only the tumbling mode is reproduced. In the tumbling mode, the radial position of the particle’s centre of mass is continuously oscillating periodically. The magnitude of the oscillation increases with the density ratio $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ . It is seen that when $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=0.3$ , the average $r^{\ast }=0.442$ with $r^{\ast }\in (0.412,0.470)$ , the oscillation magnitude is $\unicode[STIX]{x1D6FF}r=0.470-0.412=0.058$ . When $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ and $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=3.0$ , it oscillates with $\unicode[STIX]{x1D6FF}r=0.084$ and $\unicode[STIX]{x1D6FF}r=0.163$ , respectively, in the radial direction. The oscillation in the radial direction is attributed to the particle’s inertia and the strong interaction between the particle and the tube wall. For a specific narrow tube, the magnitude of the radial oscillation increases with the inertia of the particle.

Figure 8. (a) Phase diagram in a three-dimensional parameter space. (b,c) Phase diagram in the planes $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=3.0$ and $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=0.3$ , respectively.

Figure 9. Phase diagram for an oblate ellipsoid inside Poiseuille flow ( $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ ).

Hence, for $D/A\geqslant 1.9$ , the motion mode is not affected by the particle’s inertia ( $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}\leqslant 3$ ). Figure 8 shows the phase diagram in a three-dimensional parameter space. As we do not intend to provide accurate borders between the modes, it is only a schematic diagram. From figure 8, it is seen that the tumbling mode is still in the right part of the space ( $D/A\geqslant 1.9$ ), which is independent of both $Re$  and $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ under our considered parameter space.

Groups III and IV in table 3 are the cases in narrower tubes but with lower and higher $Re$ , respectively. For the lower $Re$  cases (Group III), in the narrower tubes ( $D/A=1.5$ ) the particle adopts the log-rolling mode when $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=0.3$ and 0.6, while $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1.0$ , 2.0 and 3.0, the state is kayaking. For higher $Re$  cases (Group IV), the particle adopts the log-rolling mode at $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=0.3$ , 0.6 and 1.0, but when $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=2.0$ and 3.0, it still takes the kayaking mode. Hence, for a constant $Re$ , the kayaking mode is preferred when the particle’s inertia increases.

Hence, for $D/A<1.9$ , the mode distribution would be affected by the particle’s inertia ( $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ ). The border separating the kayaking and log-rolling modes may move upwards or downwards depending on the particle’s inertia. Figure 8 shows the approximate result. If the planes with constant $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ are viewed from front to back, the border moves upwards. It means that on a plane with higher $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ , the area of the kayaking mode increases. In other words, due to the inertial effect, the kayaking mode becomes more common.

6.1 The particle inertial effect

The approximate phase diagram in a three-dimensional parameter space for an oblate ellipsoid is shown in figure 14. The parabolic cylinder consists of a cluster of parabolic curves. Inside the parabolic cylinder there is the intermediate mode. It is seen that the mode distribution would be affected by the particle’s inertia ( $\unicode[STIX]{x1D70C}_{p}$ ). For wider tubes, the log-rolling mode is dominant in the phase diagram. For narrower tubes, the border separating the intermediate and the other two modes may move upwards or downwards depending on the particle’s inertia. If the planes with constant $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ are viewed from front to back, the border moves upwards. It means that on a plane with higher $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ , the area of the log-rolling mode increases. For $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=3$ , the inclined mode almost disappears under our considered parameter space and the log-rolling mode occupies a greater area in the $D/A-Re$ plane. In other words, at a slightly higher $Re$ , the inertia of the oblate particle helps it maintain the log-rolling state.

Figure 14. Phase diagram for an oblate ellipsoid. (a) Phase diagram in a three-dimensional parameter space. The vertical plane ABCD separates the inclined and log-rolling modes. The rotational mode appears inside the parabolic cylinder is the intermediate mode. The region of the inclined mode is bounded by the plane ABCD, the curved plane BCEF, and the ‘wall’ of the box. (b,c) Phase diagram in the planes $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=3.0$ and $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=0.3$ , respectively.

6.2 The inclined mode

In the following, the reasons for the existence of the inclined mode are further discussed. The streamlines around the particle and pressure contours on an oblate ellipsoid are illustrated in figure 15. The pressure is normalized by $p=(p-c_{s}^{2}\unicode[STIX]{x1D70C}_{f})/\unicode[STIX]{x1D70C}_{f}U_{m}^{2}$ . In the case for $D/A=2$ and $Re=354$ , the dimensional parameters are $\unicode[STIX]{x1D70F}=0.565$ , tube length $L=1920lu$ and diameter $D=160lu$ . The oblate ellipsoid has dimensions $a=20lu$ , $b=c=40lu$ . Under the pressure difference $\unicode[STIX]{x0394}p=2.49\times 10^{-3}~mu/(lu\cdot ts^{2})$ between the two ends of the tube, the maximum velocity in the Poiseuille flow is $U_{m}=0.0958lu/ts$ and terminally the particle reaches the inclined mode with a constant velocity $U_{p}=-0.0559lu/ts$ in the $z$ -direction without rotation.

Figure 15. Streamlines and pressure contours on the surface of an oblate ellipsoid in Poiseuille flow. The case is $D/A=2$ , $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=1$ , and $Re=354$ . The position of the centre of the particle is at $x=0$ and $y=0.468$ . The streamlines are drawn in the inertial frame moving with $U_{p}$ , which is the particle’s velocity in the $z$ -direction in the space-fixed coordinate system.

Figure 15(a) shows the streamlines inside the axi-symmetric plane (plane ABCD, i.e.  $(x=0)$ -plane) where the $x^{\prime }$ -axis stays. The angle between the $x^{\prime }$ -axis and $z$ -axis is approximately $67^{\circ }$ because $\cos (\unicode[STIX]{x1D6FE})\approx 0.38$ . The streamlines are drawn in an inertial frame moving with the same velocity as that of the particle. Figure 15(b) shows the pressure contours and surrounding streamlines viewed from the $y$ -direction. It is seen that the contours and streamlines are almost symmetric about the $(x=0)$ -plane.

We can qualitatively understand why the particle is able to maintain the steady state. In figure 15(a), on the lower and upper surface of the particle, the particle experienced lower and higher pressure, respectively. That may generate a positive torque (counter-clockwise). On the other hand, from the directions of the streamlines, it is seen that on the whole the shear stress acting on the surface would generate a negative torque (clockwise). Hence the torque due to the pressure difference will be balanced by the negative torque due to the viscous force acting on the surface. In this way, the inclined mode is able to be an equilibrium state.

Because both the log-rolling mode and the inclined mode are stable modes, an intermediate state combining the characteristics of both these stable modes should also possibly exist.