On high Taylor number Taylor vortices

Taylor-Couette flow (TCF) is perhaps among the most studied flows in fluid mechanics. In the 100 years since Taylor’s monumental work (Taylor 1923), it has provided an excellent testing ground for theoretical, experimental and numerical studies of rotating shear flows. How shear and Coriolis forces alter flow characteristics is important in various applications, and TCF was designed so that they can easily be adjusted by changing the rotation speed of the inner and outer cylinders. Researchers have long been fascinated by the numerous metastable flow patterns observed in the relatively low Taylor number regime (e.g. Andereck et al. 1986). On the other hand, it was only a decade ago that the study of high Taylor number flows became active. Great efforts were made to investigate the nature of turbulence in the parameter space by means of high Taylor number experiments (e.g. Paoletti & Lathrop 2011; Dennis et al. 2011; van Gils et al. 2011, 2012; Huisman et al. 2012) and direct numerical simulations (e.g. Ostilla et al. 2013; Brauckmann & Eckhardt 2013, 2017; Ostilla-Mónico et al. 2014ab). As summarised in the review paper by Grossmann et al. (2016), and in fact seen in the pioneering experiments by Lathrop et al. (1992), Lewis & Swinney (1999), fully developed turbulence has a surprisingly clean asymptotic character, while there seems to be no definitive Navier-Stokes-based theory to explain it.

This study aims to reveal the asymptotic properties of steady axisymmetric solutions at high Taylor numbers and to compare them with the experimental and numerical results. The analysis of such solutions, known as Taylor vortex solutions, goes back to the weakly nonlinear analysis, for example, by Davey et al. (1962). With modern computational power, it is possible to calculate solutions up to the Taylor number used in the experiments. Of course, the use of Newton’s method is essential as the solution is unstable in the high Taylor number regime.

The idea that an unstable solution with a relatively simple structure with respect to time can capture some characteristics of turbulence is not as absurd as one might think. It is well-known in dynamical systems theory that chaotic dynamics can be approximated by a sufficiently large number of periodic orbits (see Cvitanović et al. 2012, for example). For moderate Reynolds number shear flows, it was repeatedly reported that a good approximation of the turbulent dynamics can be obtained with a small number of periodic orbits (Kawahara et al. 2012; Willis et al. 2013; Krygier et al. 2021). In phase space, these periodic orbits usually appear with stationary or travelling wave solutions in their vicinity. The advantage of focusing on simple solutions is that their asymptotic nature may be justified theoretically by mathematical analyses. Over the past decade, it has been established that the matched asymptotic expansion is a powerful tool in understanding the behaviour of steady or travelling wave solutions in shear flows (Hall & Sherwin 2010, Deguchi et al. 2013, Deguchi & Hall 2014ab, Deguchi 2015, Dempsey et al. 2017, Deguchi & Walton 2018). Therefore if we are allowed to assume that there is a simple solution that roughly captures the scaling properties of turbulence within the vast phase space, then there is hope for a logical explanation of the scaling from first principles.

In the high Taylor number numerical and experimental studies, the parameter dependence of angular momentum transport was a major focus. The driving force behind those studies was the ‘analogy’ between turbulent Rayleigh-Bénard convection (RBC) and TCF. This analogy was well known at least in the 1960s and has risen and fallen throughout the history of turbulence research (Bradshaw 1969, Dubrulle & Hersant 2002). The recent studies have been influenced by Eckhardt et al. (2007) who argued that the phenomenology of RBC turbulence proposed by Grossmann & Lohse (2000) can be applied to TCF as well. Shortly later the aforementioned high Taylor number TCF experiments confirmed that the scaling of the Nusselt number $Nu$ is similar to that observed for RBC by He et al. (2012) ($Nu$ for TCF is defined as the torque on the cylinder wall normalised by its laminar value). It should be remarked that despite the analogy that has been believed, this similarity in the ultimate scaling is actually not at all obvious. As pointed out, for example, by Chandrasekhar (1961), Robinson (1967), Veronis (1970), and Lezius & Johnston (1976), for the two flows to be equivalent they must be at least axisymmetric. Moreover, the exact equivalence of the two flows requires an infinitesimally narrow cylinder gap, co-rotating cylinders, and a Prandtl number of unity.

It is an interesting question, then, to forget the latter three conditions and ask whether the analogy in the sense of the $Nu$ scaling holds for an axisymmetric Taylor vortex and a two-dimensional roll cell. The equations governing both flows are not the same, but they certainly have a similar structure. Many researchers have theoretically studied the large Rayleigh number nature of roll cells in RBC over the years; see Pillow (1952), Robinson (1967), Wesseling (1969), Chini & Cox (2009), Hepworth (2014) for constant Prandtl number flows, Roberts (1979), Jimenez & Zufiria (1987), and Vynnycky & Masuda (2013) for asymptotically large Prandtl number flows. Waleffe et al. (2015) and Sondak et al. (2015) shortened the wavelength of the RBC roll cell and found that there is a special wavelength at which the Nusselt number reaches a maximum value. Interestingly, this optimised Nusselt number is close to that obtained experimentally, although the Rayleigh number they used is much lower than the one used by He et al. (2012). More recently, Wen et al. (2022) continued the same solution branch to higher Rayleigh numbers and claimed that the maximum Nusselt number corresponds to the so-called classical scaling, where the Nusselt number is proportional to the Rayleigh number to the power of one-third (Malkus et al. 1954; Priestley 1954; Grossmann & Lohse 2000; Kawano et al. 2021). The study by Kooloth et al. (2021) confirmed that roll cells with different wavelengths play important roles in two-dimensionally restricted RBC turbulence. This paper is motivated by all the above RBC studies.

In the classical turbulence regime of TCF, where Taylor vortices are robustly observed in experiments, the symmetry restriction of the flow may not be necessary for a good agreement between the steady solution and turbulence, and at least a better agreement than in RBC can be expected. The situation is different in the ultimate turbulence regime, where eddies of different sizes and wavelengths are present, making the structure far more complex than classical turbulence. However, as long as the cylinders are not strongly counter-rotating, the experimental observations indicate that vortices of about the scale of the gap are still present, suggesting that Taylor vortices may play some role in the dynamics.

It should also be noted that previous theoretical studies have shown that analytical approximations can be obtained for vortices with extremely short wavelengths driven by thermal or Coriolis forces (Hall & Lakin 1988; Bassom & Hall 1989; Bassom & Blennerhassett 1992; Denier 1992; Blennerhassett & Bassom 1994). In this type of asymptotic theory, the mean flow varies by a finite amount from the base flow and is therefore called a strongly nonlinear theory. However, the scaling of momentum and heat transport of the nonlinear state is not so different from that of the basic flow, and in this sense, the character of fully developed turbulence is not well captured. Attempts have been made to construct asymptotic solutions with larger amplitudes, but a complete understanding is still lacking. This paper proposes a solution to this long-standing problem.

In the next section, we begin by formulating our problem. Comparisons of the Taylor vortex solutions with previous experiments and simulations are then carried out in section 3. In section 4, a matched asymptotic expansion analysis is performed for the case of Taylor vortices with an aspect ratio about unity. We will see in section 5 that the asymptotic structure of the solution dramatically changes when the axial period becomes asymptotically short. In section 6, how the short-period vortices develop from the laminar solution is investigated in detail using a matched asymptotic expansion. Finally, in section 7, the main findings are summarised and discussed.

Taylor-Couette flow can be described by the Navier-Stokes equations in the cylindrical coordinates $(r,\theta,z)$. If the flow is axisymmetric, the governing equations are written as

The operators $D$ and $\triangle$ are defined as $D=\partial_{t}+u\partial_{r}+w\partial_{z}$ and $\triangle=\partial_{r}^{2}+r^{-1}\partial_{r}+\partial_{z}^{2}.$ The length and velocity scales are chosen so that the cylinder gap is unity and the no-slip conditions on the cylinder walls are described as

using the Reynolds numbers associated with the rotation of the inner and outer cylinders, $R_{i}$ and $R_{o}$. Note that our non-dimensionalisation implies that using the radius ratio $\eta=r_{i}/r_{o}<1$ the inner and outer radii are specified as

respectively. The circular Couette flow solution is written as

For other nontrivial solutions of (1) periodicity is imposed in the interval $z\in[0,2\pi/k]$ using the axial wavenumber $k$. The momentum equations (1a)-(1c) can be simplified as

Steady solutions of the above system can be computed without regarding their stability by using the numerical code used in our previous studies (Deguchi & Altmeyer 2013; Deguchi et al. 2014). The code is based on the Newton-Raphson method applied to the Chebyshev-Fourier discretised system. To calculate the Taylor vortex with an aspect ratio of about unity, we typically used up to 250th Chebyshev polynomials and 250th Fourier harmonics. This spatial resolution is more than sufficient for $Ta=O(10^{10})$. We shall see that when $k$ is large, we can reach higher Taylor numbers. In this case, the highest degree of Chebyshev polynomials was increased to 450 to fully resolve the very thin near-wall boundary layer structures.

Instead of the two Reynolds numbers, the majority of high Taylor number TCF studies summarised in Grossmann et al. (2016) used the Taylor number $Ta$ and the rotation rate $a$. They are easily found by the standard parameters as

The former parameter is similar to the shear Reynolds number used in Dubrulle et al. (2005) and is zero for the rigid body rotation case (their second parameter, the rotation number, is a function of $a$ and $\eta$). The Nusselt number is defined by

where

is the mean azimuthal velocity.

In the limit $\eta\rightarrow 1$, the gap becomes much narrower than the cylinder radius, and the local flow can be represented in Cartesian coordinates. As noted by Deguchi (2016) there are several variations on the narrow gap limit, two of which are relevant to this paper. One is of course the rotating plane Couette flow, where the system is in perfect agreement with RBC if the Prandtl number is unity (see Chandrasekhar 1961, Robinson 1967, Veronis 1970, and Lezius & Johnston 1976, Brauckmann & Eckhardt 2017, for example). The other utilises an argument similar to the derivation of the Görtler vortex (Hall 1983), and is used, for example, by Denier (1992). For completeness, the difference between the two limits is highlighted in Appendix A.

Here we use the radius ratio $\eta=5/7$ and fix the outer cylinder ($a=0$) to compute the Taylor vortex solutions. Significant deviations from the narrow gap approximation can be observed at this radius ratio. That parameter choice is frequently used in experiments and numerical simulations (see section 3 of Grossman et al. 2016) and is hence convenient for the comparison. In experiments, it has been confirmed that the effect of endwalls on torque/mean flow measurements is negligible (see van Gils et al. (2012), for example).

The red curve in figure 1 shows the Taylor vortex solution calculated with the fixed axial wavenumber $k$. According to the simulations, the Taylor cells are relatively robust even in the numerically generated classical turbulent flows, with a cell aspect ratio of approximately unity. More specifically, direct numerical simulations (DNS) by Ostilla et al. (2013) imposed the axial periodicity $2\pi$ and typically observed three vortex pairs (i.e. $k=3$ modes). This motivated the choice of $k=3$ for our Taylor vortex computation.

The solution branch bifurcates from the circular Couette flow at $Ta\approx 10^{4}$, and as the Taylor number increases, it loses stability with respect to three-dimensional perturbations. The onset of turbulence is $Ta\approx 3\times 10^{6}$. The blue circles in the figure are the DNS by Ostilla et al. (2013) and Ostilla-Mónico et al. (2014a). Despite the relatively short axial periodicity imposed, the simulations are known to agree with the experimental results (the squares and triangles in figure 1). Both the experiments and simulations indicate the existence of a transition point $Ta\approx 10^{8}$ at which the behaviour of $Nu$ changes. The turbulence below/above the transition point is referred to as the classical/ultimate regime. The Nusselt number of the ultimate turbulence has the scaling $Nu\propto Ta^{\beta}$ with the exponent $\beta\approx 0.38$, which is also seen in the RBC experiments (He et al. 2012). The empirical asymptotic prediction $Nu=0.009\,Ta^{0.38}$ sits well below the theoretical upper bound of $Nu$ derived by Ding & Marensi (2019). The exponent 1/2 is typical for upper bounds using the energy method. Similar asymptotic behaviour of $Nu$ but with a logarithmic correction has been proposed using the phenomenology of the log law of turbulent boundary layers (see Grossmann et al. 2016). It is, however not known whether such scaling can be clearly observed in experiments.

As long as the solution has an aspect ratio of about unity, it captures the nature of classical turbulence surprisingly well. This is best illustrated by a comparison of mean flows shown in figure 2a. Here

is the shifted and normalised mean angular velocity (denoted by $\langle\overline{\omega}\rangle_{z}$ in Ostilla et al. 2013). At sufficiently high Taylor numbers, the boundary layer formation near the cylinder walls is clearly visible. As seen in figure 1, the solution gives a reasonable estimate of $Nu$ in the classical turbulence regime. The cross-section of the Taylor vortex (see figure 3b) reveals that the boundary layer actually surrounds the vortex core, where the azimuthal velocity appears to be rather uniform. The dynamics of the boundary layers in the classical turbulence are known to be somewhat quiescent (see Ostilla et al. 2013, for example), and their qualitative structure is reminiscent of figure 2b. The literature summarised in Grossman et al. (2016) often implies that the Prandtl-Blasius theory could be applied to the boundary layer, and we shall see in section 4 that this is, in fact, true for the asymptotic limit of the steady solution. Moreover, the theory to be presented in section 4 yields the scaling of $Nu\propto Ta^{0.25}$ and $Re_{w}\propto Ta^{0.5}$, which are well agreeing with the turbulent observations. Here $Re_{w}$ is the wind Reynolds number, i.e. the typical roll cell circulation speed normalised by the viscous velocity scale $\nu/d$ ($\nu$ is the kinematic viscosity of the fluid and $d$ is the cylinder gap). The precise definition of the wind Reynolds number varies according to the literature. Huisman et al. (2012) measured the standard deviation of the radial velocity $u$, while Ostilla et al. (2013) computed the average of $u^{2}+w^{2}$ and then square rooted it. Both definitions give the same scaling, but with different prefactors.

In the ultimate turbulence, on the other hand, the situation appears to be more intricate. The snapshots of turbulence are typically accompanied by eddies of a vast scale, whereby the large-scale Taylor vortex is blurred in the time-averaged field. In particular, small turbulent eddies appearing in the boundary layer have been pointed out as a critical qualitative difference between the two turbulent regimes separated by the transition point seen in figure 1 (Dong 2007, Ostilla-Mónico et al. 2014a). Therefore solutions with larger wavenumbers may play a more critical role in the turbulent dynamics. Our Nusselt number results also support this speculation. As seen in figure 1, the Taylor vortex with fixed $k$ underestimates experimental $Nu$ in the ultimate regime. However, if $k$ is optimised to maximise $Nu$ (the crosses in figure 1), it can reach the experimental values for $Ta\lesssim 10^{11}$.

Figure 3a shows how $Nu$ changes as the wavenumber $k$ of the Taylor vortex is varied. The $Nu$ curve has two local maxima. The bimodal distribution of $Nu$ and the scaling of the two peaks are remarkably similar to those seen in the RBC roll cell computation by Waleffe et al. (2015) and Sondak et al. (2015). Based on their calculations at several Rayleigh numbers $Ra$, the first (second) peak has the wavenumber scaling $k\propto Ra^{0.217}$ ($Ra^{0.256}$) with the Nusselt number scaling exponent $\beta$ slightly larger (smaller) than 0.31. We shall provide a theoretical rationale for the scaling in sections 5 and 6.

For the ultimate turbulence regime, the approximation of the turbulent mean flow by the $k=3$ state is slightly worse (figure 2b). However, the mean flow of course varies with $k$. As the value of $k$ is increased from 3, the mean flow of the solution in the core region approaches a turbulent profile in the vicinity of the second peak. Figures 3c,d show the flow field at each peak, which is also examined in detail in the following two sections.

In figure 2a, we saw that some kind of homogenisation occurs in the core region of the Taylor vortex, and a thin boundary layer emerges around it. Such a flow structure looks very much like the high Rayleigh number RBC roll cell, which was studied intensively from the 1950’s to the 1970’s, for example, by Pillow (1952), Robinson (1967), and Wesseling (1969). More recently, it was reported that when the slip walls are imposed, the asymptotic solution can be calculated semi-analytically and is in good agreement with the numerical solutions (Chini & Cox 2009; Hepworth 2014). However, as noted by Robinson (1967), for no-slip walls, the scaling and structure of the boundary layer have to be modified, which makes analytical progress difficult. The situation in the Taylor vortex is close to the latter problem, as we shall see below.

First, let us examine the core region. For RBC, the roll cell vorticity in this region is known to become constant due to the Prandtl-Batchelor theorem, which states that the homogenisation of the vorticity occurs in the region where the streamlines are closed in two-dimensional inviscid flows (Prandtl 1904; Feynman & Lagerstrom 1956; Batchelor 1956). Moreover, the temperature becomes constant as well because the temperature equation has a similar structure to that of the vorticity equation (see Moore & Weiss 1973, for example). To better see what precise physical quantities are homogenised in the core of the Taylor vortex, it is desirable to choose large gaps. Figure 4 shows the results for $\eta=0.5$. The colour map shown in figure 4b clearly indicates that it is actually angular momentum $rv$ that becomes uniform in the region where the streamlines are closed. Also, figure 4c shows that it is not the azimuthal vorticity $\omega$ that homogenises, but $\omega/r$.

Building on the above observations, now we introduce the new variables $\Gamma=Ta^{-1/2}rv$ and $\Omega=\omega/r$, whereby the governing equations (6) are rewritten as

Here, the subscript $r$ or $z$ denotes a partial differentiation. The no-slip conditions become

If $\Gamma$ is considered as temperature, $\Omega$ as roll cell vorticity and $Ta$ as Rayleigh number, one may notice that the structure of the equations is very similar to that of the two-dimensional RBC. The last term in the right hand side of (4.1b) corresponds to the Coriolis force when the RPCF limit is taken, and it plays a role of buoyancy. (This term is not exactly the Coriolis force unless the system is in the RPCF limit, but for simplicity we will call it ‘Coriolis force’ hereafter; see Appendix A also.) This analogy in the sense of the structure of the equations would explain why $\Gamma$ and $\Omega$ become constant in the Taylor vortex core, as seen in figure 4, at least on an intuitive level.

In fact, following Batchelor (1956), a mathematical argument can be developed. First, assuming the typical roll cell circulation strength $Re_{w}$ (the wind Reynolds number) is asymptotically large, we expand $\Psi=Re_{w}\Psi_{0}+\cdots,\Gamma=\Gamma_{0}+\cdots,\Omega=Re_{w}\Omega_{0}+\cdots$. Then the leading order part of (11a) becomes $\Psi_{0r}\Gamma_{0z}-\Psi_{0z}\Gamma_{0r}=0$, which suggests that the function $\Gamma_{0}$ only depends on $\Psi_{0}$. Now, in the $r$-$z$ plane, take a region $\mathcal{A}$ enclosed by a streamline $\mathcal{C}$ oriented counterclockwise (thus $\Psi$ is constant along $\mathcal{C}$). Integrating (11a) over $\mathcal{A}$, we obtain

The left hand side vanishes because upon using the Stokes theorem it becomes

The second equality comes from the fact that the gradient of $\Psi$ and the line element on $\mathcal{C}$, $d\mathbf{l}$, are orthogonal on the streamline. While the leading order part of the right hand side can be transformed by again applying the Stokes theorem:

Here $(u_{0},w_{0})=(-r^{-1}\Psi_{0z},r^{-1}\Psi_{0r})$ is the leading order part of the roll velocity. The integral in the last line should not vanish because we are assuming the existence of a strong circulation due to the swirling motion of the Taylor roll. Thus the equation implies $\frac{d\Gamma_{0}}{d\Psi_{0}}=0$ at the value of $\Psi_{0}$ on ${\mathcal{C}}$ we choose. The above argument holds for any region enclosed by a streamline, and hence the value of $\Gamma_{0}$ must be constant in the core region. The argument above does not change when an arbitrary constant is added to $\Gamma$. This implies that if power series asymptotic expansion of $\Gamma$ is considered the homogenisation also occurs in all the higher-order terms as long as they have a spatial scale of $O(1)$. Therefore for steady solutions, the non-uniform component in $\Gamma$ is exponentially small.

Likewise we can show that the value of $\Omega_{0}$ is constant as well in the core. By integrating (11b) over ${\mathcal{A}}$,

and of course the left hand side should vanish. Since we already know that $\Gamma_{0}$ is a constant plus an exponentially small fluctuation, the last term in the right hand side of (16) can be neglected. From (11b) we see that $\Omega_{0}$ is a function of $\Psi_{0}$, and thus the leading order part of integral (16)

yields $\frac{d\Omega_{0}}{d\Psi_{0}}=0$.

In the above argument, we have assumed that the swirling speed of the rolls is sufficiently large, but to see exactly how large it is, we need to analyse the boundary layer. Let us now take the region $\mathcal{A}$ as large as possible and assume that a viscous boundary layer appears around its boundary $\mathcal{C}$. In this core region (see figure 5a), we have the estimation $Re_{w}=O(\Psi|_{c})=O(\Omega|_{c})$, where the subscript $c$ indicates that the physical quantities are measured in the core. At the roll cell perimeter $\mathcal{C}$ without loss of generality, we can set $\Psi=0$. Hence, if the thickness of the boundary layer is written as $\epsilon$, the size of the streamfunction there can be estimated as $O(\Psi|_{b})=O(\epsilon\Psi|_{c})=O(\epsilon Re_{w})$ (the subscript $b$ stands for boundary layer). In order to ensure the viscous-convective balance of (11b) in this thin layer we further require $O(\Psi|_{b})=O(\epsilon^{-1})$, and therefore $Re_{w}=O(\epsilon^{-2})$.

As seen in figure 4c, parts of the boundary layer are attached to the walls, where the flow has to fulfil the no-slip conditions. In order for the streamfunction to be modified in the near-wall boundary layer, both sides of (11c) must balance, so $O(\Omega|_{b})=O(\epsilon^{-1}Re_{w})$ using $\partial_{r}=O(\epsilon^{-1})$ and $O(\Psi|_{b})=O(\epsilon^{-1})$. Another essential part of the boundary layer is the plume, where the layer is detached from the wall. When the boundary layer leaves the wall, it typically forms a sharp corner. Similar to the RBC cases, as we round the corner, the sizes of $\Omega$ and $\Gamma$ are unchanged. This can be justified by considering a streamline within the boundary layer. The streamline passes through the $O(\epsilon^{1/2})$ neighbourhood of the corner, where the flow is inviscid, and hence $\Omega$ and $\Gamma$ are functions of the streamfunction; see also the RBC literature introduced at the beginning of this section.

The plume layer is no longer parallel to the wall, and the Coriolis force acting there drives the whole Taylor cell. This plume balance allows us to completely fix the flow scaling in terms of $Ta$. Assuming $\partial_{r}=O(1)$ and $\partial_{z}=O(\epsilon^{-1})$ in (11b), the viscous-convective terms of $O(\epsilon^{-2}\Omega|_{b})=O(\epsilon^{-5})$ counter-balances the Coriolis term of $O(\epsilon^{-1}Ta)$ when $\epsilon=Ta^{-1/4}$. Here we used the fact that within the near-wall boundary layer, the size of $\Gamma$ is $O(1)$ from the boundary conditions, and the same size must be used in the plume.

The asymptotic structure can be summarised as follows. Within the core region we use the expansions

where $\omega_{0}$ and $\gamma_{0}$ are constants. The expansion is consistent with the Prandtl-Batchelor structure to leading order, and $\Psi_{0}$ must be determined by

in the aforementioned region $\mathcal{A}$ with the boundary condition $\Psi_{0}=0$ on $\mathcal{C}$.

Let $l$ and $n$ be the distance along $\mathcal{C}$ and its inward normal, respectively. Then using the rescaled normal variable $N=\epsilon^{-1}n$, the flow within the boundary layer can be expanded as

The leading order equations in the boundary layer are

where $\varphi$ is the angle between the $z$-axis and the $n$-axis.

When the boundary layer is in contact with the cylinder wall, the following boundary conditions must be imposed at $N=0$:

For the near-wall boundary layer, $\cos\varphi=0$ so (22) is none other than Prandtl’s boundary layer equation. While for large $N$, the flow must match the core solution. Let $U(l)$ be the value of $\Psi_{0n}$ on $\mathcal{C}$, which can be found by the core flow problem (19). Then as $N\rightarrow\infty$, the boundary layer solution must satisfy

When the boundary layer is detached from the wall, the similar far-field conditions should be applied on both sides $N\rightarrow\pm\infty$.

The theoretical boundary layer thickness $\epsilon=Ta^{-1/4}$ implies the Nusselt number scaling $Nu\propto Ta^{\beta}$ with $\beta=1/4$. We can check if this is consistent with the angular momentum transport. By averaging (11a) with respect to $z$ and further radially integrating from $r_{i}$ to $r$, the transport balance can be obtained as

Here we used overline to denote the average with respect to $z$ and wrote $\Gamma=\overline{\Gamma}(r)+\widetilde{\Gamma}(r,z)$. Let us consider this balance at a sufficient distance from both walls. In the core region, $\widetilde{\Gamma}$ would be exponentially small because of the Prandtl-Batchelor theorem, so the transport is extremely inefficient. The plume region of thickness $O(\epsilon)$ is hence the main contributor to the left-hand side, which can be estimated as $O(\epsilon(\partial_{z}\Psi|_{b})\widetilde{\Gamma}|_{b})=O(\epsilon^{-1})$ using the scaling $\widetilde{\Gamma}|_{b}=O(1)$, $\Psi|_{b}=O(\epsilon^{-1})$. This term indeed balances with the second term on the right-hand side.

The asymptotic structure derived here is consistent with the behaviour of the numerical Taylor vortex solution. Figure 6 summarises the numerical results for $\eta=0.5$, $a=-1/8$. The red solid curve in the figure shows the scaling $Nu\propto Ta^{1/4}$, which can also be seen in the data presented in figure 1 ($\eta=5/7,a=0$). The green dashed and the blue dotted curves are $Ta^{-1/2}rv$ and $Ta^{-1/2}\omega/r$ measured at the centre of the roll cell. According to the theory, they converge towards the constants $\gamma_{0}$ and $\omega_{0}$, respectively. The scaling of $\omega/r$ corresponds to the scaling of $u$, $w$ in the core, and hence the scaling of $Re_{w}$.

One may have noticed that the exponent $\beta=1/4$ of the Nusselt number differs from the exponent $\beta=1/3$ deduced in the asymptotic analysis of Chini & Cox (2009) and Hepworth (2014). This discrepancy is not due to the differences in the flow driving mechanism but to the boundary conditions at the walls. If the zero stress condition is imposed, the linear extrapolation of the core streamfunction already satisfies the boundary condition to leading order. This means that the balance $O(\Omega|_{b})=O(\epsilon^{-1}Re_{w})$ we assumed for the no-slip case is not necessary. For the slip wall case, the magnitude of the vorticity does not change in the core and the boundary layer, so the strong vortex layer seen in figure 4c does not appear. If we use the balance $O(\Omega|_{b})=O(\Omega|_{c})=O(\epsilon^{-2})$ instead for the scaling argument of the plume we have $\epsilon=Ta^{-1/3}$, as expected.

For the slip wall RBC, further analytical progress has been made using the fact that the boundary layer equations become linear and the roll cells are rectangular (Chini & Cox 2009; Hepworth 2014). In our case, however, we have to rely on numerical calculations because the boundary layer equations are fully nonlinear. Moreover, the shape of the core region is nontrivial due to the small vortices appearing near the corners (see figure 4c). This means that the core and boundary layer equations (4.8), (4.10), (4.11) need to be iteratively solved by updating the core shapes and constants $\gamma_{0}$ and $\omega_{0}$. Such numerical calculations are too challenging and out of the scope of this paper. Differences in the structure of the boundary layer also affect the $Pr$ dependence of the asymptotic solution of RBC. For the slip wall case, asymptotic solutions for different $Pr$ can be obtained by rescaling the $Pr=1$ solution. However, this is possible because the boundary layer is linear, which is not the case for the no-slip case.

Finally, we show that the above analysis provides some insights into the structure of the mean flow. As already seen in figure 4b, the angular momentum $rv$ is almost constant in the core region. Therefore, the mean angular momentum $r\overline{v}$ is expected to become a constant away from the wall. This is indeed the case for the numerical Taylor vortex solution (figure 7a). In the studies of TCF turbulence, on the other hand, the mean angular velocity $q$ is usually plotted, and it is often noted that the profile is linear in the core. At first glance, this appears to be the case, for example, when looking at figure 2a, but this is because the cylinder gap ($\eta\approx 0.714$) is too narrow to clearly see the radial dependence of the profile (see figure 7a for the $q$ profile for a wide gap case $\eta=0.5$). If the numerical results by Ostilla-Mónico et al. (2014b) are summarised in terms of $rv$, as shown in figure 7b, they clearly show the constant angular momentum property in the core. Note that when $a$ becomes too large, the Taylor vortex appears to favour the vicinity of the inner cylinder, and the homogenisation is only observed there. In the long history of TCF studies, some researchers have also pointed out that the turbulent mean flows might have a uniform angular momentum zone (Wattendorf 1935, Taylor 1935, Smith & Townsend 1982, Lewis & Swinney 1999, Dong 2007, Brauckmann & Eckhardt 2017). However, this fact is not widely known, probably because the mathematical reasons behind it has not been elucidated.

The aim of this section is to examine the asymptotic behaviour of the solutions at the first and second peaks seen in figure 3a. To this end, in figure 8, we performed similar calculations at two higher Taylor numbers. The results are summarised using the theoretical scaling to be derived in this section. Essentially, the scaling of $k$ and $Nu$ represent the width of the roll cell in the $z$ direction and the thickness of the boundary layer adjacent to the cylinder wall, respectively. The computations of the solutions are not easy, especially around the first peak, where very thin boundary layers need to be resolved. Although the convergence of the numerical solutions to the asymptotic states is still not perfect, the theoretical scalings are consistent with the overall features of the numerical data.

The structure of the flow field in both peaks can be roughly divided into a boundary layer near the wall and a core region in the middle of the gap, as we have seen in figures 2 and 3. On closer inspection, one further notices that the asymptotic structures of the two flows are quite different. For example, in the first peak solution, $\overline{\Gamma}$ has a flat profile in the core region (figure 9a), but this is not the case in the second peak solution (figure 9b). Figure 10 examines how the core flows develop from the near wall region adjacent to the inner cylinder. In the first peak solution, the near-wall structure is somewhat similar to the $k=O(1)$ case, with the wall boundary layer becoming a sharp plume as rounding the corner (top panel). On the other hand, in the second peak solution, no apparent plume can be recognised away from the wall, and the flow only varies slowly in the $z$-direction (bottom panel). The core flow inherits this property, as seen from figure 11, where the axial structure of the fluctuation fields is plotted at the midgap $r=r_{m}=(r_{o}+r_{i})/2$. For the second peak solution, only a single Fourier mode plays a major role in the core, while for the first peak, a number of harmonics participate in forming the plume. A theoretical explanation of those differences will be deduced below, together with the detailed scalings of the flows.

The goal of this section is to derive a matched asymptotic expansion consistent with the behaviour of the second peak solutions. This asymptotic state can be reached by decreasing the wavenumber from the linear critical point, as seen in figure 8b. However, for $k/Ta^{1/4}$ ranging from about 0.6 to 0.8, the appropriate scaling of the Nusselt number is $O(Ta^{0})$, which is different from that observed for the second peak state; see figure 14a. We shall study this transitional state in section 6.1, followed by the analysis of the second peak in section 6.2.

The most important previous work throughout this section is due to Denier (1992), who applied the Hall & Lakin (1988) type high wavenumber asymptotic theory to Taylor-Couette flow in the narrow gap limit. As briefly commented in that paper, the extension of the theory to wide-gap cases should be straightforward. In section 6.1, we derive an analytic expression of the Nusselt number and compare it with the numerical solution for the first time.

As the solution moves away from the linear critical point in the transitional state, the vortex grows from the vicinity of the inner cylinder (figure 14c). The asymptotic state corresponding to the second peak appears when the vortex fills the entire gap. A similar scenario was suggested in Denier (1992) using a matched asymptotic expansion, but the scaling obtained is unfortunately not consistent with the behaviour of the numerical solutions. The main reason for this is that the matching is not actually possible in the two-layer boundary layer structure assumed in Denier (1992). In section 6.2, we will resolve that difficulty by modifying the structure of the near wall-zone to three layers.

We choose $\delta=k^{-1}$ as a small perturbation parameter and introduce the scaled axial variable $Z=\delta^{-1}z$. According to Hall (1982), the linear critical point of curved flow problems typically has the Taylor number expansion $Ta=\delta^{-4}T_{0}+\cdots$. This is also true for Taylor Couette flow from the linear critical point to the second peak.