2 Numerical set-up

The fluid inside the torus is considered viscous and incompressible and therefore satisfies the incompressible Navier–Stokes equations

where $\boldsymbol{u}(\boldsymbol{r},t)$ is the velocity vector, with components $\boldsymbol{u}=(u_{r},u_{{\it\theta}},u_{z})$ in cylindrical coordinates, $p$ is the (dimensionless) pressure which has been rescaled by dividing by the Reynolds number $\mathit{Re}$ and $\boldsymbol{f}$ represents a volume force. The Reynolds number is based on the diameter of the pipe and the bulk velocity, i.e.  $\mathit{Re}=2R_{p}U_{b}/{\it\nu}$ .

The steady solutions and the stability results presented in this paper have been computed with PaStA, an in-house developed software package. PaStA is a Fortran90 code, written in primitive variables and based on the finite element method (FEM) employing triangular $\mathbb{P}_{2}{-}\mathbb{P}_{1}$ elements; it relies on the MUMPS library (Amestoy, Duff & L’Excellent Reference Amestoy, Duff and L’Excellent2000) for the solution of linear systems (for details, see Canton Reference Canton2013). The time-dependent simulations, on the other hand, have been performed by means of the spectral element code nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008). Convergence and accuracy tests have been performed to a more than satisfactory degree by comparing results to reference ones, whenever possible, or by comparing the two codes when results were new. These tests are omitted from this article for conciseness; some can be found in Canton (Reference Canton2013), Noorani et al. (Reference Noorani, El Khoury and Schlatter2013). However, we demonstrate excellent agreement between the linear analysis and the nonlinear DNS performed with two different codes, employing two different methods, in § 5.

The stability of the flow is studied in a classical framework: as a first step, a fixed point for the system, i.e. a steady solution to (2.1), is computed with the nonlinear solver, which employs Newton’s method. Subsequently, its stability properties are investigated by linearising the system about this steady state, and computing the leading eigenvalues of the linearised Navier–Stokes operator, $\mathscr{L}$ . By employing an expansion in normal modes, the perturbation fields are defined as

where $k\in \mathbb{R}$ is the streamwise wavenumber, ${\it\omega}={\it\omega}_{\mathit{r}}+\text{i}{\it\omega}_{\mathit{i}}\in \mathbb{C}$ is the eigenvalue and $\hat{\boldsymbol{u}}\in [\mathbb{C}]_{2}^{3}$ and $\hat{p}\in [\mathbb{C}]_{2}$ form the eigenfunction $\hat{\boldsymbol{x}}=(\hat{\boldsymbol{u}},\hat{p})$ . Note that here, and in the following, the symbols $[\mathbb{A}]_{d}^{c}$ are used to denote a $c$ -dimensional vector space $\mathbb{A}$ in a $d$ -dimensional physical space; it is clear that when $c$ is omitted $\mathbb{A}$ represents a scalar space. Although $k$ is in principle a real number, the $2{\rm\pi}$ -periodicity of the torus restricts it to $\mathbb{Z}$ .

The linearised Navier–Stokes equations, discretised via the finite element method, are thus reduced to a generalised eigenvalue problem of the form

where $\hat{\boldsymbol{x}}$ is the discretised eigenvector, $\unicode[STIX]{x1D648}\in \mathbb{R}^{\mathit{dof}\times \mathit{dof}}$ is the (singular) generalised mass matrix, $\unicode[STIX]{x1D647}_{k}\in \mathbb{C}^{\mathit{dof}\times \mathit{dof}}$ represents the discretisation of the linearised Navier–Stokes operator $\mathscr{L}$ , parametrised by the streamwise wavenumber $k$ . Finally, ‘ $\text{dof}$ ’ is a shorthand for the total number of degrees of freedom, the number of unknowns of the problem, which accounts for both velocity and pressure, i.e.  $\text{dof}=c(\mathbb{P}_{2}~\text{nodes})+(\mathbb{P}_{1}~\text{nodes})$ . The eigensolutions are computed to machine precision through an interface to the ARPACK library (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998), where the shift-invert spectral transformation is used to accelerate convergence to the desired portion of the spectrum. The eigenmodes determine whether the fixed point is stable or not, to locate the non-hyperbolic points of the system and, in the event of a bifurcation point, to determine the type of the bifurcation.

Once a bifurcation is identified and its type determined, its neutral curve can be traced in parameter space, separating the stable and unstable regions. In order to track a bifurcation, an augmented set of equations describing the system at the bifurcation point, i.e. where the critical eigenvalue has zero growth rate, needs to be solved. Equation (2.4a–c ) presents this system in the case of a Hopf bifurcation, which is the type of bifurcation encountered by this flow,

where $\unicode[STIX]{x1D649}=\unicode[STIX]{x1D649}(\mathit{Re},{\it\delta},\boldsymbol{x})\in \mathbb{R}^{\mathit{dof}\times \mathit{dof}}$ , $\unicode[STIX]{x1D648}=\unicode[STIX]{x1D648}({\it\delta})\in \mathbb{R}^{\mathit{dof}\times \mathit{dof}}$ , and $\unicode[STIX]{x1D647}_{k}=\unicode[STIX]{x1D647}_{k}(\mathit{Re},{\it\delta},\boldsymbol{x})\in \mathbb{C}^{\mathit{dof}\times \mathit{dof}}$ . The unknowns for this system are: the base flow, i.e. the steady solution to (2.1), comprising both velocity and pressure, $\boldsymbol{x}=(\boldsymbol{u},\boldsymbol{p})\in \mathbb{R}^{\mathit{dof}}$ ; the frequency of the critical eigenmode, ${\it\omega}_{\mathit{r}}\in \mathbb{R}$ ; and the corresponding eigenvector, $\hat{\boldsymbol{x}}=(\hat{\boldsymbol{u}},\hat{\boldsymbol{p}})\in \mathbb{C}^{\mathit{dof}}$ . The pair $\mathit{Re},{\it\delta}\in \mathbb{R}$ can act as unknown/parameter or, vice versa, as parameter/unknown, where one is fixed while the other is obtained as part of the solution. Equation (2.4a ) is a real equation which determines the base flow: $\unicode[STIX]{x1D649}$ constitutes the discretised steady-state Navier–Stokes equations, i.e. it is shorthand for a steady solution to (2.1). This is the nonlinear equation that is solved when computing a steady state, regardless of whether it is on the neutral curve or not. Equation (2.4b ) is a complex equation which represents the eigenvalue problem (2.3) in the particular case where a pair of complex conjugate eigenvalues has zero growth rate. This equation forces the solution of the system to be on the neutral curve. Equation (2.4c ) is a complex equation as well, it fixes the phase and amplitude of the eigenvector. The constant vector ${\bf\phi}$ is chosen as the real part of the initial guess for $\hat{\boldsymbol{x}}$ , so that this linear equation mimics an $L^{2}$  norm.

Newton’s method is employed for this large nonlinear system as well. In contrast with the standard procedure, though, the full Jacobian matrix for (2.4a–c ) is not assembled, as its size would be $3\text{dof}+2\times 3\text{dof}+2$ , and its sparsity pattern would need to be created. Instead, a block Gauss factorisation and linear algebra are used to split the resulting system into five $\text{dof}\times \text{dof}$ linear systems, two real and three complex. In addition to reducing memory requirements, these systems have matrices with the same sparsity pattern as those already used for the solution of the steady state and the eigenvalue problem. The method is similar to those presented in Salinger et al. (Reference Salinger, Bou Rabee, Pawlowski, Wilkes, Borroughs, Lehoucq and Romero2002), but here the Jacobian matrix for the steady state is not equivalent to the linearised Navier–Stokes operator $\unicode[STIX]{x1D647}_{k}$ . The tolerance for the iterative solution of (2.4a–c ) was chosen as to have an uncertainty of $\mathit{Re}$ on the neutral curve of $\pm 10^{-4}\,\%$ on a mesh with 98 870 elements, corresponding to $\text{dof}=646\,327$ .

3 Base flow

The fluid, flowing inside the pipe at constant mass rate, is maintained in motion by a volume force aligned with the centreline of the torus and defined as $\boldsymbol{f}=\boldsymbol{e}_{{\it\theta}}F/r$ . Here $F$ is a scalar constant and $r$ represents the distance from the axis of the torus. The inverse proportionality to $r$ is required so that the integral of $\boldsymbol{f}$ along the curvilinear coordinate $\boldsymbol{s}=r{\bf\theta}$ is independent of $r$ , i.e.  $\int \boldsymbol{f}\boldsymbol{\cdot }\text{d}\boldsymbol{s}=\int (F/r)r\boldsymbol{e}_{{\it\theta}}\boldsymbol{\cdot }\text{d}{\bf\theta}=F({\it\theta}-{\it\theta}_{0})~\forall r$ . This property confers on $\boldsymbol{f}$ a role similar to that of the pressure in straight pipes (or channels) and it ensures continuity with a straight pipe for a torus in the limit of vanishing curvature. In fact, when ${\it\delta}\rightarrow 0$ , $r\rightarrow \infty$ on the section of the toroidal pipe, then $\boldsymbol{f}$ tends to a value which is uniform on the section of the pipe, as the pressure in a straight pipe.

The value of $F$ is determined via the secant method such as to have unitary bulk velocity. This cannot be done analytically since the relationship between $\boldsymbol{f}$ and $U_{b}$ is nonlinear and a priori unknown (see Ito Reference Ito1959).

The steady-state solution inherits both the rotational invariance and the mirror symmetry of the toroidal geometry and it is thus independent of ${\it\theta}$ and symmetric with respect to the equatorial plane of the torus. In addition to being evident from geometrical considerations, this has been confirmed by three-dimensional DNS performed for values of ${\it\delta}$ and $\mathit{Re}$ below criticality. The invariance with respect to ${\it\theta}$ allows the solutions to be computed on a two-dimensional section (retaining three velocity components), sensibly reducing the computational cost.

The base flow is characterised by the presence of two symmetric, counter-rotating Dean vortices. The side with the counterclockwise rotating vortex is visible in the contours of the streamfunction in figure 2(a,c), and the centres of the vortices are indicated with plus symbols. These vortices develop as a result of the centrifugal force arising from breaking the central symmetry characterising a straight pipe and are thus present at every Reynolds number and for any non-zero curvature. The question of a critical Dean number below which no recirculation occurs has been discussed in textbooks (see, for instance, Isachenko, Osipova & Sukomel Reference Isachenko, Osipova and Sukomel1974, p. 218). However, even for ${\it\delta}=10^{-11},~\mathit{Re}=1$ , corresponding to $De=\mathit{Re}\sqrt{{\it\delta}}\approx 3\times 10^{-6}$ , the vortices are present, albeit exhibiting low intensity; their presence for any non-zero curvature can also be shown analytically, cf. Dean (Reference Dean1927).

As can be observed in figure 2, which depicts two steady states for the same Dean number, $De=759$ , the curvature of the torus plays a primary role in determining the behaviour of the base flow, starting at least from ${\it\delta}=0.1$ . It greatly influences, among other features, the shape of the vortices, the position of the maximum velocity and the magnitude of the in-plane flow. An example can be observed in figure 2: the maximum of $u_{{\it\theta}}$ moves from the outer portion of the pipe for low curvature, figure 2(a), to the centre of the vortices, 2(c), for high curvature. Moreover, as discussed in the next section, the curvature is of primary importance for the characteristics of the first instability. We thus observe that the curvature of the pipe should be considered as a separate parameter with as much importance as the Reynolds number and that the Dean number alone is insufficient to describe the flow fully.

Figure 2.

Steady state solutions for $De=759$ at different curvatures and Reynolds number: (a,b) ${\it\delta}=0.1$ , $\mathit{Re}=2400$ ; (c,d) ${\it\delta}=0.9$ , $\mathit{Re}=800$ . Panels (a,c) depict the magnitude of the streamwise velocity component $u_{{\it\theta}}$ , and contours of the in-plane streamfunction ${\it\psi}$ for positive values; the plus symbols are located at the position of the maximum and minimum of ${\it\psi}$ . Panels (b,d) show the magnitude of the in-plane velocity. Quantities are normalised with the bulk velocity. The ‘I/O’ labels indicate the inner and outer portions of the pipe (see figure 1).

Figure 3.

Eigenvalue spectra for the flow in a toroidal pipe. (a) Presents the spectrum for ${\it\delta}=0.01,\mathit{Re}=4300$ , corresponding to $De=430$ , while (b) is for ${\it\delta}=0.3,\mathit{Re}=3400$ , corresponding to $De\approx 1862$ . Both cases are moderately above criticality: the flow for ${\it\delta}=0.01$ becomes unstable for $\mathit{Re}>4257$ , while for ${\it\delta}=0.3$ $\mathit{Re}_{\mathit{c}}=3379$ . The least-stable eigenvalues are highlighted with plus symbols. The wavenumber of the critical mode for ${\it\delta}=0.01$ is 74, while for ${\it\delta}=0.3$ $k=7$ . The corresponding eigenvectors are depicted in figure 4. Note the striking difference in the structure of the spectra.

4 Stability analysis

The first step of the stability analysis was the computation of eigenvalue spectra on a coarse grid of points in parameter space, $({\it\delta},\mathit{Re})\in \left([0.005,1],[500,10^{4}]\right)$ . All curvatures in the grid revealed spectra with unstable eigenvalues for $\mathit{Re}<6\times 10^{3}$ . This result confirms the numerical observations by Di Piazza & Ciofalo (Reference Di Piazza and Ciofalo2011) as well as the experiments by Kühnen et al. (Reference Kühnen, Holzner, Hof and Kuhlmann2014, Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) and extends these studies, which are limited to the low curvature range. Most importantly, it proves that this flow is linearly unstable, in obvious contrast to the flow inside a straight pipe (see Schmid & Henningson (Reference Schmid and Henningson1994) and citations therein). Each spectrum in the grid exhibits a pair of complex conjugate eigenvalues as most unstable mode, indicating that the cause of the modal instability is a supercritical Hopf bifurcation. In fact, DNS confirm that the flow transitions from a steady state to a periodic state, characteristic of the stable limit cycle identified by this bifurcation. For a thorough discussion, see § 5.

Figure 3 presents, as examples, two supercritical eigenvalue spectra. These are for very low and moderate curvature, ${\it\delta}=0.01$ and $0.3$ , respectively, and each has a complex conjugate eigenpair with positive growth rate. Although the instability is caused by a Hopf bifurcation for every curvature, the associated critical mode changes with the curvature of the pipe. This can be seen in figure 4, which shows the different spatial structure of the two critical modes presented in figure 3. The eigenmodes differ in their wavenumbers, temporal frequencies, spatial structure and symmetry properties. Figure 3 demonstrates a striking difference in the structure of the two spectra and the variation of the frequency of the critical mode. The set of spectra computed on the grid revealed a non-monotonic relationship between critical Reynolds number and curvature; even the variation of eigenmode characteristics, e.g. wavenumber, as a function of $({\it\delta},\mathit{Re})$ does not have a continuous trend. It is clear that a detailed analysis of the whole parameter space cannot be performed employing a coarse grid of points and interpolation between them.

Figure 4.

Critical modes for (a) ${\it\delta}=0.01,~\mathit{Re}=4257$ and (b) ${\it\delta}=0.3$ , $\mathit{Re}=3379$ . Both modes are antisymmetric with respect to the equatorial plane of the torus. The wavenumbers are $k=74$ and $k=7$ respectively. Isocontours of opposite values of the streamwise velocity component are depicted in blue and red, $u_{{\it\theta}}=\pm 0.15$ in (a) and $u_{{\it\theta}}=\pm 2.2$ in (b). The vertical cuts depict the streamwise velocity and in-plane streamlines of the corresponding base flow. The fluid is flowing clockwise, as indicated by the arrows.

Employing the method described in § 2, the neutral curve of every critical mode has thus been traced in the parameter space as a function of curvature and Reynolds number. The result is a complex picture, depicted in figure 5. It presents five ‘families’ and three ‘isolated’ modes; their envelope constitutes the global neutral curve for the flow. All eigenmodes on this curve share the characteristic of representing a travelling wave, while all of the other properties are modified along the neutral curve. As anticipated in § 3, the curvature of the pipe plays a key role in determining the identity of the critical mode and its characteristics. The neutral curves cannot, in fact, be collapsed onto one another as a function of a single non-dimensional parameter, thereby confirming the need for two parameters ( ${\it\delta}$ and $\mathit{Re}$ ) to describe this flow.

A ‘family’ comprises eigenmodes that share common characteristics, while the eigenvalues designated as ‘isolated’ are those which contribute to the global neutral curve while being dissimilar to the modes belonging to the neighbouring neutral curves. In more detail, all modes belonging to a given family display the same spatial structure and symmetry properties, they have approximately equal phase speed and comparable wavelength. In addition, while the wavenumber does not vary monotonically along the envelope of neutral curves, it does so inside families. Furthermore, eigenvalues in the same family lie on the same branch. This is very clear for low curvatures: in figure 3(a) the two uppermost branches comprise eigenvalues belonging to the same family, the same holds for the two, nearly straight-line branches just below them, and so on. Possibly even more characteristic is the fact that the eigenmodes forming a family have neutral curves with a very similar trend $({\it\delta},\mathit{Re}({\it\delta}))$ . This can readily be observed in figure 5 where each neutral curve is purposely plotted beyond the envelope line to illustrate this feature.

Figure 5.

Neutral curve in the ${\it\delta}{-}\mathit{Re}$ plane. ${\it\delta}\in [0.002,1]$ . Each line corresponds to the neutral curve of one mode. The neutral curve for the flow is formed by the envelope of the lines. Five families (black and blue) and three isolated modes (green) are marked by labels. Symmetric modes are indicated with continuous lines while antisymmetric modes are represented with dashed lines. Note that the curves are not interpolated, i.e. they are segments connecting computed solutions with ${\rm\Delta}{\it\delta}=\mathit{O}(10^{-3})$ . The uncertainty on the Reynolds number is $\pm 10^{-4}\,\%$ .

5 Comparison with the nonlinear behaviour

Nonlinear DNS have been performed to confirm the results of the linear stability analysis. The values of ${\it\delta}$ and $\mathit{Re}$ of these simulations are reported in table 1. A randomly perturbed Poiseuille profile was used as the initial condition for all simulations. Each value of curvature investigated has a transition from a steady to a periodic regime, confirming the presence of an Hopf bifurcation. Figure 7 confirms it by depicting the amplitude of saturated flow oscillations as a function of $\mathit{Re}$ ; the result is $A\propto \sqrt{\mathit{Re }-\mathit{ Re }_{c}}$ , characteristic of a Hopf bifurcation (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983). Amongst others, the most relevant quantities used for the comparison are: the critical Reynolds number, the wavenumber of the travelling wave, its symmetry properties, growth rate and temporal frequency. The simulations show excellent agreement with the results of the linear analysis, confirming all properties related to the critical modes. While quantities such as the wavenumber, symmetry, growth rate and temporal frequency could be measured directly, the critical Reynolds number needed to be derived from the other quantities. To do so, the amplitude of the oscillating flow has been measured in simulations across the neutral curve. The envelopes of these oscillations have then been extracted from the time intervals where the signals exhibited exponentially growing, respectively decaying behaviour (depending on $\mathit{Re}$ , of course). The time-averaged growth rate and the Reynolds numbers of the two simulations have then been employed in a two-point, linear interpolation to estimate the Reynolds number corresponding to zero growth rate for the nonlinear flow. An example of this procedure is shown in figure 7 for ${\it\delta}=0.05$ . For this curvature the interpolation between the DNS for $\mathit{Re}=3700$ and $3750$ gives $\mathit{Re}_{\mathit{c}}\approx 3713.2$ , the linear analysis predicts a critical Reynolds number $\mathit{Re}_{\mathit{c}}=3713.36$ . All other quantities compared between nonlinear DNS and the linear analysis are affected by differences of the same order of magnitude.

Figure 7.

Saturation amplitude, $A$ , measured as a function of Reynolds number for curvature ${\it\delta}=0.05$ . The critical Reynolds number is $\mathit{Re}_{c}=3713$ . The round (black) marker corresponds to the critical $\mathit{Re}$ computed through the linear stability analysis. The dashed line represents a second-order polynomial fit to the DNS data; its intercept with the horizontal axis was computed by linearly interpolating the growth rate, respectively decay, of the oscillations measured in the simulations for $\mathit{Re}=3750$ , respectively $\mathit{Re}=3700$ .

Table 1.

Values of curvature and Reynolds number at which direct numerical simulations were performed. All DNS were performed in a full, three-dimensional torus.

6 Discussion and conclusions

The present results are compared with the literature in figure 8 (in the range where previous data are available). It can be noted that, in spite of the different methods employed, all previous works indicate subcritical transition for curvatures below 0.02–0.04. This phenomenon was expected since ${\it\delta}$ is nearing its lower limit and the effect that the centrifugal force has on the flow becomes less relevant. Moreover, the diameter of the torus is so large in comparison to that of the pipe that the feedback-like mechanism caused by the $2{\rm\pi}$ -periodicity in the streamwise direction becomes almost irrelevant to the flow.

Figure 8.

Comparison of the neutral curve with experimental and numerical data from the literature. The continuous (black) line represent the envelope of families F4S and F5A (the intersection is at ${\it\delta}=0.016$ , where the two dotted lines highlight the continuation of the two families). The data by Cioncolini & Santini (Reference Cioncolini and Santini2006) refers to the first discontinuity in their friction measurements, while the reported data by Sreenivasan & Strykowski (Reference Sreenivasan and Strykowski1983) is the curve they refer to as the ‘conservative lower critical limit’, and that by Noorani & Schlatter (Reference Noorani and Schlatter2015) is for sublaminar friction drag. The interpolants provided by Kühnen et al. (Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) are also depicted: with a dashed (blue) line the fit for the supercritical transition and with a dotted (blue) line the fit for subcritical transition. All results agree on subcritical transition for very low curvature, with the data by Kühnen et al. (Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) falling on top of that by Sreenivasan & Strykowski (Reference Sreenivasan and Strykowski1983).

Higher curvatures, on the other hand, present a different picture. Both Sreenivasan & Strykowski (Reference Sreenivasan and Strykowski1983) and Kühnen et al. (Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) appear to have identified family F4S. The former selected the appearance of the first ‘bursts’ near the outer wall as a definition for the critical Reynolds number. This explains the higher estimate of $\mathit{Re}$ where supercritical transition occurs, but also answers their question on the cause behind the non-monotonic behaviour of their data: the switch from subcritical to supercritical transition. Kühnen and coworkers, on the other hand, aware of the supercritical transition identified by Webster & Humphrey (Reference Webster and Humphrey1993), Di Piazza & Ciofalo (Reference Di Piazza and Ciofalo2011), report the Reynolds number corresponding to the onset of the periodic state. These measurements are in fair agreement with F4S. They confirm the relevance of the linear instability and that the first instability is caused by a Hopf bifurcation. The discrepancy in Reynolds number, and the trend of their interpolant may well be due to an optimistic estimate of the accuracy of their measurements, as well as the absence of data at higher curvatures. In fact, the DNS performed at the rightmost point of their dataset, ${\it\delta}=0.1,~\mathit{Re}=3202$ , confirms the absence of oscillations.

In conclusion, the toroidal geometry allows the isolation of the effect that ${\it\delta}$ has on the first instability. It will now be possible to separate these effects from others when investigating differently curved pipes. Results show that a very little curvature is sufficient to make this flow linearly unstable via a Hopf bifurcation that leads the flow to a periodic regime. The neutral curve associated with this bifurcation has been traced in parameter space exploring, for the first time, the whole range of possible curvatures. The Reynolds number and the wavelength of the critical modes diverge when approaching ${\it\delta}=0$ , suggesting a continuous transition to the straight pipe, which has modes with high ${\it\lambda}$ , i.e. quasi-streamwise independent, as the least-stable ones.