1. Introduction
Deformable particles, such as fluid drops, elastic particles, capsules, vesicles and red blood cells, suspended in channel flows have long been known to migrate away from the channel wall at low Reynolds numbers (
${Re}$
), with few exceptions for fluid drops with a certain range of the viscosity ratio (Karnis, Goldsmith & Mason Reference Karnis, Goldsmith and Mason1963; Goldsmith Reference Goldsmith1971; Leal Reference Leal1980; Sekimoto & Leibler Reference Sekimoto and Leibler1993; Mortazavi & Tryggvason Reference Mortazavi and Tryggvason2000; Hodges, Jensen & Rallison Reference Hodges, Jensen and Rallison2004; Shin & Sung Reference Shin and Sung2011; Villone & Maffettone Reference Villone and Maffettone2019; Rezghi, Li & Zhang Reference Rezghi, Li and Zhang2022; Takeishi et al. Reference Takeishi, Ishimoto, Yokoyama and Rosti2025). The migration of deformable particles away from the wall, sometimes called the axial accumulation, is understood to be caused by nonlinear effects between particle deformation and flow field due to the presence of the wall (Leal Reference Leal1980; Takemura, Magnaudet & Dimitrakopoulos Reference Takemura, Magnaudet and Dimitrakopoulos2009; Sugiyama & Takemura Reference Sugiyama and Takemura2010). Along with the wall effect, the curvature or shear gradient of the velocity profile also produces a tendency for migration towards the centreline, where the shear is lower. Kaoui et al. (Reference Kaoui, Ristow, Cantat, Misbah and Zimmermann2008) demonstrated by a two-dimensional (2-D) numerical simulation in the low-
${Re}$
limit that a vesicle in unbounded Poiseuille flow migrates towards the centre of the flow, as a result of asymmetric deformation between higher and lower shear sides. In the Stokes flow limit, Villone et al. (Reference Villone, Greco, Hulsen and Maffettone2016) and Villone (Reference Villone2019) showed by numerical analyses that a neo-Hookean elastic particle always migrates towards the channel centreline both in circular and square channel flows. They also found that the migration velocity is larger for particles that are more deformable and located closer to the channel wall, due to larger deformation.
In recent years, inertial migration of particles in rectangular channel flows at finite
${Re}$
has been intensively studied in the field of microfluidics, as it is associated with a promising technology for the passive manipulation, focusing and sorting of biological cells and particles (Di Carlo Reference Di Carlo2009; Martel & Toner Reference Martel and Toner2014; Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2019; Kalyan et al. Reference Kalyan, Torabi, Khoo, Sung, Choi, Wang, Treutler, Kim and Hur2021; Lee, Kim & Yang Reference Lee, Kim and Yang2023; Zhang et al. Reference Zhang2024). Currently, most microfluidic devices for the separation and sorting of biological cells are designed to work mainly based on the difference in cell size. However, the operation based on their deformability is useful and of practical importance for improving diagnostic, therapeutic, biological and other performance, since many target cells, such as circulating tumour cells and malaria-infected cells, have a different stiffness from surrounding cells, and various diseases, such as cancer, blood diseases and inflammation, often accompany cell deformability alterations (Hou et al. Reference Hou, Bhagat, Chong, Mao, Tan, Han and Lim2010; Hur et al. Reference Hur, Henderson-MacLennan, McCabe and Di Carlo2011; McFaul, Lin & Ma Reference McFaul, Lin and Ma2012; Yang et al. Reference Yang, Lee, Ahn, Kang, Shim, Lee, Hyun and Kim2012; Preira et al. Reference Preira, Grandne, Forel, Gabriele, Camara and Theodoly2013; Wang et al. Reference Wang, Mao, Byler, Patel, Henegar, Alexeev and Sulchek2013; Holmes et al. Reference Holmes, Whyte, Bailey, Vergara-Irigaray, Ekpenyong, Guck and Duke2014; Park et al. Reference Park2016; Krüger et al. Reference Krüger, Holmes and Coveney2014a
; Guo et al. Reference Guo, Duffy, Matthews, Islamzada and Ma2017; Connolly, McGourty & Newport Reference Connolly, McGourty and Newport2021; Stathoulopoulos et al. Reference Stathoulopoulos, Passos, Kaliviotis and Balabani2024). Despite the importance of the deformability of suspended particles, the fundamentals of its effect on inertial migration have not been fully understood.
The present study aims to investigate the impact of the deformability of suspended particles on their inertial migration experimentally and numerically from fluid-dynamical points of view. To simplify the problem, we treat the flow of dilute suspensions of deformable particles (hydrogel microspheres in the experiment and viscous hyperelastic particles in the numerical simulation) in a Newtonian fluid flowing through straight square channels.
Inertial migration of particles suspended in channel flows was first observed more than 60 years ago by Segré & Silberberg (Reference Segré and Silberberg1961). They reported that neutrally buoyant spherical particles in circular tube flows migrate laterally towards an annulus with a radius 0.6 times the tube at relatively low
${Re}$
(Segré & Silberberg Reference Segré and Silberberg1961, Reference Segré and Silberberg1962). This phenomenon is referred to as the Segre–Silberberg (SS) effect and the particle focusing annulus is called the SS annulus. On the other hand, in rectangular channel flows, which are commonly used in microfluidics, suspended particles were found to first focus on a ring, called pseudo-SS ring (pSS ring), and then migrate towards several discrete points in downstream cross-sections (Chun & Ladd Reference Chun and Ladd2006; Choi, Seo & Lee Reference Choi, Seo and Lee2011; Nakagawa et al. Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015; Shichi et al. Reference Shichi, Yamashita, Seki, Itano and Sugihara-Seki2017). In particular, spherical particles in square channel flows eventually focus on four points located near the centre of the channel faces at relatively low
${Re}$
(Di Carlo et al. Reference Di Carlo, Irimia, Tompkins and Toner2007, Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Zhou & Papautsky Reference Zhou and Papautsky2013; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014; Miura, Itano & Sugihara-Seki Reference Miura, Itano and Sugihara-Seki2014). This focusing position is called the channel face equilibrium position or midline equilibrium position (MEP). If the focusing positions and cross-stream migration velocities can be controlled by the size, shape and deformability of suspended particles, the SS effect can be used for continuous particle separation and sorting in a label-free, external field-free manner using microchannels of simple geometry. To this end, extensive studies have been performed to develop microfluidic devices applying the SS effect to suspension flows of living cells and particles (see the review articles such as those by Di Carlo Reference Di Carlo2009; Bhagat et al. Reference Bhagat, Bow, Hou, Tan, Han and Lim2010; Karimi, Yazdi & Ardekani Reference Karimi, Yazdi and Ardekani2013; Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2019; Razavi Bazaz et al. Reference Razavi Bazaz, Mashhadian, Ehsani, Saha, Krüger and Ebrahimi Warkiani2020; Tang et al. Reference Tang, Zhu, Jiang, Zhu, Yang and Xiang2020; Kalyan et al. Reference Kalyan, Torabi, Khoo, Sung, Choi, Wang, Treutler, Kim and Hur2021; Lee et al. Reference Lee, Kim and Yang2023; Zhang et al. Reference Zhang2024, and the references therein).
The inertial focusing of spherical particles on the midlines (MEP) can be explained by the inertial lift exerted on them in square channel flows. The inertial lift is known to consist mainly of the shear gradient-induced lift, acting in the direction towards higher shear, and the wall effect, acting in the direction away from the channel wall (Ho & Leal Reference Ho and Leal1974; Schonberg & Hinch Reference Schonberg and Hinch1989; Asmolov Reference Asmolov1999; Matas, Morris & Guazzelli Reference Matas, Morris and Guazzelli2004). In circular tube flows, the balance between the outward lift due to the shear gradient and the inward lift due to the wall effect determines the radius of the SS annulus (Matas, Morris & Guazzelli Reference Matas, Morris and Guazzelli2009). In the pressure-driven flow through a square channel, on the other hand, the shear rate is higher closer to the channel walls with the maximum at the centre of the channel faces, and it is lower near the centre of the cross-section and the corners. Thus, the direction of the inertial lift due to the shear gradient is mostly outwards in the radial direction and towards the midlines, or away from the diagonal, in the azimuthal direction. The radial component, which is more pronounced than the azimuthal component, has a tendency to balance with the wall effect, so that particles approach the pSS ring first, and then migrate towards the
$y$
- or
$z$
-axis (MEP) due to the presence of the azimuthal component of the shear gradient-induced lift towards these axes.
In addition to these inertial lifts, deformable particles in channel flows experience lifts due to particle deformation (deformation-induced lift) as already noted. Thus, the inertial focusing of deformable particles is expected to be different from that of rigid particles. In fact, polyethylene glycol and alginate hydrogel particles suspended in a rectangular channel flow were demonstrated to focus closer to the channel centreline as they become more deformable (Ding et al. Reference Ding, Xiao, Zhang, Zhang, Zhu, Wang and Xue2025). Recent numerical studies have reported several features of inertial migration specific to deformable particles in square channel flows. Schaaf & Stark (Reference Schaaf and Stark2017) showed that most of the elastic capsules migrate towards equilibrium positions on the diagonal of the channel cross-section and their focusing positions become closer to the centre for more deformable particles. The focusing position on the diagonal is referred to here as the diagonal equilibrium position (DEP) and that on the channel centre as the centre equilibrium position (CEP). Raffiee, Dabiri & Ardekani (Reference Raffiee, Dabiri and Ardekani2017) demonstrated that deformable capsules focus on the DEP in a relatively wide range of
${Re}$
, in contrast to rigid particles, which focus on the MEP in the same
${Re}$
range. Using hyperelastic particles, Esposito et al. (Reference Esposito, Romano, Hulsen, D’Avino and Villone2022) have recently reported the existence of another equilibrium position, which is located at an intermediate position between the MEP and DEP. In certain ranges of parameters, this intermediate equilibrium position (IEP) is stable whereas both MEP and DEP are unstable.
All these numerical studies agree that equilibrium positions of deformable particles can exist on the diagonal in square channel flows at certain parameter values, but to the authors’ knowledge, no experimental evidence for the presence of the DEP has been provided, except for our recent study on human red blood cells (RBCs) (Tanaka & Sugihara-Seki Reference Tanaka and Sugihara-Seki2022). Human RBCs have a biconcave discoid shape at rest, with a diameter, thickness and volume of approximately 8,
$2{-}3\,\unicode{x03BC} \textrm{m}$
and
$94\,\unicode{x03BC} \textrm{m}^3(=V)$
, respectively (Hochmuth Reference Hochmuth1987). Their equivalent diameter is
$d=2\sqrt [3]{3V/4\pi } \sim 5.6\,\unicode{x03BC} \textrm {m}$
. RBCs are surrounded by a thin membrane, inside which is a haemoglobin solution. Due to this structure and the excess surface area relative to the volume, they are highly deformable. Tanaka & Sugihara-Seki (Reference Tanaka and Sugihara-Seki2022) demonstrated that at low
${Re}$
with negligible inertia, RBCs migrate towards the channel centre (CEP) in square channel flows, in accord with previous studies (Goldsmith Reference Goldsmith1971; Mchedlishvili & Maeda Reference Mchedlishvili and Maeda2001; Sasaki et al. Reference Sasaki, Seki, Itano and Sugihara-Seki2017; Losserand, Coupier & Podgorski Reference Losserand, Coupier and Podgorski2019). In contrast, at finite
${Re}$
, they focus near four points on the diagonals, indicating the presence of the DEP. Since RBCs have a rather specific shape, biconcave discoid, our first aim of the present study is to demonstrate experimentally the presence of the DEP using initially spherical deformable particles. To this end, using hydrogel spheres with
$5.2\,\unicode{x03BC} \textrm {m}$
diameter and square channels with
$50\,\unicode{x03BC} \textrm {m}$
width, we experimentally investigate their inertial focusing in square channel flows. The obtained particle distributions in the channel cross-section are compared with those of RBCs or rigid spherical particles with similar sizes to explore the effect of particle deformation on their migration in square channel flows.
These hydrogel particles with a diameter of
$5.2\,\unicode{x03BC} \textrm {m}$
are the largest of their kind that can currently be synthesised and the square channels used in this study,
$50\,\unicode{x03BC} \textrm {m}$
wide (
$600\,\textrm {mm}$
long), are the smallest commercially available. Thus, it is difficult to perform further experiments with different parameters, since suspended particles need to have a finite size
$({\gtrsim}0.07)$
relative to the channel width for inertial focusing (Di Carlo et al. Reference Di Carlo, Irimia, Tompkins and Toner2007; Bhagat et al. Reference Bhagat, Kuntaegowdanahalli and Papautsky2008). Instead, we conduct numerical simulations for a viscous hyperelastic particle suspended in square channel flows, using a full Eulerian method (Sugiyama et al. Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011). By comparing the numerical and experimental results, we estimate the stiffness of the hydrogel particles and compare the obtained value with reported Young’s moduli for similar particles. Next, the impact of the deformability of the particles on their inertial focusing is investigated numerically by systematically changing the Reynolds number
${Re}$
and the capillary number
$Ca$
(or the Laplace number,
$La$
, defined as the ratio of the particle Reynolds number and the capillary number; see § 2.1). Several types of focusing patterns of the particles in the cross-section are obtained. As the particles become more deformable, the focusing pattern shifts from four-point focusing on the MEP, to eight-point focusing on the MEP and DEP, to four-point focusing on the DEP, and to single-point focusing on the CEP. These transitions of the particle focusing pattern are discussed in terms of the interplay between the inertial lift and the particle deformation-induced lift. In contrast to a previous numerical study (Esposito et al. Reference Esposito, Romano, Hulsen, D’Avino and Villone2022), we observe a bistable state of the MEP and DEP, but no single focusing on the IEP.
The present numerical computations indicate that the stability of the channel centre (CEP) is changed at
$La\sim 8$
for the blockage ratio of 0.2; it is stable for
$La\lt 8$
and unstable for
$La\gt 8$
. This numerical result is accounted for by the lateral force exerted on the hyperelastic particle near the channel centre under the assumption of small
${Re}$
and
$Ca$
. The discussion is based on the asymptotic expansion of the lateral force in terms of
${Re}$
and
$Ca$
, together with the numerical results for rigid spheres at
$Ca=0$
and for hyperelastic particles at
${Re}=0$
.
In this manuscript, the experimental and numerical methods are described in § 2. The experimental results using hydrogel particles are shown in § 3.1. The numerical results, including the comparison with the experimental results and discussion, are presented in § 3.2. Section 3.3 is devoted to the discussion of the stability of the CEP based on the inertial lift and deformation-induced lift. Concluding remarks are drawn in § 4. In Appendix A, the details of the numerical method for a rigid sphere using the immersed boundary method are described.
2. Methods
2.1. Dimensionless parameters
If suspended particles are neutrally buoyant rigid spherical particles with diameter
$d$
and the suspending fluid is a Newtonian fluid with density
$\rho$
and viscosity
$\mu$
, then important dimensionless parameters for the inertial migration in square channel flows are the Reynolds number
${Re}=\rho UD/\mu$
and the blockage ratio
$\kappa =d/D$
, where
$D$
and
$U$
represent the channel width and the average flow velocity, respectively. The particle Reynolds number is defined as
${\textit{Re}}_p={Re}\boldsymbol{\cdot }\kappa ^2$
. We consider here the case in which the concentration of particles is low enough that particle–particle interactions can be neglected. In the case of elastic particles with elastic modulus
$G$
, the capillary number
$Ca=\mu U/DG$
, qualifying the ratio between the viscous stress and the elastic stress, is added as a control parameter, and
$d$
represents the undeformed diameter of the particles. The ratio between
${\textit{Re}}_p$
and
$Ca$
gives the Laplace number
$La={\textit{Re}}_p/Ca=\rho d^2G/\mu ^2$
, which is expressed in terms of the fluid and particle properties, independent of the flow velocity and channel width.
2.2. Materials and experimental methods
Figure 1. (a) Chemical structure of hydrogel microspheres and (b) optical microscopy image of packed hydrogel microspheres (diameter
$d=5.2\,\unicode{x03BC} \textrm {m}$
). NIPAm,
$N$
-isopropylacrylamide; BIS,
$N,{}N'$
-methylenebis(acrylamide); AAc, acrylic acid.
As deformable suspended particles, we used poly(
$N$
-isopropylacrylamide)-based hydrogel microspheres fluorescently labelled by 5-aminofluorescein (figure 1) (Kawamoto et al. Reference Kawamoto, Yanagi, Nishizawa, Minato and Suzuki2023, Reference Kawamoto, Minato and Suzuki2024). These hydrogel particles were synthesised by a modified aqueous precipitation polymerisation method under the condition of
$N$
-isopropylacrylamide (NIPAm) 70 mol %,
$N,N'$
-methylenebis (acrylamide) (BIS) 1 mol % and acrylic acid (AAc) 29 mol %. As shown in figure 1(b), the particle size is almost uniform, with an average diameter
$d= 5.2\,\unicode{x03BC} \textrm {m}$
. Details of these hydrogel microspheres and their deformation at the air/water interface are described by Minato et al. (Reference Minato, Murai, Watanabe, Matsui, Takizawa, Kureha and Suzuki2018).
Figure 2. (a) Experimental set-up and (b) method of image analyses: an example of obtained images (top), distribution of particle centroids (middle) and probability density function (p.d.f.) in the channel cross-section (bottom).
The hydrogel particles were suspended in glycerol aqueous solution at volume fractions of
${0.003}{-}0.01\,\%$
. The density and viscosity of the solution are
$1.05 \times 10^3 {\,\rm kg\,m}^{-3}$
and 1.72 mPa s, respectively, at
$22\,{}^{\circ }\textrm {C}$
. As shown in figure 2(a), a syringe pump (Nexus 6000, ISIS) was used to infuse the suspension into a straight glass tube with a square cross-section of width
$D= 50\,\unicode{x03BC} \textrm {m}$
and length
$L=50\!-\!600\,\textrm {mm}$
(VitroCom). Face-on fluorescence images of the tube cross-section near the outlet (1–2 mm upstream of the outlet) were taken along the tube axis from the downstream side using a high-speed camera with a built-in image intensifier (SV200i, Photron) under the illumination of a 3 W blue laser (Kentech), equipped with a longpass filter (SCHOTT OG-530, Edmund Optics) and an ultralong working distance objective (SLMPLN
$50\times$
,
$100\times$
, Olympus). An LED light was also used to detect the tube wall. When using a
$50\times$
objective, the pixel size is
$0.268\times \, 0.268\,\unicode{x03BC} \textrm {m}^2$
. The images obtained were analysed using the public domain software ImageJ (NIH) to detect the position of the centroid of each particle in the cross-section. In each experiment, we counted more than 300 particles to determine their distribution and the probability densities (figure 2
b).
2.3. Numerical methods
Figure 3. Configuration for the numerical simulation
$(L=2D)$
.
In our numerical simulations, we considered a single deformable particle flowing in a channel of length
$L$
with a square cross-section with side
$D$
as shown in figure 3. The channel is filled with a Newtonian fluid and the pressure-driven Poiseuille flow is directed along the
$x$
-coordinate. We assume that both the fluid and the particle are incompressible and have the same density
$\rho$
and dynamic viscosity
$\mu$
, which means the particle is neutrally buoyant. The deformable hydrogel particle is modelled as a viscous hyperelastic solid rather than a capsule or a fluid drop. We specifically employ the neo-Hookean model (Mooney Reference Mooney1940; Rivlin Reference Rivlin1948; van Hoogstraten, Slaats & Baaijens Reference van Hoogstraten, Slaats and Baaijens1994), which has the simplest constitutive law for such hyperelastic materials. The particle is initially spherical with a diameter
$d$
under unstressed conditions.
This study explores a broad parameter space to achieve more comprehensive understandings rather than to reproduce specific phenomena. We employed two dimensionless numbers
${Re}$
and
$Ca$
as the governing parameters.
Our simulation employs the fully Eulerian method developed by Sugiyama et al. (Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011). This method is a finite difference method based on the VOF method by Hirt & Nichols (Reference Hirt and Nichols1981), solving the fluid–structure coupling problem using the VOF function. To derive the nondimensional governing equations, we introduce the following non-dimensional variables,
\begin{equation} \frac {\partial }{\partial t}=\frac {U}{D}\frac {\partial }{\partial t^*},\quad \boldsymbol{u}=U\boldsymbol{u}^*,\quad p=\rho U^2p^*,\quad \boldsymbol{\nabla }=\frac {1}{D}\nabla^{*}. \end{equation}
Then, we obtain the non-dimensional governing equations:
\begin{align}\nabla^{*}\boldsymbol{\cdot }\boldsymbol{u}^*&=0,\end{align}
\begin{align}\frac {\partial \boldsymbol{u}^*}{\partial t^*}+\big(\boldsymbol{u}^*-\boldsymbol{e}_x{V}_x^*\big)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^*&=-\nabla^{*}p^*+\frac {1}{\textit{Re}}{\nabla }^{*2}\boldsymbol{u}^*+\frac {1}{{Re} Ca}\nabla^{*}\boldsymbol{\cdot }\left (\alpha \boldsymbol{B}^{\prime*}\right )-\boldsymbol{e}_x\frac {\mathrm{d}P^*}{\mathrm{d}x^*}_{\textit{drive}},\end{align}
\begin{align}\frac {\partial \alpha }{\partial t^*}+(\boldsymbol{u}^*-\boldsymbol{e}_x{V}_x^*)\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha &=0,\end{align}
\begin{align}\frac {\partial \boldsymbol{B}^*}{\partial t^*}+(\boldsymbol{u}^*-\boldsymbol{e}_x{V}_x^*)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B}^*&=\boldsymbol{L}^*\boldsymbol{\cdot }\boldsymbol{B}^*+\boldsymbol{B}^*\boldsymbol{\cdot }\boldsymbol{L}^{*\top },\end{align}
where
$V_x$
denotes the
$x$
-component of the particle’s velocity,
$\alpha$
denotes the volume fraction of solid,
$\boldsymbol{L}=\partial\boldsymbol{u}/\partial\boldsymbol{x}$
denotes the velocity gradient tensor and
$\boldsymbol{B}=\boldsymbol{F\boldsymbol{\cdot }F}^{\top }$
denotes the left Cauchy–Green deformation tensor, where
$\boldsymbol{F}=\partial \boldsymbol{x}/\partial \boldsymbol{X}$
is the deformation gradient,
$\boldsymbol{x}$
is the current coordinates and
$\boldsymbol{X}$
is the reference coordinates (Bonet & Wood Reference Bonet and Wood2008).
Following Sugiyama et al. (Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011), we monolithically describe the pressure over the entire incompressible media. Further, we write the mixture stress as
$\boldsymbol{\sigma }=-p\boldsymbol{I}+(1-\alpha ){}\boldsymbol{\sigma}^{\prime}_{f}+\alpha \boldsymbol{\sigma}^{\prime}_{s}$
, where
$\boldsymbol{\sigma }^{\prime}$
denotes the deviatoric stress, and the subscript
$f$
and
$s$
refer to the fluid and solid phases, respectively. The deviatoric tensor is defined as
$\boldsymbol{T}^{\prime}=\boldsymbol{T}-\mathrm{tr} (\boldsymbol{T} )\boldsymbol{I}/3$
. The term
$\mathrm{d}P^*/\mathrm{d}x^*_{\textit{drive}}$
indicates the constant driving pressure, which is analytically derived from the solution of duct flow without particles (Cornish Reference Cornish1928):
\begin{equation} -\frac {\mathrm{d}P^*}{\mathrm{d}x^*}_{\textit{drive}}=\frac {4}{\textit{K{Re}}},\quad K\approx {0.160577}\ldots . \end{equation}
Although the fully Eulerian method (Sugiyama et al. Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011) has been demonstrated to facilitate solving various problems such as flows bounded by compliant walls (Rosti & Brandt Reference Rosti and Brandt2017; Esteghamatian, Katz & Zaki Reference Esteghamatian, Katz and Zaki2022) and flows including soft particles (Rosti & Brandt Reference Rosti and Brandt2018; Rosti et al. Reference Rosti, Brandt and Mitra2018; Prasad, Sharma & Kulkarni Reference Prasad, Sharma and Kulkarni2022), the numerical solution is likely to be unstable especially when long-time computation is required, as in the present study. Such a numerical instability is more pronounced with the higher advection speed of the fluid–structure interface.
In this study, we investigate the motion of particles moving downstream in the
$x$
-direction. The particle’s
$x$
-coordinate can be fixed within the computational domain, as the computational system is homogeneous in the
$x$
-direction. This technique is implied by the term including
$\boldsymbol{e}_xV_x^*$
in (2.3)–(2.5), minimising the advection velocity. Note that this formulation is based on a coordinate transformation, in which the
$x$
coordinate in the computational space attaches to the particle centroid, and holds even if
$V_x^*$
varies with time. For the time and spatial discretisation of the advection terms with the constant velocity (
$\partial _{t^*}-\boldsymbol{e}_xV_x^*\boldsymbol{\cdot }\boldsymbol{\nabla }$
) in (2.3)–(2.5), a discrete Fourier interpolation (Gerz, Schumann & Elghobashi Reference Gerz, Schumann and Elghobashi1989) is applied to reduce numerical diffusion and instability.
Furthermore, the MTHINC scheme (Ii et al. Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012b
) is employed for the convection of the volume fraction (
$\partial _{t^*}\alpha +({\boldsymbol{u}}^*\boldsymbol{\cdot }\boldsymbol{\nabla })\alpha$
) in (2.4). The MTHINC scheme approximates the VOF function distribution across the elastic-fluid interface using a hyperbolic tangent function, effectively suppressing interface smearing and numerical instability. The sharpness parameter relative to the grid is set to
$\beta =2$
, following the rigid particle simulation in Appendix A.
Other spatial derivative terms are discretised using the second-order central scheme except for the advection term of
$\boldsymbol{B}$
. The advection term
${\boldsymbol{u}}^*\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B}^*$
in (2.5) is spatially discretised by the fifth-order WENO scheme (Liu, Osher & Chan Reference Liu, Osher and Chan1994; Jiang & Shu Reference Jiang and Shu1996).
For the time procedures, the second-order Adams–Bashforth scheme (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2012) is applied to the convection term of
$\boldsymbol{u}^*\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^*$
in (2.3),
$\boldsymbol{u}^*\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B}^*$
and the right-hand side of (2.5). The Crank–Nicolson scheme (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2012) is applied to the viscous and elastic terms in (2.3). To implicitly treat the pressure together with (2.2), the Simplified Marker-and-Cell (SMAC) method (Amsden & Harlow Reference Amsden and Harlow1970) is employed. To solve the Poisson equation directly, we use a fast Fourier transform (FFT) and a tri-diagonal matrix algorithm (TDMA).
For the boundary conditions, we apply periodic conditions along the channel and no-slip conditions for the walls. The computational domain is discretised with a grid size of
$256\times 128\times 128$
and
$128\times 64\times 64$
for
$\kappa =0.1$
and
$\kappa =0.2$
, respectively. The channel length and width are defined as
$L=2$
and
$D=1$
, respectively, and the grid resolution is given by
$\Delta x=\Delta y=\Delta z=1/128\text{ or }1/64$
. The effect of channel length
$L$
associated with the periodic boundary condition is discussed and the present choice of
$L/D=2$
is validated in figure S1 of the Supplementary material, which is available at https://doi.org/10.1017/jfm.2025.10574.
The Courant–Friedrichs–Lewy (CFL) number is fixed at a constant value of 0.0625, and the time step
$\Delta t$
is determined based on the maximum propagation velocity
$U_{\textit{prop}}$
. Here,
$U_{\textit{prop}}$
is determined by comparing the maximum advection velocity
$U_{\textit{adv}}= \max (|u_x-V_x|,|u_y|,|u_z|)$
and the shear/transverse wave speed in the elastic particle
$U_{\textit{els}}=\sqrt {G/\rho }$
.
\begin{align} U_{\textit{prop}}&=\max \left (U_{\textit{adv}},U_{\textit{els}}\right ), \end{align}
\begin{align} \Delta t&=\frac {\mathrm{CFL}\boldsymbol{\cdot }\Delta x}{U_{\textit{prop}}}. \end{align}
In particular, for the advection terms, we employ the advective interpolation scheme (Kajishima & Taira Reference Kajishima and Taira2016, § 3.5.1), which satisfies the relations
$\partial _j(u_ju_i)=u_i\partial _ju_j+u_j\partial _ju_i$
in a discretised form and ensures that the momentum and the quadratic quantity are highly conserved.
We also employ the immersed boundary method (Kajishima et al. Reference Kajishima, Takiguchi, Hamasaki and Miyake2001) for a calculation of a rigid particle, as detailed in Appendix A, and compare the results with those of the deformable particles.
2.4. Taylor deformation parameter
The Taylor deformation parameter
$\mathcal{T}$
is introduced as a measure of the deformation of elastic particles. This parameter is mathematically defined by the lengths of the semi-major axes
$l_m$
and semi-minor axes
$l_n$
of an ellipse in a two-dimensional plane (Taylor Reference Taylor1934),
\begin{equation} \mathcal{T}=\frac {l_m-l_n}{l_m+l_n}. \end{equation}
Here,
$\mathcal{T}$
decreases as the shape of an object approaches a circle, whereas greater deformations result in higher
$\mathcal{T}$
values. As this study involves three-dimensional particles, it is necessary to determine quantities corresponding to the semi-major and semi-minor axes in a spherical geometry. In this study, we evaluate
$\mathcal{T}$
using another well-known approach by approximating the shape of a deformed sphere as a triaxial ellipsoid with the same moment of inertia to capture the deformation in three dimensions (Ramanujan & Pozrikidis Reference Ramanujan and Pozrikidis1998; Takeishi et al. Reference Takeishi, Rosti, Imai, Wada and Brandt2019). The principal moments of inertia
$I_{a,b,c}$
relative to the centre of mass of the triaxial ellipsoid can be determined using the VOF function
$\alpha$
. The semi-axes of the triaxial ellipsoid
$a,b,c$
are obtained from the following equations:
\begin{equation} a=\sqrt {\frac {5}{2M}(-I_a+I_b+I_c)},\quad b=\sqrt {\frac {5}{2M}(I_a-I_b+I_c)},\quad c=\sqrt {\frac {5}{2M}(I_a+I_b-I_c)}. \end{equation}
where
$M$
denotes the mass of the particle. Finally,
$l_m$
and
$l_n$
are determined as
$l_m=\max (a,b,c), l_n=\min (a,b,c)$
.
3. Results and discussion
3.1. Experimental results
Figures 4(a)–4(e) show representative distributions of hydrogel particles in the channel cross-section for
$L/D=1000{-}12\,000$
at
${Re}=0.1{-}10$
, and figure 4(f) shows probability density functions (p.d.f.s) of particles for
$L/D=12\,000$
. At
${Re} = 0.1$
shown in figure 4(a), particles gradually approach the centre of the channel cross-section, and focus near the centre at the most downstream (see figure 4
f). This axial accumulation is consistent with previous studies for deformable particles at low
${Re}$
(Karnis et al. Reference Karnis, Goldsmith and Mason1963; Goldsmith Reference Goldsmith1971). At
${Re}=0.5$
and 1, particles migrate inwards until they focus along a small ring together with the channel centre. At
${Re} = 5$
and 10, particles are aligned along a ring (pSS ring) upstream and focus near four points on the diagonals downstream. This focusing position corresponds to the DEP. It is noteworthy that the particle focusing positions move away from the channel centre with increasing
${Re}$
.
Figure 4. (a)–(e) Distributions of hydrogel particles
$(d= 5.2\,\unicode{x03BC} \textrm {m},D=50\,\unicode{x03BC} \textrm {m})$
and (f) probability density functions (p.d.f.s) at
$L/D=12\,000\, (L = 600\,\textrm {mm})$
.
Figure 5. (a) Average distances
$\langle r\rangle /(D/2)$
of hydrogel particles from the channel centre and (b) number fraction Pd of particles located within
$\pm 10^\circ$
from the diagonals, at
${Re} = 0.1$
(triangles), 0.5 (open squares), 1 (closed squares), 5 (open circles) and 10 (closed circles).
For the particle distributions shown in figures 4(a)–4(e), figures 5(a) and 5(b) plot the average distance of the particle centroid from the channel centre,
$\langle r\rangle /(D/2)$
, and the number fraction of particles located within
$\pm 10^\circ$
from the diagonals,
$P_d$
(see the inset of figure 5
b), as a function of
$L/D$
. The error bar represents the standard deviation (only the upper or lower side is plotted). At low
${Re}\,(\leqslant 1)$
, the average distances first decrease with
$L/D$
up to 6000 and become almost constant further downstream, while
$P_d$
keeps nearly constant values
${\sim} 20^\circ /90^\circ =0.22$
, that is, for random distribution, for all
$L/D$
. This result indicates that, at
${Re}\leqslant 1$
, particles have reached equilibrium positions in the radial direction within a distance of
$L/D={6000}$
from the channel inlet, with no preferential distribution in the azimuthal direction. In contrast, the average distances at
${Re}\geqslant 5$
remain nearly constant independent of
$L/D\ (\geqslant {1000})$
and
$P_d$
increases almost monotonically with
$L/D$
. This result implies that, at
${Re}\geqslant 5$
, particles have reached an equilibrium radial position (pSS ring) in a short distance from the channel inlet
$(L/D \lt {1000})$
and migrate circumferentially towards the diagonal further downstream. The first part of this migration, mostly in the radial direction towards the pSS ring, is known as the first phase of the inertial migration, which represents a rather rapid motion, and the following slow migration along the pSS ring towards the DEP is the second phase. The present result demonstrates that hydrogel particles also exhibit a two-phase property of the inertial migration in square channel flows.
Figure 6. Distributions of (a) hydrogel particles in glycerol aqueous solutions, (b) rigid spherical particles
$(d=5\,\unicode{x03BC} \textrm {m})$
in glycerol aqueous solutions and (c) human red blood cells in blood plasma (Tanaka & Sugihara-Seki Reference Tanaka and Sugihara-Seki2022)
$(D=50\,\unicode{x03BC} \textrm {m}, L=600\,\textrm {mm})$
.
Figure 6 illustrates a comparison of the distribution of hydrogel particles with those of rigid spherical particles
$(d=5\,\unicode{x03BC} \textrm {m})$
and human red blood cells (RBCs) reported in our previous study (Tanaka & Sugihara-Seki Reference Tanaka and Sugihara-Seki2022). RBCs have a biconcave discoid shape at rest, with an equivalent diameter
$d\sim 5.6\,\unicode{x03BC} \textrm {m}$
, and they are highly deformable, as noted in § 1. In the case of small inertia
$({Re}=0.1)$
, deformable particles (hydrogel particles and RBCs) focus near the channel centre, whereas rigid spherical particles are dispersed widely in the channel cross-section, as inferred from the motion of rigid spherical particles along the channel axis in the Stokes flow. On the other hand, in the case of finite inertia at
${Re}=10$
, both hydrogel particles and RBCs focus near the DEP, whereas rigid spherical particles focus near the MEP. This result indicates that the focusing on the DEP results from the effect of particle deformation. At
${Re}=1$
, rigid particles focus along the pSS ring with higher concentrations near the MEP. With regards to deformable particles, hydrogel particles are aligned along a small ring together with the channel centre, whereas RBCs focus on four points on the diagonals together with the centre. Although these two focusing patterns of deformable particles may appear different, they are essentially similar and the difference presumably arises from differences in the size of the pSS ring. Possibly due to differences in the deformability and shape of the suspended particles and/or the rheological properties of the suspending media, the pSS ring for the hydrogel particles is much smaller than that for RBCs. Since the flow field near the channel centre is nearly axisymmetric, the azimuthal component of the lift exerted on the hydrogel particles on the pSS ring is much weaker than that on the RBCs. As a result, the focusing in the second phase of the inertial migration is much slower for the hydrogel particles than that for RBCs.
3.2. Numerical results
Figure 7. Comparisons between the experimental and numerical results for the particle distribution in the channel cross-section at corresponding distances from the inlet
$(\kappa = 0.1)$
. The upper panels: trajectories of the centroid of hyperelastic particles from initial positions (open circles) to final positions (closed circles) during the dimensionless time (a)–(c)
$t^*(\equiv tU/D)\sim {2000}$
and (d,e) 1000, in the first quadrant of the cross-section. The lower panels: final positions in the entire cross-section (dots) and corresponding particle distributions obtained experimentally at (a)–(c)
$L/D = {2000}$
, (d,e) 1000.
The elastic modulus
$G$
or the Laplace number
$La$
for hyperelastic particles was searched for the one that best reproduces the distribution of hydrogel particles experimentally observed, since the elastic moduli of these particles are unknown. First, numerical simulations with
$\kappa =0.1$
and
${Re}=10$
corresponding to figure 4(e) were performed for various
$La$
values. With
$La$
set to 5, 10 and 50, trajectories of the centroid of hyperelastic particles projected onto the channel cross-section are shown in the upper panels of figure 7(a–c), starting from various initial positions (open circles) located in the lower right half of the first quadrant to final positions (closed circles) during a dimensionless time
$t^*(\equiv tU/D)\sim {2000}$
. The lower panels illustrate their final positions in the entire cross-section, obtained by taking into account the symmetry with respect to the
$y$
-,
$z$
-axes (midlines) and diagonals, along with the distribution of hydrogel particles observed experimentally at
$L/D = {2000}$
(figure 4
e). Comparisons between experimental and numerical results in figure 7(a–c) indicate that the experimentally obtained distribution is best reproduced when adopting
$La = 10$
. To confirm this value, the numerical simulations with
$La = 10$
for different
${Re}$
and different distances
$L/D$
are performed and compared with corresponding distributions of hydrogel particles. Figures 7(d) and 7(e) are two examples, in which the experimental and numerical results are in good agreement in favour of
$La\sim 10$
.
The Laplace number
$La = 10$
gives the elastic modulus
$G=1.04\,\textrm {kPa}$
, using the values for the density and viscosity of the fluid and the particle diameter. Assuming a Poisson’s ratio of 0.5, Young’s modulus of the hydrogel particles is estimated to be
${\sim}3.13\,\textrm {kPa}$
. Banquy et al. (Reference Banquy, Suarez, Argaw, Rabanel, Grutter, Bouchard, Hildgen and Giasson2009) used atomic force microscopy to measure Young’s modulus of similar but much smaller hydrogel particles and reported that it ranges from 18 kPa for particles with a cross-linker content of 1.7 mol % to 211 kPa for particles with a cross-linker content of 15 mol %. The present estimate of Young’s modulus
${\sim}3.1\,\textrm {kPa}$
for particles with a cross-linker content of 1 mol % is slightly smaller but comparable to extrapolations from these reported values. Although the distributions of hydrogel particles experimentally obtained at
${Re}\geqslant 3$
are successfully reproduced by the present numerical simulation, comparisons at lower
${Re}(\lesssim 1)$
are difficult because of the large computational time required for small particles with
$\kappa =0.1$
to reach specific focusing positions at lower
${Re}$
.
Figure 8. (a) Trajectories of the centroid of hyperelastic particles in the first quadrant of the channel cross-section, (b) the final positions in the entire cross-section and (c) snapshots of particles at the final position, for
$\kappa =0.2$
at
${Re}=40$
.
To explore the effect of particle deformation on the inertial migration in wider ranges of
${Re}$
and
$La$
, we computed the motion and deformation of hyperelastic particles with larger blockage ratios. Figure 8 shows trajectories of the particle centroid starting from two initial positions
$(y_0^*,z_0^*)=(0.35,0.07)$
and
$(0.35,0.28)$
, their final positions in the entire cross-section obtained from symmetry considerations and final particle shapes, for
$\kappa =0.2$
at
${Re}=40$
in the range of
$La$
from 10 to
$\infty$
(rigid spherical particles). The computations for rigid spherical particles were conducted using a computer scheme reported previously (Nakagawa et al. Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015). In all cases, particles first migrate inwards until they reach a certain distance from the channel centre and then move circumferentially towards the
$y$
-axis or diagonal, except in the case of
$La=10$
, where they continue to move towards the channel centre. This feature clearly shows the two-phase property of the inertial migration in square channel flows, and the circumferential motion in the second phase is performed along the pSS ring towards a stable equilibrium position. The final position of rigid spherical particles is located on the
$y$
-axis (MEP), in agreement with previous studies (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Miura et al. Reference Miura, Itano and Sugihara-Seki2014). The hyperelastic particles with relatively large
$La$
(
$=400$
and 200) also exhibit focusing on the MEP, and the focusing positions move closer to the channel centre with decreasing
$La$
. On the other hand, for smaller
$La(=50)$
, the final focusing positions are located on the diagonals (DEP). In the case of
$La=120$
between these two cases, hyperelastic particles focus on eight points on the midlines and diagonals in the entire cross-section, implying a bistable state of the MEP and DEP. At the smallest value of
$La(=10)$
, a stable equilibrium position is located on the channel centre even at moderate
${Re}(=40)$
.
Figure 8(c) depicts snapshots of hyperelastic particles located at the final position, when viewed along the
$x$
-,
$y$
- and
$z$
-axes. Note that these are shapes at a certain instant in time. It is seen from figure 8(c) that, in general, a smaller
$La$
results in a larger particle deformation, as expected. However, the shift of the final position towards the channel centre with decreasing
$La$
suppresses that increase (see figures 10
b and 10
c). Notably, for
$La=10$
, the deformation is rather small and the particle has an almost axisymmetric shape, rounded at the front and nearly flat at the back, since its final position is very close to the channel centre, where the shear rate is small and axisymmetric.
Figure 9. (a) Trajectories (above the diagonal) are shown with dashed lines for
$t^*\leqslant 20$
and solid lines for
$t^*\gt 20$
, while the magnitude of the in-plane velocity (below the diagonal) is represented by the colour bar. (b) Time evolution of the dimensionless
$y$
-coordinate
$y^*$
(top) and
$z$
-coordinate
$z^*$
(bottom), with the inset focusing on the
$t^*$
-axis ranges from 0 to 100. The initial positions are
$(y_0^*,z_0^*)=(0.35,0.07),(0.35,0.28)$
. The blockage ratio
$\kappa =0.2$
and
${Re}=40$
.
To demonstrate the two-phase property of inertial migration in more detail, figure 9 shows the lateral speed and the movement of the particle centroid. The particle is seeded at the same initial positions as in figure 8 or symmetric positions relative to the diagonal for
$La=50\!-\!400$
. In the upper left half of figure 9(a), the trajectories of the particle are drawn with dashed lines for
$t^*\leqslant 20$
, corresponding to a roughly estimated dimensionless time for the first phase, and the trajectories for
$t^*\gt 20$
are drawn with solid lines. In the lower right half of figure 9(a), the lateral velocity magnitude of the particle is indicated by colour along the trajectories. The variations of the position of the particle centroid are plotted as a function of time in figure 9(b). It is seen from figure 9(a) that the lateral velocity in the first phase is
$O(10^{-2})$
relative to the mean flow velocity, whereas that in the second phase is
$O(10^{-3})$
for
$La=50$
and
$O(10^{-4})$
for
$La=120{-}400$
. The rapid migration in the first phase and the slow migration in the second phase can be also confirmed from figure 9(b) together with the insets. A distinct feature observed in figure 9(b) is that the centroid of hyperelastic particles fluctuates even after approaching a stable position.
In each case shown in figures 8 and 9, the focusing pattern of particles can be easily identified. However, for some parameter values including
$La$
between 10 and 50 at
${Re} = 40$
, the process became highly challenging. In some cases, the migration velocity of the particle becomes very small before reaching the midline, diagonal or centre. In other cases, excessive deformation caused numerical instability. The former case suggests the presence of another type of equilibrium position, presumably located between the MEP and DEP on the pSS ring. This focusing position may correspond to the equilibrium position called the intermediate equilibrium position (IEP) by Esposito et al. (Reference Esposito, Romano, Hulsen, D’Avino and Villone2022). Their numerical analyses showed that hyperelastic particles suspended in square channel flows migrate to the IEP between the midline and the diagonal (eight equilibrium positions in the whole cross-section) in a rather wide range of parameters. In the present study, it is practically difficult to determine whether the particle’s intermediate position is a focusing position or a transient one. This is mainly because the lateral migration in the second phase is generally slow as seen in figure 9, and it becomes extremely slow near the transitions between different focusing positions.
To avoid ambiguity, we performed the computation in fixed time steps of
$4\times 10^6$
at a constant
$\mathrm{CFL}=0.0625$
, starting from several initial positions of the particle for various
${Re}$
and
$Ca$
at
$\kappa =0.2$
. Each final position is classified into four types as follows: CEP when the final position is located within
$0.05D$
from the channel centre, MEP or DEP when the azimuthal angle of the final position from the midline or diagonal is within
$\pm \pi /30$
, and IEP for the rest. When the final position is categorised as the IEP, the extra computation of
$3\times 10^6$
steps is conducted and its final position is checked again according to the above criteria. This procedure is repeated as far as the final position is categorised as the IEP, up to the total time steps of
$2.2\times 10^7$
. The computational time
$t^*=tU/D$
for
$2.2\times 10^7$
time steps depends on
$U_{\textit{prop}}$
(which is primarily determined by
${Re}$
and
$Ca$
) under the constant CFL number, with the maximum
$t^*\sim 1.029\times 10^4$
for
${Re}=50$
and
$Ca=0.02$
, and the minimum
$t^*\sim 1.359\times 10^3$
for
${Re}=1$
and
$Ca=0.004$
.
Figure 10. (a) Particle focusing patterns, (b) distance of the final position from the channel centre and (c) time average of
$\mathcal{T}$
at the final position for
$\kappa =0.2$
. Each dataset is calculated for each
${Re}$
and
$Ca$
from three initial positions:
$(y_0^*,z_0^*)=(0.0625,0.03125),(0.3125,0.03125),(0.3125,0.28125)$
. Each symbol represents a focusing position: circles for CEP, squares for MEP, diamonds for DEP and triangles for IEP. Each coloured region corresponds to a focusing pattern: (A) CEP, (B) DEP, (C) MEP and (D) MEP + DEP.
The results of such a classification with three initial positions
$(y_0^*,z_0^*)=(0.0625,{}0.03125)$
,
$(0.3125,0.03125)$
and
$(0.3125,0.28125)$
are shown in figure 10(a), in which three symbols representing the type of the final positions are superimposed at each point in the parameter space for
$1\leqslant {Re}\leqslant 10^2$
and
$4\times 10^{-3}\leqslant Ca\leqslant 3\times 10^{-2}$
. When three symbols are identical (in most of these cases, the three final positions are very close to each other), it is highly suggested that the final position represents a stable equilibrium position or a focusing position observed in the experiment. Such cases consist of a single-point focusing at the CEP (pattern A, area coloured yellow in figure 10
a), four-point focusing at the DEP (pattern B, area coloured green) and four-point focusing at the MEP (pattern C, area coloured blue). In the parameter space between patterns (B) and (C) for
$Ca\geqslant 10^{-2}$
highlighted in orange, a bistable state of the DEP and MEP is suggested. This focusing pattern, corresponding to the pattern for
$La=120$
in figure 8(b), is referred to as pattern (D). In uncoloured areas near the edge of coloured areas, where transitions between different focusing patterns are expected, the IEP (triangles) is observed together with other symbols. This indicates that all three final positions do not converge to a single point in the current computational time. Although single or multiple stable states of equilibrium positions (CEP, DEP, MEP and IEP) may be reached after much longer computations, it is currently hard to identify the focusing pattern in these areas. As far as examined in the present study, we do not observe a sole focusing pattern at the IEP. This is a significant difference from the results of Esposito et al. (Reference Esposito, Romano, Hulsen, D’Avino and Villone2022), as discussed later.
The effect of particle deformability is dominant and that of inertia is small in the upper left of figure 10(a), and vice versa in the lower right. From top left to bottom right, coloured areas sequentially represent patterns (A), (B) and (C). There is an area of pattern (D) between patterns (B) and (C). A notable feature in figure 10(a) is that the boundary between patterns (A) and (B) lies approximately on a straight line with
$La=8$
, which is related to the stability of the CEP. Similarly, the stability of the DEP may be changed on another line with
$La\sim 200$
. The implication of the former is discussed in § 3.3. The latter is beyond the scope of the present study and will be left for future work.
We consider the transitions of the particle focusing pattern shown in figure 10(a), based on the effects of inertia and particle deformation. As noted in § 1, the inertial lift consists mainly of the shear gradient-induced lift, acting in the direction towards higher shear, and the wall effect, acting in the direction away from the channel wall. In the square cross-section, the direction of the shear gradient-induced lift is mostly outwards in the radial direction and towards the
$y$
- or
$z$
-axis (midlines) in the azimuthal direction. Thus, rigid particles are focused on the midlines (MEP) at finite
${Re}$
and small
$Ca$
(pattern C). On the other hand, if the particle deformation-induced lift is exerted on hyperelastic particles in the direction towards lower shear, as noted already, then its direction is nearly opposite to the shear gradient-induced lift, i.e. mostly inwards and towards the diagonal in the square cross-section. As a result, more deformable particles tend to focus on the diagonal (DEP) rather than on the midline (MEP), and their focusing positions approach the channel centre (CEP), as seen in figure 10(a). Thus, an increase in the deformation-induced lift directed towards lower shear accounts for the transition of the focusing pattern from pattern (C) through (B) to (A) shown in figure 10
(a).
Figures 10(b) and 10(c) depict the distance of the final position from the channel centre and the average of the Taylor deformation parameter
$\mathcal{T}$
during the last 50 000 steps, respectively. Figure 10(b) clearly shows that, at constant
${Re}$
, all of MEP, DEP and IEP approach the channel centre with decreasing
$La$
. This indicates that the final positions of more deformable particles approach closer to the channel centre. All curves for the DEP roughly collapse onto a single curve, which is similar to the feature observed for capsules (Schaaf & Stark Reference Schaaf and Stark2017). Figure 10(b) shows that this is also almost true for the MEP, and that the master curve for the MEP is located below the DEP master curve. Additionally, the DEP appears only for
$La\gtrsim 8$
, below which the focusing position is changed to the channel centre (CEP), as expected from figure 10(a).
Figure 10(c) shows that
$\mathcal{T}$
at the final position increases with increasing
$La$
up to several tens, and then gradually decreases at each
${Re}$
. Considering that
$La$
provides a measure of the impact of inertia relative to particle deformation, a decrease in
$\mathcal{T}$
with increasing
$La$
seems reasonable. The inverse trend for small
$La$
is explained by the fact that more deformable particles tend to approach closer to the channel centre, where smaller shear rates cause smaller deformation.
Recently, Esposito et al. (Reference Esposito, Romano, Hulsen, D’Avino and Villone2022) reported focusing patterns of hyperelastic particles for
$\kappa =0.2$
in wider ranges of
${Re}$
and
$Ca$
. They found four types of particle focusing patterns, among which patterns (A) and (C) are consistent with the present study. Each of patterns (A) and (C) appears in a similar range of
${Re}$
and
$Ca$
in figure 10(a), whereas in their study, the focusing on the IEP emerges in the range where patterns (B) and (D) appear in the present study, and pattern (B) does not appear in the present range of
$1\leqslant {Re}\leqslant 10^2$
and
$4\times 10^{-3}\leqslant Ca\leqslant 3\times 10^{-2}$
. Esposito et al. (Reference Esposito, Romano, Hulsen, D’Avino and Villone2022) never found pattern (D) in their whole range of parameters. These discrepancies may be due to the different models used for the particles: the particles in the present study are viscoelastic, but those in their models include no viscosity. Differences in the criteria used to categorise the focusing positions may also be possible causes of the discrepancies.
Figure 11. Distance of the final position from the channel centre. Coloured plots show results for particles in the present study, and black markers represent capsule data from Krüger et al. (Reference Krüger, Kaoui and Harting2014b
)
$(\kappa =0.2)$
.
As deformable particles, capsule models have been often adopted in numerical simulations to study inertial migration. Among these studies, Krüger et al. (Reference Krüger, Kaoui and Harting2014b
) investigated the rheological properties of a capsule suspension in the pressure-driven flow between two infinite parallel plates (Misbah Reference Misbah2014). At a low volume fraction of particles, they demonstrated that the lateral equilibrium position of the capsules increases with increasing
${Re}$
up to
$45$
(corresponding to
${\sim} 60$
in the present definition of Reynolds number), followed by a decrease, as shown in figure 11. Figure 11 replots the present results of the final lateral position
$r^*_e$
shown in figure 10(b) as a function of
${Re}$
for various
$Ca$
, superimposing the results of Krüger et al. (Reference Krüger, Kaoui and Harting2014b
). In their study, the capillary number
$Ca'$
is defined as the ratio of a typical shear stress magnitude to a characteristic elastic particle stress
$\kappa _s/a$
, where
$\kappa _s$
and
$a$
represent the shear elasticity of the capsule membrane and the radius of the undeformed capsule, respectively. Figure 11 demonstrates that each plot at a constant
$Ca$
or
$Ca'$
exhibits a similar upward convex property, although the hyperelastic particles show a gradual decrease in
$r_e^*$
at
${Re} \gt 40$
only for cases with larger
$Ca (\geqslant 0.02)$
. For the hyperelastic particles with smaller
$Ca$
,
$r^*_e$
drops sharply at
${Re}\sim 40$
due to the transition of stable equilibrium position from the DEP to the MEP (see figure 10
a). Additionally, both particles exhibit the trend that
$r^*_e$
at constant
${Re}$
is smaller for larger
$Ca$
, implying that more deformable particles tend to focus closer to the channel centre.
For capsules in a confined square channel, Schaaf & Stark (Reference Schaaf and Stark2017) reported capsule migration behaviour in the same geometry. Although strict quantitative comparisons are not possible due to differences in the deformation metrics used to characterise particle stiffness (in terms of
$La$
and
$Ca$
), their results show behaviour broadly similar to our findings for viscous hyperelastic particles. As shown in figures 3(b) and 3(d) of Schaaf & Stark (Reference Schaaf and Stark2017), which correspond to our results at
${Re}\sim 5$
and
$50$
, respectively, the changes in focusing patterns exhibit nearly the same trends. Capsules migrate towards the CEP at low
$La$
, and the focusing positions transition from the DEP to the MEP with increasing
$La$
. Regarding the migration velocity, a two-phase migration property is also observed, with softer capsules migrating faster, consistent with the elastic particle results in figure 9. The primary difference lies in the bistable state between the DEP and MEP. This bistability is likely absent in the results of Schaaf & Stark (Reference Schaaf and Stark2017) because they computed trajectories from a single initial position, whereas capturing bistability requires simulations from at least two distinct initial positions. Another difference is the dependence of
$r^*_e$
and
$\mathcal{T}$
on
$La$
, shown in figure 10(b, c). Schaaf & Stark (Reference Schaaf and Stark2017) reported a critical Laplace number
$La_c \sim 1$
for capsules with
$\kappa =0.2$
to move away from the CEP and noted that
$La_c$
changes sensitively with
${Re}$
. In contrast, our results indicate
$La_c \sim 8$
without clear dependence on
${Re}$
within the low-
${Re}$
regime as shown in figure 10(b). Additional details are provided in figures S2 and S3 in the Supplementary material.
In some cases in the present study, identifying the equilibrium position of the particle is very difficult, mainly due to their extremely slow migration in the second phase. This difficulty is unavoidable when the equilibrium position is determined from the calculation of trajectories. Instead of calculating trajectories, Schaaf & Stark (Reference Schaaf and Stark2017) estimated the lift exerted on a capsule by applying an adjustable force evenly distributed over all the membrane vertices to hold the capsule in place. By evaluating the external force to compensate the lift for various prescribed positions of the capsule, they obtained the lift map and determined the equilibrium position where the lift vanishes. In the present study, however, it is difficult to apply this method to the hyperelastic particles and thus to obtain the lift map.
3.3. Discussion of the stability of the channel centre
The phase diagram (figure 10
a) suggests that the boundaries of focusing patterns appear to be formed along lines of constant
$La$
. Here,
$La$
is defined by the ratio of
${\textit{Re}}_p$
to
$Ca$
, indicating that both inertia
$({\textit{Re}}_p)$
and deformation
$(Ca)$
influence the transition of the focusing positions. This trend becomes even more apparent in the phase diagram plotted using
$({Re}, La)$
coordinates, as shown in figure S4 of the Supplementary material. In the following, we first explain why
$La$
serves as the characteristic scaling parameter for this phenomenon. Then, with the aid of additional numerical simulations, we present a semi-analytical discussion of the transition between CEP and DEP near the line of
$La=8$
using an asymptotic expansion of the forces acting on the particle, with
${\textit{Re}}_p$
and
$Ca$
as perturbed parameters.
3.3.1. Asymptotic expansion
Cox & Brenner (Reference Cox and Brenner1968) applied matched asymptotic expansions to predict the migration of rigid spheres in Poiseuille flow with low
${Re}$
. This approach was extended to include inertial effects on the lateral migration of a neutrally buoyant rigid sphere in a Newtonian fluid, through an asymptotic expansion in terms of
${Re}$
(Ho & Leal Reference Ho and Leal1974; Hood, Lee & Roper Reference Hood, Lee and Roper2015). For deformable drops, Chan & Leal (Reference Chan and Leal1979, Reference Chan and Leal1981) applied perturbation methods and the reciprocal theorem to analyse deformation-induced migration. The effects of inertia and deformation on the lateral migration of bubbles near a wall were experimentally investigated by Takemura et al. (Reference Takemura, Takagi, Magnaudet and Matsumoto2002), and theoretically analysed by Magnaudet, Takagi & Legendre (Reference Magnaudet, Takagi and Legendre2003) through a perturbation expansion in terms of the Galileo and Bond numbers. Sugiyama & Takemura (Reference Sugiyama and Takemura2010) also performed an asymptotic expansion in terms of
$Ca$
to analyse the lateral migration of slightly deformed bubbles rising near a wall, complementing the previous work by Magnaudet et al. (Reference Magnaudet, Takagi and Legendre2003).
Using such an asymptotic expansion, the force acting on an slightly deformable sphere can be expressed as
\begin{equation} \boldsymbol{F}=\mu{d}Uf(y,z)\boldsymbol{e}_x-\mu{d}\boldsymbol{D}\boldsymbol{\cdot }\boldsymbol{u}+{\textit{Re}}_p\mu{d}U\boldsymbol{f}_{\!I}+\textit{Ca}\mu{d}U\boldsymbol{f}_{\!D}+o\left ({\textit{Re}}_p\right )+o\left (\textit{Ca}\right ). \end{equation}
The first and second terms represent the force on a rigid sphere in a creeping flow. The third and fourth terms correspond to the first-order terms of inertial lift and deformation-induced lift, respectively, expanded in terms of
${\textit{Re}}_p$
and
$Ca$
. Here,
$f(y,z)$
is the force coefficient dependent on the cross-sectional position,
$\boldsymbol{D}$
is the drag coefficient matrix, and
$\boldsymbol{f}_{\!I}$
and
$\boldsymbol{f}_{\!D}$
are the dimensionless vectors representing the inertia and deformation-induced forces, respectively.
For small
${\textit{Re}}_p$
and
$Ca$
, the cross-sectional components (denoted by the
$\perp$
symbol) of the force can be expressed as
\begin{equation} \boldsymbol{F}^{\perp }=-\mu{d}\boldsymbol{D}^{\perp }\boldsymbol{\cdot }\boldsymbol{u}^{\perp }+{\textit{Re}}_p\mu{d}U\boldsymbol{f}_{\!I}^{\perp }+\textit{Ca}\mu{d}U\boldsymbol{f}_{\!D}^{\perp }. \end{equation}
The particle is at the equilibrium position, provided that
\begin{equation} \boldsymbol{F}^{\perp }=\boldsymbol{0},\quad \boldsymbol{u}^{\perp }=\boldsymbol{0}, \end{equation}
from which, (3.2) becomes
\begin{equation} La\boldsymbol{f}_{\!I}^{\perp }+\boldsymbol{f}_{\!D}^{\perp }=0, \end{equation}
where
$La={\textit{Re}}_p/Ca$
is the Laplace number. Since
$\boldsymbol{f}_{\!I}$
and
$\boldsymbol{f}_{\!D}$
are functions of the cross-sectional coordinates and the blockage ratio
$\kappa$
, (3.4) indicates that
$La$
determines the particle’s equilibrium positions. Once
$\boldsymbol{f}_{\!I}$
and
$\boldsymbol{f}_{\!D}$
are known, we may evaluate
$La$
at the equilibrium position. Subsequently, we shall demonstrate that (3.4) accounts for the change of the focusing pattern from the CEP to the DEP beyond
$La\sim 8$
in figure 10(a). Upon decomposing (3.2) on the polar coordinates with the origin at the channel centre, we consider the following relation between the radial force
$F_r$
and velocity
$u_r$
of the particle near the origin:
\begin{equation} F_r=-\mu{d}D_{rr}u_r+{\textit{Re}}_p\mu{d}Uf_{\textit{Ir}}+\textit{Ca}\mu{d}\textit{U f}_{\textit{Dr}}. \end{equation}
In § 3.3.2, we estimate
${f}_{\textit{Ir}}$
from numerical simulations for a rigid sphere
$(Ca=0)$
at finite
${Re}$
. In § 3.3.3, we estimate
${f}_{Dr}$
from simulations for a hyperelastic particle (finite
$Ca$
) at
${Re}=0$
. In § 3.3.4, we examine the stability of the CEP.
3.3.2. Rigid sphere fixed within the cross-section
Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) computed the inertial lift force by simulating a fixed rigid particle within the cross-section using COMSOL Multi-physics. They expressed it in a dimensional form as
$F_I=\rho U^2d^2\kappa f_I=\mu{d}U{\textit{Re}}_pf_I$
, where
$f_I$
is a function of the particle position. Nakagawa et al. (Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015) also investigated the lateral migration of rigid spheres using the immersed boundary method and clarified the stability of equilibrium positions through simulations of both fixed and moving particles. Following the approach of Nakagawa et al. (Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015), we employ the immersed boundary method (detailed in Appendix A) for the numerical simulation of a rigid sphere, which freely moves in the streamwise direction and rotates in all the directions but rests at the cross-sectional position.
For a rigid sphere fixed within the cross-section
$(u_r=0,f_{Dr}=0)$
, (3.5) becomes
\begin{equation} F_r={\textit{Re}}_p\mu{d}\textit{U f}_{\textit{Ir}}. \end{equation}
Computing
$F_r$
along the straight lines across the centre within the cross-section, we estimate radial profiles of
$f_{\textit{Ir}}$
based on (3.6), as shown in figure 12. Each coloured symbol represents the results along the midline
$(z_p/y_p=0)$
, the diagonal line
$(z_p/y_p=1.0)$
and the intermediate line
$(z_p/y_p=0.5)$
, where
$(y_p,z_p)$
is the particle’s centre of mass. The grey marks represent the results of Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) and Nakagawa et al. (Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015) at
$z_p/y_p=0$
, showing good agreements with the present results. We also simulated the moving particles and verified the particle motion towards the equilibrium point of
$f_{\textit{Ir}}=0$
. In all the cases, the particles near the channel centre migrated outwards, implying the instability of the CEP pattern as expected from figure 12. The solid lines represent odd polynomial fits (up to the fifth degree) of the discrete points near the origin
$(0\leqslant r^*\leqslant 0.25)$
. The dashed lines indicate linear fits with slopes
$a_1$
corresponding to the first-order coefficients of the fitted functions. The force approximately increases along these linear lines, which are nearly identical. Thus, we obtain
\begin{equation} f_{\textit{Ir}}=\frac {F_r}{\rho U^2d^2\kappa }\approx a_1r^*,\quad a_1\approx (1.238\pm 0.003)\ \text{at}\ r^*\ll 1. \end{equation}
Figure 12. Force profile acting on a particle fixed at
$(y_p,z_p)$
with
${Re} =40, Ca=0$
and
$\kappa =0.22$
, along three lines: (red squares) the midline
$z_p/y_p=0$
, (green triangles) the diagonal line
$z_p/y_p=1.0$
and (blue triangles) the intermediate line
$z_p/y_p=0.5$
. The case of
$z_p/y_p=0$
is compared with the previous studies: (grey circles) Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) and (grey diamonds) Nakagawa et al. (Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015). (Inset) Fitting functions of a fifth-degree odd polynomial and linear fits with slopes of the first-degree coefficients.
3.3.3. Deformable particle in Stokes flow
The deformation-induced force
$f_{Dr}$
in (3.5) is estimated by computing a hyperelastic particle moving in the Stokes flow. Assuming the quasi-steady motion of the particle (namely,
$F_r=0$
) and employing Stokes’s law
$(D_{rr}=3\pi )$
in an infinite domain with
${Re} =0$
, (3.5) becomes
\begin{equation} 0=-3\pi\!\mu{d}u_r+Ca\mu{d}\textit{U f}_{Dr}. \end{equation}
From (3.8),
$f_{Dr}$
is determined using the particle velocity
$u_r^*$
, namely
\begin{equation} f_{Dr}=\frac {3\pi }{\textit{CaU}}u_r=\frac {3\pi }{Ca}u_r^*\quad \left (\because \ u_r =Uu_r^*\right ). \end{equation}
Figure 13. (a) Trajectories of the centroid of hyperelastic particles for
${Re}=0, Ca=0.01$
and
$\kappa =0.2$
with initial positions
$y_0^*=0.375,z_0^*=0.03125,0.125,0.25,0.375$
. Circles indicate the initial positions, and crosses mark the end points of the calculation. (b) Velocity distributions of deformable particles in Stokes flow under the same conditions of panel (a). Dash-dotted lines indicate fifth-degree odd polynomial fits, while dotted lines represent linear approximations using the first-degree coefficients. The fitting range is
$r^*\in [0:0.3]$
.
The distribution of
$u_r^*$
is found as the numerical solution to the equation set (2.2)–(2.5), where the left-hand side of (2.3) is omitted. The numerical simulations were conducted for seven initial positions given by
$y_0^*=0.375,\ z_0^*=0.03125,\ 0.0625i\ (i=1,\ 2,\ldots ,\ 6)$
and three capillary numbers
$Ca=0.005,\ 0.01,\ 0.015$
. Figure 13 shows typical results of the particle trajectory and velocity. In all the cases, the particles migrate towards the channel centre, exhibiting the CEP pattern as expected (figure 13
a). It is seen from the inset of figure 13(b) that the plots for different initial positions are likely to collapse onto a single line near the origin. By fitting the plots of
$u_r^*$
versus
$r^*$
with a fifth-order polynomial (odd polynomial up to degree 5) and approximating it near the centre
$(r^*\ll 1)$
, the dimensionless velocity can be expressed as
\begin{equation} u_r^*\approx b_1^{'}r^*, \end{equation}
with
$b_1^{'}=1.07\times 10^{-2}$
for
$Ca=0.01$
. The asymptotic expansion in § 3.3.1 supposes that
$f_{Dr}$
in (3.9) is independent of
$Ca$
and dependent only on
$r^*$
near the origin for given
$\kappa$
. Introducing
$b_1=b_1'/Ca$
gives
\begin{equation} f_{Dr}=\frac {3\pi }{\textit{Ca}}b_1^{'}r^*=3\pi b_1r^*. \end{equation}
We find
$b_1=1.07$
for
$Ca=0.01$
,
$b_1=0.863$
for
$Ca=0.005$
, and
$b_1=1.12$
for
$Ca=0.015$
. Thus,
$b_1$
remains nearly constant, specifically
\begin{equation} b_1=-(1.02\pm 0.14). \end{equation}
3.3.4. Stability of CEP
Substituting (3.7) and (3.11) into (3.5) with
$u_r=0$
gives
\begin{equation} F_r\approx \textit{Ca}\mu{d} U\left (\textit{Laa}_1r^*+3\pi b_1r^*\right ). \end{equation}
Note that
$a_1$
has a positive value, indicating that the inertial force drives the particle away from the channel centre, while
$b_1$
has a negative value, resulting in the particle migrating towards the centre. Equations (3.3) and (3.13) indicate that the particle is at the equilibrium position if
$r^*=0$
or
$La a_1r^*+3\pi b_1r^*=0$
. Here, we consider a particle located slightly offset from the channel centre. In the case where
$La a_1r^*+3\pi b_1r^*\lt 0$
, the radial force
$F_r$
is negative and, therefore, the particle migrates towards the centre
$r^*=0$
, corresponding to the equilibrium position. By contrast, in the case where
$La a_1r^*+3\pi b_1r^*\gt 0$
(namely
$F_r\gt 0$
), the particle migrates away from the centre. Therefore, together with
$F_r=0$
, we find the critical Laplace number
$La_c$
:
\begin{align} \textit{La}_c&=-\frac {3\pi b_1}{a_1}=7.75\pm 1.05, \end{align}
beyond which the CEP pattern becomes unstable. This critical value in (3.14) shows a good agreement with the boundary at
$La\sim 8$
in figure 10(a), indicating that the asymptotic expansion explains the transition of the particle focusing between the CEP and DEP regions.
4. Conclusions
In the present study, we used hydrogel microspheres as deformable particles to investigate experimentally the inertial migration of deformable particles flowing through square channels. The experiments have shown that at finite
${Re}(\geqslant 5)$
, hydrogel particles with blockage ratio of 0.1 focus near four points located on the diagonals, which indicates the presence of stable equilibrium positions on the diagonal (DEP). As
${Re}$
decreases, the particle focusing position moves inwards until the particles focus on the channel centre (CEP), i.e. axial accumulation, at low
${Re}\,(=0.1)$
. Corresponding numerical simulations using viscous hyperelastic particles predicted that Young’s modulus of these hydrogel particles would be
${\sim}3.1\,\textrm {kPa}$
, slightly smaller but comparable to extrapolations from previously reported values. Additionally, two-phase migration property observed for hydrogel particles was reproduced by the numerical simulation for hyperelastic particles at finite
${Re}$
; particles move rather quickly mostly in the radial direction towards the pSS ring in the first phase, followed by a slow migration along the pSS ring towards stable equilibrium positions in the second phase.
Furthermore, the numerical simulations for hyperelastic particles with a larger blockage ratio
$(=0.2)$
showed that, as the particles become more deformable, their focusing pattern shifts from four-point focusing on the MEP (pattern C), to eight-point focusing on the MEP and DEP (pattern D), to four-point focusing on the DEP (pattern B), and to single-point focusing on the channel centre (pattern A). The particle focusing positions approach the channel centre with increasing particle deformability. These transitions of the particle focusing pattern can be accounted for by an increase in the particle deformation-induced lift, acting towards the lower shear.
In addition, we confirmed that the migration behaviour observed for hydrogel particles qualitatively agrees with previously reported capsule migration (Krüger et al. Reference Krüger, Kaoui and Harting2014b
; Schaaf & Stark Reference Schaaf and Stark2017). While quantitative comparison is complicated by differences in deformation metrics and particle models, key features such as the two-phase migration process and the transition of focusing patterns with particle deformability show remarkable qualitative similarity. Notably, both studies observe particle focusing moving towards the channel centre with increasing deformability, and a corresponding variation in equilibrium positions depending on flow conditions. In the present study, we have constructed a phase diagram showing the particle focusing patterns as a function of
${Re}$
and
$Ca$
(and
${Re}$
and
$La$
in the Supplementary material).
Moreover, our study extends these findings by performing detailed stability analyses and theoretical modelling through asymptotic expansions with two perturbed parameters,
${\textit{Re}}_p$
and
$Ca$
, suggesting that the equilibrium position is governed by the balance between inertial lift and deformation-induced lift. Extra numerical simulations were performed to examine the stability of the CEP, and the predicted critical value from the theoretical analysis showed good agreement with the boundary of CEP and DEP regions in the phase diagram. The boundary between the DEP and MEP at high
${Re}$
(or higher
$La\sim 200$
in figure 10
a) is difficult to analyse using asymptotic expansions, as the assumptions of small
${\textit{Re}}_p$
and
$Ca$
no longer hold. Additionally, further investigation is required to clarify why deformable particles tend to migrate towards the DEP at low
${Re}$
, whereas rigid particles focus at the MEP under similar conditions. These two points/issues are beyond the scope of the present study and are left for future work.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.10574.
Funding
This research was partially supported by JSPS KAKENHI Grant Numbers JP20H02072, and JP24K00809. YH was supported by JST SPRING, Grant Number JPMJSP2138. Part of the results was obtained by supercomputer Fugaku at the RIKEN R-CCS through the HPCI System Research Project (Project ID: hp220106), SQUID at D3 centre, Osaka University, and GENKAI at Research Institute for Information Technology, Kyushu University.
Declaration of interests
The authors report no conflict of interest.
Appendix A.
The numerical simulations of neutrally buoyant rigid particle were conducted using the immersed boundary method (IBM) based on Kajishima et al. (Reference Kajishima, Takiguchi, Hamasaki and Miyake2001). The rigid density is the same as the fluid one. Here, we will outline the numerical algorithm applicable only when the particle is spherical and the density
$\rho$
is spatially uniform and time-invariant.
The governing equations are
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0, \end{align}
\begin{align} \rho \frac {\mathrm{D} \boldsymbol{u}}{\mathrm{D} t}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }+\boldsymbol{f}_p, \end{align}
where
$\mathrm{D}/\mathrm{D}t=\partial /\partial t+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }$
is the material derivative, and
$\boldsymbol{f}_p$
is the body force for coupling the rigid and fluid motions. Note that although we use the method to fix the
$x$
-coordinate of the particle centroid in the computational space as described in § 2.3, we here write equations of the physical space in the inertial coordinates (where the particle freely moves) for simplicity of explanation. Additionally, the stress tensor
$\boldsymbol{\sigma }$
includes the influence of driving pressure
$\mathrm{d}p/\mathrm{d}x_{\textit{drive}}$
. The IBM consists of two sequential procedures in the progression from the
$n$
th time step to the
$(n+1)$
th one.
In the first procedure, the entire domain is regarded as a fluid. To satisfy (A1) and (A2) with
$\boldsymbol{f}_p=0$
, the velocity
$\hat {\boldsymbol{u}}$
is found by the SMAC algorithm
\begin{align} \tilde {\boldsymbol{u}}&=\boldsymbol{u}^n+\Delta t\left (-\frac {3}{2}\bigl (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\bigr )^n+\frac {1}{2}\bigl (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\bigr )^{n-1}-\frac {1}{\rho }\boldsymbol{\nabla } p^n \right .\nonumber\\&\quad\left . +\,\frac {1}{\rho }\left \{\frac {3}{2}\left (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{\prime}\right )^n-\frac {1}{2}\big(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{\prime}\big)^{n-1}\right \}\right ), \end{align}
\begin{align}\nabla^{2}\delta p&=\frac {\rho }{\Delta t}\boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}},\end{align}
\begin{align}\hat {\boldsymbol{u}}&=\tilde {\boldsymbol{u}}-\frac {\Delta t}{\rho }\boldsymbol{\nabla }\delta p,\end{align}
where
$\tilde {\boldsymbol{u}}$
is the unprojection velocity,
$\delta p=p^{n+1}-p^n$
is the pressure increment and
$\boldsymbol{\sigma }^{\prime}=\boldsymbol{\sigma }-\mathrm{tr}(\boldsymbol{\sigma })\boldsymbol{I}/3$
is the deviatoric stress tensor. In (A3), the second-order Adams–Bashforth (AB2) scheme is applied to
$\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}$
and
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{\prime}$
.
In the second procedure, the rigid and fluid motions are coupled. The AB2 scheme updates the particle centroid position
$\boldsymbol{x}_p$
in a way,
\begin{equation} \boldsymbol{x}_p^{n+1}=\boldsymbol{x}_p^n+\Delta t\left (\frac {3}{2}\boldsymbol{v}_p^n-\frac {1}{2}\boldsymbol{v}_p^{n-1}\right ), \end{equation}
where
$\boldsymbol{v}_p$
denotes the translational velocity of the particle. To have the same smoothness across the interface as the present fluid-structure interaction simulation (see § 2.3), the rigid volume fraction
$\alpha$
at the location
$\boldsymbol{x}$
is given by
\begin{equation} \alpha =\frac {1}{2}\left (1-\textrm {tanh}\frac {\beta (|\boldsymbol{r}_p|-d/2)}{\Delta x}\right ), \end{equation}
where
$\boldsymbol{r}_p = \boldsymbol{x}-\boldsymbol{x}_p$
is the relative position vector to the particle centroid,
$d$
is the sphere diameter and
$\beta (=2)$
is the sharpness parameter relative to the grid width
$\Delta x$
. From (A2)–(A5),
$\boldsymbol{u}^{n+1}$
is determined as
\begin{equation} \boldsymbol{u}^{n+1}=\hat {\boldsymbol{u}}+\frac {\Delta t}{\rho }\boldsymbol{f}_p. \end{equation}
To impose the rigid motion inside the particle, the coupling force
$\boldsymbol{f}_p$
in (A8) is
\begin{equation} \boldsymbol{f\!}_p=\alpha ^{n+1}\rho \frac {\boldsymbol{u}_p^{n+1}-\hat {\boldsymbol{u}}}{\Delta t}, \end{equation}
where
$\boldsymbol{u}_p^{n+1}$
is the particle velocity at the
$(n+1)$
th time step given by
\begin{equation} \boldsymbol{u}_p^{n+1}=\boldsymbol{v}_p^{n+1}+\boldsymbol{\omega }_p^{n+1}\times \boldsymbol{r}_p^{n+1}, \end{equation}
where
$\boldsymbol{\omega }_p$
denotes the rotational velocity of the particle. In the subsequent development, we will clarify how to identify
$\boldsymbol{v}_p^{n+1}$
and
$\boldsymbol{\omega }_p^{n+1}$
. The equations of motion for the spherical particle are
\begin{align} m_p\frac {\mathrm{d}\boldsymbol{v}_p}{\mathrm{d} t}&=\oint _{S_p}\boldsymbol{n}_p\boldsymbol{\cdot }\boldsymbol{\sigma }\mathrm{d} S, \end{align}
\begin{align} \boldsymbol{I}_{p}\boldsymbol{\cdot }\frac {\mathrm{d} \boldsymbol{\omega }_p}{\mathrm{d} t}&=\oint _{S_p}\boldsymbol{r}_{p}\times (\boldsymbol{n}_p\boldsymbol{\cdot }\boldsymbol{\sigma })\mathrm{d} S, \end{align}
where
$m_p (=\pi \rho d^3/6)$
is the particle mass,
$S_p$
is the particle surface,
$\boldsymbol{n}_p$
is the outward unit normal vector on the surface and
$\boldsymbol{I}_p (=\pi \rho d^5 \boldsymbol{I}/60)$
is the inertia tensor. As long as the grid resolution is sufficiently fine, the rigid volume fraction
$\alpha$
possesses the nature of the Heaviside function, and thus
$m_p$
and
$\boldsymbol{v}_p$
may be written as
\begin{equation} m_p=\rho \int _V \alpha \, \mathrm{d}V,\quad \boldsymbol{v}_p=\frac {\int _V\alpha \boldsymbol{u}\mathrm{d}V}{\int _V\alpha\mathrm{d}V}=\frac {\int _V\alpha \boldsymbol{u}_p\mathrm{d}V}{\int _V\alpha\mathrm{d}V}, \end{equation}
where
$V$
is the volume of the entire domain. To be treated in the IBM, (A11) is approximated with the grid scale smoothness. Introducing a smoothed delta function
$\delta =\left |\boldsymbol{\nabla }\alpha \right |$
, which is non-zero only near the interface, we express the normal vector
$\boldsymbol{n}_p$
and the surface integral as
\begin{equation} \boldsymbol{n}_p=-\frac {\boldsymbol{\nabla }\alpha }{\delta },\quad \oint _{S_p}\ldots\mathrm{d}S=\int _V\ldots\delta \, \mathrm{d}V. \end{equation}
From (A14), an arbitrary tensor
$\boldsymbol{T}$
satisfies
\begin{equation} \oint _{S_p}\boldsymbol{n}_p\boldsymbol{\cdot }\boldsymbol{T}\mathrm{d}S=-\int _V\boldsymbol{\nabla }\alpha \boldsymbol{\cdot }\boldsymbol{T}\mathrm{d}V=-\int _V\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\alpha \boldsymbol{T}\right )\mathrm{d}V+\int _V\alpha \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{T}\mathrm{d}V. \end{equation}
Applying the divergence theorem to the first term on the right-hand side of (A15) gives
\begin{equation} \int _V\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\alpha \boldsymbol{T}\right )\mathrm{d}V=\oint _{\partial V}\boldsymbol{n}\boldsymbol{\cdot }(\alpha \boldsymbol{T})\mathrm{d}S=0. \end{equation}
Here, we supposed no particle (i.e.
$\alpha =0$
) on the boundary
$\partial V$
of the entire domain. From (A2), (A15) and (A16), the right-hand side of (A11) is rewritten in a volume integral form
\begin{equation} \oint _{S_p}\boldsymbol{n}_p\boldsymbol{\cdot }\boldsymbol{\sigma }\, \mathrm{d} S =\int _{V}\alpha \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }\, \mathrm{d}V =\int _{V}\alpha \rho \frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D}t}\, \mathrm{d}V -\int _{V}\alpha \boldsymbol{f}_p\, \mathrm{d}V. \end{equation}
Using the conservation of the volume fraction
$\mathrm{D}\alpha /\mathrm{D}t=0$
and the Lagrangian time derivative
$\mathrm{d}/\mathrm{d}t=\partial /\partial t+\boldsymbol{v}_p\boldsymbol{\cdot }\boldsymbol{\nabla }$
along the particle centroid together with (A1), the first term on the right-hand side of (A17) becomes
\begin{equation} \int _V\alpha \rho \frac {\mathrm{D}\boldsymbol{u}}{\mathrm{D}t}\,\mathrm{d}V =\rho \int _V \frac {\mathrm{d}(\alpha \boldsymbol{u})}{\mathrm{d}t}\, \mathrm{d}V +\underbrace {\rho \int _V \boldsymbol{\nabla }\boldsymbol{\cdot }\big(\alpha (\boldsymbol{u}-\boldsymbol{v}_p)\boldsymbol{u}\big)\, \mathrm{d}V}_{=0}. \end{equation}
On the right-hand side of (A18), the first term may be written as
\begin{equation} \rho \int _V \frac {\mathrm{d}(\alpha \boldsymbol{u})}{\mathrm{d}t}\, \mathrm{d}V =\rho \frac {\mathrm{d}}{\mathrm{d}t}\int _V \alpha \boldsymbol{u}\, \mathrm{d}V \end{equation}
owing to the time-invariance of the domain in the inertial coordinates. Further, the second term vanishes because of (A16). The relations
$\boldsymbol{f}_p \ne 0$
and
$\alpha =1$
are valid only inside the particle. Therefore, the translational motion equation (A11) with (A13) and (A17)–(A19) reduces into
\begin{equation} \int _V\boldsymbol{f}_p\, \mathrm{d}V=0. \end{equation}
From (A9), (A13) and (A20), we arrive at the expression of the translational velocity at the
$(n+1)$
th time step
\begin{equation} \boldsymbol{v}_p^{n+1}=\frac {\int _V\alpha ^{n+1}\boldsymbol{u}_p^{n+1}\,\mathrm{d}V}{\int _V\alpha ^{n+1}\,\mathrm{d}V}=\frac {\int _V\alpha ^{n+1}\hat {\boldsymbol{u}}\mathrm{d}V}{\int _V\alpha ^{n+1}\,\mathrm{d}V}. \end{equation}
Through a similar process to the derivation of (A20), the rotational motion equation (A12) reduces into
\begin{equation} \int _V\boldsymbol{r}_p^{n+1}\times \boldsymbol{f}_p\,\mathrm{d}V=0. \end{equation}
From (A9), (A10) and (A22), we can express the
$i$
th component of the rotational velocity at the
$(n+1)$
th time step in Einstein’s index notation
\begin{equation} \omega _{p,i}^{n+1}=\frac {\int _V\alpha ^{n+1}\big(\varepsilon _{ijk}r^{n+1}_{p,j}\hat {u}_k\big)\,\mathrm{d}V}{\int _V\alpha ^{n+1}\big (\big |\boldsymbol{r}_{p}^{n+1}\big |^2-\big (r_{p,i}^{n+1}\big )^2\big)\,\mathrm{d}V}\quad \left (\text{no sum. with respect to }i\right ), \end{equation}
where
$\varepsilon _{ijk}$
denotes the Levi–Civita symbol. The expressions (A10), (A21) and (A23) indicate that the velocity
$\boldsymbol{u}_p^{n+1}$
of the neutrally buoyant particle depends only on
$\hat {\boldsymbol{u}}$
,
$\boldsymbol{x}_p^{n+1}$
and
$\alpha ^{n+1}$
(which are given by (A5), (A6) and (A7), respectively), and is identified without time-integral operation for
$\boldsymbol{v}_p$
and
$\boldsymbol{\omega }_p$
. Once these quantities are known, we can determine the coupling force
$\boldsymbol{f}_p$
by (A9) and update
$\boldsymbol{u}^{n+1}$
by (A8).
In the numerical simulation with the rigid sphere fixed within the cross-section in § 3.3.2, the particle obeys (A22) and freely rotates. Further, it freely translates in the
$x$
-direction despite
$\boldsymbol{e}_y\boldsymbol{\cdot }\boldsymbol{v}_p = 0$
and
$\boldsymbol{e}_z\boldsymbol{\cdot }\boldsymbol{v}_p = 0$
. For such a translational motion, the reduced equation is
\begin{equation} \boldsymbol{e}_x\boldsymbol{\cdot }\int _V\boldsymbol{f}_p\, \mathrm{d}V=0. \end{equation}
Thus,
$\boldsymbol{x}_p^{n+1}$
,
$\alpha ^{n+1}$
,
$\boldsymbol{v}_p^{n+1}$
,
$\boldsymbol{\omega }_p^{n+1}$
,
$\boldsymbol{f}_p$
and
$\boldsymbol{u}^{n+1}$
are determined in the above-mentioned manner except that (A21) is replaced with
\begin{equation} \boldsymbol{v}_p^{n+1}=\boldsymbol{e}_x\frac {\int _V\alpha ^{n+1}\boldsymbol{e}_x\boldsymbol{\cdot }\hat {\boldsymbol{u}}\,\mathrm{d}V}{\int _V\alpha ^{n+1}\mathrm{d}V}. \end{equation}
In the steady state with
${\mathrm{d}\boldsymbol{v}_p}/{\mathrm{d} t}=0$
, from (A13), (A17)–(A19) and (A24), we may estimate the cross-sectional force
$\boldsymbol{F}^{\perp }$
(see § 3.3.2) as
\begin{equation} \boldsymbol{F}^{\perp }=\boldsymbol{P}^{\perp }\boldsymbol{\cdot }\oint _{S_p}\boldsymbol{n}_p\boldsymbol{\cdot }\boldsymbol{\sigma }\, \mathrm{d} S =-\boldsymbol{P}^{\perp }\boldsymbol{\cdot }\int _V\boldsymbol{f}_p\, \mathrm{d}V=-\int _V\boldsymbol{f}_p\, \mathrm{d}V, \end{equation}
where
$\boldsymbol{P}^{\perp }=\boldsymbol{I}-\boldsymbol{e}_x\boldsymbol{e}_x$
is the projection tensor onto the cross-section.
Open access
