Linear stability analysis of the flow between rotating cylinders of wide gap

This study investigated by an analytical method the flow that develops in the gap between concentric rotatingcylinderswhentheTaylornumber Ta exceedsthefirstcriticalvalue.Concentriccylindersrotating at the speed ratio µ = 0 are investigated over the radius ratio range 0 . 20 ≤ η ≤ 0 . 95. This range includes configurations characterised by a larger annular gap width d than classical journal bearing test cases and by a Taylor number beyond the first critical Taylor number at which Taylor vortices develop. The analysis focuses on determining the parameters for the direct transition from axisymmetric Couette flow to wavy Taylor vortex flow. The results show a marked change in trend as the radius ratio η reduces below 0.49 and 0.63 for the azimuthal wave-numbers m = 2 and 3 respectively. The axial wavenumber increases so that the resulting wavy Taylor vortex flow is characterised by vortex structures elongated in the radial direction, with a meridional cross-section that is significantly elliptical. The linear stability analysis of the perturbation equations suggests this instability pattern is neutrally stable. Whereas a direct transition from axisymmetric Couette flow is not necessarily the only route for the onset of wavy Taylor vortex flow, the significant difference between the predicted pattern at large gap widths and classical wavy Taylor


Introduction
The flow in the gap between concentric independently rotating cylinders has attracted great attention over the years.This dates back to 1888 and 1890, when Mallock [1] and Couette [2] conducted independent experiments using concentric rotating cylinders.The work of Taylor [3] provides experimental and analytical explanations for the appearance and for the development of flow instabilities between rotating concentric cylinders.Taylor's experiment [3] showed that, when the angular velocity of the inner cylinder is increased above a certain threshold, the steady and axisymmetric Couette flow becomes unstable.This leads to a new secondary steady flow state, the axisymmetric Taylor vortex flow, which is characterised by the formation of regularly spaced vortices spanning the gap between the inner cylinder and the outer cylinder in the axial direction.From the time Taylor [3] first established the threshold at which the axisymmetric Couette flow becomes unstable, further investigations have been conducted both experimentally and analytically to determine the onset of this rotational instability in terms of the first critical Taylor number Ta c .
These include the research contributions by Chandrasekhar [4], Donnelly [5], Brewster et al. [6], Jeffreys [7], Taylor [3], Donnelly and Schwarz [8], and Di Prima [9].This first critical Taylor number Ta c is described by Deshmukh et al. [10] and by Batten et al. [11] as the value of the Taylor number Ta at which the first transition occurs from axisymmetric Couette flow (CF) to axisymmetric Taylor vortex flow (TVF).
The hydrodynamic instability and the transition between different regimes of the flow between a non-rotating (Ω 2 = 0) outer cylinder and a rotating (Ω ̸ = 0) inner cylinder, at the speed ratio µ = Ω 2 Ω −1 = 0, can be determined as functions of different non-dimensional parameters.These parameters are the ratio of the  As the Taylor number increases, the steady and axisymmetric Taylor vortices that develop in the annular region of concentric cylinders (Fig. 1a) change to unsteady, asymmetric and wavy Taylor vortices [3,[13][14][15] (Fig. 1b).These flow patterns are illustrated by the flow visualisation with aluminium powder from Koschmieder [12] reported in Fig. 1.Cole [16] and Zati et al. [17], as quoted in Di Prima et al. [18], and Coles [14] and Debler et al. [19] found by experiment that wavy Taylor vortex flow could be obtained by a transition from axisymmetric Couette flow as well as by a transition from axisymmetric Taylor vortex flow.It is of interest to consider the specific case in which a wavy Taylor vortex flow is obtained by a sudden increase in the Taylor number, starting from an axisymmetric Couette flow.
The wavy Taylor vortex flow is characterised by travelling azimuthal waves that are superimposed on the steady and axisymmetric Taylor vortices [20].The azimuthal waves rotate around the inner cylinder at some wave angular speed ω [21,22].The azimuthal waves have a defined integer azimuthal wave-number, m, and move at a finite wave velocity ωR o in the azimuthal direction [23].The wave velocity is observed to be approximately 30% to 50% of the surface speed of the inner cylinder, depending on the Taylor number and other conditions [14,15,[23][24][25].The azimuthal wave-number m, the axial wave-number a, and the Taylor number Ta all depend upon the value of η and on the aspect ratio Γ [25,26].
The asymmetric, unsteady and wavy Taylor vortices have been studied experimentally by Coles [14], Schwarz et al. [15], Nissan et al. [13], Snyder [27], Debler et al. [19], and Castle et al. [28].Di Prima [23] studied the asymmetric disturbances analytically by using the Galerkin method to solve the resulting eigen-value problem for the case of a small gap (η → 1) at µ = 0 and at µ = 0.5.The results of his investigation show that the critical speed, and hence the Ta at which a mode is cut on, increases with increasing azimuthal wave-number m.Krueger et al. [22] studied analytically the asymmetric disturbances using both a small-gap approximation and the full linear perturbation equations over the ranges 0.6 ≤ η ≤ 0.95 and −1 ≤ µ ≤ 0. Davey et al. [29] employed a double amplitude expansion for their analytical investigation.Eagles [30] extended the method employed by Davey et al. [29] using fifth-order amplitude expansions.He used the full linear disturbances equations with and without the smallgap approximation by Krueger et al. [22].Bisshopp [31] provided a numerical solution for the narrow gap approximation for the case where the cylinders rotate in opposite directions and the ratio of the azimuthal wave-number to the axial wave-number is sufficiently small.Roberts [32] used direct numerical calculations for µ = 0 over a gap width ratio range 0.75 ≤ η ≤ 0.95 and confirmed the results by Di Prima [23].Roberts [32] departed from the small gap approximation, documenting the gap width ratios η = 0.95, 0.90, 0.85 and 0.75.For various values of m at different η, Roberts [32] predicted the corresponding values of the critical axial wave-number a c , the critical angular wave velocity ω c , and the critical value of the Taylor number at which there is a transition from axisymmetric Couette flow to wavy Taylor vortex flow.This is herein referred to as the second critical Taylor number Ta c2 .His results showed that the Ta c2 for the onset of wavy Taylor vortices is slightly higher than Ta c for axisymmetric Taylor vortices.In all the cases reported by Roberts [32], an increment in the azimuthal wave-number m at constant radius ratio η is associated to a higher second critical Taylor number and critical axial wave-number but to a lower critical angular wave velocity ω c .As η decreases, the value of ω c also decreases.Eagles [30] reported that the Taylor number for the onset of wavy Taylor vortex flow was about 11% above the first critical Taylor number Ta c , for the results obtained without the small gap approximation.The result obtained by Eagles for η = 0.95 is in qualitative agreement with the results obtained by Davey et al. [29].
The analytical investigation of Roberts [32] for the stability of axisymmetric Couette flow to asymmetric disturbances does not provide results for the second critical Taylor number, the critical axial wave-number, and critical angular wave velocity for η < 0.75 at µ = 0.This information is only available for just a few combinations of η and µ from a handful of publications, like Di Prima et al. [18].
It is therefore of interest to further explore the flow regime in an idealised geometry where the radius ratio η ≤ 0.65.This work aims to explore a subset of the non-dimensional parameter space of the asymmetric vortex flow for the case of an infinitelength coaxial cylinder assembly in which the inner cylinder is rotating and the outer cylinder is fixed.This study wishes to describe the bulk flow parameters of the wavy Taylor vortex flow regime, including the azimuthal wave-number, the axial wavenumber, and the azimuthal wave speed, based on a new numerical solution of the analytical model by Roberts [32] for a low radius ratio configuration.
The main body of this paper is structured into seven sections.Section 2 summarises the analytical model of Roberts [32] for the linear stability analysis of axisymmetric Couette flow and proposes its application over the extended radius ratio range 0.28 ≤ η ≤ 0.95.Section 3 presents two new numerical procedures for solving the coupled ordinary differential equations in the analytical model by Roberts [32].In Section 4, the new numerical procedures are cross-validated with one another and are further validated by comparing their results against the corresponding values published in Roberts [32], in Schwarz et al. [15], and in Krueger et al. [22], over the radius ratio range 0.75 ≤ η ≤ 0.95, for asymmetric disturbances, and by Roberts [32], Krueger et al. [22], Di Prima et al. [18], and Donnelly and Schwarz [8] over the range 0.20 ≤ η ≤ 0.95 for axisymmetric disturbances.Section 5 presents the main novel contribution from this work, namely the characterisation of the asymmetric and wavy Taylor vortex flows that develop in the annular region of coaxial cylinders at wider gaps, over the range 0.28 ≤ η ≤ 0.75.This section emphasises the substantial qualitative and quantitative differences that the resulting wavy Taylor vortex flow pattern is predicted to have compared to the narrower gap width assemblies considered by Roberts [32] and similarly by other authors.Section 5 then compares the results reported in Section 4 with other numerical results from the literature for four different radius ratios.This shows that the methods of Section 3 can satisfactorily predict the non-dimensional parameters (Ta c2 , a c , ω c ) reported in the literature.It then provides new predictions for η = 0.50 and for η = 0.53, which was tested by the authors in past work [33].Section 6 defines the scope boundaries of the current work and clarifies its limitations with respect to practical applications.Conclusions from this research work are presented in Section 7, specifically, on what new knowledge has been exposed by this study.

Linear stability condition of the axisymmetric Couette flow
In this investigation, the stability equations of Roberts [32] are solved with two numerical schemes different than in the 1965 seminal work [32].The two numerical procedures provide results for values of azimuthal wave-number m greater than those presented by Roberts and by other investigators, for µ = 0, over the range η ≤ 0.75.
The linear stability analysis of the axisymmetric Couette flow used in the present work is detailed in section '(b) solutions to the non-symmetric problem' of the Appendix by Roberts [32] to a journal article by Donnelly and Schwarz [8].Roberts [32] considers the inviscid, incompressible, steady Couette flow developing between co-axial circular cylinders of radii R i and R o assuming that the flow is unrestricted in the axial dimension (L → ∞).The problem is stated in non-dimensional form, with R o normalising length and 2π /Ω normalising time.Taylor [3] provides a steady, spanwise uniform azimuthal velocity profile that satisfies the radial equilibrium and the solid wall conditions at µ = 0, which is referred to as the Couette flow in Roberts [32], and is where Ta = 1 2 η −2 (1 + η) (1 − η) 5 T ′ and η ≤ r ≤ 1. Roberts [32] examines the stability of this primary flow by solving the characteristic value problem of the perturbation equations, (A19) in Roberts [32], that are six first-order differential equations with complex coefficients In Eqs. ( 3)-( 8), the velocity perturbations u, v, w in the radial, azimuthal, and axial directions are assumed to be harmonic and to vary as exp i[m The non-dimensional axial wavenumber a ′ = aR o , in which a is the dimensional axial wavenumber in metres, the azimuthal wave-number m is integer as the domain is 2π periodic, the azimuthal wave non-dimensional , and X , Y , Z are auxiliary variables used to unroll the problem to coupled firstorder differential equations.In particular, X is related to the perturbation in pressure through the continuity equation.A solution to Eqs. ( 3)-( 8) is sought under the wall boundary condition u = v = w = 0 at r = η and at r = 1.This is achieved by generating a base of three independent solutions, denoted by subscripts 1, 2, 3, which satisfy the wall boundary condition at r = η, and seeking a linear combination of these three solutions that also satisfies the wall boundary condition at r = 1.The base of three independent solutions is generated by setting at r = η the starting conditions The parabolic set of Eqs. ( 3)-( 8) is then integrated in r numerically up to r = 1, where the complex valued determinant is evaluated in order to locate the complex zeros of ∆.These zeros correspond to linear combinations of the base solutions that satisfy the wall boundary condition at r = 1.

Numerical methods for solving the linear stability system
To perform the numerical integration of Eqs. ( 3)-( 8), two different numerical methods, hereafter referred to as method 1 and method 2, were used.Using two methods for solving the same linear stability system provided cross-validated predictions.This is important for the wide gap configurations for which comparative data from the literature is not yet available.
In method 1, the continuous derivatives of Eqs. ( 3)-( 8) were approximated by first-order forward differences on a staggered grid.This staggered grid discretisation differs from the uniform grid discretisation used in Roberts and it is a novel element of this work.In the staggered grid scheme, the velocity perturbations u, v, w are defined on nodes located at i∆r, while the auxiliary variables X , Y , Z are defined on the staggered grid with nodes located at (i + 1 2 )∆r, 1 ≤ i ≤ 1 + (1 − η)/∆r, with i integer.The staggered grid enables the use of stencil-centred values for u, v, w in the computation of the derivatives for X , Y , Z and vice-versa, eliminating the requirement for a local interpolation.A half-step width integration is used to prime the first internal nodes on the staggered grid from the boundary values of X , Y , Z .A sensitivity analysis was performed on the dependence of the value of the determinant on ∆r that led to selecting a step length of ∆r = 0.00015625, that is lower than ∆r = 0.004 used in Roberts [32].
In Eqs. ( 3)-( 8), T ′ , ω ′ , a ′ , and m appear as free parameters.For a given azimuthal wave-number m, locating the complex zeros involves a search in (T ′ , ω ′ , a ′ ), starting from initial guesses for (T ′ , ω ′ , a ′ ), for which the values from Table A4 in Roberts [32] are used.The way this search is conducted is novel for this problem.It involves searching for the minimum of the modulus of the determinant |∆| using the non-linear Nelder-Mead method as described by Nelder and Mead [34], Powell [35], and Rosenbrock [36].The Matlab implementation of the Nelder-Mead algorithm in double precision is used to obtain the zeros to machine precision.The complex valued determinant provides only two constraints, through its real and imaginary values, for the three free parameters (T ′ , ω ′ , a ′ ).
A sweep in a ′ is performed in the neighbourhood of a ′ from Table A4 in Roberts [32], with a resolution ∆a ′ = 0.001R o , solving for (T ′ , ω ′ ) by the Nelder-Mead method.This sweep identifies the minimum in T ′ , which is Ta ′ c and is the lowest Taylor number of azimuthal wave-number m.This is the critical Taylor number at the azimuthal wave-number m with its associated non-dimensional critical angular velocity ω c ′ and non-dimensional critical axial wave-number a c ′ , defined by the curve in (T ′ , ω ′ , a ′ ).
In method 2, the continuous derivatives of Eqs. ( 3)-( 8) were approximated by the improved Euler method [37] on a uniform mesh ∆r = 0.0000488.This finite differences central scheme is second-order accurate.By this, the complex valued determinant ∆ of Eq. ( 10) is evaluated explicitly, given T ′ , ω ′ , a ′ , and m.For a given integer-valued azimuthal wave-number m, a search in (T ′ , ω ′ , a ′ ) is performed using the non-linear Nelder-Mead method for the minima of |∆| and then the lowest T ′ among the search outputs is selected.This gives a second independent estimate of Ta ′ c , ω c ′ , and a c ′ , for a given m.

Numerical predictions of asymmetric instability modes for narrow gaps d
The (Ta ′ c , ω ′ c , a ′ c ) results obtained from method 1 and method 2 match one another to within five significant figures over the range 1 ≤ m ≤ 3 and 0.28 < η ≤ 0.95, at µ = 0.This crossvalidation provides good confidence in these results.The new numerical solutions of the perturbation equations are compared in Fig. 2 and in Fig. 3 with the results reported by Roberts [32], by Schwarz et al. [15], and by Kruger et al. [22], for different values of η over the range 0.75 ≤ η ≤ 0.95, at µ = 0. Fig. 2 shows the second critical Taylor number, Ta c2 , against the radius ratio of the cylinders, η.Fig. 3 shows the critical axial wave-number, a c , against the radius ratio of the cylinders, η.The critical axial wave-number a c is normalised by the cylinder gap width d, so that a c In the ordinate scale shown, the results from methods 1 and 2 overlap.The points corresponding to the same m value are outlined in Fig. 2 and in Fig. 3 by a unique line pattern.In Fig. 2 and in Fig. 3, the dashed line designates the Ta c2 and a c values at m = 1, the dashed dot line designates the Ta c2 and a c values at m = 2, while the dotted line designates the Ta c2 and a c values at m = 3, over the range 0.75 ≤ η ≤ 0.95.
The current numerical methods provide Ta c2 and a c values for m ≥ 2 in Figs. 2 and 3, extending the range of values reported both by Roberts [32] and by Krueger et al. [22] at η = 0.75.The analytical results from Krueger et al. [22] were based on a small gap approximation, while the results from the present work and from Roberts [32] are based on the full linear perturbation equations.This produces some scatter among the diamond symbols at η = 0.95 in Fig. 2 and in Fig. 3.This difference notwithstanding, Figs.  and 2 (full symbols), Roberts [32], Krueger et al. [22], and experiment by Schwarz et al. [15] (open symbols).Fig. 3. Variation of a c for µ = 0 over the range 0.75 ≤ η ≤ 0.95.Methods 1 and 2 (full symbols), Roberts [32] and Krueger et al. [22] (open symbols).and 3 show a good agreement between the results from methods 1 and 2 and the ones from Roberts [32] and from Krueger et al. [22].
The new Ta c2 and a c values for m = 2 and m = 3 follow a trend similar to the one for the corresponding values at m = 1 from Roberts [32], over the range 0.75 ≤ η ≤ 0.95.This gives further confidence in the current numerical methods and forms the basis for testing them at η < 0.75.Fig. 4 shows the change in the normalised critical angular velocity for the first three non-zero wave-numbers m = 1, 2, 3 over the same gap width ratio range 0.75 ≤ η ≤ 0.95 as for Figs. 2 and  3.In Fig. 4, ω c is normalised by Ω.The predictions from methods 1 and 2 are compared against the corresponding published results by Roberts [32].At small gap widths, the normalised angular velocity assumes a value close to 0.5 of the rotational speed of the inner cylinder, indicating that the asymmetric disturbances rotate at about the average of the rotational speed between the inner and the outer cylinders.This is consistent with the representation of the wavy Taylor vortices as perturbations to an otherwise axially uniform Couette flow that, at a sufficiently small gap, would feature a constant angular velocity gradient from the inner cylinder to the outer cylinder.Therefore, the bulk angular velocity towards the (full symbols), Roberts [32] (open symbols).middle of the gap, where the centres of the wavy Taylor vortices are likely to be, is expected to be approximately 0.5 ΩR 1 .As the cylinder gap increases and η decreases, boundary layers form over the surfaces of the cylinder [33] and the angular velocity of the bulk flow at the middle of the gap reduces.Accordingly, Fig. 4 shows a monotonic reduction in ω c at increasing gap widths.There is a good overlap between ω c predicted by methods 1 and 2, which are shown by overlapping filled symbols, and the corresponding results from the literature [32].The new predictions for ω c at η = 0.75 show a consistent trend to that of the remainder of the dataset, that is, a monotonic decrease of ω c at decreasing η.
The good agreement among the predicted values of Ta c2 , a c , and ω c for asymmetric disturbances over the radius ratio range 0.75 ≤ η ≤ 0.95 indicates that the new numerical procedures of methods 1 and 2 for determining the second critical Taylor number for the first few modes provide satisfactory and consistent results with past published work.
It is of interest to extend the validation of the new numerical methods for determining the stability parameters of the axisymmetric Couette [2] flow under an axisymmetric flow constraint, for which published results are available over a wider radius ratio range.This validation is presented in the next section.

Numerical predictions of symmetric instability modes for wider gaps d
Figs. 5 and 6 report similar results to the ones presented in Fig. 2 and in Fig. 3 on extended axes, for the azimuthal wave-number m = 0 at µ = 0.These results cover the wider radius ratio range 0.20 ≤ η ≤ 0.95 and are compared against published results by Roberts [32], Krueger et al. [22], Di Prima et al. [18], and Donnelly and Schwarz [8].Fig. 5 shows the first critical Taylor number, Ta c , against the radius ratio of the cylinders, η, while Fig. 6 shows the corresponding critical axial wave-number, a c , against the radius ratio of the cylinders, η.The critical axial wave-number a c is normalised as in Fig. 3 [22,32].The overlap between the predictions by Roberts [32] (open symbols), Di Prima et al. [18] (open symbols), and the ones from methods 1 and 2 (filled symbols) shows that there is a good agreement between the m = 0 results over the range 0.20 ≤ η ≤ 0.95.This further confirms the validity of the present numerical methods as unified schemes for axisymmetric and asymmetric Taylor instability modes, at µ = 0. Figs. 5 and 6 have the same trend as the ones in Fig. 2 and in Fig. 3, respectively, in that Ta c and a c increase with decreasing η.This indicates that the Methods 1 and 2 (full symbols), Roberts [32], Krueger et al. [22], Di Prima et al. [18], and experimental data from Donnelly and Schwarz [8] (open symbols).transition from the axisymmetric Couette flow to the axisymmetric Taylor vortex flow occurs at higher Taylor numbers in cylinders of progressively larger gap configurations.Furthermore, the axial dimension of the resulting Taylor instabilities (axisymmetric Taylor vortices) progressively reduces compared to the cylinder gap d.
Comparing the results of Figs. 2 and 5, as well as of Figs. 3  and 6, it is observed that the second critical Taylor number and the critical axial wave-number of the asymmetric perturbations are above the first critical Taylor number and the critical axial wave-number of the axisymmetric perturbations.Furthermore, the increments are positively correlated to the azimuthal wavenumber m, in agreement with Koschmieder [12].
The new numerical methods for determining the Taylor vortex flow characteristics upon departure from the axisymmetric Couette flow, by the model of Roberts [32], give results consistent with one another and with past published work.This is for both axisymmetric and for asymmetric perturbations.This motivates the application of these newly validated numerical tools to exploring the departure from the axisymmetric Couette flow by a symmetric perturbations, where the cylinder gap width d is wider.

Numerical predictions of asymmetric modes at larger gaps d
Figs. 7 and 8 show the variation of Ta c2 and of a c with η at µ = 0 over the ranges 1 ≤ m ≤ 3 and 0.28 ≤ η ≤ 0.95, using line patterns consistent to the ones of Figs. 2, 3, and 4. Fig. 7 reports the same trend of Ta c2 increasing with decreasing η, which denotes configurations with a progressively larger gap d between the cylinders.The trends of Ta c2 versus η in Fig. 7 were found to be adequately described by the exponential function + c exp (dη), for m = 1, 2, 3, with (a, b, c, d) determined by a least squares regression to a residual value R ≥ 0.998.This multivariate regression was performed using the MAT-LAB open curve fitting tool ''cftool''.Fig. 7 gives a graphical view of the quality of this curve fit by showing the three curve fits for m = 1, 2, 3 by the lines and the values for Ta c2 sampled at different η by methods 1 and 2, by the symbols.Each symbol type clusters appreciably close to its corresponding curve fit, of which Table 1 reports the regressed coefficients.
At small gaps, over the range 0.75 ≤ η ≤ 0.95, Ta c2 and η are inversely related so that a steady and monotonic increase in Ta c2 is predicted at decreasing values of η.As the gap width becomes Methods 1 and 2 (full symbols), Roberts [32], Krueger et al. [22], and Di Prima et al. [18] (open symbols).wave-number a c with the radius ratio η for the azimuthal modes m = 1 to 3, over the same η range as Fig. 7.The trends of a c versus η in Fig. 8 were found to be adequately described by third-order polynomials of the form a c = aη + bη 2 + cη 3 + d, over two were determined by the same least squares regression technique as for the curve fits in Fig. 7, to residual values R ≥ 0.999.Table 2 reports the regressed coefficients obtained by ''cftool'' in MATLAB.
At the higher radius ratios, over the range 0.75 ≤ η ≤ 0.95, with a narrower gap d between the cylinders, Fig. 8 shows a decrease in a c with increasing radius ratio η, with a c from all three modes m = 1, 2, 3 converging towards π as η → 0.95.This a c value corresponds to toroidal Taylor vortices of circular cross-section in the meridional plane of the cylinder assembly, arranged in pairs of contra-rotating vortices, which repeat along the length of the cylinder assembly with an axial periodicity of approximately 2d.
This gives an axial wave-number for the vortex pair of 2π /2d and a non-dimensional axial wave-number a c ∼ π.At η < 0.95, there is a monotonic increase in the non-dimensional axial wavenumber a c with increasing azimuthal wave-number m.This can be explained by the radial displacement of adjacent vortices at m > 0 which, by having different azimuthal phase angles, can pack more tightly in the axial direction.At the lower radius ratios 0.41 ≤ η ≤ 0.75, with a wider gap d between the cylinders, the nondimensional critical axial wave-number is predicted to increase significantly and non-linearly with decreasing radius ratio η and increasing azimuthal wave-number m.There is a noticeable change in the gradient da c /dη for m = 3 at the radius ratio η ≃ 0.63, below which a steep rise in non-dimensional critical axial wave-number a c is predicted.This points to a new behaviour for the m = 3 mode, for which a c becomes closer to 2π than to π.This may indicate either significant vortex packing or the presence of Taylor vortices of elliptical cross-section in the meridional plane.This change in behaviour motivated the use of two different polynomial curve fits either side of η = 0.63, for which two set of regressed coefficients are reported in Table 2. Similarly, the azimuthal wave-number m = 2 is predicted to exhibit a rapid increase in the axial wave-number a c at radius ratios η ≤ 0.49, for which a c also increases to values ≥ 2π .As for m = 3, the predicted a c versus η trend of Fig. 8 for m = 2 indicates a change in behaviour at radius ratios η ≤ 0.49, which appears not to be documented in the literature.The higher values of a c over the range 0.41 ≤ η ≤ 0.49 are likewise compatible with either significant vortex packing, enabled by the variation in the Taylor vortex centre position in the azimuthal direction, or with Taylor vortices of elliptical cross-section, with the ellipse major axis in the radial direction, or possibly with a non-uniform radial displacement of the vortex cores.The lowest azimuthal wave-number m = 1 for non-symmetric Taylor instabilities does not exhibit evidence of any change of behaviour over the radius ratio range 0.28 ≤ η ≤ 0.95 like the one for the higher azimuthal wave-numbers m = 2 and m = 3.
Such difference could reflect a substantially axially constant phase angle of the m = 1 wave-number, similar to the one displayed by Fig. 1(b).This would prevent neighbouring Taylor vortices from stacking azimuthally, by having different phase angles, which would result in a higher axial wave-number behaviour possibly similar to that predicted for m = 2 and m = 3.However, Fig. 1(b) suggests that a substantially constant phase angle can also occur at higher azimuthal wave-numbers, since the number azimuthally periodic undulations suggests m > 3 and possibly m = 7 in the flow visualisation of Fig. 1(b).Alternatively, m = 1 may exhibit a similar change of behaviour for η < 0.28, which is outside the range tested by the authors.The possibility of a substantial difference in behaviour between m = 1 and m = 2 and 3 remains an open matter for further investigation.
In Fig. 9, the variation of the non-dimensional critical angular wave speed ω c with the radius ratio η is presented for the same three azimuthal wave-numbers, m = 1, 2, 3, using the same third degree polynomial function that was used to describe the trends in Fig. 8.The regressed coefficients are presented in Table 3 for different ranges of the cylinder radius ratio η.The multi-variable least squares regression, performed using ''cftool'' in MATLAB, gave residual values R ≥ 0.999 for the three curves in Fig. 9.The good overlap between the curves and symbols in Fig. 9 visually confirms a good quality fit to the ω c predictions at discrete values of η from methods 1 and 2 of Section 3.
Over the higher radius ratio range 0.75 ≤ η ≤ 0.95, ω c increases almost linearly with η, assuming slightly higher values as m increases from m = 1 to 3. As the radius ratio η increases, ω c increases monotonically toward ω c = 0.5 at η = 0.95, for all m.This is consistent with ω c = 0.5 obtained from the small gap approximation analysis by Krueger et al. [22].It may also be explained intuitively, as in Roberts [32], that the angular speed of the pattern may be estimated by the average of the rotating speeds of the bounding cylinders.The progressive reduction in critical angular wave speed at decreasing η may indicate a displacement of the centres of the wavy Taylor vortices towards the outer fixed cylinder, where the angular velocity of the local flow is lower than 0.5 of the inner cylinder angular speed, due to the outer cylinder no-slip wall effect.
Over the lower radius ratio range 0.28 ≤ η ≤ 0.75, with a wider gap d between the cylinders, Fig. 9 shows a progressive non-linear decrease in ω c with reducing η for the azimuthal wavenumber m = 1, essentially following the same trend as over the narrower gap ratio range 0.75 ≤ η ≤ 0.95 but for a second-order non-linear type effect.The trend for the azimuthal wave-numbers m = 2 and m = 3 is instead quite different.Following an initial decrement in ω c with reducing η up to η = 0.63 for m = 3 and up to η = 0.49 for m = 2, there is a trend reversal with the critical angular speed ω c increasing rapidly with reducing η.This trend seems to be peculiar to this wider gap range and is a novel result from the present numerical investigation.
This result has practical implications on the design of any instrumentation aimed at monitoring the angular wave speed either in a m = 2 or in a m = 3 asymmetric Taylor vortex flow at low tip radius ratios.This configuration is predicted to onset from an axisymmetric Couette flow above a relatively high second critical Taylor number Ta c2 , from Fig. 7, and be characterised by a high nondimensional critical angular wave speed ω c , from Fig. 9, therefore the sampling rate for resolving azimuthally this wavy Taylor vortex flow requires careful planning and the use of hardware capable of a suitably high acquisition rate.
In this work, the current numerical model is not tested for azimuthal wave-numbers higher than 3.
Table 4  Comparison with experiment is limited to modes m = 0 and m = 1 by the available experimental data.The lack of experimental data for comparison at small radius ratios prompted the authors not to evaluate radius ratios η < 0.28 in the current investigation, albeit methods 1 and 2 of Section 3 enable such computation to be performed by the interested reader.

Scope and limitations of the results
The results presented in this paper are limited to predicting the transition from axisymmetric Couette flow to either axisymmetric Taylor vortex flow or to wavy Taylor vortex flow in concentric cylinders of infinite length.In experiment, the finite length of the cylinders is likely to affect the model predictions, more so at the wider gap d that this paper specifically addresses.This is because, for the same cylinder length L, the normalised length L/d that characterises the flow aspect ratio is lower in a wide gap configuration 0.28 ≤ η < 0.75 than for the 0.75 ≤ η ≤ 0.95 geometries of narrower gap studied by Roberts [32].The analytical model of Section 2 considers the direct transition from axisymmetric Couette flow to wavy Taylor vortex flow of azimuthal wave-number m ∈ [1,3].An alternative mechanism for the onset of wavy Taylor vortex flows of azimuthal wavenumber m > 0 is by staging through the axisymmetric Taylor vortex flow regime, as discussed in Di Prima et al. [18], and/or, more generally, through a combination of modes of azimuthal wave-numbers less than m.Applications in which the inner cylinder target rotational speed is reached through a modest acceleration are likely to be more prone to exhibit such staging behaviour.Therefore, the current predictions are likely to be limited to applications with a high spool up rate of the inner cylinder.In this context, Table 4 indicates a steady increase in the second critical Taylor number with increasing m.This makes more difficult implementing a rapid rotational acceleration to Ta > Ta c2 for the higher m numbers, in a practical engineering application, while maintaining a laminar flow.This motivated limiting the current investigation to m ≤ 3.
Section 2 identifies the conditions for neutral stability of the axisymmetric Couette flow system.A cylindrical assembly rotating at Ta > Ta c2 for m = 3 also meets the neutral stability criteria for m < 3. The analysis therefore provides only the sufficient condition for a given wavy Taylor vortex mode to be cut on.It does not predict the mode energy content, nor does it details the interaction among the cut-on modes, namely, whether the cut-on modes are mutually supportive or competitive and, if so, whether there is dominance by any given mode.
The regressions for the analytical estimation of the non-dimensional parameters (Ta c2 , a c , ω c ) for different azimuthal wave-numbers m are simple curve fits.The specific expressions were determined by testing several options among commonly used algebraic functions.The selection of the specific expressions is not meant to reflect either a simplified model or an interpretation of the flow physics.Therefore, demonstrating any appropriateness of the selected algebraic functions beyond the fact that they provide numerically good residual values of R ≥ 0.99, when fitted to the data, over the stated η range, is outside the scope of this paper.

Conclusions
The majority of the past investigations on the Taylor vortex flow focussed on high radius ratio (small gap width) configurations, relevant to journal bearing flows.The emphasis of this paper is instead on low radius ratio (wide gap width) configurations.These configurations are relevant to engineering applications of high societal and economic impact, such as the annular flow in oil wells, the flow in submerged pumps for water wells, and the flow in vertical turbine pumps.
The idealised geometry of two concentric cylinders of infinite length was modelled analytically as in Roberts [32].Numerical solutions of the perturbation equations for the axisymmetric Couette flow enabled the estimation of the sufficient conditions for the onset of axisymmetric and of wavy Taylor vortex flows, over a wider range of radius ratios 0.28 ≤ η ≤ 0.95 than in Roberts [32], for integer azimuthal wave-numbers m ∈ [0, 3], at the speed ratio µ = 0. Two iterative numerical procedures were used to obtain such solutions.Each procedure provides, by itself, a unified solution method for both axisymmetric and wavy Taylor instability modes.
Simpler analytical expressions for the non-dimensional parameters (Ta c2 , a c , ω c ) are also provided, based on curve fitting, that give a first approximation that is non-iterative and easily computable.These results give the first approximation to use in the two iterative procedures of Section 3.
The linear stability analysis using the small perturbation model by Roberts [32] predicted a change in behaviour in the wider gap configurations 0.41 ≤ η ≤ 0.75 for the azimuthal wave-numbers m = 2 and 3.The variation of the critical axial wave-number a c and of the critical azimuthal wave speed ω c with η shows a sudden change in the a c and ω c gradients at the radius ratios η = 0.49, for m = 2, and η = 0.63, for m = 3.This new behaviour leads to m = 2 and 3 azimuthal modes characterised by tightly packed vortices (a c ∼ 2π ) in the axial direction and by a higher azimuthal wave speed than over the larger cylinder radius ratio range 0.75 ≤ η ≤ 0.95.This new behaviour was not observed at any azimuthal wavenumber m < 2. It appears to be unreported in the literature, possibly due to the preponderance of work at smaller gap ratios, using radius ratios in the range 0.75 ≤ η ≤ 0.95.The new behaviour may have relevance to the pattern of wear in the outer pipe of wells, due to abrasion by suspended particulate, as well as to the input energy required for pumping, due to the mixing loss produced by the wavy Taylor vortex flows.
The current linear stability analysis provides limited insight into the modal energy distribution and into the mode dominance for configurations in which more than one azimuthal mode is cut on.It is also possible that the flow arrives at the higher azimuthal wave-number modes by a process of staging through lower m values, which the current model does not account for, although this behaviour is reasonably less likely with a fast acceleration to the target angular speed Ω of the inner cylinder.In view of the aforementioned significant economic and societal relevance of these flows, it would be of interest testing for this new behaviour experimentally at different gap widths.

2η 2 d 4 1 − η 2
and the Reynolds number Re = ΩR i d/ν or the Taylor number Ta =
lists the second critical Taylor number Ta c2 , the critical axial wave-number a c , and the critical angular wave velocity ω c characterising the transition from axisymmetric Couette flow to wavy Taylor vortex flow of azimuthal wave-number m ≤ 3, for discrete values of η.Predictions of Ta c2 , a c , and ω c are listed alongside published reference data [15,22,32].Table 4 and Figs. 5 and 7 show a monotonic increase in Ta c and in Ta c2 with decreasing η.The difference Ta c2 -Ta c also increases with decreasing η, indicating a wider separation in Taylor number of the transition from axisymmetric Couette flow to Taylor vortex flows of different azimuthal wave numbers.

Table 1
Regressed coefficients for the estimation of Ta c2 over the gap width range 0.28 ≤ η ≤ 0.95.

Table 2
Regressed coefficients for the estimation of a c over the gap width range 0.28 ≤ η ≤ 0.95.

Table 3
Regressed coefficients for the estimation of ω c over the gap width range 0.28 ≤ η ≤ 0.95.