A bar-joint model based on the corrected resistive force theory for artificial flagellated micro-swimmers propelled by acoustic waves

An artificial micro-swimmer normally refers to an object synthesized at the microscopic scale [1, 2] and capable of interacting with a fluidic environment to achieve locomotion. While natural micro-swimmers, i.e. microorganisms, have been studied since the invention of Leeuwenhoek's microscope in the 17th century, artificial micro-swimmers were not thought possible until Richard Feynman's famous lecture [3] on nanotechnology in 1959. Since then, people have fantasized about tiny untethered robots [4] to kill bacteria and implement targeted therapy in vivo. For example, a swarm of micro-swimmers actuated by acoustic power can swim along the intestinal tract to kill the cancer cells by generating heat under an alternated magnetic field (i.e. magnetic thermotherapy), as illustrated in figure 1. Nevertheless, these tasks are still very challenging, and most of the achievements are limited to laboratory studies [5, 6] where the micro-robots appear in diverse structures, such as bubbles [7, 8], nanorods [9, 10], helical filaments [11, 12], and flagella [13, 14]. In this work, we focus on artificial flagellated micro-swimmers (AFMSs).

Zoom In Zoom Out Reset image size

An AFMS usually consists of a relatively rigid head and a flexible tail, as shown in the inset of figure 1, where the tail is customarily called a flagellum in microbiology [15]. We focus on the flagellated structure, especially the sperm-like one, because it is easy to fabricate through 3D printing [14] and is suitable for medical applications due to the potential of cargo transportation by its head [16]. Various actuation strategies have been proposed for artificial micro-swimmers, which can be classified into two categories [17]: the self-propulsion approaches, e.g. the catalytic reactions [18], and external actuation approaches, such as optical [19], thermal [20], magnetic [21], and acoustic actuation [22]. For AFMSs and considering medical applications, the acoustic and magnetic actuations are the most suitable because they are tissue- and circulation-friendly and have been employed in diagnosis and therapies [23]. Here we choose acoustic excitations for propulsion and navigation and also consider magnetic power to assist turning.

It is noted that the inertia of flagella and surrounding fluids can be ignored because the Reynolds number is sufficiently low at the microscale and when the dynamic response is away from resonance [24]. In a low Reynolds number (LRN) regime, the Navier–Stokes equation can be recast into the linear Stokes equation [25] and the linearity of the Stokes equation leads to the so-called Purcell's scallop theorem [26] which forbids objects with a reciprocal motion from moving. Instead, a traveling wave in a flagellum can generate net propulsion [27]; hence, whipping the flagellum of a micro-swimmer can propel itself forward. One may note that the strong acoustic streaming around the sharp edges is a well-established mechanism for propelling acoustically actuated swimmers [28, 29]. To simplify the swimming problem, as described by Liu and Ruan [24], the effect of actuation (e.g. an acoustic power) is simplified to actuation forces and torques on the head of a flagellated micro-swimmer, which leads to flagellum wiggling and propulsion. However, the expressions for fluid forces on a solid body are normally so complicated that numerical simulations are necessary to predict the locomotion [30]. Fortunately, if a body is slender (i.e. the ratio of the characteristic width to length is less than 0.1 [31]), an asymptotic expression in terms of the slenderness can be obtained, known as the slender body theory (SBT) [32–36]. However, the SBT is still too complex to be applied in most cases, because the interactions between different parts of the body have to be considered based on the SBT. If omitting these interactions and only considering the drags, i.e. each segment of a slender flagellum independently experiences resistance to its motion, the SBT is simplified into the resistive-force theory (RFT) [35, 37, 38].

Based on the RFT and the balance of the drag and elastic forces, Machin [39] proposed the first-order governing equation for a uniform flagellum. His work was followed by Wiggins and Goldstein [40] who established the so-called hyper-diffusion equation to describe an AFMS. Then, the propulsive force generated by whipping an elastic filament was experimentally verified based on the head movement [41]. The partial differential equations (PDEs) of a variable cross-section flagellum were proposed by Singh and Yadava [42]. Kaynak et al [14] first reported the fabrication of such sperm-like micro-swimmers and observed their swimming under the actuation of acoustic waves. Liu and Ruan [24] investigated the corresponding swimming mechanisms by assuming that the swing of the head arising from the acoustic actuation leads to the whipping of the flagellum and that the interaction between the flagellum and its surrounding fluid can generate propulsion. They employed the Galerkin method [43] to solve the governing equations and exhibited the results consistent with the experimental ones. It should be noted that the models of the AFMS are usually one-dimensional (1D) [6], which confines the swimmer's motion to a straight line. Considering the need for manipulation or steering in applications, such a 1D model is insufficient. Establishing at least a two-dimensional (2D) AFMS model to address plane motion is therefore indispensable.

One may derive the 2D governing PDEs of a flagellum via Hamilton's principle [44] if the dynamic profile of the flagellum is a known function [45]. Unfortunately, in swimming the wiggling profile of a flagellum is not ad hoc. In addition, a flagellum becomes curved in turning [46], bringing about geometric nonlinearly, and the influence of the rigid-body motion further increases the nonlinearity of problem. Therefore, a discrete flagellum model was developed to solve the 2D swimming problem. Inspired by Purcell's three-link swimmer [26], a flexible flagellum has been modeled by using multiple rigid bars linked by flexible joints (aka. the bar-joint models) [47], which intrinsically fulfills the inextensible constraints of a flexible tail [48, 49]. Moreau et al [50] compared the classical continuum model with the bar-joint model and demonstrated that the latter could achieve better numerical performance. However, these models were not verified by simulations or experiments.

It is generally nontrivial to verify and correct micro-swimmer models experimentally because geometries and material properties are often uncertain (or undermined in literature). In addition, the wiggling profile of a flagellum is difficult to capture under a microscope. Hence, numerical simulations dealing with fluid-structure interactions (FSIs), which are usually based on the finite element method (FEM), were more desirable to verify a simplified micro-swimmer model. Rorai et al [31] studied a sperm-like micro-swimmer in a 3D infinite domain by FEM, and demonstrated an excellent agreement between the SBT and the FSI simulations in terms of the propulsive matrix coefficients. Nevertheless, the dynamic motion of the flagellum was specified a priori, and the swimmer for each time was fixed, i.e. the elasticity of the flagellum and the resultant rigid-body motion were not investigated. Curatolo and Teresi [51] employed the automatic remeshing technique [52, 53] to couple the structure's rigid-body motion with FSI simulations for a fish-like swimmer. However, they did not compare their simulation results with those of RFT or SBT.

Undoubtedly, differences exist between the results of FSI simulations and the RFT or SBT models due to various reasons. For example, the RFT and SBT can hardly address fluid problems around a flagellum tip or a head-flagellum joint [31]. The fluid force exerted on a head is also difficult to estimate because the head is normally neither a sphere nor a slender body. Most importantly, the RFT and SBT are initially established for a slender cylindrical body, but the flagellum in 2D simulations (considering that 3D simulations are computationally too demanding) approximates a slender body with rectangular cross-sections, which is often the case in experiments (e.g. [14]). In addition, the most prevailing AFMS fabrication is based on layer-by-layer photocuring (i.e. 3D printing) [54], which leads to flagella with rectangular cross-sections and hence, the need for modifying SBT and RFT for them.

The SBT for arbitrary cross-section has been developed by adding a dimensionless coefficient tensor K [55], depending only on the cross-sectional geometry, to the fluid velocity field. Borker and Koch [56] extended this method to particle dynamics in shear flows and determined the magnitudes of components of K through FSI simulations. However, such a geometry-dependent K is insufficient to predict the locomotion of an AFMS because of the aforementioned complexities. In particular, when dealing with AFMS turning, the corrected SBT is still very difficult to solve [57]. Hence, we argue that a correction to RFT, considering both the effect of non-circular cross-section and the match of the linear terminal speed with the corresponding FSI simulation, may be a better strategy of approximation.

In the subsequent sections, we will follow the above thought to develop the governing equations of an AFMS based on a bar-joint model and a corrected resistive-force theory (CRFT) for AFMS with rectangular cross-sections. In section 2, we will derive the expressions of CRFT and the governing equations of a bar-joint model. In section 3, we will first state the swimming problem and corresponding FSI simulations, then examples of convergence analyses of FSI simulation will be presented. After that, the correction factor K is determined in terms of different geometrical, materials, and acoustic-actuation parameters. Comparisons of swimmer's trajectories between the CRFT-based bar-joint model and the FSI simulation will then be presented as a verification of our theory, and a demonstration of the acoustic or magneto-acoustic steering strategy for the AFMS will be shown as the examples of applications. Finally, we will conclude the results of this research in section 4.

The correction of SBT for the case of non-circular cross-section has been developed in [55]. Here we briefly introduce the concept and then simplify it to a more convenient form, i.e. the CRFT. In a LRN regime, the Navier–Stokes equations of fluid reduce to the linear Stokes equations [25] given by:

$$\begin{align} - \nabla p + \mu {\nabla ^2}{\mathbf{u}} = 0{\text{, }}\nabla \cdot {\mathbf{u}} = 0,\end{align} \tag{ 1a,b }$$

where u and p are the velocity and pressure fields of a fluid, respectively, and µ the dynamic viscosity. Considering a slender body of unit length, let ${\mathbf{S}}\left( s \right) = {\left[ {{S_\parallel }\left( s \right),{S_ \bot }\left( s \right)} \right]^{\text{T}}}$ be the distribution of Stokeslets [58], where the superscript T denotes transpose and 0 ⩽ s ⩽ 1 is the arclength along the body centerline. Applying the linear superposition method (equation (1a,b ) are linear) to solve the fluid field and pressure distribution, the total fluid force F _f acting on a flagellum of unit length can be obtained as:

$$\begin{align}\,{{\mathbf{F}}_f} = - 8\pi \mu \int_0^1 {{\mathbf{S}}\left( s \right){\text{d}}s} .\end{align} \tag{ 2 }$$

Let P(s) and κ(s) be the coordinates and the half-width of the body at the position s along the body centerline, where κ(s) also indicates how slender the body is. If κ is sufficiently small and the cross-section is circular, one can obtain the classical expression of S based on the SBT [35]:

$$\begin{align} {\mathbf{S}}\left( s \right) &= \frac{{2{\mathbf{I}} - {\hat{\mathbf{x}}}\left( s \right){\hat{\mathbf{x}}}\left( s \right)}}{{4\ln \kappa \left( s \right)}} \cdot \left\{ {{\mathbf{u}}\left( s \right) - {\mathbf{V}}\left( s \right) + {S_ \bot }\left( s \right){\boldsymbol{\hat{\mathbf{y}}}}\left( s \right)} \right.\nonumber \\ &\quad+ \ln \left[ {4s\left( {1\! -\! s} \right)} \right]{\mathbf{S}}\left( s \right) \!+\! {S_\parallel }\left( s \right)\left\{ {\ln \left[ {4s\left( {1 \!- \! s} \right)} \right] - 2} \right\}{\hat{\mathbf{x}}}\left( s \right) \nonumber \\ &\quad + \int_{ - s}^{1 - s} \left[ \frac{{{\mathbf{S}}\left( {s + \delta } \right)}}{{\left| {{{\mathbf{P}}_{s\delta }}} \right|}} - \frac{{{\mathbf{S}}\left( s \right)}}{{\left| \delta \right|}} + \frac{{{{\mathbf{P}}_{s\delta }}{{\mathbf{P}}_{s\delta }} \cdot {\mathbf{S}}\left( {s + \delta } \right)}}{{{{\left| {{{\mathbf{P}}_{s\delta }}} \right|}^3}}}\right.\nonumber\\ &\qquad\qquad\quad\left.\left. -\,\, \frac{{{S_\parallel }\left( s \right){\delta ^2}{\hat{\mathbf{x}}}\left( s \right)}}{{{{\left| \delta \right|}^3}}} \right]{\text{d}}\delta \right\} \end{align} \tag{ 3 }$$

where ${\hat{\mathbf{x}}}\left( s \right)$ and ${\hat{\mathbf{y}}}\left( s \right)$ are orthogonal unit vectors that are tangential and normal to the body centerline, respectively, I the identity matrix, P _sδ = P(s) – P(s + δ) the distance vector, and V(s) the local translational velocity of the body.

To deal with arbitrary flagellum cross-sections or 2D problems (i.e. with rectangular cross-sections), we follow the method proposed in [55] to introduce a dimensionless coefficient tensor K to the expression of fluid velocity field. Based on the SBT for circular flagellum cross-section, the relative fluid velocity field can be expressed as: ${\mathbf{U}} = {\mathbf{C}}\cdot{\mathbf{S}}$, where U(s) = u(s) – V(s) represents the relative velocity of the fluid with respect to the flagellum part at point s. For non-circular ones, it is assumed that the relative fluid velocity for the non-circular case is: ${\mathbf{U}} = \left( {{\mathbf{C}} + {\mathbf{K}}} \right)\cdot{\mathbf{S}}$, where C is the resistive tensor, and the Stokeslets S is still based on the expression for a slender body with a circular cross-section. Note that the revised expression of U can still meet the zero-divergence condition for fluid velocity (i.e. equation (1b)) as long as K is a symmetric tensor with zero divergence.

In [55], K depends only on the shape of the cross-section, which can be directly estimated for some special shapes without performing FSI simulations. Nevertheless, in most cases, FSI simulations are required to determine K because the dynamics of a solid body (e.g. the wiggling profile of a flagellum), which affects U and K, cannot be determined a priori, as reported by Borker and Koch [56]. Hence, we assume that K = [K_ij ] (i, j = 1, 2 for 2D problems) is the only correction tensor to determine and that it must be determined based on the comparison between FSI simulation and the SBT model. In the following, we reduce the corrected SBT to the corrected RFT by neglecting interactions between different parts of a flagellum. This is generally an acceptable simplification in the study of micro-swimmers [35] considering the limited accuracy in experimental observation. In RFT, the fluid forces F _f acting on a slender body is assumed to be proportional to the local velocity of the body. We shall then include K in the coefficients of proportionality to establish the CRFT.

For shorthand, we recast equation (3) as:

$$\begin{align}{\mathbf{S}}\left( s \right) = {\boldsymbol{\mit\kappa }}\left( s \right)\cdot\left\{ {{\mathbf{U}}\left( s \right) + \psi \left[ {{\mathbf{S}}\left( s \right)} \right]} \right\},\end{align} \tag{ 4 }$$

where ${\boldsymbol{\mit\kappa }}\left( s \right) = {{\left[ {2{\mathbf{I}} - {\hat{\mathbf{x}}}\left( s \right){\hat{\mathbf{x}}}\left( s \right)} \right]} \mathord{\left/ \right. } {\left[ {4\ln \kappa \left( s \right)} \right]}}$, and ψ[S(s)] is the functional of Stokeslets S(s) and ψ(0) = 0. Taking the correction term ${\mathbf{K}}\cdot{\mathbf{S}}$ into account, equation (4) is rewritten as:

$$\begin{align}{\mathbf{S}}\left( s \right) = {\boldsymbol{\mit\kappa }}\left( s \right)\cdot\left\{ {{\mathbf{U}}\left( s \right) + {\mathbf{K}}\cdot{\mathbf{S}}\left( s \right) + \psi \left[ {{\mathbf{S}}\left( s \right)} \right]} \right\}.\end{align} \tag{ 5 }$$

The approximation of equation (5) can be made based on the scenario of a straight flagellum of a uniform cross-section. In this scenario, U and κ are constants, and the direction vectors ${\hat{\mathbf{x}}}$ and ${\hat{\mathbf{y}}}$ are independent of position s, which greatly simplifies the solution. An iterative procedure can then be adopted. First, the zeroth-order approximation S⁽⁰⁾(s) can be determined by applying S = 0 to the right-hand side of equation (5) and is expressed as:

$$\begin{align}{{\mathbf{S}}^{\left( 0 \right)}}\left( s \right) = {\boldsymbol{\mit\kappa }}\cdot{\mathbf{U}}.\end{align} \tag{ 6 }$$

Substituting equation (6) into the right-hand side of equation (5) again results in the first-order approximation S⁽¹⁾(s). The second-order approximation S⁽²⁾(s) can be obtained similarly and is normally accurate enough for the RFT [35]. Eventually, we can obtain the total fluid force F _f on a uniform straight flagellum of unit length with a unidirectional motion by substituting S⁽²⁾(s) into equation (2).

Denoted by ${U_\parallel }\left( s \right)$ and ${U_ \bot }\left( s \right)$ the components of U(s) tangential and normal to the local centerline of a flagellum, respectively, the corresponding fluid forces in the RFT are formulated as:

$$\begin{align}{F_f}_\parallel = {c_\parallel }{U_\parallel }{\text{, }}{F_f}_ \bot = {c_ \bot }{U_ \bot }.\end{align} \tag{ 7a,b }$$

Based on the second-order approximation S⁽²⁾(s) and equation (2), the two coefficients of resistive forces in equation (7a, b ) in the CRFT can be derived (cf appendix A) as:

$$\begin{align}{c_\parallel } &= 2\pi \mu \left[ \frac{1}{{\ln \left( {1/\kappa } \right)}} + \frac{{\frac{3}{2} - \ln 2 - \frac{K}{4}}}{{{{\ln }^2}\left( {1/\kappa } \right)}} \right.\nonumber\\ &\left.\quad+\,\, \frac{{{{\left( {\frac{3}{2} - \ln 2 - \frac{K}{4}} \right)}^2} + 1 - \frac{{{\pi ^2}}}{{12}}}}{{{{\ln }^3}\left( {1/\kappa } \right)}} \right],\end{align} \tag{ 8 }$$

and

$$\begin{align}{c_ \bot } &= 4\pi \mu \left[ \frac{1}{{\ln \left( {1/\kappa } \right)}} + \frac{{\frac{1}{2} - \ln 2 + \frac{K}{2}}}{{{{\ln }^2}\left( {1/\kappa } \right)}} \right.\nonumber\\ &\left.\quad+\,\, \frac{{{{\left( {\frac{1}{2} - \ln 2 + \frac{K}{2}} \right)}^2} + 1 - \frac{{{\pi ^2}}}{{12}}}}{{{{\ln }^3}\left( {1/\kappa } \right)}} \right],\end{align} \tag{ 9 }$$

where the correction factor K = K₁₁ = −K₂₂ due to the property of zero divergence of the tensor K. Note that when K = 0, equations (8) and (9) reduce to the RFT expressions in [35]. The proposed CRFT is then applicable to formulate the governing equations of an AFMS based on the bar-joint model.

Figure 2 demonstrates the bar-joint model of the AFMS. The flagellum is simplified into rigid bars with a constant length 2 l. They are indexed (denoted by i) from 1 to n and connected with linear torsion springs [47]. The ith spring joint connects the bars i and (i+ 1) and the first (i = 1) bar is assumed to be rigidly connected to the head. The motion of the ith bar and the head (i = 0) can be represented by the vector of time-dependent variables X _i (t) = (x_i, y_i, θ_i )^T, where x_i and y_i are coordinates of the midpoint of the ith bar (or the center of the head), and θ_i the angle of rotation with respect to the x-axis (anticlockwise rotation is taken as positive). Correspondingly, the velocity vector of the head (i = 0) and bars (i > 0) can be represented by ${{\dot{\mathbf{X}}}_i} = {\left( {{{\mathop {{\text{ }}x}\limits^ \cdot }_i},{{\mathop {{\text{ }}y}\limits^ \cdot }_i},{{\mathop \theta \limits^ \cdot }_i}} \right)^{\text{T}}}$, where the overhead dot denotes time derivative.

Zoom In Zoom Out Reset image size

Based on the CRFT, one can easily obtain the three equilibrium equations of external forces for the bar-joint model, given by:

$$\begin{align}\left[ {\begin{array}{*{20}{c}} {{{\mathbf{A}}_0}}&{\mathbf{A}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathop {\mathbf{X}}\limits^\cdot }_0}} \\ {\mathop {\mathbf{X}}\limits^\cdot } \end{array}} \right] + {{\mathbf{F}}_{\text{0}}} = {\boldsymbol{0}},\end{align} \tag{ 10 }$$

where A₀ and A are the resistive matrices of the head and bars (detailed expressions can be seen in equations (B16) and (B19) in appendix B), and F₀ = [F_extx , F_exty, τ_ext + x₀ F_{ext x} − y₀ F_extx ]^T is the vector of the external actuation, containing horizontal and horizontal forces (arising from the acoustic actuation) and a torque (about the origin, arising from the magnetic power) on the head.

At the ith joint, the velocities determined from the two connected bars should be equal, i.e. V _i (2l) = V _{i +1}(0). These constraints ensure the inextensibility conditions of a flexible flagellum [48, 49]. Since the first bar and the head are assumed to be rigidly connected, one can obtain below (2n + 1) equations in terms of the constraints:

$$\begin{align}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}}&{{{\mathbf{T}}_V}} \\ { - {\mathbf{I}}}&{\mathbf{B}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathop {\mathbf{X}}\limits^\cdot }_0}} \\ {\mathop {\mathbf{X}}\limits^\cdot } \end{array}} \right] = {\boldsymbol{0}},\end{align} \tag{ 11 }$$

where I is a 3 × 3 identity matrix, T _V the matrix of kinematic constraints of neighboring bars (see equation (B20)), and B the kinematic constraints between the first bar and the head (see equation (B23)).

The polymeric (viscoelastic) flagellum is described based on the discrete beam theory [50] in which the inherent damping in most polymers is nonnegligible [24]. The material of the micro-swimmer is considered to be a kind of isotropic polymer described by the Kelvin–Voigt viscoelastic model [59]. The moment induced by the torsion spring at the ith joint M_i is expressed in a forward difference formula as:

$$\begin{align}{M_i} = E\frac{{{J_i} + {J_{i + 1}}}}{2}\frac{{{\theta _{i + 1}} - {\theta _i}}}{{2l}} + \eta \frac{{{J_i} + {J_{i + 1}}}}{2}\frac{{{{\mathop \theta \limits^\cdot }_{i + 1}} - {{\mathop \theta \limits^\cdot }_i}}}{{2l}},\end{align} \tag{ 12 }$$

where E and η are the Young's modulus and viscosity coefficient, respectively [60–62], and J_i denotes the second moment of area of the ith bar, where an average value of J between two neighboring bars is adopted to indicate the second moment of area at the ith joint in case the flagellum is non-uniform. Then, the moment equilibrium of the structure leads to the below (n − 1) equations (cf appendix B):

$$\begin{align}\left[ {\begin{array}{*{20}{c}} {{{\mathbf{D}}_0}}&{\mathbf{D}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathop {\mathbf{X}}\limits^\cdot }_0}} \\ {\mathop {\mathbf{X}}\limits^\cdot } \end{array}} \right] = {{\mathbf{F}}_E},\end{align} \tag{ 13 }$$

where D and D₀ are the matrix of moment balance related to the bars and head, respectively, which are expressed in equations (B25) and (B29), and F _E is a list of torques, which is expressed in equation (B31) in appendix B.

Finally, combining equations (10)–(13), the governing ODEs of the bar-joint model based on the CRFT can be expressed as:

$$\begin{align}\left[ {\begin{array}{*{20}{l}} {{{\mathbf{A}}_0}}&{\mathbf{A}} \\ {\boldsymbol{0}}&{{{\mathbf{T}}_V}} \\ { - {\mathbf{I}}}&{\mathbf{B}} \\ {{{\mathbf{D}}_0}}&{\mathbf{D}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathop {\mathbf{X}}\limits^\cdot }_0}} \\ {\mathop {\mathbf{X}}\limits^\cdot } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {{\mathbf{F}}_{\text{0}}}} \\ {\boldsymbol{0}} \\ {\boldsymbol{0}} \\ {{{\mathbf{F}}_E}} \end{array}} \right].\end{align} \tag{ 14 }$$

The nondimensionalization of equation (14) can be taken by choosing the half-length of a bar l, the period of the acoustic actuation T = 1/f_r , and the fluid viscosity μ as the reference variables based on the Buckingham π theorem [63]. The dimensionless quantities are detailed in appendix B. Equation (14) can be solved by using the built-in ODE solver ode15s in MATLAB [64].

We consider an AFMS requiring acoustic excitations to propel. The concerned AFMS is composed of a rigid head and a flexible slender flagellum, as exhibited in figure 3, where the structure has been meshed for the FSI simulation. The AFMS suspends in an aqueous solution and swims under an acoustic actuation that leads to the oscillations of the head and the whipping of the flagellum. According to the study of Liu and Ruan [24], the effect of periodic acoustic excitation for an AFMS can be assumed to be a periodic force F_ext acting on the head with the actuation frequency f_r , that is ${F_{{\text{ext}}}}\left( t \right) = {F_a}\cos \omega t$, where ω= 2πf_r is the angular frequency and F_a is the amplitude. The head oscillation induced by F_ext whips the flagellum, giving rise to the advancement of the AFMS at an average speed of V_ave. We focus firstly on the swimming problem that the average displacement of the micro-swimmer during several actuation periods keeps a straight line. It has been demonstrated both experimentally [14, 65] and theoretically [24, 40, 47] that the micro-swimmer tends to align itself perpendicular to the direction of F_ext. Hence, in the following FSI simulations, the only boundary condition is that the centroid of the head (or the end of the flagellum if the micro-swimmer has no head) is subject to a transverse harmonic force F_ext. Consequently, the head oscillates transversely with an amplitude A_amp.

Zoom In Zoom Out Reset image size

Figure 3 shows the mesh of the FSI model that was simulated using a built-in fully coupled 2D FSI solver in COMSOL Multiphysics [66]. The whole geometry can be divided into two domains: a rectangular fluid domain of 2000 × 500 µm² (capacious enough for a micro-swimmer) and a solid one for the AFMS in the center of the fluid. The radius of the circular head and the length of the flagellum are denoted by R_H and L, respectively. For a 2D simulation, the cross-section of the flagellum is rectangular. The flagellar width W_y can be either uniform or tapered with a maximum of W_{y 0} at the root (the joint between the head and the flagellum). The micro-swimmer is assumed to have a uniform thickness of W_z .

We investigate the steady-state swimming performances using the trajectory of the center of the head (or the right end of the flagellum). In the numerical simulations, the incompressible creeping flow module without the inertial term (Stokes flow) is employed for the fluid domain, while the solid mechanics module with plane strain approximation for the solid domain [67]. The micro-swimmer consists of two material models: a viscoelastic material for the flagellum and a rigid domain for the head [59]. Considering the expected medical applications and the current experiments, the flagellar material is assumed to be a Kelvin–Voigt viscoelastic material for representing an organic polymer, such as the polypyrrole employed in [13] and the polyethylene glycol employed in [14]. F_ext is applied along the y-direction to the centroid of the head. If headless, F_ext is applied to the right end of the flagellum. The automatic remeshing technique is employed to tackle the problems induced by large rigid-body motions, where the fluid domain (i.e. the gray area in figure 3) is set as the moving mesh region based on the Winslow smoothing method [68], in which the mesh distortion parameter is set to unity.

We first exemplify 2D simulations of two AFMS models as shown in figures 4(a) and (b). The first (figure 4(a)) one is a uniform flagellum (without head) with dimensions: L = 200 μm, W_z = 40 µm, and W_y = 12 µm, and swims in a fluid with the dynamic viscosity: µ = 0.1 Pa·s. The second one (figure 4(b)) is a sperm-like micro-swimmer referring to the AFMS described in [14]. It has a thickness, W_z = 20 µm, and a cylindrical head with a radius: R_H = 20 µm, and a tapered flagellum: L = 200 μm and W_y = W_{y 0}(1 – λX), where 0 ⩽ X⩽ 1 is the length normalized by L, W_{y 0} = 20 μm, and the tapering index λ = 0.9. Both micro-swimmers have the material properties: E = 10 MPa and η = 10 Pa·s and in the static fluids (u = 0).

Zoom In Zoom Out Reset image size

Figures 4(a) and (b) illustrate the fluid pressure distribution at t = 0.005 s. An vertical acoustic force F_ext was applied at the right with a frequency of f_r = 1000 Hz and an amplitude of 2.7 μN and 4 µN, for the uniform flagellum (figure 4(a)) and the sperm-like micro-swimmer (figure 4(b)), respectively, to maintain an amplitude of A_amp = 20 μm at the actuated end. It is observed that the average trajectory within a period of acoustic excitation is nearly a straight line parallel to the x-axis. The average (over one period) velocity in the first two periods varies (transient motion) and afterwards stabilizes at 383.4 μm·s⁻¹ (figure 4(a)) or 5537 μm·s⁻¹ (figure 4(b)). For the micro-swimmer herein, considering the fluid density ρ ∼ 10³ kg·m⁻³, which leads to the Reynolds number Re = ρLV/μ = 7.668 × 10⁻⁴, which is sufficiently small to warrant the employment of the RFT. The stabilized terminal velocity V_ave is obtained with the simulation domain having around 10⁴ elements (figure 4(a)) or 2 × 10⁴ elements (figure 4(b)). Figure 4(c) illustrates the relative errors for the cases without and with head, respectively, in terms of the number of elements used in the FSI simulations. It is noted that the simulation results have converged (i.e. the deviation is within 5%) when the number of elements is more than 10⁴, indicating that the FSI simulation is convergent and stable.

The correction factor K is determined by comparing the terminal average velocities V_ave obtained from the FSI simulations and the CRFT-based bar-joint model (the convergence analysis of the bar-joint model can be seen in appendix C). As an example, figure 5(a) shows the results of CRFT-based bar-joint model when K sweeps from 0.1 to 100 for the case of a uniform flagellum without head. In solving the bar-joint model, we also set A_amp = 20 μm to determine V_ave (blue dotted curve in figure 5(a)), and the deviations (red curve in figure 5(a)) from the simulation result (383.4 μm·s⁻¹) has a sharp minimum at K = 3.2, which is less than 1%.

Zoom In Zoom Out Reset image size

In principle, adjusting K can lead to an excellent agreement between simulation and CRFT results. In the following estimate of K for different cases, we set the criterion of acceptable deviation in V_ave to be 5% to balance the computational efficiency and accuracy and demonstrate several relations between K and the flagellum slenderness α (α = W_y /L) in figures 5(b)–(d). In order to ensure that the flagellum is still a slender body, we set α in the range of 0.05–0.07 [31]. Figures 5(b) and (c) show that K varies almost linearly with α, and the slope depends on the thickness W_z, for both the uniform flagellum without a head and the sperm-like micro-swimmer with a rigid-cylindrical head, respectively. The detailed expressions of K in terms of α obtained through linear regression are listed in table 1.

In table 1, one can note that the existence of a head can change the slope of the linear K–α relation from a negative value to a positive one. Besides, as observed in figure 5(c), the variation of R_H can hardly influence the K–α curve. In addition, K doubles when the depth W_z increases from 20 to 40 µm, as shown in figures 5(b) and (c), and also doubles when L increases from 200 to 400 µm, as shown in figure 5(d). These results suggest that, roughly speaking, K is linearly dependent on the volume of the flagellum. In figure 5(d), α is sampled from 0.05 to 0.1, and the fitting equation obtained in the α range of 0.05–0.07 is used. In the 0.07–0.1 range, the difference between the fitting equation and the actual results of K are still below 3%.

When a tapered flagellum is discretized into be a bar-joint model, these bars have uniform dimensions. In this case, the fitting equations of K for a uniform flagellum may be directly employed. For example, consider the micro-swimmer with L = 200 µm and R_H = W_z = W_{y 0} = 20 µm, assume that the width of the ith bar, i.e. W_yi , follows the tapering equation (equation (C1)) and that the slenderness of the ith bar, i.e. α_i , equals W_yi /L, the correction factor of the ith bar, i.e. K_i , can then be estimated according to the fitting curve of K in figure 5(d). In this way we can directly obtain the terminal average velocity V_ave based on the CRFT model without implementing the FSI simulation. We herein designate this method as the direct-K method. Straightforward though this method is, its outcomes can deviate from those of the FSI simulation, and the corresponding error in terms of the tapering index λ is illustrated in blue dotted line in figure 6(a). For the classical RFT model, i.e. K_i = 0 for each bar, the error exceeds 90%. For the direct-K method, it is noted that the error is around 12% for λ > 0.2, which indicates that if we directly apply the fitting expression of K arising from the case of a uniform flagellum to the case of a tapered one, the error can decrease by about 80% from the results of RFT. However, the deviation is still considerable.

Zoom In Zoom Out Reset image size

More accurate results for a tapered flagellum can be achieved by assuming all bars having identical K, i.e. K_i = K. In this way, we treat a tapered flagellum to be an equivalent uniform one with an effective K, and thus this method is regarded as the effective-K method. The red circles in figure 6(a) represent the effective K at different λ, and the red solid line is the corresponding linear fitting: K = −1.29λ + 6.43. It is noted that in figure 6(a) for λ ⩽ 0.1, magnitude of the effective K is nearly identical to that of the case of λ = 0, and the error of the effective-K method is close to 0, indicating that the effect of tapering on hydrodynamics only works when λ > 0.1.

Note that the above results of K are obtained with the specific material and actuation parameters. We find that A_amp and η hardly affect the determination of K, while f_r, E and μ can have considerable influence. For example, for the micro-swimmer as shown in figure 4(b) in section 3.1, we varied f_r, E and μ to check their influence on K. It is noted that the combined term β = μf_r /E is an important non-dimensional parameter in the governing equations of flagellum's hydrodynamics. For example, β enters into the elastohydrodynamic penetration length [40]: $\sqrt[4]{{\left( {EJ} \right)/\left( {2\pi {c_ \bot }{f_r}} \right)}}$, proportional to$\sqrt[4]{{{1 \mathord{\left/ \right. } \beta }}}$, and the sperm number [69]: $L\sqrt[4]{{\left( {2\pi {c_ \bot }{f_r}} \right)/\left( {EJ} \right)}}$, proportional to $\sqrt[4]{\beta }$. Hence, we plot the relations between K and β in figure 6(b) with the corresponding fitting equations listed in table 1. It is noted that when β is between 5 and 20, the variations of K in terms of f_r, E and μ collapse to the same curve. When β >20, the effects of E and μ are still the same and tends to level off; nevertheless, the effect of head excitation frequency f_r tends to lean upwards. The difference in K increases with the increase in f_r and μ (or decrease in E) when β> 20 because raising f_r can considerably enhance V_ave (which leads to a larger K) as it can induce the change of the flagellum's wiggling profile (e.g. bringing about smaller wavelength). Instead, increasing μ or reducing E only leads to larger resistive forces under the same wiggling profile, which results in a small enhancement in V_ave.

In order to explicate the validity of the method for correcting RFT, we compare three head trajectories (or the trajectories of the right end of the flagellum for the case without head) calculated by FSI simulations and by CRFT. Figure 7 exhibits the trajectories and their errors from FSI simulations and the CRFT models for the acoustic-actuation period from 20 to 25. The corresponding geometries of AFMSs at four instances which equally divide the last period in the FSI simulation (outline) and the CRFT model (centerline) are given as the insets. The supplementary videos of the corresponding FSI simulations and relevant codes of the CRFT model for verification are provided in the supplementary documents.

Zoom In Zoom Out Reset image size

Figures 7(a), (c) and (e) show the trajectories of FSI simulations and CRFT models for the case of a uniform flagellum without a head, a uniform flagellum with a head, and a tapered flagellum with a head, respectively. To approximate the trajectories based on the CRFT model, the corresponding K are determined to be 3.2, 6.3, and 5.2, respectively, based on the fitting equations suggested in the last section. It is noted that each of the trajectories obtained based on the CRFT model has an offset from the FSI-simulated one, reflecting the slight difference in the propulsive forces. This is due to the determination of K, in which we allowed a 5% tolerance in the motility difference. The distance between the positions of the AFMS head obtained from the FSI simulation and those from the CRFT solution, after dividing by the flagellum length L, is defined as the position error of CRFT solution, which are shown in figures 7(b), (d) and (f) (from 20 to 25 periods) and their insets (for first 40 periods). One may note that in figures 7(b) and (f), there are some small spikes, whereas in figure 7(d) the errors are simply oscillating. Several extreme error points in figures 7(b), (d) and (f) are indicated in figures 7(a), (c) and (e), respectively, where the corresponding differences between the two solutions are indicated by dashed lines. It is noted that the spikes of errors are due to the oscillatory motions, resembling the interference of waves.

The insets in figures 7(b), (d), and (f) show that the position errors in 40 periods are well below 4%. Because of the small motility difference (i.e. the tolerance in determining K), the trajectory deviations are caused; nevertheless, these deviations are well bounded and remain small when normalized by the length of the flagellum. Such small errors indicate that the CRFT-based bar-joint model, as described in section 3, renders an excellent approximation to numerical simulations. It should be noted that such position errors can be further reduced if we set a more stringent requirement in accuracy in determining K.

We expect that the proposed CRFT-based bar-joint model can ultimately be employed for the path planning and navigation control of AFMSs because of the excellent agreement between the FSI and the simplified model has been revealed. Here, as the final demonstration, we employ CRFT to study the turning strategy for an AFMS, which is essential and inevitable in practical applications. It is noted that FSI simulation for turning is unconverted in our trials due to excessively distorted elements; hence, we only provide the results of CRFT model.

In the following demonstration, the AFMS is identical to the one shown in figure 7(e). We expect that the external force F_ext (F_ext = [F_extx , F_exty ]^T) can be applied to the head of an AFMS through applying acoustic waves. According to [24], under an acoustic planewave rendering a 2D fluid displacement field A_acou with the amplitude A_acou, the oscillation amplitude of the head (R_H = 20 μm) of the AFMS A_amp is approximately 1.67A_acou if the sound wavelength is far greater than the dimensions of the head. In this case, the oscillatory velocity in the fluid is nonzero, i.e. u = ∂A_acou/∂t.

One may note from the trajectories in figure 7 that an AFMS tends to move along the direction perpendicular to F_ext that induces head oscillations, which is also the case in the experiments of acoustic actuation [65]. Hence, steering an acoustically propelled AFMS might be achieved by rotating F_ext and we can design such a turning method using only acoustic power based on the CRFT model. Considering an acoustic planewave that can gradually change the direction of propagation, hence changing the direction of F_ext on the head of an AFMS. As shown in figure 8(a), we let the angle of the acoustic wavevector rotates from 90° to 60° (counterclockwise from positive direction of x-axis) in 10 periods. After the change, the angle of the micro-swimmer is around 5°. Maintaining F_ext at 60° for 100 periods, the angle of the AFMS gradually drift to be close to 60°. Further turning F_ext to 150° in 10 periods leads the angle of the micro-swimmer to move along the straight line of 60°. The detailed trajectory of the AFMS is shown in figure 8(b).

Zoom In Zoom Out Reset image size

The above example shows that using only acoustic power to control AFMSs is possible in theory, although the orientation relation between an AFMS and acoustic wavevector is nontrivial, depending on the parameters of the actuation and the swimmer [70]. If the direction of wavevectors can be accurately adjusted, these waves can induce head oscillation and flagellum wiggling for propulsion and steering. In these attempts, our AFMS model can provide the formalism for the trajectory prediction and manipulation, which will lay the foundation for controlling micro-swimmers to carry out precision therapies for cancers and other diseases requiring non-invasive and localized treatments.

If a torque can be applied to the head of an AFMS, turning will be easier and more controllable. This can be achieved by introducing magnetic particles (MPs) into the head [71], resulting in a magnetized head along an easy axis (M), and applying an external magnetic field (B) unparallel with the head's easy axis, resulting in a torque τ_ext = M× B on the head, as sketched in the inset in figure 9(b). Note that the acoustic power is still used for propulsion but not for steering. The external magnetic field used in hyperthermia-based therapy is usually around 10 mT [72] and that the volumetric magnetic susceptibility is in the order of 10⁵ for ferromagnetic MPs [73], accordingly, the magnetic-induced torque by single MP of 1 μm is around 4 × 10⁻¹² N·m. With five MPs in the head, the total magnetic torque on the head is 2 × 10⁻¹¹ N·m, which is assumed in the following calculation. For such a magneto-acoustic method, we design the variations of F_extx , F_exty and τ_ext as shown in figure 9(a), which results in the trajectory of the AFMS as shown in figure 9(b). Hence, the proposed magneto-acoustic maneuver is feasible based on the existing technologies and our theoretical analysis of magneto-acoustic steering can be employed for targeted therapies which require the accurate and rapid navigation.

Zoom In Zoom Out Reset image size