1. Introduction
Flow through curved piping networks is ever-present with implications throughout engineering and physiology (Berger, Talbot & Yao Reference Berger, Talbot and Yao1983). Significant to the flow through curved pipes is the presence of secondary motions in the form of pairs of counter-rotating vortices (Kühnen et al. Reference Kühnen, Holzner, Hof and Kuhlmann2014; Lupi, Canton & Schlatter Reference Lupi, Canton and Schlatter2020) arising due to centrifugal forces. The introduction of simple curvature introduces new structures that are not present in the straight pipe, leading to different instability mechanisms (Canton, Schlatter & Örlü Reference Canton, Schlatter and Örlü2016; Gelfgat Reference Gelfgat2020; Lupi et al. Reference Lupi, Canton and Schlatter2020). Understanding these mechanisms as well as the nature of the secondary structures is important for comprehending such flows and manipulating their effects (Berger et al. Reference Berger, Talbot and Yao1983; Canton et al. Reference Canton, Schlatter and Örlü2016).
Transition to turbulence in unidirectional pipe flow is a nonlinear problem and is known to be linearly stable at all Reynolds numbers (Darbyshire & Mullin Reference Darbyshire and Mullin1995; Grossmann Reference Grossmann2000; Drazin & Reid Reference Drazin and Reid2004). The onset of instabilities in unidirectional straight pipe flow was first investigated by Reynolds (Reference Reynolds1883) and has been studied extensively since. The critical Reynolds number at which transition from the laminar regime occurs, and is sustained, is reported as
${Re} \approx 2300$
(Kerswell Reference Kerswell2005), or
${Re} \approx 2024$
by Avila et al. (Reference Avila, Moxey, de Lozar, Avila, Barkley and Hof2011) where a cross-over between transient and sustained turbulence is defined. The transition point can be significantly higher or lower depending on the inlet conditions (Eckhardt Reference Eckhardt2007). In the study by Darbyshire & Mullin (Reference Darbyshire and Mullin1995), turbulent growth was observed at
${Re} \gt 1800$
with sufficiently sized inlet perturbations. Below
${Re} \approx 1760$
, the same study does not report the maintenance of any growth despite attempts to disrupt the laminar state by injecting fluid into the flow and stirring the supply tank. The onset of instabilities in unidirectional straight pipe flow is highly sensitive to initial conditions and under controlled experiments, the transition has been delayed to a Reynolds number of
$10^5$
(Pfenninger Reference Pfenninger1961). In practice, this is unlikely but it highlights the dependence of transition on external disturbances, indicating an increasing critical Reynolds number with a reduction in finite-amplitude perturbations (Kerswell Reference Kerswell2005).
Unidirectional, Newtonian flow through straight pipes is governed by the Reynolds number only. Reciprocating straight pipe flows also require consideration of the frequency and amplitude of the oscillation (Xu et al. Reference Xu, Warnecke, Song, Ma and Hof2017). For reciprocating flows, the experimental study by Hino, Sawamoto & Takasu (Reference Hino, Sawamoto and Takasu1976) reported the nature of turbulent onset in the Reynolds number –Stokes parameter space. Their study, which defined the Reynolds number in terms of the Stokes layer thickness,
${Re}_{\delta _\omega } = U \delta _\omega / \nu$
, where
$U$
is the velocity amplitude,
$\delta _\omega = \sqrt {2\nu /\omega }$
is the Stokes layer thickness,
$\omega$
is the angular reciprocating frequency and
$\nu$
is the fluid kinematic viscosity, and the Stokes parameter,
$\lambda = d/(2\delta _\omega)$
, where
$d$
is the pipe diameter, found three regimes: weakly, conditionally and fully turbulent (Hino et al. Reference Hino, Sawamoto and Takasu1976). In the conditionally turbulent regime growth was observed during the deceleration phase of the cycle before it decayed back to a laminar state during the acceleration phase (Hino et al. Reference Hino, Sawamoto and Takasu1976). In this regime the critical Reynolds number was
${Re}_{\delta _\omega } \simeq 550$
for a Stokes parameter
$\lambda \gt 1.6$
(Hino et al. Reference Hino, Sawamoto and Takasu1976). The study by Xu et al. (Reference Xu, Warnecke, Song, Ma and Hof2017) also reported similar regimes by which pulsating pipe flows could be categorised. In the high-frequency regime for the Womersley parameter,
$\alpha = d \sqrt {2 \pi f / \nu }/2 \gt 12$
, where
$f$
is the pulsation frequency, with pulsation amplitude,
$A = U_o / U_s = 0.4$
, where
$U_o$
is the oscillating component and
$U_s$
is the steady component of the velocity, Xu et al. (Reference Xu, Warnecke, Song, Ma and Hof2017) reported that the transition threshold appeared agnostic to the frequency of the flow due to the sustainment of turbulent puffs appearing identical to those observed in steady flow experiments. At low frequencies (
$\alpha \lesssim 2.5$
), Xu et al. (Reference Xu, Warnecke, Song, Ma and Hof2017) reported the transition could be approximated by the quasi-steady assumption. The experimental study by Trip et al. (Reference Trip, Kuik, Westerweel and Poelma2012) on pulsatile pipe flow also found little influence from the Womersley parameter on transition in the high-frequency regime (
$\alpha \gt 10$
), with the nature of the turbulent structure such that growth was also observed during the deceleration phase and re-laminarisation occurring during the acceleration phase (Trip et al. Reference Trip, Kuik, Westerweel and Poelma2012).
In the straight pipe limit, the flow is stable to infinitesimal perturbations requiring finite-amplitude perturbations to initiate turbulence when the Reynolds number exceeds a critical threshold (Grossmann Reference Grossmann2000; Avila et al. Reference Avila, Moxey, de Lozar, Avila, Barkley and Hof2011). As the curvature increases, the base flow is fundamentally altered due to the appearance of secondary Dean flows (Dean Reference Dean1928a , Reference Deanb ), and subsequently, finite-amplitude perturbations are no longer required to initiate a change to the laminar state. When the effects of the curvature are large, the laminar flow becomes susceptible to linear transition mechanisms when the Reynolds number is above the critical threshold.
In constant curvature pipes, secondary vortices are present in the basic flow unlike the straight pipe. The vortices arise due to a centrifugal imbalance associated with forcing the fluid through the curvature, appearing as two counter-rotating vortices (Dean Reference Dean1928a
,
Reference Deanb
) oriented directly above and below the plane of symmetry of the pipe with their axes aligned with the direction of the flow. Early studies by Eustice (Reference Eustice1910) showed experimentally the presence of secondary flows using die filaments in a series of curved tubes before Dean (Reference Dean1928b
) proposed a mathematical solution, limited to small curvatures, for the appearance of the vortices which were later given his name. Canton, Örlü & Schlatter (Reference Canton, Örlü and Schlatter2017) studied the steady flow characteristics through a torus and observed that, for any Reynolds number and curvature, there always existed secondary Dean vortices (Canton et al. Reference Canton, Örlü and Schlatter2017). Canton et al. (Reference Canton, Örlü and Schlatter2017) reported the features of the secondary flow as highly dependent on the curvature of the pipe, in particular the shape and location of the core of the vortices, and the position of the maximum streamwise velocity (Canton et al. Reference Canton, Örlü and Schlatter2017). Lupi et al. (Reference Lupi, Canton and Schlatter2020) studied the unidirectional flow through a pipe with localised curvature,
$\delta = 1/3$
, where
$\delta$
is the ratio of the radius of curvature of the centreline of the pipe to the radius of the pipe cross-section, with a
$90$
-degree bend. They observed two counter-rotating vortices downstream of the curved region (Lupi et al. Reference Lupi, Canton and Schlatter2020), similar to what is observed in the torus and helical pipe.
The instability mechanisms of the unidirectional flow through the closed-system torus and open helical pipes were experimentally studied by Kühnen et al. (Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) for low values of curvature. Kühnen et al. (Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) investigated the curvature range
$\delta \in [0.01, 0.1]$
reporting a cross-over from subcritical to supercritical transition mechanisms at
$\delta = 0.028$
(Kühnen et al. Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015). Above curvatures of
$0.028$
, Kühnen et al. (2015) observed a supercritical mechanism. Below this, the transition mechanism was observed as subcritical (Kühnen et al. Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015), similar to what occurs in the straight pipe.
Canton et al. (Reference Canton, Schlatter and Örlü2016) extended on the work of Kühnen et al. (Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) in the low-curvature range with a numerical study covering the entire curvature range,
$\delta \in [0.002, 1.0]$
, in the closed-system torus. They found all curvatures to become unstable to a linear mechanism at Reynolds numbers, based on the diameter of the pipe, in the range
${Re} \in [\approx 2000, \approx 6500]$
(Canton et al. Reference Canton, Schlatter and Örlü2016). Lupi et al. (Reference Lupi, Canton and Schlatter2020) also reported a linear instability mechanism in the unidirectional flow through the bent pipe at
${Re} \approx 2531$
. In both studies on the torus and bent pipe, the linear stability analyses were able to predict the transition closely, matching with the observations from direct numerical simulation and confirming a linear instability mechanism in both geometries.
Flows through more complex curved pipes were studied by Jacob et al. (Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023). They investigated the turbulent onset of oscillatory flows in both single and double bifurcating tubes; the same geometry that is used in this study, shown in figure 1. Similar to the flow through the torus, helical and bent pipes, Jacob et al. (Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023) also reported the appearance of vortices as counter-rotating pairs downstream of the bifurcating joints, where one pipe transitions into two. In the opposite direction of flow, when two pipes merge into one, two counter-rotating vortex pairs were observed. Jacob et al. (Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023) also observed the appearance of conditional turbulence at peak flow rates, in contrast to the straight pipe where conditional turbulence is reported during the peak deceleration phase of the cycle. Jacob et al. (Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023) suggested the turbulence grows on the unstable vortices. The frequencies at which the turbulent growth appeared were orders of magnitude higher than the oscillatory frequency of the flow (Jacob, Tingay & Leontini Reference Jacob, Tingay and Leontini2021).
Figure 1. Single-bifurcation
$(1:2)$
and double-bifurcation
$(1:2:4)$
geometries used in this study. The lengths of the primary, secondary (single-bifurcation only) and tertiary pipes (double-bifurcation only) are not shown to scale. The red markers indicate the origin (
$x, y, z = 0$
). The blue markers indicate the primary pipe free ends. The yellow markers indicate the secondary pipe free ends (single-bifurcation only). The green markers indicate the tertiary pipe free ends (double-bifurcation only).
It is clear that the instability mechanisms in straight pipes and pipes with curvature are different. The straight pipe loses stability due to a nonlinear process in contrast to pipes with curvature which become unstable due to a linear mechanism associated with the secondary Dean vortices. A quantitative understanding of the onset of instabilities, and subsequent turbulence, in complex curved pipe networks, constant curvature pipes and the bent pipe has implications for vascular and respiratory flows. Previous studies investigating flows through helical and bent pipes as well as the torus have been motivated by understanding the effects of curvature on transition (e.g. Canton et al. Reference Canton, Schlatter and Örlü2016; Kühnen et al. Reference Kühnen, Braunshier, Schwegel, Kuhlmann and Hof2015) or to obtain accurate friction factor data from pressure drop measurements (e.g. Cioncolini & Santini Reference Cioncolini and Santini2006; Kühnen et al. Reference Kühnen, Holzner, Hof and Kuhlmann2014). We are motivated by trying to understand the fundamental transition processes occurring in airway-like geometries. Even though respiratory flows are typically reciprocating, previous work (Jacob et al. Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023) has indicated that the turbulent transition occurs at frequencies much faster than any reciprocation with little impact from the history of the flow, even for conditions relevant to high-frequency medical ventilation (Hibberd et al. Reference Hibberd, Leontini, Scott, Pillow, Miedema, Rimensberger and Tingay2024; Scott et al. Reference Scott, Jacob, Tingay and Leontini2024), suggesting a study of the unidirectional flow is warranted.
This study investigates the linear instabilities of the unidirectional steady flows through the single and double bifurcations using linear stability analysis and direct numerical simulation. We find that the critical Reynolds number, structure, wavelength and frequency of all modes depend on the formation of the secondary Dean flows which themselves are dependent on the flow direction through, and configuration of (i.e. the number of generations) the geometry. We start by investigating the formation of the secondary Dean vortices in the basic flow which drive the instabilities and then conduct linear stability analyses, in two geometries, of the steady flows in both directions through each domain. Our results indicate that the linear stability analyses accurately predict the transition Reynolds number and the frequencies at which the leading linear modes oscillate. The direct numerical simulation (DNS) confirms the linear stability analyses to accurately predict the dynamics of the transition from laminar to time-dependent flows and the linear stability analyses confirms that the onset of turbulence in the flows through complex curved pipe networks is generated by linear instability mechanisms of the unstable Dean vortices.
Note that throughout this paper the term ‘bifurcation’ has two meanings. The primary meaning of the term describes one pipe splitting into two. The secondary meaning of the term relates to dynamical systems and is used to describe a change in behaviour, i.e. a spatial symmetry-breaking bifurcation (SSBB) or change in flow state. To avoid confusion, we use the acronym SSBB to describe the dynamical systems behaviour where applicable.
The remainder of this paper is organised as follows: § 2 describes the details of the geometries, the extent of the study and the computational set-up for performing DNS, calculating the base flows and performing the linear stability analyses. Section 3 details the characteristics of the base flows and the results from the linear stability. Section 4 compares the results from the linear stability with DNSs. Section 5 briefly discusses the relevance of this research for respiratory flows and the limitations of the study, and Section 6 provides a short discussion of the results and concluding remarks.
2. Methods
2.1. Fluid domain
The two geometries – a single and a double bifurcation – used in this study were developed based on the approximately self-similar dimensions of the human airway detailed in Grotberg (Reference Grotberg1994). The planar-symmetric model is based on the fundamental self-similar features that are observed in the first sixteen generations of the human respiratory system where a generation represents a bifurcation of an airway vessel into two smaller vessels. Each generation, or bifurcation, is defined by its primary vessel or pipe diameter, its secondary pipe diameter, the angle between the two secondary pipe centrelines and the length of each straight pipe section. In this study, these relationships are as follows:
$L_n/D_n=3.5$
for pipe lengths without free ends,
$D_{n+1}/D_n=0.79$
,
$\phi =64^{\circ }$
and
$R_n=5D_n$
, where
$L_n$
is the respective pipe length,
$D_n$
is the primary pipe diameter,
$\phi$
is the angle between the secondary pipe centrelines,
$R_n$
is the radius of curvature of the pipe centreline and
$n$
is the generation index. Figure 1 details a plan view of both geometries used in this study. From top to bottom, to establish definitions, throughout this study we refer to the top pipe as the primary pipe or generation. Where the primary pipe splits into two, these are referred to as the secondary pipes and, in the double-bifurcation model where the secondary pipes split again, these are referred to as the tertiary pipes. We have used the same index,
$n$
, as detailed in Grotberg (Reference Grotberg1994), where
$1$
, and
$2$
refer to the primary and secondary generations, respectively. For this study we have chosen an angle of
$64^{\circ }$
for all bifurcating joints to maintain geometric similarity with previous studies (Jacob et al. Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023; Wanigasekara et al. Reference Wanigasekara, Jacob, Manasseh and Leontini2024). We have extended the entry and exit pipe lengths to remove the effects imposed by the boundaries. The relationship between the primary and secondary generations for both geometries used in this study are detailed in table 1.
Table 1. Geometric relationships between generations for the bifurcation models used in this study.
The planar geometries lie in the
$x{-}z$
plane and the primary pipe cross-sections are parallel to the
$x{-}y$
plane aligning the primary pipes with the
$z$
-axis of the global coordinate system. Each bifurcation, where pipe
$n$
splits into pipe
$n+1$
, transitions through a large curved region of the domain following the centrelines shown in figure 1, noting that the centrelines run through the centre of the pipe sections. All bifurcating transitions feature a median value of curvature,
$\delta = a/R_n = 0.0895$
, where
$a$
is the pipe radius equal to
$(D_n/2+D_{n+1}/2)/2$
.
Throughout this study we refer to the axial and transverse components of velocity, designated as
$U_{ax}$
and
$U_{tr}$
, respectively. The definitions of these terms are local in nature. The axial component of a velocity vector, at any point in the domain, is the component of the vector that is tangent to the local pipe centreline. The transverse component of a velocity vector, at any point in the domain, is the component of the vector that is projected on to a plane which lies normal to the local pipe centreline. The origin of both geometries where
$x, y, z = 0$
, is located where the straight pipe section,
$L_1$
, terminates. The origins are shown as red markers in figure 1.
2.2. Boundary conditions
At the free ends of the primary pipes (indicated by blue markers in figure 1) a parabolic Hagen–Poiseuille profile was applied
\begin{equation} \boldsymbol {u}_z(r)=U^*\left [1-\frac {4r^2}{D_1^2}\right]\!, \end{equation}
where
$U^*$
is the peak centreline velocity at the primary pipe free end and
$r$
is the radial coordinate of the pipe. At all rigid wall boundaries, zero relative motion was implemented using a no-slip Dirichlet condition, i.e.
$\boldsymbol {u} = 0$
, and the pressure gradient normal to the wall (
$\partial p / \partial \boldsymbol {n}$
) is set such that mass conservation is enforced at the wall where
$\boldsymbol {u}$
is the non-dimensional velocity,
$p$
is the non-dimensional pressure and
$\boldsymbol {n}$
is the vector normal to the wall. The non-dimensionalisation is detailed in § 2.3. At the free ends of the secondary pipes (in the single-bifurcation geometry indicated by yellow markers in figure 1) and tertiary pipes (in the double-bifurcation geometry indicated by green markers in figure 1) a zero pressure open Dirichlet condition was applied for the pressure, i.e.
$p = 0$
, and zero normal velocity gradient,
$\partial \boldsymbol {u} / \partial \boldsymbol {n} = 0$
. For the direction of flow from the tertiary or secondary pipes to the primary pipe, the profile defined in (2.1) was reversed at the primary pipe free end. For this direction of flow, the boundary condition for the velocity at the secondary or tertiary pipe free ends was modified to preserve numerical stability. The modified boundary condition adds a positive, mean velocity divergence to the single layer of elements at the exits of the domain. The strength of the spatially varying artificial divergence,
$D_f(x) = C_f[1-(x_n/E_t)^2]$
, eventually reduces to zero at some location downstream of the boundaries. Here,
$C_f$
is a user-defined constant controlling the strength of the artificial divergence,
$x_n$
is the normal distance to the boundary and
$E_t$
is the thickness of the element layer at the boundary (Fischer et al. Reference Fischer, Loth, Lee, Lee, Smith and Bassiouny2007). The modified boundary condition enforces some local divergence while maintaining it globally. The modification, detailed in the study by Fischer et al. (Reference Fischer, Loth, Lee, Lee, Smith and Bassiouny2007), was used with the same source code used in this study for high-Reynolds-number vascular flows where the effects were shown to be contained near the boundary. Other studies using the same geometry and source code to that used here reported the impact of the modified condition to be contained to regions within
$3{-}4$
diameters of the boundary (Jacob et al. Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023), well clear of the region in the domain of interest.
2.3. Governing equations and solver
To calculate the steady-state, three-dimensional base flows we solve the non-dimensional, incompressible Navier–Stokes equations
\begin{equation} \left . \begin{array}{ll} \displaystyle \frac {\partial \boldsymbol {u}}{\partial \tau } + (\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u} = - \boldsymbol{\nabla }p + \frac {1}{{Re}} {\nabla} ^2 \boldsymbol {u} + \boldsymbol {f}, \\[8pt] \displaystyle \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u} = 0, \end{array}\right \} \end{equation}
where
$\boldsymbol {u} = \boldsymbol {u}^* / U^*$
is the non-dimensional velocity field and
$\boldsymbol {u}^*$
is the velocity,
$p = p^* / (\rho U^{*2})$
is the non-dimensional pressure,
$p^*$
is the pressure and
$\rho$
is the fluid density,
$\tau = t^* U^* / D_1$
is the non-dimensional time,
$D_1$
is the characteristic length equal to the primary pipe diameter (detailed in § 2.1) and
$t^*$
is the time. The Reynolds number is
${Re} = U^* D_1 / \nu$
, where
$\nu$
is the kinematic viscosity. The acceleration term,
$\boldsymbol {f}$
, is a force introduced to implement the selective frequency damping framework detailed in § 2.5 which was used to find the steady base flows. For the DNSs,
$\boldsymbol {f} = 0$
. Aside from the fluid domain and the direction of flow, the Reynolds number is the only parameter which governs the flow conditions.
The linear stability analysis requires the calculation of the linear evolution of perturbations of the steady base flow. By assuming the state vectors can be decomposed into a base flow and perturbation, i.e.
$\boldsymbol {u} = \boldsymbol{U} + \boldsymbol {u}^\prime$
and
$p = P + p^\prime$
, where
$\boldsymbol{U}$
is the velocity base flow,
$\boldsymbol {u}^\prime$
is the velocity perturbation,
$P$
is the pressure base flow and
$p^\prime$
is the pressure perturbation, equations for the evolution of the perturbations are formed as
\begin{equation} \left . \begin{array}{ll} \displaystyle \frac {\partial \boldsymbol {u}^\prime }{\partial \tau } + (\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^\prime + (\boldsymbol {u}^\prime \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U} = - \boldsymbol{\nabla }p^\prime + \frac {1}{{Re}} {\nabla}^{2} \boldsymbol {u}^\prime , \\[8pt] \displaystyle \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u}^\prime = 0. \end{array}\right \} \end{equation}
The equations for the flow (2.2) and the perturbations (2.3) were solved using the spectral element method, detailed by Patera (Reference Patera1984), with the open-source code, Nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008). The mesh for constructing the fluid domain consisted of
$20\,800$
hexahedral elements for the cases investigating the single-bifurcation geometry, and
$49\,920$
hexahedral elements for the cases investigating the double-bifurcation geometry. All simulations used twelfth-order tensor-product Lagrange polynomials, in a Galerkin formulation, using Gauss–Legendre–Lobatto quadrature points for calculating the integrals of the equations in weak form.
The temporal domain was discretised using three-way time splitting to solve for the terms of the Navier–Stokes equations. A third-order Adams–Bashforth scheme was used to establish the advective substep, the Crank–Nicholson scheme was used for the diffusive substep (Fischer Reference Fischer2003) and the pressure correction was achieved by forming a Poisson equation, formed by taking the divergence of the second substep equation and enforcing the divergence-free condition at the end of this substep.
Grid resolution studies were performed on the base flow and perturbation calculations to confirm convergence of global and local parameters by increasing the interpolation polynomial order of each spectral element (
$p$
-refinement). These studies are detailed in Appendix A. Subsequently, all simulations (base flow and perturbation) were run with total node counts of
$27\,684\,800$
and
$66\,443\,520$
for the single- and double-bifurcation cases, respectively. The complete system was computed on parallel architecture across typically
$2048$
CPUs requiring a total wall time of
${\approx}\, 10$
days.
2.4. Extent of the study
This study considers four cases across two geometries. We refer to inflow as the flow direction from the primary pipe to the secondary and tertiary pipes. Outflow is the opposite of this, with the flow direction from the secondary or tertiary pipes towards the primary pipe. The flow directions are shown in figure 1. For simplicity, we refer to the inflow cases as
$12in$
and
$124in$
for the single- and double-bifurcation geometries, respectively. We refer to the outflow cases as
$12out$
and
$124out$
for the single- and double-bifurcation geometries, respectively. These are summarised in table 2.
Table 2. Nomenclature for the four cases studied and their definition.
2.5. Direct numerical simulation
To establish the base flows for the linear stability analysis each case was initiated with a DNS, driven by a parabolic Hagen–Poiseuille boundary condition (2.1) at the free end of the primary pipe. Each case was run for
$\tau \ge 10^3$
time units ensuring a statistical steady state was achieved using a variable time step no greater than
$10^{-3}$
time units. The ratio of the run time to achieve a statistical steady state to the advective time scale,
$\tau _{adv} = L_{FD} / U^*$
, in the double-bifurcation geometry, where
$L_{FD}$
is the length of the fluid domain and is equal to
$\tau /\tau _{adv} = 1000/L_{FD} / U^* \approx 15$
, confirming sufficient time for any transients to wash through the domain.
For the linear stability analysis to proceed, a steady base state or base flow is required about which the governing equations can be linearised (Jordi et al. Reference Jordi, Cotter and Sherwin2014, Reference Jordi, Cotter and Sherwin2015). Steady solutions were found using the selective frequency damping (SFD) framework developed by Åkervik et al. (Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006) that applies a forcing term to the momentum equation, expressed as
\begin{equation} \dot {\boldsymbol{q}} = \mathcal{N}(\boldsymbol{q}) - \chi (\boldsymbol{q} - \boldsymbol{w}),\end{equation}
where
\begin{equation} \dot {\boldsymbol{q}} = \mathcal{N}(\boldsymbol{q}) \end{equation}
is the Navier–Stokes equation expressed in operator form, with
$\boldsymbol{q} = [\boldsymbol {u}\; p]^T$
representing the state vector consisting of the velocity and the pressure, and
\begin{equation} - \chi (\boldsymbol{q} - \boldsymbol{w}) = \boldsymbol {f}, \end{equation}
which is the forcing term. This form of the forcing term implements a first-order Butterworth low-pass causal filter, applied to the momentum equation (2.2). Here,
$\chi$
is the filter gain,
$\boldsymbol{q}$
is the flow field solution and
$\boldsymbol{w}$
is the low-pass filtered solution, expressed as
$\boldsymbol{w} = T \times \boldsymbol{q}$
, where
$T$
is a temporal filter kernel with filter width,
$\varDelta$
. When expressed in differential form, detailed in Pruett et al. (Reference Pruett, Gatski, Grosch and Thacker2003), the filter is defined as
\begin{equation} \frac {\partial \boldsymbol{w}}{\partial \tau } = \frac {\boldsymbol{q} - \boldsymbol{w}}{\varDelta }. \end{equation}
The forcing term drives the flow towards the steady-state solution defined as a modification of
$\boldsymbol{q}$
with reduced temporal oscillations (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006). The filter parameters,
$\chi$
and
$\varDelta$
, were selected based on a priori knowledge of the statistically stable flows. The filter width was selected to suppress the lowest unstable oscillations, observed from the DNS. The filter width relates to the cutoff frequency,
$\omega _c$
, from the relation
$\varDelta = 1 / (\omega _c 2 \pi)$
. The gain (sometimes referred to as the control coefficient) was set as low as possible, with typical starting values of
$\chi = 0.2$
. This enabled the system to converge to the steady state in the least amount of time, avoiding an over damped and slowly evolving system. The SFD was switched on for times
$\tau \geq 10^3$
, and generally required around
$1.5 \times 10^3$
time units to reach the converged steady state. The base flow was considered converged when the global
$L_2$
norm reached
${\leq} 10^{-9}$
. This criterion was established from a test of the SFD framework detailed in Appendix B.
2.6. Linear stability analysis
The linear modes of each base flow were calculated using an Arnoldi iteration technique (Arnoldi Reference Arnoldi1951). The technique is well established and only an outline of the main points is provided here. If the equations of motion for the perturbation (equations (2.3)) are expressed in operator form
\begin{equation} \dot {\boldsymbol{q}}^\prime = \mathcal{L}\boldsymbol{q}^\prime , \end{equation}
where
$\boldsymbol{q}^\prime = [\boldsymbol {u}^\prime \; p^\prime]^T$
, the aim is to find the eigenvalues and the eigenvectors (or linear modes) of
$\mathcal{L}$
. Rather than solving for the eigenvalues of this full system, it can be related to the linear map
\begin{equation} \boldsymbol{q}^{\prime }(\tau +T) = L \boldsymbol{q}^{\prime }(\tau), \end{equation}
which states that the disturbance,
$\boldsymbol{q}^{\prime }$
, at time,
$\tau +T$
, is equal to the action of the linear operator,
$L$
, acting on the current state of the flow,
$\boldsymbol{q}^{\prime }$
, at time,
$\tau$
. In this expression
$T$
is the time between perturbation snapshot vectors that are formed by the action of the linear operator (i.e. integrating the linearised Navier–Stokes equations (2.3) forward from time,
$\tau$
to
$\tau +T$
) which constructs the Krylov subspace. The formation of the complete matrix, and solving for the eigenvalues and eigenvectors, is computationally difficult and not required. Using the Krylov method, the eigenvalues and eigenvectors of the lower-dimensional subspace (strictly referred to as Ritz pairs) are accurate approximations of the eigenvalues and eigenvectors of the full matrix.
The Krylov sequence of the operator
$L$
is expressed as
\begin{equation} \mathcal{K} = [\boldsymbol{q}, L \boldsymbol{q}, L^2 \boldsymbol{q}, L^3 \boldsymbol{q},\ldots]. \end{equation}
Our iterative Arnoldi technique constructs an orthonormal basis,
$\boldsymbol{Q}$
, and a truncated upper-Hessenberg matrix
$\boldsymbol{H}$
from the Krylov matrix
$\mathcal{K}$
. The matrix
$\boldsymbol{H}$
is related to the action of the operator
$L$
via the similarity transformation
\begin{equation} L=\boldsymbol{QHQ}^{*} \quad \text{or} \quad L\!\boldsymbol{Q}=\boldsymbol{QH}, \end{equation}
where
$\boldsymbol{Q}^*$
is the conjugate transpose. The linear modes for each base flow reported in this study are constructed from the eigenvectors of the upper-Hessenberg matrix. The eigenvalues of
$\boldsymbol{H}$
are the multipliers
$\mu$
of
$L$
, which are related to the eigenvalues
$\sigma$
of the original system
$\mathcal{L}$
via
$\mu = e^{\sigma T}$
.
Specifically, the SFD-stabilised velocity field,
$\boldsymbol{U}$
, calculated from (2.2) was loaded as the base flow, for the time-stepping linearised Navier–Stokes solver in (2.3). Randomised low-amplitude noise of the order
$10^{-5}$
was set as the initial condition across the entire computational domain to initiate the growth of the perturbations for the velocity and the pressure fields. The perturbations were integrated forward in time to establish the action of the linear operator,
$L$
, which was captured across successive perturbation snapshot vectors that were output at intervals of
$T = 0.1$
time units to establish the Krylov subspace.
The algorithm used to construct the orthonormal basis and upper-Hessenberg approximation and to solve for the eigenvalues and eigenvectors of the linear modes was executed outside of the time-stepping linearised Navier–Stokes solver, similar to the technique outlined by Mamun & Tuckerman (Reference Mamun and Tuckerman1995). While this does not allow the common implicitly restarted Arnoldi method to be used, it allows for the size of the Krylov sequence to be adjusted after the calculation of the snapshots which can be beneficial. A Krylov subspace of typically
$10{-}14$
vectors was used to calculate the linear modes for each base flow by iterating through a much larger vector ensemble of
${\approx}150$
, i.e. we calculate the eigenvalues and eigenvectors associated with the Krylov subspace spanned by the vectors
$1$
to
$14$
, then
$2$
to
$15$
, through to
$137$
to
$150$
to identify the consistently recurring linear modes.
Modes with multipliers (eigenvalues of
$L$
) such that
$|\mu | \gt 1$
grow over time and therefore are unstable. If the mode has an associated time scale, i.e. there is a imaginary component of the eigenvalue, the frequency at which the mode oscillates is established from the relation
\begin{equation} f=\frac {1}{2 \pi T}\tan ^{-1}\left (\frac {\mu _i}{\mu _r}\right)\!, \end{equation}
where
$\mu _i$
is the imaginary component and
$\mu _r$
is the real component of the eigenvalue. If the mode is purely real, i.e. there is no imaginary component associated with the eigenvalue, then there is no periodic behaviour associated with the growth of the mode.
The in-house Arnoldi algorithm used in this study has been verified against the work by Lupi et al. (Reference Lupi, Canton and Schlatter2020) in the bent pipe with excellent results. The verification study is detailed in Appendix C.
3. Results
3.1. Observations of flow transitions and states from direct numerical simulations
We start by investigating the characteristics of the base flows. In all four cases we ran the DNS computations at increments in Reynolds number of
$250$
to isolate where the flows transition from a steady to a time-dependent regime. All DNS computations were run for
$\tau \geq 1000$
time units to ensure a statistical steady state was reached (detailed in § 2.5). From the DNS data we have plotted the standard deviation of the velocity magnitude,
$\boldsymbol {u}_{mag} = \sqrt {\boldsymbol {u}_x^2 + \boldsymbol {u}_y^2 + \boldsymbol {u}_z^2}$
, as a function of the Reynolds number, shown in figure 2 where
$\boldsymbol {u}_x$
,
$\boldsymbol {u}_y$
and
$\boldsymbol {u}_z$
are the velocity components in the
$x$
,
$y$
and
$z$
directions, respectively. From this relationship we have fitted a power-law curve and extrapolated the curve to the zero crossing to approximate only the transition Reynolds number for each case, also shown in figure 2 and listed in table 3. Note in all plots we have used three data points, except for the case of the double-bifurcation outflow,
$124out$
, where we have used five data points due to the outlier at
${Re} = 6500$
.
Figure 2. Results from DNS computations to approximate where each case transitions from steady to time-dependent. Each panel displays the standard deviation ((SD), red triangles) calculated from DNS of the velocity magnitude (
$y$
-axis) plotted as a function of
${Re}$
(
$x$
-axis), a power-law curve fit (red dashed lines) and an interpolated value of
${Re}$
(black filled circles) calculated from the power-law curve fit. Panel (a) is from the single-bifurcation inflow. The probe point is located
${\approx}5.0D_1$
along the
$z$
-axis, from the origin, in the secondary pipe. Panel (b) is from the double-bifurcation inflow. The probe point is located
${\approx}\, 5.0D_1$
along the
$z$
-axis, from the origin, in the secondary pipe. Panel (c) is from the single-bifurcation outflow. The probe point is located
${\approx}-4.0D_1$
along the
$z$
-axis, from the origin, in the primary pipe. Panel (d) is from the double-bifurcation outflow. The probe point is located
${\approx}\, -6.0D_1$
along the
$z$
-axis, from the origin, in the primary pipe.
Table 3. Results from DNS computations to isolate where each case transitions from laminar to time-dependent and where symmetry breaks occur and are resurrected. The values of
${Re}$
were calculated from Figure 2. Each case number is detailed in Table 2.
In the double-bifurcation outflow study,
$124out$
, we observe two SSBBs in the temporally stable base flow that is unique to this case. The first SSBB occurs between
${Re} = 4250$
and
$4500$
where the base flow symmetry is altered. The symmetry is then partially recovered due to a second SSBB arising between
${Re} = 5750$
and
$6000$
. The range of
${Re}$
where the symmetry breaks occur are detailed in the third and fourth columns of table 3 and discussed in § 3.2.3.
3.2. Steady base flows obtained using selective frequency damping
Steady flows were obtained using SFD. For all cases, the SFD was started from a DNS case that was statistically stationary. We implemented this process regardless of whether the DNS case showed any temporal variation; this way we could guarantee that the SFD did not introduce any spurious flow features, and the eventual steady states found were reliable. In all cases, the stabilised base flow solution calculated with the SFD framework, at a value of Reynolds number just above the critical value, appeared visually indistinguishable to the base flow solution just below the critical value (obtained with the SFD framework switched on and the unforced naturally steady solution obtained with the SFD switched off). In all four comparisons the difference in the Reynolds number was
$250$
. This confirmed that the base flow solution at a Reynolds number just above the critical value was representative of the base flow solution at the critical Reynolds number.
3.2.1. Inflow topology
The stable base flows all exhibit pairs of counter-rotating Dean vortices in regions of the flow downstream of the curvature. For the inflow cases,
$12in$
and
$124in$
illustrated in figures 3 and 4, respectively, the vortices are observed as a counter-rotating pair that are mirror symmetric about the
$x{-}z$
plane. The single vortex pairs appear in both secondary pipes and are mirror symmetric about the
$y{-}z$
plane. In both cases, it is not until the flow enters the straight section of the secondary pipes that the vortices start to form clear structures. The double-bifurcation inflow case,
$124in$
, consists of two counter-rotating pairs (four vortices) appearing in the tertiary pipe, that are mirror symmetric about the
$x{-}z$
and
$y{-}z$
planes, before recovering the two-vortex system (one counter-rotating pair).
Figure 3 shows the structure of the stable base flow for the single-bifurcation inflow case,
$12in$
, at
${Re} = 7500$
. At pipe section
$3D$
, upon exiting the curved region of the geometry, the vortices start to form clear structures. The structure of the vortex pairs assumes a kidney-like profile where a deficit in the secondary flow is evident at the inside of the secondary pipe (left-hand side of the pipe section with respect to the image on the page). Between sections
$3D$
and
$4D$
the vortex centres start to shift towards the centre of the pipe as the deficit in the secondary flow is reduced. By pipe section
$5D$
the structure of the two vortices occupy almost the full pipe cross-section. Only small changes to the location of the vortex centres are observed as they move horizontally towards the inner pipe wall as the flow progresses further down the secondary pipe. Downstream of pipe section
$6D$
the secondary flow structure is fully formed. Two small vortices can be seen forming at approximately pipe section
$5D$
where a deficit in the secondary flow grows on the outside of the pipe, increasing in size as the primary vortices shift across the pipe section.
Figure 3. Single-bifurcation inflow case,
$12in$
, secondary pipe base flow topology. The pipe section views display the spatially evolving stable base flow at
${Re} = 7500$
. The local axial velocity (
$U_{ax}$
) is overlaid by the transverse velocity (
$U_{tr}$
) streamlines. The left-hand side of each pipe section view pertains to the inside of the secondary pipe from the image on the right-hand side. The direction of flow through the pipe section views is into the page. The image on the right-hand side shows the location of each pipe section in the secondary pipe. The vertical arrow indicates the direction of flow through the domain. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the image on the right-hand side.
Figure 4 shows the structure of the stable base flow for the double-bifurcation inflow case,
$124in$
, at
${Re} = 5000$
. The formation of the vortices in the secondary pipe (panel (a)) undergo a similar evolution to what is observed in the single-bifurcation case,
$12in$
, with the core of the vortices following an almost identical path as the flow progresses through the domain. In case
$124in$
, by pipe section
$5D$
the structure of the two vortices occupies almost the full pipe cross-section, and is fully formed downstream of pipe section
$6D$
. The spatially evolved base flow topologies for both inflow cases, shown in figures 3 and 4, appear qualitatively similar despite the difference in Reynolds number. Panel (b) of figure 4 shows the stable base flow in the tertiary pipe of the double-bifurcation inflow case,
$124in$
, at
${Re} = 5000$
. The vortices emerge from the second bifurcation as two counter-rotating pairs before eventually recovering the two-vortex system.
Figure 4. Double-bifurcation inflow case,
$124in$
, base flow topology. The pipe section views display the spatially evolving stable base flow at
${Re} = 5000$
. The local axial velocity
$(U_{ax})$
is overlaid by the transverse velocity
$(U_{tr})$
streamlines. The left-hand side of each pipe section view pertains to the inside of the pipe from the image on the right-hand side. The direction of flow through the pipe section views is into the page. The images on the right-hand side show the location of each pipe section in the secondary (a) and tertiary (b) pipes. The vertical arrows indicate the direction of flow through the domain. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the images on the right-hand side in both panels.
3.2.2. Outflow topology
When the direction of flow is reversed for the outflow cases,
$12out$
and
$124out$
illustrated in figures 5 and 6, respectively, the vortices are observed as two counter-rotating pairs downstream of the first bifurcation (the bifurcation separating the tertiary and secondary pipes in the case of
$124out$
). In the unique case of
$124out$
, in the primary pipe, four counter-rotating pairs are observed before eventually decaying back to a four-vortex system (two counter-rotating pairs). The structure of the four-vortex system observed in the primary pipe of the outflow case,
$12out$
, and the secondary pipe of case,
$124out$
, is characterised by symmetry about the vertical and horizontal planes running through the centre of the pipe sections (until pipe section
$3D$
, in
$124out$
, where the symmetry is lost about the horizontal plane). The topology of the flow in the secondary pipe, in case
$124out$
, also inherits symmetry about the
$y{-}z$
plane. However, the eight-vortex system, in the primary pipe of case
$124out$
, only displays symmetry about a horizontal plane running through the centre of the pipe section.
Figure 5 shows the structure of the stable base flow for the single-bifurcation outflow case,
$12out$
, at
${Re} = 4750$
. The structure of the vortices exhibit symmetry about the
$x{-}z$
and
$y{-}z$
planes. Upon exiting the bifurcation (section
$0D$
) the vortices are essentially established with only subtle changes in the location of the vortex centres, moving towards the centre of the pipe section, as the flow progresses downstream of the bifurcation.
Panel (b) of figure 6 shows the structure of the stable base flow in the secondary pipe of case,
$124out$
, at
${Re} = 6250$
which appears qualitatively similar to the features observed in figure 5, after pipe section
$3D$
.
Panel (a) of figure 6 shows the structure of the stable base flow for case,
$124out$
, at
${Re} = 6250$
in the primary pipe. Upon exiting the bifurcation, the secondary flow exhibits four significant counter-rotating vortex pairs orientated around the circumference of the pipe section. A symmetry is apparent about the
$x{-}z$
plane; however, the symmetry about the
$y{-}z$
plane is lost. This vortex structure is preserved until pipe section
$-2D$
, after which the two vortex pairs on the left-hand side begin to merge and relocate towards the centre of the pipe. The two vortex pairs on the opposite (right-hand) side survive until pipe section
$-5D$
, after which the four-vortex system is recovered. However, the symmetry about the
$y{-}z$
plane is never recovered.
Figure 5. Single-bifurcation outflow case,
$12out$
, primary pipe base flow topology. The pipe section views display the spatially evolving stable base flow at
${Re} = 4750$
. The local axial velocity
$(U_{ax})$
is overlaid by the transverse velocity
$(U_{tr})$
streamlines. The left-hand side of each pipe section view pertains to the left-hand side of the primary pipe from the image on the right-hand side. The direction of flow through the pipe section views is out of the page. The image on the right-hand side shows the location of each pipe section in the primary pipe. The vertical arrow indicates the direction of flow through the domain. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the image on the right-hand side.
Figure 6. Double-bifurcation outflow case,
$124out$
, base flow topology. The pipe section views display the spatially evolving stable base flow at
${Re} = 6250$
. The local axial velocity
$(U_{ax})$
is overlaid by the transverse velocity
$(U_{tr})$
streamlines. The left-hand side of each pipe section view pertains to the left-hand side (panel (a)) and inside (panel (b)) of the pipe from the image on the right-hand side. The direction of flow through the pipe section views is out of the page. The image on the right-hand side shows the location of each pipe section in the primary (top) and secondary (bottom) pipes. The vertical arrows indicate the direction of flow through the domain. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the images on the right-hand side in both panels.
3.2.3. Spatial symmetry-breaking bifurcation in the double-bifurcation outflow
In the (natural) temporally stable
$124out$
base flow we observe a loss of vortex symmetry in the primary pipe of the basic flow, shown in panel (a) of figure 6. With respect to the orientation of the pipe section views in the plane of the page, a symmetry about the horizontal plane is clear but the symmetry about the vertical plane is lost. This asymmetry in the temporally stable secondary flows is not a feature of all Reynolds numbers and we first detect this modification between
${Re} = 4250$
and
$4500$
, shown in figure 7. At
${Re} = 4250$
, shown in the top rows of panels (a) and (b), up to eight clear vortex structures are present with symmetry about the vertical and horizontal planes. By
${Re} = 4500$
the symmetry of the basic flow is altered as the vortices appear to distort around the circumference of the pipe, detailed in the bottom rows of panels (a) and (b). The temporally stable base flow undergoes a SSBB at Reynolds number
${Re}_{g_{1}} \approx 4375$
, where the flow structure is altered.
Figure 7. Double-bifurcation outflow case,
$124out$
, primary pipe base flow topology, detailing the first SSBB where the geometric symmetry of the Dean vortices is lost. The pipe section views display the spatially evolving, temporally stable, base flows. The top row of panel (a) details the pipe sections at locations
$0D$
to
$-3D$
, along the
$z$
-axis, of the primary pipe at
${Re} = 4250$
where the symmetry is maintained. The bottom row of panel (a) details the pipe sections at locations
$0D$
to
$-3D$
, along the
$z$
-axis, of the primary pipe at
${Re} = 4500$
where the symmetry is lost. Panel (b) details the same information for pipe sections
$-4D$
to
$-7D$
. Both panels are displayed as the local axial velocity
$(U_{ax})$
overlaid by the transverse velocity
$(U_{tr})$
streamlines. The direction of flow is out of the page. The curved arrows indicate the direction of the vortices. The location of each pipe section, in the primary pipe, is shown in the image on the right-hand side of panel (a) in figure 6.
Figure 8. Evolution of the temporally stable base flow for the double-bifurcation outflow case,
$124out$
. The local axial velocity
$(U_{ax})$
is overlaid by the transverse velocity
$(U_{tr})$
streamlines. Each pipe section view is shown at increments of
${Re} = 250$
starting at
${Re} = 4250$
to
${Re} = 6000$
. All section views are from the plane located
$-3 D$
downstream of the origin in the primary pipe shown in the image on the right-hand side of panel (a) in figure 6. The direction of flow is out of the page. The lines used to illustrate the orientation of the vortex centres are denoted as
$S_1$
,
$S_2$
,
$S_3$
and
$S_4$
. The four quadrants are denoted as
$Q_1$
,
$Q_2$
,
$Q_3$
and
$Q_4$
.
The evolution of the temporally stable flow, on a representative selected plane
$(-3D)$
, in the primary pipe is shown in figure 8. Starting at
${Re} = 4250$
, vortex symmetry about the horizontal and vertical planes is well defined. Lines,
$S_1$
and
$S_2$
, used to illustrate the orientation of the vortex centres, show the four main vortices are both orientated at an angle of
$39^{\circ }$
to the horizontal plane. By
${Re} = 4500$
, the symmetry about the horizontal and vertical planes is lost, as the vortices appear to stretch in a clockwise direction around the perimeter of the pipe section. Line
$S_1$
, running through the centres of the vortices lying in the second and fourth quadrants, now makes an angle of
$47^{\circ }$
with the horizontal plane, while line
$S_2$
, running through the centre of the vortices lying in the first and third quadrants, is reduced to an angle of
$26^{\circ }$
with the horizontal plane. The position (the angle between
$S_1$
and the horizontal plane) of the vortex centres lying in quadrants two and four remain fixed from
${Re} = 4500$
to
$5750$
despite the changes in their topology. The two main vortices, lying in quadrants one and three, gradually move closer to the horizontal plane as the Reynolds number is increased throughout the same range (
${Re} = 4500$
to
$5750$
). At
${Re} = 5000$
, the vortices in quadrants two and four continue to stretch around the perimeter of the pipe before two smaller vortices eventually detach by
${Re} = 5500$
. At
${Re} = 5750$
, the two smaller detached vortices have drifted to quadrants one and three where they begin to merge with the two main vortices occupying the same quadrants.
Throughout the Reynolds-number range spanning
${Re} = 4250$
to
$5750$
, where we observe a loss of vortex symmetry about the horizontal and vertical planes, we also observe a
$180^{\circ }$
vortex reflection about both lines,
$S_1$
and
$S_2$
. However, the base flow undergoes a second SSBB, at
${Re}_{g_{2}} \approx 5875$
, where the point symmetry through the origin is lost but a symmetry about the horizontal plane is recovered, shown at
${Re} = 6000$
. We now observe two significant vortices in quadrants two and three. Lines
$S_1$
and
$S_2$
both make an angle of
$28^{\circ }$
with the horizontal plane noting that the lines
$(S_1$
and
$S_2)$
pierce the vortex centres in the adjacent quadrants. That is,
$S_1$
pierces the centre of the vortices lying in
$Q_1$
and
$Q_3$
and
$S_2$
pierces the centre of the vortices lying in
$Q_2$
and
$Q_4$
. The remaining two vortices in quadrants two and three make an angle of
$22^{\circ }$
with the vertical plane, shown as lines
$S_3$
and
$S_4$
, in red. However, they do not pass through a vortex centre in the adjacent quadrants.
The angles of the lines
$S_1$
,
$S_2$
,
$S_3$
and
$S_4$
, with the horizontal plane, illustrated in figure 8, are shown plotted as a function of the Reynolds number in figure 9. The two SSBBs are summarised in table 4.
Another way to summarise the SSBBs is to consider the symmetry group of the resulting flows. Prior to SSBB 1, there are two reflection symmetries, one about the horizontal axis, and one about the vertical axis. After SSBB 1, these reflection symmetries are lost; the remaining symmetry is the combination of a reflection about the horizontal, then a reflection about the vertical, axis. Finally, after SSBB 2, the only remaining symmetry is the reflection about the horizontal axis. These changes in symmetry group are outlined in table 5.
The basic flow condition shown at
${Re} = 6000$
(
${Re} \gt {Re}_{g_{2}}$
) in figure 8 is close to the (Hopf) bifurcation to a time-dependent flow. We note there exists a true basic flow which is mirror symmetric about the
$x$
-axis
$(y = 0)$
and the
$y$
-axis (
$x = 0$
) which undergoes a stationary symmetry-breaking bifurcation leading to the condition shown in figure 8. Therefore, the flows shown in figure 8 could also occur mirrored about either axis, with the choice of direction set by the initial conditions. These changes in the base flow structure affects the mechanism by which the base flow loses stability, in comparison with the first three cases, and is discussed in § 3.3.3.
Table 4. Summary of SSBBs detected from the temporally stable base flows.
Table 5. The symmetry groups of the flow with progression of
${Re}$
for the
$124out$
case. Reflection about the horizontal and vertical axes are designated
$R_H$
and
$R_V$
, respectively.
Figure 9. Evolution of the location of the Dean vortices in the temporally stable base flow for the double-bifurcation outflow case,
$124out$
. The markers show the angular locations of the main vortex centres with respect to the horizontal plane shown in figure 8. On the
$y$
-axis is the absolute value of the angle (in degrees) between the lines
$S_1$
,
$S_2$
,
$S_3$
and
$S_4$
and the horizontal plane as a function of
${Re}$
(
$x$
-axis). The black filled markers show the angle between
$S_1$
and the horizontal plane. The black unfilled markers show the angle between
$S_2$
and the horizontal plane. The red filled marker shows the angle between the lines
$S_3$
and
$S_4$
, and the horizontal plane. The values of
${Re}$
at angles
$39^{\circ }$
and
$28^{\circ }$
(
$4250$
and
$6000$
, respectively) between lines
$S_1$
and
$S_2$
and the horizontal plane are the same. They are shown offset by
${Re} = 20$
for visibility.
3.3. Linear stability analyses
3.3.1. General
In this section we report the results of the linear stability analyses of the steady base flows reported in § 3.2, which are either naturally stable or obtained by stabilising the flow using SFD. The primary result reported here is that all four cases transition to a time-dependent flow via a linear instability mechanism when the Reynolds number exceeds a critical threshold. The leading modes in all cases exhibit complex conjugate eigenvalue pairs, indicating a Hopf bifurcation.
We perform the stability analysis on each base flow, calculating the leading modes, where ‘leading’ refers to the fastest growing, or slowest decaying modes, i.e. those with the largest magnitude multiplier,
$|\mu |$
. We first do this at one value of
${Re}$
, repeating the process at the next value of
${Re}$
, continuing over values of
${Re}$
for which we have calculated base flows. We identify the linear modes belonging to the same solution branch, at each value of
${Re}$
, via their values of
$|\mu |$
, frequency and spatial structure. The transition values of
${Re}$
(those at marginal stability where
$|\mu | = 1$
) for each mode are found by interpolating between the values of
$|\mu |$
. This process allows us to determine the most unstable mode from each case. The most unstable mode transitions at the lowest, or critical, Reynolds number.
3.3.2. Inflow
For inflow case
$12in$
, we predict the most unstable mode to transition at a critical Reynolds number,
${Re}_{12in_{1}} \approx 7369$
, growing with a frequency of
$f \approx 0.52$
. Specifically, to predict the critical Reynolds number (i.e. when
$|\mu | = 1$
), we established a power-law relationship of the multiplier (
$|\mu |$
) as a function of the Reynolds number from the data in table 6. The most unstable mode is referred to as mode 1 (table 6). We detect two more leading modes throughout the range of Reynolds numbers that become unstable very close to
${Re}_{12in_{1}}$
, growing with higher frequencies of
$f \approx 1.46$
and
$f \approx 1.93$
. These two modes become unstable at Reynolds numbers
${Re}_{12in_{2}}\approx 7480$
and
${Re}_{12in_{3}}\approx 7580$
, and are referred to as modes
$2$
and
$3$
, respectively (table 6).
Table 6. Eigenvalues for the leading modes in the single-bifurcation inflow case,
$12in$
, as a function of
${Re}$
. An asterisk indicates the interpolated value of
${Re}$
at which the mode becomes unstable using a power-law fit. Mode 1 becomes unstable at the critical
${Re}$
.
The second inflow case,
$124in$
, becomes unstable at a lower critical Reynolds number than case
$12in$
, while growing with a higher frequency of
$f \approx 0.89$
. The most unstable mode transitions at critical Reynolds number
${Re}_{124in_{1}} \approx 4844$
, and is referred to as mode 1 (table 7). At a slightly higher Reynolds number a second mode becomes unstable, growing with a higher frequency of
$f \approx 1.18$
. The second leading mode becomes unstable at Reynolds number
${Re}_{124in_{2}} \approx 4850$
, and is referred to as mode
$2$
(table 7). The leading modes for both inflow cases, from the data in tables 6 and 7, are plotted in terms of the multiplier and frequency as a function of the Reynolds number, displayed in figure 10.
Table 7. Eigenvalues for the leading modes in the double-bifurcation inflow case,
$124in$
, as a function of
${Re}$
. An asterisk indicates the interpolated value of
${Re}$
at which the mode becomes unstable using a power-law fit. Mode 1 becomes unstable at the critical
${Re}$
.
Figure 10. Multipliers,
$|\mu |$
, and frequencies,
$f$
, of the leading modes from inflow studies
$12in$
and
$124in$
plotted as a function of
${Re}$
. Panel (a) shows the data from case,
$12in$
, from table 6. Panel (b) shows the data from case,
$124in$
, from table 7. The top row in each panel shows the multiplier,
$|\mu |$
, of each mode as a function of
${Re}$
. The dashed line indicates where the multiplier is equal to unity. The bottom row in each panel shows the frequency,
$f$
, of each mode as a function of
${Re}$
. Filled circles indicate the mode is stable. Unfilled circles indicate the mode is unstable.
The structures of the leading modes from both inflow cases in the secondary pipes are qualitatively similar, as shown in figures 11 and 12. In both geometries, when visualised as the axial component of the velocity
$(U_{ax})$
, the modes grow on the inside leg of both the secondary pipes (the outer region of the pipe with respect to the curvature), downstream of the first bifurcation, appearing as narrow crescent-like profiles truncated by small lobes near the centre of the pipe. The crescent-like profile in the double bifurcation,
$124in$
, almost forms a half-circle in comparison with a quarter-circle in the single bifurcation,
$12in$
. It is not until the curvature of the first bifurcation is cleared that the structure of the modes begin to appear. The location where growth occurs in both cases is
${\approx}\, 5D_1$
downstream of the origin in case,
$12in$
(mode 1), and
${\approx}\, 4.1D_1$
downstream of the origin in case,
$124in$
(mode 1). Due to the altered conditions in the double bifurcation, the lobe-like structure of the modes near the centre of the pipe is affected as it approaches the next generation. In case
$12in$
, the growth of each mode is preserved for
${\approx}\, 4{-}5D_1$
, where a clear wavelength is sustained, shown in the inset detail views of figure 11. In case
$124in$
, growth is observable in the outermost tertiary pipes and both secondary pipes. For the most unstable mode (mode 1), growth is most evident in the outermost tertiary pipes with a faint trace in the secondary pipes, displayed in the top rows of panels (a) and (b) of figure 12. The second leading mode, in case
$124in$
, exhibits a slightly more defined structure in the secondary pipes as well as clear growth in the outermost tertiary pipes. Both modes have clear streamwise wavelengths in both generations that survive for
${\approx}\, 2{-}2.4D_1$
in the secondary, and
${\approx}\, 7.4D_1$
in the tertiary, pipes. There is no observable structure in the inner tertiary pipes, for either mode, of case
$124in$
.
Despite the subtle differences observed in the mode structures of the inflow cases in the secondary pipes, and the difference in the critical Reynolds number, the topology of the base flow that triggers the linear instability appears almost identical in both cases. The similarity in the secondary pipe base flows is shown on the left-hand side pipe section views in figures 11 and 12. The leading modes subsequently arise in similar locations, growing in the region of high strain rates on the inside of the secondary pipes where the two base flow vortices start to diverge before wrapping around the perimeter of the pipe. In both inflow cases, the spatial position in the geometry where the growth is observed coincides with this typical base flow topology when the Reynolds number exceeds the critical threshold. Figure 13 contrasts the base flow topology for cases
$124in$
(top half) and
$12in$
(bottom half) illustrating the similarities in the vortex size and orientation.
Figure 11. Leading unstable modes for the single-bifurcation inflow case,
$12in$
. The image on the left shows the region of growth of all unstable modes in the secondary pipe. The inset detail views show the wavelength, from top to bottom, of modes
$1$
,
$2$
and
$3$
, respectively, as the axial component of the velocity,
$U_{ax}$
. The pipe section views show the typical base flow condition as the streamlines of the transverse component of the velocity,
$U_{tr}$
(left), the mode structure as the axial component of the velocity,
$U_{ax}$
(centre), and the mode structure as the horizontal transverse component of the velocity,
$U_{tr}$
(right). The velocity perturbations are normalised by the maximum magnitude. The left-hand side of each pipe section view (side A) pertains to the inside of the secondary pipe (side A) of the far left image. The right-hand side of each pipe section view (side B) pertains to the outside of the secondary pipe (side B) of the far left image. The direction of flow through the pipe sections is into the page. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the image on the left-hand side.
Figure 12. Leading unstable modes in the double-bifurcation inflow case,
$124in$
. The top row of images are from the secondary pipe. The bottom row of images are from the tertiary pipe. The image on the left of each row shows the region of growth of all unstable modes in the secondary and tertiary pipes. The inset detail views show the wavelength, from top to bottom, of modes
$1$
and
$2$
, respectively, as the axial component of the velocity,
$U_{ax}$
. The pipe section views show the typical base flow condition as the streamlines of the transverse component of the velocity,
$U_{tr}$
(left), the mode structure as the axial component of the velocity,
$U_{ax}$
(centre), and the mode structure as the horizontal transverse component of the velocity,
$U_{tr}$
(right). The velocity perturbations are normalised by the maximum magnitude. The left-hand side of each pipe section view (side A) pertains to the inside of the secondary or tertiary pipe (side A) of the far left image. The right-hand side of each pipe section view (side B) pertains to the outside of the secondary or tertiary pipe (side B) of the far left image. The direction of flow through the pipe sections is into the page. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the images on the left-hand side.
3.3.3. Outflow
For the outflow case,
$12out$
, we detect a distinct leading mode that becomes unstable at critical Reynolds number
${Re}_{12out_{1}} \approx 4734$
. The mode grows with a frequency of
$f \approx 3.21$
, and is referred to as mode
$1$
(table 8). No other modes approaching transition are observed. The flow in the second outflow case,
$124out$
, becomes unstable at a higher critical Reynolds number than the single-bifurcation case,
$12out$
, while growing with a lower frequency of
$f \approx 0.23$
. The most unstable mode transitions at critical Reynolds number,
${Re}_{124out_{1}} \approx 6000$
, and is referred to as mode 1 (table 9). At a higher Reynolds number a second mode becomes unstable, growing with a higher frequency of
$f \approx 2.11$
. The second mode becomes unstable at Reynolds number,
${Re}_{124out_{2}} \approx 6312$
, and is referred to as mode
$2$
(table 9). The leading modes for both outflow cases, from the data in tables 8 and 9, are plotted in terms of the multiplier and frequency as a function of the Reynolds number, displayed in figure 14.
The structure of the leading unstable mode for case
$12out$
is shown in figure 15. The mode grows on a fully symmetric four-vortex base flow in the high strain regions where the two counter-rotating vortex pairs diverge. The mode appears at
${\approx}\, 0.25D_1$
downstream of the origin and is preserved for
${\approx}\, 2.95D_1$
with a clear wavelength. In case
$124out$
, we observe a completely different mode structure due to the changes in the base flow, shown in figure 16. The most unstable mode (mode 1) grows in the primary and secondary pipes. In the primary pipe, growth is apparent across the entire pipe section, concentrated in regions around the cores of the main vortices. It appears at
${\approx}\, 1.0D_1$
and is preserved for
${\approx}\, 5.3D_1$
. In the secondary pipe, a wavelength is evident, however, a relationship between the regions of growth, shown in the pipe section, and the vortices is not apparent. The second leading mode is concentrated in the primary pipe, appearing at
${\approx}\, 1.0D_1$
and preserved for
${\approx}\, 6.2D_1$
. There is no observable structure in the secondary pipe.
Similar to both inflow cases, the distinct unstable mode observed in outflow case
$12out$
grows in regions of high strain in the centre of the primary pipe. In contrast, in the outflow case
$124out$
, the modes no longer display clear wavelengths and the growth does not emanate from the high strain regions associated with the diverging base flow vortices. In the primary pipe, the growth now appears to concentrate in the regions around the main vortex cores, occupying large areas of the pipe section. Further, we note that, in the
$124out$
case, the unstable mode bifurcates from a base flow that is not reflection symmetric about the
$y$
-axis, and therefore the unstable mode can grow from the base flow shown in figure 8, or its mirror image. In the secondary pipe, in the case of mode 1, there is no clear relationship between the location of the growth and the base flow vortices.
Table 8. Eigenvalues for the leading mode in the single-bifurcation outflow case,
$12out$
, as a function of
${Re}$
. An asterisk indicates the interpolated value of
${Re}$
at which the mode becomes unstable using a power-law fit. Mode 1 becomes unstable at the critical
${Re}$
.
Figure 13. Comparison of the typical base flow topology present in the single- and double-bifurcation inflow cases,
$12in$
and
$124in$
. The image contrasts the size and orientation of the vortices in the single-bifurcation (bottom half) and double-bifurcation (top half) studies which trigger the linear instability after
${Re}$
has exceeded the critical threshold. The red crosses indicate the location of the vortex centres. The curved arrows indicate the direction of the vortices.
Table 9. Eigenvalues for the leading modes in the double-bifurcation outflow case,
$124out$
, as a function of
${Re}$
. An asterisk indicates the interpolated value of
${Re}$
at which the mode becomes unstable using a power-law fit. Mode 1 becomes unstable at the critical
${Re}$
.
Figure 14. Multipliers,
$|\mu |$
, and frequencies,
$f$
, of the leading modes from outflow studies
$12out$
and
$124out$
plotted as a function of
${Re}$
. Panel (a) shows the data from case,
$12out$
, from table 8. Panel (b) shows the data from case,
$124out$
, from table 9. The top row in each panel shows the multiplier,
$|\mu |$
, of each mode as a function of
${Re}$
. The dashed line indicates where the multiplier is equal to unity. The bottom row in each panel shows the frequency,
$f$
, of each mode as a function of
${Re}$
. Filled circles indicate the mode is stable. Unfilled circles indicate the mode is unstable.
Figure 15. Leading unstable mode in the single-bifurcation outflow case,
$12out$
. The image on the left shows the region of growth of the unstable mode in the primary pipe. The inset detail view shows the wavelength as the axial component of the velocity,
$U_{ax}$
. The pipe section views show the typical base flow condition as the streamlines of the transverse component of the velocity,
$U_{tr}$
(left), the mode structure as the axial component of the velocity,
$U_{ax}$
(centre), and the mode structure as the horizontal transverse component of the velocity,
$U_{tr}$
(right). The velocity perturbations are normalised by the maximum magnitude. The left-hand side of each pipe section view (side A) pertains to the left-hand side of the primary pipe (side A) of the far left image. The right-hand side of each pipe section view (side B) pertains to the right-hand side of the primary pipe (side B) of the far left image. The direction of flow through the pipe sections is out of the page. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the image on the left-hand side.
Figure 16. Leading unstable modes in the double-bifurcation outflow case,
$124out$
. The top row of images are from the primary pipe. The bottom row of images are from the secondary pipe. The image on the left of each row shows the region of growth of all unstable modes in the primary and secondary pipes. The inset detail views show the wave structure, from top to bottom, of modes
$1$
and
$2$
, respectively, as the axial component of the velocity,
$U_{ax}$
. The pipe section views show the typical base flow condition as the streamlines of the transverse component of the velocity,
$U_{tr}$
(left), the mode structure as the axial component of the velocity,
$U_{ax}$
(centre), and the mode structure as the horizontal transverse component of the velocity,
$U_{tr}$
(right). The velocity perturbations are normalised by the maximum magnitude. The left-hand side of each pipe section view (side A) pertains to the left-hand side of the primary pipe and the inside of the secondary pipe (side A) of the far left image. The right-hand side of each pipe section view (side B) pertains to the right-hand side of the primary pipe and the outside of the secondary pipe (side B) of the far left image. The direction of flow through the pipe sections is out of the page. The curved arrows indicate the direction of the vortices. The location of each pipe section is with respect to the position along the
$z$
-axis, from the origin. The orientation of the coordinate system relates to the images on the left-hand side.
4. Comparison with direct numerical simulation
We have carried out a series of DNS studies to compare the results from the linear stability analyses with the nonlinear solutions. We compare the critical Reynolds numbers and the frequencies, calculated from the time histories obtained from DNS, with the linear stability results. We also compare the mode structures obtained from the linear stability results with the mode structures from the DNS solutions as well as compare the dominant wavenumbers of the linear and nonlinear cases obtained from spatial analyses of the streamwise flows.
For the frequency comparisons, the DNS solutions were restarted from the same stable base flows about which the nonlinear Navier–Stokes equations were linearised (detailed in § 2.5), with the SFD switched off. Each case was run until a statistical steady state was reached, enabling inspection of the frequencies from both the early and saturated time histories. The frequencies were obtained by fast Fourier transform within a sliding Hanning window. In both inflow cases, we observe long transient periods in the early time history that differ to what is observed in the saturated DNS solution. The linear stability analyses, from the inflow studies, detect the early transient frequencies and match closely with the saturated nonlinear frequencies. In comparison with the inflow studies, the outflow cases reach a saturated solution almost instantaneously in the DNS time history, which matches well with the linear stability analyses from the outflow studies.
For the mode structure comparisons, we visualise the DNS flow topology at similar pipe section locations to where the linear modes were visualised in § 3.3. For the cases in the single-bifurcation geometry (
$12in$
and
$12out$
), we observe an almost identical match between the mode structures. For the cases in the double-bifurcation geometry (
$124in$
and
$124out$
), we observe very similar structures in the pipes where growth is dominant (the tertiary pipe in case
$124in$
and the primary pipe in case
$124out$
). In the pipes where growth is faint, the similarities are fewer. Similarly, for the comparisons of the spatial analyses, we illustrate the results from identical streamwise locations to where the corresponding linear modes were visualised in § 3.3. The velocity data for the spatial analyses were obtained from a line orientated through regions of prominent mode growth in the streamwise direction, parallel to the local pipe centreline. The instantaneous wavenumbers were obtained by fast Fourier transform within a sliding Hanning window. In all cases, we observe similar dominant wavenumbers from the linear and nonlinear flows over similar spatial locations.
In §§ 4.1 and 4.2 we compare the nonlinear characteristics with the stability results for the inflow and outflow studies, respectively.
4.1. Inflow
Figures 17–20 show samples of the time history and frequency spectrum from cases,
$12in$
and
$124in$
, respectively. In both cases, the time history shown in figures 17 and 19 have been measured at a Reynolds number slightly higher than when the leading modes are all unstable. In both figures, the DNS has been restarted from the stabilised base flow, with the SFD switched off, and run until the solution saturates to a statistical steady state.
For the first DNS inflow case,
$12in$
, we detect a dominant low-frequency signal in the spectrum of
$f \approx 0.70$
, at time
$\tau \lt 1308$
, shown in panel (a) of figure 18. At time
$1309 \lt \tau \lt 1315$
, in the same figure, we observe two higher-frequency components of
$f \approx 1.52$
and
$f \approx 1.99$
. Panel (b) in figure 18 shows the changes to the spectrum over the time interval
$1315 \lt \tau \lt 1450$
. The highest-frequency mode (
$f \approx 1.99$
) becomes less dominant with only a light trace observable at time
$\tau \gt 1380$
. The lower-frequency mode (
$f = 1.52$
) is observed briefly at time
$\tau \approx 1340$
, before it reduces and the spectrum starts to break up (
$\tau \gt 1400$
). At time
$1450 \lt \tau \lt 1600$
, in panel (c) there is still a light trace of the high-frequency signal (
$f \approx 1.99$
) but the dominant frequency is much lower (
$f \approx 0.55$
). At times
$\tau \gt 1850$
, in panel (d), there is a single dominant frequency in the spectrum. The most unstable mode, detected from the linear stability, oscillates with a frequency of
$f \approx 0.52$
, matching closely with the initial frequency from the early time history in panel (a). Further, it is almost identical with the saturated frequency from panel (d) of figure 18. The two higher transient frequencies of
$f \approx 1.52$
and
$f \approx 1.99$
, that appear most dominant at time
$1309 \lt \tau \lt 1317$
, in panels (a) and (b) are close to the second and third unstable modes detected from the linear stability analysis with frequencies of
$f \approx 1.46$
and
$f \approx 1.93$
, respectively.
Figure 17. Velocity time series from the single-bifurcation inflow study,
$12in$
, at
${Re} = 7750$
. The sample has been restarted from the stabilised base flow with the SFD switched off at time
$\tau = 1295$
and run until a statistical steady state was reached. Panel (a) shows the full time series,
$1295 \lt \tau \lt 2220$
. Panel (b) is at time
$1295 \lt \tau \lt 1315$
, from the sample between the first two vertical dashed lines in panel (a). Panel (c) is at time
$1315 \lt \tau \lt 1450$
, from the sample between the second and third vertical dashed lines in panel (a). Panel (d) is at time
$2170 \lt \tau \lt 2220$
, from the sample between the fourth and fifth vertical dashed lines in panel (a). The time series is the
$x$
-component of velocity from a probe point located
${\approx}\, 5.0D_1$
along the
$z$
-axis, from the origin, in the secondary pipe.
Figure 19 shows the time history samples from case
$124in$
. In panel (a) of figure 20 at times
$2435 \lt \tau \lt 2447$
, we detect a strong frequency signal of
$f \approx 0.88$
. This signal starts to break up over the time interval
$2448 \lt \tau \lt 2600$
, shown in panels (a) and (b). At time
$2600 \lt \tau \lt 2700$
, the signal appears to saturate, shown in panel (c), which is confirmed in panel (d) over the time interval
$2700 \lt \tau \lt 2800$
, where the signal appears almost identical. The saturated frequency spectrum appears to consists of three components; a high, low and intermediate frequencies of
$f \approx 1.12$
,
$0.88$
and
$1.0$
, respectively. In both panels (c) and (d) of figure 19 a time-dependent switching is apparent where we observe intermittent mode dominance, potentially due to pattern competition (Ciliberto & Gollub Reference Ciliberto and Gollub1984). The almost periodic switching in the nonlinear solution sees the frequency oscillate between the two frequencies similar to those detected from the linear stability analysis. The initial frequency of
$f \approx 0.88$
, in the early time history
$(\tau \lt 2447)$
is almost identical with the frequency of the most unstable mode with a frequency
$f \approx 0.89$
, which also forms the lower bound of the switching in the saturated flow in panels (c) and (d). The second unstable mode detected from the linear stability analysis with a frequency of
$f \approx 1.18$
appears slightly higher than the upper bound of the switching in the saturated flow with a peak frequency of
$f \approx 1.12$
.
Throughout the frequency spectrum of case
$124in$
, we also observe a very low-frequency signal of
$f \approx 0.1$
that is most prominent in panel (a). Interestingly, in the saturated solutions shown in panels (c) and (d), we observe this signal, albeit with a very low intensity, appearing with similar periodicity; essentially at a similar time to when the main signal jumps to the upper bound, i.e. at times
$\tau \approx 2620$
and
$2665$
, in panel (c) and times
$\tau \approx 2715$
and
$2760$
, in panel (d). However, this frequency is not detected in the linear stability analysis.
Figure 18. Frequency spectrum for the single-bifurcation inflow study,
$12in$
, calculated from the time series shown in figure 17. Panel (a) is the inital spectrum from the sample restarted from the stabilised base flow over the time interval
$1295 \lt \tau \lt 1315$
. Panel (b) is the spectrum over the time interval
$1315 \lt \tau \lt 1450$
. Panel (c) is the spectrum from the time interval
$1450 \lt \tau \lt 1850$
. Panel (d) is the spectrum over the time interval
$1850 \lt \tau \lt 2220$
. The leading linear modes in each panel are shown as solid, dashed and dotted lines.
Figure 19. Velocity time series from the double-bifurcation inflow study,
$124in$
, at
${Re} = 5000$
. The sample has been restarted from the stabilised base flow with the SFD switched off at time
$\tau = 2435$
and run until a statistical steady state was reached. Panel (a) shows the full time series
$2435 \lt \tau \lt 2750$
. Panel (b) is at time
$2435 \lt \tau \lt 2455$
, from the sample between the first two vertical dashed lines in panel (a). Panel (c) is at time
$2455 \lt \tau \lt 2700$
, from the sample between the second and third vertical dashed lines in panel (a). Panel (d) is at time
$2700 \lt \tau \lt 2750$
, from the sample between the third and fourth vertical dashed lines in panel (a). The time series is the
$x$
-component of velocity from a probe point located
${\approx}\, 5.0D_1$
along the
$z$
-axis, from the origin, in the secondary pipe.
Figure 20. Frequency spectrum for the double-bifurcation inflow study,
$124in$
, calculated from the time series shown in figure 19. Panel (a) is the initial spectrum from the sample restarted from the stabilised base flow over the time interval
$2435 \lt \tau \lt 2455$
. Panel (b) is the spectrum over the time interval
$2455 \lt \tau \lt 2600$
. Panel (c) is the spectrum from the time interval
$2600 \lt \tau \lt 2700$
. Panel (d) is the spectrum from the time interval
$2700 \lt \tau \lt 2800$
. The leading linear modes in each panel are shown as solid and dashed lines.
Shown in figures 21 and 22 is a visual comparison of the mode structures from the linear stability analyses and the DNS solutions for the inflow cases
$12in$
and
$124in$
, respectively. The structures are compared in the pipe section views at similar locations. For the nonlinear DNS flows we illustrate the structure at two values of the Reynolds number to see if the mechanisms observed from the linear stability analyses are likely to be preserved at values above the critical threshold. To obtain the nonlinear mode structures the DNS solution was decomposed into a mean velocity component and a temporally fluctuating velocity component,
$\boldsymbol {u} = \overline {\boldsymbol {u}} + \hat {\boldsymbol {u}}$
. The displayed fluctuating velocity component,
$\hat {\boldsymbol {u}}$
, was calculated by subtracting the mean velocity field,
$\overline {\boldsymbol {u}}$
, from the full velocity field,
$\boldsymbol {u}$
.
For the single-bifurcation inflow case,
$12in$
in figure 21, the nonlinear modes from the DNS are displayed at
${Re} = 7750$
and
$8250$
. At
${Re} = 7750$
the structure is almost identical to what is observed from the linear stability. At
${Re} = 8250$
the profile of the crescent is slightly reduced but the basic structure is preserved.
For the double-bifurcation inflow case,
$124in$
in figure 22, the nonlinear modes from the DNS are displayed at
${Re} = 5250$
and
$5500$
. At
${Re} = 5250$
in the secondary pipe, the cresent on the left-hand side of the pipe section in the nonlinear mode is smaller than what is observed in the linear mode. However, at
${Re} = 5500$
these profiles are much closer in appearance. In the top row of figure 22 the structure from the secondary pipe is barely visible at
${Re} = 5250$
, but is slightly more defined at
${Re} = 5500$
. At both values of Reynolds number the structures are less prevalent in the secondary pipes (compared with the tertiary pipe), similar to what is observed in the linear modes. In the tertiary pipe, the nonlinear mode is clearly visible at both values of Reynolds number. At
${Re} = 5250$
the structure in the tertiary pipe is visible in the outer pipe only, the same as what is observed from the linear stability analysis.
Figure 21. Comparison of linear (a) and nonlinear (b) mode structures for the single-bifurcation inflow study,
$12in$
. The linear modes
$1$
,
$2$
and
$3$
are identical to the images displayed in figure 11. The nonlinear mode structures from DNS are at
${Re} = 7750$
and
${Re} = 8250$
. All images are displayed as the local axial component of velocity,
$U_{ax}$
. Nonlinear mode at
${Re} = 7750$
colour range,
$U_{ax} \pm 10^{-4}$
. Nonlinear mode at
${Re} = 8250$
colour range,
$U_{ax} \pm 3.0 \times 10^{-4}$
.
Figure 22. Comparison of linear (a) and nonlinear (b) mode structures for the double-bifurcation inflow study,
$124in$
. The linear modes
$1$
and
$2$
are identical to the images displayed in figure 12. The nonlinear mode structures from DNS are at
${Re} = 5250$
and
${Re} = 5500$
. The images in the top row are from the secondary pipe. The images in the bottom row are from the tertiary pipe. All images are displayed as the local axial component of velocity,
$U_{ax}$
. Nonlinear mode at
${Re} = 5250$
colour range,
$U_{ax} \pm 10^{-2}$
. Nonlinear mode at
${Re} = 5500$
colour range,
$U_{ax} \pm 10^{-1 }$
.
Shown in figures 23, 24 and 25 are the spectra from the streamwise spatial analyses for both inflow cases. Figure 23 displays the wavenumbers detected in the secondary pipe from case
$12in$
, and figures 24 and 25 display the wavenumbers detected in the secondary and outer tertiary pipes, respectively, from case
$124in$
. In all plots, the spectra of the linear modes are calculated for data extracted over an identical streamwise location as the corresponding mode displayed in § 3.3. The results from DNS are calculated from the temporally fluctuating component,
$\hat {\boldsymbol {u}}$
, of the decomposed flows shown in figures 21 and 22 and are displayed over the entire range.
Mode
$1$
, from case
$12in$
, displays a dominant wavenumber of
$k \approx 1.0$
that remains largely unchanged, shown in panel (b) of figure 23, before the signal breaks up at location
${\approx}\, 9.0D_1$
. Mode
$2$
in panel (c) and mode
$3$
in panel (d) show clear wavenumbers which gradually start to increase near the end of the domain analysed. Mode
$2$
has a dominant wavenumber of
$k \approx 2.46$
that increases to
$k \approx 2.6$
, by location
${\approx}\, 10.6D_1$
before it eventually breaks up. Mode
$3$
displays a slightly higher dominant wavenumber of
$k \approx 3.61$
. The wavenumber (k
${\approx}\, 3.61$
) gradually starts increasing after location
${\approx}\, 8.5D_1$
to
$k \approx 4.9$
by location
$10.5D_1$
. From the DNS spectrum in panel (a) a dominant wavenumber of
$k \approx 1.0$
is detected that is almost identical to mode
$1$
.
Figure 24 displays the wavenumbers in the secondary pipe for inflow case,
$124in$
. Mode
$1$
in panel (b) displays a dominant wavenumber of
$k \approx 0.7$
that remains unchanged until location
${\approx}\, 5.6D_1$
. The signal then breaks up and a higher wavenumber of up to
$k \approx 8.5$
can be observed. Mode
$2$
, in panel (c) displays a dominant wavenumber of
$k \approx 2.51$
up to location
${\approx}\, 4.9D_1$
before it increases to
$k \approx 3.8$
. The signal is observed to increase again to
$k \approx 7.0$
between locations
${\approx}\, 5.7D_1$
and
${\approx}\, 6.2D_1$
. The DNS spectrum in panel (a) displays a clear wavenumber of
$k \approx 2.0$
between locations
${\approx}\, 4.2D_1$
and
${\approx}\, 5.9D_1$
. From location
$5.9D_1$
to
$6.5D_1$
a lower wavenumber of
$k \approx 0.7$
is detected. In the outer tertiary pipe, in figure 25, mode
$1$
displays a dominant wavenumber of
$k \approx 2.05$
in panel (b) before the signal starts to change after location
${\approx}\, 22.5D_1$
. Mode
$2$
displays a higher dominant wavenumber of
$k \approx 2.97$
up to location
${\approx}\, 22.0D_1$
before the signal starts to increase to
$k \approx 5.0$
. The DNS spectrum in panel (a) displays a dominant wavenumber of
$k \approx 1.92$
up to location
${\approx}\, 22.5D_1$
, slightly less than the wavenumber observed from mode
$1$
(
$k \approx 2.05$
). Thereafter, at least two short-lived signals of
$k \approx 3.8$
and
$1.7$
are observed.
Results from the linear stability analyses and the DNS, for inflow case
$124in$
, show the growth is predominantly in the outer tertiary pipes. Only faint traces are observable in the secondary generation. Subsequently, the streamwise spatial analyses detect no obvious similarities between the linear and nonlinear structures, unlike what is observable in the outer tertiary pipe. The tertiary generation shows a close match between the DNS and linear mode
$1$
where the growth is well established. For the first case,
$12in$
, the streamwise spatial analysis also indicates a close match between linear mode
$1$
and the saturated DNS flow, similar to what is observed from the frequency spectrum in figure 18.
Figure 23. Single-bifurcation inflow case,
$12in$
, streamwise spatial analysis. In all plots the spectrum is obtained from the fast fourier transform (FFT) of the local axial component of the velocity,
$(U_{ax})$
, over a line orientated in the secondary pipe. The line spans the coordinates
$x = 1.35$
,
$y = 0$
,
$z = 4.15$
to
$x = 5.88$
,
$y = 0$
,
$z = 11.4$
, orientating it parallel to the centreline near the inner wall. The
$x$
-axis is the streamwise position along the line in the secondary pipe. The
$y$
-axis is the wavenumber
$(k)$
. Each linear mode is displayed over an identical streamwise position and length as the corresponding mode displayed in the inset detail view in figure 11. The nonlinear plot is displayed over the entire range (
$5.0D$
to
$11.5D$
). Panel (a) is the nonlinear solution from DNS at
${Re} = 7750$
. Panel (b) is linear mode
$1$
. Panel (c) is linear mode
$2$
. Panel (d) is linear mode
$3$
.
Figure 24. Double-bifurcation inflow case,
$124in$
, streamwise spatial analysis. In all plots the spectrum is obtained from the FFT of the local axial component of the velocity,
$(U_{ax})$
, over a line orientated in the secondary pipe. The line spans the coordinates
$x = 0.92$
,
$y = 0$
,
$z = 3.5$
to
$x = 2.795$
,
$y = 0$
,
$z = 6.5$
, orientating it parallel to the centreline near the inner wall. The
$x$
-axis is the streamwise position along the line in the secondary pipe. The
$y$
-axis is the wavenumber
$(k)$
. Each linear mode is displayed over an identical streamwise position and length as the corresponding mode displayed in the inset detail view in figure 12. The nonlinear plot is displayed over the entire range (
$4.1D$
to
$6.5D$
). Panel (a) is the nonlinear solution from DNS at
${Re} = 5500$
. Panel (b) is linear mode
$1$
. Panel (c) is linear mode
$2$
.
The linear stability analysis predicts the inflow cases become unstable at critical Reynolds numbers
${Re} \approx 7369$
and
${Re} \approx 4844$
for cases
$12in$
and
$124in$
, respectively. The DNS results in § 3.1 approximated critical Reynolds numbers of
${Re} \approx 7381$
and
${Re} \approx 5051$
for cases,
$12in$
and
$124in$
, respectively. The inflow comparisons are summarised in the first five rows of table 10.
Table 10. Summary of transition (neutrally stable)
${Re}$
and frequencies detected from linear stability and DNS. An asterisk indicates the critical
${Re}$
for each case.
Figure 25. Double-bifurcation inflow case,
$124in$
, streamwise spatial analysis. In all plots the spectrum is obtained from the FFT of the local axial component of the velocity,
$(U_{ax})$
, over a line orientated in the tertiary pipe. The line spans the coordinates
$x = 5.0$
,
$y = 0$
,
$z = 7.35$
to
$x = 11.05$
,
$y = 0$
,
$z = 10.3$
, orientating it parallel to the centreline approximately through the centre of the pipe. The
$x$
-axis is the streamwise position along the line in the tertiary pipe. The
$y$
-axis is the wavenumber
$(k)$
. Each linear mode is displayed over an identical streamwise position and length as the corresponding mode displayed in the inset detail view in figure 12. The nonlinear plot is displayed over the same range (
$17.6D_1$
to
$25.0D_1$
). Panel (a) is the nonlinear solution from DNS at
${Re} = 5500$
. Panel (b) is linear mode
$1$
. Panel (c) is linear mode
$2$
.
4.2. Outflow
Figures 26–29 show the samples of the time history and frequency spectrum from cases,
$12out$
and
$124out$
, respectively. In both cases the time histories shown in figures 26 and 28 have been measured at a Reynolds number slightly higher than when the leading modes are all unstable. In both figures, the DNS has been restarted from the stabilised base flow, with the SFD switched off, and run until the solution saturated to a statistical steady state.
Figure 26 shows the sample of the time history from the first DNS outflow case,
$12out$
. We detect three short-lived low-frequency signals in the spectrum (
$f \approx 1.2$
,
$f \approx 1.0$
and
$f \approx 0.1$
) at time
$\tau \lt 1815$
, shown in panel (a) of figure 27. At times
$\tau \gt 1815$
, there is a single dominant frequency in the spectrum of
$f \approx 3.21$
. From the linear stability, we detect a single unstable mode with an identical frequency (
$f \approx 3.21$
).
Figure 28 shows the sample of the time history from case,
$124out$
. In panel (a) of figure 29 at time,
$\tau \gt 3935$
, we detect a signal almost immediately of
$f \approx 0.25$
. At time
$3936 \lt \tau \lt 3938$
, we observe a short-lived component in the spectrum with a low intensity peak frequency of
$f \approx 1.95$
. At time
$3940 \lt \tau \lt 3965$
, the spectrum begins to saturate with a dominant frequency plus low intensity peaks of up to
$f \approx 2.1$
, shown in panels (a) and (b). At time
$\tau \gt 4025$
, in panels (c) and (d), there is a clear low-frequency signal. The first frequency of
$f \approx 0.25$
closely matches with the most unstable mode detected from the linear stability analysis of
$f \approx 0.23$
. The short-lived component in the spectrum at time
$3936 \lt \tau \lt 3938$
, with the low intensity peak frequency of
$f \approx 1.95$
, is close to our second leading mode with a frequency of
$f \approx 2.11$
.
Figure 26. Velocity time series from the single-bifurcation outflow study,
$12out$
, at
${Re} = 4750$
. The sample has been restarted from the stabilised base flow with the SFD switched off at time
$\tau = 1785$
and run until a statistical steady state was reached. Panel (a) shows the full time series,
$1785 \lt \tau \lt 2360$
. Panel (b) is at time
$1785 \lt \tau \lt 1800$
, from the sample between the first two vertical dashed lines in panel (a). Panel (c) is at time
$1800 \lt \tau \lt 1825$
, from the sample between the second and third vertical dashed lines in panel (a). Panel (d) is at time
$2340 \lt \tau \lt 2350$
, from the sample between the fourth and fifth vertical dashed lines in panel (a). The time series is the
$x$
-component of velocity from a probe point located
${\approx}\, -4.0D_1$
along the
$z$
-axis, from the origin, in the primary pipe.
Figure 27. Frequency spectrum for the single-bifurcation outflow study,
$12out$
, calculated from the time series shown in figure 26. Panel (a) is the initial spectrum from the sample restarted from the stabilised base flow over the time interval
$1785 \lt \tau \lt 2015$
. Panel (b) is the spectrum over the time interval
$2015 \lt \tau \lt 2350$
. The linear mode in each panel is shown as a solid line.
Figure 28. Velocity time series from the double-bifurcation outflow study,
$124out$
, at
${Re} = 6500$
. The sample has been restarted from the stabilised base flow with the SFD switched off at time
$\tau = 3925$
and run until a statistical steady state was reached. Panel (a) shows the full time series,
$3925 \lt \tau \lt 4660$
. Panel (b) is at time
$3925 \lt \tau \lt 3945$
, from the sample between the first two vertical dashed lines in panel (a). Panel (c) is at time
$3945 \lt \tau \lt 4025$
, from the sample between the second and third vertical dashed lines in panel (a). Panel (d) is at time
$4400 \lt \tau \lt 4650$
, from the sample between the fourth and fifth vertical dashed lines in panel (a). The time series is the
$x$
-component of velocity from a probe point located
${\approx}\, -6.0D_1$
along the
$z$
-axis, from the origin, in the primary pipe.
Figure 29. Frequency spectrum for the double-bifurcation outflow study,
$124out$
, calculated from the time series shown in figure 28. Panel (a) is the initial spectrum from the sample restarted from the stabilised base flow over the time interval
$3935 \lt \tau \lt 3950$
. Panel (b) is the spectrum over the time interval
$3950 \lt \tau \lt 4025$
. Panel (c) is the spectrum from the time interval
$4025 \lt \tau \lt 4400$
. Panel (d) is the spectrum over the time interval
$4400 \lt \tau \lt 4650$
. The leading linear modes in each panel are shown as solid and dashed lines.
Shown in figures 30 and 31 is a visual comparison of the mode structures from the linear stability analyses and the DNS solutions for outflow cases
$12out$
and
$124out$
, respectively, obtained in the same way as the inflow comparisons. For the single-bifurcation outflow case,
$12out$
in figure 30, the nonlinear modes from the DNS are displayed at
${Re} = 4750$
and
$5250$
. At
${Re} = 4750$
the structure is almost identical to what is observed from the linear stability. At
${Re} = 5250$
the structure is slightly modified but the basic structure is preserved, similar to the single-bifurcation inflow study when the Reynolds number is increased.
For the double-bifurcation outflow case,
$124out$
in figure 31, the nonlinear modes from the DNS are displayed at
${Re} = 6250$
and
$6500$
. In the top row is the mode structure from the primary pipe and in the bottom row is the mode structure from the secondary pipe. At
${Re} = 6250$
the growth of the nonlinear mode is concentrated in the primary pipe and almost undetectable in the secondary pipe. The structure in the primary pipe is very similar, but not identical, to mode 1 detected from the linear stability. The nonlinear mode at
${Re} = 6250$
appears predominantly in the top and bottom regions of the pipe (with respect to its orientation on the page) around the periphery of where the vortex centres appear, similar to both the linear mode structures. At
${Re} = 6500$
the nonlinear mode structure in the primary pipe has altered significantly with only small remnants of the structure that are observable at
${Re} = 6250$
, like the small lobe-like profiles around the periphery of the pipe section, being preserved when the Reynolds number is increased above the critical threshold. At
${Re} = 6500$
, structures in the secondary pipe are observable, and similar to the linear modes although dissimilar in topology, appear relatively independent of the base flow vortex structures. Unlike the previous three cases which appear to maintain the basic mode structure detected from the linear stability analysis when the Reynolds number is increased (over the values tested), the comparison of the linear and nonlinear mode structures for the double-bifurcation outflow case,
$124out$
, suggests that this mechanism may only be present near the critical threshold.
Figure 30. Comparison of linear (a) and nonlinear (b) mode structures for the single-bifurcation outflow study,
$12out$
. The linear mode is identical to the image displayed in figure 15. The nonlinear mode structures from DNS are at
${Re} = 4750$
and
${Re} = 5250$
. All images are displayed as the local axial component of velocity,
$U_{ax}$
. Nonlinear mode at
${Re} = 4750$
colour range,
$U_{ax} \pm 5.0 \times 10^{-2}$
. Nonlinear mode at
${Re} = 5250$
colour range,
$U_{ax} \pm 7.0 \times 10^{-2}$
.
Figure 31. Comparison of linear (a) and nonlinear (b) mode structures for the double-bifurcation outflow study,
$124out$
. The linear modes
$1$
and
$2$
are identical to the images displayed in figure 16. The nonlinear mode structure from DNS are at
${Re} = 6250$
and
${Re} = 6500$
. The images in the top row are from the primary pipe. The images in the bottom row are from the secondary pipe. All images are displayed as the local axial component of velocity,
$U_{ax}$
. Nonlinear mode at
${Re} = 6250$
colour range,
$U_{ax} \pm 2.0 \times 10^{-2}$
. Nonlinear mode at
${Re} = 6500$
colour range,
$U_{ax} \pm 1.5 \times 10^{-1}$
.
Shown in figures 32, 33 and 34 are the streamwise spatial analyses for both outflow cases, obtained in the same way as the inflow studies. Figure 32 displays the wavenumbers in the primary pipe from case
$12out$
. Mode
$1$
displays a dominant wavenumber of
$k \approx 6.1$
, that is unchanged throughout the spectrum, shown in panel (b). From the DNS spectrum, in panel (a), a dominant wavenumber of
$k \approx 6.15$
is detected that is almost identical to the linear mode.
Figures 33 and 34 display the wavenumbers in the primary and secondary pipes respectively, from case
$124out$
. In the primary pipe, mode
$1$
in panel (b), has a dominant wavenumber of
$k \approx 0.32$
from location
$1.0D_1$
to
$-1.5D_1$
, before briefly increasing, followed immediately by a decrease to
$k \approx 0.25$
at location
$-4.3D_1$
. Mode
$2$
, in panel (c) displays a clear wavenumber of
$k \approx 1.2$
from spatial position
$1.0D_1$
to
$-2.2D_1$
before the signal breaks up and at least two (higher and lower) wavenumbers are detected. From the DNS spectrum in panel (a), a dominant wavenumber of
$k \approx 0.29$
is observed between locations
$0.5D_1$
to
$-5.2D_1$
. Slightly higher short-lived signals are also detected throughout the spectrum of up to
$k \approx 1.2$
.
In the secondary pipe, for case
$124out$
, a single wavenumber of
$k \approx 0.81$
is detected from linear mode
$1$
over the spectrum. From the DNS case, a brief, length-varying signal is detected between locations
$3.2D_1$
and
$4.5D_1$
ranging from
$k \approx 1.2$
to
$k \approx 4.0$
before the signal settles to
$k \approx 0.77$
. In both the linear and nonlinear cases, growth in the secondary pipes is far less prominent than what is observed in the primary pipe. For linear mode
$2$
, growth in the secondary pipe is undetectable (highlighted in figure 31).
From the DNS in the primary pipe a clear wavenumber of
$k \approx 0.29$
is apparent throughout the majority of the spectrum, slightly below what is detected from linear mode
$1$
(
$k \approx 0.32$
). A significantly higher wavenumber of
$k \approx 1.2$
is dominant in linear mode
$2$
, however, this is only detected intermittently in the DNS spectrum. In the secondary pipe, we also detect similar wavenumbers between the DNS and linear mode
$1$
. For case
$124out$
, where multiple unstable linear modes are detected, the streamwise spatial analysis indicates a close match between linear mode
$1$
and the saturated DNS flow, similar to what is observed from the frequency spectrum in figure 29.
Figure 32. Single-bifurcation outflow case,
$12out$
, streamwise spatial analysis. In all plots the spectrum is obtained from the FFT of the local axial component of the velocity,
$(U_{ax})$
, over a line orientated in the primary pipe. The line spans the coordinates
$x = 0.015$
,
$y = 0$
,
$z = -0.25$
to
$x = 0.15$
,
$y = 0$
,
$z = -5.25$
, orientating it parallel to the centreline and offset
$0.015$
in the
$x$
-direction from the centre of the pipe section. The
$x$
-axis is the streamwise position along the line in the primary pipe. The
$y$
-axis is the wavenumber
$(k)$
. The linear mode is displayed over an identical streamwise position and length as the corresponding mode displayed in the inset detail view in figure 15. The nonlinear plot is displayed over the same range (
$-0.25D$
to
$-3.2D$
). Panel (a) is the nonlinear solution from DNS at
${Re} = 4750$
. Panel (b) is linear mode
$1$
.
Figure 33. Double-bifurcation outflow case,
$124out$
, streamwise spatial analysis. In all plots the spectrum is obtained from the FFT of the local axial component of the velocity,
$(U_{ax})$
, over a line orientated in the primary pipe. The line spans the coordinates
$x = 0.3$
,
$y = 0$
,
$z = 1.0$
to
$x = 0.3$
,
$y = 0$
,
$z = -14.0$
, orientating it parallel to the centreline near the wall. The
$x$
-axis is the streamwise position along the line in the primary pipe. The
$y$
-axis is the wavenumber
$(k)$
. Each linear mode is displayed over an identical streamwise position and length as the corresponding mode displayed in the inset detail view in figure 16. The nonlinear plot is displayed over the entire range (
$1.0D$
to
$-5.2D$
). Panel (a) is the nonlinear solution from DNS at
${Re} = 6500$
. Panel (b) is linear mode
$1$
. Panel (c) is linear mode
$2$
.
The linear stability predicts that the outflow cases become unstable at critical Reynolds numbers of
${Re} \approx 4734$
and
${Re} \approx 6000$
for cases
$12out$
and
$124out$
, respectively. The DNS results in § 3.1 approximated critical Reynolds numbers of
${Re} \approx 4601$
and
${Re} \approx 5924$
for cases,
$12out$
and
$124out$
, respectively. The outflow comparisons are summarised in the last three rows of table 10.
5. Implications for respiratory flows
5.1. Medical ventilation
Respiratory flow applications are the motivation for this research. The frequencies at which our unstable modes oscillate are much higher than the oscillatory frequency of the flows encountered during regular breathing, and the frequencies currently used during high-frequency ventilation, which suggests that the current results of this paper using steady flows are potentially relevant. Further evidence for this comes from the studies of Jacob et al. (Reference Jacob, Tingay and Leontini2021, Reference Jacob, Tingay and Leontini2023), in which turbulent bursts were observed in reciprocating flows through geometries identical to those used here. In the double bifurcation, Jacob, Tingay & Leontini (Reference Jacob, Tingay and Leontini2023) reported frequencies of
$f \approx 0.23$
. In all cases, the turbulent bursts were shown to be a function of the instantaneous flow rate, suggestive of a quasi-steady instability mechanism.
In the double-bifurcation study by Jacob et al. (Reference Jacob, Tingay and Leontini2023), during the outflow portion of the cycle, the onset of turbulence was apparent near an instantaneous Reynolds number
${Re} \approx 6500$
(based on the instantaneous flow rate) in the primary pipe. Whereas during the inflow portion of the period, in the double bifurcation, the onset of turbulence was apparent near
${Re} \approx 7000$
. In both directions, the onset was localised within the primary and secondary pipes (during outflow), and the secondary and tertiary pipes (during inflow), before turbulence propagated through the domain.
In the clinical application of high-frequency ventilation, for the treatment of an infant patient, the peak instantaneous Reynolds number reaches
$\mathcal{O}(10^4)$
with an oscillatory frequency of ventilation in the range
$10{-}15\,\rm Hz$
(Scott et al. Reference Scott, Jacob, Tingay and Leontini2024; Hibberd et al. Reference Hibberd, Leontini, Scott, Pillow, Miedema, Rimensberger and Tingay2024). Scaling the frequencies that are observed from the leading unstable modes to an infant respiratory system with a tracheal diameter of
$5.5\,$
mm (Griscom 1982), using the definition of the Reynolds number in § 2.3, a non-dimensional frequency defined as
$\boldsymbol {f}=fD_1 / U^*$
and the kinematic viscosity of air
$\nu = 1.51 \times 10^{-5}$
, gives the dimensional frequencies listed in table 11.
Table 11. Dimensional frequencies of the leading unstable modes scaled to an infant patient respiratory system with a trachea diameter of
$5.5\,$
mm.
Figure 34. Double-bifurcation outflow case,
$124out$
, streamwise spatial analysis. In both plots the spectrum is obtained from the FFT of the local axial component of the velocity,
$(U_{ax})$
, over a line orientated in the secondary pipe. The line spans the coordinates
$x = 0.83$
,
$y = 0$
,
$z = 2.7$
to
$x = 3.52$
,
$y = 0$
,
$z = 7.0$
, orientating it parallel to the centreline approximately through the centre of the pipe. The
$x$
-axis is the streamwise position along the line in the secondary pipe. The
$y$
-axis is the wavenumber
$(k)$
. The linear mode is displayed over an identical streamwise position and length as the corresponding mode displayed in the inset detail view in figure 16. The nonlinear plot is displayed over the same range (
$3.2D$
to
$7.0D$
). Panel (a) is the nonlinear solution from DNS at
${Re} = 6500$
. Panel (b) is linear mode
$1$
.
Turbulent fluctuation, and subsequent mixing, is often considered to be of importance to gas transport in the airway, especially during high-frequency ventilation (Hibberd et al. Reference Hibberd, Leontini, Scott, Pillow, Miedema, Rimensberger and Tingay2024; Scott et al. Reference Scott, Jacob, Tingay and Leontini2024). Our results suggest that the mechanisms we observe from the linear stability analysis are potentially preserved above the critical threshold, before fully turbulent flows prevail, in at least three of the four cases which potentially provide a gas transport mechanism through the upper airway generations.
Some computational models of the airway network have attempted to include the impacts of turbulence (Herrmann, Tawhai & Kaczka Reference Herrmann, Tawhai and Kaczka2016), with parameterised models based on the turbulent transition in the reciprocating flow in a straight pipe. However, this study shows that turbulent onset is significantly different in the bifurcation geometry inherent to the airway with likely ramifications for turbulent mixing, as well as other transport mechanisms including mean streaming (Chang Reference Chang1984; Wanigasekara et al. Reference Wanigasekara, Jacob, Manasseh and Leontini2024).
We note that what is referred to as turbulent flow throughout this study, is not turbulence in the classic sense, i.e. large velocity fluctuations with multiple scale structures and an energy cascade. What is observed from DNS, just above the stability threshold, are in fact small fluctuations at very high frequencies (detailed in table 11). Despite these flows not being characterised by fully developed turbulence, there appears to be capacity for the flows to conditionally generate non-negligible Reynolds stresses, and by extension eddy viscosity and potentially enhanced mixing. Although the induced transport may only be a fraction of the speed of the underlying convective flow, this may be important over longer time scales for medical ventilation.
5.2. Limitations
The focus of this study is the linear instability of flows through idealised airway models using planar-symmetric geometry with a Hagen–Poiseuille velocity profile at the inlet to simulate the unidirectional flows. Our results indicate there is a fundamental linear instability mechanism that results in the transition from laminar to time-dependent, and subsequently turbulent, conditions when the Reynolds number exceeds the critical threshold and it is likely that this mechanism translates, to some extent, to the respiratory flows through more complex human systems. However, the linear mechanisms reported are likely not the entire picture when considering the complex dynamics of oscillatory flows through the human respiratory system. There are potentially other external perturbations which may contribute to instability mechanisms in the human airway. This study has not considered the presence of finite-amplitude perturbations from the upper or lower airway generations, or convective instabilities and non-modal mechanisms that may result in large transient growth of disturbances.
6. Discussion and conclusion
Our results confirm that the flows through complex curved pipe networks transition from laminar to temporally varying flows, and subsequently turbulence, via a linear mechanism associated with the longitudinal vortices, driven by a centrifugal imbalance when the flow passes through the curved regions of the geometry. This is in contrast to the nonlinear process that occurs in flow domains approaching the straight pipe limit. We verify the results from the linear stability, with DNS, and observe the frequencies of these slightly supercritical flows to be approximately equal to the frequencies detected from the linear stability analysis. The structures and frequencies of the unstable modes are completely dependent on the formation of the vortices, which themselves depend on the direction of the flow through the geometry and the number of bifurcations in the domain. For inflows, we observe a decrease in the critical Reynolds number with the addition of the tertiary bifurcation. However, the instability mechanism and mode structure in the secondary pipes are qualitatively similar, growing on two-vortex base flow topologies. In the tertiary pipes, in case
$124in$
, the base flow evolves to a four-vortex system where similarities in the mode structure with that seen in outflow case
$12out$
are apparent. The modes grow near the centre of the pipe, in regions of high strain rates, on the four vortices. For outflows, we observe an increase in the critical Reynolds number with the addition of the tertiary bifurcation. More importantly, the tertiary bifurcation triggers a different instability mechanism which arises due to a modification of the secondary base flow structure. We track this modification to two regular SSBBs in the (natural) temporally stable base flow. These results are potentially important for understanding and modelling gas transport through the upper generations of the human respiratory system.
Acknowledgements
T.J.A.S. acknowledges the Australian Government and Swinburne University of Technology for the provision of a Research Training Program Scholarship. T.J.A.S. and J.S.L. acknowledge the support of resources and services from the National Computational Infrastructure (NCI) and the Pawsey Supercomputing Research Centre’s Setonix Supercomputer (https://doi.org/10.48569/18sb-8s43), via competitive grant IZ4 which is supported by the Australian Government as well as the OzSTAR national facility at Swinburne University of Technology which is funded by Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government, and from the Victorian Higher Education State Investment Fund (VHESIF) provided by the Victorian Government.
Declaration of interests
The authors report no conflict of interest.
Appendix A
A.1 Grid resolution study
A grid resolution study was carried out to confirm the accuracy of the results. The grid resolution study was performed by monitoring both global and local parameters. Globally, the energy norm based on the average of the root-mean-square of the velocity across the solution domain was tracked and the linear modes from the stability analysis eigendata were inspected. Locally, the structure of the flow at a section downstream of the bifurcation (into the secondary pipe) that is known to be susceptible to the onset of instabilities was inspected in the base flow and the structure of an unstable mode was inspected where growth occurred.
For the single-bifurcation geometry, the first global and local resolution studies, relating to the base flow, were run across polynomial orders (p-order) of the basis functions of the spectral elements of
$\mathrm{8,10,12,14}$
and
$\mathrm{16}$
at
${Re} = 9600$
. The results of the global resolution study showed a decreasing global energy norm as the p-order was increased. An exponential curve, expressed as a function of the p-order (
$E_{norm} = 0.33e^{-0.2011{p-order}} + 0.5663$
), was fitted to the measured global data to estimate a converged energy norm. Using this approximation for convergence (
$\mathrm{0.5663}$
), the errors associated with each p-order were estimated (detailed in table 12).
Table 12. Results from the first global grid resolution study relating to the global energy norm. Results show the error associated with p-orders
$8,10,12,14$
and
$16$
at
${Re} = 9600$
under inflow conditions. The reported values of
$E_{norm}$
were obtained at time,
$\mathrm{\tau = 115}$
. The last column, CPU hours, indicates the total wall clock time required to run each case from the study.
For the local indicator of convergence, the flow structure located
$\mathrm{5.5}$
diameters downstream of the bifurcation at time,
$\mathrm{\tau = 100}$
was monitored. By p-order
$\mathrm{12}$
the main features of the flow were resolved and changes in the structure appeared qualitatively diminished with respect to subsequent increases in the p-order (shown in figure 35).
Figure 35. Results from the first local grid resolution study. All pipe sections shown are located downstream of the bifurcation in the secondary pipe. Results from all p-orders are shown as the velocity magnitude at
${Re} = 9600$
under inflow conditions.
Based on the energy norm convergence criterion, the error in the solution for p-order
$\mathrm{12}$
was less than
$\mathrm{6 \,\%}$
. The maximum critical value for the Reynolds numbers investigated in this study was less than
$\mathrm{7750}$
, equating to
$\mathrm{\approx 20 \,\%}$
less than the Reynolds number used in the grid convergence study, suggesting the error associated with p-order
$\mathrm{12}$
is likely to be less than the reported estimation. The local indicator of convergence showed that the structure of the flow was resolved by p-order
$\mathrm{12}$
, giving confidence that the important features were captured. Refining the grid resolution further required significantly greater computational resources making this difficult. To decrease the estimated error to below
$2\,\%$
required changing the p-order from 12 to 16, resulting in three times more wall clock time (based on this study), rendering the improvement in accuracy not worth the computational expense. Subsequently, interpolation polynomial order
$12$
was used to establish all base flows reported in this study.
The second global and local resolution studies considered the linear modes outputted from the Arnoldi algorithm and the local structure in the flow where the growth was evolving. Based on the global and local resolution results, relating to the base flow solution, the base flow for calculating the leading modes was p-order
$12$
. From this base flow, the perturbations were evolved forward in time using p-orders
$12, 14$
and
$16$
for the linear solver. Three Arnoldi calculations were performed, for each linearised run, to confirm the leading modes being tracked were the same (detailed in table 13). Based on the value of the multipliers (relating to the growth rate), the error across the three interpolation polynomial orders was
$\lt\!1 \,\%$
.
Table 13. Results of the global grid resolution study relating to the unstable modes from interpolation polynomial orders
$12, 14$
and
$16$
calculated from the interpolation polynomial order
$12$
base flow at
${Re} = 9600$
.
Locally, the growth of the leading modes was compared to confirm the structure across the three resolutions relating to the linear solver. The mode structures were fully resolved by p-order
$\mathrm{12}$
. Subsequently, all simulations (base flow and perturbation) were run with p-order
$12$
, giving a total node count of
$27\,684\,800$
.
Figure 36. Results from the second local grid resolution study. All pipe sections are located downstream of the bifurcation in the secondary pipe. Results from all p-orders are shown as the velocity magnitude at
${Re} = 9600$
under inflow conditions.
For the double-bifurcation geometry, we have used a p-order of
$12$
for all base flow and perturbation calculations, based on the resolution study detailed in Jacob et al. (Reference Jacob, Tingay and Leontini2023), who used an identical macro-mesh. Subsequently, all simulations (base flow and perturbation) were run with a total node count of
$66\,443\,520$
.
Appendix B
B.1 Selective frequency damping framework convergence study
Figures 37 and 38 show the results from a convergence test of the SFD framework from the single-bifurcation inflow case,
$12in$
, at
${Re} = 7250$
. The test tracked the difference in the velocity magnitude time history from two set-ups: a DNS solution with the SFD framework switched off, and a DNS solution with the SFD framework switched on. The difference was correlated with the global
$L2$
norm outputted from the SFD framework. In figure 37, the undamped velocity magnitude time history is plotted from a probe point located
${\approx}\, 15$
diameters downstream of the origin in the secondary pipe of the DNS solution in black, which stabilises naturally given sufficient time. Plotted in red is the solution after the SFD framework has been switched on, from the same probe, with all other parameters remaining identical. In this validation case, the SFD framework was switched on at time
$\tau = 75$
, just prior to the epoch of time where the turbulence appeared most prominent, and run until a selected convergence criterion was reached. The final values of the dimensionless velocity magnitude shown in figure 37 were both equal to
$0.609119$
correlating with an root-mean-square error (
$L2$
) norm of
$10^{-9}$
confirming a convergence criterion accurate to six or more significant figures. Figure 38 shows the truncated convergence history plotting the
$L2$
norm, outputted by the SFD algorithm, correlating with the red trace as the solution approaches the true steady state.
This test case was known to stabilise naturally given enough time making it suitable for testing the accuracy of the algorithm and establishing an adequate convergence criterion. The steady-state base flow solutions were considered converged when the global
$L2$
norm
$\leq 10^-9$
. These converged base flow solutions formed the basis around which the governing equations were linearised and subsequently used for all linear stability calculations.
Figure 37. SFD framework convergence test at
${Re} = 7250$
from the single-bifurcation inflow case,
$12in$
. The black line tracks the velocity magnitude of the undamped DNS solution. The red line tracks the velocity magnitude of the same case after the SFD framework has been switched on at
$\tau =75$
. All other parameters are identical. The probe point from which the data was collected is located
${\approx}\, 15$
diameters downstream of the origin in the centre of the secondary pipe. The time history of the SFD solution is correlated with an
$L2$
norm time history shown in figure 38. The SFD solution used a filter gain of
$0.15$
and filter width of
$3.2$
. The initial ramping up velocity data (
$\tau \lt 35$
) from the DNS trace has been omitted for clarity.
Figure 38. The SFD framework convergence test at
${Re} = 7250$
from the single-bifurcation inflow case,
$12in$
. The three traces are outputted from the SFD algorithm plotting the global
$L2$
norm from each grid point after the forcing was switched on at
$\tau = 75$
. The time history of the
$L2$
norm is correlated with the time history of the SFD solution (red trace) in figure 37. The SFD solution used a filter gain of
$0.15$
and filter width of
$3.2$
.
Appendix C
C.1 In-house Arnoldi algorithm validation
A validation study was carried out to verify the accuracy of the results obtained from our in-house Arnoldi algorithm. We compare the results generated by our in-house code with the study detailed in Lupi et al. (Reference Lupi, Canton and Schlatter2020). We develop an identical numerical set-up to replicate the conditions for unidirectional flow through a
$90$
degree bent pipe and compare the eigenvalues and frequency of the unstable eigenmode, as well as the mode structure. The data reported from our in-house code were run at
${Re} = 2550$
, the same as that reported in Lupi et al. (Reference Lupi, Canton and Schlatter2020).
Our algorithm successfully identifies a pair of unstable complex conjugate eigenvalues which are in excellent agreement with the values reported by Lupi et al. (Reference Lupi, Canton and Schlatter2020), detailed in table 14.
Figure 39 presents a direct comparison between the computed eigenmode structure from our in-house code (left) and the reference structure reported by Lupi et al. (Reference Lupi, Canton and Schlatter2020) (right). The spatial structure of the eigenmode computed from our code reproduces the key features observed by Lupi et al. (Reference Lupi, Canton and Schlatter2020).
Table 14. Results of the validation study comparing the eigenvalues generated from our in-house Arnoldi algorithm with the results reported in Lupi et al. (Reference Lupi, Canton and Schlatter2020) from the uni-directional flow through the
$90$
degree bent pipe.
Figure 39. Comparison of the unstable eigenmode structure obtained from our in-house Arnoldi algorithm (panel (a)) and the reference results reported by Lupi et al. (Reference Lupi, Canton and Schlatter2020) (panel (b)). Both eigenmode structures are shown at
${Re} = 2550$
. In the plane of symmetry, pseudocolours of the outward velocity component are shown. In the pipe cross-sections, in-plane streamlines and pseudocolours of the streamwise velocity component are shown. The heat map is arbitrary units. The arrow indicates direction of flow through the domain. Inset is a magnified view of the mode structure inside the bend.
Open access
