3.1. Fluid pumping

The motion of the artificial cilium reaches a steady state after a few periods. When $\theta _0 = 0$, the cilium beating patterns are periodic and symmetric over one period (figure 2a). The particle follows a symmetric ‘figure-8’ trajectory. As a result, the time-averaged disturbance flow field shows no net fluid pumping along the $y$ direction (figure 2b). However, when $\theta _0 < 0$, we observe a distinct power ($F_0 < 0$) and recovery stroke ($F_0 > 0$). A typical example is shown in figure 2(c). The filament is stretched out during the power stroke (red). As $\boldsymbol {F}_0$ reverses, the filament is bent with large deformation and moves towards the $+y$ direction during the recovery stroke (blue). The particle traces out an asymmetric trajectory and the disturbance flow field shows a clear net flux along the $-y$ direction (figure 2d). The centre of mass of the filament also traces out an asymmetric trajectory, similar to the natural cilium (Brumley et al. Reference Brumley, Wan, Polin and Goldstein2014).

Figure 2.

Time lapse of the artificial cilium over one forcing period in the steady state and the time-averaged disturbance flow field with $\eta =3.5$ and $\beta =0.08$ for (a,b) $\theta _0=0$ (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.436), and (c,d) $\theta _0=-1.16$ (supplementary movie 2). The force amplitude $A=65$. Time runs from red to blue. The cyan lines in (a,c) indicate the initial orientations. The dark dashed lines show the particle trajectories and the green dashed lines show the centre-of-mass position of the filament. The colour bars in (c,d) show the magnitude of the flow speed and arrows indicate flow directions.

The pumping performance can be characterized by the flux of the disturbance flow field obtained by integrating the Blake tensor over the $x$–$z$ plane perpendicular to the pumping direction (Liron Reference Liron1978). The resulting instantaneous flux due to a point force of unit strength along the $+y$ direction located at a distance $h$ from the no-slip wall is given by $h/{\rm \pi} \mu$. The net flux $Q$ of the composite cilium consists of the flux generated by the filament $Q_{{f}}$ and the flux generated by the particle $Q_{{p}}$

(3.1)\begin{equation} Q = Q_{{f}} + Q_{{p}} = \frac{1}{{\rm \pi} \mu} \left\langle\int_{{-}L/2}^{L/2}x(s)f_y(s)\,{\rm d}s\right\rangle + 6b \langle \chi x_{{p}}\boldsymbol{v}_{{p}}\boldsymbol{\cdot} \hat{\boldsymbol{y}}\rangle, \end{equation}

where $\langle \rangle$ denotes time average over a period in steady state and the distance of the particle from the wall $x_{{p}} = (x+b\cos \theta )|_{s={L}/{2}}$. The particle flux $Q_{{p}} \approx 6b\langle x_p \boldsymbol {v}_p\boldsymbol{\cdot} \hat {\boldsymbol {y}} \rangle = 6b\langle x_p \dot {y}_p \rangle = 6bS/\tau$, where $S$ is the area swept by the particle (Osterman & Vilfan Reference Osterman and Vilfan2011). This states that $Q_{{p}}$ is proportional to the area enclosed by the non-reciprocal trajectory of the particle.

Figure 3(a) shows the effect of $\theta _0$ with fixed force amplitude $A=65$. As $|\theta _0|$ increases (tilting towards the $-y$ direction), both $|Q_{{f}}|$ and $|Q_{{p}}|$ increase. Since the area swept by the particle is much larger than the area swept by the centre-of-mass position of the filament (figure 2c), $|Q_{{p}}|$ is always larger than $|Q_{{f}}|$. Surprisingly, for sufficiently large $|\theta _0|$, $|Q_{{f}}|$ and $|Q_{{p}}|$ drop abruptly to smaller values. The time lapse of the filament deformation shows that the abrupt change of $Q_{{p}}$ is accompanied by a reversal of filament bending direction. As shown by the two insets in figure 3(a), the filament is bent downwards during the recovery stroke for $\theta _0 = -1.12$ and upwards for $\theta _0 = -0.96$. This difference in the filament deformation leads to a difference in the particle trajectories, and therefore an abrupt change in $Q_{{p}}$. A similar discontinuity is observed when varying the force amplitude $A$ with fixed $\theta _0=-0.9$ (figure 3b). As $A$ is increased, the flux increases. When $A$ becomes sufficiently large, both $|Q_{{f}}|$ and $|Q_{{p}}|$ jump to larger values with a reversal of the bending direction.

Figure 3.

(a,b) Time-averaged flux $|Q_{{f}}|$ (left axis, dark triangles), $|Q_{{p}}|$ (left axis, dark squares) and the total flux $|Q|$ (left axis, dark circles) as a function of (a) the initial tilt angle $\theta _0$ and (b) the force amplitude $A$ with $\beta =0.08$ and $\eta =3.5$. The flux is scaled by $B/\mu L$. In (a) $A = 65$, and in (b) $\theta _0 = -0.9$. The deviations $|\Delta \theta |=|\bar {\theta }_{{m}}-\theta _0|$ are shown on the right axes by the blue diamonds. The insets show the time lapse of the artificial cilium for different values of (a) $\theta _0$ (supplementary movies 3 and 4) and (b) $A$ corresponding to the filled red and green squares. The initial orientations are marked by the cyan lines. Panel (c) shows $|Q_{{p}}|$ as a function of $\theta _0$ and $A$.

The filament bending direction is likely determined by deviation of cilium orientation at the start of the recovery stroke from the initial orientation, which is also the natural orientation at rest with no external force. We compute the deviation as $\Delta \theta = \bar {\theta }_{{m}}(t=n)-\theta _0$, where the average orientation in steady state $\bar {\theta }_{{m}} = \int _{-1/2}^{1/2}\theta (s)\,{\rm d}s$ and $n$ is an integer. The deviation $|\Delta \theta |$ is largest for the symmetric case with $\theta _0 = 0$ (figure 2a). Figure 3(a,b) shows that $|\Delta \theta |$ decreases as the cilium is more tilted (with fixed $A$) or as the driving force $A$ decreases (with fixed $\theta _0$). The filament changes from downward bending to upward bending as $|\Delta \theta |$ exceeds a critical value around 0.15 in figure 3(a) and 0.18 in figure 3(b). As an example, the case of $\theta _0 = -1.12$ (downward bending) in figure 3(a) shows smaller deviation than the case of $\theta _0 = -0.96$ (upward bending). A similar observation can be made by comparing the two insets in figure 3(b) with $A = 50.9$ and $A = 69.8$.

We then compute the density plot of $Q_{{p}}$ as functions of $\theta _0$ and $A$. The values of $A$ and $\theta _0$ are limited to avoid the cilium touching the no-slip wall. As shown in figure 3(c), two regions with large negative flux are identified at large $A$ and $\theta _0$ corresponding to upward and downward bending. The sharp boundary separating these two regions spans a wide range for $\theta _0 \lessapprox -0.8$ and $A \gtrapprox 55$.

3.2. Linear stability analysis

When a filament is under compression at its tip along the tangential direction, the filament buckles as the compression exceeds a critical value (Landau & Lifshitz Reference Landau and Lifshitz1986). To investigate whether the tangent component of $\boldsymbol {F}_0$ is large enough to induce buckling in our system, especially at large values of $|\theta _0|$, we perform a linear stability analysis on the composite cilium with a tangential compression force $\varGamma$ acting on the particle.

Consider a small deformation from an initially straight filament with $\theta _0 = 0$, then $x\sim s$, $\boldsymbol {p} \sim (1, y_x)$, and $\boldsymbol {p}^{\perp } \sim (-y_x, 1)$. We use non-dimensional equations and ignore the non-local interaction. The linearized tension equation is $T_{ss} = 0$ with $T_{s} = 0$ at $s=-1/2$ and $T+3/2c\beta \chi T_{s} = -\varGamma$ at $s=1/2$. This leads to $T = -\varGamma$. We linearize (2.1) as

(3.2)\begin{equation} -\alpha \varGamma y_{xx}-\alpha y_{xxxx} = y_{t}, \end{equation}

where $\alpha = \eta ^{-4}(c+2)$. The boundary conditions are

(3.3)\begin{equation} y = 0,\quad y_{x} = 0, \end{equation}

at $s=-1/2$, and

(3.4a,b)\begin{equation} y_{xx} + \beta^3\eta^{4} y_{xt} = 0,\quad y_{xxx} - 3/4\eta^{4}\beta\chi (y_{t}+\beta y_{xt}) = 0, \end{equation}

at $s=1/2$. We consider perturbations of the form $y(x,t) = \phi (x)e^{\lambda t}$, and solve the resulting eigenvalue problem numerically using centred finite differences in the bulk and sided differences at the boundaries. Figure 4 shows that the real part of the largest eigenvalue ${Re}(\lambda )$ decreases first and then increases as $\varGamma$ is increased. The system becomes unstable if ${Re}(\lambda ) > 0$. For $\beta = 0$, the critical value $\varGamma ^{\ast } \approx 37.6$, which agrees with the result given in De Canio et al. (Reference De Canio, Lauga and Goldstein2017). The critical value becomes smaller as $\beta$ is increased. For $\beta = 0.08$, $\varGamma ^{\ast } \approx 22.6$.

Figure 4.

Linear stability analysis. (a) The real part of the largest eigenvalue ${Re}(\gamma )$ as a function of $\varGamma$ for different values of $\beta$. (b) In steady state, the $y$-component position of the particle (dark curve, left axis), actuation force $F_0$ (red solid curve, right axis) and filament tangent force $F_{{tang}}$ at $s=1/2$ (red dashed curve, right axis) as a function of time with $A = 65$, $\eta = 3.5$, $\beta = 0.08$ and $\theta _0 = -1.12$. The green dashed lines indicate $F_0^{\ast }$, the value of $F_0$ at the maximum of $|F_{{tang}}|$ around $t \approx 7.1$.

We observe a signature of buckling instability in our simulations. As shown in figure 5(b), computed with $\theta _0=-1.12$ (corresponding to the left inset in figure 3a), although $F_0$ increases rapidly first during the recovery stroke starting at $t=7.0$, the particle remains almost fixed for a period of time with little change in its $y$-component position. Meanwhile, the magnitude of the filament tangent force at $s=1/2$, expressed as $F_{{tang}}|_{s=1/2} = [(-\boldsymbol {r}_{sss}+T\boldsymbol {r}_s)\boldsymbol{\cdot} \boldsymbol {p}]_{s=1/2}$, increases and reaches a maximum at $t\approx 7.1$. The filament then buckles and the particle moves towards the $+y$ direction with $F_{{tang}}$ released. Therefore, we use $F^{\ast } = F_0(t\approx 7.1)$ to estimate the compression force acting on the composite cilia as $\varGamma = |F^{\ast }\sin (\theta _0)| \approx 37.9$. We apply the same approximation to the case of $\theta _0 = -0.96$ (the right inset in figure 3a) and obtain $\varGamma \approx 41.0$. Both estimates are larger than $\varGamma ^{\ast }$, indicating that buckling instability indeed occurs at large $|\theta _0|$ and the abrupt change in the net flux is caused by a reversal of the filament buckling direction.

Figure 5.

Results of the segmental model with $\beta =0.08$, $\eta = 3.5$, and $\gamma _j=1/3$. (a,b) The time-averaged flux generated by the particle $|Q_{p}|$ (dark squares, left axis) as a function of (a) $\theta _0$ and (b) $A$. In (a) $A = 18.0$, and in (b) $\theta _0 = -1.0$. The flux is scale by $K/\mu$. Insets show the time lapse of the artificial cilia for different values of (a) $\theta _0$ (supplementary movies 5 and 6) and (b) $A$ corresponding to the red and green squares. Cyan lines indicate initial orientations. The deviations $|\Delta \theta |=|\bar {\theta }_{{m}}-\theta _0|$ are shown on the right axes by the blue circles. (c) Relative deflection $\theta _3-\theta _2$ under a constant driving force $F_0 = 18.0$ as a function of time $t$ for different values of $\delta \theta$ with $\theta _0 = -1.0$.

3.3. Segmental model

In the segmental model, the generated flux is also dominated by the contribution from the particle, especially at large $|\theta _0|$. Similar to the elastic model, abrupt change in $Q_{{p}}$ is also observed when varying $\theta _0$. Figure 5(a) shows that $|Q_{{p}}|$ first increases monotonically as $|\theta _0|$ increases from 0. For sufficiently large $|\theta _0|$, $|Q_{{p}}|$ suddenly jumps to a smaller value, along with a reversal of the buckling direction during the recovery stroke: segment 2 and segment 3 buckle upwards ($\theta _2 > \theta _3$, see inset with $\theta _0 = -1.0$) when $\theta _0 \gtrsim -1.0$ and downwards ($\theta _2 < \theta _3$, see inset with $\theta _0 = -1.2$) when $\theta _0 \lesssim -1.0$. An abrupt change in $Q_{{p}}$ is also found when the force amplitude $A$ is changed (figure 5b). As $A$ exceeds a critical value around 22.5, segments 2 and 3 switch from a downward buckling to an upward buckling with a significant increase in $|Q_{{p}}|$. Comparing the two cases shown by the insets of figure 5(b) with $A=13.6$ and $20.6$, $|Q_{{p}}|$ is tripled with an apparent increase in the enclosed area by the particle trajectory. The average orientation of the cilium is $\bar {\theta }_{{m}} = (\theta _1+\theta _2+\theta _3)/3$, and the deviation from the natural orientation $\Delta \theta = \bar {\theta }_{{m}}-\theta _0$. The transition from downward buckling to upward buckling is accompanied with a pronounced increase in $|\Delta \theta |$.

Finally, to verify that the reversal of the buckling direction is indeed caused by the deviation from the natural orientation, we perform the following numerical experiment. We apply a constant driving force $F_0=18.0$ along the $+y$ direction and evolve the system for a time duration of $0.5$. The deviation from a straight line is varied: we set $\theta _1(t=0) = \theta _0$, $\theta _2(t=0) = \theta _0-\delta \theta$ and $\theta _3 (t=0)= \theta _0-2\delta \theta$, where $\delta \theta$ is the magnitude of the deviation. The buckling direction is characterized by the relative deflection between segments 2 and 3. As shown in figure 5(c), segments 2 and 3 buckle downwards with $\theta _3-\theta _2 > 0$ for $\delta \theta \lessapprox 0.13$. An abrupt change occurs when $\delta \theta \gtrapprox 0.13$, segments 2 and 3 buckle upwards with $\theta _3-\theta _2 < 0$ for most of the time and reach equilibrium positions as $t\to 0.5$.