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 Reference Segalini and Garrett2017). 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 (Reference Howarth1951), Banks (Reference Banks1965) and, more recently, Segalini & Garrett (Reference Segalini and Garrett2017), 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 (Reference Stewartson1958) describes the flow qualitatively well (Calabretto et al. Reference Calabretto, Denier and Levy2019), hence, following Stewartson (Reference Stewartson1958), 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 Reference Segalini and Garrett2017), (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 (Reference Stewartson1958) 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 (Reference Stewartson1958) 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 (Reference Segalini and Garrett2017) but with the introduction of the Neumann conditions for the equator, following Dennis et al. (Reference Dennis, Ingham and Singh1981) 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 (Reference Stewartson1958), 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 (Reference Segalini and Garrett2017), 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 (Reference Segalini and Garrett2017). 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 (Reference Segalini and Garrett2017) use the Banks (Reference Banks1965) 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 Reference Henson2003). 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 (Reference Smith2023), 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 (Reference Smith2023).

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

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. (Reference Dennis, Ingham and Singh1981) 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 (Reference Segalini and Garrett2017); 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 Reference Smith2023). In order to test the numerical method outlined above and in Smith (Reference Smith2023), 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 (Reference Ghia, Ghia and Shin1982). 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 (Reference Smith2023) 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. (Reference Ghia, Ghia and Shin1982) 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 (Reference Smith2023).

Lastly, in order to verify that the Stewartson (Reference Stewartson1958) 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. (Reference Blackburn, Lee, Albrecht and Singh2019) 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. (Reference Calabretto, Levy, Denier and Mattner2015, Reference Calabretto, Denier and Levy2019) 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. (Reference Calabretto, Denier and Levy2019), 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. (Reference Calabretto, Levy, Denier and Mattner2015) and Smith (Reference Smith2023) 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. (Reference Calabretto, Levy, Denier and Mattner2015), 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. (Reference Calabretto, Denier and Levy2019) and Dennis et al. (Reference Dennis, Ingham and Singh1981).