3.1. Numerical method to obtain the stream function

For a given body $B$ , we solve the biharmonic equation (2.2) for the streamfunction $\psi (x,y)$ with boundary conditions (2.3)–(2.5) using MATLAB’s PDE toolbox. The computational domain is the rectangle $R=[-L,L]\times [-1,1]$ with the body $B$ removed, where $L\gg 1$ is chosen to be sufficiently large such that doubling $L$ results in a change in the drag of less than $5 \times 10^{-5}$ . For bodies with area of order unity, this condition is generally satisfied by using $L=20$ . Results for the flow around some elliptical bodies are shown in figure 2 (not all the computational domain is shown); these bodies are drag-minimising over the family of elliptical bodies. Results for the flow around some optimal bodies are shown in figures 3 and 4.

3.2. Shape optimisation

The minimisation of the drag under the constraint of a given fixed area (2.8) is achieved by parametrising the body shape using (Zabarankin Reference Zabarankin2013)

(3.1) \begin{equation} \zeta (t) = x(t) + i y(t) = C \sum _{j=0}^N \left ( a_j \cos \frac {\pi t}{2} + i b_j \sin \frac {\pi t}{2} \right ) T_j(2t-1), \end{equation}

where $t \in [0,1]$ and $T_j(\boldsymbol{\cdot })$ are the Chebyshev polynomials of the first kind and $C$ is a normalisation constant that fixes the body area. Equation (3.1) gives the body shape in the positive quadrant and the rest of the shape is obtained via reflection.

For any integer $N \geq 0$ , solution to the optimisation problem (2.8) involves finding the values of $a_j$ and $b_j$ that minimise the drag (note that there are $2N+2-1=2N+1$ free parameters owing to the fixed area condition, which manifests via determining $C$ ). For $N=0$ , the shape (3.1) is simply an ellipse. The representation (3.1) allows for any sharp angle at the noses provided that $N \gt 0$ ; we will use this explicitly in § 5.3.

Figure 3. Flow past the general optimal body with areas (a,c) $A=0.1$ and (b,d) $A=1$ : (a,b) streamfunction and streamlines; (c,d) vorticity.

Figure 4. Flow past the general optimal body with area $A=4$ : (a) streamfunction and streamlines; (b) vorticity.

The half-extent of the shape in the $x$ and $y$ directions is given by

(3.2) \begin{equation} a= {\textrm{max}} \left (x\right ) = C \sum _{j=0}^N a_j (-1)^j,\qquad b= {\textrm{max}} \left (y\right ) = C\sum _{j=0}^N b_j. \end{equation}

To obtain the drag-minimising shape, the strategy for a given area $A$ uses alternating descent and proceeds by first finding the optimal elliptical shape ( $N=0$ ), which is a one-parameter optimisation problem. We then increase $N\to N+2$ , and vary each parameter $a_j$ and $b_j$ in turn seeking the smallest drag and using the elliptical body as the initial guess. This iteration continues until a desired tolerance is reached and we then increase $N$ again. We generally continued until $N=6$ , beyond which the reduction in the drag is typically less than $0.001\,\%$ . Results for the optimal shapes of different areas are shown in figures 5 and 6 (black dots show the results of the simple ‘unconstrained’ optimisation described here).

Figure 5. Optimal shapes with areas (a,b) $A=0.025$ , (c,d) $A=0.1$ , (e,f) $A=0.5$ : result for a parameterisation (3.1) with (a,c,e) $N=2$ ; (b,d,f) $N=4$ . The black dots are the shapes from the original optimisation approach (see § 3.2) and the blue lines are from the improved, constrained optimisation (see § 5.3). The red dashed lines in panels (a)–(d) are Richardson’s shape (6.1), accurate for $A\ll 1$ . The shapes are not smooth at the two noses and the apices have internal angles less that $\pi$ .

Figure 6. Optimal shape with area $A=1$ using a parametrisation (3.1) with (a) $N=2$ and (b) $N=4$ . (c) Optimal shape with area $A=2$ and $N=4$ . The black dots are the shapes from the original optimisation approach (see § 3.2) and the blue lines are from the improved, constrained optimisation (see § 5.3). In panel (c), the magenta dash-dotted line is the leading-order asymptotic prediction (6.17), whilst the red dashes show the minimal drag shape calculated from lubrication theory; see also figure 12.

4.1. Small elliptical bodies

In this subsection, we calculate the optimal elliptical shape for $A \ll 1$ ; see also figure 2(a,b) for the corresponding numerical result. The optimal shape (and the associated flow field) are obtained via asymptotic analysis. The strategy is to construct the outer flow (away from the small body), which contains an undetermined coefficient in the streamfunction before determining the inner flow (which depends on the body shape) and then asymptotically matching the two solutions to determine the coefficients.

4.1.1. Outer solution

In the outer problem, the effect of the small body arises as a Stokeslet at the origin, without having to resolve its detailed shape. We write the streamfunction as

(4.1) \begin{equation} \psi _{\textit{out}} = y-\frac {y^3}{3} + \hat {\psi }, \end{equation}

where the first two terms are associated with the far-field flow and the third term accounts for a Stokeslet at the origin superposed with its image terms to satisfy the boundary conditions on the walls (Happel & Brenner Reference Happel and Brenner1983; Crowdy & Davis Reference Crowdy and Davis2013). The strength of the Stokeslet (which will end up being proportional to the leading-order drag on the small body) is shape-dependent and will be determined by matching to an appropriate inner region.

Away from the origin, we retain the biharmonic equation (2.2) ${\nabla} ^4 \hat {\psi }=0$ . We seek a solution for $\hat {\psi }$ that satisfies the boundary conditions as $x \to \pm \infty$ and on the two plates, $y=\pm 1$ (ignoring the matching conditions as $x,y \to 0$ for now);

(4.2) \begin{align} \hat {\psi } & \to 0 \,\,\qquad \qquad \qquad {\textrm{as}} \qquad\qquad |x| \to \infty , \end{align}

(4.3) \begin{align} \hat {\psi } &= \frac {\partial \hat {\psi }}{\partial y}=0 \qquad \qquad {\textrm{on}}\qquad\qquad y =\pm 1. \end{align}

To obtain the solution, we superpose the images of the Stokeslet and take a Fourier cosine transform in $x$ ,

(4.4) \begin{equation} \varPsi ( k, y) = \int _0^{\infty } \hat {\psi }(x,y) \cos kx \, {\textrm{d}} x, \qquad \hat {\psi }(x,y)=\frac {2}{\pi } \int _0^{\infty } \varPsi ( k, y) \cos kx \, {\textrm{d}} k. \end{equation}

The biharmonic equation (2.2) is then recast as

(4.5) \begin{equation} \frac {\partial ^4 \varPsi }{\partial y^4} - 2 k^2 \frac {\partial ^2 \varPsi }{\partial y^2} + k^4 \varPsi = 0. \end{equation}

We apply the boundary conditions at $y=\pm 1$ and in the far-field, and then invert the Fourier transform to obtain (Happel & Brenner Reference Happel and Brenner1983; Pozrikidis Reference Pozrikidis1992; Crowdy & Davis Reference Crowdy and Davis2013)

(4.6) \begin{equation} \hat {\psi } = \sum _{n=0}^\infty \left [B_{2n} f_{2n}(x,y) +C_{2n} g_{2n}(x,y)\right ], \end{equation}

where

(4.7) \begin{align} f_{2n}(x,y) =& \int _0^{\infty } k^{2n} \left [y e^{- k |y|} + \frac {(2 k-1 +e^{-2 k}) y \cosh k y - 2 k \sinh k y}{\sinh 2 k - 2 k} \right ] \frac {\cos k x}{ k} \, {\textrm{d}} k, \nonumber \\ g_{2n}(x,y)=& \int _0^{\infty } k^{2n} \left [\frac {y}{|y|} e^{- k |y|} + \frac {2 k y \cosh k y -(1+ 2 k +e^{-2 k})\sinh k y}{\sinh 2 k - 2 k} \right ] \cos k x \, {\textrm{d}} k, \end{align}

where we have exploited the symmetry properties of the flow in both the $x$ and $y$ directions.

To understand which terms in (4.6) are relevant to our analysis, it is helpful to understand the behaviour of the outer solution (4.6) for $r=\sqrt {x^2 +y^2} \ll 1$ . The component functions (4.7) have the behaviour

(4.8) \begin{equation} f_{2n} \sim r^{1-2n} \quad {\textrm{and}} \quad g_{2n} \sim r^{-1-2n} \quad {\textrm{as}} \quad r \to 0. \end{equation}

We will later show that a small body of length scale $\varepsilon \ll 1$ (and inner radial coordinate $R = r/\varepsilon = O(1)$ ) has inner solution that grows as $O(R \log R)$ in its far-field. Hence, treating logarithms as $O(1)$ and matching between regions yields that the coefficients in (4.6) have asymptotic sizes of $B_{2n} = \mathcal{O}(\varepsilon ^{2n})$ and $C_{2n}=\mathcal{O}(\varepsilon ^{2n+2})$ . Therefore, the leading-order version of (4.6) is

(4.9) \begin{align} \hat {\psi } \sim B_0 f_0(x,y). \end{align}

Using (4.9) in the outer stream function (4.1), we can deduce that the dominant behaviour of (4.1) for small $r$ is

(4.10) \begin{equation} \psi _{\textit{out}} \sim \left [1 +B_0 \left (\log \frac {1}{r} -\alpha _1 \right )\right ] y, \end{equation}

where the constant $\alpha _1$ is given by the following expression:

(4.11) \begin{equation} \alpha _1 = \int _0^{\infty } \frac {1-2k+2k^2-e^{-2k}}{k(\sinh 2k -2k)}-\frac {e^{-k}}{k} \, {\textrm{d}} k\approx 0.4157. \end{equation}

With the outer solution in hand, we next obtain the inner solution for an ellipse and formally match this to the outer solution.

4.1.2. Inner solution for an elliptical body

For the inner problem in the vicinity of the ellipse, we write

(4.12) \begin{equation} x=\varepsilon X, \quad y= \varepsilon Y, \quad r=\varepsilon R, \end{equation}

where

(4.13) \begin{equation} \varepsilon =\sqrt {A/\pi } \ll 1, \end{equation}

and so $\varepsilon$ represents the radius of a circle of equivalent area to the ellipse under consideration. The ellipse boundary is defined by

(4.14) \begin{equation} (X,Y) = (\tilde {a} \cos \eta , \tilde {b} \sin \eta ), \qquad \eta \in [0,2\pi ), \end{equation}

and so $ \tilde {a} \tilde {b}=1$ , noting that these are related to the semi-minor and semi-major axes in $(x,y)$ coordinates via $(a,b) = \sqrt {A/\pi }\,(\tilde {a},\tilde {b})$ . We also define elliptical coordinates

(4.15) \begin{equation} X = \tilde {c} \cosh \xi \cos \eta , \qquad Y = \tilde {c} \sinh \xi \sin \eta , \end{equation}

where

(4.16) \begin{equation} \tilde {c}=\sqrt {\tilde {a}^2-\tilde {b}^2}, \quad \frac {\tilde {a}}{\tilde {c}} = \cosh \xi _0, \quad \frac {\tilde {b}}{\tilde {c}} = \sinh \xi _0, \quad \xi _0 = \frac {1}{2} \log \frac {\tilde {a}+\tilde {b}}{\tilde {a}-\tilde {b}}, \end{equation}

and so the body boundary corresponds to $\xi =\xi _0$ . In elliptical coordinates, the biharmonic equation becomes

(4.17) \begin{equation} {\nabla} ^2 \left ( \frac {1}{\cosh ^2 \xi - \cos ^2 \eta } {\nabla} ^2 \psi \right ) =0, \quad {\textrm{where}} \quad {\nabla} ^2 = \frac {\partial ^2}{\partial \xi ^2}+\frac {\partial ^2}{\partial \eta ^2}. \end{equation}

The solution with no-slip, no-flux on the body and the correct symmetry properties is determined by Shintani, Umemura & Takano (Reference Shintani, Umemura and Takano1983), and given as

(4.18) \begin{equation} \psi_{\textit{in}} = E \big [ (\xi -\xi _0) \sinh \xi - \sinh \xi _0 \cosh \xi _0 \sinh \xi + \sinh ^2 \xi _0 \cosh \xi \big ] \sin \eta , \end{equation}

for some constant $E$ that is to be determined via matching with the outer solution.

4.1.3. Matching the inner and outer solutions

It is useful to note the large $R=\sqrt {X^2+Y^2}$ behaviour of the elliptical coordinates,

(4.19) \begin{equation} \xi \sim \log \frac {2R}{\tilde {c}}, \qquad \eta \sim \theta . \end{equation}

We expand the inner solution (4.18) in terms of the outer coordinates and find that the leading-order contributions are

(4.20) \begin{equation} \psi_{\textit{in}} \sim E \left [ \log \frac {2r}{\tilde {c} \varepsilon } - \xi _0 - \sinh \xi _0 \cosh \xi _0 + \sinh ^2 \xi _0 \right ] \frac {r}{\tilde {c} \varepsilon } \sin \theta . \end{equation}

Then, matching with the near-field outer solution (4.10) and comparing like terms furnishes the following two conditions:

(4.21) \begin{align}-\frac {E}{\tilde {c} \varepsilon } & = B_0,\end{align}

(4.22) \begin{align}\frac {E}{\tilde {c} \varepsilon } \left ( \log \frac {2}{\tilde {c} \varepsilon } - \xi _0 - \sinh \xi _0 \cosh \xi _0 + \sinh ^2 \xi _0 \right ) & = 1-B_0 \alpha _1 ,\end{align}

from which we find

(4.23) \begin{equation} -\frac {E}{\tilde {c} \varepsilon } = B_0 = \frac {-1}{\log \dfrac {1}{\varepsilon } -\alpha _1 - \dfrac {\tilde {b}}{\tilde {a}+\tilde {b}}+ \log \dfrac {2}{\tilde {a}+\tilde {b}} }. \end{equation}

The drag on the body (2.7) may then be obtained,

(4.24) \begin{equation} D=\frac {4 \pi }{\log \dfrac {1}{\varepsilon } -\alpha _1 - \dfrac {\tilde {b}}{\tilde {a}+\tilde {b}}+ \log \dfrac {2}{\tilde {a}+\tilde {b}} }, \end{equation}

and by setting $\tilde {a}=\tilde {b}=1$ , we can recover the result for the drag on a circle (Takaisi Reference Takaisi1955, Reference Takaisi1956),

(4.25) \begin{equation} D_{\textit{circle}}=\frac {4\pi }{\log \dfrac {1}{\varepsilon }-\alpha _1-\dfrac {1}{2}} \approx \frac {4\pi }{\log \dfrac {1}{\varepsilon }-0.9157}. \end{equation}

The optimal ellipse (for $A \ll 1$ ) is obtained by minimising (4.24) subject to $\tilde {a}\tilde {b}=1$ , which furnishes

(4.26) \begin{align} (a,b) = (\tilde {a},\tilde {b}) \sqrt {\frac {A}{\pi }} &=\left (\sqrt {\sqrt {2}+1}, \, \sqrt {\sqrt {2}-1}\,\right )\sqrt {\frac {A}{\pi }} \nonumber\\ &\approx \left (1.5538,\,0.6436\right )\sqrt {\frac {A}{\pi }}. \end{align}

This shape is shown in $(X,Y)$ coordinates in figure 8.

Figure 8. Comparison of the optimal elliptical shape (4.26) and the optimal general body shape (6.1) for small areas. The colour indicates the square of the vorticity on the body surface (relative to its mean), which is unity everywhere on the surface of the general optimal body and given by (4.29) for the ellipse.

The predictions (4.26) are compared with the numerical results for the optimal ellipse in figures 7(c) and 7(d). The leading-order aspect ratio of the optimal ellipse is $\sqrt {2}-1 \approx 0.4142$ , which is independent of $\varepsilon$ , i.e. the area; see also figure 7(bii). The minimal drag can then be calculated from (4.24) and (4.26), yielding

(4.27) \begin{equation} D_{\textit{ellipse}}=\frac {4\pi }{\log \dfrac {1}{\varepsilon }-\alpha _1-1+\dfrac {1}{\sqrt {2}}+\dfrac {1}{2}\log \big (2\sqrt {2}-2 \big )} \approx \frac {4\pi }{\log \dfrac {1}{\varepsilon }-0.8027}, \end{equation}

which shows good agreement with the numerical results for $A \ll 1$ ; see figure 7(a). The ratio of the drag on this optimal ellipse to the drag on a circle of the same area is

(4.28) \begin{equation} \frac {D_{\textit{ellipse}}}{D_{\textit{circle}}} \approx \frac {\log \dfrac {1}{\varepsilon }-0.9157}{\log \dfrac {1}{\varepsilon }-0.8027}\approx \frac {\log \dfrac {1}{A}-0.6867}{\log \dfrac {1}{A}-0.4607}. \end{equation}

Although this ratio approaches unity as $A\to 0$ , the convergence is logarithmically slow. Indeed, the drag on the optimal ellipse is 12 % less than the optimal circle for $A=0.1$ , 5 % less for $A=0.01$ and 3.5 % less for $A=0.001$ . Equation (4.28) can also be interpreted as the small circle having the same drag as the small optimal elliptical body of $25.4\,\%$ greater area.

The vorticity on the boundary of the optimal elliptical body is

(4.29) \begin{equation} \omega =\frac {2\tilde {a} \sin \eta }{\varepsilon \big(\tilde {a} ^2-\tilde {c}^2 \cos ^2 \eta \big ) \left(\log \dfrac {1}{\varepsilon } -0.8027 \right)}, \end{equation}

which is compared with the numerical result for $A=0.025$ in figure 2(e).

The body surface vorticity (4.29) vanishes at the upstream and downstream stagnation points ( $\eta =0^{{\circ}},180^{{\circ}}$ ), suggesting that the optimal general body requires a sharper apex angle. We formalise this observation later via variational analysis in § 5.1 and via apex angle analysis in § 5.2. The surface vorticity on the ellipse has local maxima at $\eta =\pm 22.5^{{\circ}}, \pm 157.5^{{\circ}}$ ; the vorticity is indicated by the colour of the ellipse in figure 8. This figure also shows the small general optimal body and comparison suggests that at the locations on the ellipse near the local maxima of surface vorticity, the body must be ‘shrunk’ to obtain the optimal general body; see also the variational argument in § 5. Likewise, the body should be ‘grown’ near the local minima of surface vorticity at the two ends of the ellipse.

4.2. Large elliptical bodies

In the regime of a large body, it is convenient to denote the body by

(4.30) \begin{equation} -a\leqslant x \leqslant a, \quad -h(x)\leqslant y \leqslant h(x), \end{equation}

where $h(x)$ is the shape of the boundary. We rescale the streamwise coordinate with

(4.31) \begin{equation} \chi = \frac {x}{a}, \end{equation}

with the body lying in $-1\leqslant \chi \leqslant 1$ . In the regime of large bodies ( $A\gg 1$ ), the numerical results of figures 7(c) and 7(d) indicate that the half-length of the optimal ellipse is proportional to the area, $a\propto A$ , whilst its maximum thickness, $h(0)$ , converges to a constant. Hence, the ratio of the channel thickness to the body length is proportional to $1/A\ll 1$ in this regime. As the two fluid gaps between the body and the channel walls are long and relatively thin, the flow is then predominantly in the horizontal direction with $v \ll u$ . The exceptions are the two adjustment regions in the vicinity of the two apices in which the streamwise velocity in the centre of the channel transitions from unity to zero (or from zero to unity at the trailing end) over a distance $x =\mathcal{O}(1)$ or equivalently $\chi = \mathcal{O}(1/a)$ (e.g. Graebel Reference Graebel1965). We neglect the adjustment behaviour in these regions and assume that, away from the body, $u=1-y^2$ , $v=0$ in $|\chi |\gt 1$ . Above and below the body, we deploy the lubrication approximation in $|\chi |\lt 1$ with

(4.32) \begin{align} u &= u_0\left (y;h(\chi )\right ) + a^{-2} u_2\left (y;h(\chi )\right ) + \mathcal{O} (a^{-4} ), \\[-9pt] \nonumber \end{align}

(4.33) \begin{align} v &= a^{-1} v_0\left (y;h(\chi )\right ) + a^{-3} v_2\left (y;h(\chi )\right ) + \mathcal{O} (a^{-5} ). \end{align}

We use the higher-order lubrication equations to calculate velocities up to and including terms of order $a^{-2}$ in the upper half of the channel, $h\lt y\lt 1$ , which are given in Appendix A (see also Tavakol et al. Reference Tavakol, Froehlicher, Holmes and Stone2017).

The drag on the body (2.10) may be written as

(4.34) \begin{equation} D = 2 \int _{-a}^a\left ( \frac {\partial u}{\partial y}- \frac {\partial ^2 u}{\partial y^2} \right )\bigg \rvert _{y=1} \, {\textrm{d}} x ={2a \int _{-1}^1\left ( \frac {\partial u}{\partial y}- \frac {\partial ^2 u}{\partial y^2} \right )\bigg \rvert _{y=1} \, {\textrm{d}} \chi ,} \end{equation}

since by assumption, $u=1-y^2$ in $|x|\gt a$ . This can be used to obtain the leading-order lubrication approximation for the drag,

(4.35) \begin{align} D= 8 a \int _{-1}^1 \frac {1+h}{(1-h)^3} \, {\textrm{d}} \chi \quad {\textrm{where}} \quad h(\chi )=b\sqrt {1-\chi ^2}. \end{align}

The vorticity is given by

(4.36) \begin{equation} {\omega } =-\frac {\partial u_0}{\partial y}+a^{-2} \left ( \frac {\partial v_0}{\partial \chi }-\frac {\partial u_2}{\partial y} \right ) + \mathcal{O}(a^{-4}), \end{equation}

and on the boundary of the body, this becomes

(4.37) \begin{align} {\omega }(y=h) = -\frac {4}{(1-h)^2}+a^{-2} \left [\frac {8}{5(1-h)^2} \left ( \frac {\partial h}{\partial \chi }\right )^2 + \frac {16}{15(1-h)} \frac {\partial ^2 h}{\partial \chi ^2} \right ] +\mathcal{O}(a^{-4}).\\[-18pt]\nonumber \end{align}

For a large elliptical body, we seek the semi-major and semi-minor axes, $a$ and $b$ , which minimise the leading-order drag (4.35) subject to fixed area $a b\pi =A$ . By writing $a=A/(b\pi )$ , this becomes a one-parameter minimisation problem and the optimal choice, obtained via numerical quadrature, is

(4.38) \begin{equation} b\approx 0.2553, \qquad a\approx 1.2469 A, \end{equation}

which gives a minimal drag of

(4.39) \begin{equation} D\approx 38.9526, \qquad a \approx 48.5705 A. \end{equation}

The optimal large elliptical body has fixed semi-minor axis and is simply extended in the streamwise direction at larger areas; see figure 7(b). The leading-order asymptotic predictions for the drag on the optimal ellipse (4.39) and the semi-minor and semi-major axes of the optimal ellipse (4.38) are compared with the numerical results in figure 7, and there is good agreement for $A \geqslant 2$ . We note that the next term in the drag on the ellipse is proportional to $\partial h/\partial \chi$ and $\partial ^2 h/\partial \chi ^2$ , and so is technically singular as $\chi \to \pm 1$ . Hence, a full and formal asymptotic treatment of this behaviour (which is beyond the scope of the present study) would require matching with a transition region near the noses when the lubrication equations start to break down, in which the flow is governed by the full Stokes equations.

The vorticity on the body (4.37) with $h=b\sqrt {1-\chi ^2}$ is compared with the numerical result for the vorticity on the optimal ellipse for $A=4$ in figure 2(f). The vorticity is singular at the apices, as expected, where the lubrication assumption breaks down (the vorticity must vanish at $\chi =\pm 1$ because each apex is a stagnation point).

5.1. Vorticity on the surface of an optimal body

In the case of uniform flow in an unbounded domain, the optimal body is the body that has constant vorticity on its boundary (Pironneau Reference Pironneau1973; Richardson Reference Richardson1995; Zabarankin Reference Zabarankin2013). Here, we show that an optimal body in Poiseuille flow has approximately (but not exactly) constant surface vorticity,

(5.1) \begin{equation} \omega ^2= {\textrm{const.}} +\mathcal{O} (b^2)\quad {\textrm{on}} \quad \partial B, \end{equation}

where $b$ is the greatest half-thickness of the body and $\partial B$ is the boundary of the body. There is a change of sign of $\omega$ across the centreline. This result is derived as follows.

First, from no-slip on the body boundary, $u_i=0$ on $\partial B$ , we find that

(5.2) \begin{equation} {0=\int _{\partial B}u_i \sigma _{\textit{ij}} n_j \,{\textrm{d}}S.} \end{equation}

We then consider a large rectangle $R=[-L,L]\times [-1,1]$ with $L\gg 1$ where the boundary of $R$ (which does not include the body boundary) is denoted by $\partial R$ and the fluid domain is denoted by $R\setminus B$ . Upon applying the divergence theorem, we obtain the relation

(5.3) \begin{equation} {0=-\iint _{R\setminus B} \frac {\partial u_i}{\partial x_j}\sigma _{\textit{ij}} \, {\textrm{d}}x\,{\textrm{d}}y - \int _{\partial R}u_i \sigma _{\textit{ij}} n_j \,{\textrm{d}}S,} \end{equation}

where the unit normal points into the fluid domain. We use mass continuity and symmetry in the indices to recast the first integral as the viscous dissipation,

(5.4) \begin{equation} {E_t \equiv \frac {1}{2} \iint _{R\setminus B} \left ( \frac {\partial u_i}{\partial x_j}+ \frac {\partial u_j}{\partial x_i} \right )^2 \, {\textrm{d}}x\,{\textrm{d}}y=\iint _{R\setminus B} \frac {\partial u_i}{\partial x_j}\sigma _{\textit{ij}} \, {\textrm{d}}x\,{\textrm{d}}y .} \end{equation}

The integrand in the second integral in (5.3) vanishes on $y=\pm 1$ owing to no-slip, which leaves the far-field contribution from the ends of the rectangle at $x=\pm L$ , where $u=1-y^2$ , $v=0$ and $\partial p/\partial y=0$ so that

(5.5) \begin{equation} {- \int _{\partial R}u_i \sigma _{\textit{ij}} n_j \,{\textrm{d}}S=\frac {4}{3} \Delta p,} \end{equation}

where $\Delta p=[-p]_{-L}^L$ is the pressure decrease. Equation (5.3) then relates the viscous dissipation in the channel to the pressure drop,

(5.6) \begin{equation} {\frac {4}{3} \Delta p = E_t.} \end{equation}

In an unbounded domain, the drag on a body in Stokes flow is equal to the viscous dissipation in the fluid (e.g. Montenegro-Johnson & Lauga Reference Montenegro-Johnson and Lauga2015), but in Poiseuille flow in a channel, the drag on the body also depends on the second moment of the $x$ traction about the centreline of the body, which we denote by $J$ . A similar manipulation to that already performed can be used to show that

(5.7) \begin{equation} J\equiv \int _{\partial B} y^2 \sigma _{1j} n_j \,{\textrm{d}}S = \frac {2}{3} \Delta p + 2 \int _{-L}^{L} \frac {\partial u}{\partial y} \bigg \rvert _{y=1} \, {\textrm{d}} x + \frac {16}{3} L. \end{equation}

Combining (5.6) and (5.7) with (2.9) furnishes

(5.8) \begin{equation} D=E_t -\frac {16}{3}L + J. \end{equation}

We now perturb the body via ${\boldsymbol{x}} \to {\boldsymbol{x}} + \delta c(s) {\boldsymbol{n}}$ with $\delta \ll 1$ , where $s$ denotes arclength, and we denote the small change in velocity by $\delta \tilde {\boldsymbol{u}}$ and in the stress tensor by $\delta \tilde {\boldsymbol{\mathsf{\mathit{\sigma }}}}$ . The perturbation in the drag (5.8) becomes

(5.9) \begin{equation} \delta \tilde {D} = \delta \tilde {E}_t + \delta \tilde {J}. \end{equation}

Here, the perturbations $\tilde {E}_t$ and $\tilde {J}$ are given by straightforward variational analysis (details of the calculations are given in Appendix B; see also Pironneau Reference Pironneau1973; Richardson Reference Richardson1995; Zabarankin Reference Zabarankin2013; Montenegro-Johnson & Lauga Reference Montenegro-Johnson and Lauga2015),

(5.10) \begin{equation} \tilde {E}_t = 2 \int _{\partial B} c(s) \omega ^2 \, {\textrm{d}}S, \qquad \tilde {J}=\int _{\partial B} y^2 \left (\tilde {\sigma }_{1j} + \sigma _{1j} c(s) \right ) n_j \, {\textrm{d}} S. \end{equation}

We introduce a cost function of the form

(5.11) \begin{equation} \mathcal{F} = D + \lambda \left (A_0-A \right ), \end{equation}

where $\lambda$ is the Lagrange multiplier enforcing fixed area $A=A_0$ , and we find that

(5.12) \begin{equation} 0= \tilde {\mathcal{F}} = \int _{\partial B} 2 c \omega ^2 + y^2 \left (\tilde {\sigma }_{1j} + \sigma _{1j} c \right ) n_j - \lambda c \, {\textrm{d}} S. \end{equation}

Hence, as $y\to 0$ (i.e. in the limit of small bodies), the vorticity will be constant to leading order. Specifically, we can rewrite (5.12) as

(5.13) \begin{equation} 0= \int _{\partial B} 2 c \omega ^2 - \lambda c \, {\textrm{d}} S + \mathcal{O}(b^2). \end{equation}

The vorticity on the body hence satisfies

(5.14) \begin{equation} \omega ^2 = \omega _0^2 + \mathcal{O}(b^2), \end{equation}

for some non-zero constant $\omega _0$ . Over all areas, we show later (see § 6.2) that the optimal body has greatest thickness $b \leqslant 0.215$ and so $b^2\leqslant 0.046$ , and that the vorticity squared satisfies $\omega _0^2 \geqslant 42.189$ . Thus, $b^2$ ends up being numerically small when compared with $\omega _0^2$ . Moreover, it can be shown that the constant $\omega _0$ is of $\mathcal{O}(1/b)$ as $A \to 0$ (i.e. as $b\to 0$ ) (Richardson Reference Richardson1995). Hence, the vorticity nowhere vanishes nor blows up on the body, which is important for analysing the apex; see § 5.2. The approximate condition (5.1) becomes exact in the limit $A \to 0$ .

A shape with smooth boundary everywhere will always have a stagnation point with vanishing vorticity at the location on the boundary where the separating streamline meets the object. A stagnation point violates the optimality condition and so the optimal shape must have a sharp corner.

These results confirm that $\omega ^2$ is approximately (but not exactly) constant over the body surface. Equation (5.12) indicates that given a small non-optimal body shape (say an ellipse), the drag can generally be decreased by shrinking the body in regions where the vorticity on the boundary is larger than its mean and expanding the body in regions where the vorticity is small; e.g. figure 8.

The drag can be written alternatively as

(5.15) \begin{equation} D= \frac {3}{2}E_t + 2 \int _{-L}^{L} \frac {\partial u}{\partial y} \bigg \rvert _{y=1} \, {\textrm{d}} x. \end{equation}

This description decomposes the drag into a contribution that is proportional to the dissipation and one that is proportional to the friction on the walls. For small bodies, the perturbation to the second term disappears i.e. $\tilde {D}$ only has a contribution from the vorticity on the body squared (first term in (5.12)) for $A \ll 1$ . Furthermore, if the channel wall boundary conditions were no flux and free slip ( $v=\partial u/\partial y=0$ ) rather than no flux and no slip ( $u=v=0$ ), then the condition of constant vorticity on the optimal body would be exact for all areas.

5.2. Optimal nose angle at the front and back is $102.5^\circ$

We now investigate the nose angle of the lowest-drag body in a channel. Whilst the argument we will use to determine the angle is based upon a local analysis, it relies upon the surface vorticity result from § 5.1, which takes into account the entire flow domain, including the channel walls.

Consider a polar coordinate system centred at the trailing edge of the body $(a,0)$ , so that

(5.16) \begin{equation} x= a +\tilde {r} \cos \tilde {\theta }, \qquad y = \tilde {r} \sin \tilde {\theta }. \end{equation}

Locally to this trailing edge, the streamfunction takes the form (Moffatt Reference Moffatt1964)

(5.17) \begin{equation} \psi =\tilde {r}^\lambda f_{\lambda }(\tilde {\theta }), \qquad f_{\lambda }(\pm \beta ) = f^{\prime}_{\lambda }(\pm \beta ) = 0. \end{equation}

where $\beta =\pi -\alpha$ and $\alpha$ is the half-angle at the leading and trailing edges; see figure 1. The vorticity on the body (locally to the trailing edge) is $-{\nabla} ^2 \psi \sim \tilde {r}^{\lambda -2}$ . For the vorticity to be non-zero and non-singular as $\tilde {r} \to 0$ (from (5.1)), we require $\lambda =2$ . The $\tilde {\theta }$ dependence is then

(5.18) \begin{equation} f_{2}(\tilde {\theta }) = B_1 \sin 2 \tilde {\theta } +B_2 \tilde {\theta }, \end{equation}

to satisfy the biharmonic equation, where $B_1$ and $B_2$ are arbitrary constants. For a non-trivial solution that satisfies the boundary conditions, we require (Richardson Reference Richardson1995)

(5.19) \begin{equation} 2 \beta = \tan 2 \beta , \end{equation}

which has the solution

(5.20) \begin{equation} \beta \approx 128.727^{\textrm{o}}, \quad \alpha \approx 51.273^{\textrm{o}}. \end{equation}

Therefore, the optimal nose angle is $2\alpha \approx 102.5^{\textrm{o}}$ . Condition (5.20) is verified numerically in § 5.3 by computing the optimal shapes without using this constraint and comparing the results to constrained optimisation with the nose angle imposed; see also figures 5 and 6.

5.3. Applying constraints to the numerical optimisation

Our previous results can be used to reduce the computational effort of the numerical method (§ 3.2) for obtaining the optimal shape for $N\geqslant 2$ . Namely, by using the results that the gradient at $x=0$ vanishes owing to symmetry and that the gradient at the nose is fixed (see § 5.2), the number of free parameters in the shape optimisation is reduced by two. The condition that the noses are non-smooth and have fixed angle will be imposed as a constraint.

For the form of the shape given by (3.1), the derivatives with respect to $t$ are given by

(5.21) \begin{align} \frac {{\textrm{d}} x}{{\textrm{d}} t} & = C \sum _{j=0}^N \left ( -a_j \frac {\pi }{2} T_j(2t-1) \sin \frac {\pi t}{2} + 2 a_j j U_{j-1} (2t-1) \cos \frac {\pi t}{2} \right )\!, \end{align}

(5.22) \begin{align} \frac {{\textrm{d}} y}{{\textrm{d}} t} & = C \sum _{j=0}^N \left ( b_j \frac {\pi }{2} T_j(2t-1) \cos \frac {\pi t}{2} + 2 b_j j U_{j-1} (2t-1) \sin \frac {\pi t}{2} \right )\!, \end{align}

where $U_j$ are Chebyshev polynomials of the second kind.

At the nose ( $t=0$ ), the gradient of the boundary is fixed and so

(5.23) \begin{equation} - \tan \alpha = \frac {{\textrm{d}} y}{{\textrm{d}} x}= \frac {\sum _{j=0}^N b_j \dfrac {\pi }{4} (-1)^{j} }{\sum _{j=1}^N a_j j^2 (-1)^{j-1} }, \end{equation}

where $\alpha$ is the interior half-angle, given by (5.20). At the midpoint ( $t=1$ ), the gradient of the boundary vanishes giving

(5.24) \begin{equation} 0= \sum _{j=1}^N b_j j^2. \end{equation}

Hence, for $N\geqslant 2$ , we impose

(5.25) \begin{equation} b_N = -\frac {1}{N^2} \sum _{j=1}^{N-1} b_j j^2 \end{equation}

as well as

(5.26) \begin{equation} a_N = \frac {\pi }{4} \cot \alpha \frac {(-1)^N}{N^2} \sum _{j=0}^N b_j (-1)^{j}- \frac {(-1)^N}{N^2} \sum _{j=1}^{N-1} a_j j^2 (-1)^{j}. \end{equation}

The constraints (5.25)–(5.26) reduce the number of free parameters to be determined in the optimisation process of § 3.2 from $2N+1$ to $2N-1$ and so the computational effort is lower. Differences between the optimisation results with and without these constraints for $N\geqslant 4$ are imperceptible; see figures 5 and 6. The drag on the optimal body for $N=2,4,6$ is compared with the drag on the optimal ellipse ( $N=0$ ) for different areas in figure 9 (using the constraints (5.25)–(5.26) for $N\geqslant 2$ ). The change in drag between $N=4$ and $N=6$ was found to be less than $0.002\,\%$ for $A \leqslant 4$ .

Figure 9. Drag reduction of optimal bodies relative to the optimal ellipse for different areas. The ellipse corresponds to the shape (3.1) with $N=0$ . Increasing $N$ is associated with more degrees of freedom in the shape optimisation. For $N\gt 0$ , the optimal body (and the associated drag) was calculated using the constrained approach described in § 5.3.

6.1. General optimal body of small area

Numerical results for the flow past the optimal body with $A=0.1$ are shown in figure 3. Richardson (Reference Richardson1995) calculated the optimal shape in an unbounded uniform Stokes flow, which is also the shape for which the vorticity has constant magnitude on its boundary. To construct the exact solution, Richardson (Reference Richardson1995) conformally mapped the optimal body onto the unit circle and rewrote the three boundary conditions on the body boundaries (no-slip, no flux and constant vorticity) for this new domain. Richardson used a Goursat representation, which leaves two analytic functions (and the body shape) that must be determined by using the boundary conditions on the body surface and in the far-field. After various manipulations, the problem reduces to determining an analytic function that satisfies the correct jump conditions.

In the limit of small area, the optimal body in a channel has exactly constant surface vorticity (see § 5.1) and so it must have the same shape as Richardson’s solution. However, the drag on the optimal body in channel flow is different owing to wall effects. This difference can be understood by considering the asymptotic regions in the mathematical problem: in both the free-space and the channel flow configurations, the inner region for the flow around the optimal body has the same form, but the outer region is different (it must either account for inertia or it must account for the channel walls); see also § 1.

To understand the regime of confined flow past the small optimal body, we consider an inner coordinate system given by (4.12), with the body having area $A=\pi \varepsilon ^2$ and where we define $Z=X+i Y$ . The body is given by

(6.1) \begin{equation} (X,Y) = \left (\mathcal{X}(s),\mathcal{Y}(s) \right ), \end{equation}

and the shape of the optimal body in the positive quadrant is (Richardson Reference Richardson1995)

(6.2) \begin{equation} \mathcal{X}(s) = \sqrt {2} \int _0^\infty \frac {f'(k/2)}{f(k)} \sin (ks) \, {\textrm{d}}k, \quad \mathcal{Y}(s) = \sqrt {2} \int _0^\infty \frac {f(k/2)}{f(k)} \cos (ks) \, {\textrm{d}}k, \end{equation}

where $0\leqslant s \lt \infty$ parametrises the body and $f(k)=k\cosh k - \sinh k$ . The shape in the other three quadrants is obtained by symmetry. The area in $(X,Y)$ coordinates is $\pi$ and so the area in the original $(x,y)$ coordinates is $A$ . Equation (6.2) corresponds to the optimal shape with (in Richardson’s notation) dimensionless length scale $C=\pi /(2\sqrt {2})$ . The half-extent and half-width of this body are

(6.3) \begin{equation} \max \left [\mathcal{X}(s)\right ] =\frac {3 \pi }{4 \sqrt {2}} \approx 1.66608, \quad \max \left [\mathcal{Y}(s)\right ] \approx 0.67065. \end{equation}

The optimal shape (6.1) is shown in figure 8, where it is compared with the optimal ellipse for $A \ll 1$ .

The solution for the outer flow away from the optimal body has the same functional form as that in § 4.1.1 for the optimal elliptic body, with the difference in shape yielding different values for the matching constant $B_0$ . To determine its value, we note that the near-field expansion of the outer flow (4.10) is equivalent to the following complexified velocity:

(6.4) \begin{equation} u+ i v = -B_0 \log r + \frac {B_0}{2} e^{2 i \theta } -\left (\alpha _1+\frac {1}{2} \right ) B_0 +1, \quad {\textrm{as}} \quad r\to 0. \end{equation}

We now use (6.4) to determine $B_0$ by matching with the inner solution for Richardson’s optimal shape.

For the inner solution, we write the velocity using the Goursat representation as (Muskhelishvili Reference Muskhelishvili1977; Richardson Reference Richardson1995)

(6.5) \begin{equation} U+i V = \phi (Z) - Z \overline {\phi '(Z)} - \overline {\chi (Z)}. \end{equation}

To agree with the outer behaviour (6.4), we write the two analytic functions as (cf. Richardson Reference Richardson1995)

(6.6) \begin{equation} \phi (Z) = \frac {-B_0}{2} \log Z + \phi _0 (Z), \qquad \chi (Z) = \frac {B_0}{2} \log Z +Q + \chi _0(Z), \end{equation}

with $\phi _0 (Z),\chi _0(Z) \to 0$ as $|Z| \to \infty$ , and so the constant $Q$ is obtained via matching (6.5) in the large- $|Z|$ limit with (6.4) to deduce

(6.7) \begin{equation} Q = -B_0 \left (\log \frac {1}{\varepsilon } -\alpha _1-\frac {1}{2} \right )-1. \end{equation}

Then, the solution of the inner flow problem (in which the magnitude of the body’s surface vorticity is constant) furnishes (see Richardson Reference Richardson1995)

(6.8) \begin{equation} Q=\frac {-B_0}{2} \varLambda , \end{equation}

where

(6.9) \begin{align} \varLambda &=2\log \frac {\pi }{4} +\frac {4}{\pi } G -1 +\int _0^{\infty } \frac {\left (s^2-\sinh s \tanh s\right )\left (s \sinh s - \cosh s +1\right )}{s^2 \sinh s \left (s \cosh s -\sinh s\right )} \, {\textrm{d}} s \nonumber\\ &\approx -0.23586, \end{align}

and

(6.10) \begin{equation} G=\sum _{n=0}^\infty \frac {(-1)^n}{(2n+1)^2}\approx 0.915965 \end{equation}

is Catalan’s constant.

Equating (6.7) and (6.8) determines the constant

(6.11) \begin{equation} B_0 =\frac {-1}{\log \dfrac {1}{\varepsilon } -\alpha _1-\dfrac {1}{2}-\dfrac {\varLambda }{2}} \approx \frac {-1}{\log \dfrac {1}{\varepsilon } -0.7978}. \end{equation}

The drag is therefore

(6.12) \begin{equation} D \approx \frac {4 \pi }{\log \dfrac {1}{\varepsilon } -0.7978}, \end{equation}

which is a slight improvement on the optimal ellipse (4.27) with the drag proportional to $D=\mathcal{O} ( (\log \frac {1}{A} )^{-1} )$ as $A \to 0$ for both shapes. The prediction (6.12) is compared with the numerical results for the drag in figure 10(a).

Figure 10. (a) Drag on the general optimal body as a function of area. Black lines show the numerical results, the blue dash-dotted line shows the small-area result (6.12) and red dashed lines shows the large-area result (6.28). (b) Aspect ratio of the general optimal body.

The aspect ratio of this small optimal shape is $0.4025$ and that of the small optimal ellipse is $0.4142$ ; see figure 8 and also figure 10(b) for a comparison with the numerical results. The ratio of the drag on this optimal body to that on the optimal ellipse with the same area is

(6.13) \begin{equation} \frac {D_{\textit{optimal}}}{D_{\textit{ellipse}}} \approx \frac {\log \dfrac {1}{\varepsilon }-0.8027}{\log \dfrac {1}{\varepsilon }-0.7978}\approx \frac {\log \dfrac {1}{A}-0.4607}{\log \dfrac {1}{A}-0.4509}. \end{equation}

Equation (6.13) can interpreted as the small optimal ellipse having the same drag as the small optimal general body of $0.98\,\%$ greater area. For $A=0.1$ , the optimal body has 0.53 % less drag than the optimal ellipse of the same area; for $A=0.01$ , it has 0.24 % less drag; and for $A=0.001$ , it has 0.15 % less drag than the optimal ellipse (see also figure 9).

Richardson’s shape (shown as red dashed lines) is compared with the numerically computed optimal shape in figure 5(a–d) and we see excellent agreement for $A=0.025$ and reasonable agreement for $A=0.1$ . Figure 5(a–d) indicates that the first-order correction to Richardson’s shape, which accounts for modest wall effects, would make the optimal shape more slender.

The magnitude of the vorticity is constant (to leading order) on this optimal small-area body and is given by

(6.14) \begin{equation} |\omega | = \frac {-B_0}{2} \frac {\pi }{C \varepsilon } = \frac {\sqrt {2}}{\varepsilon \left (\log \dfrac {1}{\varepsilon } - 0.7978 \right )}. \end{equation}

Equation (6.14) can be compared with the vorticity magnitude on a small optimal ellipse (4.29), which varies with angle $\eta$ , to estimate the four regions in which the ellipse has larger vorticity than the optimal body and so protrudes outside the optimal body (see figure 8). In the small- $\varepsilon$ limit, the $\log 1/\varepsilon$ terms in the denominators dominate the small difference in the numerical constants $0.8027$ and $0.7978$ , from which we may deduce that the regions of relatively higher body vorticity are $\varepsilon$ -independent in the limit. Using this limit, in the positive quadrant, the ellipse has higher vorticity in the region given by angles that satisfy

(6.15) \begin{equation} 14^{{\circ}} \lt \eta \lt 59^{{\circ}}, \end{equation}

with the other three regions obtained by symmetry.

6.2. Large bodies ( $A \gg 1$ )

Over most of their extent, optimal large bodies have constant thickness; see figure 4. Therefore, as a first approximation, we seek the rectangular body of length $2a$ and half-thickness $h_0$ that minimises the leading-order drag (according to the lubrication approximation; see § 4.2),

(6.16) \begin{equation} D= 8 a \int _{-1}^1 \frac {1+h}{(1-h)^3} \, {\textrm{d}} \chi = \frac {16a(1+h_0)}{(1-h_0)^3} \end{equation}

subject to fixed area, $4ah_0=A$ . This minimisation problem reduces to a single quadratic equation for $h_0$ with physical solution

(6.17) \begin{equation} h_0=\frac {\sqrt {7}-2}{3}\approx 0.2153, \end{equation}

which is independent of $A$ and is slightly thinner at $x=0$ than the optimal large ellipse, defined in (4.38). This prediction for the thickness shows excellent agreement with the numerical results away from the ends of the body; see the magenta dash-dotted line in figure 6(c). At the ends of the body, a more subtle treatment is required to account for $h=0$ and the nose angle, which is discussed as follows.

The prediction (6.17) corresponds to a half-length of

(6.18) \begin{equation} a=\frac {3A}{4(\sqrt {7}-2)}\approx 1.1614 A, \end{equation}

and a minimal drag,

(6.19) \begin{equation} D=\left (\frac {536+208\sqrt {7}}{27}\right )a \approx 40.2339 a, \qquad D=\left (\frac {632+238\sqrt {7}}{27}\right )A\approx 46.7292A, \end{equation}

Comparing the minimal drag result (6.19) for the general optimal body to the equivalent result for the optimal ellipse (4.39) suggests a drag reduction of 3.8 %, which is significantly larger than the drag reduction for small bodies. The vorticity on the surface of this optimal body (away from the ends) is predicted to be

(6.20) \begin{equation} |\omega | = \frac {4}{(1-h_0)^2} = \frac {32+10\sqrt {7}}{9}\approx 6.4953, \end{equation}

which is compared to the numerical results in figure 11.

Figure 11. Magnitude of the vorticity on the surface of the optimal body for $A=0.1, 1, 4$ . The vorticity is approximately (but not exactly) constant; see § 5.1. The red dashed line shows the leading-order prediction for large $A$ , (6.20), whilst the magenta dash-dotted line is the improved prediction from integrating (6.25) and using (4.37).

To account for the ends of the body, we introduce the scaled coordinate $\zeta =1-a\chi =1-x$ (see § 4.2) and rewrite the drag as

(6.21) \begin{align} D[h]&=\int _0^a F(h,h',h'') \, {\textrm{d}} \zeta \end{align}

(6.22) \begin{align} &= 16 \int _0^{a} \frac {1+h}{(1-h)^3} + \frac {h-4}{15(1-h)^2} \frac {\partial ^2 h}{\partial \zeta ^2} + \frac {3(h+1)}{5(1-h)^3} \left ( \frac {\partial h}{\partial \zeta }\right )^2 \, {\textrm{d}} \zeta , \\[9pt] \nonumber \end{align}

where a prime denotes differentiation with respect to $\zeta$ and we have incorporated the next order term from the lubrication approximation; see Appendix A. Combining (6.21) with the fixed area constraint,

(6.23) \begin{equation} A[h]=\int _0^{a} G(h) \, {\textrm{d}} \zeta =4\int _0^{a} h \, {\textrm{d}} \zeta , \end{equation}

and writing $J=F-\lambda G$ (where $\lambda$ is a Lagrange multiplier) furnishes the following Euler–Lagrange equation for the optimal shape

(6.24) \begin{equation} \frac {\partial J}{\partial h} -\frac {{\textrm{d}}}{{\textrm{d}} \zeta }\frac {\partial J}{\partial h'}+\frac {{\textrm{d}}^2}{{\textrm{d}} \zeta ^2}\frac {\partial J}{\partial h''}=0, \end{equation}

which must be solved subject to $h(a)=h_0$ , $h'(a)=h''(a)=0$ and $h(0)=0$ . The conditions at $\zeta =a$ determine $\lambda$ and then (6.24) becomes

(6.25) \begin{equation} 15 \left [2+h- \frac {(2+h_0)(1-h)^4}{(1-h_0)^4} \right ]= 4\left ( \frac {\partial h}{\partial \zeta }\right )^2 (2h+7) +8\frac {\partial ^2 h}{\partial \zeta ^2} (1-h)(2+h). \end{equation}

We integrate this ordinary differential equation for $h(\zeta )$ numerically, shooting from $\zeta =0$ and selecting $h'(0)$ to ensure the boundary conditions as $\zeta \to a\gg 1$ are satisfied. We find that $h'(0)=0.969$ , which is a more acute nose angle than predicted by the complete nose analysis. Of course, the present analysis is not strictly valid at the nose because the lubrication approximation breaks down. Nonetheless, the solution, $h(\zeta )$ , does an excellent job at predicting the shape; see the red dashed line in figures 6(c) and 12. It can also be used to approximate the vorticity on the body (4.37) (see figure 11 where it is compared with the numerical results for $A=4$ ).

Figure 12. Comparison between the optimal shape computed with unconstrained optimisation (blue line) and the optimal shape obtained from lubrication theory (red dashed line from (6.25)) for $A=2$ (the zoomed out version is shown in figure 6 c).

With the numerical solution to (6.25), we can calculate the area of this optimal body and its accompanying drag from (6.23) and (6.21),

(6.26) \begin{equation} A\approx 0.8610 a -0.2162, \qquad D\approx 40.2339 a -6.2621, \end{equation}

where the leading-order terms are identical to the result for the optimal rectangle and, as expected, the rounding of the noses has reduced both the area and the drag (for a fixed value of the length $a$ ). The aspect ratio of the optimal body is given by

(6.27) \begin{equation} \frac {b}{a} \approx \frac {0.1853}{A+0.2162}, \end{equation}

which is compared with the numerical results in figure 10(b). The drag can be written in terms of the area as

(6.28) \begin{equation} D \approx 46.7292A + 2.4365, \end{equation}

where the first term is as for the rectangle and the second term reflects the less efficient contribution from the noses. This prediction is compared with the numerical results in figure 10(a) and shows excellent agreement even for only moderate values of $A$ . As for the ellipse, the drag is proportional to the area with $D = \mathcal{O}(A)$ as $A \to \infty$ , but with a different prefactor.