Inertial focusing of spherical capsule in pulsatile channel flows
Naoki Takeishi Kenta Ishimoto Naoto Yokoyama Marco Edoardo Rosti

Nodeless superconductivity and topological nodal states in molybdenum carbide
    Toni Shiroka
    
=============================================================================



§ ABSTRACT 
We present numerical analysis of the lateral movement of spherical capsule in the steady and pulsatile channel flow of a Newtonian fluid, for a wide range of oscillatory frequency. Each capsule membrane satisfying strain-hardening characteristic is simulated for different Reynolds numbers Re and capillary numbers Ca.

Our numerical results showed that capsules with high Ca exhibit axial focusing at finite Re similarly to the inertialess case.
We observe that the speed of the axial focusing can be substantially accelerated by making the driving pressure gradient oscillating in time.
We also confirm the existence of an optimal frequency which maximizes the speed of axial focusing, that remains the same found in the absence of inertia.
For relatively low Ca, on the other hand, the capsule exhibits off-centre focusing, resulting in various equilibrium radial positions depending on Re.
Our numerical results further clarifies the existence of a specific Re for which the effect of the flow pulsation to the equilibrium radial position is maximum.

The roles of channel size and viscosity ratio on the lateral movements of the capsule are also addressed.



  capsule,
  hyperelastic membrane,
  inertial focusing,
  off-centre focusing,
  pulsatile channel flow,
  computational biomechanics.




§ 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 <cit.>,
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 <cit.>. 

Although the techniques allow us to reduce the complexity and costs of clinical applications by using small amount of blood samples,
archetypal inertial focusing system requires steady laminar flow through long channel distances L_f, which can be estimated as L_f = π H/(Re_p f_l), where H is the dimension of the channel (or its hydraulic diameter) and f_l is a non-dimensional lift coefficient <cit.>.
So far, various kind of geometries have been proposed to achieve the required distance for inertial focusing in a compact space, e.g., sinusoidal, spiral, and hybrid channels <cit.>.
Without increasing Re_p, the recent experimental study by <cit.> achieved inertial focusing of 0.5-μm-size particle (Re_p ≈ 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,
utilizing oscillatory flow can be an alternative and practical strategy for inertial focusing in microfluidic devices.
Recently, <cit.> experimentally investigated the effects of the Womersley number (Wo) 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 to 10%) compared to what would be required in a steady flow <cit.>.


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 <cit.>,
numerical simulations <cit.>,
and experimental observations <cit.>,
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 <cit.>.


It is now well known that a deformable particle at low Re migrates toward the channel axis under steady laminar flow <cit.>.
Hereafter, we call this phenomenon as “axial focusing”.
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 <cit.>.
For finite Re (> 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 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.


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 <cit.> and numerical simulations <cit.>.


For instance, <cit.> 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 ≤ d_0/W ≤ 0.8,
using a rectangular channel with a high aspect ratio of W/H ≈ 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 <cit.>.
The experimental results can be qualitatively described using a spherical capsule <cit.> or droplet model <cit.>.
In more recent experiments by <cit.>, the authors investigated the effect of Re (1 < Re < 40) and capillary number Ca – ratio between the fluid viscous force and the membrane elastic force – (0.1 < Ca < 1) on the lateral equilibrium of bubbles in rectangular microchannels and different bubble-to-channel size ratios with 0.48 ≤ d_0/W ≤ 0.84.
The equilibrium position of such soft particles results from the competition between Re and Ca,
because high Re induce the 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 <cit.>.
Although a number of numerical analyses regarding inertial focusing have been reported in recent years mostly for spherical particles <cit.>,
the equilibrium positions of soft particles is still debated owing to the complexity of the phenomenon.
<cit.> showed that the equilibrium position in a cross section of rectangular microchannel with d_0/H = 0.2 shifts toward the wall as Re increases from 1 to 100.
<cit.> also performed numerical simulations of spherical capsules in a square channel for 0.1 ≤ d_0/H ≤ 0.4 and 5 ≤ Re ≤ 100 without viscosity contrast,
and showed that the equilibrium position was nearly independent of Re.
In a more recent numerical analysis by <cit.>,
simulations of a spherical hyperelastic particle in a circular channel with d_0/D = 0.2 were performed with 100 ≤ Re ≤ 400 and Weber number (We) with 0.125 ≤ We ≤ 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 toward the channel centre <cit.>.
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 μm in a cylindrical microchannel with D = 2R = 20–50 μm (i.e., R/a_0 = 2.5–6.25) for a wide range of oscillatory frequency.
Each capsule membrane, following the Skalak constitutive (SK) law <cit.>, 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 <cit.>.
The remainder of this paper is organised as follows.
Section 2 gives the problem statement and numerical methods,
Section 3 presents the numerical results for single spherical capsule. Finally, a summary of the main conclusions is reported in Section 4.
A description of numerical verifications is presented in the Appendix.




§ PROBLEM STATEMENT



 §.§ Flow and capsule models and setup

We consider the motion of an initially spherical capsule with diameter d_0 (= 2a_0 = 8 μm) flowing in a circular channel diameter D (= 2R = 20–50 μm),
with a resolution of 32 fluid lattices per capsule diameter d_0.
The channel length is set to be 20a_0, following previous numerical study <cit.>.
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 Fig. 7 in <cit.>),
we further assess the effect of this length on the lateral movement of a capsule in Appendix <ref> (figure <ref>a).

The capsule consists of a Newtonian fluid enclosed by a thin elastic membrane, sketched in figure <ref>.




The membrane is modeled as an isotropic and hyperelastic material following the SK law <cit.>, in which the strain energy w_SK and principal tensions in the membrane τ_1 and τ_2 (with τ_1 ≥τ_2) are given by

    w_SK/G_s = 1/4( I_1^2 + 2I_1 - 2I_2 + C I_2^2),

and

    τ_i/G_s = η_i/η_j[ η_i^2 - 1 + C ( η_i^2 η_j^2 - 1 ) ],   for (i, j) = (1, 2)  or (2, 1).

Here, w_SK is the strain energy density function,
G_s is the membrane shear elastic modulus,
C is a coefficient representing the area incompressibility,
I_1 (= η_1^2 + η_2^2 - 2) and I_2 (= η_1^2η_2^2 - 1) are the invariants of the strain tensor,
with η_1 and η_2 being the principal extension ratios.
In the SK law (<ref>),
the area dilation modulus is K_s = G_s (1 + 2C).
In this study, we set C = 10^2 <cit.>,
which describes an almost incompressible membrane.
Bending resistance is also considered <cit.>,
with a bending modulus k_b = 5.0 × 10^-19 J <cit.>.
These values have been shown to successfully reproduce the deformation of red blood cells in shear flow <cit.> and the thickness of cell-depleted peripheral layer in circular channels (see Figure A.1 in <cit.>).
Neglecting inertial effects on the membrane deformation,
the static local equilibrium equation of the membrane is given by

    ∇_s ·τ̌ + q̌ = 0̌,

where ∇_s (= ( Ǐ - ňň) ·∇) is the surface gradient operator,
ň is the unit normal outward vector in the deformed state,
q̌ is the load on the membrane,
and τ̌ is the in-plane elastic tension that is obtained using the SK law (equation <ref>).

The fluids are modeled with the incompressible Navier–Stokes equations for the fluid velocity v̌:

    ρ( ∂v̌/∂ t + v̌·∇v̌)
    	   = ∇·σ̌^f  + ρf̌, 
    ∇·v̌   = 0,

with 

    σ̌^f = -pǏ + μ( ∇v̌ + ∇v̌^T ),

where σ̌^f is the total stress tensor of the flow,
p is the pressure, 
ρ is the fluid density,
f̌ is the body force,
and μ is the viscosity of the liquid,
expressed using a volume fraction of the inner fluid I (0 ≤I≤ 1) as:

    μ = { 1 + ( λ - 1 ) I}μ_0,

where λ (= μ_1/μ_0) is the viscosity ratio, 
μ_0 is the external fluid viscosity, and μ_1 is the internal fluid viscosity.


The dynamic condition coupling the different phases requires the load q̌ to be equal to the traction jump ( σ̌^f_out - σ̌^f_in) across the membrane:

    q̌ = ( σ̌^f_out - σ̌^f_in) ·ň,

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 ∂ p_0/∂ z (= ∂_z p_0), 
which can be related to the maximum fluid velocity in the channel by ∂_z p_0 = -4 μ_0 V_max^∞/R^2.
The pulsation is given by a superimposed sinusoidal function,
such that the total pressure gradient is

    ∂_z p (t) = ∂_z p_0 + ∂_z p_a sin( 2 π f t ).



The problem is governed by six main non-dimensional numbers, including i) the Reynolds number Re and ii) the capillary number Ca defined as:

    Re = ρ D V_max^∞/μ_0, 
       Ca = μ_0 γ̇_m a_0/G_s = μ_0 V_max^∞/G_sa_0/4 R,

where V_max^∞ (= 2V_m^∞) is the maximum fluid velocity in the absence of any cells,
V_m^∞ is the mean fluid velocity,
and γ̇_m (= V_m^∞/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 iii) the viscosity ratio λ, iv) the size ratio R/a_0, v) the non-dimensional pulsation frequency f^∗ = f/γ̇_m, and vi) the non-dimensional pulsation amplitude ∂_z p_a^∗ = ∂_z p_a/∂_z p_0. Considered the focus of this study, we decide to primarily investigate the effect of Re, R/a_0, and f^∗.
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 λ, Reynolds number Re (> 1), and Ca.
We confirmed that the flow at Re = 0.2 well approximates an almost inertialess flow for single- <cit.> and multi-cellular flow <cit.>.
Unless otherwise specified, we show the results obtained with ∂_z p_a^∗ = 2 and λ = 1.



 §.§ Numerical simulation

The governing equations for the fluid are discretised by the LBM based on the D3Q19 model <cit.>.
We track the Lagrangian points of the membrane material points x̌_m (X̌_m,t) over time,
where X̌_m is a material point on the membrane in the reference state.
Based on the virtual work principle,
the above strong-form equation (<ref>) can be rewritten in weak form as 

    ∫_S û̌·q̌ dS = ∫_S ϵ̌̂̌ : τ̌ dS,

where S is the surface area of the capsule membrane, and û̌ and ϵ̌̂̌ = ( ∇_s û̌ + ∇_s û̌^T )/2 are the virtual displacement and virtual strain, respectively.
The FEM is used to solve equation (<ref>) and obtain the load q̌ acting on the membrane <cit.>.
The velocity at the membrane node is obtained by interpolating the velocities at the fluid node using the immersed boundary method <cit.>.
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 <cit.> and front-tracking method <cit.> 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 <cit.>:
in our simulation, the relative volume error is always maintained lower than 1.0 × 10^-3%,
as tested and validated in our previous study of cell flow in circular channels <cit.>.
All procedures were fully implemented on a GPU to accelerate the numerical simulation.
More precise explanations for numerical simulations including membrane mechanics are provided in our previous works <cit.>.

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 <cit.>.



 §.§ Analysis of capsule deformation

Later, we investigate the in-plane principal tension T_i (with T_1 ≥ T_2) and the isotropic tension T_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 T_iso = (T_1 + T_2)/2 for the deformed capsule.
The averaged value of T_iso is then calculated as

    ⟨ T_iso⟩ = 1/ST∫_T∫_S T_iso (x̌_m, t) dS dt,

where T is the period of the capsule motion.
Hereafter, ⟨·⟩ 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 γ̇_m t = 200,
and we start accumulating the statistics from γ̇_m t ≥ 400 to fully cancel the influence of the initial conditions. 




§ RESULTS



 §.§ Axial focusing of the capsule under steady channel flow (Re < 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 <ref>(a) shows side views of the capsule during its axial focusing in channel of size R/a_0 = 2.5 for different Ca (= 0.05, 0.2, and 1.2).
The capsule, initially placed at r^∗_c = r_c/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 process is commonly observed for every Ca.
The time history of the radial position of the capsule centroid r_c is shown in figure <ref>(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:

    r_c^∗ = βexp(-α t^∗),

where t^∗ (= γ̇_mt) is the non-dimensional time,
and α (> 0) and β are two coefficients that can be found by a least-squares fitting to the plot.
Fitting are performed using data between the initial (r_c/R = 0.55) and final state (Δ x_LBM/R ≤ 0.01 for R/a_0 = 2.5),
defined as the time when the capsule is within one mesh size (Δ x_LBM) from the channel axis.


Performing time differentiation of equation (<ref>),
the non-dimensional velocity of the capsule centroid ṙ_c^∗ can be estimated as:

    ṙ_c^∗ = -α r_c^∗.

This linear relation (<ref>) 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 α and the lift force on the capsule are provided in Appendix <ref>.

Figure <ref>(a) shows the coefficient α as a function of Ca.
As expected from figure <ref>(b),
the value of α increases with Ca.
Since the capsule deformability is also affected by the viscosity ratio λ,
its influence on α is also investigated in figure <ref>(b).
At a fixed Ca (= 1.2),
the value of α decreases with λ.



To further proof that α 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_0/R.
Note that a case with larger channel for constant Re denotes smaller V_max^∞, resulting in smaller G_s (i.e., softer capsule) for constant Ca.
Figure <ref>(a) is one of the additional runs at Ca = 0.2,
where the capsule is initially placed at r_0/R = 0.75.
Figure <ref>(b) is the time history of the radial position of the capsule centroid r_c for different initial positions r_0/R.
We observe that the exponential fitting is still applicable for these runs, with the coefficient α reported in figure <ref>(c).
These results provide a confirmation that α is indeed independent of the initial radial position r_0/R.
Furthermore, the fitting provided in equation (<ref>) 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 <ref> (see also figure <ref>).





 §.§ 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 <ref>(a),
which are obtained with f^∗ = 0.02 and Re = 10.
The simulations are started from a off-centre radial position r_0/R = 0.4.

At the end of the migration, the least deformable capsule (Ca = 0.05) exhibits an ellipsoidal shape with an off-centred position (figure <ref>a, left), while the most deformable one (Ca = 1.2) exhibits the typical parachute shape at the channel centreline (figure <ref>a, right).
Detailed trajectories of these capsule centroids r_c/R are shown in figure <ref>(b),
where the non-dimensional oscillatory pressure gradient ∂_z p^∗ (t^∗) (= 1 + 2 sin(2 π f^∗ t^∗)) is also displayed.
The least deformable capsule (Ca = 0.05) fluctuates around the off-centre position r_c/R (≈ 0.2),
and the waveform of r_c/R lags behind ∂_z p^∗ (t^∗).
The capsule with large Ca (= 1.2), on the other hand, immediately exhibits axial focusing, reaching the centerline within one flow period (figure <ref>b).
Therefore, axial and off-centre focusing strongly depend on Ca.

Figure <ref>(c) is the time history of the isotropic tension T_12.
The major waveforms of T_iso are synchronised with ∂_z p^∗ 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 <ref> (see figure <ref>).

To clarify whether fast axial focusing depends on the phase of oscillation or not,
an antiphase pulsation (i.e., ∂_z p_a^∗ = -2) is given by ∂_z p^∗ (t^∗) = 1 - 2 sin(2 π f^∗ t^∗).
Time histories of the capsule centroid r_c/R and membrane tension T_iso under such condition are shown in figures <ref>(d) and <ref>(e),
where the case at the same Ca = 1.2 from figures <ref>(b) and <ref>(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 ∼ 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 almost by 50% in prograde pulsation (∂_z p_a^∗ = 2) compared to that in the steady flow,
the time in antiphase pulsation is decreased only by 1%.
Such small acceleration in antiphase pulsation comes from relatively small deformation in early periods (figure <ref>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 <ref>(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 comparing to 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, on the other hand, decreases with the pulsation of the flow.
Figure <ref>(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 traveled distance with pulsation.




Following the classification by <cit.>,
our problem is oscillatory dominated,
since the oscillation amplitude is one order of magnitude greater than the steady flow component (i.e., O(sω/u̅^') ∼ 10^1,
where s is the centreline displacement amplitude and u̅^' 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,

consistently with previous numerical study (see figure 4a in <cit.>).



We now focus on axial focusing (i.e., cases of relatively high Ca) at finite Re.
As discussed in figure <ref>(d),
previous study showed that the speed of the axial focusing can be accelerated by the flow pulsation <cit.>. 

An acceleration indicator of the axial focusing [1 - T/T_st] at Re = 10 is summarised in figure <ref>,
as a function of f^∗ (= f/γ̇_m),
where the results at Re = 0.2 <cit.> are also supperposed.
Although the initial radial position of the capsule r_0/R 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. <cit.>;
the estimation of the dimensional frequency is therefore not trivial and left as for future investigations.
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.





 §.§ 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 <ref>(a) shows representative time history of the capsule centroid during inertial (or off-centre) focusing at Re = 30 and f^∗ = 0.02 for different initial position of the capsule r_0/R (= 0.1 and 0.4),
where insets represent snapshots of the lateral view of deformed capsule at various time γ̇_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_0
(except when r_0 = 0 for which the capsule remains at centreline).
Hereafter, each run case is started from a slightly off-centre radial position r_0/R = 0.4 (R/a_0 = 2.5).
For the trajectory at early times (γ̇_m t ≤ 20),
fitting by equation (<ref>) still works.
At quasi-steady state (γ̇_m t > 20),
the capsule centroid fluctuates around an off-centre position r_c/R (≈ 0.3).
Thus, the trajectory of the capsule during inertial focusing can be expressed as

    r_c^∗ = 
      βexp(-α t^∗)   for t^∗≤ t_ax^∗
    
      r_e^∗ + Δ r_osci^∗   for t^∗ > t_ax^∗
      ,

where t_ax^∗ is the time period during axial focusing,
r^∗_e is the equilibrium radial position of the capsule centroid due to inertia,
and Δ  r_osci^∗ 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^∗_e = ⟨ r_c^∗⟩ = ( 1/T) ∫_t^∗^t^∗ + T r_c (t^') dt^'.



Figure <ref>(b) shows the time histories of the capsule centroid r_c/R at f^∗ = 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 ≤ 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.

Figure <ref>(c) shows the waveforms of r_c/R at the end of the migration (γ̇_m≥ 350),
where the instantaneous values are  normalised by their respective amplitudes χ_amp and and shifted so that each baseline is the mean value χ_m.
Although the delay of r_c/R from the oscillatory pressure gradient ∂_z p^∗ tends to decrease as Re increases,
the overall waveforms of r_c/R well follow that of ∂_z p^∗, as shown in figure <ref>(b).
To quantify the waveform of r_c/R and its correlation to ∂_z p^∗,
we extract the dominant (or peak) frequency f^∗_peak of r_c/R with a discrete Fourier transform,
whose principle and implementation are described in <cit.>,
and the result are shown as a function of Re in figures <ref>(d).
In the cases of Re ≤ 6, the capsule does not exhibit off-centre focusing, and thus the plots are displayed for Re ≥ 7 only.
The value of f^∗_peak collapses on the frequency of ∂_z p^∗ with f^∗ = 0.02 for Re ≥ 7 (figure <ref>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 <ref>(a) and <ref>(b) show the time average of the radial position or equilibrium position ⟨ r_c⟩/R and the isotropic tension ⟨ T_iso⟩, 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 ≤ 10, 
followed by a more moderate increases for Re > 10;
these general tendency 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 <ref>a).

The results also show that small fluctuations of the capsule radial position (SD(r_c/R) < 10^-2) are accompanied by large fluctuations of the membrane tension (SD(T_iso) > 10^-1).




 §.§ Effect of oscillatory frequency on capsule behaviour under pulsatile channel flow

Finally, we investigate the effect of the oscillatory frequency f^∗ on the equilibrium radial position ⟨ r_c⟩/R at Re = 10,
with the results summarised in Figure <ref>(a),
where those at steady flow are also displayed at the point f^∗ = 0.
The results clearly suggest that there exists a specific frequency to maximise ⟨ r_c⟩/R, independently of Re.
Interestingly, such effective frequency (f^∗ = 0.05) are close to or slightly larger than those maximising the axial focusing speed (see figure <ref>).
Comparing to steady flow,
the equilibrium radial position ⟨ r_c⟩/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^∗ = 5) collapses on that obtained with steady flow.


Figure <ref>(b) shows the time average of the isotropic tension ⟨ T_iso⟩ as a function of f^∗.
The values of ⟨ T_iso⟩ decrease as f^∗ increases because of the reduction of the shear stress when moving closer to the channel centreline (i.e., small ⟨ r_c⟩/R).
The results of large capsule deformation at relatively small frequencies are consistent with a previous numerical study by <cit.>, 
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, consistently with  predictions obtained by asymptotic theory <cit.>. Thus, capsules at low frequencies exhibit an overshoot phenomenon,
in which the peak deformation is larger than that its value in steady shear flow. 

By increasing channel diameter D (= 2R = 30 μm, 40 μm, and 50 μm),
we also investigate the effect of the size ratio R/a_0 (= 3.75, 5, and 6.25) on the equilibrium radial position ⟨ r_c⟩/R.
Figure <ref>(a) is the time history of r_c/R for different size ratios R/a_0 at Re = 30, and f^∗ = 0.02,
where the trajectories obtained with the steady flow are also displayed.
All run cases are started from r_0/R = 0.4.
The equilibrium radial positions increase with R/a_0,
while the contribution of oscillatory flow to ⟨ r_c⟩/R becomes small as well as its fluctuation.
This is quantified in figure <ref>(b),
where ⟨ r_c⟩/R is shown as a function of the size ratio R/a_0.
Although the equilibrium radial position ⟨ r_c⟩/R increases with R/a_0,
indicating that dimensional equilibrium radial position ⟨ r_c⟩ also increases with R,
the isotropic tension ⟨ T_iso⟩/G_s decreases as shown in figure <ref>(c).
This is because the distance from the capsule centroid to the wall (R - ⟨ r_c⟩) increases with R, resulting in lower shear stress.
Oscillatory-dependent off-centre focusing is summarised in figure <ref>(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 about 0.6R for Re = DV/ν = O(1) <cit.>,
and shift to larger radius for larger Re <cit.>,
where V is the average axial velocity <cit.>.
Our numerical results show that capsules with low deformability (Ca = 0.05) are still in ⟨ r_c⟩/R ∼ 0.5 even for the largest channels (R/a_0 = 6.25; R = 25 μm) and Reynolds numbers (Re = 30), both in the steady and pulsatile flows (figure <ref>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 ≤ Re ≤ 400 and 0.00125 ≤ We ≤ 4 <cit.>.
There, the authors showed that the particle radial position is ⟨ r_c⟩/R ∼ 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 <ref>b and <ref>d) because of the low shear stress acting on the membranes (figure <ref>c).
In other word, a large amplitude is required for oscillaton-induced off-centre focusing in 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 (<ref>).










We have provided insights about the coefficient α (> 0) in r_c^∗ = βexp(-α t^∗), which potentially scales the lift force and depends on shape, i.e., capillary number Ca and viscosity ratio λ.

 



§ 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^∗.
The roles of size ratio R/a_0, viscosity ratio λ, 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 traveling 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^∗≈ 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, on the other hand, the capsule exhibits off-centre focusing, resulting in an equilibrium radial position ⟨ r_c⟩/R which depends on Re.
There also exists a specific frequency to maximise ⟨ r_c⟩/R, which is independent of Re.
Interestingly, such effective frequency (f^∗≈  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 appear at O(Re) = 10^0,
and shifts to an almost perfect syncrohisation as Re increases.
Thus, there is almost no contribution of flow pulsation to ⟨ r_c⟩/R at relatively low Re (≤ 5) or large Re (≥ 30),
while the contribution of the pulsation to ⟨ r_c⟩/R is maximised at moderate Re (≈ 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 channel diameter increase,
and thus relatively large oscillatory amplitude is required in such high Re and large channels.

Given that the speed of inertial focusing can be controlled by oscillatory frequency,
the results obtained here can be utilized for label-free cell alignment/sorting/separation techniques, e.g., for circulating tumor cells in cancer patients or precious hematopoietic cells such as colony-forming cells.



§ ACKNOWLEDGEMENTS

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.
The presented study was partially funded by 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).




§ CONFLICTS OF INTEREST

The authors report no conflict of interest.






§ NUMERICAL SETUP 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 between these different channel lengths in figure <ref>, 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).








§ 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


    f_1^L = - ξ_1 ẏ.


Note that the drag coefficient ξ_1 > 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,

    γ̇= -2 V_max^∞/R^2 y.

In the presence of the background shear,
the shear-induced lift force in general appears,
and this is proportional to the shear strength <cit.>,

    f_2^L = -ξ_2 γ̇= 2 ξ_2 V_max^∞/R^2 y,

where the coefficient ξ_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, α, as

    α = 2 ξ_2/ξ_1V_max^∞/R^2,

we obtain the evolution equation for the position y as

    ẏ = - α y.

This equation is easily solved if α is constant and the result is the exponential accumulation to the channel centre, consistent with the numerical results. 



§ NEO-HOOKEAN SPHERICAL CAPSULE

The NK constitutive law is given by

    w_NH/G_s = 1/2( I_1 - 1 + 1/I_2 + 1).



Figure <ref>(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_0/R, and viscosity ratio λ) are the same as described in <ref>.
Even at relatively small Ca (= 0.2),
the NH-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 <ref>(b),
where the result at the highest Ca (= 1.2) obtained with SK law described in figure <ref>(b) is also superposed.
The results suggest that equation (<ref>) still works even for NH-spherical capsules, although the applicable ranges of Ca are relatively small compared to those described by the SK law.



§ TAYLOR PARAMETER

The SK-spherical capsule deformation is quantified by the Taylor parameter D_12, defined as

    D_12 = | a_1 - a_2 |/a_1 + a_2,

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 <cit.>.

Figure <ref> shows the time history of D_12 at Re = 10, R/a_0 = 2.5, and f^∗ = 0.02.
Differently from what observed for the isotropic tension shown in figure <ref>(c),
the off-centred capsule exhibits large D_12,
which well responds to the oscillatory pressure ∂_z p^∗. 

Thus, the magnitude of D_12 is strongly correlated with the capsule radial position (and the consequent shear gradient).



Figures <ref>(a–c) are the time average of D_12.
Overall, these results exhibit trends comparable to those of ⟨ T_iso⟩, previously shown in figures <ref>(b), <ref>(b), and <ref>(c).
Despite the similarities,
the axial-symmetric shaped capsule,
typical of large Ca,
exhibits small D_12 (figure <ref>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 T_iso as an indicator of membrane deformation.


jfm



