%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
Design optimisation is observed throughout nature, and it has long been used in the private sector for better and more efficient design. Many industries, particularly aeronautics (Jameson, Martinelli & Pierce Reference Jameson, Martinelli and Pierce1998; Alonso, LeGresley & Pereyra Reference Alonso, LeGresley and Pereyra2009) and ship design (Campana et al. Reference Campana, Liuzzi, Lucidi, Peri, Piccialli and Pinto2009), use fluid dynamical considerations when optimising shapes (Mohammadi & Pironneau Reference Mohammadi and Pironneau2004), such as minimising drag or maximising lift. These optimisations are usually subject to constraints; for instance, a ship might need to have a fixed volume, or certain shapes might be too costly to engineer. In nature, too, optimisation is subject to physical constraints; for example, the beat patterns of spermatozoa (Pironneau & Katz Reference Pironneau and Katz1974; Lauga & Eloy Reference Lauga and Eloy2013) and cilia (Osterman & Vilfan Reference Osterman and Vilfan2011; Eloy & Lauga Reference Eloy and Lauga2012) are subject to energetic constraints on the molecular motors responsible for the motion. Similarly, the shapes of squirming microorganisms correspond closely to those minimising power dissipation subject to a maximum curvature (Vilfan Reference Vilfan2012). Inspired by nature, shape optimisation can be used to improve the design of biomimetic artificial swimmers (Keaveny, Walker & Shelley Reference Keaveny, Walker and Shelley2013).
In a classical paper, Pironneau (Reference Pironneau1973) used a variational approach to show that, in the absence of inertia, the lowest drag shape of fixed volume is front–back symmetric with constant vorticity on its surface. Furthermore, the front and rear of the shape taper into cones of semi-cone angle
$60^{\circ }$
. This shape was then computed numerically by Bourot (Reference Bourot1974) by minimising the drag among a finite family of shapes. Pironneau’s optimal profile has subsequently been observed in the explosively launched spores of ascomycete fungi (Roper et al.
Reference Roper, Pepper, Brenner and Pringle2008). The constraint of fixed volume is appropriate in this biological case, as a certain amount of material is required in order to generate a new organism.
Recently, mechanical optimisations have begun to be examined for more efficient drug delivery, using material rather than chemical properties to achieve specific effects (Champion, Katare & Mitragotri Reference Champion, Katare and Mitragotri2007; Mitragotri & Lahann Reference Mitragotri and Lahann2009; Petros & DeSimone Reference Petros and DeSimone2010). Variations in the shape, size, mechanical and surface properties of drugs can significantly impact drug-delivery processes such as phagocytosis (Champion & Mitragotri Reference Champion and Mitragotri2006), whereby macrophages engulf solid particles. Gratton et al. (Reference Gratton, Ropp, Pohlhaus, Luft, Madden, Napier and DeSimone2008) showed that microscopic slender rod-like particles with aspect ratio three, corresponding to eccentricity
$e=0.9428$
, were subsumed four times faster than spherical particles of equal volume.
A potentially important constraint for designer drug delivery, and also biomimetic microrobots, is fixed particle surface area. Indeed, surface area affects diffusion rates, is associated with toxicity (Monteiller et al.
Reference Monteiller, Tran, MacNee, Faux, Jones, Miller and Donaldson2007) and is important in setting the rate of production of reactants for catalysis and chemical reaction. In this study, we determine the analogous fixed-surface-area shape to Pironneau’s fixed-volume optimal; that is, the convex axisymmetric shape that minimises drag subject to fixed surface area. These constraints restrict us to a family of simply manufactured shapes resembling Pironneau’s optimum, and specifically rule out fractal-like shapes of finite surface area folded into a very small space. We derive the fixed-surface-area shape numerically using a representation whereby its tangent angle is described by a very small number of odd Legendre modes, together with a local analysis of the equations of fluid motion. We show that this ‘other’ optimal shape is almost twice as slender as the fixed-volume shape, has a surface vorticity square proportional to the mean surface curvature (Pironneau Reference Pironneau1973; Bourot Reference Bourot1974) and near either end is locally conical about an angle of
$30.8^{\circ }$
.
2. Mathematical analysis
Two aspects of the optimisation problem are amenable to an analytical approach, namely the optimality condition and a local shape analysis near either end of the body.
2.1. Optimality condition
We begin by deriving the shape’s optimality condition, following Pironneau (Reference Pironneau1973). In Newtonian flows at microscopic scales, the fluid flow velocity
$\boldsymbol{u}$
driven by a body force per unit volume
$\boldsymbol{F}$
is governed by the Stokes flow equations,
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}p={\it\mu}{\rm\nabla}^{2}\boldsymbol{u}+\boldsymbol{F},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$
with pressure
$p$
and viscosity
${\it\mu}$
. Consider a shape
$S$
described by
$\boldsymbol{x}$
translating in this fluid at a fixed speed
$\boldsymbol{U}$
in the absence of fluid body forces (
$\boldsymbol{F}=\mathbf{0}$
). In the reference frame of the body the velocity satisfies the no-slip boundary condition as
$\boldsymbol{u}(\boldsymbol{x})=\mathbf{0}$
and
$\boldsymbol{u}=-\boldsymbol{U}$
at infinity. We seek the shape that minimises the drag. The rate of working of the drag is equal to the rate of viscous dissipation in the fluid, and so we wish to minimise
$$\begin{eqnarray}D=2{\it\mu}\int _{V}e_{ij}e_{ij}\text{d}V,\end{eqnarray}$$
where
$V$
refers to the volume outside of the shape and
$e_{ij}=(\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i})/2$
.
Perturbing
$S$
into a new shape described by
$\boldsymbol{x}+{\it\alpha}\boldsymbol{n}$
, where
$\boldsymbol{n}$
is the normal to the surface and
${\it\alpha}$
is assumed to be a small length everywhere, we denote the resulting change in the velocity field around the shape by
${\it\delta}\boldsymbol{u}$
. The change in the dissipation is given by
$$\begin{eqnarray}{\it\delta}D={\it\mu}\int _{V}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\left(\frac{\partial {\it\delta}u_{i}}{\partial x_{j}}+\frac{\partial {\it\delta}u_{j}}{\partial x_{i}}\right)\text{d}V.\end{eqnarray}$$
Using integration by parts and symmetries in the indices we obtain
$$\begin{eqnarray}{\it\delta}D=-2{\it\mu}\int _{V}{\it\delta}u_{i}\left(\frac{\partial ^{2}u_{i}}{\partial x_{j}\partial x_{j}}+\frac{\partial ^{2}u_{j}}{\partial x_{i}\partial x_{j}}\right)\text{d}V+2{\it\mu}\int _{S}{\it\delta}u_{j}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)n_{i}\text{d}S.\end{eqnarray}$$
Since
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$
, we have identically
$\partial _{j}\partial _{i}u_{j}=0$
. The first term in the first integral in (2.4) is the Laplacian of the velocity field, which we know from (2.1a,b
) in the absence of body forces is given by
${\it\mu}\partial _{j}\partial _{j}u_{i}=\partial _{i}p$
. Furthermore, a Taylor expansion of the no-slip boundary condition
$(\boldsymbol{u}+{\it\delta}\boldsymbol{u})(\boldsymbol{x}+{\it\alpha}\boldsymbol{n})=\mathbf{0}$
allows us to obtain
$$\begin{eqnarray}{\it\delta}\boldsymbol{u}(\boldsymbol{x})=-{\it\alpha}\frac{\partial \boldsymbol{u}}{\partial n}\cdot\end{eqnarray}$$
With these two results, (2.4) simplifies to
$$\begin{eqnarray}{\it\delta}D=-2\int _{V}{\it\delta}u_{i}\frac{\partial p}{\partial x_{i}}\text{d}V-2{\it\mu}\int _{S}{\it\alpha}\frac{\partial u_{j}}{\partial n}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)n_{i}\text{d}S.\end{eqnarray}$$
Using integration by parts, the first integral in (2.6) becomes
$$\begin{eqnarray}\int _{V}-p\frac{\partial {\it\delta}u_{i}}{\partial x_{i}}\text{d}V+\int _{S}p{\it\delta}u_{i}n_{i}\text{d}S.\end{eqnarray}$$
The first term in (2.7) is zero because the perturbation flow is divergence free, while the second term is zero because of (2.5) and the divergence-free condition. Because of the no-slip boundary condition for
$\boldsymbol{u}$
and the fact that it is divergence free, the last integral in (2.6) can be further simplified to
$$\begin{eqnarray}{\it\delta}D=-2{\it\mu}\int _{S}{\it\alpha}\frac{\partial u_{j}}{\partial n}\frac{\partial u_{j}}{\partial n}\text{d}S=-2{\it\mu}\int _{S}{\it\alpha}\bigg|\frac{\partial \boldsymbol{u}}{\partial n}\bigg|^{2}\text{d}S.\end{eqnarray}$$
In the fixed-surface-area case, we consider a cost function of the form
$$\begin{eqnarray}J=D+{\it\lambda}\left(\int _{S}\text{d}S-S_{0}\right)\!,\end{eqnarray}$$
where
${\it\lambda}$
is the Lagrange multiplier enforcing the constraint of fixed surface area,
$S=S_{0}$
. Variations in
$S$
arising from the perturbation to the shape are given by
$2\int {\it\kappa}{\it\alpha}\text{d}S$
, where
${\it\kappa}$
is the mean surface curvature, and thus the perturbation to
$J$
is given by
$$\begin{eqnarray}{\it\delta}J={\it\delta}D+{\it\lambda}{\it\delta}S=-2{\it\mu}\int _{S}{\it\alpha}\bigg|\frac{\partial \boldsymbol{u}}{\partial n}\bigg|^{2}\text{d}S+2{\it\lambda}\int _{S}{\it\kappa}{\it\alpha}\text{d}S=2\int _{S}{\it\alpha}\left(-{\it\mu}\bigg|\frac{\partial \boldsymbol{u}}{\partial n}\bigg|^{2}+{\it\lambda}{\it\kappa}\right)\text{d}S.\end{eqnarray}$$
By setting
${\it\delta}J=0$
, the optimality condition is then
$$\begin{eqnarray}|\partial \boldsymbol{u}/\partial n|^{2}\propto {\it\kappa},\end{eqnarray}$$
where the first term can also be seen as the norm squared of the surface vorticity. It should be noted that in the fixed-volume case, the cost function is
$J=D+{\it\lambda}(V-V_{0})$
, and since
${\it\delta}V=\int {\it\alpha}\text{d}S$
, the optimal shape has constant surface vorticity (Pironneau Reference Pironneau1973).
(a) Flow over the vertex of a cone illustrating the streamfunction for the minimal-drag fixed-volume case of
${\it\theta}_{0}=120^{\circ }$
. (b) The function
$H({\it\nu})$
in (2.17) for a selection of cone angles, highlighting the minimal-drag fixed-surface-area case of
${\it\theta}_{0}=149.179^{\circ }$
.
2.2. Local analysis
As shown by Pironneau (Reference Pironneau1973) and Bourot (Reference Bourot1974), optimal shapes have pointy American-football-like ends. The second aspect of the optimisation problem that can be treated analytically is a derivation of the optimal angle at either end of the locally conical shape.
We consider a conical shape and describe the axisymmetric flow around it using spherical coordinates
$(r,{\it\theta})$
, with
$r$
measured from the tip and
${\it\theta}$
from the axis of symmetry in the fluid (figure 1
a). The cone is assumed to be located at
${\it\theta}={\it\theta}_{0}$
, where
${\it\theta}_{0}$
may be
$90^{\circ }$
in the case of a blunt body, and the goal of the analysis is to derive the value of
${\it\theta}_{0}$
consistent with the optimality condition. Locally around the cone, the streamfunction
${\it\psi}(r,{\it\theta})$
is expected to take the asymptotic form
${\it\psi}=r^{n}f({\it\theta})$
. The equations and boundary conditions satisfied by the streamfunction are
$$\begin{eqnarray}\displaystyle E^{2}(E^{2}{\it\psi})=0,\quad E^{2}=\frac{\partial ^{2}}{\partial r^{2}}+\frac{\sin {\it\theta}}{r^{2}}\frac{\partial }{\partial {\it\theta}}\left[\frac{1}{\sin {\it\theta}}\frac{\partial }{\partial {\it\theta}}\right],\quad \left.\frac{\partial {\it\psi}}{\partial {\it\theta}}\right|_{{\it\theta}_{0}}=0,\quad {\it\psi}({\it\theta}_{0})={\it\psi}(0)=0. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
In the case of fixed volume, the optimality condition of constant surface vorticity yields
$$\begin{eqnarray}|\partial \boldsymbol{u}/\partial n|=|E^{2}{\it\psi}/r\sin {\it\theta}|=\text{const.}\Rightarrow n=3.\end{eqnarray}$$
The solution for
$f({\it\theta})$
given in Pironneau (Reference Pironneau1973) does not, however, satisfy the condition that
${\it\psi}({\it\theta}_{0})=0$
(2.12a−d
), and we rederive it here. Plugging the ansatz
${\it\psi}=r^{3}f({\it\theta})$
into (2.12a−d
) leads to the ordinary differential equation
$$\begin{eqnarray}f^{\prime \prime }-\cot {\it\theta}f^{\prime }+6f=A\cos {\it\theta}+B.\end{eqnarray}$$
Together with the particular solution, the full solution is then
${\it\psi}=cr^{3}(at+bt^{3}+1)$
for
$t=\cos {\it\theta}$
. In order to enforce the boundary conditions, the constants must be
$a=-3/2t_{0}$
,
$b=1/2t_{0}^{3}$
. Finally, since the streamline along
${\it\theta}=0$
(or
$t=1$
) is split by the cone, we require
$f(0)=0$
, i.e.
$$\begin{eqnarray}\frac{1}{2t_{0}^{3}}-\frac{3}{2t_{0}}+1=0\Rightarrow (t_{0}+1/2)(t_{0}-1)^{2}=0,\end{eqnarray}$$
and therefore
$t_{0}=-1/2$
or
${\it\theta}_{0}=120^{\circ }$
. The optimal fixed-volume body thus has locally conical ends with a semi-angle
${\it\theta}_{sc}=60^{\circ }$
, agreeing with the result of Pironneau (Reference Pironneau1973). The streamlines associated with this flow are displayed in figure 1(a).
For the fixed-surface-area shape, the optimality condition depends on the mean curvature
${\it\kappa}$
of the surface (2.11). It should be noted that this means that we expect the vorticity on the surface of our optimal shape to diverge towards its front and rear ends. At the conical tip,
${\it\kappa}\propto 1/r$
, and following the same local analysis, this condition may be written as
$$\begin{eqnarray}r|\partial \boldsymbol{u}/\partial n|^{2}=r|E^{2}{\it\psi}/r\sin {\it\theta}|^{2}=\text{const.}\Rightarrow n=5/2.\end{eqnarray}$$
In order to derive the cone angle corresponding this value
$n=5/2$
, we turn to Wakiya (Reference Wakiya1976), who considered Stokes flows near the tips of rigid bodies. In particular, Wakiya showed that, for a given cone angle
${\it\theta}_{0}$
, the exponent
$n$
of the streamfunction satisfying (2.12a−d
) is given by the first non-trivial root (i.e.
${\it\nu}\neq 0,1$
) of
$$\begin{eqnarray}H({\it\nu})=P_{{\it\nu}}(t_{0})\left[t_{0}P_{{\it\nu}}^{\prime \prime }(t_{0})+P_{{\it\nu}}^{\prime }(t_{0})\right]-t_{0}\left[P_{{\it\nu}}^{\prime }(t_{0})\right]^{2},\end{eqnarray}$$
where
${\it\nu}+2=n$
and
$P_{{\it\nu}}$
is the Legendre function of the first kind of, possibly fractional, degree
${\it\nu}$
. The function
$H({\it\nu})$
is plotted for cone angles in figure 1(b). The fixed-volume shape is the special case where
${\it\nu}=1$
is a double root of (2.17). In the case of a fixed-surface constraint, we have
${\it\nu}=1/2$
, and we can numerically invert (2.17) to obtain the corresponding value of
$t_{0}$
. This is illustrated in figure 1(b), and we obtain
${\it\theta}_{0}=149.179^{\circ }$
in our coordinate system, corresponding to a semi-cone angle of
$$\begin{eqnarray}{\it\theta}_{sc}\approx 30.8^{\circ }.\end{eqnarray}$$
We thus expect to obtain a more slender optimal shape than in the fixed-volume case.
3. Numerical shape optimisation
We now proceed with a computational approach in order to derive the optimal shape and confirm the validity of our local analysis. The fundamental solution to the Stokes flow (2.1a,b
) driven by a point force
$\boldsymbol{f}{\it\delta}(\boldsymbol{x}-\boldsymbol{y})$
located at
$\boldsymbol{y}$
is given by the Stokeslet tensor
$\unicode[STIX]{x1D61A}_{ij}(\boldsymbol{x},\boldsymbol{y})$
,
$$\begin{eqnarray}u_{i}(\boldsymbol{x})=\unicode[STIX]{x1D61A}_{ij}(\boldsymbol{x},\boldsymbol{y})f_{j},\quad \unicode[STIX]{x1D61A}_{ij}(\boldsymbol{x},\boldsymbol{y})=\frac{1}{8{\rm\pi}}\left(\frac{{\it\delta}_{ij}}{r}+\frac{r_{i}r_{j}}{r^{3}}\right)\!,\end{eqnarray}$$
where
$r_{i}=x_{i}-y_{i}$
and
$r^{2}=r_{1}^{2}+r_{2}^{2}+r_{3}^{2}$
. For axisymmetric flow, the Green’s function corresponding to a ring of point forces
$M_{{\it\alpha}{\it\beta}}(\boldsymbol{x}-\boldsymbol{x}_{0})$
can then be obtained by integrating the Stokeslet in the azimuthal direction (Pozrikidis Reference Pozrikidis1992),
$$\begin{eqnarray}\displaystyle & \displaystyle M_{xx}=2k\left(\frac{{\it\sigma}_{0}}{{\it\sigma}}\right)^{1/2}\left(F+\frac{\hat{x}^{2}}{\hat{r}^{2}}E\right)\!,\quad M_{x{\it\sigma}}=-k\frac{\hat{x}}{{\it\sigma}_{0}{\it\sigma}^{1/2}}\left[F-({\it\sigma}^{2}-{\it\sigma}_{0}^{2}+\hat{x}^{2})\frac{E}{\hat{r}^{2}}\right], & \displaystyle \nonumber\\ \displaystyle & \displaystyle M_{{\it\sigma}x}=k\frac{\hat{x}}{{\it\sigma}}\left(\frac{{\it\sigma}_{0}}{{\it\sigma}}\right)^{1/2}\left[F+({\it\sigma}^{2}-{\it\sigma}_{0}^{2}-\hat{x}^{2})\frac{E}{\hat{r}^{2}}\right], & \displaystyle \nonumber\\ \displaystyle & \displaystyle M_{{\it\sigma}{\it\sigma}}=\frac{k}{{\it\sigma}_{0}{\it\sigma}}\left(\frac{{\it\sigma}_{0}}{{\it\sigma}}\right)^{1/2}\left\{({\it\sigma}_{0}^{2}+{\it\sigma}^{2}+2\hat{x}^{2})F-\left[2\hat{x}^{4}+3\hat{x}^{2}({\it\sigma}_{0}^{2}+{\it\sigma}^{2})+({\it\sigma}^{2}-{\it\sigma}_{0}^{2})^{2}\right]\frac{E}{\hat{r}^{2}}\right\}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Here,
$F$
and
$E$
are complete elliptic integrals of the first kind of argument
$k$
,
$$\begin{eqnarray}F(k)=\int _{0}^{{\rm\pi}/2}\frac{\text{d}{\it\omega}}{(1-k^{2}\cos ^{2}{\it\omega})^{1/2}},\quad E(k)=\int _{0}^{{\rm\pi}/2}(1-k^{2}\cos ^{2}{\it\omega})^{1/2}\text{d}{\it\omega},\end{eqnarray}$$
where
$k^{2}=4{\it\sigma}{\it\sigma}_{0}/[\hat{x}^{2}+({\it\sigma}+{\it\sigma}_{0})^{2}]$
, while
$\hat{x}=x-x_{0}$
and
$r^{2}=\hat{x}^{2}+({\it\sigma}-{\it\sigma}_{0})^{2}$
are axial and radial coordinates
$(x,{\it\sigma})$
.
The Newtonian fluid flow around our axisymmetric body can then be reduced to a boundary integral of ring forces taken over the meridional plane,
$$\begin{eqnarray}u_{{\it\alpha}}(\boldsymbol{y})=\frac{1}{8{\rm\pi}}\int _{C}M_{{\it\alpha}{\it\beta}}(\boldsymbol{x},\boldsymbol{y})f_{{\it\beta}}(\boldsymbol{y})\text{d}l,\end{eqnarray}$$
where the double layer potential has been eliminated (rigid-body motion). The meridional plane of the rigid body is discretised into straight-line segments of constant force per unit length; i.e. the components
$f_{x},f_{{\it\sigma}}$
in (3.4) are constant over each element
$s_{i}$
. The velocity of the body is specified at the midpoint of each element. Numerical evaluation of each non-singular line integral is performed with four-point Gaussian quadrature, while singular integrals have the Stokeslet singularity removed and are integrated analytically. The axisymmetric Green’s function is evaluated using BEMLIB (Pozrikidis Reference Pozrikidis2002). The drag on the body is then calculated by summing the products of the area of each element by the calculated traction.
By computing the drag, we can create a series of profiles with decreasing drag relative to the equivalent sphere. For this, we adopt the approach of Bourot (Reference Bourot1974) whereby the drag is minimised within a family of parametrised shapes, a method also used by Zabarankin (Reference Zabarankin2013) to determine optimal profiles in elastic media. By the symmetry of the underlying equations, we restrict our attention to axisymmetric shapes with front–back symmetry (Pironneau Reference Pironneau1973). An efficient description of such shapes may be obtained by expanding their tangent angle,
${\it\phi}(s)$
, in odd Legendre polynomials,
$P_{2n-1}$
,
$$\begin{eqnarray}{\it\phi}(s)=\mathop{\sum }_{n=1}^{N}A_{n}P_{2n-1}(s),\end{eqnarray}$$
for
$s$
arclength along the meridional plane. The
$(x,{\it\sigma})$
positions of the shape are then
$$\begin{eqnarray}x=\int _{C}\cos [{\it\phi}(s)]\text{d}s,\quad {\it\sigma}=\int _{C}\sin [{\it\phi}(s)]\text{d}s.\end{eqnarray}$$
Integration of (3.6a,b
) is performed numerically by Simpson’s rule over an ultrafine mesh, prior to the discretisation of the boundary into straight-line elements. The volume and surface area of the resultant shape are calculated by summing conical frustra defined by the ultrafine mesh. Optimisation is performed for
$n=1\rightarrow N$
sequentially, and the converged optimal coefficients for the
$n=N-1$
shape are used as an initial guess for the
$n=N$
shape, with
$A_{N}$
initially set to zero. Optimisation is carried out using the fmincon function in Matlab, employing a sequential quadratic programming (SQP) algorithm. In constrained optimisations, the sum of the Legendre coefficients is forced to equal the analytically calculated tip angle at every iterative step. Hessians for each step are computed numerically via finite differences.
In order to validate the numerical procedure, we compare the drag on ellipsoids of varying eccentricity and fixed volume as calculated by our boundary element code with the analytical result of Chwang & Wu (Reference Chwang and Wu1975),
$$\begin{eqnarray}\text{Drag}=6{\rm\pi}{\it\mu}aUC_{f1},\quad C_{f1}=\frac{8}{3}e^{3}\left[-2e+(1+e^{2})\log \frac{1+e}{1-e}\right]^{-1}.\end{eqnarray}$$
The percentage relative errors in the drag between our work and Chwang & Wu (Reference Chwang and Wu1975) for 100, 200, 400 and 800 mesh elements are, respectively, 0.0102 %, 0.00256 %, 0.000641 % and 0.000160 %. The drag of successive optima for increasing numbers of modes converges very quickly. As a very high degree of accuracy is required to fully resolve the optimisation for higher numbers of modes, we discretise the shape into 800 elements. As a further validation to ensure that forces are sufficiently resolved around the conical ends of spindle-like bodies, the drag force on a spindle of
$120^{\circ }$
is found to be in agreement with that provided by Wakiya (Reference Wakiya1976).
As a second numerical check, we verify that expansion of the tangent angle in terms of Legendre modes provides an efficient description of the lowest drag shape of fixed volume. We consider both an optimisation scheme whereby the front and rear semi-cone angles are constrained to
$60^{\circ }$
(see above), and a fully unconstrained optimisation. Figure 2(a) shows our constrained and unconstrained optima for six Legendre modes in comparison with the shape given by Bourot (Reference Bourot1974), described by 10 coefficients of the equation
$$\begin{eqnarray}\frac{r({\it\theta})}{{\it\lambda}}=2-\frac{2}{\sqrt{3}}\sin {\it\theta}+\mathop{\sum }_{n=1}^{10}B_{n}\sin ^{n+1}{\it\theta}.\end{eqnarray}$$
The values for the coefficients given by Bourot (Reference Bourot1974) are truncated to 4 s.f., and their shape does not quite have the same volume as a sphere of unit radius. Nonetheless, it is clear that our description of the tangent angle by odd Legendre modes yields the optimal shape quickly, see e.g. the convergence of the relative drag in figure 2(b). This rapid convergence is why a large number of elements (
$N=800$
) are required in the unconstrained optimisation for the angle at the tip to converge to
$60^{\circ }$
. It should be noted that the optimal shape for only three Legendre modes is almost indiscernible from our final optimal. After six modes in the unconstrained optimisation, the calculated tip angle converges with 2 % relative error. The relative drag of the optimal shape (ratio of the drag to that of a sphere of equivalent volume) is
$D_{r}=0.9542549$
, in agreement with the result of Bourot,
$D_{r}=0.95425$
. Pironneau’s optimality condition of a constant surface vorticity is also confirmed to be valid in figure 2(c).
(a) The lowest drag shape of fixed volume as calculated by constrained (fixed front/rear angle) and unconstrained optimisation, together with the shape originally given by Bourot (Reference Bourot1974). (b) The drag relative to a sphere of unit volume of successive optima as found by our constrained optimisation. (c) Pironneau’s optimality condition over the surface of Bourot’s optimum, the optimal prolate spheroid of fixed volume and our constrained optimum.
In order to provide a baseline comparison for the numerical shapes, the drag was first minimised among the family of prolate spheroids of varying eccentricity for both fixed volume and fixed surface area using the analytical result in (3.7a,b
) and our boundary element code. The computational code achieves results identical to within
$10^{-8}$
of the analytical ones and the optimal spheroids have eccentricities of
$e=0.8588$
(aspect ratio 1.952) in the fixed-volume case and 0.9688 (aspect ratio 4.037) in the fixed-surface-area case. We show these results in figure 3, and they allow us to anticipate two key facts of drag minimisation for fixed surface areas. First, we can predict that the lowest drag shape of fixed surface area is likely to be more slender than its fixed-volume counterpart (this was already anticipated in the local analysis above). Furthermore, we can expect that the reduction in drag of the fixed-surface-area optimum relative to its equivalent sphere will be greater than that for the fixed-volume counterpart. Indeed, for spheroids, the reduction in drag is 11 % for fixed surface area but only 4.45 % for fixed volume.
(a) The drag on prolate spheroids of fixed volume (dark blue) and surface area (light yellow) relative to their equivalent spheres as a function of eccentricity. The minimum-drag prolate spheroids are marked with pluses. (b) The shapes of the minimum-drag prolate spheroids of fixed volume and surface area, normalised to unit volume.
We now compute the fixed-surface-area analogue of Pironneau’s shape using both unconstrained and constrained (i.e. fixed-tip-angle) optimisation. After six modes in the unconstrained optimisation, the calculated tip angle converges to the analytically calculated optimal of
$30.8^{\circ }$
(2.18) with a 1 % relative error. The optimal shapes found with both optimisation approaches are displayed in figure 4(a). The unconstrained and constrained results are indistinguishable. We also included in figure 4(a) the fixed-volume optimal, with all shape volumes normalised to one in order to facilitate comparison. As anticipated, the fixed-area optimal is much more slender than the best shape of fixed volume (aspect ratio of 4.162 against 2.109, which are both slightly above the aspect ratios of the optimal spheroids). As with the volume constraint, the drag converges here to its optimal value with very few modes, which is illustrated in figure 4(b). The first six unconstrained coefficients describing the tangent angle in (3.5) are
$$\begin{eqnarray}\left.\begin{array}{@{}rclrclrcl@{}}A_{1}\ & =\ & -0.47956, & A_{2}\ & =\ & -0.035208, & A_{3}\ & =\ & -0.0096958,\\ A_{4}\ & =\ & -0.0042400, & A_{5}\ & =\ & -0.0022198, & A_{6}\ & =\ & -0.0012228.\end{array}\right\}\end{eqnarray}$$
(a) The lowest drag shape of fixed surface as calculated by unconstrained (
$+$
) and constrained optimisation (i.e. fixed tip angles, light yellow line), in comparison with the fixed-volume case (dark blue line). Both shapes are normalised to have unit volume. (b) The drag relative to a sphere of unit surface area for successive optima using constrained optimisation. (c) Vorticity (
$+$
) along the surface of the optimal shape of fixed surface area shown with its mean curvature, scaled by the ratio of vorticity to curvature at
${\it\theta}=90^{\circ }$
(solid purple line).
How much better than the equivalent sphere is the optimal shape? The drag ratio we obtain is equal to
$D_{r}=0.8871581$
, meaning a drag reduction of about 11.3 %. This is significantly greater than the 4.57 % reduction for the case of fixed volume. Finally, we are able to test numerically the optimality criterion derived in (2.11). We plot in figure 4(c) both the value of the flow vorticity on the surface and the mean surface curvature,
${\it\kappa}$
. As can be seen, our computational results show a perfect proportionality between both, confirming the validity of the local analytical criterion discussed in § 2.1.
4. Conclusion
Shape optimisation at small scales is of potential significance for problems in chemistry and colloidal science. In this work, we used the boundary element method together with a mathematical analysis to determine computationally the convex axisymmetric shape that minimises drag in Stokes flow subject to the constraint of a fixed surface area. The shape is about twice as slender as the fixed-volume case, and the front and rear ends of the shape are tangent to a cone of semi-angle
${\approx}30.8^{\circ }$
. This was borne out by minimisation of the drag with respect to coefficients of an efficient representation of the shape’s tangent angle in terms of Legendre modes. Although we did not mathematically prove that the shape we obtained is globally optimal, or in fact unique, among all axisymmetric convex ones, since it quickly converges to it from all initial conditions, we conjecture that it is.
From a practical standpoint, while the optimality conditions for the fixed-volume and fixed-surface-area shapes entail the vorticity on the surfaces of these shapes, the drag turns out to be somewhat insensitive to details of the surface vorticity. In fact, the optimal prolate spheroids of fixed volume and surface area both have zero vorticity at the front and rear ends, but only exhibit 0.1 % and 0.3 % more drag than their respective full optima. This insensitivity engenders the notion of ‘nearly optimal design’ in Stokes flow; optimisation of simple shapes can be effective and helpful to experimental design.
Acknowledgements
This work was funded in part by a European Union Marie Curie CIG to E.L. and a Royal Commission for the Exhibition of 1851 Research Fellowship to T.D.M.-J.
