Swimming of a sphere in a viscous incompressible fluid with inertia

Swimming of a sphere in a viscous incompressible fluid with inertia

B. U. Felderhof ufelder@physik.rwth-aachen.de Institut für Theorie der Statistischen Physik
RWTH Aachen University
Templergraben 55
52056 Aachen
   R. B. Jones r.b.jones@qmul.ac.uk Queen Mary University of London, The School of Physics and Astronomy, Mile End Road, London E1 4NS, UK
July 9, 2019

The swimming of a sphere immersed in a viscous incompressible fluid with inertia is studied for surface modulations of small amplitude on the basis of the Navier-Stokes equations. The mean swimming velocity and the mean rate of dissipation are expressed as quadratic forms in term of the surface displacements. With a choice of a basis set of modes the quadratic forms correspond to two Hermitian matrices. Optimization of the mean swimming velocity for given rate of dissipation requires the solution of a generalized eigenvalue problem involving the two matrices. It is found for surface modulations of low multipole order that the optimal swimming efficiency depends in intricate fashion on a dimensionless scale number involving the radius of the sphere, the period of the cycle, and the kinematic viscosity of the fluid.

47.15.G-, 47.63.mf, 47.63.Gd, 87.17.Jj

I Introduction

The swimming of fish and the flying of birds still pose problems to theory 1 (). The analysis can be based on the Navier-Stokes equations for the flow of a viscous fluid with a no-slip boundary condition at the surface of the body with periodically changing shape. For simplicity the fluid may be taken to be incompressible. The fluid is then characterized by its shear viscosity and mass density.

Most of the theoretical work has been concerned with either of two limiting situations. The swimming of microorganisms is well described by the time-independent Stokes equations of low Reynolds number hydrodynamics 2 ()3 (). The work in this area has been reviewed by Lauga and Powers 4 (). In the opposite limit of inviscid flow the analysis is based on the Euler equations with the effect of viscosity relegated to a boundary layer. The flow is predominantly irrotational, apart from the boundary layer and a wake of vorticity. The work in this field was reviewed by Sparenberg 5 ()6 () and by Wu 7 ()8 (). The modeling of bird flight was reviewed by Pennycuick 9 () and by Shyy et al. 10 (). The problem has also been addressed in computer simulation 11 ().

It is important to have a model covering the full range of kinematic viscosity. The seminal work of Taylor 12 () on the swimming of a sheet in the Stokes limit was extended to a fluid with inertia by Reynolds 13 () and by Tuck 14 (). The swimming of a planar slab has also been studied in the full range 15 (). The disadvantage of these models is the infinite length of the system which precludes study of the finite size effects which are believed to be crucial in the inviscid limit.

In earlier work we have studied small amplitude swimming of a body in a viscous fluid with inertia from a general point of view 16 (). As an example we studied the swimming of a sphere by means of potential flow 17 (). Later we showed that in the Stokes limit the addition of viscous modes leads to a significantly enhanced optimal efficiency 18 (). In the following we study the swimming of a sphere in the full range of kinematic viscosity. For simplicity we assume axisymmetric flow.

The effect of Reynolds stress turns out to be quite important. At small values of the kinematic viscosity it largely cancels the effect of viscous stress and pressure in the optimal efficiency. As a consequence the efficiency varies little with kinematic viscosity in a wide range. Similar behavior was found for an assembly of interacting rigid spheres 19 () and for a planar slab 15 (). We find that in swimming of a single sphere with surface distortions consisting of a running wave of three or five low order multipolar modes the optimal efficiency in the inviscid limit tends to the value for potential swimming. We expect that this feature holds more generally.

Ii Flow equations

We consider a flexible sphere of radius immersed in a viscous incompressible fluid of shear viscosity and mass density . In the laboratory frame, where the fluid is at rest at infinity, the flow velocity and the pressure satisfy the Navier-Stokes equations


The fluid is set in motion by time-dependent distortions of the sphere. We shall study axisymmetric periodic distortions which lead to a translational swimming motion of the sphere in the direction in a Cartesian system of coordinates. The surface displacement is defined as the vector distance


of a point on the displaced surface from the point on the sphere with surface . The fluid velocity in the rest frame is required to satisfy


This amounts to a no-slip boundary condition. The instantaneous translational swimming velocity and the flow pattern follow from the condition that no net force is exerted on the fluid.

We perform a perturbation expansion in powers of the displacement . The first order flow velocity and pressure satisfy the linearized Navier-Stokes equations. The translational swimming velocity of the sphere, averaged over a period, is denoted by . To first order in displacements the mean swimming velocity vanishes. We have shown previously 16 () that to second order the mean swimming velocity may be calculated as the sum of a surface and a bulk contribution


In spherical coordinates the surface contribution may be expressed as an integral of a mean surface velocity , defined by


where the overhead bar indicates the time average over a period. The surface contribution to the swimming velocity is given by the spherical average 17 ()


The bulk contribution corresponds to the term in the Navier-Stokes equations. The time-averaged second order flow velocity and pressure satisfy the inhomogeneous Stokes equations 16 ()


with boundary condition


The right hand side in Eq. (2.7) represents a force density . The bulk part of the second order flow satisfies Eq. (2.7) with the no-slip boundary condition . The contribution to the mean swimming velocity follows from the flow velocity at infinity in the rest frame and the condition that no net force is exerted on the fluid. The integral of the force density is canceled by the surface integral of an induced force density on the sphere at rest. As we shall show later, the bulk contribution may be calculated with the aid of an antenna theorem 20 ().

In the following we consider in particular the case of harmonic time variation. It is then convenient to introduce complex notation. The surface displacement is written as


with complex amplitude . The corresponding first order flow velocity and pressure are given by


The time-averaged surface velocity may be expressed as


The time-averaged second order rate of energy dissipation is given by 16 ()


where is the first order stress tensor with Cartesian components


The efficiency of swimming is defined as


This quantity may be expressed as a ratio of two forms which are quadratic in the surface displacement, and involve two Hermitian matrices. Earlier we have studied the zero frequency limit of the problem where inertia plays no role 18 (). It has been shown by Shapere and Wilczek 18A () that in this limit the above definition of efficiency is preferable to that of Lighthill 18B (). The efficiency defined in Eq. (2.14) is essentially the ratio of speed and power and is relevant in the whole range of scale number.

Iii Matrix formulation

The explicit calculation requires the choice of a basis set of solutions of the linearized Navier-Stokes equations. After removal of the exponential time-dependent factor the equations read


with the variable


In our previous work 17 () we have chosen a set of modes identical to those used earlier in a hydrodynamic scattering theory of flow about a sphere 21 (). In the present axisymmetric case we can use a reduced set of solutions. Moreover, it turns out that in the numerical work it is advantageous to use a different normalization. Thus we use the modes


with modified spherical Bessel functions 22 () and vector spherical harmonics defined by


with Legendre polynomials and associated Legendre functions in the notation of Edmonds 23 (). The notation is identical to that used in Ref. 26. In particular and . The solutions are associated with vanishing pressure variation. The notation for the flow field is identical to that in our previous work for zero frequency 18 (). At zero frequency there is no pressure variation associated with these irrotational flow fields. We remark here that the above basis set is not suitable at low frequency owing to the singularity of the solutions at . At low frequency we must use a modified set, as discussed later.

We expand the first order flow velocity and pressure in the modes, given by Eq. (3.3), as


with complex coefficients , which can be calculated as moments of the function on the surface . Correspondingly the displacement function is expanded as


We define the complex multipole moment vector as the one-dimensional array


Then can be expressed as


with a dimensionless Hermitian matrix and the notation


where the subscript takes the two values with and . The subscripts correspond to notation used in earlier work 17 (),21 () for modes proportional to those in Eq. (3.3). We impose the constraint that the force exerted on the fluid vanishes at any time. This requires . We implement the constraint by dropping the first element of and erasing the first row and column of the matrix . We denote the corresponding modified vector as and the modified matrix as .

The time-averaged rate of dissipation can be expressed as


with a dimensionless Hermitian matrix . We denote the modified matrix obtained by dropping the first row and column by .

With the constraint the mean swimming velocity and the mean rate of dissipation can be expressed as


Optimization of the mean swimming velocity for given mean rate of dissipation, taken into account with a Lagrange multiplier , leads to the eigenvalue problem


Both matrices and are hermitian, so that the eigenvalues are real. With truncation at maximum -value the truncated matrices and are -dimensional. The maximum positive eigenvalue is of particular interest. Its corresponding eigenvector provides the swimming mode of maximal efficiency. The truncated matrices correspond to swimmers obeying the constraint that all multipole coefficients for vanish.

With use of Eq. (2.6) we find the contribution to the matrix . The remainder follows from the contribution to the mean swimming velocity. The elements of the matrices and are complex numbers which can be calculated by substitution of the expansions in Eq. (3.5) and (3.6) into the expressions (2.11) and (2.12). We have calculated the elements of the matrices and in our earlier work 17 (). In order to make contact with our subsequent work on the axisymmetric case in the limit of zero frequency 18 () we have performed an independent calculation for axial symmetry. In the calculation we use angular integrals which are detailed in Appendix A.

The matrix is diagonal in , and the matrices and are tridiagonal in . The matrices are frequency-dependent via the variable . We write


where is the dimensionless viscosity 25 (). The maximum eigenvalue depends on the variable . We call the scale number. It is related to the Roshko number , where , by , if in we use the sphere diameter as the characteristic length .

From Eq. (3.13) we see that for given fluid properties the frequency must decrease with increasing radius as in order to keep the scale number constant. For water the kinematic viscosity takes the value , and in air it takes the value . Hence for air we have with in cm and in Hz.

Iv Effect of Reynolds stress

The force density in Eq. (2.7) can be written alternatively as the divergence of a Reynolds stress . In order to find the contribution to the mean swimming velocity caused by this stress we must solve Eq. (2.7) with no-slip boundary condition at . We explained the procedure below Eq. (2.8). The mean swimming velocity corresponds to the matrix as in Eq. (3.8).

According to theory developed earlier 16 () the second order swimming velocity corresponding to the mean Reynolds stress can be written as a sum of two contributions,


where corresponds to the integral of the force density


according to Stokes’ law


The remainder corresponds to a solution of the Stokes problem Eq. (2.7) corresponding to the force density and a freely moving sphere of radius .

In order to find the contribution we use an antenna theorem derived earlier 20 (). The force density is decomposed into vector spherical harmonics as


with scalar functions given by the angular integrals


and with remainder given by a sum of higher order vector spherical harmonics . According to the antenna theorem 20 () the Green function solution of the Stokes equations, assumed to be valid in all space in the absence of the sphere, corresponds to a flow velocity for given by


with flow fields given by 24 ()




and remainder given by a sum of higher order vector spherical harmonics . The Green function tends to zero at infinity.

The velocity in Eq. (4.1) follows from Faxén’s theorem as 24 ()


By construction the flow pattern tends to plus a flow which decays faster than at large distance from the origin.

V Expansion at low frequency

As we have noted in Sec. III the singular behavior of the flow fields at leads to numerical difficulties at low frequency 21 (). Thus instead of we use


with coefficient


It may be checked that at zero frequency


where is the viscous mode function


used in the zero frequency theory 18 ().

The pressure corresponding to the velocity mode in Eq. (5.1) is given by


Correspondingly the expansion in Eq. (3.5) must be replaced by


With multipole vector defined by


the mean swimming velocity and mean rate of dissipation can be expressed as


with matrices and .

The above expansion is not suitable at high frequency owing to numerical difficulties in the eigenvalue problem analogous to Eq. (3.13). Thus we must use two different expansions in the two regimes of low and high frequency. The eigenvalues for the two eigenvalue problems are of course the same, and the eigenvectors are related by the matrix transforming one basis set into the other. At zero frequency the two matrices in the representation of this section are identical to those derived earlier 18 ().

Vi Results

With a small number of modes at low multipole order the calculations can be performed analytically. With higher multipoles included the generalized eigenvalue problem can be solved numerically by use of the Eigensystem command of Mathematica. In our earlier work 17 () we have performed calculations involving just potential modes. In that case the Reynolds stress vanishes. As we have shown at zero frequency 18 () the inclusion of viscous modes leads to a significantly higher maximum efficiency.

The qualitative behavior can be seen from a simple model with only modes of orders and included. In this case there are just three modes, the dipolar potential mode at , the viscous mode at , and the quadrupolar potential mode at . The modes are given by Eq. (3.3) at high frequency, and by Eqs. (5.1) and (5.5) for the viscous modes at low frequency. It turns out that in this case the high frequency expansion works well numerically over the whole range of interest. It suffices to compare with the zero frequency results 18 ().

From Eqs. (7.11) and (7.17) of Ref. 18 the matrix at zero frequency is given by


and the matrix is given by


We denote the corresponding maximum eigenvalue by . The eigenvalue problem yields 26 ()


The corresponding eigenvector has components .

With fluid inertia included the maximum eigenvalue as a function of the scale number is given by the expression


with numerator


and denominator


The expression is derived by use of the characteristic equation from the matrices and given explicitly in Appendix B. The denominator is related to the determinant of by .

In Fig. 1 we plot the function as a function of . The low frequency expansion is given by


in accordance with Eq. (6.3). The function shows a maximum value at , corresponding to the eigenvector with components for the modes of Eq. (3.3). For large values of the numerical calculation can be performed by a method discussed in the next section. At high scale number the contribution of the Reynolds stress to the mean swimming velocity is quite important. If it were omitted, then the maximum eigenvalue would grow as . Actually, the maximum eigenvalue tends to the limiting value


the optimal value for potential swimming with just the dipolar and quadrupolar modes. The eigenvalue and the corresponding eigenvector are found from Eqs. (6.1) and (6.2) with the second row and column of the matrices deleted.

The maximum eigenvalue at given increases as more modes are included. In Fig. 2 we show the maximum eigenvalue as a function of for surface displacements given by a superposition of the five modes for . In this case the numerical calculation is more demanding than for . We use the modes of Sec. V for and the modes of Eq. (3.3) for . The maximum eigenvalue at zero frequency is and the corresponding eigenvector has components . There is a weak maximum at . The corresponding eigenvector has components in terms of the modes of Sec. V. As increases beyond the function decays to , the optimal value for potential swimming with modes with corresponding eigenvector .

Vii Behavior for large scale number

The behavior shown in Figs. 1 and 2 requires further discussion of the asymptotic properties for large scale number. The function shown in Fig. 2 is calculated from a complicated analytic expression involving exponential integrals, like the expression Eq. (6.4) for . The expression is derived from the explicit expressions for the matrices and given in Appendix B, by use of the characteristic equation.

The exponential integrals appear multiplied by exponentials. It is therefore useful to define


In the expression for we need the value of at and at for large positive . Clearly the function is analytic in the complex plane apart from a branch cut along the negative real axis. The values of which are needed in Fig. 2 can be found accurately by numerical integration in Eq. (7.1). The expression for involves powers of up to and down to , so that in numerical calculations for large we need integer programming.

To scrutinize the behavior for large values of it is useful to derive series expansions. By integration by parts we derive


where is given by a sum of terms,


and the remainder is given by


The sum corresponds to the sum of the first terms in the asymptotic expansion of the exponential integral 22 ().

By replacing by the sum for low values of in the expressions for and we can derive series expansions in inverse powers of . For we find


At the sum of the first four terms shown agrees with the exact value to four decimal places. For we find by the same method


At the first two terms of the expansion in Eq. (7.6) yield a value accurate to one percent and at the value is accurate to five decimal places.

The eigenvector corresponding to the maximum eigenvalue can also be found by the same numerical method. For example, at we find the eigenvector for the modes of Sec. III. This shows that for this scale number all five modes contribute significantly to the eigenvector. The contributions from and cancel to a large extent. We find the ratios and . This shows again the importance of the Reynolds stress.

We note that for these modes of low order the whole calculation can be performed analytically, though at the expense of quite complicated expressions. It follows from the expressions in Appendix B that the elements of the matrix remain finite in the limit and apparently there is a cancellation from the matrix leading to the limiting behavior shown in Figs. 1 and 2.

Viii Simple model

It is worthwhile to discuss the properties of the simple model introduced at the beginning of Sec. VI in some more detail. We denote the coefficients of the zero frequency modes of Ref. 18 as . These are the amplitudes of the dipolar mode, the stresslet, and the quadrupolar mode, respectively. By comparison of the amplitudes of vector spherical harmonics on the spherical surface we find that the coefficients are related to the coefficients of the modes of Sec. III by the relations


These coefficients give the same surface displacement as the zero frequency modes for the set . It follows from Eq. (8.1) that if one considers surface displacements characterized by a chosen set of coefficients , then the coefficient of the boundary layer mode grows in absolute magnitude beyond all bounds as the scale number increases, whereas the amplitudes of the potential modes remain bounded. This causes an increase of dissipation as increases, owing to steep gradients in the boundary layer. Hence the efficiency decreases with increasing . As an example we show in Fig. 3 the ratio


where is calculated from the vector given below Eq. (6.3) by use of Eq. (8.1). The ratio starts at the optimal value in the limit , but eventually decreases to zero at large . In Fig. 4 we show the surface deformation at four chosen times.

The simplest way of constructing a good mode of swimming for large is to avoid the boundary layer altogether and use surface displacements corresponding to the potential mode given below Eq. (6.8). In this case , independent of . This is a mode of high efficiency, only slightly less efficient than the optimal mode given by the maximum eigenvector of Eq. (3.12).

Ix Discussion

The analysis shows that for given surface distortion the mean swimming velocity and mean power of a sphere depend in an intricate way on the dimensionless scale number , defined in Eq. (3.13) and related to the Roshko number by . Optimization of the mean swimming velocity for given power at fixed leads to a generalized eigenvalue problem. The eigenvector with largest eigenvalue within a class of strokes characterizes the optimal stroke in that class. Explicit expressions for the two Hermitian matrices characterizing the eigenvalue problem are given in Appendix B. The expressions are the pinnacle of the present analysis.

The results are of particular interest for large values of the scale number. It turns out that in this range it is crucial to take the Reynolds stress into account. Apparently on time average over a period the effect of pressure gradients is largely canceled by the effect of Reynolds stress. Quantitatively the effect is measured by a comparison of the contributions and to the mean swimming velocity. A numerical example is given at the end of Sec. VII.

In the asymptotic range of very large scale number we find that the optimal efficiency tends to the value for potential flow for the class of axisymmetric strokes involving the three modes of order , as well as for axisymmetric strokes involving the five modes of order . The same behavior was found for the swimming of a deformable slab 15 (). It may be conjectured that for a sphere this behavior occurs also for more complicated strokes. It is shown in Figs. 1 and 2 that in the asymptotic range the swimming is somewhat less effective than in the Stokes limit. There is an appreciable change in the mode of optimal swimming as is seen by a comparison of the eigenvectors of displacements at the end of Secs. VI and VII.

For the strokes considered the optimal efficiency shows a maximum in the intermediate range of scale number. In this regime the optimal efficiency varies little with the scale number. For a swimmer the actual value of the scale number is determined by its geometrical dimension and the kinematic viscosity of the fluid. For that particular value of the swimmer can optimize its stroke to second order in the amplitude from the solution of the eigenvalue problem. It will be of interest to determine how the efficiency varies as the amplitude of stroke is increased. This requires numerical solution of the Navier-Stokes equations with no-slip boundary condition for the moving surface and is beyond the scope of the present investigation.

Appendix A

In this Appendix we provide expressions for angular integrals which occur in the calculation. We consider vector-functions of the form


and scalar functions of the form


In the calculation of the matrix we need the angular integral


where in the last equation we have omitted the dependence on for brevity. Similarly


In the calculation of the matrix we need the angular integral


as well as the integral