2.1 Numerical discretisation

The present analysis is performed via DNS of the incompressible Navier–Stokes equations. The equations are discretised with the spectral-element code Nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008) using a $\mathbb{P}_{N}-\mathbb{P}_{N-2}$ formulation. After an initial mesh-dependency study, the polynomial order was set to $N=8$ for the velocity and, consequently, $N=6$ for the pressure. We consider a $90^{\circ }$ bent pipe with curvature $\unicode[STIX]{x1D6FF}=0.3$ for a Reynolds number $Re=11\,700$ , corresponding to a friction Reynolds number $Re_{\unicode[STIX]{x1D70F}}\approx 360$ (referred to the straight pipe sections). A straight pipe of length $L_{i}=7D$ precedes the bent section (see § 3.1), and a second straight segment of length $L_{o}=15D$ follows it. Further details about the mesh, including element number and size, are reported in table 2. The supplementary video movie1.m4v shows the set-up and a visualisation of the flow.

Table 2.

Details of the mesh employed in the present work. $n_{el}$ corresponds to the number of hexahedral elements, while $n_{dof}$ is the total number of degrees of freedom (velocity and pressure). Quantities indicated as $(\cdot )^{+}$ are reported in viscous scaling, and the numbers between parenthesis correspond to minimum and maximum values.

2.2 Inflow boundary and divergence-free synthetic eddy method

Since the aim of the present work is to reproduce and study swirl-switching, a periodic or quasi-periodic phenomenon, the treatment of the inflow boundary is of utmost importance. The flow field prescribed at the inflow boundary should not introduce any artificial frequency, in order to avoid the excitation of unphysical phenomena or a modification of the frequencies inherent to the swirl-switching. A recycling method, as the one used by Rütten et al. (Reference Rütten, Meinke and Schröder2001, Reference Rütten, Schröder and Meinke2005) and Carlsson et al. (Reference Carlsson, Alenius and Fuchs2015), should therefore be avoided, as highlighted in § 1.

In the present work the velocity field at the inlet boundary of the straight pipe preceding the bend is prescribed via a divergence-free synthetic eddy method (DFSEM). This method, introduced by Poletto et al. (Reference Poletto, Revell, Craft and Jarrin2011) and based on the original work by Jarrin et al. (Reference Jarrin, Benhamadouche, Laurence and Prosser2006), works by prescribing a mean flow modulated in time by fluctuations in the vorticity field. The superposition of the two reproduces, up to second order, the mean turbulent fluctuations of a reference flow and requires a short streamwise adjustment length to fully reproduce all quantities. The fluctuations are provided by a large number of randomly distributed ‘vorticity spheres’ (or ‘eddies’) which are generated and advected with the bulk velocity in a fictitious cylindrical container located around the inflow section. When a sphere exits the container, a new, randomly located sphere is created to substitute it. The cylindrical container is dimensioned such that newly created eddies do not touch the inlet plane upon their creation and they have stopped affecting it before exiting the container, i.e. the cylinder extends from $-\max (D_{eddies})$ to $\max (D_{eddies})$ (see figure 5 in Poletto, Craft & Revell Reference Poletto, Craft and Revell2013, for an illustration of the container).

The random numbers required to create the fluctuations are generated on a single processor with a pseudo random number generator (Chandler & Northrop Reference Chandler and Northrop2003) featuring an algorithmic period of $2^{1376}$ iterations, large enough to exclude any periodicity in the simulations, which feature approximately $10\,000$ synthetic eddies. The size of the spheres is selected to match the local integral turbulence length scale, and their intensity is scaled to recover the reference turbulent kinetic energy, producing isotropic but heterogeneous second-order moments. The method prescribes isotropic turbulence, instead of the correct anisotropic variant, because it was shown that the former leads to a shorter adjustment length in wall-bounded flows (see figure 11 in Jarrin et al. Reference Jarrin, Benhamadouche, Laurence and Prosser2006). In order to satisfy the continuity equation, no synthetic turbulence is created below $(R-r)^{+}<10$ ; however, this does not significantly affect the adjustment length since the dynamics of the viscous sublayer is faster than the mean and converges to a fully developed state in a shorter distance. The turbulence statistics necessary for the method, specifically the mean flow $U(r)$ , the turbulent kinetic energy $k(r)$ and the dissipation rate $\unicode[STIX]{x1D716}(r)$ , were extracted from the straight pipe DNS performed by El Khoury et al. (Reference El Khoury, Schlatter, Noorani, Fischer, Brethouwer and Johansson2013). Section 3.1 presents the validation of our implementation of the DFSEM; more details can be found in Hufnagel (Reference Hufnagel2016).

2.3 Proper orthogonal decomposition

Besides point measures, we use POD (Lumley Reference Lumley, Yaglom and Tatarski1967) to extract coherent structures from the DNS flow fields and identify the mechanism responsible for swirl-switching. More specifically, we use snapshot POD (Sirovich Reference Sirovich1987) where $n$ three-dimensional, full-domain flow fields of dimension $d$ (corresponding to the number of velocity unknowns) are stored as snapshots. POD decomposes the flow in a set of orthogonal spatial modes $\unicode[STIX]{x1D731}_{i}(\boldsymbol{x})$ and corresponding time coefficients $a_{i}(t)$ ranked by kinetic energy content, in decreasing order. The most energetic structure extracted by POD corresponds to the mean flow and will be herein named ‘zeroth mode’, while the term ‘first mode’ will be reserved for the first time-dependent structure.

A series of instantaneous flow fields (snapshots) is ordered column-wise in a matrix $\unicode[STIX]{x1D64E}\in \mathbb{R}^{d\times n}$ and decomposed as:

where $\boldsymbol{U}\in \mathbb{R}^{d\times d}$ , $\unicode[STIX]{x1D72E}=\text{diag}(\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x1D70E}_{m},0)$ , with $m=\min (d,n)$ , and $\boldsymbol{V}\in \mathbb{R}^{n\times n}$ . The decomposition in (2.1) is obtained by computing the singular value decomposition (SVD) of $\unicode[STIX]{x1D648}^{1/2}\unicode[STIX]{x1D64E}\unicode[STIX]{x1D64F}^{1/2}$ , where $\unicode[STIX]{x1D648}$ is the mass matrix and $\unicode[STIX]{x1D64F}$ is the temporal weights matrix, which results in $\tilde{\unicode[STIX]{x1D650}}\unicode[STIX]{x1D72E}\tilde{\unicode[STIX]{x1D651}}^{\intercal }$ , where $\tilde{\unicode[STIX]{x1D650}}$ and $\tilde{\unicode[STIX]{x1D651}}$ are unitary matrices ( $\tilde{\unicode[STIX]{x1D650}}^{\intercal }\tilde{\unicode[STIX]{x1D650}}=\unicode[STIX]{x1D644}$ and $\tilde{\unicode[STIX]{x1D651}}^{\intercal }\tilde{\unicode[STIX]{x1D651}}=\unicode[STIX]{x1D644}$ ); the POD modes are then obtained as $\boldsymbol{U}=\unicode[STIX]{x1D648}^{-1/2}\tilde{\unicode[STIX]{x1D650}}$ and $\boldsymbol{V}=\unicode[STIX]{x1D64F}^{-1/2}\tilde{\unicode[STIX]{x1D651}}$ . To improve the convergence of the decomposition, we exploit the symmetry of the pipe about the $I$ – $O$ plane, which results into a statistical symmetry for the flow, and store an additional mirror image for each snapshot (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993).

Figure 2.

Recovery of fully developed turbulence statistics for the divergence-free synthetic eddy method at $Re_{\unicode[STIX]{x1D70F}}=360$ , compared to the reference values by El Khoury et al. (Reference El Khoury, Schlatter, Noorani, Fischer, Brethouwer and Johansson2013). Panel (a) shows the ratio between the DFSEM and the reference data as a function of streamwise distance from the inflow plane. The grey shaded area indicates a $\pm 1\,\%$ tolerance with respect to the data by El Khoury et al. (Reference El Khoury, Schlatter, Noorani, Fischer, Brethouwer and Johansson2013). Panels (b–d) show classical statistical profiles as a function of radial position at $s_{i}=5D$ . Solid lines indicate the reference data, while symbols represent the current results (note that the number of shown points is reduced and does therefore not represent the grid resolution; see table 2).

3.1 Inflow validation

An auxiliary simulation for $Re_{\unicode[STIX]{x1D70F}}=360$ was set-up to test the performance of the DFSEM in a $25D$ long straight pipe, provided with the same mesh characteristics used for the bent pipes. Classical statistical quantities were used for the validation and compared with the reference values by El Khoury et al. (Reference El Khoury, Schlatter, Noorani, Fischer, Brethouwer and Johansson2013). The comparison is presented in figure 2(a) as a function of distance from the inflow boundary, and shows that the DFSEM approaches a fully developed turbulent state (within $\pm 1\,\%$ of error) at approximately $5D$ from the inflow boundary. Figure 2(b–d) presents the velocity, stress profiles and the turbulent kinetic energy budget at the chosen streamwise position of $s_{i}=5D$ , which confirm the recovery of fully developed turbulence by the divergence-free synthetic eddy method.

A length of $7D$ was therefore chosen for the straight pipe preceding the bent section, in order to allow for some tolerance and to account for the (weak, up to $1D$ ) upstream influence of the Dean vortices (Anwer, So & Lai Reference Anwer, So and Lai1989; Sudo, Sumida & Hibara Reference Sudo, Sumida and Hibara1998). For comparison, the more commonly used approach where random noise is prescribed at the inflow requires a development length between 50 and $110D$ (Doherty, Monty & Chong Reference Doherty, Monty and Chong2007). POD modes were also computed to further check the correctness of this method. The results, not reported here for conciseness, were in good agreement with those presented by Carlsson et al. (Reference Carlsson, Alenius and Fuchs2015) for a periodic straight pipe, that is, streamwise invariant modes with azimuthal wavenumbers between 3 and 7.

3.2 Two-dimensional POD

Two-dimensional POD, considering all three velocity components, is employed as a first step in the analysis of swirl-switching. Instantaneous velocity fields are saved at a distance of $2D$ from the end of the bent section and are used, with their mirror images, to assemble the snapshot matrices (Berkooz et al. Reference Berkooz, Holmes and Lumley1993). The 1234 velocity fields used for the decomposition were saved at a sampling frequency of $St=0.25$ , and the sampling was started only after the solution had reached a statistically steady state. As a consequence of exploiting the mirror symmetry, all modes are either symmetric or antisymmetric, a condition to which they would converge provided that a sufficient number of snapshots had been saved.

Figure 3.

Pseudo colour of the streamwise velocity component and streamlines of the in-plane velocity components for the first three POD modes (i–iii). The modes are oriented as in figure 1(b). The snapshots were extracted at $s_{o}=2D$ .

The first three modes are shown in figure 3 by means of pseudo colours of their streamwise velocity component and streamlines of the in-plane velocity components. Two out of three modes are antisymmetric: (i, ii) and are in the form of a single swirl covering the whole pipe section, (i), and a double swirl, (ii), formed by two counter-rotating vortices disposed along the inner–outer direction on the symmetry plane. The third mode, (iii), resembles a harmonic of the Dean vortices.

These findings are in agreement with previous experimental work, such as that of Sakakibara et al. (Reference Sakakibara, Sonobe, Goto, Tezuka, Tada and Tezuka2010) and Kalpakli & Örlü (Reference Kalpakli and Örlü2013), which attributed the dynamics of swirl-switching to the antisymmetric modes. The frequency content of these modes is presented in figure 4, in terms of Welch’s power spectral density estimate for the time coefficients of the first three modes, corresponding to the structures shown in figure 3. It can be observed that the spectra have a low peak-to-noise ratio and that each mode is characterised by a different spectrum and peak frequency, in agreement with previous 2-D POD studies: see, e.g. figure 8 in Hellström et al. (Reference Hellström, Zlatinov, Cao and Smits2013), which presents peaks with similar values to the present ones, although their study was for a slightly larger curvature of $\unicode[STIX]{x1D6FF}=0.5$ . This fact has caused some confusion in the past, with disagreeing authors attributing different causes to the various peaks, without being able to come to the same conclusion about the frequency, or the structure, of swirl-switching. The reason is that swirl-switching is caused by a three-dimensional wave-like structure, as will be shown by 3-D POD in § 3.3, and a two-dimensional cross-flow analysis cannot distinguish between the spatial and temporal amplitude modulations created by the passage of the wave. A simple analytical demonstration of this concept is provided in the appendix A, and shows that conclusions drawn from a flow reconstruction based on 2-D POD modes (see, e.g. Hellström et al. Reference Hellström, Zlatinov, Cao and Smits2013) are incomplete.

Figure 4.

Welch’s power spectral density (PSD) estimate for the time coefficients $a_{i}$ of the most energetic 2-D POD modes. The frequencies are scaled with pipe diameter and bulk velocity. The 2-D modes were extracted at $s_{o}=2D$ . The markers and corresponding labels report the frequency of the highest peak of each spectrum.

3.3 Three-dimensional POD

For the 3-D POD, the same snapshots as for 2-D POD were used. In order to reduce memory requirements, the snapshots were interpolated on a coarser mesh before computing the POD. This is, however, not a problem because the swirl-switching is related to large-scale fluctuations in the flow.

Figure 5.

Panel (a) shows the four most energetic three-dimensional POD modes. The four longitudinal cuts show pseudo colours of the normal velocity component $u_{n}$ , while the eight cross-sections display the in-plane streamlines and are coloured by streamwise velocity $u_{s}$ . The supplementary material includes two videos showing the reconstruction of the flow based on these modes. Panel (b) shows the swirl intensity, measured by circulation $\unicode[STIX]{x1D6E4}$ , along the streamwise axis of the two most energetic modes, $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ , and their envelope. The spatially decaying, wave-like behaviour can be appreciated.

Figure 6.

Time coefficients $a_{i}(t)$ of the two most energetic three-dimensional POD modes. Panel (a) shows the temporal signal, which allows us to observe the qualitative quarter-period phase shift of mode 2 with respect to mode 1; panel (b) shows the (colour coded) time over coefficients $a_{1}$ and $a_{2}$ , illustrating the oscillating character. The time axis is (arbitrarily) cut at $t=50D/U_{b}$ for illustration purposes, the total recorded signal is over $300D/U_{b}$ .

Figure 7.

Welch’s power spectral density estimate for the time coefficients $a_{i}$ of the most energetic 3-D POD modes. The frequencies are scaled with pipe diameter and bulk velocity. The markers and corresponding labels report the frequency corresponding to the peak of each spectrum. The range of the Strouhal number is identical to that of figure 4 to ease comparison.

The four most energetic modes are depicted in figure 5 by means of pseudo colours of normal and streamwise velocity components, as well as streamlines of the in-plane velocity. It can be observed that the modes come in pairs: 1–2 and 3–4, as is usual for POD modes and their time coefficients in a convective flow. The first coherent structure extracted by the POD is formed by modes 1 and 2 and constitutes a damped wave-like structure that is convected by the mean flow (see figure 5 for the spatial structure and figure 6 for the corresponding time coefficients; the supplementary video movie2_modes0-2.mov shows their behaviour in time). This is not a travelling wave such as those observed in transitional flows, as those of the examples mentioned in the introduction, but a coherent structure extracted by POD from a developed turbulent background that persists in fully developed turbulence, and is just a regular component of the flow on which irregular turbulent fluctuations are superimposed (see, e.g. Manhart & Wengle Reference Manhart and Wengle1993, for a similar case). Nevertheless, the present wave-like structure could be a surviving remnant of pre-existing, purely time-periodic, flow structures formed in the bent section and arising in the process of transition to turbulence past bends (see, e.g. the case of the flow past a circular cylinder by Sipp & Lebedev (Reference Sipp and Lebedev2007)). It was found that the first instability of the flow inside of a toroidal pipe is characterised by the appearance of travelling waves (Kühnen et al. Reference Kühnen, Holzner, Hof and Kuhlmann2014; Canton et al. Reference Canton, Schlatter and Örlü2016). It is therefore possible that similar waves appear in the transition to turbulence of the present flow case, and continue to modulate the large scales of the flow at high Reynolds numbers while being submerged in small-scale turbulence. To support this hypothesis, the frequencies and wavelengths of the present coherent structures are in the same range as those measured in toroidal pipes (Canton et al. Reference Canton, Schlatter and Örlü2016), and Brücker (Reference Brücker1998) observed swirl-switching even for $Re$ as low as 2000, although the measured oscillations had very low amplitude.

The present modes are, obviously, not strictly periodic in space nor in time: as can be seen in figure 5(b) showing the swirl intensity, the intensity of the modes is essentially zero upstream of the bend ( $s<0$ ), reaches a maximum at approximately $1D$ downstream of the bend end, and then decreases with the distance from the bend. Furthermore, the respective time coefficients are only quasi-periodic, as can be observed from their temporal signal, depicted in figure 6(a), and by their frequency spectra, figure 7(a). Nevertheless, it can be observed in figure 5 that the spatial structure of these modes is qualitatively sinusoidal along the streamwise direction $s_{o}$ , with a wavelength of approximately 7 pipe diameters. The figures in Brücker (Reference Brücker1998) actually already suggest the appearance of a wave-like structure in the presence of swirl-switching.

This wave-like structure is formed by two counter-rotating swirls, visible in the 2-D cross-sections in figure 5, which are advected in the streamwise direction while decaying in intensity and, at the same time, move from the inside of the bend towards the outside, as can be seen in the longitudinal cuts in figure 5 and in the supplementary video movie2_modes0-2.mov. The temporal amplitude of these modes is also qualitatively cyclic, as illustrated by the projection along the time coefficients in figure 6(b). The wave-like behaviour can be appreciated even better in the aforementioned video showing the flow reconstructed with these two modes, movie2_modes0-2.mov. Modes 1 and 2 are phase shifted by a quarter of their quasi-period: figure 5 shows that the structure of mode 2 is located approximately a quarter of a wavelength further downstream of the structure of mode 1; while, figure 6 illustrates the constant delay of the time coefficient of mode 2 with respect to that of mode 1.

The second structure, formed by modes 3 and 4, has a spatial layout that closely resembles that of the first pair, i.e. it is also a wave-like structure, and constitutes the first ‘harmonic’ of the wave formed by modes 1 and 2. The spatial structure of modes 3 and 4 has half of the main wavelength of modes 1 and 2, and the highest peak in the spectrum of the third and fourth time coefficients is at exactly twice the frequency of the peaks of $a_{1}$ and $a_{2}$ , as can be seen in figure 7(a). The video movie3_modes0-4.mov shows the reconstruction of the flow field by including modes 3 and 4. It can be observed that these modes introduce oscillations with higher frequency and smaller amplitude when compared to the reconstruction employing only modes 1 and 2.

As can be observed from figure 5, the modes do not present any connection to the straight pipe section preceding the bend. This is in direct contrast with the findings of Carlsson et al. (Reference Carlsson, Alenius and Fuchs2015), whose results were likely altered by the interference of an intrinsic frequency and wavelength on the recycling inflow boundary with the structure of the swirl-switching. Our results are, instead, in agreement with Noorani & Schlatter (Reference Noorani and Schlatter2016) who observed swirl-switching in a toroidal pipe (i.e. in the absence of a straight upstream section) confirming that these large-scale oscillations are not caused by structures formed in the straight pipe, but by an effect which is intrinsic to the bent section.

The power spectral density analysis of the time coefficients of these modes, computed as a Welch’s estimate, is presented in figure 7. One can see that, unlike the PSD of the two-dimensional POD modes (figure 4), the three-dimensional modes present two distinct peaks, one per pair of modes. The peak for the first modal pair is located at $St\approx 0.16$ , which is in the range of Strouhal numbers found by both Brücker (Reference Brücker1998) and Hellström et al. (Reference Hellström, Zlatinov, Cao and Smits2013). More importantly, this frequency is the lowest for this pair of modes and matches that given by the wavelength and propagation speed of the wave as well as that of the swirl-switching, as observed by reconstructing the flow field with the most energetic POD modes (see the online movies).

The present analyses were also performed on a pipe with curvature $\unicode[STIX]{x1D6FF}=0.1$ for the same Reynolds number. Swirl-switching was observed in this case as well, with dynamics which is qualitatively identical to the one observed for $\unicode[STIX]{x1D6FF}=0.3$ , but is characterised by lower frequencies, peaking at $St\approx 0.045$ . The lower frequencies and larger scales (wavelength of approximately $20D$ ) characterising the wave-like structure at this curvature meant that a quantitative analysis was too expensive with the present set-up. We have therefore limited this work to the study of one curvature only, but preliminary, not converged results can be found in Hufnagel (Reference Hufnagel2016).