Squirming Sphere in an Ideal Fluid

Introduction

The purpose of this work is to derive analytically and with an explicit formula the motion due to the variation of the surface of a squirming body in a fluid. Such a derivation has been done before, but only in two dimensions, and there is no literature on any attempts to generalize to three dimensions, which is considerably more involved algebraically. But utilizing software such as Mathematica, the derivation is quite straightforward.

The primary motivation behind this work is the search for innovative propulsion mechanisms. Surface variations may not lead to a very high speed, but they provide higher maneuverability and efficiency and are the primary source of motion of aquatic animals. Two types of fluid are studied - an Eulerian (potential) fluid where all viscous effects are ignored and inertia is the only factor, and a Stokesian (creeping) fluid where all inertial effects are ignored and viscous friction is the only factor. We consider only those cases because they linearize the equations governing the fluid motion, and permit analytic solutions.

Hydromechanical Connections

Consider a body immersed in a fluid. We regard its shape as a shape manifold \(M\). And the configuration manifold

\[Q=M\times G\]

represents the complete position and orientation of the body in space.

\(Q=M\times G\) can be thought of as a trivial bundle with a group structure G, and an action of \(G\) on \(Q\) given by \(\Phi_{h}(r,g)=(r,hg)\) for \((r,g)\in Q\) and \(h\in G.\) An element of the tangent space is given by \(\dot{q}=(\dot{r},\dot{g})\in TQ\).

We are interested in computing a connection on the configuration bundle \(Q\) that will allow us to compute the fiberwise translations resulting from cyclical changes in shape.

A connection is a vector valued one-form \(\Gamma:TQ\to\mathfrak{g}\) on \(Q\).

We follow [Kel98] in defining two types of connections:

  • Mechanical connection:

    \[\Gamma_{\text{mech}}:(q,\dot{q})\to\mathbb{I}^{-1}(q)J(q,\dot{q})\]
    where \(\mathbb{I}:\mathfrak{g}\to\mathfrak{g}^{*}\) is the locked inertia tensor such that
    \[\langle\mathbb{I}(q)\xi,\eta\rangle=\ll\xi_{Q}(q),\eta_{Q}(q)\gg_{KE}\]

    and

    \[\mathbb{I}(r,g)=Ad_{g^{-1}}^{*}\mathbb{I}_{\text{loc}}(r)Ad_{g^{-1}}\]
    and \(J:TQ\to\mathfrak{g}^{*}\) is the momentum map.
  • Stokes connection:

\[\Gamma_{\text{Stokes}}:(q,\dot{q})\to\mathbb{V}^{-1}(q)K(q,\dot{q})\]

where \(K:TQ\to\mathfrak{g}^{*}\) is a momentum map

\[\langle K(q,\dot{q}),\xi\rangle=\langle F(q,\dot{q}),\xi_{Q}(q)\rangle\]
associated with a dissipative force \(F\), and \(\mathbb{V}:\mathfrak{g}\to\mathfrak{g}^{*}\) is a viscosity tensor such that

\[\langle\mathbb{V}(q)\xi,\eta\rangle=\langle F(\xi_{Q}(q)),\eta_{Q}(q)\rangle\]

and

\[\mathbb{V}(r,g)=Ad_{g^{-1}}^{*}\mathbb{V}_{\text{loc}}(r)Ad_{g^{-1}}\]

Kelly uses the following local expressions

\[\begin{split}\begin{aligned} \Gamma_{\text{mech}} & = & Ad_{g}(g^{-1}\dot{g}+A_{\text{mech}}\dot{r})\\ \Gamma_{\text{mech}} & = & Ad_{g}(g^{-1}\dot{g}+A_{\text{Stokes}}\dot{r})\end{aligned}\end{split}\]

and for zero initial momentum

\[\begin{split}\begin{aligned} g^{-1}\dot{g} & = & -A_{\text{mech}}\dot{r}\\ g^{-1}\dot{g} & = & A_{\text{Stokes}}\dot{r}\end{aligned}\end{split}\]

which completely describe the evolution of a Potential and Stokesian swimmers respectability.

Squirming Sphere

Consider a three dimensional nearly spherical body whose surface is described by the following equation:

\[r=F(\theta,\varphi;s_{1},\dots,s_{n})\]

Such a body has already been studied previously numerically [Ruf05], and the current method that allows us to obtain explicit formulas for the displacement as will be later seen, are valuable tools for validating numerical simulations.

Coordinate system.

Coordinate system.

As an example we consider the propulsive movement of a roughly spherical device whose boundary shape is modulated by a small amount. We parametrize the boundary in terms of the shape variables \(s_{i}(t)\), and express it as a sum of an unperturbed base shape and a perturbed part:

\[F(\theta;s)=1+\epsilon\,\left(s_{1}\cos2\theta+s_{2}\cos3\theta\right)\]

The configuration manifold is

\[Q=\mathbb{R}^{2}\times SE(3)\]
Boundary perturbations.

Boundary perturbations.

Potential Flow

Irrotational potential flow is best described in terms of the velocity potential, which satisfies

\[\nabla^{2}\phi=0\]

with the Neumann boundary condition

\[\begin{equation}\boldsymbol{\nabla}\phi\cdot\boldsymbol{n}\big|_{r=F}=\boldsymbol{u}\cdot\boldsymbol{n}\big|_{r=F}\label{eq:boundary}\end{equation}\]

where \(\boldsymbol{u}\big|_{r=F}\) is the prescribed boundary motion.

For an axisymmetric problem, such as this one, the solution of the Laplace equation can be expressed as a series

\[\phi(r,\theta)=\sum_{k=0}^{\infty}\left(A_{k}r^{k}+\frac{B_{k}}{r^{k+1}}\right)P_{k}(\cos\theta)\]

where \(P_{k}(x)\) are the Legendre polynomials satisfying the orthogonality condition:

\[\int_{-1}^{+1}P_{m}(x)P_{n}(x)dx=\frac{2}{2n+1}\delta_{mn}\]

We require that \(\boldsymbol{u}\to0\) as \(r\to\infty\), and therefore we must take \(A_{k}=0\) for \(\forall k\).

\[\phi(r,\theta)=\sum_{k=0}^{\infty}a_{k}\frac{P_{k}(\cos\theta)}{r^{k+1}}\equiv\sum_{k=0}^{\infty}a_{k}\Phi_{k}\]

The problem of solving for \(\phi\) is reduced to finding the coefficients \(a_{k}\).

First, in order to enforce the boundary conditions on the perturbed surface, we need to compute the normal \(\boldsymbol{n}\) at each point. The two tangential components on the surface are

\[\begin{split}\begin{aligned} \boldsymbol{t_{1}} & =\frac{\partial}{\partial\theta}(F\boldsymbol{\hat{r}})=\frac{\partial F}{\partial\theta}\hat{\boldsymbol{r}}+F\boldsymbol{\hat{\theta}}\\ \boldsymbol{t_{2}} & =\frac{\partial}{\partial\varphi}(F\boldsymbol{\hat{r}})=\frac{\partial F}{\partial\varphi}\boldsymbol{\hat{r}}+F\sin\theta\,\boldsymbol{\hat{\varphi}}\end{aligned}\end{split}\]

and therefore the normal component can be evaluated as the cross product of the two

\[\begin{split}\begin{aligned} \boldsymbol{n} & =\boldsymbol{t_{1}}\times\boldsymbol{t_{2}}\\ & =F^{2}\sin\theta\:\boldsymbol{\hat{r}}-F\frac{\partial F}{\partial\theta}\sin\theta\,\boldsymbol{\hat{\theta}}-F\frac{\partial F}{\partial\varphi}\,\boldsymbol{\hat{\varphi}}\end{aligned}\end{split}\]

Note that we need not normalize the expression for \(\boldsymbol{n}\) since it occurs on both sides of the boundary condition equation Eq.\ref{eq:boundary} . Expanding the spherical coordinates \((r,\theta,\varphi)\) of \(\boldsymbol{n}\) in terms of \(\epsilon\):

\[\begin{split}\begin{aligned} \boldsymbol{n} & =\mathbf{n}^{(0)}+\epsilon\mathbf{n}^{(1)}+\epsilon\mathbf{n}^{(2)}\\ & =\left(\begin{array}{c} 1\\ 0\\ 0 \end{array}\right)+\epsilon\left(\begin{array}{c} 2\left(c[2\theta]s_{1}+c[3\theta]s_{2}\right)\\ 2s[2\theta]s_{1}+3s[3\theta]s_{2}\\ 0 \end{array}\right) \\ & \quad +\epsilon^{2}\left(\begin{array}{c} \left(c[2\theta]s_{1}+c[3\theta]s_{2}\right){}^{2}\\ \left(c[2\theta]s_{1}+c[3\theta]s_{2}\right)\left(2s[2\theta]s_{1}+3s[3\theta]s_{2}\right)\\ 0 \end{array}\right)\end{aligned}\end{split}\]

where we have used the abbreviations \(c[\theta]=\cos(\theta)\) and \(s[\theta]=\sin(\theta)\).

Next, we expand the potential in a power series in terms of the small perturbation parameter \(\epsilon\) .

\[\phi=\phi^{(0)}+\epsilon\phi^{(1)}+\epsilon^{2}\phi^{(2)}+...+\epsilon^{i}\phi^{(i)}+...\]

where

\[\phi^{(i)}=\sum_{k=0}^{\infty}a_{k}^{\phantom{k}(i)}\Phi_{k}\]

Plugging back in Eq.\ref{eq:boundary} , we expand both sides in terms of \(\epsilon\)

\[\begin{split}\begin{aligned} \nabla\phi\cdot\boldsymbol{n}\big|_{r=F} & =\sum_{i=0}^{\infty}\epsilon^{i}\sum_{k=0}^{\infty}a_{k}^{(i)}f_{k}^{(i)}\\ \boldsymbol{u}\cdot\boldsymbol{n}\big|_{r=F} & =\sum_{i=0}^{\infty}\epsilon^{i}u^{(i)}\end{aligned}\end{split}\]

and equate terms of equal power

\[\sum_{k=0}^{\infty}a_{k}^{(i)}f_{k}^{(i)}=u^{(i)}\]

where the first few terms \(f_{k}^{(0)}\) are

\[\begin{split}\begin{aligned} f_{0}^{(0)} & =-1\\ f_{1}^{(0)} & =-2\eta\\ f_{2}^{(0)} & =\frac{3}{2}\left(1-3\eta^{2}\right)\\ f_{3}^{(0)} & =2\left(3\eta-5\eta^{3}\right)\\ f_{4}^{(0)} & =-\frac{5}{8}\left(3-30\eta^{2}+35\eta^{4}\right)\\ f_{5}^{(0)} & =-\frac{3}{4}\left(15\eta-70\eta^{3}+63\eta^{5}\right)\\ f_{6}^{(0)} & =\frac{7}{16}\left(5-105\eta^{2}+315\eta^{4}-231\eta^{6}\right)\\ f_{7}^{(0)} & =\frac{1}{2}\left(35\eta-315\eta^{3}+693\eta^{5}-429\eta^{7}\right)\end{aligned}\end{split}\]

where \(\eta=\cos\theta\).

By the linearity of Laplace’s equation (using superposition) one can write following Kirchoff,

\[\phi=\phi_{z}\dot{z}+\phi_{s1}\dot{s}_{1}++\phi_{s2}\dot{s}_{2}\]

Then, up to second order in \(\epsilon\)

\[\begin{split}\begin{aligned} \phi_{z} & =-\frac{P_{1}}{2r^{2}}+\epsilon\left(\frac{9P_{1}s_{1}}{10r^{2}}+\frac{9P_{2}s_{2}}{7r^{3}}-\frac{6P_{3}s_{1}}{5r^{4}}-\frac{48P_{4}s_{2}}{35r^{5}}\right)\\ & \quad +\epsilon^{2}\bigl(\frac{P_{1}\left(-3339s_{1}^{2}-4990s_{2}^{2}\right)}{2450r^{2}}-\frac{1781P_{2}s_{1}s_{2}}{735r^{3}}+\frac{6P_{3}\left(308s_{1}^{2}-45s_{2}^{2}\right)}{1925r^{4}}\\ & \qquad \qquad +\frac{42264P_{4}s_{1}s_{2}}{13475r^{5}}-\frac{16P_{5}\left(91s_{1}^{2}-90s_{2}^{2}\right)}{637r^{6}}-\frac{64P_{6}s_{1}s_{2}}{11r^{7}}-\frac{512P_{7}s_{2}^{2}}{143r^{8}}\bigg)\\ \\ \phi_{s1} & =\epsilon\left(\frac{P_{0}}{3r}-\frac{4P_{2}}{9r^{3}}\right) + \epsilon^{2}\bigg(-\frac{14P_{0}s_{1}}{15r}-\frac{3P_{1}s_{2}}{35r^{2}} + \frac{8P_{2}s_{1}}{21r^{3}}+\frac{2P_{3}s_{2}}{5r^{4}} \\ & \qquad \qquad \qquad \qquad \qquad \qquad -\frac{512P_{4}s_{1}}{525r^{5}}-\frac{64P_{5}s_{2}}{63r^{6}}\bigg)\\ \\ \phi_{s2} & =\epsilon\left(\frac{3P_{1}}{10r^{2}}-\frac{2P_{3}}{5r^{4}}\right)+\epsilon^{2}\bigg(-\frac{34P_{0}s_{2}}{35r}-\frac{309P_{1}s_{1}}{350r^{2}}-\frac{13P_{2}s_{2}}{315r^{3}}\\ & \qquad \qquad \qquad \qquad \qquad \qquad +\frac{122P_{3}s_{1}}{225r^{4}}+\frac{1128P_{4}s_{2}}{1925r^{5}}-\frac{208P_{5}s_{1}}{189r^{6}}-\frac{1856P_{6}s_{2}}{1617r^{7}}\bigg)\end{aligned}\end{split}\]

The total kinetic energy is given by

\[\begin{split}\begin{aligned} L & =\frac{1}{2}\int_{\mathfrak{F}}dV\,|\nabla\phi|^{2}\\ & \frac{1}{2}\int_{0}^{2\pi}d\phi\,\int_{0}^{\pi}\sin\theta d\theta\int_{r=F}^{\infty}r^{2}dr\:|\nabla\phi|^{2}\\ & =\pi\int_{-1}^{1}d\eta\int_{r=F}^{\infty}r^{2}dr\:|\nabla\phi|^{2}\\ & =\frac{\pi}{5}\dot{z}^{2}-\epsilon\frac{\pi}{5}\text{ }\left(\frac{49}{15}\text{ }s_{1}\dot{z}^{2}+\frac{6}{5}\dot{s}_{2}\dot{z}\right)\\ & \quad +\epsilon^{2}\frac{\pi}{5}\bigg(\frac{230}{189}\dot{s}_{1}^{2}+\frac{241}{225}\dot{s}_{2}^{2}+2\left(\frac{1069}{525}s_{1}\dot{s}_{2}-\frac{628}{245}s_{2}\dot{s}_{1}\right)\dot{z}+\\ & \qquad \qquad + \left(\frac{576}{175}s_{1}^{2}+\frac{144857}{18865}s_{2}^{2}\right)\dot{z}^{2}\bigg)\\ & \quad +O(\epsilon^{3})\end{aligned}\end{split}\]

This Lagrangian is invariant under translation along \(z\), so its momentum map is

\[\begin{split}\begin{aligned} J & = \frac{2\pi}{5}\text{ }\dot{z}-\epsilon\frac{2\pi}{5}\left(\frac{49}{15}\text{ }s_{1}\dot{z}+\frac{3}{5}\dot{s}_{2}\right) \\ & \quad +\epsilon^{2}\frac{2\pi}{5}\left(\frac{1069}{525}s_{1}\dot{s}_{2}-\frac{628}{245}s_{2}\dot{s}_{1}+\left(\frac{576}{175}s_{1}^{2}+\frac{144857}{18865}s_{2}^{2}\right)\dot{z}\right)+O(\epsilon^{3})\end{aligned}\end{split}\]

The locked inertia tensor is given by

\[\mathbb{I}(s_{1},s_{2}):\quad\xi\mapsto\frac{2\pi}{5}\text{ }\xi\left(1-\epsilon\frac{49}{15}\text{ }s_{1}+\epsilon^{2}\left(\frac{576}{175}s_{1}^{2}+\frac{144857}{18865}s_{2}^{2}\right)\right)+O(\epsilon^{3})\]

Therefore,

\[\begin{split}\begin{aligned} \Gamma_{\text{mech}} & =\mathbb{I}^{-1}J \\ & =dz-\epsilon\frac{3}{5}ds_{2}+\epsilon^{2}\left(\frac{8}{105}s_{1}\dot{s}_{2}-\frac{628}{245}s_{2}\dot{s}_{1}\right)+O(\epsilon^{3})\end{aligned}\end{split}\]

Self-propulsion of the sphere requires that its motion remain horizontal with respect to this connection, therefore

\[\dot{z}=\epsilon\frac{3}{5}\dot{s}_{2}+\epsilon^{2}\left(\frac{628}{245}s_{2}\dot{s}_{1}-\frac{8}{105}s_{1}\dot{s}_{2}\right)\]

Stokes Flow

The solution of creeping-motion can be reduce in two dimensions to the biharmonic equation in terms of the stream function

\[\nabla^{4}\psi=0\]

which satisfies the no-slip and no-penetration boundary condition on the surface

\[u\big|_{r=F}=(u_{r},u_{\theta})=(\dot{F},0)\]

Following [Lea07] p.462 the general solution is

\[\psi=\sum_{k=1}^{\infty}(A_{k}r^{k+3}+B_{k}r^{k+1}+C_{k}r^{2-k}+D_{k}r^{-k})Q_{k}(\eta)\]

where \(Q_{k}(x)\) are polynomial functions defined in terms of the Legendre polynomials

\[Q_{k}(x)=\int_{-1}^{x}P_{k}(x)dx\]

and satisfying the orthogonality condition

\[\int_{-1}^{+1}\frac{Q_{m}(x)Q_{n}(x)}{1-x^{2}}dx=\frac{2}{n(n+1)(2n+1)}\delta_{mn}\]

Again we want the velocity to vanish at infinity, so we will neglect the first two terms

\[\psi=\sum_{k=1}^{\infty}(a_{k}r^{2-k}+b_{k}r^{-k})Q_{k}(\eta)\]
\[u_{r}=-\frac{1}{r^{2}}\frac{\partial\psi}{\partial\eta}\qquad u_{\theta}\sqrt{1-\eta^{2}}=-\frac{1}{r}\frac{\partial\psi}{\partial r}\]

Again we will exploit the linearity of the equation by using Kirchoff’s method to express the streamfunction as

\[\psi=\psi_{z}\dot{z}+\psi_{s1}\dot{s}_{1}++\psi_{s2}\dot{s}_{2}\]
\[\begin{split}\begin{aligned} \psi_{z} & =\left(\frac{1}{2r}-\frac{3r}{2}\right)Q_{1}+\\ & +\epsilon\bigg[Q_{3}\left(\frac{12s_{1}}{5r^{3}}-\frac{12s_{1}}{5r}\right)+Q_{1}\left(-\frac{9s_{1}}{10r}+\frac{9rs_{1}}{10}\right)\\ & \qquad \qquad + Q_{4}\left(\frac{24s_{2}}{7r^{4}}-\frac{24s_{2}}{7r^{2}}\right)+Q_{2}\left(\frac{27s_{2}}{14}-\frac{27s_{2}}{14r^{2}}\right)\bigg]\\ & +\epsilon^{2}\bigg[Q_{1}\left(-\frac{3r\left(448s_{1}^{2}+1045s_{2}^{2}\right)}{1225}+\frac{4683s_{1}^{2}+8125s_{2}^{2}}{2450r}\right)\\ & \qquad \qquad +Q_{2}\left(-\frac{1287}{490}s_{1}s_{2}+\frac{2197s_{1}s_{2}}{490r^{2}}\right)\\ & \qquad \qquad +Q_{3}\left(-\frac{6\left(308s_{1}^{2}-135s_{2}^{2}\right)}{1925r}-\frac{4\left(308s_{1}^{2}+45s_{2}^{2}\right)}{1925r^{3}}\right)\\ & \qquad \qquad +Q_{4}\left(-\frac{10764s_{1}s_{2}}{2695r^{4}}-\frac{1836s_{1}s_{2}}{2695r^{2}}\right)+Q_{7}\left(\frac{2048s_{2}^{2}}{143r^{7}}-\frac{1536s_{2}^{2}}{143r^{5}}\right)\\ & \qquad \qquad +Q_{5}\left(\frac{48\left(91s_{1}^{2}-54s_{2}^{2}\right)}{637r^{5}}-\frac{32\left(91s_{1}^{2}-18s_{2}^{2}\right)}{637r^{3}}\right)\\ & \qquad \qquad +Q_{6}\left(\frac{224s_{1}s_{2}}{11r^{6}}-\frac{160s_{1}s_{2}}{11r^{4}}\right)\bigg]\end{aligned}\end{split}\]

Due to the lack of time, I am unable to continue with the calculations. The Mathematica code that I use seems to have trouble computing the \(\psi_{s1}\) and \(\psi_{s2}\) components. This could be related the 3D equivalent of the Stoke’s paradox, or could simply be a bug in the code. In any case, it would be important to carry this calculation to completion at some point in the future in order to compare the displacement resulting from creeping flow to that in potential flow.

Discussion

A calculation similar to the one done above was carried by Kelly for a two dimensional circle where the algebra is less involved. Figure [Flo:comp] compares the results. In all cases the shape trajectory is a unit circle \((s_{1},s_{2})\) shape space.

Geometric phase for a particular gate. TODO: stokes 3d.

Geometric phase for a particular gate. TODO: stokes 3d.

Explicitly, we quote the results. For the the potential flow case, Kelly obtained:

\[\dot{z}=\epsilon^{2}(s_{2}\dot{s}_{1}-s_{1}\dot{s}_{2})+O(\epsilon^{3})\]

and for Stoke’s case

\[\dot{z}=-\epsilon^{2}\left(\frac{1}{4}s_{2}\dot{s}_{1}+\frac{1}{2}s_{1}\dot{s}_{2}\right)+O(\epsilon^{2})\]

If we compare those, to our result

\[\begin{equation}\dot{z}=\epsilon\frac{3}{5}\dot{s}_{2}+\epsilon^{2}\left(\frac{628}{245}s_{2}\dot{s}_{1}-\frac{8}{105}s_{1}\dot{s}_{2}\right)\label{eq:magic}\end{equation}\]

The first thing that we notice is the presence of \(O(\epsilon^{1})\) in the 3d case (our solution). However, this term has no effect on the net displacement and it only results in a more pronounced oscillation of the center of mass during each cycle. Also we notice that the body achieves higher displacement in 3d.

This simplicity with which such explicit equations such as Eq.\ref{eq:magic} can be obtained is a testament to the power of the geometric approach to mechanics.

References

[Blo07]Anthony Bloch. Nonholonomic Mechanics and Control. Springer, September 2007. ISBN 0387955356.
[Hol08]Darryl D. Holm. GEOMETRIC MECHANICS: Dynamics and Symmetry. Imperial College Press, illustrated edition edition, May 2008. ISBN 1848161956.
[Kel98]Scott D Kelly. The Mechanics and Control of Robotic Locomotion with Appliations to Aquatic Vehicles. PhD thesis, California Institute of Technology, May 1998.
[KP04]M. Kohr and I. Pop. Viscous Incompressible Flow for Low Reynolds Numbers (Advances in Boundary Elements). WIT Press (UK), illustrated edition edition, July 2004. ISBN 1853129917.
[Lea07]L. Gary Leal. Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes. Cambridge University Press, 1 edition, June 2007. ISBN 0521849101.
[MR02]Jerrold E. Marsden and Tudor S. Ratiu. Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems. Springer, 2nd edition, December 2002. ISBN 038798643X.
[MB99]R. Mason and J. Burdick. Propulsion and control of deformable bodies in an ideal fluid. In 1999 IEEE International Conference on Robotics and Automation, 1999. Proceedings, volume 1. 1999.
[MRR06]J. B. Melli, C. W. Rowley, and D. S. Rufat. Motion planning for an articulated body in a perfect planar fluid. SIAM Journal on Applied Dynamical Systems, 5:650, 2006.
[NK07]S. Nair and E. Kanso. Hydrodynamically coupled rigid bodies. Journal of Fluid Mechanics, 592:393–411, 2007.
[Ruf05]Dzhelil S Rufat. Locomotion of a three dimensional variable shape body in inviscid fluid. May 2005.