Modal stability analysis of toroidal pipe flow approaching zero curvature

Abstract The present study investigates the modal stability of the steady incompressible flow inside a toroidal pipe for values of the curvature $\delta$ (ratio between pipe and torus radii) approaching zero, i.e. the limit of a straight pipe. The global neutral stability curve for $10^{-7} \leq \delta \leq ~10^{-2}$ is traced using a continuation algorithm. Two different families of unstable eigenmodes are identified. For curvatures below $1.5 \times 10^{-6}$, the critical Reynolds number ${{Re}}_{cr}$ is proportional to $\delta ^{-1/2}$. Hence, the critical Dean number is constant, ${{De}}_{cr} = 2\,{{Re}}_{cr}\,\sqrt {\delta } \approx 113$. This behaviour confirms that the Hagen–Poiseuille flow is stable to infinitesimal perturbations for any Reynolds number and suggests that a continuous transition from the curved to the straight pipe takes place as far as it regards the stability properties. For low values of the curvature, an approximate self-similar solution for the steady base flow can be obtained at a fixed Dean number. Exploiting the proposed semi-analytic scaling in the stability analysis provides satisfactory results.


Introduction
The computations by Meseguer & Trefethen (2003) showed that the Hagen-Poiseuille flow is stable to infinitesimal disturbances even for Re = 5 × 10 6 , with Re being the Reynolds number based on the bulk velocity U b , the radius of the pipe cross-section R p , and the kinematic viscosity ν (a different definition is used in the reference).Extrapolating their numerical results, the authors inferred that the laminar flow in a straight pipe is modally stable for Re → ∞.This conclusion stems from numerical evidence.Analytical proof supporting it exists only for axisymmetric perturbations (Herron 1991), while it is missing in the case of non-axisymmetric disturbances.Nonetheless, transition to turbulence in a straight pipe has been known to occur at much lower values of Re since the seminal work by Reynolds (1883), even though a deeper physical understanding of the process has been achieved only recently (Barkley 2016;Avila, Barkley & Hof 2023).
For curved pipes, the flow becomes turbulent at higher Reynolds numbers than for a straight geometry (Taylor 1929;White 1929), and different transition mechanisms take place at the inner and outer walls, respectively (Sreenivasan & Strykowski 1983).For values of the curvature -defined as the ratio between the radius of the cross-section of the pipe and the curvature radius of the torus centreline (δ = R p /R c ) -below approximately 0.028, the flow undergoes subcritical transition (Kühnen et al. 2015), and the intermittent coexistence of laminar flow and turbulent puffs can be observed (Rinaldi, Canton & Schlatter 2019).On the other hand, for higher curvatures, the transition is initiated by a supercritical Hopf bifurcation, as shown numerically by Canton, Schlatter & Örlü (2016) and experimentally by Kühnen et al. (2014), and a bifurcation cascade takes place (Canton et al. 2020).In a restricted region of the (δ, Re) parameter space, both scenarios are possible, and the asymptotic state is dictated by the initial condition (Canton et al. 2020).
The modal stability analysis of straight pipe flow, or flows in toroidal pipes with curvature approaching zero, is not directly relevant in practice, since the transition to turbulence is triggered by finite-amplitude perturbations (Kühnen et al. 2015).Nevertheless, the study of the global neutral stability curve for these curvatures can provide further theoretical evidence of the stability of the Hagen-Poiseuille flow to infinitesimal perturbations at any Reynolds number.This conjecture is based on the fact that the flow in a torus tends to that in a straight pipe for δ → 0, and the same behaviour is therefore expected for its stability properties.To verify this hypothesis here, the global neutral stability curve computed by Canton et al. (2016) is extended to values of the curvature between 10 −7 and 10 −2 .For this purpose, we employ a tailored continuation algorithm and accurate numerical computations of both base flow and disturbances.

Base flow
The dynamics of the incompressible flow of a viscous Newtonian fluid in a toroidal pipe is described by the incompressible Navier-Stokes equations, which are made dimensionless by scaling with the bulk velocity U b , the radius of the cross-section of the pipe R p , and the constant density of the fluid ρ, and read They are formulated using an orthogonal toroidal reference system {s, r, θ} (Germano 1982), shown in figure 1 together with a sketch of the toroidal pipe.Therefore, u = (u s , u r , u θ ) T is the velocity vector, p is the non-dimensional pressure, and Re = U b R p /ν.The flow is driven by a volume force f = (F/h s )e s , where e s is the unit vector in the streamwise direction, F is a scalar constant representing the forcing amplitude, and h s = 1 + δr sin(θ ).This force mimics the effect of a streamwise pressure gradient, as described by Noorani & Schlatter (2015) and Canton, Örlü & Schlatter (2017).
Provided that the streamwise pressure gradient is set to zero since the flow is driven by the volume force defined above, the incompressible Navier-Stokes equations (2.1) written in toroidal coordinates read as follows.For the s-momentum, (2.2a) For the r-momentum, For the θ -momentum, (2.2c) The incompressibility constraint is It needs to be remarked that for δ = 0, (2.2) can be written as the equations for a straight pipe flow (Meseguer & Trefethen 2003), confirming that the flow in a torus tends to the Hagen-Poiseuille flow as δ → 0.
The steady base flow Q = (U, P) T satisfies the steady counterpart of the system (2.1) subject to no-slip and impermeability boundary conditions on the pipe wall and homogeneity in the streamwise direction s.Therefore, the problem can be solved on a two-dimensional cross-section of the toroidal pipe while retaining all three velocity components.

Modal stability analysis
The modal stability analysis is performed considering infinitesimal perturbations q = (u , p ) T evolving on top of the steady base flow Q.To this aim, the incompressible Navier-Stokes equations are linearised about the steady basic state, reading The formulation of (2.3) in toroidal coordinates is reported in Appendix A. At the wall, the same boundary conditions as for the nonlinear problem are imposed.Since the base flow is homogeneous in the streamwise direction, the following normal mode ansatz can be introduced for the perturbations where α = 2πR p /l s ∈ R is the streamwise wavenumber, with l s indicating the wavelength, and λ = ω + iσ ∈ C, with σ being the growth rate and ω the angular frequency.In this framework, the modal stability analysis is often termed BiGlobal in the literature (Theofilis 2003(Theofilis , 2011)).
Substituting (2.4) into (2.3), the following generalised eigenvalue problem for the complex-valued eigenmode qα = ( ûα , pα ) T and the corresponding eigenvalue λ is obtained for each streamwise wavenumber α: where R is a singular operator (identity for the velocity, zero for the pressure), and L α corresponds to the Fourier transform of the linearised Navier-Stokes operator, i.e.
It is clear from the ansatz (2.4) that infinitesimal perturbations grow exponentially in time if Im(λ) = σ > 0, thus the flow exhibits a modal instability.Note that the values of the streamwise wavenumber α are discrete, such that k = α/δ ∈ Z, because of the geometrical periodicity of the torus.

Numerical method
3.1.Spatial discretisation The numerical computations are carried out on a two-dimensional domain using an in-house developed MATLAB code based on a Fourier-Chebyshev spectral collocation method and written in primitive variables.Fourier basis functions are used in the azimuthal direction θ , whereas Chebyshev polynomials are employed in the radial direction r.The latter are defined for r ∈ [−R p , R p ].In this way, the nodes are clustered only close to the wall.Moreover, by using an even number of Chebyshev polynomials, i.e. the highest order being odd, there are no grid points at the pipe centreline such that any singularity is avoided.The computational domain, sketched in figure 2, consists of N p = n r × n θ collocation points, where n r corresponds to half the nodes defined on −R p ≤ r ≤ R p , and n θ represents the number of grid points in the azimuthal direction.
The one-dimensional differentiation matrices are obtained using the routines in the MATLAB DMSUITE by Weideman & Reddy (2000).The two-dimensional differentiation matrices for the first-and second-order derivatives in the azimuthal direction can then be defined using the Kronecker product ⊗ as where D (1) θ,1D and D (2) θ,1D are the one-dimensional differentiation matrices in θ , and I n stands for the identity matrix of order n.Because of the node ordering and the sign change of the derivatives across r = 0, the two-dimensional differentiation matrices in the radial direction have a more complex structure.They have different expressions for the streamwise (subscript s) and the radial/azimuthal (subscript r/θ) velocity components and are denoted by In this example, n r = 3 and n θ = 6.The node ordering is also indicated. where k is a matrix incorporating the sign changes for even and odd derivatives for each component, given by (3.5)

Steady-state solver
The numerical computation of the steady base flow consists of finding the zeros of the nonlinear function N (q) representing the steady counterpart of (2.1), where q = (u, p) T .This task is performed by employing the Newton-Raphson method, which at the iteration n + 1 reads where J | q n is the Fréchet derivative of N (q) computed in q n .The initial guess for the Newton-Raphson method is provided by either the solution of the linear Stokes equations or a previous steady state computed for different parameters δ and Re.The solution is required to have a volumetric flow rate This condition is imposed by augmenting the system (3.6) with the constraint Q s ( f ) − πR 2 p U b = 0, which provides the required forcing amplitude F. Because a collocated grid is used, the solution exhibits spurious pressure modes (Canuto et al. 1988).The only spurious pressure mode affecting the problem under investigation is the checkerboard mode.In addition, the pressure is determined up to an additive constant because it appears in the equations solely through its gradient, and only Dirichlet boundary conditions for the velocity are prescribed.However, this constant pressure mode is not referred to as parasitic in the literature (Canuto et al. 1988).Two different techniques are implemented to solve the indeterminacy of the pressure.The first approach consists of imposing that the solution has a null projection on the checkerboard and constant pressure modes.The second method implies solving the augmented system where J is the Jacobian matrix, q = q n+1 − q n , N represents the discrete version of N (q n ), n includes the checkerboard and constant pressure modes, and β is an additional variable of the augmented system.For a more detailed explanation of this method, the reader is referred to Gustavsson (2003).The algorithm is considered to reach convergence when 1 where R is the discrete form of the operator R.
For different values of δ and Re, the steady state computed with the present code is compared with that obtained by Canton et al. (2017) using the in-house developed, unstructured, low-order, finite-element code PaStA.The Euclidean norm per grid point of the relative difference between the fields computed with the two codes is at most of the order of 10 −3 , whereas the relative difference between the friction factors is found to be of the order of 10 −4 .Good agreement is obtained for both the location and magnitude of the maxima of streamwise and in-plane velocities.It is worth noting that to allow the comparison between the basic states, linear interpolation between the finite-element mesh and the spectral grid is performed, which introduces additional sources of differences.The accuracy of all computed base flows is related to R N 2 , and thus to the tolerance chosen for the Newton-Raphson method.

Tracking of the global neutral stability curve
The stability of the flow in a torus depends on three main parameters: Re, δ and α.The critical Reynolds number is computed by employing a continuation algorithm that exploits the Newton-Raphson method to find the zeros of σ (Re) for different values of the curvature δ and the streamwise wavenumber α.The eigenvalues are obtained using either the shift-and-invert Arnoldi method implemented in the MATLAB command eigs or the generalised Rayleigh quotient iteration described in algorithm 1.
The latter approach has a lower computational cost and allows tracking a family of marginally stable modes at different curvatures or streamwise wavenumbers.Neutral stability curves Re cr (δ) are computed for several fixed values of the streamwise wavenumber using this method.The envelope of these curves represents the global neutral stability curve of the flow.The critical Reynolds number is also computed as a function of α for fixed values of the curvature to obtain an improved initial guess for tracking the least stable mode.Since the steady base flow in a torus cannot be derived analytically, in contrast with other flows (e.g.channel and straight pipe flows), the continuation method is coupled with the steady-state solver (which is also based on the continuation of the previous solution) to compute the basic state at each iteration of the algorithm.
The numerical method is validated against previous results in the literature.The comparison between the spectrum computed with the present code for δ = 0.01, Algorithm 1 Generalised Rayleigh quotient iteration Given the generalised eigenvalue problem Ax = λBx Input: A, B, x 0 Normalise initial vector Normalise eigenvector Re = 2150, 0 ≤ α ≤ 1 and that obtained by Canton et al. (2016) for the same parameters is shown in figure 3. Note that the value of the Reynolds number in the present work is half of that in Canton et al. (2016) because of the different length scales used for normalisation.The present code is found able to capture correctly the unstable branch, and good agreement is observed for the remaining portion of the spectrum as well.

Global neutral stability curve
Modal stability analysis is carried out according to the theoretical framework presented in § 2.2.Preliminary computations are performed for values of the curvature one order of magnitude apart from each other to identify how the eigenvalue spectrum and the most unstable eigenmode change with δ.For all curvatures, the maximum value of the streamwise wavenumber considered at this stage is at least α = 1 (the wavelength is normalised with the pipe radius).The MATLAB built-in function eigs is used to compute the portion of the spectrum closest to marginal stability.A convergence study is also performed to determine the required spatial resolution for capturing the instability.A portion of the eigenvalue spectrum for δ = 10 −5 , Re = 14 550, 0.005 ≤ α ≤ 1.55 is shown in figure 4 for three spatial resolutions.Two different branches of eigenvalues can be observed close to criticality, i.e. σ = 0.A spectral grid with n r × n θ = 25 × 50 is sufficient for converging the branch characterised by lower angular frequency.This branch is associated with the centre modes defined in the following.On the other hand, the branch with higher values of ω, associated with wall modes, requires spatial resolution n r × n θ = 35 × 70.The spatial structure of the eigenmodes belonging to these two branches is described in more detail in the following.
At this stage of the analysis, a bisection method is used to obtain an estimate of the critical Reynolds number for a given curvature.However, computing accurately Re cr with this approach is considerably expensive because of the slow convergence rate of the method.Better and faster computation of the neutral stability curve can be obtained by employing a continuation algorithm that exploits the Newton-Raphson method, as described in § 3.3.
The global neutral stability curve obtained pursuing the latter approach for 10 −7 ≤ δ ≤ 10 −2 is shown in figure 5, together with that computed by Canton et al. (2016) for 10 −2 ≤ δ ≤ 1.For curvatures between 1.5 × 10 −4 and 10 −2 , the critical Reynolds number increases smoothly with decreasing curvature.The critical eigenmodes are antisymmetric with respect to the equatorial plane and located predominantly in the central part of the cross-section, as shown in figure 6 for the case at δ = 10 −3 , Re = 3575 and α = 0.338.These modes are referred to as centre modes.Their tracking is carried out on a spectral grid with n r × n θ = 25 × 50, according to the resolution study discussed previously.
The global neutral stability curve exhibits a kink at curvatures close to 1.5 × 10 −4 since a different family of modes is found at criticality for 1.5 × 10 −6 δ 1.5 × 10 −4 .They are symmetric with respect to the equatorial plane and concentrated mainly in the near-wall region, towards the outer side of the bend.Therefore, they are referred to as wall modes.An example of the spatial structure of these modes is shown in figure 7, which displays the real part of the critical eigenmode for δ = 10 −4 , Re = 5338 and α = 1.848.The wall modes are associated with higher absolute values of the streamwise wavenumber α and the angular frequency ω than the centre modes.For δ 1.5 × 10 −4 , the neutral stability curve associated with this family of modes shows that they become unstable for higher values of Re than the centre modes.As stated previously, a higher spatial resolution is employed for tracking wall modes.The need for finer grid spacing is likely ascribed to their spatial structure, which exhibits a thin layer close to the wall.
For curvatures below approximately 1.5 × 10 −6 , the critical Reynolds number follows the power-law trend Re cr ∝ δ −1/2 .This result implies that Re cr diverges to infinity as δ → 0 and that the Dean number De = 2 Re √ δ at criticality is constant, approximately equal to 113, for these values of δ.The critical eigenmodes correspond to the centre modes described above since the wall modes become unstable at higher values of the Reynolds number in this range of curvatures.
For curvatures below approximately 10 −4 , the centre modes exhibit a peculiar behaviour: they become unstable at a certain Reynolds number, and restabilise for higher values of Re, becoming eventually unstable as the Reynolds number is increased further.To better understand this phenomenon, the growth rate associated with this family of modes is computed for several values of α and Re by means of the generalised Rayleigh quotient iteration.In this way, only the centre modes are tracked.
Neutral stability curves in the Re-α plane are identified by considering the flow as marginally stable when the absolute value of the computed growth rate is below 10 −9 .They are shown in figure 8 for δ = 10 −3 , 10 −4 , 10 −5 .The range of considered streamwise wavenumbers varies depending on the curvature.For δ = 10 −3 , the centre modes do not restabilise after becoming unstable.On the other hand, for δ = 10 −4 , a small detached region of instability is observed for 0.135 α 0.185 and 6200 Re 7500.This island of instability arises because the growth rate of the centre modes does not exhibit a monotonic trend with Re after crossing the marginal stability limit for the considered parameters.Hence, the centre modes become unstable and then restabilise in this region of the parameter space.A similar behaviour is also observed at δ = 10 −5 for 0.0375 α 0.0725.For this value of the curvature, the analysis is also carried out using higher spatial resolution in both radial and azimuthal directions, confirming the robustness of the results.Islands of instability have not been observed by Canton et al. (2016), most likely because they are characteristic of low curvatures.Nevertheless, they appear in spiral Couette flow (Meseguer & Marques 2000) and spiral Poiseuille flow (Meseguer & Marques 2002) as a result of the competition between different instability mechanisms.However, it needs to be remarked that they arise within the same family of modes in the case under investigation.
The streamwise wavenumber α and the phase speed v p = ωR p /α of the critical modes are reported in figure 9.For curvatures at which centre modes occur at criticality, the streamwise wavenumber exhibits a clear decreasing trend as δ decreases.In particular, for δ 1.5 × 10 −6 , the critical streamwise wavenumber is proportional to the square root of the curvature.According to Trefethen, Trefethen & Schmid (1999), the least stable mode in a straight pipe occurs for α / = 0. Hence, it is expected that the critical α for toroidal pipes approaches this value smoothly as δ → 0. For δ ≈ 1.5 × 10 −6 and δ ≈ 1.5 × 10 −4 , step discontinuities can be observed in the critical streamwise wavenumber since wall modes, characterised by larger absolute values of α, are the least stable for 1.5 × 10 −6 δ 1.5 × 10 −4 .For the wall modes, the streamwise wavenumber at criticality decreases as δ is reduced, but at a lower rate than for centre modes.The wavelength of wall modes is thus less affected by the curvature.
The critical phase speed increases slightly with decreasing curvature for 1.5 × 10 −4 δ 10 −2 , and it exhibits discontinuities at δ ≈ 1.5 × 10 −4 and δ ≈ 1.5 × 10 −6 due to the switching between the two families of critical modes.For the wall modes, v p decreases as the curvature is lowered, and its value is close to the amplitude of the base flow streamwise velocity in the region close to the outer wall, where the eigenmodes are predominant.On the other hand, for curvatures below approximately 1.5 × 10 −6 , the phase speed at criticality is almost constant, approximately equal to 1.5U b .Given that the critical modes  for δ 1.5 × 10 −6 have their core close to the centre of the pipe, this value of the critical phase speed is not surprising since it is almost equal to the magnitude of the streamwise velocity component of the base flow in this region.

Influence of the base flow
The stability analysis of toroidal pipe flows requires the numerical computation of the steady base flow, which increases both the complexity and computational cost of the algorithm.This task can be avoided if a sufficiently accurate eigenvalue spectrum can be obtained using an approximate basic state.In this section, the eigenvalue spectrum computed using different approximate base flows is compared to that obtained by linearising the incompressible Navier-Stokes equations about the numerically computed steady state.

Semi-analytic base flow
An approximate base flow can be derived by expanding the steady incompressible Navier-Stokes equations in powers of δ.Considering only the leading-order terms, after some algebraic manipulations, the equations can be written as follows.For the s-momentum, De (5.1b)For the θ -momentum, The incompressibility constraint is (5.1d) Numerical evidence shows that F ∝ √ δ for a fixed value of De (see figure 10).Therefore, for a given Dean number, the field (U s , U r / √ δ, U θ / √ δ, P/δ) T is a solution of (5.1) for any combination of δ and Re.
The higher-order terms in the expansion can be neglected for curvatures approaching zero, and an approximate steady base flow for a given combination of curvature and Reynolds number can then be derived simply by scaling the solution obtained for different values of δ and Re yielding the same Dean number.This result confirms the conclusion of Canton et al. (2017), i.e. that the Dean number similarity does not hold for any value of De, even at low curvatures.It needs to be remarked that the semi-analytic scaling introduced here can be applied only for low curvatures, as in the solution proposed by Dean (1927Dean ( , 1928)), but it is valid for any Dean number.Additionally, since the critical Dean number is constant for low curvatures, the proposed scaling implies that the magnitude of the secondary flow at criticality decreases with decreasing δ.

Eigenvalue spectra
The base flow for low curvatures can be constructed synthetically from previously computed steady solutions by exploiting the semi-analytic scaling described in § 5.1, thus avoiding its explicit numerical computation.However, it is not guaranteed that the accuracy of the base flow asymptotic expansion carries over to the eigenvalue spectrum.Indeed, highly non-normal operators can be significantly affected by small perturbations (Trefethen et al. 1993).The effect on the eigenvalues of having an approximate basic state is thus assessed for δ = 10 −7 , Re = 500 000.Three different approximate basic states are employed in the stability analysis: a base flow built from the solution at δ = 10 −5 , Re = 50 000 exploiting the scaling described above, the Hagen-Poiseuille flow (retaining non-zero curvature in the linearised incompressible Navier-Stokes equations), and the numerically computed base flow with in-plane velocity components U r and U θ set to zero.The leading eigenvalues for these basic states are reported in figure 13, together with the spectrum for the numerically computed base flow.Using the Hagen-Poiseuille flow as basic state completely fails to capture the spectrum, and considerably unstable eigenvalues are observed.This discrepancy is ascribed to the substantial difference between the base flow at δ = 10 −7 , Re = 500 000 and the laminar flow in a straight pipe.Setting the in-plane velocity components to zero in the numerically computed base flow does not provide satisfactory results either.This finding implies that the secondary motion in the basic state has a crucial role in determining the stability properties of the flow, even though it has low intensity (see also Canton et al. 2017).Nevertheless, as stated above, the magnitude of the secondary flow at criticality decreases as δ → 0. Therefore, better agreement is expected for lower curvatures.
No substantial differences are observed between the eigenvalue spectra computed using the scaled basic state and the numerically computed base flow.Hence, the truncation of the asymptotic expansion does not have a significant influence on the stability properties.
It should be pointed out that it is not possible to perform the stability analysis using Dean's expansion of the base flow (Dean 1927) since it is valid only for De 37.94.Thus, it loses any meaning at criticality, where De cr ≈ 113.Indeed, at Dean numbers as high as those on the global neutral stability curve, the base flow derived by Dean exhibits reversed flow, which is not physical in the considered configuration.However, the semi-analytic scaling described in § 5.1 can be applied in these cases.

Discussion and outlook
The present work extends the modal stability analysis of toroidal pipe flows carried out by Canton et al. (2016) to values of the curvature δ approaching zero, i.e. a straight pipe.The global neutral stability curve is traced through a continuation algorithm.Because of the properties of the geometry, the critical modes are either symmetric or antisymmetric with respect to the equatorial plane, depending on the value of the curvature.Although the flow in a toroidal pipe undergoes subcritical transition for δ 0.028, the present work highlights that the critical Reynolds number for modal instability diverges to infinity as the curvature approaches zero.The Dean number at criticality is approximately equal to 113 for low curvatures.For sufficiently low values of δ, the streamwise wavenumber of the critical mode is proportional to √ δ, whereas the critical phase speed plateaus at a value of approximately 1.5U b .These findings constitute further proof of the fact that the Hagen-Poiseuille flow is indeed stable to infinitesimal disturbances for any Reynolds number, as stated by Meseguer & Trefethen (2003) based on stability calculations.They also indicate a continuous transition of the stability properties of curved pipe flows towards the straight pipe limit.
The influence on the eigenvalue spectrum of having an approximate base flow is also investigated.Asymptotic expansions, corroborated by numerical evidence, show that the field (U s , U r / √ δ, U θ / √ δ, P/δ) T is an approximate self-similar solution of the steady incompressible Navier-Stokes equations for any combination of δ and Re corresponding to a fixed Dean number.The error of the approximation is found to be negligible for very low curvatures.For a given Dean number, an approximate base flow can then be derived from a solution at different curvature and Reynolds number.For δ = 10 −7 and Re = 500 000, this approximate basic state provides an eigenvalue spectrum in good agreement with that obtained using the numerically computed base flow.On the other hand, using the Hagen-Poiseuille flow as basic state and retaining non-zero curvature in the linearised incompressible Navier-Stokes equations does not yield accurate results, implying that the secondary flow is pivotal in the stability analysis.Neither does the stability analysis for the numerically computed base flow with the in-plane velocity components set to zero correctly capture the eigenvalue spectrum, confirming this finding.
Finally, it is worth mentioning that a validation of the present results by means of direct numerical simulations or experiments is not straightforward.Indeed, assuming that one can eliminate all the disturbances leading to subcritical transition, a complete torus of length L t = 2πR p /δ needs to be investigated in order not to miss any potentially relevant wavenumber.Considering an experimental facility consisting of a pipe with cross-section diameter 1 m (the pipe facility at CICLoPE has diameter 90 cm; see Örlü et al. 2017), a complete torus with curvature δ = 10 −7 will have length approximately 31 000 km, which is of the same order of magnitude of Earth's circumference (approximately 40 000 km).

Figure 1 .
Figure1.Sketch of the toroidal pipe showing the radius of the cross-section of the pipe R p , the curvature radius of the torus centreline R c , and the reference system {s, r, θ}.

Figure 2 .
Figure2.Sketch of the polar mesh for the Fourier-Chebyshev spatial discretisation on the pipe cross-section.In this example, n r = 3 and n θ = 6.The node ordering is also indicated.

Figure 3 .
Figure 3. Portion of the eigenvalue spectrum for δ = 0.01, Re = 2150, 0 ≤ α ≤ 1.The region close to the unstable branch is shown.The black dashed line indicates marginal stability (σ = 0).Black crosses indicate results from present code; maroon circles indicate data from Canton et al. (2016).

Figure 9 .
Figure 9. (a) Streamwise wavenumber α and (b) phase speed v p of the critical modes as functions of the curvature.Solid lines indicate data for 10 −7 ≤ δ ≤ 10 −2 (computed in the present work), whereas dashed lines show data adapted from Canton et al. (2016) for 10 −2 ≤ δ ≤ 1.

Figure 10 .
Figure 10.Amplitude of the volume force F for steady base flows at De = 200 and different values of the curvature.The dashed line indicates the trend √ δ.
Numerical verificationNumerical computations of steady base flows Q = (U, P) T for different combinations of δ and Re corresponding to the same Dean number support the semi-analytic scaling derived above.Figure11shows the absolute value of the difference between the scaled base flows (U s , U r / √ δ, U θ / √ δ, P/δ) T at De = 200 computed for δ = 10 −6 and δ = 10 −8 .The maximum difference between the two cases is of the order of 10 −5 .The energy norm of the relative difference between scaled base flows separated by one decade in δ at De = 200 decays linearly as the curvature decreases, as displayed in figure12.

Figure 11 .Figure 12 .
Figure 11.Absolute value of the difference of the scaled velocity components and pressure at De = 200 between the steady base flows at δ = 10 −6 , Re = 10 5 and δ = 10 −8 , Re = 10 6 , i.e. separated by one order of magnitude in Re.The inner wall of the bend is located on the bottom, whereas the top corresponds to the outer wall.