1. Introduction
In a pipe flow at a finite channel (or particle) Reynolds number
$Re$
(
$Re_p$
), a rigid spherical particle exhibits migration perpendicular to the flow direction, originally reported by Segre & Silberberg (Reference Segre and Silberberg1962), the so-called inertial focusing or tubular pinch effect, where the particles equilibrate at a distance from the channel centreline as a consequence of the force balance between the shear-induced and wall-induced lift forces. The phenomenon is of fundamental importance in microfluidic techniques such as label-free cell alignment, sorting and separation techniques (Martel & Toner Reference Martel and Toner2014; Warkiani et al. Reference Warkiani, Khoo, Wu, Tay, Bhagat, Han and Lim2016; Zhou et al. Reference Zhou, Mukherjee, Gao, Luan and Papautsky2019). Although the techniques allow us to reduce the complexity and cost of clinical applications by using small amounts of blood samples, archetypal inertial focusing system requires steady laminar flow through long channel distances
$L_f$
, which can be estimated as
$L_f = \pi \mathcal {H}/(Re_p f_l)$
, where
$\mathcal {H}$
is the dimension of the channel (or its hydraulic diameter) and
$f_l$
is a non-dimensional lift coefficient (Di Carlo Reference Di Carlo2009). So far, various kinds of geometries have been proposed to achieve the required distance for inertial focusing in a compact space, e.g. sinusoidal, spiral and hybrid channels (Bazaz et al. Reference Bazaz, Mashhadian, Ehsani, Saha, Krüger and Warkiani2020). Without increasing
$Re_p$
, the recent experimental study by Mutlu et al. (Reference Mutlu, Edd and Toner2018) achieved inertial focusing of
$0.5$
$\mu\rm m$
particles (
$Re_p \approx 0.005$
) in short channels by using oscillatory channel flows. Since the oscillatory flows allow a suspended particle to increase its total travel distance without net displacement along the flow direction, using oscillatory flow can be an alternative and practical strategy for inertial focusing in microfluidic devices. Recently, Vishwanathan & Juarez (Reference Vishwanathan and Juarez2021) experimentally investigated the effects of the Womersley number (
$\alpha$
) on inertial focusing in planar pulsatile flows, and evaluated the lateral migration (or off-centre focusing) speed on a small and weakly inertial particle for different oscillatory frequencies. They concluded that inertial focusing is achieved in only a fraction of the channel length (
$1\,\%$
–
$10$
%) compared with what would be required in a steady flow (Vishwanathan & Juarez Reference Vishwanathan and Juarez2021). Sun et al. (Reference Sun, Yu and Shao2009) performed two-dimensional (2-D) simulations of a neutrally buoyant circular particle in oscillatory pressure-driven channel flows for
$Re \geqslant 50$
. Their results indicated that lower oscillatory frequency makes the equilibrium position closer to the channel centreline, while higher oscillatory frequency maintains the equilibrium positions similarly to the steady flow conditions. However, it remains unknown whether the equilibrium position of deformable capsules under pulsatile channel flows can be formulated in the same context as that of a rigid circular particle.
While a number of studies have analysed the off-centre focusing of rigid spherical particles under steady flow by a variety of approaches, such as analytical calculations (Ho & Leal Reference Ho and Leal1974; Schonberg & Hinch Reference Schonberg and Hinch1989; Asmolov Reference Asmolov1999), numerical simulations (Bazaz et al. Reference Bazaz, Mashhadian, Ehsani, Saha, Krüger and Warkiani2020; Feng et al. Reference Feng, Hu and Joseph1994; Shao et al. Reference Shao, Yu and Sun2008; Yang et al. Reference Yang, Wang, Joseph, Hu, Pan and Glowinski2005; Yu et al. Reference Yu, Phan-Thien and Tanner2004) and experimental observations (Di Carlo Reference Di Carlo2009; Karnis et al. Reference Karnis, Goldsmith and Mason1966; Matas et al. Reference Matas, Morris and Guazzelli2004), the inertial focusing of deformable particles such as biological cells, consisting of an internal fluid enclosed by a thin membrane, has not yet been fully described, especially under unsteady flows. Due to their deformability, the problem of inertial focusing of deformable particles is more complex than with rigid spherical particles, as originally reported by Segre & Silberberg (Reference Segre and Silberberg1962). It is now well known that a deformable particle at low
$Re$
migrates towards the channel axis under steady laminar flow (Karnis et al. Reference Karnis, Goldsmith and Mason1963). Hereafter, we call this phenomenon ‘axial focusing’. A recent numerical study showed that, in almost inertialess condition, the axial focusing of a deformable spherical capsule can be accelerated by the flow pulsation at a specific frequency (Takeishi & Rosti Reference Takeishi and Rosti2023). For finite
$Re$
(
${\gt } 1$
), however, it is still uncertain whether the flow pulsation can enhance the off-centre focusing or impede it (i.e. axial focusing). Therefore, the primary objective of this study is to clarify whether a capsule lateral movement at finite
$Re$
in a pulsatile channel flow can be altered by its deformability. The second objective is to clarify whether the
$Re$
-dependent equilibrium radial position of a capsule in a channel or travelling time are controllable by oscillatory frequency.
At least for steady channel flows, inertial focusing of deformable capsules including biological cells have been investigated in recent years both by means of experimental observations (Warkiani et al. Reference Warkiani, Khoo, Wu, Tay, Bhagat, Han and Lim2016; Zhou et al. Reference Zhou, Mukherjee, Gao, Luan and Papautsky2019) and numerical simulations (Raffiee et al. Reference Raffiee, Dabiri and Ardekani2017; Schaaf & Stark Reference Schaaf and Stark2017; Takeishi et al. Reference Takeishi, Yamashita, Omori, Yokoyama, Wada and Sugihara-Seki2022). For instance, Hur et al. (Reference Hur, Henderson-MacLennan, McCabe and Di Carlo2011) experimentally investigated the inertial focusing of various cell types (including red blood cells, leukocytes and cancer cells such as a cervical carcinoma cell line, breast carcinoma cell line and osteosarcoma cell line) with a cell-to-channel size ratio
$0.1 \leqslant { d}_0/W \leqslant 0.8$
, using a rectangular channel with a high aspect ratio of
$W/H \approx 0.5$
, where
${ d}_0$
,
$W$
and
$H$
are the cell equilibrium diameter, channel width and height, respectively. They showed that the cells can be separated according to their size and deformability (Hur et al. Reference Hur, Henderson-MacLennan, McCabe and Di Carlo2011). The experimental results can be qualitatively described using a spherical capsule (Kilimnik et al. Reference Kilimnik, Mao and Alexeev2011) or droplet model (Chen et al. Reference Chen, Xue, Zhang, Hu, Jiang and Sun2014). In more recent experiments by Hadikhani et al. (Reference Hadikhani, Hashemi, Balestra, Zhu, Modestino, Gallaire and Psaltis2018), the authors investigated the effect of
$Re$
(
$1 \lt Re \lt 40$
) and capillary number
$Ca$
– ratio between the fluid viscous force and the membrane elastic force – (
$0.1 \lt Ca \lt 1$
) on the lateral equilibrium of bubbles in rectangular microchannels and different bubble-to-channel size ratios with
$0.48 \leqslant { d}_0/W \leqslant 0.84$
. The equilibrium position of such soft particles results from the competition between
$Re$
and
$Ca$
, because high
$Re$
induces off-centre focusing, while high
$Ca$
, i.e. high deformability, allows axial focusing. However, notwithstanding these recent advancements, a comprehensive understanding of the effect on the inertial focusing of these two fundamental parameters has not been fully established yet.
Numerical analysis more clearly showed that the ‘deformation-induced lift force’ becomes stronger as the particle deformation is increased (Raffiee et al. Reference Raffiee, Dabiri and Ardekani2017; Schaaf & Stark Reference Schaaf and Stark2017). Although a number of numerical analyses regarding inertial focusing have been reported in recent years mostly for spherical particles (Bazaz et al. Reference Bazaz, Mashhadian, Ehsani, Saha, Krüger and Warkiani2020; Banerjee et al. Reference Banerjee, Rosti, Kumar, Brandt and Russom2021), the equilibrium positions of soft particles is still debated owing to the complexity of the phenomenon. Kilimnik et al. (Reference Kilimnik, Mao and Alexeev2011) showed that the equilibrium position in a cross-section of a rectangular microchannel with
${ d}_0/H = 0.2$
shifts towards the wall as
$Re$
increases from
$1$
to
$100$
. Schaaf & Stark (Reference Schaaf and Stark2017) also performed numerical simulations of spherical capsules in a square channel for
$0.1 \leqslant { d}_0/H \leqslant 0.4$
and
$5 \leqslant Re \leqslant 100$
without viscosity contrast, and showed that the equilibrium position was nearly independent of
$Re$
. In a more recent numerical analysis by Alghalibi et al. (Reference Alghalibi, Rosti and Brandt2019), simulations of a spherical hyperelastic particle in a circular channel with
${d}_0/D = 0.2$
were performed with
$100 \leqslant Re \leqslant 400$
and Weber number (
$We$
) with
$0.125 \leqslant We \leqslant 4.0$
, the latter of which is the ratio of the inertial effect to the elastic effect acting on the particles. Their numerical results showed that regardless of
$Re$
, the final equilibrium position of a deformable particle is the centreline and harder particles (i.e. with lower
$We$
) tended to rapidly migrate towards the channel centre (Alghalibi et al. Reference Alghalibi, Rosti and Brandt2019). The behaviour of a capsule subjected to pulsatile channel flow was addressed in the pioneering work by Maestre et al. (Reference Maestre, Pallaresa, Cuestaa and Scott2019), where the migration velocity during axial focusing was investigated at
$O(Re) \leqslant 10^{-2}$
and
${ d}_0/D = 0.5$
for
$Ca = 0.075$
–
$1.2$
. Despite these efforts, the inertial focusing of capsules subjected to pulsatile flow at finite inertia cannot be estimated based on these achievements.
Aiming for the precise description of the inertial focusing of spherical capsules in pulsatile channel flows, we thus perform numerical simulations of individual capsules with a major diameter of
${ d}_0 = 2a_0 = 8$
$\mu\rm m$
in a cylindrical microchannel with
$D = 2R = 20$
–
$50$
$\mu\rm m$
(i.e.
$R/a_0 = 2.5$
–
$6.25$
) for a wide range of oscillatory frequencies. Each capsule membrane, following the Skalak constitutive (SK) law (Skalak et al. Reference Skalak, Tozeren, Zarda and Chien1973), is simulated for different
$Re$
,
$Ca$
and size ratio
$R/a_0$
. Since this problem requires heavy computational resources, we resort to GPU computing, using the lattice-Boltzmann method (LBM) for the inner and outer fluids, and the finite-element method (FEM) to describe the deformation of the capsule membrane. This model has been successfully applied in the past for the analysis of the capsule flow in circular microchannels (Takeishi et al. Reference Takeishi, Yamashita, Omori, Yokoyama, Wada and Sugihara-Seki2022; Takeishi & Rosti Reference Takeishi and Rosti2023). The remainder of this paper is organised as follows. Section 2 gives the problem statement and numerical methods, and § 3 presents the numerical results for a single spherical capsule. Finally, a summary of the main conclusions is reported in § 4. A description of numerical verifications is presented in the Appendix.
2. Problem statement
2.1. Flow and capsule models and set-up
We consider the motion of an initially spherical capsule with diameter
${ d}_0$
(= 2
$a_0$
= 8
$\mu\rm m$
) flowing in a circular channel diameter
$D$
(= 2
$R$
= 20–50
$\mu\rm m$
), with a resolution of 32 fluid lattices per capsule diameter
${ d}_0$
. The channel length is set to be 20
$a_0$
, following a previous numerical study (Takeishi et al. Reference Takeishi, Yamashita, Omori, Yokoyama, Wada and Sugihara-Seki2022). Although we have investigated in the past the effect of the channel length
$L$
and the mesh resolutions on the trajectory of the capsule centroid (see figure 7 of Takeishi & Rosti Reference Takeishi and Rosti2023), we further assess the effect of this length on the lateral movement of a capsule in Appendix A (figure 12
$a$
).
Visualisation of a spherical capsule with radius
$a_0$
in a channel with radius of
$R$
(
$R/a_0$
= 2.5) under a pulsatile flow with velocity
$V^\infty$
, which can be decomposed into the steady parabolic flow
$V_0^\infty$
and the oscillatory flow
$V_{ {osci}}^\infty$
in the absence of any capsule. The capsule, initially placed at off-centre radial position
$r^\ast _{ {c0}}= r_{ {c0}}/R = 0.4$
, travels in the radial direction. In the figure, the lengths travelled by the capsule in the flow (
$z$
) direction are not to scale for illustrative purpose. Hereafter, the same modification will be applied for visualisation.
The capsule consists of a Newtonian fluid enclosed by a thin elastic membrane, sketched in figure 1. The membrane is modelled as an isotropic and hyperelastic material following the SK law (Skalak et al. Reference Skalak, Tozeren, Zarda and Chien1973), in which the strain energy
$w_{\mathrm {SK}}$
and principal tensions in the membrane
$\tau _1$
and
$\tau _2$
(with
$\tau _1 \geqslant \tau _2$
) are given by
\begin{equation} \frac {w_{\textit {SK}}}{G_s} = \frac {1}{4} \left ( I_1^2 + 2I_1 - 2I_2 + C I_2^2\right ) \end{equation}
and
\begin{equation} \frac {\tau _i}{G_s} = \frac {\eta _i}{\eta _j} \left [ \eta _i^2 - 1 + C \left ( \eta _i^2 \eta _j^2 - 1 \right ) \right ] \quad \text {for} \ (i, j) = (1, 2) \ \text {or} \ (2, 1). \end{equation}
Here,
$w_{\textit {SK}}$
is the strain energy density function,
$G_s$
is the membrane shear elastic modulus,
$C$
is a coefficient representing the area incompressibility, and
$I_1$
(
$= \eta _1^2 + \eta _2^2 - 2$
) and
$I_2$
(
$= \eta _1^2\eta _2^2 - 1$
) are the invariants of the strain tensor, with
$\eta _1$
and
$\eta _2$
being the principal extension ratios. In the SK law (2.1), the area dilation modulus is
$K_s = G_s (1 + 2C)$
. In this study, we set
$C = 10^2$
(Barthès-Biesel et al. Reference Barthès-Biesel, Diaz and Dheni2002), which describes an almost incompressible membrane. Bending resistance is also considered (Li et al. Reference Li, Dao, Lim and Suresh2005), with a bending modulus
$k_b = 5.0 \times 10^{-19}$
J (Puig-de-Morales-Marinkovic et al. Reference Puig-de-Morales-Marinkovic, Turner, Butler, Fredberg and Suresh2007). These values have been shown to successfully reproduce the deformation of red blood cells in shear flow (Takeishi et al. Reference Takeishi, Imai, Nakaaki, Yamaguchi and Ishikawa2014, Reference Takeishi, Rosti, Imai, Wada and Brandt2019) and the thickness of the cell-depleted peripheral layer in circular channels (see figure A.1 of Takeishi et al. Reference Takeishi, Imai, Nakaaki, Yamaguchi and Ishikawa2014). Neglecting inertial effects on the membrane deformation, the static local equilibrium equation of the membrane is given by
\begin{equation} \nabla _s \cdot \boldsymbol {\tau } + \boldsymbol {q} = \boldsymbol {0}, \end{equation}
where
$\nabla _s (= ( \boldsymbol {I} - \boldsymbol {n} \boldsymbol {n} ) \cdot \nabla )$
is the surface gradient operator,
$\boldsymbol {n}$
is the unit normal outward vector in the deformed state,
$\boldsymbol {q}$
is the load on the membrane and
$\boldsymbol {\tau }$
is the in-plane elastic tension that is obtained using the SK law (2.1).
The fluids are modelled with the incompressible Navier–Stokes equations for the fluid velocity
$\boldsymbol {v}$
:
\begin{equation} \rho \left ( \frac {\partial \boldsymbol {v}}{\partial t} + \boldsymbol {v} \cdot \nabla \boldsymbol {v} \right ) = \nabla \cdot \boldsymbol {\sigma }^f + \rho \boldsymbol {f}, \end{equation}
\begin{equation} \nabla \cdot \boldsymbol {v} = 0, \end{equation}
with
\begin{align} \boldsymbol {\sigma }^f = -p\boldsymbol {I} + \mu \left ( \nabla \boldsymbol {v} + \nabla \boldsymbol {v}^T \right ), \end{align}
where
$\boldsymbol {\sigma }^f$
is the total stress tensor of the flow,
$p$
is the pressure,
$\rho$
is the fluid density,
$\boldsymbol {f}$
is the body force and
$\mu$
is the viscosity of the liquid, expressed using a volume fraction of the inner fluid
$\mathcal {I}$
(0
$\leqslant \mathcal {I} \leqslant$
1) as
\begin{align} \mu = \left \{ 1 + \left ( \lambda - 1 \right ) \mathcal {I} \right \} \mu _0, \end{align}
where
$\lambda$
(=
$\mu _1/\mu _0$
) is the viscosity ratio,
$\mu _0$
is the external fluid viscosity and
$\mu _1$
is the internal fluid viscosity. No density contrast is considered; that is, the ratio of densities between the external and internal fluid is assumed to be one.
The dynamic condition coupling the different phases requires the load
$\boldsymbol {q}$
to be equal to the traction jump
$ ( \boldsymbol {\sigma }^f_{ {out}} - \boldsymbol {\sigma }^f_{{in}} )$
across the membrane:
\begin{align} \boldsymbol {q} = \left ( \boldsymbol {\sigma }^f_{{out}} - \boldsymbol {\sigma }^f_{ {in}} \right ) \cdot \boldsymbol {n}, \end{align}
where the subscripts ‘out’ and ‘in’ represent the outer and internal regions of the capsule, respectively.
The flow in the channel is sustained by a uniform pressure gradient
$\partial p_0/\partial z\,(= \partial _z p_0)$
, which can be related to the maximum fluid velocity in the channel by
$\partial _z p_0 = -4 \mu _0 V_{{max}}^\infty /R^2$
. The pulsation is given by a superimposed sinusoidal function, such that the total pressure gradient is
\begin{equation} \partial _z p (t) = \partial _z p_0 + \partial _z p_a \sin {\left ( 2 \pi f t \right )}. \end{equation}
The problem is governed by six main non-dimensional numbers, including: (
$i$
) the Reynolds number
$Re$
and (
$i\hspace {-1.0pt}i$
) the capillary number
$Ca$
defined as
\begin{align} &Re = \frac {\rho D V_{{max}}^{\infty }}{\mu _0}, \end{align}
\begin{align} &Ca = \frac {\mu _0 \dot {\gamma }_{m} a_0}{G_s} = \frac {\mu _0 V_{ {max}}^{\infty }}{G_s} \frac {a_0}{4 R}, \end{align}
where
$V_{ {max}}^{\infty }$
(
$= 2V_{ {m}}^{\infty })$
is the maximum fluid velocity in the absence of any cells,
$V_{ {m}}^{\infty }$
is the mean fluid velocity and
$\dot {\gamma }_{m}$
(
$= V_{ {m}}^{\infty }/D)$
is the mean shear rate. Note that increasing
$Re$
under constant
$Ca$
corresponds to increasing
$G_s$
, namely, a harder capsule. Furthermore, we have: (
$i\hspace {-1.0pt}i\hspace {-1.0pt}i)$
the viscosity ratio
$\lambda$
; (
$i\hspace {-1.0pt}v$
) the size ratio
$R/a_0$
; (
$v$
) the non-dimensional pulsation frequency
$f^\ast = f/\dot {\gamma }_{m}$
and (
$v\hspace {-1.0pt}i$
) the non-dimensional pulsation amplitude
$\partial _z p_a^\ast = \partial _z p_a/\partial _z p_0$
. Considering the focus of this study, we decided to primarily investigate the effect of
$Re$
,
$R/a_0$
and
$f^\ast$
. Representative rigid and largely deformable capsules are considered with
$Ca = 0.05$
and
$Ca = 1.2$
, respectively.
When presenting the results, we will initially focus on the analysis of lateral movements of the capsule in effectively inertialess condition (
$Re = 0.2$
) for
$R/a_0 = 2.5$
, and later consider variations of the size ratio
$R/a_0$
, viscosity ratio
$\lambda$
, Reynolds number
$Re$
(
${\gt } 1$
) and
$Ca$
. We confirmed that the flow at
$Re = 0.2$
well approximates an almost inertialess flow for single- (Takeishi & Rosti Reference Takeishi and Rosti2023) and multi-cellular flow (Takeishi et al. Reference Takeishi, Rosti, Imai, Wada and Brandt2019). Unless otherwise specified, we show the results obtained with
$\partial _z p_a^\ast = 2$
and
$\lambda = 1$
.
2.2. Numerical simulation
The governing equations for the fluid are discretised by the LBM based on the D3Q19 model (Chen & Doolen Reference Chen and Doolen1998). We track the Lagrangian points of the membrane material points
$\boldsymbol {x}_m (\boldsymbol {X}_m,t)$
over time, where
$\boldsymbol {X}_m$
is a material point on the membrane in the reference state. Based on the virtual work principle, the above strong-form equation (2.3) can be rewritten in weak form as
\begin{equation} \int _S \boldsymbol {\hat {u}} \cdot \boldsymbol {q} {\rm d}S = \int _S \boldsymbol {\hat {\epsilon }} : \boldsymbol {\tau } {\rm d}S, \end{equation}
where
$S$
is the surface area of the capsule membrane, and
$\boldsymbol {\hat {u}}$
and
$\boldsymbol {\hat {\epsilon }} = ( \nabla _s \boldsymbol {\hat {u}} + \nabla _s \boldsymbol {\hat {u}}^T )\big /2$
are the virtual displacement and virtual strain, respectively. The FEM is used to solve (2.12) and obtain the load
$\boldsymbol {q}$
acting on the membrane (Walter et al. Reference Walter, Salsac, Barthès-Biesel and Le Tallec2010). The velocity at the membrane node is obtained by interpolating the velocities at the fluid node using the immersed boundary method (Peskin Reference Peskin2002). The membrane node is updated by Lagrangian tracking with the no-slip condition. The explicit fourth-order Runge–Kutta method is used for the time integration. The volume-of-fluid method (Yokoi Reference Yokoi2007) and front-tracking method (Unverdi & Tryggvason Reference Unverdi and Tryggvason1992) are employed to update the viscosity in the fluid lattices. A volume constraint is implemented to counteract the accumulation of small errors in the volume of the individual cells (Freund Reference Freund2007): in our simulation, the relative volume error is always maintained lower than
$1.0 \times 10^{-3}\,\%$
, as tested and validated in our previous study of cell flow in circular channels (Takeishi et al. Reference Takeishi, Imai, Ishida, Omori, Kamm and Ishikawa2016). All procedures were fully implemented on a GPU to accelerate the numerical simulation. More precise explanations for numerical simulations including membrane mechanics is provided in our previous works (see also Takeishi et al. Reference Takeishi, Rosti, Imai, Wada and Brandt2019, Reference Takeishi, Yamashita, Omori, Yokoyama, Wada and Sugihara-Seki2022).
Periodic boundary conditions are imposed in the flow direction (
$z$
-direction). No-slip conditions are employed for the walls (radial direction). We set the mesh size of the LBM for the fluid solution to
$250$
nm and that of the finite elements describing the membrane to approximately
$250$
nm (an unstructured mesh with
$5120$
elements was used for the FEM). This resolution was shown to successfully represent single- and multi-cellular dynamics (Takeishi et al. Reference Takeishi, Rosti, Imai, Wada and Brandt2019, Reference Takeishi, Yamashita, Omori, Yokoyama, Wada and Sugihara-Seki2022).
2.3. Analysis of capsule deformation
Later, we investigate the in-plane principal tension
$\tau_i$
(with
$\tau_1 \geqslant \tau_2$
) and the isotropic tension
$\tau_{ {iso}}$
in the membrane of the capsule. In the case of a two-dimensional isotropic elastic membrane, the isotropic membrane tension can be calculated by
$\tau_{ {iso}} = (\tau_1 + \tau_2)/2$
for the deformed capsule. The averaged value of
$\tau_{ {iso}}$
is then calculated as
\begin{equation} \langle \tau_{\textit {iso}} \rangle = \frac {1}{S\mathcal {T}} \int _{\mathcal {T}} \int _S \tau_{\textit {iso}} (\boldsymbol {x}_m, t) {\rm d}S {\rm d}t, \end{equation}
where
$\mathcal {T}$
is the period of the capsule motion. Hereafter,
$\langle \cdot \rangle$
denotes a spatial-temporal average. Time average starts after the trajectory has finished the initial transient dynamics, which differs for each case. For instance, at finite
$Re$
conditions, a quasi-steady state is usually attained around the non-dimensional time of
$\dot {\gamma }_{m} t = 200$
, and we start accumulating the statistics from
$\dot {\gamma }_{m} t \geqslant 400$
to fully cancel the influence of the initial conditions.
(
$a$
) Side views of the capsule during its axial focusing under steady flow for
$Ca = 0.05$
(top),
$Ca = 0.2$
(middle) and
$Ca = 1.2$
(bottom). The capsule is initially placed at
$r^\ast _{ {c0}} = 0.55$
. The coloured dot on the membrane is shown to measure the membrane rotation. (
$b$
) Time histories of the radial position of these capsule centroids
$r_{{c}}/R$
. The dashed lines are the curves given by
$r_{{c}}^\ast = C_2 \exp {(-C_1 t^\ast )}$
, where
$r_{{c}}^\ast$
(
$= r_{{c}}/R$
) is the non-dimensional capsule centroid,
$t^\ast$
(
$\dot \gamma _{m}t$
) is the non-dimensional time, and
$C_1$
and
$C_2$
are the coefficients found by a least-squares fitting to the plot. The results in the figure are obtained for
$Re= 0.2$
,
$R/a_0 = 2.5$
and
$\lambda = 1$
.
3. Results
3.1.
Axial focusing of the capsule under steady channel flow (
$Re \lt 1$
)
We first investigate the axial focusing of a capsule under steady flow, which can be assumed to be effectively inertialess (
$Re = 0.2$
). Figure 2(
$a$
) shows side views of the capsule during its axial focusing in the channel of size
$R/a_0 = 2.5$
for different
$Ca$
(
$= 0.05$
,
$0.2$
and
$1.2$
). The capsule, initially placed at
$r^\ast _{ {c0}} = r_{ {c0}}/R = 0.55$
, migrates after the flow onsets towards the channel centreline (i.e. capsule centroid is
$r_{{c}} = 0$
) while deforming, finally reaching its equilibrium position at the centreline where it achieves an axial-symmetric shape. Although the magnitude of deformation during axial focusing depends on
$Ca$
, these processes are commonly observed for every
$Ca$
. The time history of the radial position of the capsule centroid
$r_{{c}}$
is shown in figure 2(
$b$
). The results clearly show that the speed of axial focusing grows with
$Ca$
. Interestingly, all trajectories are well fitted by the following empirical expression:
\begin{equation} r_{{c}}^\ast = C_2 \exp {(-C_1 t^\ast )}, \end{equation}
where
$t^\ast$
(
$= \dot \gamma _{m}t$
) is the non-dimensional time, and
$C_1$
(
${\gt } 0$
) and
$C_2$
are two coefficients that can be found by a least-squares fitting to the plot. Fittings are performed using data between the initial (
$r^\ast _{ {c0}} = 0.55$
) and final state (
$\Delta x_{\textit{LBM}}/R \leqslant 0.01$
for
$R/a_0 = 2.5$
), defined as the time when the capsule is within one mesh size (
$\Delta x_{\textit{LBM}}$
) from the channel axis.
Performing time differentiation of (3.1), the non-dimensional velocity of the capsule centroid
$\dot {r}_c^\ast$
can be estimated as
\begin{equation} \dot {r}_{\textit {c}}^\ast = -C_1 r_{\textit {c}}^\ast . \end{equation}
This linear relation (3.2) may be understood by a shear-induced lift force propotional to the local shear strength. A more detailed description of the relationship between the coefficient
$C_1$
and the lift force on the capsule are provided in Appendix B.
Figure 3(
$a$
) shows the coefficient
$C_1$
as a function of
$Ca$
. As expected from figure 2(
$b$
), the value of
$C_1$
increases with
$Ca$
. Since the capsule deformability is also affected by the viscosity ratio
$\lambda$
, its influence on
$C_1$
is also investigated in figure 3(
$b$
). At a fixed
$Ca$
(
$= 1.2$
), the value of
$C_1$
decreases with
$\lambda$
.
Coefficient
$C_1$
(
$a$
) as a function of
$Ca$
for
$\lambda = 1$
and (
$b$
) as a function of
$\lambda$
for
$Ca = 1.2$
. The results are obtained with
$Re = 0.2$
,
$R/a_0 = 2.5$
and
$r^\ast _{ {c0}} = 0.55$
.
To further prove that
$C_1$
is independent of the initial radial position of the capsule centroid, additional numerical simulations are performed with a larger channel (
$R/a_0 = 5$
) for different
$r^\ast _{ {c0}}$
. Note that a case with larger channel for constant
$Re$
denotes smaller
$V_{ {max}}^\infty$
, resulting in a smaller
$G_s$
(i.e. softer capsule) for constant
$Ca$
. Figure 4(
$a$
) is one of the additional runs at
$Ca = 0.2$
, where the capsule is initially placed at
$r^\ast _{ {c0}} = 0.75$
. Figure 4(
$b$
) is the time history of the radial position of the capsule centroid
$r_{{c}}$
for different initial positions
$r^\ast _{ {c0}}$
. We observe that the exponential fitting is still applicable for these runs, with the coefficient
$C_1$
reported in figure 4(
$c$
). These results provide a confirmation that
$C_1$
is indeed independent of the initial radial position
$r_{c0}^\ast$
. Furthermore, the fitting provided in (3.1) is applicable even for a different constitutive law. Discussion of these results for capsule described by the neo-Hookean model, which features strain-softening, is reported in Appendix C (see also figure 13).
(
$a$
) Side views of a capsule with
$Ca = 1.2$
during its axial focusing for
$R/a_0 = 5$
, where the capsule is initially placed at
$r^\ast _{ {c0}} = 0.75$
. (
$b$
) Time histories of the radial position of the capsule centroids
$r_{{c}}/R$
for different initial positions
$r^\ast _{ {c0}}$
. (
$c$
) Coefficient
$C_1$
as a function of the initial position
$r^\ast _{ {c0}}$
. The results are obtained with
$\lambda = 1$
.
3.2. Capsule behaviour under pulsatile channel flow
Next, we investigate inertial focusing of capsules at finite
$Re$
and investigate whether the equilibrium radial position of the capsule can be altered by pulsations of the flow. Two representative behaviours of the capsule at low
$Ca$
(
$= 0.05$
) and high
$Ca$
(
$= 1.2$
) are shown in figure 5(
$a$
), which are obtained with
$f^\ast = 0.02$
and
$Re = 10$
. The simulations are started from an off-centre radial position
$r^\ast _{ {c0}}$
. Hereafter, we consider the viscosity ratio
$\lambda = 1$
for simplicity. At the end of the migration, the least deformable capsule (
$Ca = 0.05$
) exhibits an ellipsoidal shape with an off-centred position (figure 5
$a$
, left), while the most deformable one (
$Ca = 1.2$
) exhibits the typical parachute shape at the channel centreline (figure 5
$a$
, right). Detailed trajectories of these capsule centroids
$r_{{c}}/R$
are shown in figure 5(
$b$
), where the non-dimensional oscillatory pressure gradient
$\partial _z p^\ast (t^\ast )$
(
$= 1 + 2 \sin {(2 \pi f^\ast t^\ast )}$
) is also displayed. The least deformable capsule (
$Ca = 0.05$
) fluctuates around the off-centre position
$r_{{c}}/R$
(
$\approx 0.2$
) and the waveform of
$r_{{c}}/R$
lags behind
$\partial _z p^\ast (t^\ast )$
. The capsule with large
$Ca$
(
$= 1.2$
), however, immediately exhibits axial focusing, reaching the centreline within one flow period (figure 5
$b$
). Therefore, axial and off-centre focusing strongly depend on
$Ca$
.
(
$a$
) Side views of the capsule during its migration at each time at
$f^\ast = 0.02$
for
$Ca = 0.05$
(left; see the supplementary movie 1, available at https://doi.org/10.1017/jfm.2025.184) and
$Ca = 1.2$
(right; see the supplementary movie 2). (
$b$
and
$c$
) Time histories of (
$b$
) the radial position of these capsule centroids
$r_{{c}}/R$
and (
$c$
) isotropic tensions
$\tau_{ {iso}}$
, respectively. In panels (
$a$
)–(
$c$
), the results are obtained with
$\partial _z p_a^\ast = 2$
. (
$d$
and
$e$
) Time histories of (d)
$r_{{c}}/R$
and (e)
$\tau_{ {iso}}$
for
$\partial _z p_a^\ast = -2$
, where those in steady flow are also superposed. In panels (
$b$
)–(
$e$
), non-dimensional pressure gradient
$\partial _z p^\ast$
is also displayed on the right axis. The results are obtained with
$Re = 10$
,
$R/a_0 = 2.5$
and
$r^\ast _{ {c0}} = 0.4$
.
Figure 5(
$c$
) is the time history of the isotropic tension
$\tau_{ {iso}}$
. The major waveforms of
$\tau_{ {iso}}$
are synchronised with
$\partial _z p^\ast$
in both
$Ca = 0.05$
and
$Ca = 1.2$
, thus indicating that the membrane tension spontaneously responds to the background fluid flow. The Taylor parameter, a classical index of deformation, is described in Appendix D (see figure 14).
To clarify whether fast axial focusing depends on the phase of oscillation or not, an antiphase pulsation (i.e.
$\partial _z p_a^\ast = -2$
) is given by
$\partial _z p^\ast (t^\ast ) = 1 - 2 \sin {(2 \pi f^\ast t^\ast )}$
. Time histories of the capsule centroid
$r_{{c}}/R$
and membrane tension
$\tau_{\textit {iso}}$
under such conditions are shown in figures 5(
$d$
) and 5(
$e$
), where the case at the same
$Ca = 1.2$
from figures 5(
$b$
) and 5(
$c$
) are also superposed for comparison, together with the solution for steady flow. Here, we define the focusing times
$T$
and
$T_{ {st}}$
needed by the capsule centroid to reach the centreline (within a one fluid mesh corresponding to
$\sim 6\,\%$
of its radius to account for the oscillations in the capsule trajectory) under pulsatile and steady flows, respectively. Although the focusing time is decreased by almost
$50\,\%$
in prograde pulsation (
$\partial _z p_a^\ast = 2$
) compared with that in the steady flow, the time in antiphase pulsation is decreased by only
$1\,\%$
. Such small acceleration in antiphase pulsation comes from relatively small deformation in early periods (figure 5
$e$
). We now understand that fast axial focusing relies on the large membrane tension after flow onset, and our numerical results exhibit the even faster axial focusing due to the pulsation of the flow.
Figure 6(
$a$
) is the time history of the distance travelled along the flow direction (
$z$
-axis)
$r_z/D$
. The distance to complete the axial focusing (
$Ca = 1.2$
) under pulsatile flow increases compared with that in steady flow because the capsule speed along the flow direction increases by adding flow pulsation, where the circle dots represent the points when the capsule has completed the axial focusing. The capsule speed along the flow direction at
$Ca = 0.05$
, however, decreases with the pulsation of the flow. Figure 6(
$b$
) shows again the radial position of capsule centroids
$r_{{c}}/R$
as a function of
$z/D$
. The capsule trajectories obtained for
$Ca = 1.2$
remains almost the same, while the capsule trajectory for
$Ca = 0.05$
reaches equilibrium within a shorter travelled distance with pulsation. Following the classification by Vishwanathan & Juarez (Reference Vishwanathan and Juarez2021), our problem is oscillatory dominated, since the oscillation amplitude is one order of magnitude greater than the steady flow component (i.e.
$O(s\omega /\bar {u}^\prime ) \sim 10^1$
, where
$s$
is the centreline displacement amplitude and
$\bar {u}^\prime$
is the centreline velocity in a steady flow component). Notwithstanding this, the oscillatory motion was not enough to enhance the inertial focusing, in terms of channel lengths needed for the inertial focusing, because of the capsule deformations impeding the inertial focusing, consistent with previous numerical study (see figure 4
$a$
of Takeishi et al. Reference Takeishi, Yamashita, Omori, Yokoyama, Wada and Sugihara-Seki2022).
(
$a$
) Time history of the distance travelled along the flow direction (
$z$
-axis)
$z/D$
in the case shown in figure 5, where the circle dots represent the points when the capsule has completed the axial focusing. (
$b$
) Radial position of capsule centroids
$r_{{c}}/R$
as a function of
$z/D$
.
We now focus on axial focusing (i.e. cases of relatively high
$Ca$
) at finite
$Re$
. As discussed in figure 5(
$d$
), previous study showed that the speed of the axial focusing can be accelerated by the flow pulsation (Takeishi & Rosti Reference Takeishi and Rosti2023). An acceleration indicator of the axial focusing
$[1 - T/T_{{st}}]$
at
$Re = 10$
is summarised in figure 7, as a function of
$f^\ast$
(
$= f/\dot \gamma _{m}$
), where the results at
$Re = 0.2$
(Takeishi & Rosti Reference Takeishi and Rosti2023) are also supperposed. Although the initial radial position of the capsule
$r^\ast _{ {c0}}$
is slightly different between the two
$Re$
, the focusing time is commonly minimised at a specific frequency in both cases. Note that the values of the dimensional frequency depend on the estimation of
$G_s$
, which varies with the membrane constitutive laws and which is also sensitive to different experimental methodologies, e.g. atomic force microscopy, micropipette aspiration, etc. (Bao & Suresh Reference Bao and Suresh2003); the estimation of the dimensional frequency is therefore not trivial. We hereby conclude that capsules with large
$Ca$
exhibit axial focusing even at finite
$Re$
, and that their equilibrium radial positions are not altered by the flow pulsation.
We speculate that the optimal focusing frequency of
$f^\ast \approx 0.02$
, corresponding to dimensional frequency of
$f = 20$
Hz, is the membrane resonance frequency, given a reference radius of
$a_0 = 4$
$\mu\rm m$
and the surface shear elastic modulus of
$G_s = 4$
$\mu\rm N\,m^-{^1}$
(Takeishi et al. Reference Takeishi, Imai, Nakaaki, Yamaguchi and Ishikawa2014). However, there is currently no clear theoretical framework on the resonance frequency of capsule. To provide further insights into the state of resonance, we constructed a 2-D fluid membrane model (or hydrodynamic equations of bilayer membrane), obtained by Onsager’s variation principle, wherein the fluid membrane is assumed to be an almost planar bilayer membrane (Takeishi et al. Reference Takeishi, Santo, Yokoyama and Wada2024c
). Our numerical results showed that a membrane characteristic shift from an elastic-dominant to viscous-dominant state appears within the range
$40$
Hz
$\leqslant f \leqslant 400$
Hz, almost independently of surface tensions (figure 5
c of Takeishi et al. Reference Takeishi, Santo, Yokoyama and Wada2024c
). Since the resonance frequency can be formulated with intrinsic material (membrane) properties, it is expected that the value remains the same even under multi-capsule interactions. Indeed, we discovered that cross-over frequency of the storage and loss moduli in suspension of biconcave capsules modelling red blood cells (RBCs), whose inverse is defined as a relaxation time, is almost
$40$
Hz, regardless of the volume fraction of the capsules (figure 7
f of Takeishi et al. Reference Takeishi, Rosti, Yokoyama and Brandt2024b
). However note that the critical frequency was commonly estimated in terms of order of magnitude (
$O(f) = 10$
Hz) both in single- and multi-capsule dynamics as well as theoretical principles, since its exact estimation depends on
$G_s$
. Our recent numerical–experimental estimation strategy allows to quantify
$G_s$
of intact RBCs under dynamics and derive its value as
$\sim 0.5$
$\mu\rm N\,m^-{^1}$
(Takeishi et al. Reference Takeishi, Nishiyama, Nagaishi, Nashima and Sugihara-Seki2024a
), which is one order of magnitude smaller than that obtained by the stretch test (Takeishi et al. Reference Takeishi, Imai, Nakaaki, Yamaguchi and Ishikawa2014). Consequently, the dimensional optimal focusing frequency becomes
$O(f) = 10^2$
Hz, which is still in the range of the critical frequency estimated by the 2-D fluid membrane model (Takeishi et al. Reference Takeishi, Santo, Yokoyama and Wada2024c
). These results form a fundamental basis for further studies on resonance frequency of plasma membrane.
Acceleration indicator of the axial focusing
$[1 - T/T_{{st}}]$
as a function of the oscillatory frequency
$f^\ast$
for different
$Re$
(= 0.2 and 10).
$T$
and
$T_{{st}}$
are the elapsed time needed by the capsule centroid to reach the centreline under pulsatile and steady flows, respectively. The initial radial position of the capsule is set to be
$r^\ast _{ {c0}} = 0.55$
for
$Re = 0.2$
(see also figure 4
$a$
of Takeishi & Rosti Reference Takeishi and Rosti2023) and
$r^\ast _{ {c0}} = 0.45$
for
$Re = 10$
. The results are obtained with
$C{\kern-1pt}a$
= 1.2.
(
$a$
) Time histories of the radial position of the capsule centroid
$r_{{c}}/R$
at
$Re = 30$
and
$f^\ast = 0.02$
for different initial positions
$r^\ast _{ {c0}}$
(
$= 0.1$
and
$0.4$
), where insets represent snapshots of the lateral view of the deformed capsule at
$\dot {\gamma }_{m} t = (^*1)\ 60$
, (
$^\ast$
2)
$75$
and (
$^\ast$
3)
$90$
. Dashed lines are the curves
$r_{ {c}}^\ast = C_2 \exp {(-C_1 t^\ast )}$
, and the dash-dotted line denotes the equilibrium radial position of the capsule centroid. (
$b$
) Time histories of
$r_{{c}}/R$
for different
$Re$
, where dashed lines denote those in steady flow. (
$c$
) Time histories of
$r_{{c}}/R$
and
$\partial _z p^\ast$
at
$Re = 7$
(blue) and
$Re = 40$
(red), where the values are normalised by the amplitude
$\chi _{ {amp}}$
, and are shifted so that each baseline is the corresponding mean value
$\chi _{m}$
. Data are shown after
$\dot \gamma _{m} t \geqslant 350$
. (
$d$
) Peak frequency
$f^\ast _{{peak}}$
of the capsule centroid
$r_{{c}}/R$
. The solid line in panel (
$c$
) denotes the oscillatory frequency
$f^\ast = 0.02$
. The results are obtained with
$Ca = 0.05$
,
$R/a_0 = 2.5$
and
$r^\ast _{{c0}} = 0.4$
.
3.3. Effect of Reynolds number on capsule behaviour under pulsatile channel flow
We now focus on the inertial focusing of capsules at relatively small
$Ca$
, and, unless otherwise specified, we show the results obtained for
$Ca = 0.05$
. Figure 8(
$a$
) shows representative time history of the capsule centroid during inertial (or off-centre) focusing at
$Re = 30$
and
$f^\ast = 0.02$
for different initial position of the capsule
$r^\ast _{ {c0}}$
(
$= 0.1$
and
$0.4$
), where insets represent snapshots of the lateral view of the deformed capsule at various times
$\dot {\gamma }_{m} t$
(
$= 60$
,
$75$
and
$90$
), respectively. The results clearly show that the equilibrium radial position of the capsule is independent of its initial position
$r_c0$
(except when
$r_c0 = 0$
for which the capsule remains at centreline). Hereafter, each run case is started from a slightly off-centre radial position
$r^\ast _{ {c0}} = 0.4$
(
$R/a_0 = 2.5$
). For the trajectory at early times (
$\dot {\gamma }_{m} t \leqslant$
20), fitting by (3.1) still works. At quasi-steady state (
$\dot {\gamma }_{m} t \gt$
20), the capsule centroid fluctuates around an off-centre position
$r_{{c}}/R$
(
$\approx 0.3$
). Thus, the trajectory of the capsule during inertial focusing can be expressed as
\begin{equation} r_{\textit {c}}^\ast = \begin{cases} C_2 \exp {(-C_1 t^\ast )} & \text {for} \ t^\ast \leqslant t_{ {ax}}^\ast, \\ r_{ {e}}^\ast + \Delta r_{{osci}}^\ast & \text {for} \ t^\ast \gt t_{{ax}}^\ast, \end{cases} \end{equation}
where
$t_{{ax}}^\ast$
is the time period during axial focusing,
$r^\ast _{{e}}$
is the equilibrium radial position of the capsule centroid due to inertia and
$\Delta r_{{osci}}^\ast$
is a perturbation due to the oscillatory flow. Here, the equilibrium radial position is measured numerically by time averaging the radial position of the capsule centroid as
$r^\ast _{{e}} = \langle r_{{c}}^\ast \rangle = ( 1/\mathcal {T} ) \int _{t^\ast }^{t^\ast + \mathcal {T}} r_{{c}} (t^\prime ) {\rm d}t^\prime$
.
Figure 8(
$b$
) shows the time histories of the capsule centroid
$r_{ {c}}/R$
at
$f^\ast = 0.02$
for different
$Re$
, together with those with steady flow. We observe that the radial positions are greater than those at steady flow for all
$Re$
, due to the larger values achieved by the pressure gradient during the pulsation. However, the actual contribution of the oscillatory flow to the inertial focusing depends on
$Re$
. For instance, for
$Re \leqslant 7$
, the capsule exhibits axial focusing at steady flow, but a pulsatile channel flow allows the capsule to exhibit off-centre focusing. Therefore, the pulsation itself can impede the axial focusing.
Time average of (
$a$
) the radial position of the capsule centroid
$\langle r_{{c}} \rangle /R$
and (
$b$
) isotropic tensions
$\langle \tau_{{iso}} \rangle$
as a function of
$Re$
at
$f^\ast$
= 0.02, where the error bars represent the standard deviation during a period. The error bars in panel (
$b$
) are displayed only on one side of the mean value for major clarity. The results are obtained with
$Ca$
= 0.05,
$R/a_0$
= 2.5 and
$r^\ast _{ {c0}} = 0.4$
.
Figure 8(
$c$
) shows the waveforms of
$r_{{c}}/R$
at the end of the migration (
$\dot \gamma _{m}t \geqslant 350$
), where the instantaneous values are normalised by their respective amplitudes
$\chi _{ {amp}}$
and are shifted so that each baseline is the mean value
$\chi _{m}$
. Although the delay of
$r_{{c}}/R$
from the oscillatory pressure gradient
$\partial _z p^\ast$
tends to decrease as
$Re$
increases, the overall waveforms of
$r_{{c}}/R$
well follow that of
$\partial _z p^\ast$
, as shown in figure 5(
$b$
). To quantify the waveform of
$r_{{c}}/R$
and its correlation to
$\partial _z p^\ast$
, we extract the dominant (or peak) frequency
$f^\ast _{{peak}}$
of
$r_{{c}}/R$
with a discrete Fourier transform, whose principle and implementation are described by Takeishi et al. (Reference Takeishi, Rosti, Yokoyama and Brandt2024b
), and the result are shown as a function of
$Re$
in figures 8(
$d$
). In the cases of
$Re \leqslant 6$
, the capsule does not exhibit off-centre focusing and thus the plots are displayed for
$Re \geqslant 7$
only. The value of
$f^\ast _{{peak}}$
collapses on the frequency of
$\partial _z p^\ast$
with
$f^\ast = 0.02$
for
$Re \geqslant 7$
(figure 8
$d$
). The transition from the axial focusing to the off-centre focusing thus requires a synchronisation, induced by capsule deformability, between the capsule centroid and the background pressure gradient.
Figures 9(
$a$
) and 9(
$b$
) show the time average of the radial position or equilibrium position
$\langle r_{{c}} \rangle /R$
and the isotropic tension
$\langle \tau_{{iso}} \rangle$
, respectively, as a function of
$Re$
, where the error bars represent the standard deviation (
$SD$
) during a period. Overall, both these values nonlinearly increase with
$Re$
, with the mean values in the oscillating flows always greater than those in steady flows. The curves show steep increases for
$Re \leqslant 10$
, followed by a more moderate increase for
$Re \gt 10$
; these general tendencies are the same in steady or pulsatile flows. The effect of the flow pulsation is maximised at moderate
$Re$
(
$= 7$
), in which the axial focusing is impeded by the pulsatile flow (figure 9
$a$
). The results also show that small fluctuations of the capsule radial position (
$SD(r_{{c}}/R) \lt 10^{-2}$
) are accompanied by large fluctuations of the membrane tension (
$SD(\tau_{{iso}}) \gt 10^{-1}$
).
3.4. Effect of oscillatory frequency on capsule behaviour under pulsatile channel flow
Finally, we investigate the effect of the oscillatory frequency
$f^\ast$
on the equilibrium radial position
$\langle r_{{c}} \rangle /R$
at
$Re = 10$
, with the results summarised in figure 10(
$a$
), where those at steady flow are also displayed at the point
$f^\ast = 0$
. The results clearly suggest that there exists a specific frequency to maximise
$\langle r_{{c}} \rangle /R$
, independent of
$Re$
. Interestingly, such an effective frequency (
$f^\ast$
= 0.05) is close to or slightly larger than that maximising the axial focusing speed (see figure 7). Compared with steady flow, the equilibrium radial position
$\langle r_{{c}} \rangle /R$
at the effective frequency was enhanced by 640 % at
$Re$
= 7, 40 % at
$Re$
= 10, 13 % at
$Re$
= 20 and 7.6 % at
$Re$
= 30. The contribution of the oscillatory flow to the off-centre focusing becomes negligible for higher frequencies, in which the trajectory of the capsule centroid at the highest frequency considered (
$f^\ast = 5$
) collapses on that obtained with steady flow.
Time average of (
$a$
)
$\langle r_{{c}} \rangle /R$
and (
$b$
)
$\langle \tau_{ {iso}} \rangle /R$
as a function of
$f^\ast$
. The error bars in panel (
$b$
) are not displayed to increase clarity. All results are obtained with
$R/a_0 = 2.5$
,
$Re = 10$
and
$Ca = 0.05$
.
Figure 10(
$b$
) shows the time average of the isotropic tension
$\langle \tau_{ {iso}} \rangle$
as a function of
$f^\ast$
. The values of
$\langle \tau_{ {iso}} \rangle$
decrease as
$f^\ast$
increases because of the reduction of the shear stress when moving closer to the channel centreline (i.e. small
$\langle r_{{c}} \rangle /R$
). The results of large capsule deformation at relatively small frequencies are consistent with a previous numerical study by Matsunaga et al. (Reference Matsunaga, Imai, Yamaguchi and Ishikawa2015), who showed that at high frequency, a neo-Hookean spherical capsule undergoing oscillating sinusoidal shear flow cannot adapt to the flow changes, and only slightly deforms, consistent with predictions obtained by asymptotic theory (Barthès-Biesel & Rallison Reference Barthès-Biesel and Rallison1981; Barthès-Biesel & Sgaier Reference Barthès-Biesel and Sgaier1985). Thus, capsules at low frequencies exhibit an overshoot phenomenon, in which the peak deformation is larger than that its value in steady shear flow.
(
$a$
) Time history of
$r_{{c}}/R$
for different size ratios channel sizes
$R/a_0$
. (
$b$
and
$c$
) Time average of (
$b$
)
$\langle r_{{c}} \rangle /R$
and (
$c$
)
$\langle \tau_{ {iso}} \rangle /R$
as a function of
$R/a_0$
. The error bars in panels (
$b$
) and (
$c$
) are displayed only on one side of the mean value for clarity. (
$d$
) Time average of
$\langle r_{{c}} \rangle /R$
at
$Ca$
= 0.05 as a function of
$f^\ast$
for different
$R/a_0$
. All results are obtained with
$Re$
= 30,
$Ca$
= 0.05 and
$f^\ast$
= 0.02, and data at
$Re = 10$
are superposed in panel (
$d$
).
By increasing channel diameter
$D$
(=
$2R$
= 30
$\mu\rm m$
, 40
$\mu\rm m$
and 50
$\mu\rm m$
), we also investigate the effect of the size ratio
$R/a_0$
(
$= 3.75$
,
$5$
and
$6.25$
) on the equilibrium radial position
$\langle r_{{c}} \rangle /R$
. Figure 11(
$a$
) is the time history of
$r_{{c}}/R$
for different size ratios
$R/a_0$
at
$Re = 30$
and
$f^\ast = 0.02$
, where the trajectories obtained with the steady flow are also displayed. All run cases are started from
$r^\ast _{ {c0}} = 0.4$
. The equilibrium radial positions increase with
$R/a_0$
, while the contribution of oscillatory flow to
$\langle r_{{c}} \rangle /R$
becomes small as well as its fluctuation. This is quantified in figure 11(
$b$
), where
$\langle r_{{c}} \rangle /R$
is shown as a function of the size ratio
$R/a_0$
. Although the equilibrium radial position
$\langle r_{{c}} \rangle /R$
increases with
$R/a_0$
, indicating that dimensional equilibrium radial position
$\langle r_{{c}} \rangle$
also increases with
$R$
, the isotropic tension
$\langle \tau_{ {iso}} \rangle /G_s$
decreases, as shown in figure 11(
$c$
). This is because the distance from the capsule centroid to the wall (
$R - \langle r_{{c}} \rangle$
) increases with
$R$
, resulting in lower shear stress. Oscillatory-dependent off-centre focusing is summarised in figure 11(
$d$
), where the results are obtained with different channel size
$R/a_0$
and different
$Re$
(
$= 10$
and
$30$
). The result shows that oscillatory-dependent off-centre focusing is impeded as
$Re$
increases.
It is known that rigid particles align in an annulus at a radius of approximately
$0.6R$
for
$Re = D\overline {V}/\nu = O(1)$
(Segre & Silberberg Reference Segre and Silberberg1962; Matas et al. Reference Matas, Morris and Guazzelli2004, Reference Matas, Morris and Guazzelli2009), and shift to larger radius for larger
$Re$
(Matas et al. Reference Matas, Morris and Guazzelli2004, Reference Matas, Morris and Guazzelli2009), where
$\overline {V}$
is the average axial velocity (Matas et al. Reference Matas, Morris and Guazzelli2004). Our numerical results show that capsules with low deformability (
$Ca = 0.05$
) are still in
$\langle r_{{c}} \rangle /R \sim 0.5$
even for the largest channels (
$R/a_0$
= 6.25;
$R$
= 25
$\mu\rm m$
) and Reynolds number (
$Re = 30$
), both in the steady and pulsatile flows (figure 11
$b$
). Therefore, off-centre focusing is impeded even at such small particle deformability. This result is consistent with previous numerical study about a spherical hyperelastic particle in a circular channel with
$R/a_0$
= 5 under steady flow for 100
$\leqslant Re \leqslant$
400 and 0.00125
$\leqslant We \leqslant$
4 (Alghalibi et al. Reference Alghalibi, Rosti and Brandt2019). There, the authors showed that the particle radial position is
$\langle r_{{c}} \rangle /R \sim 0.5$
at the highest
$Re$
(
$= 400$
) and lowest
$We$
(
$= 0.00125$
). Our numerical results further show that the contribution of the flow pulsation to the off-centre focusing decreases as the channel size
$R/a_0$
increases (figures 11
$b$
and 11
$d$
) because of the low shear stress acting on the membranes (figure 11
$c$
). In other words, a large amplitude is required for oscillaton-induced off-centre focusing in high-
$Re$
and large channels.
The Poiseuille flow in a rigid circular pipe subject to the action of an oscillation pressure gradient was well described by Uchida (Reference Uchida1956); Womersley (Reference Womersley1955). At low frequency, oscillatory flow in the tube is better able to keep pace with the changing pressure. In the limit of zero frequency, the relation between flow and pressure is instantaneous as in a steady Poiseuille flow (see figure 16
$a$
in Appendix E). Thus, the particle (or capsule) is subject to shear stress, which results in its lateral movement due to the shear-induced lift force at low frequency. We confirmed it in figure 10(
$a$
), which is also consistent with the 2-D numerical results of a neutrally buoyant circular particle (Sun et al. Reference Sun, Yu and Shao2009). However, since the mechanism of axial focusing of capsules is primarily attributed to their deformability, the frequency-dependent axial focusing of a rigid (circular) particle remains unclear. At high frequency, however, oscillatory flow in a channel is less able to keep pace with the changing pressure, thus reaching less than the fully developed Poiseuille flow profile (almost flat velocity profile) at the peak of each cycle (see figure 16
$b$
in Appendix E). In the limit of infinite frequency, the velocity reached at the peak of each cycle is zero, that is, the fluid does not move at all. Thus, the particle does not experience shear stress and maintains its lateral position at high frequency, consistent with the 2-D numerical analysis by Sun et al. (Reference Sun, Yu and Shao2009).
Furthermore, we showed that
$\langle r_{{c}} \rangle /R$
increased with
$Re (\leqslant 30)$
, results consistent with those of rigid spherical particles in three-dimensional (3-D) steady pipe flows (Sun et al. Reference Sun, Yu and Shao2009). Such an increase in
$\langle r_{{c}} \rangle /R$
with
$Re$
can also be found for a rigid spherical particle on the square channel face, especially for
$Re \leqslant 100$
(Nakagawa et al. Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015), and also observed experimentally by Miura et al. (Reference Miura, Itano and Sugihara-Seki2014); Choi et al. (Reference Choi, Seo and Lee2011); Abbas et al. (Reference Abbas, Magaud, Gao and Geoffroy2014). It is also known that the channel face equilibrium positions decrease with
$Re$
, in particular for
$Re \gt 100$
, while the channel corner equilibrium positions continue to increase with
$Re$
(Nakagawa et al. Reference Nakagawa, Yabu, Otomo, Kase, Makino, Itano and Sugihara-Seki2015). Although our numerical results described in figure 11(
$d$
) suggest that a large amplitude is required for oscillaton-induced off-centre focusing at high
$Re$
, it remains an open question whether the off-centre focusing of capsules can indeed be enhanced by large-amplitude pulsatile flow and whether the optimal frequency remains consistent with the value measured in this study (
$f^\ast \approx 0.05$
).
Throughout our analyses, we have quantified the radial position of the capsule in a tube based on the empirical expression (3.3). We have provided insights about the coefficient
$C_1$
(
${\gt } 0$
) in
$r_{{c}}^\ast = C_2 \exp {(-C_1 t^\ast )}$
, which potentially scales the lift force and depends on shape, i.e. capillary number
$Ca$
and viscosity ratio
$\lambda$
.
4. Conclusion
We numerically investigated the lateral movement of spherical capsules in steady and pulsatile channel flows of a Newtonian fluid for a wide range of
$Re$
and oscillatory frequency
$f^\ast$
. The roles of size ratio
$R/a_0$
, and capillary number
$Ca$
on the lateral movement of the capsule have been evaluated and discussed. The first important question we focused on is whether a capsule lateral movement at finite
$Re$
in a pulsatile channel flow can be altered by its deformability. The second question is whether equilibrium radial positions or travelling time are controllable by oscillatory frequency.
Our numerical results showed that capsules with high
$Ca$
still exhibit axial focusing even at finite
$Re$
(e.g.
$Re = 10$
) and that their equilibrium radial positions cannot be altered by flow pulsation. However, the speed of axial focusing at such high
$Ca$
is substantially accelerated by making the driving pressure gradient oscillating in time. We also confirmed that there exists a most effective frequency (
$f^\ast \approx 0.02$
), which maximises the speed of axial focusing, and that it remains the same as that in almost inertialess condition. For relatively low
$Ca$
, however, the capsule exhibits off-centre focusing, resulting in an equilibrium radial position
$\langle r_{{c}} \rangle /R$
that depends on
$Re$
. There also exists a specific frequency to maximise
$\langle r_{{c}} \rangle /R$
, which is independent of
$Re$
. Interestingly, such effective frequency (
$f^\ast \approx 0.05$
) is close to that for axial focusing.
Frequency-dependent inertial focusing requires a synchronisation between the radial centroid position of the capsule and the background pressure gradient, resulting in periodic and large membrane tension, which impedes axial focusing. Such synchronisation abruptly appears at
$O(Re) = 10^0$
and shifts to an almost perfect syncronisation as
$Re$
increases. Thus, there is almost no contribution of flow pulsation to
$\langle r_{{c}} \rangle /R$
at relatively low
$Re$
(
$\leqslant 5$
) or large
$Re$
(
$\geqslant 30$
), while the contribution of the pulsation to
$\langle r_{{c}} \rangle /R$
is maximised at moderate
$Re$
(
$\approx 7$
), allowing the capsule to exhibit axial focusing in steady flow. For constant amplitude of oscillatory pressure gradient, oscillatory-dependent inertial focusing is impeded as
$Re$
and the channel diameter increase, and thus a relatively large oscillatory amplitude is required in such high
$Re$
and large channels. Throughout our analyses, we have quantified the radial position of the capsule in a tube based on the empirical expression. We hereby conclude that the knowledge obtained under inertialess conditions (Takeishi & Rosti Reference Takeishi and Rosti2023) has been extended to cases involving finite
$Re$
(
${\gt } 1$
) conditions.
Given that the speed of inertial focusing can be controlled by oscillatory frequency, the results obtained here can be used for label-free cell alignment/sorting/separation techniques, e.g. for circulating tumour cells in cancer patients or precious haematopoietic cells such as colony-forming cells.
Supplementary movie.
Supplementary movie are available at https://doi.org/10.1017/jfm.2025.184.
Funding.
N.T. acknowledges funding from Daicel Corporation. K.I. acknowledges the Japan Society for the Promotion of Science (JSPS) KAKENHI for Transformative Research Areas A (Grant No. 21H05309) and the Japan Science and Technology Agency (JST), FOREST (Grant No. JPMJFR212N). This research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding to M.E.R. from the Cabinet Office, Government of Japan. M.E.R. also acknowledges funding from the Japan Society for the Promotion of Science (JSPS), grants 24K17210 and 24K00810.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Numerical set-up and verification
To show that the channel length is adequate for studying the behaviour of a capsule that is subject to inertial flow, we have tested the channel length
$L$
(
$= 20a_0$
and
$40a_0$
) and investigated its effect on the radial positions of the capsule centroids. The time history of the radial position of the capsule centroid
$r_{{c}}$
is compared among these different channel lengths in figure 12, where the centroid position
$r_{{c}}$
is normalised by the channel radius
$R$
. The results obtained with the channel length
$L$
used in the main work (
$= 20 a_0$
) are consistent with those obtained with twice longer channel (
$L = 40 a_0$
).
Time history of the radial position
$r/R$
for different channel lengths
$L$
(
$= 20 a_0$
and
$40 a_0$
) and different
$Re$
(
$= 30$
and
$40$
). In all runs, the capsule is initially placed at
$r^\ast _{ {c0}} = 0.4$
. The results are obtained with
$R/a_0 = 2.5$
and
$Ca = 0.05$
.
Appendix B. Lift force on a capsule in a Poiseuille flow
We consider an object immersed in a Poisseulle flow, assuming that the flow is in the (steady) Stokes regime and that the object size is much smaller than the channel size. We also neglect any boundary effects acting on the object. Let
$y$
be the position relative to the channel centre. Due to the linearity of the Stokes equation, the object experiences a hydrodynamic resistance proportional to its moving velocity, given by
\begin{equation} f_1^L = - \xi _1 \dot {y}. \end{equation}
Note that the drag coefficient
$\xi _1 \gt 0$
is only determined by the viscosity and the shape (including the orientation) of the particle. We then consider the effects of the background Poiseuille flow. We have assumed that the channel size is much larger than the particle size, and hence the background flow to the particle is well approximated by a local shear flow with its local shear strength,
\begin{equation} \dot \gamma = -2 \frac {V_{ {max}}^\infty }{R^2} y. \end{equation}
In the presence of the background shear, the shear-induced lift force in general appears, and this is proportional to the shear strength (Kim & Karrila Reference Kim and Karrila2005),
\begin{equation} f_2^L = -\xi _2 \dot \gamma = 2 \xi _2 \frac {V_{ {max}}^\infty }{R^2} y, \end{equation}
where the coefficient
$\xi _2$
is again only determined by the viscosity and the shape. The force balance equation on the
$y$
-direction therefore reads
$f_1^L + f_2^L = 0$
. If we introduce a new shape-dependent coefficient,
$C_1$
, as
\begin{equation} C_1 = 2 \frac {\xi _2}{\xi _1} \frac {V_{ {max}}^\infty }{R^2}, \end{equation}
we obtain the evolution equation for the position
$y$
as
\begin{equation} \dot {y} = - C_1 y. \end{equation}
This equation is easily solved if
$C_1$
is constant and the result is the exponential accumulation to the channel centre, consistent with the numerical results.
Appendix C. Neo-Hookean spherical capsule
In consideration of previous works by e.g. Lefebvre & Barthès-Biesel (Reference Lefebvre and Barthès-Biesel2007) and Wang et al. (Reference Wang, Merlo, Dupont, Salsac and Barthès-Biesel2021), the trajectory of capsule centroids are compared among different types of membrane constitutive laws for a comprehensive understanding of capsule motion in a tube and to verify whether our empirical expression (3.1) works independently of the membrane constitutive law. Here, let us take the NK constitutive law, which is given by
\begin{equation} \frac {w_{\textit{NH}}}{G_s} = \frac {1}{2} \left ( I_1 - 1 + \frac {1}{I_2 + 1} \right ). \end{equation}
(
$a$
) Side views of the capsule during its axial focusing under steady flow for
$Ca = 0.05$
(top),
$Ca = 0.1$
(middle) and
$Ca = 0.2$
(bottom).The capsule is initially placed at
$r^\ast _{ {c0}} = 0.55$
. (
$b$
) Time histories of the radial position of these capsule centroids
$r_{{c}}/R$
. The dashed lines are the curves
$r_{{c}}^\ast = C_2 \exp {(-C_1 t^\ast )}$
. The result at the highest
$Ca$
(
$= 1.2$
) obtained with SK law is also superposed. The results are obtained with
$Re= 0.2$
,
$R/a_0 = 2.5$
and
$\lambda = 1$
.
Figure 13(
$a$
) shows side views of the capsule during its axial focusing at each time for different
$Ca$
(
$= 0.05$
,
$0.1$
and
$0.2$
). Other numerical settings (
$Re$
, initial position
$r^\ast _{ {c0}}$
and viscosity ratio
$\lambda$
) are the same as described in § 3.1. Even at relatively small
$Ca$
(
$= 0.2$
), the neo-Hookean capsule exhibits large elongation after flow onsets, resulting in fast axial focusing. The trajectory and fitting for it at each
$Ca$
are shown in figure 13(
$b$
), where the result at the highest
$Ca$
(
$= 1.2$
) obtained with SK law described in figure 2(
$b$
) is also superposed. The results suggest that (3.1) still works even for neo-Hookean spherical capsules, although the applicable ranges of
$Ca$
are relatively small compared with those described by the SK law.
Appendix D. Taylor parameter
The SK spherical capsule deformation is quantified by the Taylor parameter
$D_{12}$
, defined as
\begin{equation} D_{12} = \frac {| a_1 - a_2 |}{a_1 + a_2}, \end{equation}
where
$a_1$
and
$a_2$
are the lengths of the semi-major and semi-minor axes of the capsule, and are obtained from the eigenvalues of the inertia tensor of an equivalent ellipsoid approximating the deformed capsule (Ramanujan & Pozrikidis Reference Ramanujan and Pozrikidis1998).
Time histories of the Taylor parameter
$D_{12}$
for different
$Ca$
(
$= 0.05$
and
$1.2$
) at
$Re = 0.2$
. The results are obtained with
$f^\ast = 0.02$
and
$R/a_0 = 2.5$
.
Figure 14 shows the time history of
$D_{12}$
at
$Re = 10$
,
$R/a_0 = 2.5$
and
$f^\ast = 0.02$
. Different from what is observed for the isotropic tension shown in figure 5(
$c$
), the off-centred capsule exhibits large
$D_{12}$
, which well responds to the oscillatory pressure
$\partial _z p^\ast$
. Thus, the magnitude of
$D_{12}$
is strongly correlated with the capsule radial position (and the consequent shear gradient).
Figure 15(
$a$
–
$c$
) shows the time average of
$D_{12}$
. Overall, these results exhibit trends comparable to those of
$\langle \tau_{ {iso}} \rangle$
, previously shown in figures 9(
$b$
), 10(
$b$
) and 11(
$c$
). Despite the similarities, the axial-symmetric shaped capsule, typical of large
$Ca$
, exhibits small
$D_{12}$
(figure 15
$a$
), and the capsule membrane state in pipe flows cannot be easily estimated from the deformed shape. This is why we use the isotropic tension
$\tau_{ {iso}}$
as an indicator of membrane deformation.
Time average of
$\langle D_{12} \rangle$
as a function of (
$a$
)
$Re$
(obtained with
$Ca = 0.05$
and
$R/a_0 = 2.5$
), (
$b$
)
$f^\ast$
(obtained with
$Ca = 0.05$
and
$R/a_0 = 2.5$
) and (
$c$
)
$R/a_0$
(obtained with
$Re = 10$
and
$Ca = 0.05$
). The error bars in panels (
$a$
) and (
$c$
) are displayed only on one side of the mean value and are not displayed in panel (
$b$
) for clarity.
Appendix E. Oscillatory velocity profile in a rigid tube
Let us consider the Poiseuille flow in a rigid tube, with the radius of
$R$
, subject to the action of an oscillation pressure gradient, as described by Uchida (Reference Uchida1956); Womersley (Reference Womersley1955). The governing equation for oscillatory flow in cylindrical coordinates
$(r, \theta , z)$
is
\begin{equation} \frac {\partial ^2 v_z}{\partial r^2} + \frac {1}{r}\frac {\partial v_z}{\partial r} - \frac {1}{\nu }\frac {\partial v_z}{\partial t} = \frac {1}{\mu }\frac {\partial p}{\partial z}, \end{equation}
where the pressure gradient can be represented by a Fourier series
\begin{equation} \frac {\partial p}{\partial z} = \sum _{k = 0}^\infty c_k e^{i k \omega t}, \end{equation}
with
$c_0$
corresponding to the time average pressure gradient producing the Poiseuille profile. The solution is sought in terms of the Fourier series
\begin{equation} v_z (r, t) = \sum _{k = 0}^\infty \hat {w}_k e^{i k \omega t}. \end{equation}
Inserting (E3) in (E1), one gets
\begin{equation} \frac {d^2 \hat {w}_k}{{\rm d} r^2} + \frac {1}{r}\frac {{\rm d} \hat {w}_k}{{\rm d} r} - \frac {i \omega k}{\nu } \hat {w}_k = \frac {c_k}{\mu }, \end{equation}
where
$i^2 = -1$
, and defining the dimensionless variable
$\zeta = r/R$
, the non-homogeneous equation (E4) becomes
\begin{equation} \frac {d^2 \hat {w}_k}{{\rm d} \zeta ^2} + \frac {1}{\zeta }\frac {{\rm d} \hat {w}_k}{{\rm d} \zeta } - i k \alpha ^2 \hat {w}_k = \frac {c_k}{\mu }, \end{equation}
where
$\alpha$
is the Womersley number,
\begin{equation} \alpha = R\sqrt {\frac {\omega }{\nu }}. \end{equation}
With the boundary conditions
\begin{equation} \left . \hat {w}_k \right |_{r = R} = 0, \quad \left . \frac {{\rm d} \hat {w}_k}{{\rm d} r} \right |_{r = 0} = 0, \end{equation}
the particular solution of (E5) is easily found to be
$-({i \omega k}/{\nu }) \hat {w}_k = ({c_k}/{\mu })$
, and thus the global solution of the equation (for
$k \gt 0$
) is given by
\begin{equation} \hat {w}_k = \frac {ic_k R^2}{\mu k \alpha ^2} + \mathcal {C}_1 J_0 \left ( \alpha \sqrt {k} \zeta i^{3/2} \right ) + \mathcal {C}_2 Y_0 \left ( \alpha \sqrt {k} \zeta i^{3/2} \right ), \end{equation}
where
$\mathcal {C}_1, \mathcal {C}_2$
are arbitrary constants, and
$J_0, Y_0$
are Bessel functions of order zero of the first and second kind, respectively. Recall that
$i^{1/2} = e^{i\pi /4} = (1 + i)/\sqrt {2}$
.
The boundary conditions that the global solution must satisfy in a tube are the no-slip at the wall and the finite velocity along the axis of the tube, i.e.
\begin{equation} \hat {w}_k (R) = 0, \quad |\hat {w}_k (0)| \lt \infty , \end{equation}
which provide the required conditions to determine the constants
$\mathcal {C}_1, \mathcal {C}_2$
. It is known from the properties of
$Y_0 (\zeta )$
that
$Y_0 \to \infty$
as
$\zeta$
(or
$r$
) goes to
$0$
. Thus, the second boundary condition in (E9) leads to
$\mathcal {C}_2 = 0$
, and the first boundary condition then gives
\begin{equation} \mathcal {C}_1 = \frac {- i c_k R^2}{\mu \alpha _k^2 J_0 (\alpha _k i^{3/2})}, \end{equation}
where
$\alpha _k = \alpha \sqrt {k} = R\sqrt {k\omega /\nu }$
. With these values of
$\mathcal {C}_1, \mathcal {C}_2$
, the solution
$\hat {w}_k$
is finally
\begin{equation} \hat {w}_k = \frac {i c_k R^2}{\mu \alpha _k^2} \left ( 1 - \frac {J_0 (\alpha _k \frac {r}{R} i^{3/2})}{J_0 (\alpha _k i^{3/2})} \right ) \end{equation}
and the velocity profile
$v_z (r, t)$
is therefore
\begin{equation} v_z (r, t) = - \frac {c_0 R^2}{4 \mu } \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) + \frac {R^2}{\mu } \sum _{k = 1}^\infty \mathfrak {R} \left ( \frac {ic_k}{\alpha _k^2} \left [ 1 - \frac {J_0 (\alpha _k \frac {r}{R} i^{3/2})}{J_0 (\alpha _k i^{3/2})} \right ] e^{i k \omega t} \right ), \end{equation}
where
$\mathfrak {R}$
means the real part of a complex expression.
Let us consider oscillatory flow at low frequency, i.e. small
$\alpha$
. The Taylor series of
$J_0 (\zeta i^{3/2})$
is
\begin{equation} J_0 (\zeta i^{3/2}) = 1\! - \frac {(\zeta /2)^4}{(2!)^2} + \frac {(\zeta /2)^8}{(4!)^2} - \frac {(\zeta /2)^{12}}{(6!)^2} \cdots + i \!\left ( \!\frac {(\zeta /2)^4}{(1!)^2} + \frac {(\zeta /2)^6}{(3!)^2} - \frac {(\zeta /2)^{10}}{(5!)^2} - \cdots \! \right )\!, \end{equation}
and taking the dominant terms into account, one obtains
\begin{equation} \frac {J_0 (\alpha _k \frac {r}{R} i^{3/2})}{J_0 (\alpha _k i^{3/2})} = 1 - i \frac {\alpha _k^2}{4} \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) - \frac {\alpha _k^4}{64} \left ( 3 - 4 \left ( \frac {r}{R} \right )^2 + \left ( \frac {r}{R} \right )^4 \right ) \sin {\omega t} + O(\alpha _k^6). \end{equation}
The velocity is thus written as
\begin{align} v_z (r, t) &= - \frac {c_0 R^2}{4 \mu } \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) \end{align}
\begin{align} &- \frac {R^2}{\mu } \sum _{k = 1}^\infty \left \{ \frac {c_k}{4} \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) \cos {(k\omega t)} + \frac {c_k \alpha _k^2}{64} \left ( 3 - 4 \left ( \frac {r}{R} \right )^2 + \left ( \frac {r}{R} \right )^4 \right ) \sin {(k\omega t)} \right \}, \end{align}
and if we set
$V_i = - ({c_i R^2}/{4 \mu })$
(
$i = 0, 1$
) and
$\alpha _1 = \alpha$
, the first mode
$k = 1$
becomes
\begin{equation} \frac {v_z}{V_1} = \frac {V_0}{V_1} \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) + \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) \cos {(\omega t)} + \frac {\alpha ^2}{16} \left ( 3 - 4 \left ( \frac {r}{R} \right )^2 + \left ( \frac {r}{R} \right )^4 \right ) \sin {(\omega t)}. \end{equation}
Figure 16(
$a$
) shows the velocity profile for
$\alpha = 1$
at each phase angle
$\omega t (= 0, \pi /2, \pi , 3\pi /4$
) when the Poiseuille component of the velocity is neglected (i.e.
$V_0/V_1 = 0$
). Starting from the Poseuille flow at time
$\omega t = 0$
, at the phase of
$\omega t = \pi /2$
, the Poiseuille is still positive while the corresponding pressure gradient vanishes. The phase difference disappears at mid cycle (
$\omega t = \pi$
) when the Poiseuille flow is recovered. The profile reaches its peak form at the peak pressure gradient (
$\omega t = 0, \pi$
). The velocity profile for
$V_0/V_1 = 0.5$
, which is the same condition discussed in the main text, is also shown in figure 16(
$c$
) for completeness.
Oscillatory velocity profiles in a rigid tube at (
$a$
and
$c$
) low frequency (
$\alpha = 1$
) and (
$b$
and
$d$
) high frequency (
$\alpha = 10$
). The continuous Poiseuille component of the maximum velocity
$V_0$
is neglected in panels (
$a$
) and (
$b$
), while the finite value of
$V_0/V_1$
(
$= 0.5$
), which is the same condition used in the main text, is shown in panels (
$c$
) and (
$d$
). The lines represent the profiles at different phase angles (
$\omega t$
) within the oscillatory cycle, starting from
$\omega t = 0$
and increasing by steps of
$\pi /2$
.
At high frequency, the oscillatory flow in a rigid tube is less able to keep pace with the changing pressure, thus reaching less than the fully developed Poiseuille flow profile at the peak of each cycle. The parameter
$\alpha r/R$
takes large values and the axis (
$r = 0$
) is excluded from the analysis. For high values of its argument, the asymptotic development of
$J_0 (\zeta )$
is such that
\begin{equation} J_0 (\zeta ) = \sqrt {\frac {2}{\pi \zeta }} \cos {\left ( \zeta - \frac {\pi }{4} \right )} + O(|\zeta |^{-1}), \quad \text {with}\ | \text{arg}(\zeta )| \lt 2 \pi . \end{equation}
Using the relation
$i^{3/2} = e^{i 3\pi /4}$
and
$s = \alpha _k r/R$
, we can perform the following algebraic calculation:
\begin{align} J_0 (e^{i 3 \pi /4}s) &= e^{- 3 \pi /8} \sqrt {\frac {2}{\pi s}} \cos {\left ( e^{i 3 \pi /4}s - \frac {\pi }{4} \right )} \end{align}
\begin{align} &= e^{- i 3 \pi /8} \sqrt {\frac {2}{\pi s}} \cosh {\left ( \frac {s}{\sqrt {2}} + i \left ( \frac {s}{\sqrt {2}} + \frac {\pi }{4} \right ) \right )}, \\[6pt] \nonumber \end{align}
and neglecting the decaying exponential in cosh, since we deal with large values of the argument, we obtain
\begin{equation} J_0 (e^{i 3 \pi /4}s) = e^{- 3 \pi /8} \frac {1}{\sqrt {2 \pi s}} e^{\frac {s}{\sqrt {2}}} e^{i \left ( \frac {s}{\sqrt {2}} + \frac {\pi }{4} \right )}, \end{equation}
which leads us to find
\begin{equation} \frac {J_0 (\alpha _k \frac {r}{R} i^{3/2})}{J_0 (\alpha _k i^{3/2})} \approx \frac {1}{\sqrt {r/R}} e^{-(1+i) \frac {\alpha _k}{\sqrt {2}} \left (1 - \frac {r}{R} \right )}. \end{equation}
Finally, the first mode of the velocity profile yields
\begin{align} \frac {v_z}{V_1} &= \frac {V_0}{V_1} \left ( 1 - \left ( \frac {r}{R} \right )^2 \right ) \nonumber \\ &+ \frac {4}{\alpha ^2} \left [ \sin {(\omega t)} - \frac {1}{\sqrt {r/R}} e^{-\frac {\alpha }{\sqrt {2}} \left ( 1 - \frac {r}{R} \right )} \sin {\left ( \omega t - \frac {\alpha }{\sqrt {2}} \left (1 - \frac {r}{R} \right ) \right )} \right ] + O\left ( \frac {1}{\alpha ^4} \right ). \end{align}
Figure 16(
$b$
) shows the velocity profile for
$\alpha = 10$
at each phase angle, when the continuous Poiseuille component of the velocity is neglected. While the velocity is everywhere close to zero, the profile reaches its peak form at the phase of
$\omega t = \pi /2$
(and
$3\pi /2$
), i.e. the resulting flow is in complete phase shift (by
$\pi /2$
) with respect to the pressure gradient. The velocity profile for
$V_0/V_1 = 0.5$
is also shown in figure 16(
$d$
) for completeness. Thus, in this case, the continuous Poiseuille component is retained at each time.
Open access
