2.1 Dimensionless system

As motivated in § 1, we are interested in a regime where fluid inertia is negligible. Thus, the fluid flow is governed by the Stokes and continuity equations. We define the aspect ratio of the rod to be $\unicode[STIX]{x1D716}=R_{i}/L\ll 1$ and non-dimensionalise the axial coordinate by $L$ and the cross-sectional coordinates by $R_{i}=\unicode[STIX]{x1D716}L$. We assume a non-zero constant flux $Q$ is applied at $z=0$ and thus, non-dimensionalise the axial velocity by characteristic velocity scale $W=Q/R_{i}^{2}$. We choose the pressure scaling $P=(W\unicode[STIX]{x1D707}/L\unicode[STIX]{x1D716}^{2})$ to balance the viscous terms in the Stokes equations and the cross-sectional force and torque are scaled by $\unicode[STIX]{x1D716}L\unicode[STIX]{x1D707}W$ and $\unicode[STIX]{x1D716}^{2}L^{2}\unicode[STIX]{x1D707}W$, respectively. The dimensionless Stokes equations are given in component form by

(2.2a)$$\begin{eqnarray}\displaystyle & \displaystyle -p_{x}+\unicode[STIX]{x1D716}^{2}(u_{xx}+u_{yy})+\unicode[STIX]{x1D716}^{4}u_{zz}=0, & \displaystyle\end{eqnarray}$$

(2.2b)$$\begin{eqnarray}\displaystyle & \displaystyle -p_{y}+\unicode[STIX]{x1D716}^{2}(v_{xx}+v_{yy})+\unicode[STIX]{x1D716}^{4}v_{zz}=0, & \displaystyle\end{eqnarray}$$

(2.2c)$$\begin{eqnarray}\displaystyle & \displaystyle -p_{z}+w_{xx}+w_{yy}+\unicode[STIX]{x1D716}^{2}w_{zz}=0, & \displaystyle\end{eqnarray}$$

(2.2d)$$\begin{eqnarray}\displaystyle & \displaystyle u_{x}+v_{y}+w_{z}=0, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\iint _{\unicode[STIX]{x1D6FA}}w\,\text{d}\unicode[STIX]{x1D6FA}=1,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is the two-dimensional region of the cross-section between the cylinders.

No-slip conditions at the surfaces of the inner and outer cylinders give the dimensionless boundary conditions

(2.4a-c)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u={\dot{X}}-\dot{\unicode[STIX]{x1D6E9}}(y-Y),\\ v={\dot{Y}}+\dot{\unicode[STIX]{x1D6E9}}(x-X),\\ w=0,\end{array}\right\}\quad \text{on}~\unicode[STIX]{x1D6E4}_{i}\end{eqnarray}$$

and

(2.5)$$\begin{eqnarray}\boldsymbol{u}=\mathbf{0},\quad \text{on}~\unicode[STIX]{x1D6E4}_{o}.\end{eqnarray}$$

To calculate $\dot{\boldsymbol{X}}$ and $\dot{\unicode[STIX]{x1D6E9}}$, where dots denote derivatives with respect to time, we consider conservation of linear and angular momentum for the inner rod which gives, in dimensionless form,

(2.6a,b)$$\begin{eqnarray}\ddot{\boldsymbol{X}}=\unicode[STIX]{x1D6FC}\boldsymbol{F}(\boldsymbol{X},\unicode[STIX]{x1D6E9},\dot{\boldsymbol{X}},\dot{\unicode[STIX]{x1D6E9}}),\quad \ddot{\unicode[STIX]{x1D6E9}}=\unicode[STIX]{x1D6FC}^{\prime }\unicode[STIX]{x1D70F}(\boldsymbol{X},\unicode[STIX]{x1D6E9},\dot{\boldsymbol{X}},\dot{\unicode[STIX]{x1D6E9}}),\end{eqnarray}$$

where $\boldsymbol{F}$ and $\unicode[STIX]{x1D70F}$ are the hydrodynamic force in the cross-section and $z$-component of the torque exerted on the rod due to the motion of the fluid, obtained by integrating contributions from the viscous pressure forces. Because the cylinders remain coaxial, $\unicode[STIX]{x1D70F}$ and $\boldsymbol{F}$ are the only non-zero force and torque components. The dimensionless constants in (2.6a,b) are defined as

(2.7a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D707}L^{2}/mW,\quad \unicode[STIX]{x1D6FC}^{\prime }=R_{i}^{2}\unicode[STIX]{x1D707}L^{2}/IW,\end{eqnarray}$$

where $m$ and $I$ are the rod mass and the moment of inertia about the rod axis, respectively, and $\unicode[STIX]{x1D6FC}^{\prime }$ is the same size as $\unicode[STIX]{x1D6FC}$ as $I$ scales with $R_{i}^{2}m$. We note that $\unicode[STIX]{x1D6FC}$ can be rearranged to relate the fluid density to the material density of the inner rod multiplied by the reduced Reynolds number. As discussed in the introduction, we will consider the regime where this parameter is $O(1)$. As dimensionless initial conditions for (2.6) we prescribe

(2.8a-d)$$\begin{eqnarray}\boldsymbol{X}(0)=\boldsymbol{{\mathcal{X}}},\quad \unicode[STIX]{x1D6E9}(0)=\unicode[STIX]{x1D6F7},\quad \dot{\boldsymbol{X}}(0)=U\boldsymbol{\ell },\quad \dot{\unicode[STIX]{x1D6E9}}(0)=\unicode[STIX]{x1D714},\end{eqnarray}$$

where $\boldsymbol{{\mathcal{X}}}$ is the initial location of the rod’s axis in the $(x,y)$ plane, $\unicode[STIX]{x1D6F7}$ is its initial orientation angle, $U$ is its initial velocity in a direction given by $\boldsymbol{\ell }$, a unit vector in the $(x,y)$-plane and $\unicode[STIX]{x1D714}$ is its initial angular velocity about the $z$-axis. We assume $\boldsymbol{{\mathcal{X}}}$, $U$, $\unicode[STIX]{x1D6F7}$, and $\unicode[STIX]{x1D714}$ are all $O(1)$; due to the chosen non-dimensionalisations, this requires the initial imposed rod velocities to be $\unicode[STIX]{x1D716}$ times smaller than the axial velocity.

2.2 Asymptotic analysis

We take $\unicode[STIX]{x1D716}^{2}\ll 1$ and seek expansions to our dimensionless variables of the forms

(2.9a,b)$$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{u}_{0}+\unicode[STIX]{x1D716}^{2}\boldsymbol{u}_{1}+\cdots ,\quad p=p_{0}+\unicode[STIX]{x1D716}^{2}p_{1}+\cdots ,\end{eqnarray}$$

(2.9c,d)$$\begin{eqnarray}\boldsymbol{X}=\boldsymbol{X}_{0}+\unicode[STIX]{x1D716}^{2}\boldsymbol{X}_{1}+\cdots ,\quad \unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D6E9}_{0}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6E9}_{1}+\cdots .\end{eqnarray}$$

$O(1)$

(2.10a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{\bot }^{2}w_{0}=\frac{\text{d}p_{0}}{\text{d}z},\quad \unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{u}_{0}^{\bot }+w_{0z}=0,\end{eqnarray}$$

where $\boldsymbol{u}_{0}^{\bot }=u_{0}\boldsymbol{i}+v_{0}\,\boldsymbol{j}$, and that $\text{d}p_{0}/\text{d}z$ is independent of $x$ and $y$ follows from (2.2a) and (2.2b) at leading order. To determine the cross-sectional velocity components, $u_{0}$ and $v_{0}$, we consider (2.2a) and (2.2b) at $O(\unicode[STIX]{x1D716}^{2})$, and obtain the two-dimensional Stokes equations

(2.11a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{\bot }^{2}u_{0}=p_{1x},\quad \unicode[STIX]{x1D6FB}_{\bot }^{2}v_{0}=p_{1y}.\end{eqnarray}$$

The leading-order flux condition is

(2.12)$$\begin{eqnarray}\iint _{\unicode[STIX]{x1D6FA}}w_{0}\,\text{d}\unicode[STIX]{x1D6FA}=1,\end{eqnarray}$$

and leading-order boundary conditions are

(2.13a-c)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u_{0}={\dot{X}}_{0}-\dot{\unicode[STIX]{x1D6E9}}_{0}(y-Y_{0}),\\ v_{0}={\dot{Y}}_{0}+\dot{\unicode[STIX]{x1D6E9}}_{0}(x-X_{0}),\\ w_{0}=0,\end{array}\right\}\quad \text{on}~\unicode[STIX]{x1D6E4}_{i}\end{eqnarray}$$

and

(2.14a,b)$$\begin{eqnarray}\boldsymbol{u}_{0}^{\bot }=\mathbf{0}\quad \text{and}\quad w_{0}=0,\quad \text{on }~\unicode[STIX]{x1D6E4}_{o}.\end{eqnarray}$$

The equations of motion for the rod, equation (2.6), read at leading order

(2.15a,b)$$\begin{eqnarray}\ddot{\boldsymbol{X}}_{0}=\unicode[STIX]{x1D6FC}\boldsymbol{F}_{0}(\boldsymbol{X}_{0},\unicode[STIX]{x1D6E9}_{0},\dot{\boldsymbol{X}}_{0},\dot{\unicode[STIX]{x1D6E9}}_{0}),\quad \ddot{\unicode[STIX]{x1D6E9}}_{0}=\unicode[STIX]{x1D6FC}^{\prime }\unicode[STIX]{x1D70F}_{0}(\boldsymbol{X}_{0},\unicode[STIX]{x1D6E9}_{0},\dot{\boldsymbol{X}}_{0},\dot{\unicode[STIX]{x1D6E9}}_{0}),\end{eqnarray}$$

where the leading-order expressions for the hydrodynamic forces and torque on the rod are

(2.16a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{0}=\int _{0}^{1}\left[\oint _{\unicode[STIX]{x1D6E4}_{i}}\unicode[STIX]{x1D748}_{0}\boldsymbol{n}_{i}\,\text{d}s\right]\text{d}z, & \displaystyle\end{eqnarray}$$

(2.16b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{0}=\int _{0}^{1}\left[\oint _{\unicode[STIX]{x1D6E4}_{i}}(x-X_{0})(\unicode[STIX]{x1D748}_{0}\boldsymbol{n}_{i})\boldsymbol{\cdot }\boldsymbol{j}-(y-Y_{0})(\unicode[STIX]{x1D748}_{0}\boldsymbol{n}_{i})\boldsymbol{\cdot }\boldsymbol{i}\,\text{d}s\right]\text{d}z, & \displaystyle\end{eqnarray}$$

(2.17)$$\begin{eqnarray}\unicode[STIX]{x1D748}_{0}=[-(p_{0}/\unicode[STIX]{x1D716}^{2}+p_{1})\boldsymbol{I}+(\unicode[STIX]{x1D735}\boldsymbol{u}_{0}^{\bot }+\unicode[STIX]{x1D735}\boldsymbol{u}_{0}^{\bot })^{\text{T}}].\end{eqnarray}$$

The leading-order initial conditions for the rod motion are

(2.18a-d)$$\begin{eqnarray}\boldsymbol{X}_{0}(0)=\boldsymbol{{\mathcal{X}}},\quad \unicode[STIX]{x1D6E9}_{0}(0)=\unicode[STIX]{x1D6F7},\quad \dot{\boldsymbol{X}}(0)=\boldsymbol{U},\quad \dot{\unicode[STIX]{x1D6E9}}=\unicode[STIX]{x1D714}.\end{eqnarray}$$

2.3 Leading-order solution

We seek a separable solution for $w_{0}$ of the form

(2.19)$$\begin{eqnarray}w_{0}=f(x,y)\frac{\text{d}p_{0}}{\text{d}z}.\end{eqnarray}$$

Integrating (2.10b) over a cross-section of the fluid domain between the two cylinders, using (2.13) and (2.14), and applying the divergence theorem gives

(2.20)$$\begin{eqnarray}\iint _{\unicode[STIX]{x1D6FA}}w_{0z}\,\text{d}\unicode[STIX]{x1D6FA}=0.\end{eqnarray}$$

Inserting the form (2.19) and using the imposed flux condition (2.12), we conclude that

(2.21)$$\begin{eqnarray}\frac{\text{d}^{2}p_{0}}{\text{d}z^{2}}=0\rightarrow \frac{\text{d}p_{0}}{\text{d}z}=\text{const.}\end{eqnarray}$$

and it follows that

(2.22)$$\begin{eqnarray}w_{0z}=0.\end{eqnarray}$$

Hence (2.10b) becomes

(2.23)$$\begin{eqnarray}\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{u}_{0}^{\bot }=0\end{eqnarray}$$

and the governing equation for the axial fluid flow, (2.10a), is

(2.24)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{\bot }^{2}f(x,y)=1,\end{eqnarray}$$

with no-slip boundary conditions, equations (2.13c) and (2.14b),

(2.25)$$\begin{eqnarray}f=0,\quad \text{on}~\unicode[STIX]{x1D6E4}_{i},\unicode[STIX]{x1D6E4}_{o}.\end{eqnarray}$$

Additionally, using the flux condition (2.3) we can solve for the constant pressure gradient

(2.26)$$\begin{eqnarray}\frac{\text{d}p_{0}}{\text{d}z}=\left[\iint _{\unicode[STIX]{x1D6FA}}f(x,y)\,\text{d}\unicode[STIX]{x1D6FA}\right]^{-1},\end{eqnarray}$$

and thus

(2.27)$$\begin{eqnarray}w_{0}=f(x,y)\left[\iint _{\unicode[STIX]{x1D6FA}}f(x,y)\,\text{d}\unicode[STIX]{x1D6FA}\right]^{-1},\end{eqnarray}$$

combining (2.19) and (2.26).

As $p_{0}$ is constant within a cross-section it has no contribution to the leading-order hydrodynamic force and torque given by (2.16). Additionally, the bracketed terms in (2.16) are independent of $z$, and hence (2.16) and (2.17) to calculate the force and torque reduce to

(2.28a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{0}=\oint _{\unicode[STIX]{x1D6E4}_{i}}\unicode[STIX]{x1D748}_{0}\boldsymbol{n}_{i}\,\text{d}s, & \displaystyle\end{eqnarray}$$

(2.28b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{0}=\oint _{\unicode[STIX]{x1D6E4}_{i}}(x-X_{0})(\unicode[STIX]{x1D748}_{0}\boldsymbol{n}_{i})\boldsymbol{\cdot }\boldsymbol{j}-(y-Y_{0})(\unicode[STIX]{x1D748}_{0}\boldsymbol{n}_{i})\boldsymbol{\cdot }\boldsymbol{i}\,\text{d}s, & \displaystyle\end{eqnarray}$$

(2.29)$$\begin{eqnarray}\unicode[STIX]{x1D748}_{0}=[-p_{1}\boldsymbol{I}+(\unicode[STIX]{x1D735}\boldsymbol{u}_{0}^{\bot }+\unicode[STIX]{x1D735}\boldsymbol{u}_{0}^{\bot })^{\text{T}}].\end{eqnarray}$$

Henceforth, we will drop leading-order subscripts. We will retain the subscript on the only first-order term that appears in our leading-order system, $p_{1}$, to differentiate the first-order pressure that drives flow in the cross-section, equations (2.11), from the leading-order pressure that drives the axial flow, equation (2.10a).

2.4 Model summary and computational approach

Equations (2.24), (2.25) and (2.27) give rise to Poiseuille flow, where the axial flow resistance, ${\mathcal{R}}(\boldsymbol{X},\unicode[STIX]{x1D6E9})$, is equal to the constant pressure gradient required to maintain the axial flow through the annular domain at unit flux. This is given by

(2.30)$$\begin{eqnarray}{\mathcal{R}}(\boldsymbol{X},\unicode[STIX]{x1D6E9})=\left[\iint _{\unicode[STIX]{x1D6FA}}f(x,y)\,\text{d}\unicode[STIX]{x1D6FA}\right]^{-1},\end{eqnarray}$$

a function of domain geometry. We note that for any given flux

(2.31)$$\begin{eqnarray}Q(\boldsymbol{X},\unicode[STIX]{x1D6E9})=\frac{1}{{\mathcal{R}}(\boldsymbol{X},\unicode[STIX]{x1D6E9})}\frac{\text{d}p}{\text{d}z},\end{eqnarray}$$

and so, once ${\mathcal{R}}(\boldsymbol{X},\unicode[STIX]{x1D6E9})$ is determined via (2.30), we can compute $Q(\boldsymbol{X},\unicode[STIX]{x1D6E9})$ for any specified pressure gradient. Motivated by the urological application, in which it is most natural to consider optimising flow rate for a given pressure drop, our approach in §§ 3 and 3.4 will be to fix $\text{d}p/\text{d}z=-1$ as a model input and compute $Q(\boldsymbol{X},\unicode[STIX]{x1D6E9})$ and the associated flow profile as primary outputs. In this view, minimising the resistance is equivalent to maximising the flux.

The equations governing the cross-sectional flow are the Stokes equations (2.11) and incompressibility (2.23), with no-slip boundary conditions (2.13a–b) and (2.14a). As discussed in § 1, the forces and torque on the inner rod, equations (2.28a,b), are linearly related to the imposed velocities via resistance coefficients

(2.32a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}(\boldsymbol{X},\unicode[STIX]{x1D6E9},\dot{\boldsymbol{X}},\dot{\unicode[STIX]{x1D6E9}})=-[\unicode[STIX]{x1D646}(\boldsymbol{X},\unicode[STIX]{x1D6E9})\dot{\boldsymbol{X}}+\unicode[STIX]{x1D63E}(\boldsymbol{X},\unicode[STIX]{x1D6E9})\dot{\unicode[STIX]{x1D6E9}}], & \displaystyle\end{eqnarray}$$

(2.32b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}(\boldsymbol{X},\unicode[STIX]{x1D6E9},\dot{\boldsymbol{X}},\dot{\unicode[STIX]{x1D6E9}})=-[\unicode[STIX]{x1D608}_{zz}(\boldsymbol{X},\unicode[STIX]{x1D6E9})\dot{\unicode[STIX]{x1D6E9}}+\unicode[STIX]{x1D60A}_{x}(\boldsymbol{X},\unicode[STIX]{x1D6E9}){\dot{X}}+C_{y}(\boldsymbol{X},\unicode[STIX]{x1D6E9}){\dot{Y}}], & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D646}(\boldsymbol{X},\unicode[STIX]{x1D6E9})$

$\unicode[STIX]{x1D63E}(\boldsymbol{X},\unicode[STIX]{x1D6E9})$

(2.33a,b)$$\begin{eqnarray}\unicode[STIX]{x1D646}(\boldsymbol{X},\unicode[STIX]{x1D6E9})=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D612}_{xx} & \unicode[STIX]{x1D612}_{xy}\\ \unicode[STIX]{x1D612}_{xy} & \unicode[STIX]{x1D612}_{yy}\end{array}\right],\quad \unicode[STIX]{x1D63E}(\boldsymbol{X},\unicode[STIX]{x1D6E9})=\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D60A}_{x}\\ C_{y}\end{array}\right].\end{eqnarray}$$

The resistance coefficients in (2.33a,b) are also functions of domain geometry. We will compute these in § 4. To explore the effect of geometry on both the axial and cross-sectional resistances we consider configurations (i), (ii) and (iii), as introduced in § 1. In each configuration, we will fix the cross-sectional areas of the cylinders to constrain the space available for fluid. We do this by setting the characteristic radius of the outer cylinder (scaled by $R_{i}$) which we denote $R_{o}$. This reduces the parameter space, allowing for full interrogation of the axial and transverse flows on $\boldsymbol{X}$, $\unicode[STIX]{x1D6E9}$, and the shape of the elliptical cross-sections of the inner and outer cylinders. In § 3.4, we will consider a more complete optimisation problem for axial flow resistance, relaxing our previous restrictions on the considered elliptical geometries.

3.1 (i) Circle in ellipse

In figure 3(a) we consider a circular rod (of dimensionless unit radius) inside a cylinder of elliptical cross-section (with characteristic radius $R_{o}=2$). (Figure 3a has previously appeared in Williams et al. (Reference Williams, Turney, Rauniyar, Harrah, Waters and Moulton2019a), an endourological publication to translate these relevant findings for a clinical audience.) We vary the eccentricity of the outer cylinder’s cross-section while maintaining the cross-sectional area, using (3.2) to determine the corresponding lengths of the major and minor axes. When, for each value of $e_{o}$, the rod is located at the position that maximises flux (see appendix B), we find that the flux initially increases with $e_{o}$. The maximum flux over all $e_{o}$ (data point (ii)) is nearly 50 % higher than the flux for a circular outer cylinder of the same cross-sectional area (data point (i)). The eccentricity at which the maximum flux is achieved is $e_{o}\approx 0.84$ for $R_{o}=2$. We might hypothesise that this coincides with the value for $e_{o}$ where the rod and bounding cylinder match curvature at the cylinder’s vertex, which we denote $\text{e}^{\star }$ (appendix B). However, for $R_{o}=2$, $\text{e}^{\star }\approx 0.78$, which is less than the $e_{o}$ value where maximum flux is achieved. In fact, in configuration (ii) in figure 3(a), although hard to discern from the colour map, the circular rod tangentially touches the elliptic cylinder in two locations. The colour bar in figure 3(a) gives the magnitude of the axial velocity, and demonstrates that the maximum velocity within the cross-section is also larger in configuration (ii), than in either (i) or (iii).

Figure 3.

Dimensionless flow rate, $Q$ with $\boldsymbol{X}=(X,0)$, as a function of (a) $e_{o}$, for $e_{i}=0$, $R_{o}=2$, and (b) $X$, for $e_{o}=0$, $R_{o}=3$. In (a), for each $e_{o}$, the inner circle is located at the position that maximises the axial flux (see appendix B). When $e_{o}=\text{e}^{\star }$, the ellipse matches curvature with the circle at the ellipse vertex. The colour bar provides the magnitude of the axial velocity. (b) Plots $Q$ as a function of $X$ for $e_{i}=0$ (blue) $e_{i}=0.7$ (red) $e_{i}=0.9$ (yellow) oriented at $\unicode[STIX]{x03C0}/2$. The inset plot provides a zoomed-in view from $X=0$ to $X=0.2$. The dashed black line gives the analytical solution by Piercy et al. (Reference Piercy, Hooper and Winny1933) and the black cross the solution by MacDonald (Reference MacDonald1982) (see table 2 for more details).

3.2 (ii) Ellipse in circle

When considering an elliptical inner rod and circular outer cylinder ($e_{o}=0$), we can, due to the rotational symmetry of the outer domain, fix $\boldsymbol{X}=(X,0)$ and vary the inner rod’s position along the $x$-axis only, without loss of generality. Here, we fix $R_{o}=3$, a larger cross-section to allow for a wider range of feasible positions and orientations for the inner rod, and vary $X$, $\unicode[STIX]{x1D6E9}$, and $e_{i}$. We find, as anticipated, that $Q$ increases with $X$ for all values of $e_{i}$ and $\unicode[STIX]{x1D6E9}$. Increasing $e_{i}$ can either increase or decrease $Q$ depending on the position and orientation of the inner rod. For example, when the rod is centred in the outer cylinder – and the flux is independent of $\unicode[STIX]{x1D6E9}$ due to the geometry of the domain – the largest flux is obtained when the rod is circular, and decreases with $e_{i}$. This effect can be seen for $e_{i}=0$, $e_{i}=0.7$ and $e_{i}=0.9$ oriented at $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$ in the inset plot in figure 3(b). In contrast, when the rod is sufficiently offset, the more eccentric rods at $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$ cause less obstruction near the centre of the channel, allowing for higher flow (see figure 3b).

For a circular rod ($e_{i}=0$), we validate our numerical solution against the analytical solution for offset circles (Piercy et al. Reference Piercy, Hooper and Winny1933) (dashed black line in figure 3b) in addition to confirming that our numerical solution approaches the limiting case of touching circles (MacDonald Reference MacDonald1982) as $X$ approaches 2 (black cross in figure 3b). (As Piercy et al. (Reference Piercy, Hooper and Winny1933) presents a solution for the flow rate in terms of an infinite sum, we truncate this sum to determine the flow rate (see § A.3).)

In §§ 3.1 and 3.2, we have considered the particular cases where one of the cylinders has a circular cross-section. As anticipated, in both configurations, the maximum flux was achieved for a rod tangent to the outer boundary. Perhaps unexpectedly, figure 3(a) demonstrated that the optimal configuration may not be one where the inner rod is tangent at a single point. We will explore this in more detail in § 3.4, where we will determine the true optimal geometry that maximises axial flux (for fixed area), relaxing all previous geometric restrictions. However, we will first consider another reduced parameter space, namely elliptical rods in elliptical cylinders of the same eccentricity and orientation; a configuration that allows for clear specification of the inner rod’s position, and a discussion of the effects of eccentricity as a single parameter.

3.3 (iii) Ellipse in ellipse

We now consider cylinders with elliptical cross-sections of the same eccentricities $e=e_{i}=e_{o}$ and orientation ($\unicode[STIX]{x1D6E9}=0$). We fix $R_{o}=2$ and again seek configurations that maximise the flux. (It is worthwhile to note, when comparing flux values between different geometries, that in this section and in § 3.4, $R_{o}=2$, so the space available for fluid flow is the same as in figure 3a (a circular rod in an elliptical cylinder) but less than in figure 3b (an elliptical rod in a circular cylinder) where $R_{o}=3$.)

To describe the position of the inner ellipse relative to the outer ellipse we first define $\unicode[STIX]{x1D703}$ to be the angle between the major axis of the outer ellipse and a line connecting the centres of the two ellipses, see figure 4. The effective radii of the inner and outer ellipses are then defined as

(3.3a,b)$$\begin{eqnarray}r_{i}=a_{i}b_{i}/\sqrt{a_{i}^{2}\sin ^{2}\unicode[STIX]{x1D703}+b_{i}^{2}\cos ^{2}\unicode[STIX]{x1D703}},\quad r_{o}=a_{o}b_{o}/\sqrt{a_{o}^{2}\sin ^{2}\unicode[STIX]{x1D703}+b_{o}^{2}\cos ^{2}\unicode[STIX]{x1D703}}.\end{eqnarray}$$

With these definitions, the maximum distance, $d_{max}$ between the centres at an angle $\unicode[STIX]{x1D703}$ (when the rod and the outer cylinder touch) is given by

(3.4)$$\begin{eqnarray}d_{max}=r_{o}-r_{i}.\end{eqnarray}$$

From this we define an offset parameter, $\unicode[STIX]{x1D719}$, that is the ratio of the distance between the centres of the two ellipses at an angle $\unicode[STIX]{x1D703}$ and $d_{max}$:

(3.5)$$\begin{eqnarray}\unicode[STIX]{x1D719}=d/d_{max}.\end{eqnarray}$$

Thus, the domain is characterised by three parameters: the eccentricity, $e\in [0,1)$; the angle of offset, $\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}/2]$; and the relative offset, $\unicode[STIX]{x1D719}\in [0,1]$ (see figure 4).

Figure 4.

Schematic showing a sample geometry for $e=0.8$, $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/4$ and (a) $\unicode[STIX]{x1D719}=0$, (b) $\unicode[STIX]{x1D719}=0.5$, (c) $\unicode[STIX]{x1D719}=1$. The shaded ellipse is the cross-section of the rod.

Figure 5.

Dimensionless velocity colourmaps with flux values, $Q$, for ellipses of equal eccentricities and orientations. Here $R_{o}=2$. Axes are in $(x,y)$ coordinates and these vary with eccentricity so that the available space for fluid flow is constant. Colour bars reflect different velocity values within the domains. The direction of the axis of minimal resistance, $\boldsymbol{k}_{min}$, is indicated by the dashed white line on each diagram. This will be discussed in § 4. Panels (a–d) show the effect of offset, $\unicode[STIX]{x1D719}$, for $e=0.8$ and $\unicode[STIX]{x1D703}=0$. Offset values are (a) $\unicode[STIX]{x1D719}=0.01$, (b) $\unicode[STIX]{x1D719}=0.35$, (c) $\unicode[STIX]{x1D719}=0.7$ and (d) $\unicode[STIX]{x1D719}=0.99$. Panels (e–h) show the effect of angular position, $\unicode[STIX]{x1D703}$, for $e=0.8$ and $\unicode[STIX]{x1D719}=0.99$. Angular position values are (a) $\unicode[STIX]{x1D703}=0$, (b) $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/6$, (c) $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/3$ and (d) $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$. Panels (i–l) show the effect of eccentricity, $e$, for $\unicode[STIX]{x1D719}=0.99$ and $\unicode[STIX]{x1D703}=0$. Eccentricity values are (a) $e=0$, (b) $e=0.3$, (c) $e=0.6$ and (d) $e=0.9$.

In figure 5(a–d) we fix $e=0.8$ and $\unicode[STIX]{x1D703}=0$, and plot the velocity profiles (with corresponding flux values indicated) for $\unicode[STIX]{x1D719}=0.01$, $\unicode[STIX]{x1D719}=0.35$, $\unicode[STIX]{x1D719}=0.7$ and $\unicode[STIX]{x1D719}=0.99$. We observe that $Q$ increases with $\unicode[STIX]{x1D719}$, a result previously known for $e=0$ (Piercy et al. Reference Piercy, Hooper and Winny1933; Redberger Reference Redberger1962). The effects of $\unicode[STIX]{x1D703}$ and $e$ on $Q$ are less intuitive. In figure 5(e–h) we fix $e=0.8$ and $\unicode[STIX]{x1D719}=0.99$ and vary $\unicode[STIX]{x1D703}$, plotting the associated velocity colour maps for $\unicode[STIX]{x1D703}=0$, $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/6$, $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/3$ and $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$. The flux is largest for $\unicode[STIX]{x1D703}=0$ and smallest for $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/3$ demonstrating, for the chosen parameters, non-monotonicity of $Q$ with $\unicode[STIX]{x1D703}$. We vary $e$ in figure 5(i–l) for fixed $\unicode[STIX]{x1D703}=0$ and $\unicode[STIX]{x1D719}=0.99$. The flux is largest for $e=0.6$ and smallest for $e=0.9$, also demonstrating, for the chosen parameters, non-monotonicity of $Q$ with $e$. For completeness, line plots of $Q$ as a function of independent variation of $e$, $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D703}$ can be found in appendix C.

In figures 6(a)–6(c), we show surface plots of $Q$ in $(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})$ space for three different eccentricities: $e=0$, $e=0.7$ and $e=0.9$, respectively. Note that, as $Q$ varies most significantly with $\unicode[STIX]{x1D719}$, we have restricted to $0.9\leqslant \unicode[STIX]{x1D719}\leqslant 0.99$ to isolate the effects of $\unicode[STIX]{x1D703}$ and $e$. When $e=0$ (figure 6a), $Q$ is independent of $\unicode[STIX]{x1D703}$ due to the rotational symmetry of the domain, and increases monotonically with $\unicode[STIX]{x1D719}$. For eccentric domains, $e=0.7$ and $e=0.9$ (figures 6b and 6c, respectively), $Q$ is minimal at an intermediate value of $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}=0.9$ (the smallest $\unicode[STIX]{x1D719}$ value plotted). The maximum $Q$ is seen in figure 6(b) for $e=0.7$, further validating the existence of a non-zero eccentricity value that maximises flux, $\text{e}^{max}$. This existence of a non-zero eccentricity at which flux is maximised may be of particular interest from an engineering design point of view (Williams et al. Reference Williams, Turney, Rauniyar, Harrah, Waters and Moulton2019a). However, it is noteworthy that the optimal eccentricity itself will change based on the outer geometry. We plot $\text{e}^{max}$ as a function of $R_{o}$ in figure 6(d). As $R_{o}\rightarrow \infty$, $\text{e}^{max}\rightarrow 0$. As $R_{o}$ characterises the ratio between outer and inner cross-sectional areas, the limit as $R_{o}\rightarrow \infty$ corresponds to a negligible obstruction within the channel, and hence, the maximum flux will be attained for a circular cross-section (Williams et al. Reference Williams, Turney, Rauniyar, Harrah, Waters and Moulton2019a).

Figure 6.

Surface plots (a–c) give the dimensionless flow rate, $Q$, as a function of both offset, $\unicode[STIX]{x1D719}$, and angular position $\unicode[STIX]{x1D703}$ of the inner ellipse for three different eccentricity values: (a) $e=0$, (b) $e=0.7$ and (c) $e=0.9$. Here $R_{o}=2$. Panel (d) gives $\text{e}^{max}$ as a function of the characteristic radius $R_{o}$. The $\text{e}^{max}$ value is calculated as the one that produces maximum flow over 100 values from $0.01\leqslant e\leqslant 0.99$ for each radius ratio.

In this section we have determined that for equal eccentricities and orientations of the outer and inner ellipses, there is a position for the inner ellipse and a non-zero eccentricity value that maximises flux. We note that this may not be the configuration that maximises flux for all values of $e_{i}$, $e_{o}$, $X$, $Y$ and $\unicode[STIX]{x1D6E9}$, and we consider this global optimisation problem in more detail in the following section.

3.4 Optimal geometry for axial flow

When solving for axial flow in § 3 we considered three reduced geometric parameter spaces: (i) circle in ellipse, (ii) ellipse in circle and (iii) ellipse in ellipse of the same eccentricity and orientation. We found, in each case, eccentricity values and positions for the inner rod that maximised $Q$. In this section, we will relax all previous assumptions on ellipse eccentricity and explore finding a global optimal geometry for ellipses of fixed cross sectional area, $R_{o}=2$.

For all configurations previously considered, we demonstrated that axial flow resistance was lowered by increasing the distance between the centres of the inner and outer ellipses. Intuitively, we might expect that the optimal configuration is one where the inner ellipse touches the boundary of the outer ellipse and with matching curvature, i.e. when the inner ellipse ‘hugs’ the outer (although this was shown not to be the case for configuration I in § 3.1). This assumption enables a reduced parameter space that can be approached analytically (see appendix D). In figure 7(a) we show a surface plot of $Q$ as a function of $e_{o}$ and $\unicode[STIX]{x1D703}$ produced by this calculation. We see that a global maximum occurs at $\unicode[STIX]{x1D703}=0$ and $e_{o}\approx 0.71$ with corresponding $Q=2.62$ (indicated by the red circle and corresponding flow profile figure 7b). (At this maximum, we also have $e_{i}\approx 0.45$, $\unicode[STIX]{x1D6E9}=0$.)

Figure 7.

A plot of $Q$ as a function of $e_{o}$ and $\unicode[STIX]{x1D703}$ for $R_{o}=2$. Here the inner ellipse touches the outer boundary tangentially at a single point, matching curvature. The red circle indicates the presence of a global maximum.

The optimal configuration in figure 7(a) shows a nearly 40 % improvement in flux over the optimal configuration computed in § 3.3. However, it is still not sufficient to claim this as the global optimum for ellipse geometries. In fact, the results in figure 3(a) suggested that a ‘curve-hugging’ configuration may not be optimal. Motivated by this, we relax assumptions about how the ellipses meet at the boundary and formulate an optimisation problem to traverse the full parameter space. To this end, we define a vector, $\boldsymbol{g}$, containing the five parameters governing the geometry of the domain:

(3.6)$$\begin{eqnarray}\boldsymbol{g}=(e_{i},e_{o},X,Y,\unicode[STIX]{x1D6E9})^{\text{T}}.\end{eqnarray}$$

The boundaries of the inner and outer ellipses, $\unicode[STIX]{x1D6E4}_{i}$ and $\unicode[STIX]{x1D6E4}_{o}$, are

(3.7a,b)$$\begin{eqnarray}\boldsymbol{x}^{\text{T}}\unicode[STIX]{x1D63C}\boldsymbol{x}=0,\quad \boldsymbol{x}^{\text{T}}\unicode[STIX]{x1D63D}\boldsymbol{x}=0,\end{eqnarray}$$

respectively, where $\boldsymbol{x}^{\text{T}}=(x,y,1)$ and matrices $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ are functions of $\boldsymbol{g}$. Thus, we seek the solution to the optimisation problem

(3.8)$$\begin{eqnarray}\min _{\boldsymbol{g}}-Q,\quad \text{subject to}~\boldsymbol{c}\geqslant 0,\end{eqnarray}$$

where $\boldsymbol{c}$ is a vector of six nonlinear constraints, dependent on $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ and hence $\boldsymbol{g}$, that ensure $\unicode[STIX]{x1D6E4}_{i}$ is enclosed by $\unicode[STIX]{x1D6E4}_{o}$ (see appendix E for more details). We solve (2.24) in FEniCS (Alnaes et al. Reference Alnaes, Blechta, Hake, Johansson, Kehlet, Logg, Richardson, Ring, Rognes and Wells2015), and implement MATLAB’s optimisation routine, fmincon, to solve the optimisation problem (3.8) with an interior-point method from eight randomly generated starting points in the feasible five-dimensional parameter space. As the optimisation routine seeks local minima, a multi-start approach allows for global optimisation. The optimal geometry found through this method is shown in figure 7(c) with a corresponding flux, $Q\approx 2.80$. The geometric parameters, see (3.6), are $\boldsymbol{g}\approx (0.68,0.74,1.49,0,\unicode[STIX]{x03C0}/2)^{\text{T}}$. This global optimum gives a non-negligible (nearly 7 %) increase in flux compared with the ‘curve-hugging’ optimum. Indeed, rather than perfectly matching the geometries at the boundary, the global optimal geometry is tangent to the outer boundary in two locations (this was also true in the results in § 3.1). It is interesting to note that this divides the fluid domain into two separate regions, which may be connected to previous results (Ranger Reference Ranger1994, Reference Ranger1996) that demonstrate flux enhancement through creating multiply connected regions. Ranger (Reference Ranger1996) provides a short discussion on the mechanism whereby disconnected regions can lead to increased total flux, discussing the distribution of flux per unit area for flow through a circular pipe. As flux per unit area is low at the edge, removing these regions and replacing them with disconnected areas of higher flux may lead to enhanced flow.

4.1 (i) Circle in ellipse

The analytical solution for Stokes flow between concentric and offset circular cylinders demonstrates a dependence of the forces and torque experienced by the rod on the geometry of the domain (Finn & Cox Reference Finn and Cox2001). We first consider how these results change when the bounding cylinder is elliptical. To limit the candidate geometries, we position the rod along the $x$-axis, i.e. $\boldsymbol{X}=(X,0)$. We fix the cross-sectional area and set $R_{o}=2$.

For a circular rod translating within a circular cylinder, $\unicode[STIX]{x1D612}_{xx}$ and $\unicode[STIX]{x1D612}_{yy}$ increase as the rod approaches the edge of the bounding cylinder, and $\unicode[STIX]{x1D612}_{xy}=0$, due to the symmetry of the domain (Slezkin Reference Slezkin1955). As expected, we find that for $e_{o}>0$, the same behaviour occurs. The effect of $e_{o}$ on the resistance coefficients is not quite so straightforward as the effect of position – as $e_{o}$ increases for fixed $X$, the proximity of the rod to the boundary in the vertical direction decreases while the distance to the boundary in the horizontal direction increases, and at non-equal rates. Due to these competing features, $\unicode[STIX]{x1D612}_{yy}$ exhibits an interesting dependence on $e_{o}$ for $X=0$. In figure 8(a), we plot $\unicode[STIX]{x1D612}_{yy}$ as a function of $e_{o}$ for a centred rod, and we observe a non-monotonic behaviour, highlighting the complexity of boundary interactions. Comparing resistance coefficients for $X\neq 0$ adds yet another complexity: as the range of feasible $X$ values increases with $e_{o}$ (as long as the rod can still fit tangent to the vertex of the outer cylinder), we can either compare values at absolute, $X$, or relative, $X/(a_{o}-1)$, distance from the centre, which qualitatively affects the observed trends. We find that, depending on the metric used, $\unicode[STIX]{x1D612}_{xx}$, $\unicode[STIX]{x1D60A}_{y}$, and $\unicode[STIX]{x1D608}_{zz}$ can also be non-monotonically related to $e_{o}$ when the cylinder is offset (results shown in appendix F). We note that $\unicode[STIX]{x1D60A}_{x}=0$ for a rod centred on the $x$-axis, due to arguments of symmetry.

If, instead of fixing cross-sectional area, we fix the gap between the rod and the cylinder, the domain approaches one of a rod between parallel plates as $e_{o}\rightarrow 1$. An asymptotic approximation to the Stokes force and torque on a cylinder rotating between parallel plates was obtained by Yang et al. (Reference Yang, Wolgemuth and Huber2013), where the ratio of gap width to cylinder radius defines the small parameter, $\unicode[STIX]{x1D716}$. To test whether the resistance coefficients for a rod in an elliptic cylinder approach those for a rotating rod between parallel plates (Yang et al. Reference Yang, Wolgemuth and Huber2013), we consider a rod at $X=0$ in an elliptic cylinder of fixed minor axis, $b_{o}=1.025$, which sets the gaps between the rod and the cylinder at $x=0$, to be 0.025. In figure 8(b) we plot $\unicode[STIX]{x1D608}_{zz}$ as a function of outer eccentricity (line ii), and we see that as $e_{o}\rightarrow 1$ the resistance approaches the value for parallel plates with $\unicode[STIX]{x1D716}=0.025$ (Yang et al. Reference Yang, Wolgemuth and Huber2013) (line i). As we increase $e_{o}$ with fixed minor axis, the cross-sectional area of the elliptic cylinder is also increasing. Thus, as a further geometric comparison, line iii plots $\unicode[STIX]{x1D608}_{zz}$ for a rod centred in a circular cylinder with increasing and equivalent cross-sectional area (i.e. for each $e_{o}$, the bounding circle has cross-sectional area equal to that of an ellipse of eccentricity $e_{o}$ and minor axis $b_{o}=1.025$). Although $\unicode[STIX]{x1D608}_{zz}$ decreases with $e_{o}$ due to the increase in cross-sectional area available for the fluid, $\unicode[STIX]{x1D608}_{zz}$ is larger when the outer boundary is an ellipse compared to a circle of the same cross-sectional area.

Figure 8.

Results for $R_{o}=2$. (a) Rod translating with unit velocity in the $y$-direction, centred at $(0,0)$. $\unicode[STIX]{x1D612}_{yy}$ is plotted as a function of $e_{o}$. Inset schematics are for (left to right) $e_{o}=0$, $e_{o}=0.63$ and $e_{o}=0.9$. (b) Plot of $\unicode[STIX]{x1D608}_{zz}$ for a circular rotating rod centred at $(0,0)$. The rod is rotating: (iv) between parallel plates with gap size $\unicode[STIX]{x1D716}=0.025$ (Yang et al. Reference Yang, Wolgemuth and Huber2013) (v) within an outer elliptic cylinder with $b_{o}=1.025$ ($R_{o}$ not fixed) (vi) within an outer circular cylinder of radius $R_{o}=1.025(1-e_{o}^{2})^{-1/4}$. Solution obtained from (1.3c).

4.2 (ii) Ellipse in circle

The next domain we consider is an ellipse in a bounding circle. As in § 3.2, due to the rotational symmetry of the outer domain, the number of unique combinations of position and orientation of the rod can be reduced; it is sufficient to consider the rod centred at $\boldsymbol{X}=(X,0)$ and to vary $X$, $\unicode[STIX]{x1D6E9}$ and $e_{i}$. In figure 9 we investigate how the resistance coefficients are influenced by the shape of the rod’s cross-section and its orientation and proximity to the outer boundary. In figure 9(a–c), results are plotted as functions of absolute offset, and the orientation of the rod, $\unicode[STIX]{x1D6E9}$, as well as its position, $X$, determine its proximity to the outer cylinder’s boundary, and in a non-trivial way.

Figure 9.

Panels (a) and (b) give $\unicode[STIX]{x1D608}_{zz}$ and $\unicode[STIX]{x1D612}_{xy}$, respectively, as functions of $X$ for $e_{i}=0$, $e_{i}=0.7$ and $e_{i}=0.9$ and $\unicode[STIX]{x1D6E9}=0$, $\unicode[STIX]{x1D6E9}_{0}=\unicode[STIX]{x03C0}/4$ and $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$. In (b), $\unicode[STIX]{x1D612}_{xy}$ is identically zero for $\unicode[STIX]{x1D6E9}=0$ and $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$, so these are not shown as they coincide with the line for $e_{i}=0$. Panel (c) displays $K_{min}$ for $\unicode[STIX]{x1D6E9}=0$ and $e_{i}=0$, $e_{i}=0.7$ and $e_{i}=0.9$. Panel (d) gives $C_{y}$ as a function of $e_{i}$ for $\unicode[STIX]{x1D6E9}=0$ and $X=1.45$. Colour maps of pressure $p_{1}$ near the bounding wall, due to an imposed angular velocity, are shown for $e_{i}=0$ and $e_{i}=0.7$ with the colour bar providing the magnitude of $p_{1}$, centred around $p_{1}=0$. All plots have a dimensionless outer radius of $R_{o}=3$.

We first consider resistance to rotation. We find that $\unicode[STIX]{x1D608}_{zz}$ increases as a function of offset for fixed $e_{i}$ and $\unicode[STIX]{x1D6E9}$. Results are shown in figure 9(a) for three eccentricities, $e_{i}=0$, $e_{i}=0.7$ and $e_{i}=0.9$, oriented at $\unicode[STIX]{x1D6E9}=0$, $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/4$, and $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$. For fixed position and orientation, rotational resistance increases with eccentricity. Without boundary interactions, $\unicode[STIX]{x1D608}_{zz}$ would be independent of $\unicode[STIX]{x1D6E9}$; however, due to the presence of the bounding circular cylinder, for moderate $X\neq 0$, $\unicode[STIX]{x1D608}_{zz}$ is largest when $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$ and smallest when $\unicode[STIX]{x1D6E9}=0$. As evidenced by the crossing lines for $e_{i}=0.9$ in figure 9(a), when $X$ is sufficiently large, a rod oriented at $\unicode[STIX]{x1D6E9}=0$ receives the largest opposing torque, which is likely due to the proximity of the vertex of the major axis to the boundary. We note that, if we increase the $X$ range in figure 9(a), this effect is also seen for $e_{i}=0.7$. (Results shown in appendix F.)

Turning now to translational resistance, our results validate the intuition that when the major or minor axis of the rod’s cross-section is aligned both with the direction of motion and a radial line of the outer cylinder’s cross-section, the force on the rod acts parallel to its direction of motion, and the lift force is zero. Brenner (Reference Brenner1962a) proved this to be true for the motion of a bounded arbitrary particle.

Of the three translational resistances, we first consider the lift force, $\unicode[STIX]{x1D612}_{xy}$. We find, that for fixed $e_{i}$ and $X$ away from the boundary, $\unicode[STIX]{x1D612}_{xy}$ is maximised for $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/4$. Thus the outer cylinder does not seem to change classic results for an elliptic cylinder in an unbounded fluid (Lee & Leal Reference Lee and Leal1986). To see how the lift force changes with inner eccentricity, in figure 9(b) we fix the orientation at the maximal value $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/4$ for eccentricities $e_{i}=0$, $e_{i}=0.7$ and $e_{i}=0.9$, and plot $\unicode[STIX]{x1D612}_{xy}$ as a function of $X$. We find a significant increase in lift force with eccentricity, and also a non-monotonic relation with position in the case of $e_{i}=0.9$.

Next, we analyse $\unicode[STIX]{x1D612}_{xx}$, $\unicode[STIX]{x1D612}_{yy}$ and $K_{min}$. For relatively small offset, $\unicode[STIX]{x1D612}_{xx}$ is smallest (and $\unicode[STIX]{x1D612}_{yy}$ largest) when $\unicode[STIX]{x1D6E9}=0$. We restrict ourselves to this orientation, and plot $\unicode[STIX]{x1D612}_{xx}$ (dotted line) and $\unicode[STIX]{x1D612}_{yy}$ (dashed line) as functions of $X$ for $e_{i}=0.7$ and $e_{i}=0.9$ in figure 9(c). An interesting feature is the non-monotonic behaviour of $\unicode[STIX]{x1D612}_{yy}$ with $X$. The solid coloured lines give the value of the minimum resistance, $K_{min}$. Looking particularly at $e_{i}=0.9$ we see that for small values of $X$, $K_{min}$ coincides with $\unicode[STIX]{x1D612}_{xx}$ (and thus the eigenvector $\boldsymbol{k}_{min}=\boldsymbol{i}$). However, near $X=1.2$ the lines $\unicode[STIX]{x1D612}_{xx}$ and $\unicode[STIX]{x1D612}_{yy}$ cross, and the direction of minimum resistance switches to be aligned with the $y$-axis. The same behaviour is seen for $e_{i}=0.7$, with the switching occurring at a smaller offset (near $X=1$). The solid lines denoting $K_{min}$ for each eccentricity also cross. When $X$ is small, $K_{min}$ is largest for $e_{i}=0$ and smallest for $e_{i}=0.9$ and vice versa for large $X$. Again, this can be explained by the significance of boundary interactions: far from the boundary the more eccentric ellipse has a direction of small resistance – its major axis – but closer to the boundary the major vertex is in closer proximity to the boundary, compared to ellipses of smaller eccentricity, hence the resistance in both the $x$- and $y$-directions rapidly increases. We have also demonstrated that for $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$, $\boldsymbol{k}_{min}$ is aligned with the $y$-axis for $e_{i}=0$, $e_{i}=0.7$ and $e_{i}=0.9$ (at least for $0\leqslant X\leqslant 1.45$). (Results shown in appendix F.)

Finally, we consider the coupling between translation and rotation. We find, for fixed $e_{i}$ and $X$ away from the boundary, $\unicode[STIX]{x1D60A}_{x}$ is largest for $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/4$ and $C_{y}$ is largest for $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$. (Results shown in appendix F.) When the rod is oriented at $\unicode[STIX]{x1D6E9}=0$, there is coupling between rotation and translation in the $y$-direction ($C_{y}\neq 0$). This is known analytically for $e_{i}=0$ (Slezkin Reference Slezkin1955), and in this classic case, $C_{y}$ increases monotonically with $X$. Interestingly, for more eccentric ellipses, this is no longer true; in fact, $C_{y}$ can change sign as a function of position along the $X$-axis. Alternatively, this behaviour can be illustrated by plotting $C_{y}$ as a function of eccentricity for fixed $X$ sufficiently close to the outer boundary. This is shown in figure 9 for $X=1.45$ (the furthest-right position in figure 9a–c). We observe that $C_{y}$ decreases as a function of $e_{i}$, changing sign from positive to negative. To explain this behaviour, consider a rotating rod with imposed angular velocity, in which case $C_{y}$ characterises the induced translational force in the $y$-direction. Example pressure profiles for $e_{i}=0$ and $e_{i}=0.775$ near the proximal boundary are shown in the insets, and these show that the pressure drop in the lubrication layer between the rod and the outer boundary changes direction as a function of $e_{i}$. Due to the proximity to the boundary, we can appeal to Reynolds’ lubrication equation to understand this inversion in terms of the behaviour of the fluid in the lubrication layer: for a circular surface rotating near a fixed plane, the pressure drop is in the opposite direction compared to the pressure drop for a highly curved surface moving parallel to the plane (Wannier Reference Wannier1950). For our system, rotating an ellipse of high eccentricity is analogous to moving the vertex (and point closest to the boundary) parallel to the boundary, whereas rotating an ellipse of small eccentricity does not change significantly the minimal distance to the boundary and hence behaves more like a rotating circle.

4.3 (iii) Ellipse in ellipse

In this section we revisit the configuration considered in § 3.3, with both the rod and bounding cylinder having elliptical cross-sections of the same eccentricity, $e$, and orientation ($\unicode[STIX]{x1D6E9}=0$). We set $R_{o}=2$, and parameterise $\boldsymbol{X}$ by offset $\unicode[STIX]{x1D719}$ and angular position $\unicode[STIX]{x1D703}$, illustrated in figure 4. Returning to figure 5, the direction of minimum resistance is marked with a dashed white line. We see that when the rod is centred in the outer cylinder, $\boldsymbol{k}_{min}$ lies along the major axis of the ellipse. As we move towards the wall, this direction switches to the minor axis, due to the build-up of pressure between the rod and the proximal edge of the outer cylinder. As we break symmetry and vary $\unicode[STIX]{x1D703}$, $\boldsymbol{k}_{min}$ no longer lies along an axis of symmetry of the inner ellipse and for $\unicode[STIX]{x1D719}=0.99$ it varies from vertical for $\unicode[STIX]{x1D703}=0$ to horizontal for $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$, see figure 5(e–g).

Figure 10.

Surface plots (a–c) give the value of $K_{min}$ as a function of both offset, $\unicode[STIX]{x1D719}$, and angular position $\unicode[STIX]{x1D703}$ (where these values are defined as in § 3) of the inner ellipse for three different eccentricity values: (a) $e=0$, (b) $e=0.7$ and (c) $e=0.9$. Panels (d) and (e), respectively, show the maximum and minimum value of $K_{min}$ gives the minimum value of $K_{min}$ for a given equal inner and outer ellipse eccentricity value. The geometries that provide the $\min (K_{min})$, (d), and $\max (K_{min})$, (e), for $e=0$, $e=0.4$, $e=0.7$ and $e=0.9$ are demonstrated by diagrams and the white-dashed line on each of these shows the direction of minimum resistance. $R_{o}=2$.

Figure 6(a–c) shows the axial flow rate in the $(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})$ parameter space for three different eccentricity values, and in figure 10(a–c), we plot $K_{min}$ for the same domains. For both $e=0.7$ and $e=0.9$, the smallest value of $K_{min}$ occurs at $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ and an intermediate $\unicode[STIX]{x1D719}$ value $0<\unicode[STIX]{x1D719}<1$. Here, $K_{min}$ is in the $x$-direction, aligned with the major axis of the rod. This is similar to the behaviour determined by Dvinsky & Popel (Reference Dvinsky and Popel1987) of a circular cylinder translating between two plates, parallel to the boundaries, where the lowest force was incurred at a position between the centreline and the wall, with the exact location being dependent upon the ratio of the radius of the cylinder to the distance between the walls.

In figure 10(d), the minimum value of $K_{min}$ over $(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})$-space is plotted as a function of $e$. Several points have been illustrated to show the domain that produces the value, and the $\boldsymbol{k}_{min}$ axis is demonstrated on these by a dashed white line. We see that $K_{min}$ has a local minimum near $e=0.7$. Note also that the position where the minimum $K_{min}$ is achieved changes with $e$, moving from the centre of the outer cylinder to the edge with increasing eccentricity. Figure 10(e) shows the maximum value of $K_{min}$ for increasing eccentricity, and again, the locations along with the $\boldsymbol{k}_{min}$ axis that produces these values is demonstrated at four points. Here there is a monotonic increase with eccentricity, and the location of the maximum $K_{min}$ always occurs on the major axis of the outer cylinder, closest to the wall ($\unicode[STIX]{x1D703}=0$, $\unicode[STIX]{x1D719}=0.99$).

Connecting these results with those in § 3.3, we find that for each $e$, the location where $K_{min}$ is maximised coincides with the position where $Q$ is maximised. This can be seen by comparing figure 10(a–c) with figure 6(a–c). Contour plots of $Q$, $\log (K_{min})$ and $\log (K_{max})$ for $e_{o}=e_{i}=0.7$, $\unicode[STIX]{x1D6E9}=0$ are shown in figure 11 as functions of $(X,Y)$. The plots are computed for $0\leqslant \unicode[STIX]{x1D719}\leqslant 0.99$, $0\leqslant \unicode[STIX]{x1D703}\leqslant \unicode[STIX]{x03C0}/2$, and the black boundary indicates the edge of the outer ellipse. The coincident positions of $\max (Q)$ and $\max (K_{min})$ are indicated in figures 11(a) and 11(b), respectively, by the red circles.

We note that, by parametrising the position of the rod by $\unicode[STIX]{x1D719}$, the minimum gap between the rod and the bounding cylinder depends on the angular position $\unicode[STIX]{x1D703}$. To test the sensitivity – of positions of highest minimum resistance and maximum flux – to the minimum fixed distance between the rod and the cylinder wall, we consider, in appendix G, an absolute minimum distance between the rod and bounding cylinder, rather than one that depends on $\unicode[STIX]{x1D703}$, as well as different maximum $\unicode[STIX]{x1D719}$ values. We find, for all considered minimum distances between the rod and cylinder, that the correspondence between rod positions of highest minimum resistance and maximum flux persists.

The resistance contours, figures 11(b) and 11(c), show increasingly narrow contours with proximity to the edge of the outer ellipse, demonstrating a rapid increase in the gradient of $K_{min}$ and $K_{max}$, respectively. In contrast, the distances between neighbouring $Q$ contour lines in figure 11(a) remain relatively constant. Thus, it is interesting to note that small changes to the position of the inner rod when it is near the outer boundary may incur vast changes to the resistance to cross-sectional motion, but only a small change to the corresponding axial flux. The arrows in figures 11(b) and 11(c) demonstrate the directions of minimum and maximum resistance, respectively. Along the major axis of the outer ellipse, and close to its centre, the direction of minimum resistance is towards the vertex of the inner rod, demonstrating the effect of the rod’s streamlined elliptical cross-section. However, the directions of the arrows change once the rod is sufficiently far from the centre, and the competing effects of the bounding walls become significant. This is in agreement with the results in figure 9(c). As anticipated, when the rod is near the outer boundary, the direction of minimum resistance is tangent to the wall and the direction of maximum resistance perpendicular to it. The arrows in figures 11(b) and 11(c) allow some conjecturing as to local trajectories of the inner rod; however, it is important to recall that the results in figure 11 are computed for fixed orientation of the inner rod ($\unicode[STIX]{x1D6E9}=0$). In actuality, torque will be induced by the translational motion of the rod, thus changing the resistance contours. Although similar trends exist between $Q$, $K_{min}$ and $K_{max}$, as shown in figure 11, the non-monotonic contour shapes in figure 11(b), as well as the difference in gradients, indicate that the apparent correspondence between axial and cross-sectional resistance is complicated, even in the particular reduced parameter space (of equal ellipse eccentricity and orientation) explored in figures 10 and 11.

Figure 11.

Panels (a), (b) and (c), show $Q$, $\log (K_{min})$ and $\log (K_{max})$ as functions of position for ellipses of equal eccentricity $e_{o}=e_{i}=0.7$, and equal orientation ($\unicode[STIX]{x1D6E9}=0$). The arrows in  (b) and (c) show the directions of $\boldsymbol{k}_{min}$ and $\boldsymbol{k}_{max}$, respectively. Red circles in  (a) and (b) indicate the coincident locations of maximum $Q$ and maximum $K_{min}$.