Regularized Stokeslet rings - an efficient method for axisymmetric Stokes flow, with application to the growing pollen tube
James Tyrrell, David J Smith, Rosemary J Dyson
TL;DR
The paper introduces the regularized Stokeslet ring method, an efficient axisymmetric flow solver that improves accuracy and convergence for biological fluid mechanics applications, demonstrated through various examples including pollen tube flow.
Contribution
It derives the regularized ringlet solution for axisymmetric Stokes flow, enabling more accurate and efficient modeling of biological systems compared to existing methods.
Findings
Regularized ringlets provide accurate velocity and drag calculations.
The method shows better convergence than singular solutions.
Application to pollen tube flow demonstrates biological relevance.
Abstract
The method of regularized Stokeslets, based on the divergence-free exact solution to the equations of highly viscous flow due to a spatially-smoothed concentrated force, is widely employed in biological fluid mechanics. Many problems of interest are axisymmetric, motivating the study of the azimuthally-integrated form of the Stokeslet which physically corresponds to a ring of smoothed forces. The regularized fundamental solution for the velocity (single layer potential) and stress (double layer potential) due to an axisymmetric ring of smoothed point forces, the `regularized ringlet', is derived in terms of complete elliptic integrals of the first and second kind. The relative errors in the total drag and surrounding fluid velocity for the resistance problem on the translating, rotating unit sphere, as well as the condition number of the underlying resistance matrix, are calculated; the…
| 51 | 101 | 201 | 401 | ||
|---|---|---|---|---|---|
| 0.01 | |||||
| 0.005 | |||||
| 0.001 |
| 51 | 101 | 201 | 401 | ||
|---|---|---|---|---|---|
| 0.01 | |||||
| 0.005 | |||||
| 0.001 |
| 51 | 101 | 201 | 401 | ||
|---|---|---|---|---|---|
| 0.01 | |||||
| 0.005 | |||||
| 0.001 |
| 51 | 101 | 201 | 401 | ||
|---|---|---|---|---|---|
| 0.01 | |||||
| 0.005 | |||||
| 0.001 |
| Cond | ||||
|---|---|---|---|---|
| Cond | ||||
|---|---|---|---|---|
| number of nodes | ||||||
|---|---|---|---|---|---|---|
| number of nodes | ||||||
|---|---|---|---|---|---|---|
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
MnLargeSymbols’164 MnLargeSymbols’171
Regularized Stokeslet rings – an efficient method for axisymmetric Stokes flow, with application to the growing pollen tube
J. Tyrrell
D.J. Smith
R.J. Dyson
University of Birmingham, School of Mathematics, Edgbaston, Birmingham, B15 2TT, United Kingdom
Abstract
The method of regularized Stokeslets, based on the divergence-free exact solution to the equations of highly viscous flow due to a spatially-smoothed concentrated force, is widely employed in biological fluid mechanics. Many problems of interest are axisymmetric, motivating the study of the azimuthally-integrated form of the Stokeslet which physically corresponds to a ring of smoothed forces. The regularized fundamental solution for the velocity (single layer potential) and stress (double layer potential) due to an axisymmetric ring of smoothed point forces, the ‘regularized ringlet’, is derived in terms of complete elliptic integrals of the first and second kind. The relative errors in the total drag and surrounding fluid velocity for the resistance problem on the translating, rotating unit sphere, as well as the condition number of the underlying resistance matrix, are calculated; the regularized method is also compared to 3D regularized Stokeslets, and the singular method of fundamental solutions. The velocity of Purcell’s toroidal swimmer is calculated; regularized ringlets enable accurate evaluation of surface forces and propulsion speeds for non-slender tori. The benefits of regularization are illustrated by a model of the internal cytosolic fluid velocity profile in the rapidly-growing pollen tube. Actomyosin transport of vesicles in the tube is modelled using forces immersed in the fluid, from which it is found that transport along the central actin bundle is essential for experimentally-observed flow speeds to be attained. The effect of tube growth speed on the internal cytosolic velocity is also considered. For axisymmetric problems, the regularized ringlet method exhibits a comparable accuracy to the method of fundamental solutions whilst also allowing for the placement of forces inside of the fluid domain and having more satisfactory convergence properties.
I Introduction
The Stokes equations for incompressible flow at zero Reynolds number are used extensively to model the viscous-dominated regime of microscale flows, particularly biological flows associated with cilia-driven transport, and the motility and feeding of flagellated cells such as bacteria, spermatozoa, algae and choanoflagellates. For an overview, see Lauga and Powers Lauga and Powers (2009). The fundamental solution of the Stokes flow equation, which corresponds to the flow driven by a single spatially-concentrated force is often referred to as the Oseen tensor or Stokeslet. The linearity of the Stokes flow equations enables the construction of solutions to problems involving moving boundaries with complex geometry through integral sums of Stokeslets, forming the basis for the method of fundamental solutions, slender body theory and boundary integral methods. The latter numerical method has the principle major advantage of avoiding the need to mesh the fluid volume, which has enabled highly accurate and efficient simulation of biological flow systems for several decades Phan-Thien et al. (1987); Ramia et al. (1993); Shum et al. (2010). Indeed more approximate methods based on line distributions of Stokeslets and higher order singularities also enabled major progress in this area before the present era of computationally-intensive research. For review see the earlier work of Chwang, Wu and co-authors Chwang and Wu (1975), who also explore a wide range of applications as part of a series of papers on low-Reynolds number flow Chwang and Wu (1974); Chwang (1975); Chwang and Wu (1976); Johnson and Wu (1979); Huang and Chwang (1986).
Nevertheless, two implementational issues arise with methods based on singular solutions. The first is that boundary integrals of solutions with a -type singularity can be technically complex to evaluate on or near the boundary. Line integrals associated with models of slender bodies such as cilia and flagellar are ‘more singular’, and can require careful distinction between the inside and outside of the body. Moreover, there are cases in which immersed forces due to, e.g. many suspended moving particles, are desired to be modelled by an immersed volumetric force. Cortez et al. developed the method of regularized Stokeslets Cortez (2001); Cortez et al. (2005) based on the exact, divergence-free solution to the Stokes flow equations due to a concentrated but spatially-smoothed (regularized) force. This approach has enabled the use of Stokeslet methods in a wider range of applications, such as those in which an inducing force is present in the interior of the fluid domain (as either a point in or a point/curve in ).
While conceptually elegant, the standard implementation of the method of regularized Stokeslets is computationally expensive, motivating the development of boundary element discretization Smith (2009), line integration Cortez (2018), and meshless interpolation Smith (2018) among other approaches. Many diverse biological flow problems of interest exhibit rotational symmetry, examples including spherical ‘squirmer’ swimmers Blake (1971), the conceptual toroidal swimmer of Purcell Purcell (1977), and cytosolic flow in elongating pollen tubes Chebli et al. (2013). Thus motivated, in this paper we study axisymmetric Stokes flows in which the singular and regularized Stokeslets can be integrated azimuthally to yield an axisymmetric ring of point forces. The singular solution to this problem is already known Pozrikidis (1992); the regularized solution, which we term the ‘regularized ringlet’, is newly derived. This solution forms the basis for an efficient axisymmetric method of regularized Stokeslets.
We begin in Section II by introducing the singular and regularized Stokeslets, and review their application in solving the resistance problem for a rigid body translating in a viscous fluid. In Section II.1 we derive the regularized ringlet, whose analytical solution is given in Section II.2 and is applied to the resistance problem for the translating and rotating sphere in Sections III.1 – III.2. The double layer potential (DLP), relevant to bodies undergoing volume-changing deformation, is considered in Section II.3, with analytical evaluation of the azimuthal integral in the DLP given in Appendix D. In Section III.3 we consider the case of Purcell’s toroidal swimmer Purcell (1977); Leshansky and Kenneth (2008), in which the method of regularized ringlets enables the calculation of propulsion speeds which are in excellent agreement with analytical results for both slender and non-slender tori. We then present a study of fluid flow in the angiosperm pollen tube in Section IV, an illustrative example highlighting the benefits of being able to apply smoothed forces inside the fluid (two prominent features of regularization). In Section V we conclude with a summary of results and a discussion of future work.
II Singular and regularized Stokeslet solutions
For the viscous–dominated very low Reynolds number flow associated with microscopic length scales and slow velocities, incompressible Newtonian flow is well–approximated by the steady Stokes flow equations
[TABLE]
where is dynamic viscosity, the pressure, the velocity, and the applied force per unit volume. In the case of a singular force of the form for arbitrary point force , arbitrary point at which the singularity is located, and where is the Dirac delta function, the fundamental solution Fairweather and Karageorghis (1998); Young et al. (2006) for is given by
[TABLE]
where
[TABLE]
is known as the Stokeslet.
The singularity in the Stokeslet solution can be eliminated without loss of incompressibility by regularization of the force , as described by Cortez and colleagues Cortez (2001); Cortez et al. (2005). The Dirac delta function is replaced with where is a radially symmetric, smooth ‘cutoff’ function with the property . This is in essence applying the force over a small ball, varying smoothly from a maximum at its centre to sufficiently far away, instead of using an infinite point force as in the classical Stokeslet solution. The numerical parameter dictates the radius of support of the force, and as the classical solution is recovered. Solutions for using regularized Stokeslets differ from those found using the singular Stokeslet only near the point where the force is applied. Following Cortez et al. Cortez et al. (2005) we take
[TABLE]
which yields
[TABLE]
By considering a solid body moving through the fluid, it can be shown that
[TABLE]
where is the force per unit area exerted by the body surface (denoted ) on the surrounding fluid Cortez et al. (2005). The above equation is exact; replacing the left hand side with the velocity introduces an error where on or near the body surface, and sufficiently far away.
Discretising Equation (7) using Stokeslets on the surface of the solid body enables the approximation of the fluid velocity at any point via a numerical quadrature formula
[TABLE]
where denotes the th component of the force per unit area applied at the point (a Stokeslet location) and is the quadrature weight associated with the th particle. The value of is dependent on the geometry of the body surface and in Cortez and colleagues’ work has units of area.
II.1 Derivation of the regularized ringlet
Consider a specific case of Equation (7) using a cylindrical coordinate system in which the body exhibits rotational symmetry about the axis. This symmetry enables analytical integration azimuthally, reducing the surface discretisation to a line discretisation and increasing accuracy. In doing so, we are effectively placing ‘rings’ of regularized Stokeslets at positions for and (see Figure 1 for a diagram of a single ring). This is analogous to covering the surface of the body in ‘strips’ instead of the patches used in a standard D Cartesian discretisation. With a surface parameterization where denotes arclength and , the boundary integral equation (7) reads
[TABLE]
Converting to cylindrical polar coordinates, we introduce the transformation matrix
[TABLE]
Letting indices and correspond to Cartesian and cylindrical polar bases respectively (with Einstein summation convention employed for both sets), it follows that . Assuming that velocity is evaluated at fluid point in cylindrical polars, it further follows that Recognizing that , substitution of the cylindrical forms of the velocity and force vectors into Equation (9) thus yields
[TABLE]
Under the assumption of axisymmetry, it is sufficient to only consider cases . This results in the Cartesian and the cylindrical polar planes coinciding, such that the transformation matrix with simply yields the identity matrix and Equation (11) reduces to
[TABLE]
in which the ‘ringlet’ kernel
[TABLE]
is implicitly defined. Here, is the Kronecker delta function in which if and [math] otherwise. Unlike the Stokeslet which has the symmetric property , ordering of arguments in the ringlet is important; the first and second arguments in denote the fluid point and ring location respectively, with the ring radius being the crucial non-symmetric term.
The various terms in the Stokeslet can be evaluated in cylindrical polar coordinates via
[TABLE]
The resulting form of each is given in Appendix A and used in the evaluation of in Equation (14). These integrals yield the regularized fundamental solution for an axisymmetric ring of concentrated force (the regularized ringlet). Analytical evaluation reveals that such that the rotational problem for decouples from and . The remaining nonzero yield the equations
[TABLE]
and
[TABLE]
Utilising Equations (19) – (20) in tandem models axisymmetric problems with or without azimuthal rotation, in which the fluid experiences a constant force in each principal direction at points along which the ringlet is located. This could be used, for example, to model the flow around a mobile axisymmetric body rotating about an axis defined by its direction of translation.
Approximating the integrals in Equations (19) – (20) numerically using a series of rings yields a system of equations of the form
[TABLE]
where are the radial, azimuthal, and axial components of the forces per unit area applied at ringlet location respectively, and is the quadrature weight associated with for numerical integration over . The quantity has units of length unlike its counterpart in Cortez and colleagues’ work (units of length squared). It is also possible to recompile the force per unit area and quadrature weight into a force per unit length , such that Equation (27) can alternatively represent the fluid velocity induced by a series of rings.
By considering the fluid velocity at each individual ringlet location, an invertible system can be produced. In the irrotational case (), this takes the block matrix form
[TABLE]
where
[TABLE]
and
[TABLE]
Hence the forces needed to induce a given prescribed velocity in the fluid may be found (the resistance problem). This works for both a series of translating rings and a translating axisymmetric body, with the force per unit length essentially ‘absorbing’ both the force per unit area and the quadrature weight in the latter case. Inclusion of rotation involves the formulation of a similar invertible system for Equation (28) which can be solved separately.
II.2 Analytical evaluation of the regularized ringlet
The nonzero elements of the regularized ringlet can be expressed in the form
[TABLE]
in which and
[TABLE]
with . Equation (38) is found by using the double angle formula for as well as symmetry arguments about . Following the example of Pozrikidis Pozrikidis (1992), the integrals can be computed by first expanding the numerator of the integrand in Equation (38) to obtain a series of polynomial integrals with respect to . Letting
[TABLE]
it follows that
[TABLE]
The individual integrals can be expressed in terms of complete elliptic integrals of the first and second kind, which are respectively defined
[TABLE]
The solutions for each (as can be found in Section 2.58 of Gradshteyn and Ryzhik Gradshteyn and Ryzhik (2014)) are given by
[TABLE]
from which it follows that
[TABLE]
Substitution of Equations (45) – (47) back into Equations (32) – (36) yields the complete solution for the regularized ringlet
[TABLE]
The solutions given by Equations (48) – (52) can be readily evaluated except when or . In the limit as , all . In the limit as , both and tend to finite values while This behaviour is described in detail in Appendix B. Streamlines for the flow induced by the regularized ringlet with associated unit forces in both and directions are given in Appendix C.
The form of the ringlet solutions is similar to those for the ring of singular Stokeslets (as detailed by Pozrikidis in Pozrikidis (1992) with the exception of the newly derived )
[TABLE]
and in the limit as our solutions are equivalent to their singular counterparts. This can be verified by substitution of into and is a result of the cutoff function approaching a delta distribution as .
Equations (19) – (20) and (48) – (52) provide the exact analytical solution for the fluid velocity at any point due to the drag force per unit area on the surface of a generalized axisymmetric body. Using Equations (27) – (28) in place of (19) – (20) yields the numerical solution based on discretization over the arclength . In the case of a single ring, removing the integral over and replacing the force per unit area with a force per unit length yields the exact analytical solution for the fluid velocity induced by the force acting along the ring in D space.
II.3 On the double layer potential
A more complete formulation of Equation (7) for the fluid velocity induced by a translating body is given by
[TABLE]
where the first and second integrals on the right hand side are known as the single layer potential (SLP) and the double layer potential (DLP) respectively. The stress tensor present in the DLP is given by
[TABLE]
The DLP can be neglected for problems in which the condition is satisfied (Pozrikidis, 1992, Chapter 2.3), which is the case throughout this manuscript. However we provide the expressions for the DLP in Appendix D, which may be of value for future studies of systems in which the condition is violated (e.g. bubbles).
III Simple examples and test cases
In Section II the expression for the regularized fundamental solution for an axisymmetric ring of concentrated forces, the regularized ringlet, was derived. In the following, we demonstrate the validity of the method through application to simple cases of motion.
The first two cases concern the translation and rotation of the unit sphere in a Stokesian fluid, treated independently in Sections III.1 and III.2 respectively. In Section III.3 a more complicated example is considered: the propulsion of ‘Purcell’s toroidal swimmer’ Purcell (1977); Leshansky and Kenneth (2008), powered by tank treading of the torus surface. In considering these different cases, it is shown that the method of regularized ringlets can be used to model the surface motion of axisymmetric bodies in each principal direction in a cylindrical coordinate system.
III.1 Translating unit sphere
The validity of the regularized ringlet method is illustrated by solving the resistance problem for the translating unit sphere. Given a prescribed surface velocity (), Equation (13) yields a Fredholm first kind integral equation for the unknown force distribution Smith (2018). The method of regularized ringlets (implemented here via Matlab) can be used to solve this problem.
The sphere is parametrized in the plane by for , then discretized as
[TABLE]
The velocity boundary condition is prescribed at each , and the resulting linear system is solved to yield the required force densities at each of these locations. The drag is then calculated as
[TABLE]
which is compared with the Stokes law value of . The relative errors are given for various values of and regularization parameter in Table 1 alongside the condition number of the resistance matrix in Table 2.
For given , excessively small results in the drag error becoming non–monotonic. For given , increasing eventually ceases to result in a further reduction in the relative error. This is often the case with regularized Stokeslet methods (eg. see Figure 15, Figure 16, and references Cortez et al. Cortez et al. (2005), Gallagher et al. Gallagher et al. (2018)).
A thorough comparison of our results for the translating unit sphere with those of Cortez et al. Cortez et al. (2005) can be found in Appendix E, in which it is found that using just ringlets yields consistently more accurate solutions than using regularized Stokeslets in a standard D discretization of the sphere surface. A comparison to results obtainable using the axisymmetric method of fundamental solutions (a singular Stokeslet method) is given in Appendix F, in which the relative error in the fluid velocity is also discussed. In general, although the singular method can be tuned to give smaller relative errors in either the fluid velocity or total drag separately, it cannot do so simultaneously; regularized ringlets display more satisfactory convergence properties and are the more effective method to minimize errors in both fluid velocity and total drag.
III.2 Rotating unit sphere
The solution for the steady motion of a Stokesian fluid surrounding a solid sphere rotating uniformly about a central axis is well known and can be found in e.g. Hydrodynamics by Lamb (Lamb, 1945, Chapter XI). If the sphere rotates around its axis in an cylindrical coordinate system with an angular velocity , the resulting angular velocity of the fluid is given by where . This can be written in terms of the linear velocity (more readily usable in the Stokeslet formulae) over the entire domain as,
[TABLE]
where corresponds to the surrounding fluid velocity and to the solid body rotation of the sphere respectively. The zero Reynolds number torque on this sphere is given by,
[TABLE]
the derivation of which can be found in Rubinow and Keller (1961). This torque is associated with a drag force per unit area on the surface of the sphere given by , as detailed in Appendix H. The force per unit length used in the method of regularized ringlets thus takes form
[TABLE]
although for the resistance problem this is not prescribed.
The sphere surface is again parametrized in the plane by for , discretized using ringlets at locations . Letting , the velocity is prescribed at each and the resistance matrix is constructed to yield the required force densities at each of these locations. The torque is then calculated as
[TABLE]
in which is squared since the torque is a moment, the product of distance and force. Comparing to the value , relative errors are given for various values of and regularization parameter in Table 3 alongside the condition number of the resistance matrix in Table 4.
These results are similar to those of the translating unit sphere, although for given the relative error in the drag is generally slightly larger and the condition number of the resistance matrix slightly smaller than for the results using the same on the translating unit sphere.
III.3 Purcell’s toroidal swimmer
The torus is the simplest geometry capable of describing self–propelled organisms Thaokar et al. (2007). Purcell’s toroidal swimmer Purcell (1977); Taylor (1952) describes one such organism, the geometry of which can be seen in Figure 2. Inward rotation of the torus surface produces a net force in the direction of motion of the outermost surface (against which the torus is propelled). The magnitude of this net force (and resultant propulsion speed of the torus) is dependent on the speed with which the surface of the torus rotates as well as the slenderness of the torus. It has been suggested that this mechanism could describe how a DNA miniplasmid could be turned into a self–propelled nanomachine Kulić et al. (2005). Three modes of locomotion are considered by Leshansky and Kenneth Leshansky and Kenneth (2008), corresponding to tank treading of (a) an incompressible surface, in which the tangential surface velocity is largest on the inner surface, (b) a weakly compressible surface, in which the tangential surface velocity is constant, (c) a highly compressible surface, in which the tangential surface velocity is largest on the outer surface. In the following, we restrict ourselves to looking at the case of constant tangential surface velocity.
The torus geometry is reduced to a single slenderness parameter , where and refer to the major and minor radii of the torus respectively. The torus surface is parametrized in terms of angle in the plane such that and corresponds to the outermost radial point on the surface of the torus, traversed in an anticlockwise direction. For the free swimming torus, the rigid body translation and rotation of the torus surface are related to the force per unit area exerted by the torus on the surrounding fluid by,
[TABLE]
subject to the condition of zero net force in the direction,
[TABLE]
where is the arclength parameterization of the cross section of the torus surface in the plane. We note that the additional free swimming conditions of zero net force in and zero total moment (as outlined by Phan–Thien et al. Phan-Thien et al. (1987)) are automatically satisfied by axisymmetry and respectively.
The propulsion speed of the rotating torus for any given value of and rotation speed can be determined by considering two separate situations: one in which motion is purely translational in the direction (the toroidal glider with ) and one in which motion is purely rotational (the anchored toroidal pump with where is the unit vector whose direction varies over , pointing tangential to the surface in the anticlockwise direction at all points). The glider and pump have associated force distributions and respectively. Using the regularized ringlet, surface motions can be prescribed (Figures 3a & 4a) in order to compute the associated force distributions (Figures 3b & 4b) responsible for producing each motion. For the toroidal glider the radial force is assumed to be zero (), hence
[TABLE]
which, after solving for the unknown force distribution , yields the net axial force in the direction
[TABLE]
which is nonzero. For the anchored toroidal pump the surface velocity is given by in coordinates. It follows that
[TABLE]
from which it can be determined that the net axial force is given by
[TABLE]
By linearity of the Stokes flow equations we may subtract the gliding solution from the pump solution to rewrite the system in the form
[TABLE]
where and the net force is equal to
[TABLE]
as required for the free–swimmer. The propulsion velocity of the swimming torus is thus given by , opposing the direction of outer surface motion. The speed is dependent on both the rotation speed and the slenderness ratio , so both and must be recomputed whenever one of these parameters is changed. However, Leshansky and Kenneth Leshansky and Kenneth (2008) were able to show that the propulsion speed depends linearly on the rotation speed, such that by considering the scaled propulsion speed it is only necessary vary to be able to consider all possible propulsion speeds resulting from a given constant rotational surface velocity.
Our results for the scaled propulsion velocity found using a discretization of the torus surface using regularized ringlets are compared with those obtained by Leshansky and Kenneth Leshansky and Kenneth (2008), who tackled the same problem using a line distribution of rotlets at the torus centreline (inaccurate as ) and an exact series solution via expansion in toroidal harmonics (Figure 5). Of particular note, in the limit as it is found that using rings with in the regularized ringlet method yields a scaled propulsion velocity of , representing just a error when compared to the series solution value of This error can be reduced to by using rings, at which point the value of the scaled propulsion velocity as calculated by the method of regularized ringlets is . This is a significant improvement over the solution found using a line distribution of rotlets, in which the error is for all .
In addition to giving values for the scaled propulsion speed that are in excellent agreement with the series solution of Leshansky and Kenneth Leshansky and Kenneth (2008), the regularized ringlet solution also provides the force required at each point on the torus surface to produce the swimming motion. The series and rotlet solutions do not yield this information, with the propulsion speed instead being calculated according to the net drag force on the toroidal glider in these methods. This, in combination with Figures 3a & 3b, highlights why the centreline rotlet solution is inaccurate for small values of ; as slenderness decreases, the difference between the magnitude of the force required at the inner- and outer-most surfaces to produce rigid body translation of the torus grows large. In Figure 3b (in which ), the drag force on the outer surface is approximately times larger than that on the inner surface. The constant centreline force associated with the rotlet cannot account for this discrepancy, whereas the full discretization of the torus surface using regularized ringlets can.
Figures 6a & 6b show the streamlines and magnitude of the fluid velocity in a region near the force–free toroidal swimmer with slenderness parameter , undergoing uniform counterclockwise surface rotation with unit angular velocity. This results in propagation of the torus in the direction . Fluid passing through the central hole of the torus is caught in closed streamlines, in agreement with the results of Leshansky and Kenneth Leshansky and Kenneth (2008).
IV Cytosolic flow in the pollen tube - an illustrative example of the benefits of regularization
In Section III, the method of regularized ringlets was applied to simple cases of motion to illustrate validity and applicability. These canonical examples demonstrate the use of ringlets in situations with moving boundaries. However, perhaps the most important feature of the regularized method is the ability to place source points directly in the fluid, in which case the regularized ring represents a curve in D space. The regularization parameter can then be used to control the spreading of the force (the size of the region over which it is applied). This is useful in situations in which the exact location at which a force is applied in a fluid is unknown or when a force is applied over a large area. In the following, these features of regularization are explored in more detail by considering the example of fluid flow in the angiosperm pollen tube. A more extensive study of the biomechanics of pollen tube growth (with a particular focus on the transport and distribution of elements of the cytoplasm) will follow in a future paper.
IV.1 A brief overview of the pollen tube
The pollen tube, a cellular protuberance originating from the pollen grain, is a vital component of the fertilisation process in plants. Responsible for the delivery of sperm cells from the pollen grain to the ovule, the tube may have to grow over cm in length at speeds in excess of Bove et al. (2008); Booy et al. (1992). This rapid elongation of the tube occurs via tip growth, with expansion localised close to the apical region of the cell Fayant et al. (2010). To prevent rupture, new cell membrane and cell wall material must be targeted to the sites of highest expansion. Spherical vesicles in the cytosolic region act as the delivery vectors for this material Chebli et al. (2013) and are carried by motor myosin proteins along actin filaments in a process referred to as actomyosin transport. These actin filaments, arranged in bundles in the periphery and centre of the tube, are oriented in such a way that ‘full’ peripheral vesicles travel towards the growing apex where they secrete their contents into the wall. In order to maintain the correct ratio of structural components, ‘empty’ vesicles (comprised mostly of membrane) are also secreted by the wall back into the cytoplasm, where it is hypothesised that they are picked up by myosin on the central actin bundle and travel away from the apex. The combined movement of these vesicles induces a flow in the cytosol, known as ‘cytoplasmic streaming’ or ‘cyclosis’, which further aids in the cycling of vesicles towards and away from the apical region. To our knowledge, no comprehensive attempt has been made to produce a complete model of cytosolic flow in the pollen tube in a manner deriving directly from physical principles. In this section, we show how this is easily achievable using the method of regularized ringlets.
IV.2 Mathematical model
IV.2.1 Geometry and boundary conditions
The typical geometry of the tube (in D cylindrical coordinates) along with the two vesicle populations and actin bundles can be seen in Figure 7, in which the static cylindrical shank is joined to a growing apical hemispherical cap. The boundary of the domain is split into four sections: the symmetry boundary running down the centreline of the tube, the static, impermeable peripheral wall , the growing apical hemisphere and the artificial basal boundary (where the computational domain is truncated).
The Reynolds number for fluid flow in the pollen tube can be evaluated using the typical flow speed (1\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}), tube radius ($\approx$8.13\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}) and the kinematic viscosity of water ({10}^{6}\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2},\mathrm{s}^{-1}) to find $\text{Re}\approx${10}^{-5}\text{\,}, firmly in the regime of Stokes flow. Although pollen tube growth is typically oscillatory, acceleration and deceleration are small compared to growth speed itself. This enables the use of the steady Stokes equations.
Since the cytosolic flow is induced by the actomyosin transport of vesicles along cytoskeletal actin bundles, forces must be applied inside the computational domain. Singular stokeslet methods typically do not allow for this unless the surfaces of individual vesicles are discretized, but their small size (100\text{,}\mathrm{n}\mathrm{m}$$ radius) compared to the typical length scale of the problem and their large number density make this impractical. One solution is to use the method of regularized ringlets instead, placing rings in series along the centreline of the peripheral actin bundle with the parameter being used to control the bundle thickness.
The boundary conditions for the fluid velocity are given by,
[TABLE]
in accordance with the assumptions of tip growth and axisymmetry. No restriction is placed on on the artificial boundary . For the growth velocity in the hemispherical apex, the normal displacement growth assumption of Dumais Dumais et al. (2006) is employed to define
[TABLE]
where is the angle between the outward-pointing surface normal and the positive axis (varying from [math] at the point where the hemisphere joins the shank to at the extreme apex) and is the growth speed of the tube. An example of how the growth of the boundary varies over the apical hemisphere can be seen in Figure 8, where the maximum growth speed is equal to m s*-1* and the wall velocity in the adjacent shank (not pictured) is equal to zero.
IV.2.2 Actin bundle model
A section of the pollen tube is modelled ranging from the extreme apex to a point equal to six tube radii distal where the domain is artificially truncated. In dimensionless values, the extreme apex is thus given by with being the distal truncation line. The central line of axisymmetry is given by , with denoting the peripheral boundary in the shank. The hemispherical apex is the upper-right quarter circle of radius , centred at . The peripheral actin bundle (PAB) is considered adjacent to the pollen tube wall, with its width equal to one fifth of the pollen tube radius () and extending at an angle into the apical hemisphere (see Figure 9) in accordance with the confocal microscopy imaging of Lovy-Wheeler et al. Lovy-Wheeler et al. (2005). The central actin bundle (CAB) is located further away from the apex than the PAB, and is modelled with a reduced thickness (half that of the PAB) in cylindrical coordinates to account for the central axisymmetry of the tube.
By placing rings of regularized Stokeslets in series along the centre line of the PAB, it is possible to carefully select an appropriate value for the regularization parameter such that the region over which the majority of the force distribution is applied is roughly the same as the bundle thickness. The same procedure is employed for the CAB, only with individual regularized Stokeslets rather than ringlets since the centreline coincides with . Using it is found that of the total applied force is contained within a region of radius which corresponds well to the approximate thickness of the actin bundles. This is shown in Figure 10 where the maximal value of has been scaled to for the sake of clarity. Smaller values of result in an even larger percentage of the total force being contained within but are increasingly skewed towards .
We note that our modeling ignores the presence of a third f-actin structure, the short actin bundles observed in the extreme apex Fu et al. (2001); these bundles are transient and significantly less dense than the thick peripheral and central actin bundles. Further, due to their proximity to the growing tube wall, any effect that these short actin bundles may have on the fluid velocity is likely insignificant compared to the effect of the growth of the wall.
IV.2.3 Frames of reference
In each of the velocity profiles to be produced, cytosolic velocities are prescribed along the actin bundles in accordance with values taken from the STICS analysis of Bove et al. Bove et al. (2008). This is done under the assumption that the large-scale directed movement of vesicles must be a consequence of cytosolic flow (induced by the smaller-scale movement of individual vesicles along actin). The fluid velocity is defined in terms of two different coordinate systems: the static lab frame where and the moving tip frame with . The lab frame represents the fluid velocity relative to a stationary observation point (in which on due to boundary growth), whereas the tip frame represents the fluid velocity relative to the growing tip (resulting in a static domain with at all boundary points except on ). The steady geometry of the domain for the advancing tip frame means that satisfies the steady Stokes equations, from which it follows that is also a solution (since the two differ by the uniform flow field only).
IV.2.4 Two model scenarios for comparison
In Section IV.3.1, we investigate the role of the CAB in vesicle transport. Fluid velocity is prescribed on the PAB only (see Figure 11a) in order to compare the resultant central fluid velocity to experimental results. The centreline of the PAB is modelled using the union of the straight line extending from to and the curve for . The regularized ringlet placement is extended to to ensure the velocity profile at is consistent with the rest of the tube. The fluid velocity on the straight line is given by with the fluid velocity on the curve being tangential and of constant magnitude, that is, for . On the peripheral wall (again extended to ), the velocity is in the shank and prescribed according to the growth speed of the tube in the apical hemisphere using the velocity function for .
In Section IV.3.2, the velocity of the growing boundary is varied in order to determine the effect of the growth speed of the tube on the cytosolic velocity. Fluid velocity is prescribed on each of the PAB, the CAB, and the growing boundary accordingly (see Figure 11b). Fluid velocity along the PAB is kept constant at based on the observation that cytoplasmic streaming rates are typically independent of pollen tube growth rates Vidali et al. (2001). Fluid velocity on the centreline of the central actin bundle, given by for , is also prescribed. It is assumed that vesicles are sufficiently closely packed on this bundle that the fluid velocity at its centre can be approximated by the speed of vesicles on actin (), giving . Since these node locations are at , standard regularized Stokeslets must be used here.
A regularisation parameter of is chosen in all cases and the ringlets are linearly spaced a distance approximately 0.025 apart, resulting in a smooth velocity profile. The flow velocity for (not pictured) always matches the flow velocity at almost exactly, with no further change occurring in the direction.
IV.3 Results
IV.3.1 Role of the central actin bundle
The CAB has long been hypothesized to aid in the removal of vesicles from the apical region, but direct observation of vesicle transport along actin is often hindered by the small size of vesicles (typically below the resolution limit of conventional confocal microscopes). Higher resolution imaging methods such as evanescent wave microscopy have been used to observe long-range vesicle movement (presumably a result of actomyosin transport) in the periphery of the tube Wang et al. (2006), but the limited penetration depths available in these methods (400\text{,}\mathrm{n}\mathrm{m}) do not allow for imaging of the central region. STICS analysis suggests a way of determining whether the central bundle participates in vesicle transport. Directed vesicle movement in the periphery of the tube is not seen to exceed speeds of $\approx$0.5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\,\mathrm{s}^{-1}, smaller than the 0.8\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$ observed in the centre Bove et al. (2008). By using the method of regularized ringlets and prescribing fluid velocity only along,
- •
the PAB, where 0.5\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$,
- •
the tube boundary with apical growth speed 0.1\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$,
the resulting fluid velocity in the centre should provide further insight insight into whether the central bundle participates in vesicle transport.
The results of this investigation can be seen in Figure 12, with both lab (Figure 12a) and tip (Figure 12b) frames shown. The shape of the velocity profile in Figure 12a is in excellent agreement with the STICS analysis of Bove et al. Bove et al. (2008), with a wider band of basal flow through the center than apical flow in the periphery. However, despite an arguably exaggerated prescribed peripheral fluid velocity of , fluid velocity in the centre does not achieve speeds of . This is a strong indication that the CAB must also participate in the transport of vesicles, particularly considering that that our current implementation of the Stokes equations does not account for variations in local fluid viscosity (known to be larger in the presence of filamentous actin networks Wong et al. (2004), reducing the fluid velocity induced by any given force). The increased flow speed through the centre of the tube can be easily accounted for by inclusion of additional drag from actomyosin vesicle transport (known to reach speeds of up to Vidali et al. (2001); De Win (1997)) along the CAB, with the largest velocities being observed at the very centre as a result of the reduced cytosolic volume in this region.
IV.3.2 The influence of growth speed
We now turn our attention to the effect of tube growth speed on the cytosolic velocity profile. A further three new velocity profiles are produced, based on three different growth speeds for the tube (\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$). These speeds (approximately) correspond to that of a static tube, the typical growth rate cited in Bove et al. Bove et al. (2008), and the average growth rate measured by Vidali et al. Vidali et al. (2001) for the Lilium longiflorum species. Fluid velocity is prescribed along
- •
the PAB (0.5\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$),
- •
the CAB (1\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$),
- •
and the tube boundary (0.2\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m},\mathrm{s}^{-1}$$),
in each case.
Figures 13a – 13c show the velocity profiles for each of these three tube growth speeds in the lab frame. Significant differences can be seen between the three profiles, both in the magnitude and direction of the fluid velocity. In particular, at larger growth speeds there is a wider band of cytosolic flow in the positive direction in the peripheral region, and the central band of basal flow is both narrower and of a reduced magnitude. This is an expected consequence of mass conservation, since in a tube with a faster growth speed more fluid must flow towards the apical region to fill the increasing space.
Figures 13d – 13f show the velocity profiles for each of these three tube growth speeds in the tip frame. Here, the growth speed of the tube has been subtracted from the component of the fluid velocity for each corresponding velocity profile in the lab frame. Interestingly, the differences between the velocity profiles are significantly less pronounced in this frame of reference. Differences in the magnitude of the fluid velocity still persist, but the overall shape of the profiles bear a striking similarity. The persistent shape of the velocity profiles seen in Figures 13d – 13f could help explain the observed similarities in the distribution of the apical vesicle population and other elements of the cytoplasm across multiple pollen tube species and throughout different phases of oscillatory growth.
V Conclusions and future work
In this paper, the regularized fundamental solution for the velocity (single layer potential) and stress (double layer potential) due to an axisymmetric ring of smoothed point forces (the ‘regularized ringlet’) was derived, expanding on the work of Cortez and colleagues Cortez (2001); Cortez et al. (2005). The velocity solution, written in the form of complete elliptic integrals of the first and second kind, tends to the singular solution of Pozrikidis Pozrikidis (1992) in the limit as the regularization parameter tends to [math].
For the resistance problem on the translating and rotating unit sphere, the method of regularized ringlets was shown to produce small relative errors in drag calculations and small condition numbers for the underlying resistance matrices. Further testing in the case of the translating unit sphere shows that the regularized ringlet performs favourably compared to the traditional method of regularized Stokeslets both in terms of accuracy (Appendix E) and speed (Appendix G) and that it possesses more satisfactory convergence properties than the singular ringlet (Appendix F).
The applicability of the regularized ringlet to fluid flow problems involving body surface motion in the direction (perpendicular to the line of axisymmetry) was established in the example of Purcell’s toroidal swimmer. Using the regularized ringlet method, we were able to reproduce Leshansky and Kenneth’s Leshansky and Kenneth (2008) results for the scaled propulsion velocity of the toroidal swimmer, propelled by surface tank treading against the direction of motion of its outer surface. Our results show a significant improvement over the asymptotic solution found by integrating a centre-line distribution of rotlets in the limit as slenderness decreases (). The use of regularized ringlets also yields the drag force at all points on the torus surface for the toroidal swimmer, information that is not readily available using the series or rotlet solutions.
In order to elucidate the further benefits of regularization, we studied the case of cytosolic flow in the growing pollen tube. Here the inducing drag force is the result of the directed actomyosin transport of vesicles inside the fluid, making application of singular Stokeslet methods particularly challenging. Since the area over which the force is applied is also relatively large (and somewhat indeterminate), being able to use the regularization parameter to control the spreading of the force is also vital to the solution.
Using the regularized ringlet solution for the steady Stokes equations, we were able to show that it is highly likely that the central actin bundle plays a role in the removal of vesicles from the apical pool. This has long been hypothesised, with vesicle movement through the central region seen to reach speeds of up to Bove et al. (2008) (similar to those of actomyosin transport), but direct evidence has been lacking. Fluid velocity profiles based on drag induced by actomyosin transport of vesicles along the peripheral bundle alone shows that cytosolic flow in the central region does not reach these speeds, strongly suggesting that actomyosin transport takes place on the central actin bundle.
During further investigation, we were able to show that the shape of the cytosolic velocity profile relative to the moving tip is largely independent of the growth speed of the tube (under the assumption that the actin profile is able to keep up with the advancing tip). This provides some insight into how the tube is able to consistently orient its internal architecture across multiple pollen tube species and in spite of changes in growth speed during different phases of oscillatory growth. In a future paper, we will continue our work on the pollen tube by using the cytosolic velocity profiles we have produced in tandem with an advection-diffusion-reaction equation for the spatio-temporal distribution of vesicles in the pollen tube. Investigating the conditions necessary to produce the ‘inverted vesicle cone’ almost universally observed in the apex of the angiosperm pollen tube yields new results concerning the appropriate values for different parameters pertaining to tube growth. We believe that the regularized Stokeslet ring method will be of benefit to further studies in the fluid dynamics of axisymmetric growth and beyond.
Acknowledgements.
The authors acknowledge funding from the Engineering and Physical Sciences Research Council (EPSRC) in the form of a Doctoral Training Award (EP/M508202/1 to JT) and grants (EP/N021096/1 to DJS; EP/M00015X/1 to RJD). The authors would also like to thank Prof Anja Geitmann and Dr Youssef Chebli (both of the Department of Plant Science, Faculty of Agriculture and Environmental Sciences, McGill University, Sainte-Anne-de-Bellevue, Canada) for their valuable insights on pollen tube physiology and growth.
Appendix A Regularized Stokeslets in cylindrical coordinates for ringlet evaluation
The expressions for fluid point and ring point in cylindrical coordinates are
[TABLE]
in which
Appendix B Evaluating the limiting behaviour of
To understand the behaviour of as , it is easiest to consider in the form
[TABLE]
in which
[TABLE]
with and . Evaluating for each and substituting these back into Equations (80) – (84) yields the desired results for and .
The first step in evaluating these limits is to observe that from the definition of , it follows that as either or (or both). Further, noting that
[TABLE]
in which the value of as depends on whether or (or both), it is observed that contain common finite term outside of the larger brackets. Evaluating the remaining parts of in the limit as requires employing the power series expansions of the complete elliptic integrals Gradshteyn and Ryzhik (2014), such that
[TABLE]
Letting it follows that
[TABLE]
Similarly, letting and it can be found that
[TABLE]
Finally, using and yields
[TABLE]
Compiling all of the above gives
[TABLE]
which upon substitution into the expressions for yields
[TABLE]
Appendix C Fluid flow induced by unit ringlet forces
Figure 14 shows the fluid velocity induced by a ringlet located at with associated unit force in and directions, both under free–space conditions ((a) & (c)) and in the presence of a cylindrical wall at ((b) & (d), with the wall also represented by Stokes ringlets). In the presence of the bounding wall, closed streamlines and toroidal eddies are observed.
Appendix D Evaluating the double layer potential
Recall the form of the double layer potential
[TABLE]
in which the stress tensor is given by
[TABLE]
As with the Stokeslet in the SLP, the stress tensor can be expressed in cylindrical coordinates via
[TABLE]
and the Cartesian and cylindrical forms of the flow vector are related via in which
[TABLE]
Recalling that the azimuthal component of the normal to an axisymmetric body is zero, the transformation from polar to Cartesian coordinates is where
[TABLE]
Let
[TABLE]
such that
[TABLE]
For fixed , nonzero elements of are given by
[TABLE]
such that each is a linear sum of terms of the form
[TABLE]
Letting
[TABLE]
in which for odd, each of the necessary can be written as:
[TABLE]
Using the double angle formulae for and we are able to express purely in terms of even powers of , such that
[TABLE]
where . Note that the upper limit of integration was first reduced by application of the double angle formulae, followed by the even parity of the resulting integrand about . If we further define
[TABLE]
then expanding Equation (113) for the relevant values of yields
[TABLE]
in which the integrals can be expressed in terms of complete elliptic integrals of the first and second kind ( and respecitvely) with elliptic modulus as
[TABLE]
Appendix E Translating unit sphere
Following the example of Cortez et al. Cortez et al. (2005), consider the translating unit sphere with velocity The sphere is parametrized and discretized in the same manner as in Section III.1. Assuming the sphere experiences zero azimuthal spin, the fluid velocity at any point can be approximated using the ringlets via Equation (27). Using intervals of equal size, the quadrature weight associated with numerical integration over is simply . Now considering the velocity evaluated at the location of each ring in the plane yields a system of equations which can be written in matrix form
[TABLE]
as outlined in Equations (29) - (31). By setting
[TABLE]
where are the fluid viscosity, sphere radius, and sphere speed (in ) respectively, the classical solution for Stokes flow
[TABLE]
should follow. Note that this calculation is independent of the value of by cancellation (besides the implicit requirement that for the Stokes equations to be valid), and that we use , as described in the outline of the problem.
The first test involves using regularized rings so that the grid size is given by . This is chosen such that the minimum distance between adjacent rings in the axisymmetric discretization of the sphere surface is approximately the same as the distance between adjacent points in Cortez and colleagues’ discretization using regularized Stokeslets. The regularization parameter is varied between and , and the error in the norm for the component of the flow field is recorded in each case. This error is defined as
[TABLE]
Division by is necessary for the sake of comparison of errors with later tests where the value of will change in order to alter the grid size. Initial results using regularized ringlets are shown in Figure 15 and are favourable compared to those using regularized Stokeslets. The regularized ringlet method appears to be optimal for a lower value of than the regularized Stokeslet method, with a minimal error found at with our method and with that of the regularized Stokeslet, as well as being slightly more accurate for almost all values of tested. Interestingly, the magnitude of the errors using the two methods briefly appear to coincide near the point at which the regularized Stokeslet error is minimized.
The second test involves varying the grid size for a fixed value of Regularized ringlet results are shown in Figure 16 and again compare very favourably to results using regularized Stokeslets. For larger grid sizes our errors are significantly reduced compared to those found using regularized Stokeslets. This is at least in part a result of the ringlet method being better suited to handling the small value of . With both methods, the error eventually stops decreasing as the grid size tends towards [math] since in this regime the regularization error dominates.
The final test again looks at the effect of varying the value of on the magnitude of the numerical error, this time using ringlets for a grid size approximately equal to (note that we do not use ringlets to avoid placing a ring at the point , which would result in a singular velocity when ). The velocity error is compared at two distinct points: lying on the surface of the sphere, and lying a distance of half the sphere radius away. Ringlet results are shown in Figures 17a & 17b, and once more match regularized Stokeslet results very closely. The magnitude of the error is linear with respect to on the surface of the sphere and quadratic a sufficient distance away.
Appendix F Comparison to singular solutions
Regularized solutions for Stokes flow have the advantage of being simple to implement and readily usable without needing to worry about the presence of singularities in the computational domain. This does not mean that singular solutions cannot be used, however; so long as the appropriate care is taken to deal with the singularities in some way, singular solutions can also yield excellent results. The method of fundamental solutions Young et al. (2006); Alves and Silvestre (2004) is a popular choice for implementing singular Stokeslet solutions, in which a fictitious boundary is placed outside of the computational domain and adjacent to the physical boundary of the problem considered. Stokeslets (or ‘source points’) are placed on this fictitious boundary and are associated with collocation points (typically of an equal number) on the physical boundary, with the force density for each Stokeslet being calculated using the resistance matrix such that the physical boundary conditions are satisfied. What the appropriate distance between the fictitious and physical boundaries should be is difficult to determine a priori, and in some sense this distance can be regarded as a regularization parameter for the singular problem Barrero-Gil (2012). If the separation distance is too small, the proximity between the Stokeslet singularities and the physical boundary may lead to inaccurate solutions, whereas if the distance is too large the resistance matrix may become ill-conditioned Chen et al. (2005). In some cases, placement of the fictitious boundary may also be constrained by the geometry of the problem itself, leading to solutions that are far from optimal. The difficulty in balancing all of these factors is one of the reasons for the popularity of regularized methods.
In the case of the axisymmetric ring of singular Stokeslets, some interesting behaviour occurs in the limit as the source and collocation points coincide. In this limit, tends to unity from which it follows that and in Equations (53) – (57). By employing the asymptotic expansion in which , we find that which all tend to infinity as but at a significantly slower rate than the individual Stokeslet (). The matrix elements and are similarly divergent but typically take values in the range . Figure 19 shows a comparision of the magnitude of the singularities as for (red circles) versus for (blue squares). A log-log plot must be employed due to the speed with which the Stokeslet singularity increases for small .
As a result of the slow rate with which the singularity in the axisymmetric Stokeslet ring approaches infinity, the method of fundamental solutions can be employed with a separation distance between the fictious and physical boundaries that is several orders of magnitude smaller than is possible using individual Stokeslets. This allows for the discovery of accurate solutions with a well-conditioned underlying resistance matrix. We illustrate this by once again conducting Matlab simulations for the resistance problem on the translating unit sphere (velocity ) using both regularized and singular Stokeslet rings.
As before, the sphere is parametrized in D using the arc segment for with points placed along its length at locations . For the method regularized ringlets, these denote ringlet locations. The velocity is prescribed at all ringlet locations and the resistance matrix is constructed to yield the required force densities . The condition number of the resistance matrix is determined using the \mlttfamilycond function in Matlab. The analytical solution for the drag on the unit sphere with unit velocity is known to be equal to , so that the mean relative error for the calculation of the force term resulting from the resistance matrix is given by
[TABLE]
in which we recall that is the force per unit area exerted by the sphere on the fluid and so has opposite sign to the drag. The fluid velocity at the near– and far–field locations is evaluated using the force densities according to Equation (27). The relative error in these velocities is given by
[TABLE]
based on the analytical solution given by Cortez et al. Cortez et al. (2005). Table 5 shows a summary of the results obtained for each of these four metrics for various values of regularisation parameter .
For the method of fundamental solutions, let denote collocation points on the sphere surface. source points are placed at positions where is a parameter that represents the separation distance between the physical and fictitious boundaries (Figure 18). The velocity is prescribed at all collocation points and the resistance matrix is constructed to yield the required force densities at each source point. The accuracy and applicability of the method is measured using the same four metrics as for the method of regularized ringlets. Table 6 shows a summary of results for various values of separation distance , analogous to the results in the regularized case from Table 5.
Similar results are achievable using both methods, although the singular method generally yields more well-conditioned resistance matrices and can produce relative errors of a smaller magnitude for the drag calculation and fluid velocities. It should however be noted that it is not possible to minimize each error in the singular calculation concurrently (as the smallest values for and are generally associated with larger values for the condition number and ). Despite axisymmetry enabling a drastic reduction of the separation distance for the singular problem, the question of what distance is considered ‘optimal’ still persists. This is not an issue for the regularized ringlet, in which case each error achieves a minimal value for similar values of ( to ). For excessively small (), the error becomes non–monotonic as is often the case with regularized Stokeslet methods (eg. see Figure 15, Cortez et al. Cortez et al. (2005), Gallagher et al. Gallagher et al. (2018)). The same is also true for small in the singular case. For both regularized and singular methods, reduction of regularization parameter or always appears to result in a reduction of the condition number of the resistance matrix.
Although the singular method can be tuned to give smaller relative errors in either the fluid velocity or total drag separately, it cannot do so simultaneously; regularized ringlets display more satisfactory convergence properties and are the more effective method to minimize errors in both fluid velocity and total drag.
Appendix G Computational speed
Table 7 shows a comparison of the computational times and (measured in seconds), associated with constructing the ringlet matrix (in the irrotational fluid case) and the Stokeslet matrix respectively. The increase in computational time for computing is a result of needing to compute the complete elliptic integrals and for all combinations of ring locations , which requires the construction of two further matrices and . The computational time needed for this isolated operation, is also listed in Table 7. The total additional time needed to construct is modest, typically between and of the time needed for .
The computational time associated with and (and by extension, ) can be reduced by evaluating and to a lower degree of accuracy; the Matlab function \mlttfamilyellipke(k,TOL) calculates and to the accuracy defined by TOL, which has a default value of (double-precision accuracy). This is a far greater accuracy than we are typically able to achieve using regularized Stokeslet methods, suggesting a larger value of TOL will suffice.
We note that by consideration of the size of the matrices involved and the typically small value of , we should hypothetically be able to achieve for any given value of . In practice, this is not the case. The elements for share a common form that enables them to be encoded in matrix form very efficiently. The same is not true of the elements for hence why in Table 7.
Inclusion of fluid rotation (such that an additional matrix must be constructed for incurs an additional cost of small value, approximately equal to for any given value of . We divide by under the assumption that the cost associated with constructing the additional matrix is one quarter of the cost associated with constructing the irrotational matrix .
The true value of working with ringlets in D can be seen by by evaluating the cost of inverting the matrices and versus (as is necessary in the resistance problem for evaluating the force associated with a given boundary velocity). These matrices are of size , and respectively. We investigate the cost of solving a linear system in which and are a matrix and a vector of normally distributed random data respectively, with We denote the time taken to solve this system by and in each case. The results are given in Table 8, in which it is clear to see that is significantly larger than the sum of and for all values of tested. In practice, this means that using a D Stokeslet implementation with for solving an axisymmetric resistance problem will always be significantly more costly than our D ringlet implementation with
As well as producing smaller resistance matrices with a reduced associated computational cost for a given number of nodes, our axisymmetric ringlet method also requires far fewer nodes in order to achieve the same level of accuracy as the traditional regularized Stokeslet method. In Appendix E it was shown that for the case of the translating unit sphere the axisymmetric discretization of the sphere surface with rings produces results that are consistently more accurate than the traditional D patch discretization using nodes. To put this into context, the result using ringlets in D (computational time 0.1\text{,}\mathrm{s}) thus corresponds to using $1,500,000$ regularized Stokeslets in $3$D. A lower bound estimate for the associated computational time (using $t_{3N}\propto 3N$ with data from Table [8](#A7.T8)) for method of regularized Stokeslets using this number of nodes yields a time in excess of $500$ days. Assuming each element of the resulting matrix $\mathsf{S}^{\varepsilon}$ is stored as a double-precision floating-point number, this would also require $\approx$1\text{\,}\mathrm{P}\mathrm{B} of storage.
Appendix H Force calculation for the rotating sphere
The leading order expression (in terms of Re) for the torque on the rotating sphere is given by,
[TABLE]
the derivation of which can be found in Rubinow and Keller (1961). This torque is associated with a drag force per unit area on the surface of the sphere given by , which can be verified by considering the identity,
[TABLE]
where denotes the sphere surface. Multiplication by yields,
[TABLE]
In cylindrical coordinates, we have
[TABLE]
Converting to a spherical system in which denotes the azimuthal angle and the polar angle, we substitute and such that,
[TABLE]
as required.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Lauga and Powers (2009) E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, Reports on Progress in Physics 72 , 096601 (2009).
- 2Phan-Thien et al. (1987) N. Phan-Thien, T. Tran-Cong, and M. Ramia, A boundary-element analysis of flagellar propulsion, Journal of Fluid Mechanics 184 , 533 (1987).
- 3Ramia et al. (1993) M. Ramia, D. Tullock, and N. Phan-Thien, The role of hydrodynamic interaction in the locomotion of microorganisms, Biophysical journal 65 , 755 (1993).
- 4Shum et al. (2010) H. Shum, E. A. Gaffney, and D. J. Smith, Modelling bacterial behaviour close to a no-slip plane boundary: the influence of bacterial geometry, Proceedings of the Royal Society of London A 466 , 1725 (2010).
- 5Chwang and Wu (1975) A. T. Chwang and T. Y.-T. Wu, Hydromechanics of low-Reynolds-number flow. Part 2. Singularity method for Stokes flows, Journal of Fluid Mechanics 67 , 787 (1975).
- 6Chwang and Wu (1974) A. T. Chwang and T. Y.-T. Wu, Hydromechanics of low-Reynolds-number flow. Part 1. Rotation of axisymmetric prolate bodies, Journal of Fluid Mechanics 63 , 607 (1974).
- 7Chwang (1975) A. T. Chwang, Hydromechanics of low-Reynolds-number flow. Part 3. Motion of a spheroidal particle in quadratic flows, Journal of Fluid Mechanics 72 , 17 (1975).
- 8Chwang and Wu (1976) A. T. Chwang and T. Y. Wu, Hydromechanics of low-Reynolds-number flow. Part 4. Translation of spheroids, Journal of Fluid Mechanics 75 , 677 (1976).