Laminar boundary layer separation and reattachment on a rotating sphere

Abstract

A new model of the steady boundary layer flow around a rotating sphere is developed that includes the widely observed collision and subsequent eruption of boundary layers at the equator. This is derived following the Segalini & Garrett (J. Fluid Mech., vol. 818, 2017, pp. 288–318) asymptotic approach for large Reynolds numbers but replacing the Smith & Duck (Q. J. Mech. Appl. Maths, vol. 30, issue 2, 1977, pp. 143–156) correction with a higher-order version of the Stewartson (Grenzschichtforschung/Boundary Layer Research, 1958, pp. 59–71. Springer) model of the equatorial flow. The Stewartson model is then numerically solved, for the first time, via a geometric multigrid method that solves the steady planar Navier–Stokes equations in streamfunction-vorticity form on large rectangular domains in a quick and efficient manner. The results are then compared with a direct numerical simulation of the full unsteady problem using the Semtex software package where it is found that there is broad qualitative agreement, namely the separation and reattachment of the boundary layer at the equator. However, the presence of unobserved behaviour such as a large area of reverse flow seen at lower Reynolds numbers than those observed in other studies, and that the absolute error increases with Reynolds number suggest the model needs improvement to better capture the physical dynamics.

1. Introduction

The flow around a rotating sphere in a quiescent fluid provides a useful paradigm for the study of fundamental fluid dynamical problems, such as boundary layer collisions and separations, and in other scientific fields including astrophysics. Experimental studies by Bowden & Lord (1963), Hollerbach et al. (2002), Calabretto et al. (2015) and Calabretto, Denier & Levy (2019) all observe the same flow characteristics, specifically the presence of inflow at the polar regions causing boundary layers to form on the sphere surface that are convected towards the equator. Due to the equatorial symmetry, these are formed on both hemispheres resulting in a boundary layer collision at the equator that then evolves into a radial jet creating a toroidal vortex. This has also been observed in recent numerical studies such as Calabretto et al. (2015, 2019) and Calabretto & Denier (2019), however, these do not provide insight into the underlying physics.

The flow around a rotating sphere was first studied by Stokes (1845) who described the behaviour of a slowly rotating sphere for small Reynolds numbers $\textit {Re}=a^2\varOmega /\nu$, where $a$ and $\varOmega$ are the radius and angular velocity of the sphere, respectively, and $\nu$ is the kinematic viscosity. Further theoretical results for small $\textit {Re}$ were obtained by Lamb (1924), Bickley (1938), Collins (1955), Thomas & Walters (1964) and Takagi (1977); whilst for $\textit {Re}\sim O(100)$, Dennis, Singh & Ingham (1980) calculated series solutions of Gegenbauer functions, but these become more difficult to obtain as $\textit {Re}$ increases due to the nonlinearity of terms in the series. For large $\textit {Re}$, boundary layer theory provides a suitable model of the flow near the surface where the sphere imparts the fluid with angular momentum (Segalini & Garrett 2017). The governing equations that describe the boundary layer flow were initially derived by Howarth (1951) who proposed two solution methods: a solution based on a Pohlhausen technique and a series solution based on the latitudinal angle where the latter, at leading order, recovered the von Kármán equations for the rotating disk flow, emphasising the strong connection between the flow at the poles with that of the rotating disk. Although Howarth (1951) derived the series solution, it was first solved and utilised by Banks (1965) to obtain solutions for the boundary layer. The series solution approximates the full numerical solution well at small latitudes (Manohar 1967; Banks 1976), however, as the equator is approached there is divergent behaviour, meaning a full numerical solution of the equations is mandated (Garrett & Peake 2002). Furthermore, the parabolic structure of the boundary layer equations implies that there cannot be any latitudinal stagnation points at the equator, suggesting a collision of two boundary layers, formed at both poles, is unavoidable, which must be described by a new elliptic structure (Simpson & Stewartson 1982).

The boundary layer equations of Howarth (1951) model the flow well near the sphere surface until the equator is approached where the boundary layers collide and erupt into a radial jet. In order to model this area, Stewartson (1958) proposed that this region should be a thin viscous structure described by the planar Navier–Stokes equations with an inviscid outer flow at the equator-sphere junction where he further hypothesised the existence of a small zone of recirculation. In contrast, Smith & Duck (1977), based on an analysis of a dual layer structure of overall size $O(\textit {Re}^{-3/7})\times O(\textit {Re}^{-1/2})$, conjectured a more extensive recirculation region of $O(\textit {Re}^{-3/14})$. This is opposed to the numerical study of Dennis, Ingham & Singh (1981) who found that this interaction zone is of $O(\textit {Re}^{-1/4})$; however, no recirculation has been observed in any prior experiment or numerical simulation (Segalini & Garrett 2017), until a recent numerical study by Calabretto et al. (2019) observed a small area of reverse flow at significantly large Reynolds numbers, $\textit {Re}\geq 8\times 10^4$. Despite this, no large structures of the kind Smith & Duck (1977) proposed were seen suggesting that the Stewartson (1958) model of the equatorial flow may be more qualitatively correct. Thus, around the equatorial region, the physical mechanisms and behaviour remain ambiguous. Furthermore, work concerning the stability of the flow around a rotating sphere by Segalini & Garrett (2017) features a model of the steady flow that incorporates the Smith & Duck (1977) model via a correction to the velocity profiles as the equator is approached; and by incorporating $1/\textit {Re}$ corrections with an inviscid outer flow, the resultant stability calculations agreed remarkably with experiments, despite not including the post-collisional flow. Hence, the Smith & Duck (1977) model could be negligible and have no meaningful effect on the stability of the flow above the equator, and it is therefore necessary to replace this part of the analysis with the model developed by Stewartson (1958). He formulated the equations of motion for this area, which seem to be equivalent to a streamfunction-vorticity formulation, but did not solve them; it is believed that this is due to the presence of the small zone of recirculation where the vorticity is unknown (Segalini & Garrett 2017; Calabretto et al. 2019).

In summary, the flow around the rotating sphere is well studied, but there still remains some significant uncertainty. As the boundary layer approaches the equator, the parabolic structure of the governing equations break down and do not accurately model the separation and subsequent reattachment of the boundary layer along the equator. There are currently two competing models that describe this behaviour: Stewartson (1958) and Smith & Duck (1977). Both assume the flow can be described by the planar Navier–Stokes equations, however, the higher-order analysis of Smith & Duck (1977) suggests behaviour not seen by any experiment or numerical study. The presence of a small recirculation zone at large $\textit {Re}$ does suggest that the model of Stewartson (1958) may be qualitatively correct; however, the equations have not been solved (Calabretto et al. 2019). This is believed to be due to the unknown form of the vorticity, although this can be overcome by coupling in the vorticity that results in an additional equation that can be solved. Unfortunately, this requires solving nonlinear Partial Differential Equations (PDEs) on large domains, which is not trivial. Solving these PDEs is the first aim of this work and once solutions are obtained, then this will allow further understanding of the physical mechanisms of the flow, specifically, the behaviour around the equator.

A brief summary of the boundary layer equations and Stewartson (1958) model of the equatorial flow is given in § 2 that is analogous to the planar Navier–Stokes equations. In order to solve these, a numerical method utilising the geometric multigrid method is outlined in § 3 and the results can be seen in § 4, which are then discussed in § 5.

2. Boundary layer analysis

Consider at the origin of a fixed reference frame, a rotating sphere of radius $a$ with angular velocity $\varOmega$ immersed in an otherwise quiescent fluid. Furthermore, let the problem be described by a spherical coordinate system, where $(r,\theta,\phi )$ denotes the radial, latitudinal and azimuthal coordinates, respectively. Let $(W,U,V)$ denote the velocities in the radial, latitudinal and azimuthal directions, respectively, and let $P$ be the pressure. Assuming that the flow is steady and axisymmetric, then $(U,V,W,P)$ is independent of the time $t$ and the azimuthal angle $\phi$, i.e. $\partial _t\,\cdot\, =\partial _\phi \,\cdot\, =0$, where $\partial _t = \partial /\partial t$ and $\partial _\phi = \partial /\partial \phi$. Non-dimensionalising all physical quantities by $a$, $\varOmega$ and the fluid density $\rho$, the steady, incompressible, axisymmetric Navier–Stokes problem in spherical coordinates becomes

(2.1a)\begin{gather} W\partial_r W+ \frac{U}{r}\partial_\theta W - \frac{U^2 + V^2}{r} ={-}\partial_r P + \frac{1}{\textit{Re}}\left(\nabla^2 W-\frac{2}{r^2}\partial_\theta U - \frac{2W}{r^2}-\frac{2U}{r^2}\cot{\theta}\right), \end{gather}

(2.1b)\begin{gather}W\partial_r U+ \frac{U}{r}\partial_\theta U + \frac{UW}{r} - \frac{V^2}{r}\cot{\theta} ={-}\frac{1}{r}\partial_\theta P + \frac{1}{\textit{Re}}\left(\nabla^2U+\frac{2}{r^2}\partial_\theta W - \frac{U}{r^2 \sin^2{\theta}}\right), \end{gather}

(2.1c)\begin{gather}W\partial_r V+ \frac{U}{r}\partial_\theta V + \frac{WV}{r} + \frac{UV}{r}\cot{\theta} = \frac{1}{\textit{Re}}\left(\nabla^2V - \frac{V}{r^2 \sin^2{\theta}}\right), \end{gather}

(2.1d)\begin{gather}\partial_r W + \frac{2W}{r} + \frac{1}{r}\partial_\theta U + \frac{U}{r}\cot{\theta} = 0, \end{gather}

where the derivatives are expressed as before (e.g. $\partial _r = \partial /\partial r$) and

(2.2)\begin{equation} \nabla^2\,\cdot\,= \frac{1}{r^2}\partial_r(r^2\,\cdot\,) + \frac{1}{r^2\sin^2\theta}\partial_\theta(\sin\theta\partial_\theta\,\cdot\,). \end{equation}

The boundary conditions of the system are

(2.3a)\begin{gather} U = V - \sin{\theta} = W = 0\quad \text{at } r=1, \end{gather}

(2.3b)\begin{gather}U = V = W = 0 \quad\text{as } r \to \infty. \end{gather}

The first set of conditions (2.3a) refer to the no-slip condition on the sphere surface and the second set (2.3b) specifies that the fluid is at rest far away from the sphere.

Equation (2.1) is highly coupled and nonlinear, and in order to solve it, a Direct Numerical Simulation (DNS) is needed that can require huge amounts of computational resources. To reduce this, boundary layer approximations will be made as it is expected that the bulk behaviour of the flow is located close to the sphere surface (Segalini & Garrett 2017). Furthermore, by exploiting the spherical symmetry of the geometry, the domain is truncated to one quadrant that then necessitates the homogenous Neumann boundary conditions

(2.3c)\begin{equation} \partial_\theta U = \partial_\theta V = \partial_\theta W = 0\quad \text{at}\ \theta=0 \text{ and } \frac{\rm \pi}{2}, \end{equation}

which represents symmetry conditions along the pole and equator.

The equations that govern the boundary layer above the equator are discussed in many other works including Howarth (1951), Banks (1965) and, more recently, Segalini & Garrett (2017), and so are only briefly discussed here. For convenience, the boundary layer equations are given by

(2.4a)\begin{gather} \partial_\eta \bar{W} + \partial_\theta U + U\cot{\theta} = 0, \end{gather}

(2.4b)\begin{gather}\bar{W}\partial_\eta U + U\partial_\theta U - V^2\cot{\theta} = \partial_\eta^2U, \end{gather}

(2.4c)\begin{gather}\bar{W}\partial_\eta V + U\partial_\theta V +UV\cot{\theta} = \partial_\eta^2V, \end{gather}

with

(2.5)\begin{equation} \bar{P}(\eta,\theta) ={-}\int_{\eta}^{\infty}U^2(\xi,\theta)+V^2(\xi,\theta) \,{\rm d}\xi, \end{equation}

where $\bar {W}=W/\epsilon$, $\bar {P}=P/\epsilon$,

(2.6)\begin{equation} \eta=\frac{r-1}{\epsilon} \end{equation}

is the stretched radial variable that localises the flow close to the sphere surface, and $\epsilon =\delta /a\ll 1$ is a small perturbation parameter that can be interpreted as the non-dimensional boundary layer thickness where the quantity $\delta =a/\sqrt {\textit {Re}}$ is the characteristic boundary layer thickness, i.e. $\epsilon =1/\sqrt {\textit {Re}}$. The boundary conditions are given by

(2.7a)\begin{gather} U = V - \sin{\theta} = \bar{W} = 0\quad \text{at } \eta=0, \end{gather}

(2.7b)\begin{gather}U = V = 0 \quad\text{as } \eta \to \infty, \end{gather}

(2.7c)\begin{gather}\partial_\theta U = \partial_\theta V = \partial_\theta W = 0\quad \text{at } \theta=0. \end{gather}

As discussed in other works, there is radial inflow at the poles, due to the similarity of the geometry with the rotating disk flow, as the centrifugal effect forces the fluid outwards along the sphere. The parabolic structure of (2.4) then drives the fluid towards the equator where it collides with the boundary layer emanating from the other hemisphere. However, once the equator is approached, the solutions of (2.4) cannot accurately describe the flow effectively because information about the solution is required from both upstream and along the equator, necessitating an elliptic structure of the equations.

Recent numerical studies suggest that the elliptic model of Stewartson (1958) describes the flow qualitatively well (Calabretto et al. 2019), hence, following Stewartson (1958), introduce the scaled coordinate

(2.8)\begin{equation} \beta=\frac{\theta-{\rm \pi}/2}{\epsilon}, \end{equation}

which localises the flow to the equator. The flow around the equator can be found by considering that $(U,V,W)\sim O(1)$, in order to facilitate the transfer of momentum from the sphere to along the equator, and substituting $P=\epsilon \bar {P}$ (as it is expected that $P\sim O(\epsilon )$ following Segalini & Garrett 2017), (2.8) and (2.6) into (2.1). Subsequently, by only considering leading-order terms, akin to taking the limit $\textit {Re}\rightarrow \infty$, and using the expansions $\cot \theta \approx -\epsilon \beta +O(\epsilon ^3)$ and $1/\sin ^2\theta \approx 1+O(\epsilon ^2)$, then (2.1) reduces to

(2.9a)\begin{gather} \partial_\eta W + \partial_\beta U = 0, \end{gather}

(2.9b)\begin{gather}W\partial_\eta U + U\partial_\beta U ={-}\epsilon\partial_\beta \bar{P} + \epsilon\left(\partial_\eta^2U + \partial_\beta^2U\right), \end{gather}

(2.9c)\begin{gather}W\partial_\eta V + U\partial_\beta V = \epsilon\left(\partial_\eta^2V + \partial_\beta^2V\right), \end{gather}

(2.9d)\begin{gather}W\partial_\eta W + U\partial_\beta W ={-}\epsilon\partial_\eta \bar{P} + \epsilon\left(\partial_\eta^2 W + \partial_\beta^2 W\right). \end{gather}

Viscous terms with $\epsilon$ have been retained as it is expected that they generate internal boundary layers along the sphere surface and equatorial plane of size $O(\epsilon \sqrt {\epsilon })$. At first glance, this may seem unintuitive, however, by neglecting these viscous terms Stewartson (1958) found that $U(\eta =0,\beta )\propto f(\beta )$, where $f(\beta )$ is any function such that $f(\beta )=0$ for $\beta =0$ and as $\beta \rightarrow -\infty$, which can only resolve the condition $U(0,\beta )=0$ if it is neutralised by an internal boundary layer. The pressure gradients have also been retained in order to calculate the pressure and to keep the number of unknowns and equations the same, although, they can be easily discarded as seen later. Note that Stewartson (1958) discarded these terms but they have been retained for the reasons outlined prior; hence, it is considered that (2.9) is a higher-order model. The boundary conditions are

(2.10a)\begin{gather} U = V - 1 = W = 0 \quad\text{at } \eta=0, \end{gather}

(2.10b)\begin{gather}U = \partial_{\eta} V = \partial_{\eta} W = 0 \quad\text{as } \eta \to \infty, \end{gather}

(2.10c)\begin{gather}U = \partial_{\beta} V = \partial_{\beta} W = 0 \quad\text{at } \beta=0, \end{gather}

(2.10d)\begin{gather}U - U_{BL}(\eta) = V - V_{BL}(\eta) = W = 0 \quad\text{as } \beta \to -\infty, \end{gather}

where $U_{BL}$ and $V_{BL}$ are the inflow velocities from the boundary layer solution and $\sin \theta \approx 1+O(\epsilon ^2)$. The condition (2.10a) is the no-slip condition on the sphere surface, (2.10d) refers to the inlet whilst (2.10b) is a Neumann condition for the outflow, as the radial jet will be established, and (2.10c) is a symmetry condition due to the similarity of the two hemispheres where $U=0$ is imposed to facilitate the pure radial planar flow of the boundary layer. These are the same set of equations proposed in Segalini & Garrett (2017) but with the introduction of the Neumann conditions for the equator, following Dennis et al. (1981) who found that the radial velocity does not vanish along the equator. The reason for introducing (2.10b) is because the flow should only vary noticeably at the $r\sim O(1)$ scale.

First, it is useful to note that the azimuthal component (2.9c) is uncoupled from (2.9) reducing it to a linear system of one less dimension. Furthermore, as $\eta$ and $\beta$ localise the flow to a small region at the equator-sphere junction, it can be assumed that the sphere curvature is locally flat. These two observations effectively reduce (2.9) to solving the planar incompressible Navier–Stokes equations, the difference being that ${\textit {Re}^{-1}\rightarrow \epsilon = 1/\sqrt {\textit {Re}}}$, and a linear advection–diffusion equation (2.9c) for the azimuthal velocity $V$. A streamfunction $\psi$ can be introduced via the relations

(2.11a,b)\begin{equation} W=\partial_\beta \psi \quad\text{and}\quad U={-}\partial_\eta \psi, \end{equation}

eliminating the continuity equation (2.9a). By introducing the local vorticity,

(2.12)\begin{equation} \omega=\partial_\eta U - \partial_\beta W, \end{equation}

(2.9) can be rewritten in streamfunction-vorticity form by calculating $\partial _\eta$[(2.9d)]$- \partial _\beta$[(2.9b)] to obtain

(2.13a)\begin{gather} -\omega = \Delta\psi, \end{gather}

(2.13b)\begin{gather}(\partial_\beta\psi)\partial_\eta\omega - (\partial_\eta\psi)\partial_\beta \omega = \epsilon\Delta\omega, \end{gather}

(2.13c)\begin{gather}(\partial_\beta\psi)\partial_\eta V - (\partial_\eta\psi)\partial_\beta V = \epsilon\Delta V, \end{gather}

with the pressure determined by solving

(2.14)\begin{equation} \epsilon\Delta\bar{P} = 2\left((\partial_\eta W)\partial_\beta U-(\partial_\beta W)\partial_\eta U\right), \end{equation}

where $\varDelta \,\cdot\, =\partial _\eta ^2\,\cdot\, +\partial _\beta ^2\,\cdot$ is the local Laplacian operator. The boundary conditions for the pressure are

(2.15) \begin{align} \bar{P}(\eta,\beta\rightarrow-\infty) \!= \bar{P}_{BL}(\eta),\quad \partial_\beta\bar{P} \!= 0 \quad\text{at } \beta\!=0,\quad \partial_\eta\bar{P} \!= 0 \quad \text{at } \eta\!=0 \text{ and as } \eta\rightarrow\infty, \end{align}

where $P_{BL}(\eta )$ is the boundary layer solution (2.5). The first condition refers to a matching condition from the upstream boundary layer, the second refers to the symmetry condition along the equator, the third is a Neumann condition for solid walls and the final condition represents that the solution should be constant at this scale, i.e. $O(\epsilon )$.

A similar form to (2.13) was discussed by Stewartson (1958), namely the separation of the azimuthal component and the vorticity equation (2.13a), but not the inclusion of the higher-order viscous terms. Consequently, Stewartson's model represents an inviscid model of the equatorial flow, while broadly true for the flow outside the boundary layer as investigated by Segalini & Garrett (2017), the necessity of introducing an internal boundary layer and (2.13b) allows a solution to be obtained. This explains why no solutions, analytical or numerical, have been procured as the vorticity, $\omega$, was unknown. Furthermore, as (2.9) has been transformed to (2.13), the boundary conditions (2.10) also need to be transformed, which can be seen in Appendix A.

3. Numerical methods

First, it is pertinent to note that the boundary data $U_{BL}$, $V_{BL}$ and $\bar {P}_{BL}$ for the inlet condition (2.10d) is obtained by solving (2.4) subject to (2.7) for the boundary layer upstream of the equator via the Newton space-marching method outlined in Segalini & Garrett (2017). It should also be noted that the symmetry condition (2.7c) at the pole is abandoned as the flow at this point is obtained by a disk-based model whereas Segalini & Garrett (2017) use the Banks (1965) series solution upto fifteenth order to initialise the solution; it is deemed that the first-order solution, which is equivalent to the rotating disk, is adequate enough as the correction terms of the series solution are negligible at the poles.

In order to solve (2.13) subject to the boundary conditions (A2)–(A10), nonlinear methods need to be considered. Typically, a Newton method is used although the resulting Jacobian would be too vast (${\sim }O(4N^2)$) to be realistically applicable, hence, another strategy is required that is more computationally efficient. An excellent candidate is the multigrid class of methods that do not require huge amounts of computational resources as the main idea of these methods is to spend the majority of computing time on smaller/coarser grids. This uses fewer nodes and theoretically speeds up the solution process by relying on error correction and obtaining accurate solutions on the coarsest grids (Henson 2003). The key aspect is that it is assumed that these solutions can be approximated at this level very accurately via a ‘smoothing’ stage that, in general, is a nonlinear iterative scheme where large errors are quickly damped. As the solution is restricted to a coarser grid, the small errors are amplified and the smoother quickly damps them. Continuing in this fashion allows the computation of an accurate solution on the coarsest grid that can be interpolated back to the finest grid.

The Full Multigrid–Full Approximation Storage (FMG-FAS) algorithm, outlined in Smith (2023), will be used to obtain solutions to (2.13), but first, (2.13a) and (2.13b) need to be discretised to generate a finite system of equations. Finite difference schemes will be implemented in order to be consistent with the Newton space-marching method used to obtain the boundary data (2.10d). Second-order centred differences will be utilised for all terms except the convective terms where a second-order upwind scheme will be employed of the form

(3.1a)\begin{align} {}[\partial_\eta\omega]^h_{i,j} &= \begin{cases} {\dfrac{3\omega_{i,j}-4\omega_{i-1,j}+\omega_{i-2,j}}{2h}}, & {\dfrac{\psi_{i,j+1}-\psi_{i,j-1}}{2h}}\geq 0,\\[10pt] -{\dfrac{3\omega_{i,j}-4\omega_{i+1,j}+\omega_{i+2,j}}{2h}}, & {\dfrac{\psi_{i,j+1}-\psi_{i,j-1}}{2h}}< 0, \end{cases} \end{align}

(3.1b)\begin{align}{}[\partial_\beta\omega]^h_{i,j} &= \begin{cases} {\dfrac{3\omega_{i,j}-4\omega_{i,j-1}+\omega_{i,j-1}}{2h}}, & -{\dfrac{\psi_{i+1,j}-\psi_{i-1,j}}{2h}}\geq 0, \\[10pt] -{\dfrac{3\omega_{i,j}-4\omega_{i,j+1}+\omega_{i,j+2}}{2h}}, & -{\dfrac{\psi_{i+1,j}-\psi_{i-1,j}}{2h}}< 0, \end{cases} \end{align}

and the nonlinear function $A^h_{i,j}:\mathbb {R}^2\rightarrow \mathbb {R}^2$ can be defined as

(3.2) \begin{equation} A^h_{i,j} = A^h(\psi_{i,j},\omega_{i,j}) = \begin{pmatrix} \omega_{i,j} + [\Delta\psi]^h_{i,j} \\ [\partial_\beta\psi]^h_{i,j}[\partial_\eta\omega]^h_{i,j} - [\partial_\eta\psi]^h_{i,j}[\partial_\beta\omega]^h_{i,j} - \epsilon[\varDelta\omega]^h_{i,j} \end{pmatrix}, \end{equation}

where $[\cdot ]^h_{i,j}$ denotes the finite difference approximation with grid spacing $h$ at the point $(\eta _i,\beta _j)$. A solution of $A^h_{i,j}=\boldsymbol {0}$ is sought and due to the nature of the upwind scheme and the finite size of the domain, points outside the boundary are needed. This is achieved using a three point Lagrange interpolant to obtain the following extrapolation formula:

(3.3)\begin{equation} f_{{-}1,j} = 3f_{0,j} - 3f_{1,j} + f_{2,j}. \end{equation}

It should be noted that if a Neumann condition is used then the extrapolated point can be determined by $f_{-1,j}=f_{1,j}$ via a centred difference.

The iterative smoother utilises both the Gauss–Seidel and Newton methods to produce a small, invertible, linear system that can be solved inexpensively and iterated over the whole grid multiple times. Note that as the viscosity is constant, then the Newton method will reduce to the Picard method as the chosen discretisation of the $\psi$ derivatives yield linear terms on any given node; this is in preparation for future work where viscous nonlinearities will be present. Thus, the following linear system is solved,

(3.4)\begin{equation} \begin{pmatrix} -4/h^2 & 1 \\ 0 & X^h_{i,j} \end{pmatrix}\boldsymbol{s}_{i,j} = A^h_{i,j}, \end{equation}

where the $2\times 2$ matrix is the Jacobian of $A^h_{i,j}$, and

(3.5)\begin{equation} X^h_{i,j} = \frac{3}{2h}\left(\left|[\partial_\beta\psi]^h_{i,j}\right|+\left|[\partial_\eta\psi]^h_{i,j}\right|\right)+\frac{4\epsilon}{h^2}>0. \end{equation}

At each node (3.4) is solved for the corrections $\boldsymbol {s}_{i,j}$ and the solutions $(\psi _{i,j},\omega _{i,j})$ are then updated via

(3.6)\begin{equation} \begin{pmatrix} \psi_{i,j} \\ \omega_{i,j}\end{pmatrix}\leftarrow\begin{pmatrix} \psi_{i,j} \\ \omega_{i,j}\end{pmatrix}+\alpha\boldsymbol{s}_{i,j}, \end{equation}

where $\alpha \in (0,1]$ is an under-relaxation parameter; in this study $\alpha =0.9$ unless otherwise stated. The boundary conditions and extrapolated points are then updated; it should be noted that on the outflow boundary ($\eta \rightarrow \infty$), (A9) and (A10) are solved instead of (3.4).

The restriction operator $\boldsymbol{\mathsf{I}}^{2h}_h$ is derived from a linear interpolant and can be interpreted as the weighted average of the surrounding points. A point on a coarser grid can be obtained by

(3.7)\begin{align} v^{2h}_{i,j} &=\frac{1}{16}(v^{h}_{2i-1,2j-1}+v^{h}_{2i+1,2j-1}+v^{h}_{2i-1,2j+1}+v^{h}_{2i+1,2j+1}) \nonumber\\ &\quad + \frac{1}{8}(v^{h}_{2i,2j-1}+v^{h}_{2i-1,2j}+v^{h}_{2i+1,2j}+v^{h}_{2i,2j+1}) + \frac{v^{h}_{2i,2j}}{4}, \end{align}

whilst an ‘injection’ is used on the boundary, i.e. $v^{2h}_{i,j}=v^{h}_{2i,2j}$. Similarly, the interpolation operator $\boldsymbol{\mathsf{I}}^{h}_{2h}$ is derived likewise,

(3.8) \begin{equation} \begin{aligned} v^{h}_{2i,2j} &= v^{2h}_{i,j}, \\ v^{h}_{2i+1,2j} &= \tfrac{1}{2}\left(v^{2h}_{i,j}+v^{2h}_{i+1,j}\right), \\ v^{h}_{2i,2j+1} &= \tfrac{1}{2}\left(v^{2h}_{i,j}+v^{2h}_{i,j+1}\right), \\ v^{h}_{2i+1,2j+1} &= \tfrac{1}{4}\left(v^{2h}_{i,j}+v^{2h}_{i+1,j}+v^{2h}_{i,j+1} +v^{2h}_{i+1,j+1}\right). \end{aligned} \end{equation}

The whole algorithm used to obtain solutions of (2.13a) and (2.13b) for the flow in the impinging region is displayed in figure 1; for more details about the FMG-FAS method, see Smith (2023).

Figure 1. Flow chart of the FMG-FAS algorithm. (Adapted from Smith 2023).

The solution is initialised on the fine grid via an initial guess; typically zero but a solution at a lower $\textit {Re}$ could be used to provide a better initial guess to aid convergence, as it should already capture the main behaviours of the solution. The $\eta$ and $\beta$ domains are discretised on a grid spacing $h_N$ with $\beta \rightarrow -\infty$ set to $-\textrm {ceil}(R_e^{1/4})$ (i.e. the smallest integer bigger than $R_e^{1/4}$), following the findings of Dennis et al. (1981) on the size of the interaction zone. The boundary at infinity, $\eta \rightarrow \infty$, is set to 30 in order to match the boundary layer flow of Segalini & Garrett (2017); and the boundary conditions (A2)–(A10) and extrapolated points (3.3) are assigned. On the coarse grid ($h=h_1$), the solution ($\boldsymbol {u}^{h}=(\psi ^h_{i,j},\omega ^h_{i,j})$) is obtained by using the smoothing operator (30 iterations) with $\alpha =0.95$. The residual is defined as $\boldsymbol {r}^{h} = \boldsymbol {f}^h-A^h(\boldsymbol {v}^h)$, where $\boldsymbol {v}^{h}=(\psi ^h_{i,j},\omega ^h_{i,j})$ denotes the current guess for the solution $\boldsymbol {u}^{h}$. If ${r}^h_{max}=\max \{\boldsymbol {r}^h\}>\textrm {tol}$, where $\textrm {tol}=10^{-5}$ is a prescribed tolerance, then the FAS part of the algorithm is initiated where the new problem $A^{2h}(\boldsymbol {u}^{2h}) = A^{2h}(\boldsymbol {v}^{2h})+\boldsymbol {r}^{2h}$ is solved instead. The simulation parameters, such as the number of iterations and the under-relaxation parameter, were chosen after much experimentation. Once the solutions $\psi ^h_{i,j}$ and $\omega ^h_{i,j}$ are obtained, the original flow variables are recovered via (2.11a,b) using second-order central finite differences except on the boundaries where second-order backwards differences are used.

It should be noted that this algorithm is not always guaranteed to converge, in the sense that ${r}^h_{max}<\textrm {tol}$, due to the typical problems that accompany nonlinear problems such as relying on a good initial guess, and over/under-shooting resulting in either cyclic or divergent behaviour (Smith 2023). In order to test the numerical method outlined above and in Smith (2023), recall that (2.13a) and (2.13b) are analogous to the planar Navier–Stokes equations in streamfunction-vorticity form but with $\textit {Re}^c=\epsilon =\textit {Re}^{-1/2}$, where $\textit {Re}^c$ denotes the computational Reynolds number. Hence, the lid driven cavity test case was used to probe the accuracy and robustness of the method in Appendix B where results were compared with similar solutions by Ghia, Ghia & Shin (1982). It is deemed that the numerical method provides accurate and consistent solutions reasonably quickly although there seems to be issues with convergence once $\textit {Re}^c>700$, in contrast to Smith (2023) who found difficulties arose when $\textit {Re}^c\sim 1000$. However, it is stated that this issue is problem dependent, meaning this value changes with the physical problem considered, for example, at higher $\textit {Re}^c$ the number of nodes on coarse grids may be insufficient to accurately capture the solution behaviour when interpolating to a finer grid Ghia et al. (1982) such as is the case when trying to resolve thin boundary layers.

Once the flow variables $U_{i,j}$ and $W_{i,j}$ are obtained, then the azimuthal velocity $V$ and pressure $\bar {P}$ can be determined. The uncoupled, linear equation (2.9c) for $V$ is solved similarly with the same interpolation and restriction operators but the smoothing stage is replaced with only the Gauss–Seidel method to solve linear equations, and utilising recursive V-cycles instead of the FMG-FAS algorithm. (2.9c) is discretised in the same manner using centred differences for the diffusive terms and backwards differences of the form (3.1) for the convective terms. The solution is smoothed (5 iterations) before being restricted to the coarse grid and is solved (20 iterations of Gauss–Seidel) before it is interpolated upwards to a finer grid where it is corrected and post-smoothed (3 iterations) repeating until the residual is less than $10^{-5}$. The pressure can be found by solving the Poisson equation (2.14) that is discretised using second-order centred differences creating a large sparse linear system. These types of problems can be easily handled by specialist linear algebra software such as MATLAB through the use of sparse matrices. For more information concerning these methods, see Smith (2023).

Lastly, in order to verify that the Stewartson (1958) model of the impinging region is consistent with other studies, the numerical simulation of the full Navier–Stokes problem (2.1) is needed. This is performed via a DNS using the Semtex software package developed by Blackburn et al. (2019) that can solve the incompressible Navier–Stokes equations in cartesian and cylindrical coordinate systems using the spectral element method, a combination of finite element and spectral methods, to exploit the geometric flexibility and higher accuracy of both respective methods. In two dimensions, Semtex uses standard finite element techniques to map the domain to two-dimensional quadrilateral elements, the Gauss–Lobatto–Legendre nodal shape function basis to achieve high spectral accuracy and continuous Galerkin projections. If the solution is sought in three dimensions, then this can be achieved by using Fourier expansions in an orthogonal direction but only if the solution is expected to be periodic, i.e. $2\frac {1}{2}$ dimensional. These qualities make Semtex an ideal choice as the flexibility of geometry combined with a robust, accurate DNS solver allows the symmetry of the sphere problem to be exploited to obtain reasonably high resolution solutions at considerably larger Reynolds numbers.

Following Calabretto et al. (2015, 2019) the computational domain consists of a large cylinder filled with fluid except at the origin (the centre of the cylinder) where a sphere, of radius one, rotates in the azimuthal direction with angular velocity $\varOmega$. This choice of the domain allows the exploitation of the axial symmetry of the sphere, and reduces the problem to a cylindrical coordinate system that is readily achieved in Semtex. Furthermore, as Semtex only solves the unsteady Navier–Stokes problem via backwards time integration schemes, then results may feature turbulent behaviour especially as $\textit {Re}$ increases that may make comparisons difficult. However, following the findings of Calabretto et al. (2019), using a spin-up time of $t_s=0.05$ and a time step $\Delta t=2.5\times 10^{-4}$, it is expected that the times $t\in [10,15]$ present steady flows close to the sphere surface, as the temporal variance of the velocity components vanished within this range. The computational model is described in more detail in Calabretto et al. (2015) and Smith (2023) where the cylindrical domain is discretised into 3872 quadrilateral elements each consisting of $10\times 10$ Lagrange knot points where more elements are located close towards the sphere surface and along the equator in order to accurately capture the boundary layer dynamics. It should be noted that although the computational model is the same as described by Calabretto et al. (2015), these simulations were independently computed using the ALICE supercomputer based at the University of Leicester, and it is believed these results build upon both their work and the work of Calabretto et al. (2019) and Dennis et al. (1981).

4. Results

The numerical methods outlined in § 3 enable the analysis of the solutions of (2.9), and subsequently the flow around the equator. Results for $\textit {Re}=10^4$ can be seen in figure 2 on a fine grid of spacing $h=2^{-4}$ corresponding to $N_\eta =481$ and ${N_\beta =1+\textrm {ceil}(\textit {Re}^{1/4})/h=161}$ nodes in each respective direction. Recall that the ${\beta \rightarrow -\infty}$ boundary data for the inflow velocities $U_{BL}$ and $V_{BL}$, and the pressure $\bar {P}_{BL}$ is obtained by solving (2.4) for the boundary layer upstream of the equator first. Rearranging (2.8) for $\theta$ and recalling that the $\beta \rightarrow -\infty$ boundary is set to $-\textrm {ceil}(\textit {Re}^{1/4})$, then ${\theta _{I}=-\epsilon \textrm {ceil}(\textit {Re}^{1/4})+{\rm \pi} /2}$, and the inflow boundary conditions can be easily set, for example, $U_{BL}(\eta )=U_{BL}(\eta,\theta _I)$.

Figure 2. The higher-order Stewartson (1958) model of the equatorial flow at $\textit {Re}=10^4$. Results are shown for (a) $\sqrt {W^2+U^2}$, (b) $\psi (\eta,\beta )$, (c) $P=\epsilon \bar {P}(\eta,\beta )$, (d) $V(\eta,\beta )$.

There is a smooth transfer of momentum seen in figures 2(a) and 2(b) corresponding to the separation and reattachment of the boundary layer along the equator. This is driven by an increase in pressure seen in figure 2(c) that creates an unfavourable pressure gradient that forces the boundary layer separation from the sphere. Note that this increase is not seen at the sphere surface in the DNS plot in figure 3(a) suggesting that the homogenous Neumann condition at $\eta =0$ in (2.15) is not appropriate. This could have been inferred from (2.5), which implies that above the equator $\partial _r P=\partial _\eta \bar {P}=V^2=\sin ^2\theta \sim 1$ at $\eta =0$. This is further supported by figure 3(b) that also shows that the behaviour of the pressure at the sphere surface is more complicated as it does not tend to 1 like $V^2$ but to a small variation of this; although setting the condition $\partial _\eta \bar {P}=1$ at $\eta =0$ seems like a good approximation that improves as $\textit {Re}$ increases. The azimuthal velocity $V$ can be seen in figure 2(d) where it is advected downstream along the equator by the reattached boundary layer forming the widely observed toroidal vortex around the equator. This is not surprising as this is the core behaviour of the advection–diffusion equation (2.9c), since it is uncoupled and there is no momentum transfer to facilitate any other behaviour.

Figure 3. Semtex results for the pressure $P$ at the equator. (a) Plot of $P$ at $\textit {Re}=10^4$. (b) Profiles of the pressure wall gradient.

By plotting only the positive radial flow, i.e. $W>0$, pockets of reverse flow can be seen in figure 4 where white areas denote where $W<0$. The Semtex simulations can be seen in figure 4(b) and how the boundary layer reattaches along the equator after the separation experienced upstream. It is clear that the boundary layer experiences an acceleration that turns it into a radial jet that traverses along the equator forming a toroidal vortex that surrounds the equatorial region of the sphere. As $\textit {Re}$ increases, a small region of reverse flow can be seen at $\beta \sim -2$ on the sphere surface that seems to grow in size with $\textit {Re}$. The pressure distribution for $\textit {Re}=10^4$ can be seen in figure 3(a) where a large increase in pressure generates an unfavourable pressure gradient forcing the boundary layer to separate from the sphere in a similar manner to figure 2(c). The Stewartson (1958) model can be seen in figure 4(a) where notably at the equator-sphere junction reverse flow is observed at $\textit {Re}=10^4$ that expands with increasing $\textit {Re}$. This is not observed in the Semtex simulations in figure 4(b) and has not been observed in any numerical studies for such a small $\textit {Re}$ suggesting that the Stewartson (1958) model is inconsistent compared with current numerics. This is amplified by noting that the shape of the boundary layer ‘tail’ is, qualitatively speaking, not the same compared with the DNS plots in figure 4(b). Additionally, a key trait of the equatorial flow is that once the boundary layer reattaches it accelerates and forms a radial jet. However, this acceleration is not observed in figure 4(a) for the Stewartson (1958) model. This is not surprising because in (2.9) there is no external forcing or coupling to the azimuthal momentum equation (2.9c), thus, momentum cannot be gained or lost resulting in the lack of speed experienced by the boundary layer. Nevertheless, this transfer of momentum is the core physical behaviour hypothesised by Stewartson (1958) and observed in Calabretto et al. (2019) before the boundary layer transitions into the radial jet.

Figure 4. Heatmaps of the positive radial velocity $W>0$. White areas denote areas where $W<0$, i.e. radial inflow. (a) Stewartson (1958) model. (b) Semtex simulations.

Furthermore, figure 4 seems to suggest that in the flow dictated by the higher-order Stewartson (1958) model, the boundary layer reattaches further downstream compared with the Semtex simulations. This is supported by figure 5(a) where at $\textit {Re}=10^4$ the boundary layer in the DNS has reattached by $\eta =5$ whereas in the Stewartson (1958) model it does not reattach until $\eta \approx 10$ such that reattachment is defined as the $\eta$ point where $W\sim 0.15$, as this is the maximum velocity attained by the boundary layer before separation, in the neighbourhood of $\beta =0$. Similarly, at $\textit {Re}=10^5$, the boundary layer from the Semtex simulations has reattached by $\eta =10$ but the boundary layer dictated by the Stewartson (1958) model has not yet reattached by this point. This is most likely due to the presence of the reverse flow at the equator-sphere junction in the Stewartson (1958) model that is absent in the DNS results. Finally, it is interesting to note that figure 5(a) displays another flaw of the model: the maximum of $W$ does not lie on the line $\beta =0$ in contrast to the Semtex simulations and other numerical studies (Dennis et al. 1981; Calabretto et al. 2019).

Figure 5. Profiles of the radial velocity $W$ close to the sphere surface ($\eta =0$). (a) Profiles along $\beta$ at stations of $\eta$. (b) Profiles along $\beta$ at stations of $\eta$.

As $\textit {Re}$ increases, it is expected that the area of reverse flow seen at the equator-sphere junction enlarges and in figure 4(a) plots of the positive radial velocity ($W>0$) can be seen for $\textit {Re}=8\times 10^4$ and $10^5$, corresponding to the same $\textit {Re}$ that were simulated by Calabretto et al. (2019). As $\textit {Re}$ increases, the pocket of reverse flow at the equator-sphere junction grows considerably larger, although this is unobserved in any of the simulations in figure 4(b) further suggesting that this model is unsuitable. However, a smaller pocket of reverse flow can be noticed to emerge from the sphere surface at $\beta \sim -6$ in figures 4(a ii) and 4(a iii) consistent with Calabretto et al. (2019); although the results in figures 4(b) and 5 place this pocket at $\beta \sim -1.5$. It does suggest that Stewartson (1958) is modelling an important characteristic of the flow that has been observed in simulations, but other issues remain including the qualitatively different shape of the ‘tail’ of the reattached boundary layer and that the maximum of $W$ is not always on the line $\beta =0$. These are clues that the current iteration of the model is not describing the flow reasonably accurately.

The vorticity at $\textit {Re}=10^5$ can be seen in figure 6(a) that unsurprisingly varies greatly along the boundary layer trajectory due to the variation of the streamwise velocity causing local rotation into the boundary layer. Interestingly, on the sphere surface a small region of positive vorticity is observed. This suggests anti-clockwise flow, hinting at the presence of a small vortex at the equator-sphere junction that was hypothesised to occur by Stewartson (1958). This is further implied by figure 6(b) depicting the azimuthal velocity $V$ that is advected further along the equator by the reattached boundary layer as $\textit {Re}$ increases. Hence, expanding the reach and strength of the toroidal vortex where $V\sim 0.7$ as $\eta \rightarrow \infty$ compared with $V\sim 0.5$ for $\textit {Re}=10^4$ as seen in figure 7(a). Around the equator, there is a small area $([0,5]\times [0,5])$, where $V$ is smaller in magnitude than expected, meaning that one would expect a strict monotonic decrease of $V$ in a similar fashion as the $\textit {Re}=10^4$ case as depicted in figures 2(d) and 7. However, as can be seen in figure 7(a), there is a small domain where $V$ is constant suggesting a small vortex where the higher strength swirl is advected upwards whereas the relatively weaker swirl is advected downwards creating this disparity around the equator.

Figure 6. The higher-order Stewartson (1958) model of the equatorial flow at $\textit {Re}=10^5$.

Figure 7. Profiles of the azimuthal velocity $V$ along $\eta$ at stations of $\beta$.

The presence of a vortex at the equator-sphere junction can be confirmed by plotting the radial reverse flow ($W<0$) and vector field in figure 8(a) and the latitudinal skin friction, $\textit {Re}^{-1/2}\partial _r U = \partial _\eta U$ at $\eta =0$ in figure 9. Due to the difference in magnitude of the reverse flow, as it is much smaller than the streamwise flow, it is difficult to visualise the vector fields in figure 8. However, as the latitudinal skin friction is negative, then this must imply $U<0$ near the sphere surface due to the no-slip condition $U=0$. Hence, there must be a vortex at the equator-sphere junction for the Stewartson (1958) model as hypothesised. However, it is clear that this argument does not apply to the Semtex simulations where $\partial _\eta U>0$, meaning that $U>0$ even at larger $\textit {Re}$ and, thus, no vortex can exist. Therefore, one must conclude that this vortex is not physical.

Figure 8. Heatmaps of the negative radial velocity $W<0$ and the vector field of the planar motion at $\textit {Re}=10^5$. The blue region refers to positive radial flow ($W>0$) with respect to figure 4. (a) Stewartson (1958) model. (b) Semtex simulations.

Figure 9. Profiles of the latitudinal skin friction $\textit {Re}^{-1/2}\partial _r U = \partial _\eta U$.

5. Discussion

A new model of the steady flow around the equator for the rotating sphere has been obtained that features a higher-order version of the Stewartson (1958) model. The main characteristics of the flow are modelled well, namely the separation and subsequent reattachment of the boundary layer through the smooth transfer of momentum driven by an increase in pressure at the equator leading to an unfavourable pressure gradient. Also, there is a small pocket of reverse flow that grows with increasing $\textit {Re}$ that eventually turns into a small vortex at the equator-sphere junction. This is not seen in other numerical studies or the simulations computed here and, thus, it is deemed that this vortex is not physical, but there is a smaller pocket of reverse flow that has been observed at corresponding $\textit {Re}$ by Calabretto et al. (2019) of approximately the same magnitude. On the other hand, there are significant issues: the qualitative shape of the ‘tail’ of the radial velocity, $W$, does not match with those seen in the DNS solutions in figure 4(b); the boundary layer reattaches along the equator much earlier, with respect to $\eta$; the maximum of $W$ is not on the line $\beta =0$ as was found by Dennis et al. (1981); and most importantly, the presence of behaviour not seen in other DNS studies, i.e. a large pocket of reverse flow that forms a small vortex with increasing $\textit {Re}$, suggesting that the model is not reasonably accurate. These problems are exemplified by the obvious differences with the DNS velocity profiles seen throughout, for example, in figure 5, both the positions of reattachment and reverse flow are different as well as the fact that the maximum of $W$ does not lie on $\beta =0$. This is confirmed in figure 10(b) where the absolute errors of each component are presented for $\textit {Re}=10^4$ and $\textit {Re}=10^5$. It is expected that the error decreases with increasing $\textit {Re}$ due to the asymptotic nature of the derivation of the model, however, disconcertingly, the error for the $U$ and $V$ components increases, whilst only slightly decreasing for the $W$ component. Additionally, the magnitudes of the errors are of the same order as the velocity components themselves, which is extremely concerning. The error of the radial velocity, $W$, can be attributed to the lack of speed gained by the boundary layer. This is most likely due to the dismissal of apparent negligible terms in the derivation of (2.9) that are in fact important as discussed later. This again suggests that the model is not tremendously accurate, especially considering the significant size of the error of the azimuthal velocity. Nevertheless, a small pocket of reverse flow is seen at larger Reynolds numbers that does suggest that the model needs slight modifications.

Figure 10. Absolute errors of the velocity components. Here $U_{DNS}$ denotes the DNS solution and $U_{IMP}$ denotes the higher-order Stewartson (1958) solution for the equatorial flow. Results are shown for (a) $\textit {Re}=10^4$. (b) $\textit {Re}=10^5$.

It should be noted that as the Semtex simulations are unsteady, and the flow is intrinsically unsteady, then to open discussion to unsteady behaviour is natural. As first discussed by Banks & Zaturska (1979), then subsequently by Stewartson & Simpson (1982), Dennis & Duck (1988), Van Dommelen (1990) and Calabretto et al. (2015, 2019), the presence of a finite-time singularity in the unsteady boundary layer equations provides the mechanism for the initial separation of the boundary layer. Physically, this singularity is observed as the radial jet (Calabretto et al. 2015). As this study is primarily concerned with the post-collisonal flow, in other words the flow post-singularity, in order to investigate the dynamics at large Reynolds numbers, then comparison to prior unsteady studies shall not be conducted. Nevertheless, it would be interesting to compare the early time dynamics of the numerical simulations to these other studies, in particular Van Dommelen (1990) to test their more accurate Lagrangian approach.

It is clear that the model of the flow at the equator-sphere junction models the core qualitative behaviour well: the separation and subsequent reattachment of the boundary layer, and presence of a small region of reverse radial flow. However, it also possesses various issues including the size and location of said reverse flow; the development of a small unphysical vortex at the equator-sphere junction, hypothesised by Stewartson (1958); and the shape and size of the magnitude of the radial velocity. These can be observed by the DNS simulations where the velocity profiles behave differently compared with the Stewartson (1958) model, and the absolute errors of each velocity component seem to increase with $\textit {Re}$ instead of decreasing. Nevertheless, the equations proposed by Stewartson (1958) have been solved for the first time despite these problems using the novel multigrid method of Smith (2023). It is also clear that the model, despite its problems, must be similar to a more consummate model due to the nature of its derivation and the appearance of a small pocket of reverse flow seen at the correct Reynolds numbers. It is hypothesised that terms dismissed in the derivation must be more influential than first thought. All current models of the equatorial flow assume a totally flat space, i.e. no curvature; however, at the equator one can visualise the geometry akin to a cylinder. Hence, it is postulated that perhaps azimuthal curvature terms in (2.1) are significant.

This is readily seen by considering that the planar velocity components, $U$ and $W$, are of equal order $\epsilon ^\gamma$, to facilitate the separation and reattachment of the boundary layer. By assuming $U=\epsilon ^\gamma \bar {U}$ and $W=\epsilon ^\gamma \bar {W}$, and keeping the same scaling assumptions, i.e. $P\sim O(\epsilon )$ and $V\sim O(1)$, then the radial momentum equation (2.1a) at leading order reduces to

(5.1)\begin{equation} \bar{W}\partial_\eta\bar{W}+\bar{U}\partial_\beta\bar{W}-\epsilon^{1-2\gamma} V^2={-}\epsilon^{1-2\gamma}\partial_\eta\bar{P}+\epsilon^{1-\gamma}(\partial_\eta^2\bar{W}+\partial_\beta^2\bar{W}). \end{equation}

Noticeably, the additional term $-V^2$ appears, and by increasing $\gamma$ from 0 it is clear that this additional term and the pressure gradient become more influential, both becoming of $O(1)$ at $\gamma =1/2$. This adds coupling to (2.9c) and provides a mechanism to allow momentum transfer allowing the reattached boundary layer to gain speed and transition into a radial jet. It is interesting to note that the additional term is also present in the Navier–Stokes equations when expressed in a cylindrical coordinate system, further suggesting that azimuthal curvature is important. This aspect will be one of the focus points of another paper that is currently in preparation that builds much on the content presented here where results are greatly improved and many of the discrepancies of the Stewartson (1958) model discussed disappear.

Lastly, a major discussion point in Dennis et al. (1981) was the apparent discrepancy between the size of the interaction zone at the equator-sphere junction hypothesised by Smith & Duck (1977). Smith & Duck (1977) posit a large region of size $O(\textit {Re}^{-3/14})$ corresponding to an azimuthal skin friction, $\partial _r V/\sqrt {\textit {Re}}\sim O(\textit {Re}^{-2/7})$; whilst Dennis et al. (1981), through numerical simulations, find a correlation of ${\sim }O(\textit {Re}^{-1/4})$. Both scalings are relatively close, thus, results at higher $\textit {Re}$ are needed to observe any reasonable deviation. By calculating $\partial _r V/\sqrt {R_e}$ at the equator for various $\textit {Re}$ and combining with results of Dennis et al. (1981) provide the figures in table 1. As $\textit {Re}$ increases, the magnitude of the skin friction decreases; this is not unexpected as the influence of viscosity decreases with an increase in $\textit {Re}$. Dividing the figures in table 1 by the proposed scaling and taking the average yields an approximation for the coefficients to obtain $\partial _r V/\sqrt {\textit {Re}}\sim 1.0883\textit {Re}^{-2/7}$ and $\partial _r V/\sqrt {\textit {Re}}\sim 0.779\textit {Re}^{-1/4}$. Note that the coefficient for the $\textit {Re}^{-1/4}$ scaling is the same calculated by Dennis et al. (1981) up to two decimal places. A $\log$-$\log$ plot of these lines with both the results of the higher-order Stewartson (1958) model and DNS can be seen in figure 11(a). The DNS results correlate well with the $\textit {Re}^{-1/4}$ scaling compared with the $\textit {Re}^{-2/7}$ scaling, although it is still rather close between the two, whereas the higher-order Stewartson (1958) model is wildly off correlation, increasing for larger $\textit {Re}$ and further highlighting the lack of reasonable accuracy of the model. Additionally, the relative error between the DNS results and the proposed scalings yields figure 11(b). The error of the scaling $\textit {Re}^{-2/7}$ increases with $\textit {Re}$ whereas the error of the scaling $\textit {Re}^{-1/4}$ seems to decrease, albeit slowly. Furthermore, the error of the $\textit {Re}^{-1/4}$ scaling is less than $5\,\%$ whereas, for the $\textit {Re}^{-2/7}$ scaling, the error is rarely below this threshold. Note that the relative error for the higher-order Stewartson (1958) model was $>100\,\%$ for $\textit {Re}>8\times 10^4$. Hence, figure 11 suggests that $\partial _rV/\sqrt {\textit {Re}}\sim 0.779\textit {Re}^{-1/4}$, consistent with Dennis et al. (1981), and it is deemed that the interaction zone is of size $O(\textit {Re}^{-1/4})$.

Table 1. Direct numerical simulation values of the azimuthal skin friction at $\theta ={{{\rm \pi} }/{2}}$ for various $\textit {Re}$.

Figure 11. Values of the azimuthal skin friction $-{{\partial _r V}/{\sqrt {\textit {Re}}}}$ at $\theta ={{{\rm \pi} }/{2}}$. (a) A $\log$-$\log$ plot of $-{{\partial _r V}/{\sqrt {\textit {Re}}}}$. (b) Relative error (%).