4.1. EIT in a single-phase flow

Figure 2 summarises the flow states in the ${Re}$ – ${\textit{Wi}}$ space. At a low value of ${Re}=10$ , the single-phase flow remains fully laminar for up to ${\textit{Wi}}=1000$ . Two points are selected in the computational domain to collect all the components of velocity and viscoelastic stresses at each time step to monitor their evolution to ascertain the attainment of a statistically steady-state solution. One of the ‘numerical probes’ is located at the centre of the channel ( $x^*$ = $1$ , $y^*$ = $6$ , $z^*$ = $1$ ), while the other one at $25\,\%$ distance ( $x^*$ = $0.5$ , $y^*$ = $6$ , $z^*$ = $1$ ) from one of the walls. Once the flow reaches a statistically steady state, the simulations are continued at least for another $3\lambda$ time units, during which, full flow fields are stored for a further statistical analysis.

Figure 2. Various flow states in the ${Re}$ – ${\textit{Wi}}$ space ( $ Ca=0.01, \beta =0.1$ ).

The simulations are then repeated for ${Re}=100$ . While gradually increasing ${\textit{Wi}}$ , it is observed that the flow remains laminar until ${\textit{Wi}}=100$ . However, for $Wi \geqslant 100$ , instabilities start to appear in this low- ${Re}$ flow. These instabilities start to grow with time and fluctuations are observed ultimately in all three components of velocity akin to a typical turbulent signal. Viscoelastic stresses also start to fluctuate in the time domain, as collected by the numerical probes at two different locations of the channel. It is important to emphasise that no explicit perturbations are provided from outside the system to initialise any instabilities in the flow. Only for $Wi\geqslant 100$ , the viscoelastic stresses are high enough to trigger an instability that ultimately leads to a chaotic flow regime. As the flow rate is kept constant for a fixed ${Re}$ , a sudden increase in the applied pressure gradient is also observed once the flow is transitioned from laminar to this chaotic state, indicating a rapid increase in the drag. The characteristics of this chaotic flow regime for $Wi \geqslant 100$ will be discussed in detail in the subsequent paragraphs.

When ${Re}$ is increased further to $1000$ , the transition to a chaotic flow regime starts to appear at a lower value of $Wi \geqslant 50$ . This ${Re}$ is still smaller than the minimum ${Re}_{cr} \approx 1100$ required for sustaining the inertial turbulence in a channel flow (Owolabi et al. Reference Owolabi, Poole and Dennis2016). As ${\textit{Wi}}$ is increased further, the turbulent intensity of this chaotic regime increases with a change in statistical quantities as well.

Figure 3(a) shows the velocity vectors in a vertical cutting $xz$ -plane of the viscoelastic duct flow once it reaches a statistically steady state for the ${Re}=1000$ and ${\textit{Wi}}=1000$ case. The instantaneous and time-averaged flow fields are plotted in the left and right panels of figure 3(a), respectively. As seen, the instantaneous flow field exhibits a chaotic flow regime, while the time-averaged velocity field shows an almost similar secondary flow pattern as also observed in the viscoelastic laminar duct flow with the shear-thinning effect (Li, McKinley & Ardekani Reference Li, McKinley and Ardekani2015; Naseer et al. Reference Naseer, Izbassarov, Ahmed and Muradoglu2024). This secondary flow field is developed in the duct flow primarily due to the presence of a second normal stress difference $N_2 = \bar \tau _{zz} - \bar \tau _{xx}$ . The instantaneous and time-averaged contour plots of $N_2$ are shown in the left and right panels of figure 3(b), respectively. The instantaneous plot confirms the chaotic flow pattern, while the time-averaged one resembles the typical distribution of $N_2$ in a laminar duct flow. As the EIT regime is associated with long-elongated sheets of the first normal difference ( $N_1$ ), the distribution of $N_1$ is also shown in figure 3(c) in the $yz$ -plane at the centre of the duct for the various combinations of ${Re}$ and ${\textit{Wi}}$ for which the single-phase flow becomes chaotic. Dubief et al. (Reference Dubief, Terrapon and Hof2023) categorised different EIT regimes based on the structures of $N_1$ for a 2-D channel flow. For a low value of ${Re}=100$ , a similar ‘chaotic regime’ is observed in the present 3-D duct flow as shown in figure 3(c). At a higher value of ${Re}=1000$ , but with a moderate value of ${\textit{Wi}}=100$ , this regime is shifted to ‘chaotic arrowhead’ and once the Weissenberg number is increased to a very high value of ${\textit{Wi}}=1000$ , the long-elongated sheets of $N_1$ break down completely and the flow starts to resemble a typical inertial turbulence. The iso-surfaces of the second invariant of the velocity gradient tensor ( $Q$ -criterion) is depicted in figure 3(d) to show the corresponding change in the near-wall vortices for the same combinations of ${Re}$ and ${\textit{Wi}}$ . At a very high value of ${Re}$ and ${\textit{Wi}}$ , the positive and negative contours of the $Q$ -criterion show a chaotic pattern near the wall. With decreasing ${Re}$ or ${\textit{Wi}}$ , the spanwise elongated pattern of positive and negative contours of the $Q$ -criterion starts to appear indicating that the signatures of elastic turbulence become increasingly more apparent in the EIT regime once the elastic effects dominate over the inertial effects. In the absence of any perturbations from outside the system, the flow is laminar at a lower value of ${\textit{Wi}}$ in the present scenario, as shown in figure 2. Therefore, a clear pattern of alternating contours of $Q$ -criterion is not visible here, as also observed by Dubief et al. (Reference Dubief, Terrapon and Hof2023).

Figure 3. (a) Instantaneous (left) and averaged (right) snapshots of velocity vectors showing the secondary flow field at ${Re}=1000$ and ${\textit{Wi}}=1000$ . (b) Instantaneous (left) and time-averaged (right) snapshots of second normal stress difference ( $N_2$ ) at ${Re}=1000$ and ${\textit{Wi}}=1000$ . (c) Contours of first normal stress difference ( $N_1$ ) are shown in a statistically steady EIT regime for different values of ${Re}$ and ${\textit{Wi}}$ . (d) Iso-surfaces of $Q$ -criterion are shown close to the bottom wall of the channel in a statistically steady EIT regime for different values of ${Re}$ and ${\textit{Wi}}$ . The red and grey colours represent the positive and the negative values, respectively. The contours are plotted at $\pm {5}$ , $\pm {0.05}$ and $\pm {0.005}$ for the ( ${Re}, Wi$ ) = ( $1000, 1000$ ), ( $1000, 100$ ) and ( $100, 1000$ ) cases, respectively.

Figure 4 shows various quantities used to evaluate turbulent characteristics of this single-phase chaotic flow regime for various combinations of ${Re}$ and ${\textit{Wi}}$ . The first one is the 1-D energy spectrum (figure 4 a). The energy axis is normalised by $({1}/{2})\rho _o u_o^2$ and the frequency axis by $u_o/h$ . At the higher values of $({\textit{Re}},{\textit{Wi}})= (1000, 1000)$ , the kinetic energy spectrum shows a typical $-5/3$ scaling at the large scales but it shifts to $-4$ at the smaller scales (dissipation range). These values are consistent with the results obtained in a single-phase channel flow at lower values of ${Re}$ by Mukherjee et al. (Reference Mukherjee, Singh, James and Ray2023) and Foggi Rota et al. (Reference Foggi Rota, Amor, Le Clainche and Rosti2024). The elasticity number, $El = Wi/{Re}$ , provides a measure of the relative strength of elastic to inertial effects. At $({\textit{Re}},{\textit{Wi}}) = (1000,1000)$ , $El = 1$ , indicating a balanced elasto-inertial regime. The transition to a $-4$ slope at smaller scales indicates that viscoelastic effects become significant as the length scale decreases. The steep slope suggests strong suppression of small-scale fluctuations due to elastic stresses in this highly concentrated polymer solution ( $\beta = 0.1$ ). At large scales, the flow behaves like classical turbulence, driven by inertial forces, as seen in figure 3(c). At smaller scales, high ${\textit{Wi}}$ and low $\beta$ cause elastic stresses to dominate, altering the cascade. This dual scaling suggests a cross-over frequency where the flow transitions from inertia-dominated to elasticity-dominated dynamics. The exact frequency depends on the balance among ${Re}$ , ${\textit{Wi}}$ and the polymer properties. Once ${\textit{Wi}}$ is decreased to ${\textit{Wi}}=100$ while keeping ${Re}=1000$ , the scaling of energy spectrum shifts from $-4$ to $-3$ across the majority of the spectrum (figure 4 a). It indicates that lowering ${\textit{Wi}}$ reduces the elastic effects, allowing inertial forces to distribute energy more effectively to smaller scales, thus resulting in a less steep spectrum. The $-3$ slope also suggests a transitional EIT regime where the inertial and the elastic effects are not so balanced ( $El = 0.1$ ), with less suppression of small-scale fluctuations compared with the $-4$ slope at $Wi = 1000$ . This is characterised by the ‘chaotic arrowhead’ regime in figure 3(c). Once ${Re}$ is reduced to $100$ while keeping ${\textit{Wi}}=1000$ , the energy spectrum again yields a slope of $-3$ . Here, low ${Re}$ weakens the inertial effects, increasing the relative importance of viscous and elastic forces ( $El = 10$ ). The $-3$ slope suggests a regime closer to ET, where elastic instabilities drive the flow, modulated by viscous dissipation. The similarity in the $-3$ slope between $({\textit{Re}},{\textit{Wi}}) = (1000,100)$ and $({\textit{Re}},{\textit{Wi}}) = (100,1000)$ indicates that different combinations of inertial and elastic effects can produce similar spectral characteristics. It can be inferred from these scales that once the turbulent flow is entirely governed by inertia (classical inertial turbulence), the energy spectrum assumes a slope of $-5/3$ , whereas once the turbulent regime is dominated by the elastic effects with negligible inertia, as in ET, the spectrum exhibits a slope of $-4$ . Any combination of ${Re}$ and ${\textit{Wi}}$ in between these two regimes assumes a slope in between these two limits.

Figure 4. (a) One-dimensional energy spectrum for the $({\textit{Re}},{\textit{Wi}})=(1000,1000)$ , $(1000,100)$ and $(100,1000)$ cases. (b) Distribution of the components of the average shear stress from the channel wall towards the centre in the mid-plane for $({\textit{Re}},{\textit{Wi}})=(1000,1000)$ (top panel) and $({\textit{Re}},{\textit{Wi}})=(1000,100)$ (bottom panel). (c, d, e) Distribution of turbulent kinetic energy in the mid-plane of the channel for different values of ${Re}$ and ${\textit{Wi}}$ . (f, g, h) Average elongation of polymer molecules represented by $tr(\bar B)$ shown in the same mid-plane.

The total mean shear stress in the mid-plane comprises three components: the viscous shear due to mean flow $\bar \tau _{{N}}$ , the Reynolds stress due to velocity fluctuations $\bar \tau _{{R}}$ and a polymeric stress $\bar \tau _{{P}}$ . Once the flow reaches a statistically steady state, the applied pressure gradient must balance the total shear stress at the four walls of the channel. The three components of shear stress ( $\bar \tau _{{N}} = \mu _o ({\partial \bar v}/{\partial z})$ , $\bar \tau _{{R}} = -\rho _o \overline {v'w'}$ , $\bar \tau _{{P}} = \bar \tau _{yz}$ ) are plotted for the central plane of the duct ( $x^*=1$ ) for two values of ${\textit{Wi}}=1000$ and ${\textit{Wi}}=100$ at ${Re}=1000$ (figure 4 b). The fluctuating and the time-averaged quantities are denoted by prime $(.')$ and overbar $(\bar {.})$ , respectively. Note that the averaging is performed both in time and in streamwise direction since the flow is expected to be homogeneous in the streamwise direction once a statistically steady state is reached. The stress components are normalised by the local shear stress at the wall of the respective plane. For the high value of ${\textit{Wi}}=1000$ , it is observed that the stress balance in the mid-plane of the channel resembles that of a classical inertial turbulence. The viscous stress is maximum at the wall and decays to zero at the channel centre, whereas the Reynolds stress is zero at the wall, gradually increases to its peak value and then goes to zero at the channel centre. The viscoelastic stress also follows the viscous stress profile, and becomes maximum at the wall and zero at the channel centre. As a result, the total shear stress decays linearly from the wall towards the channel centre (figure 4 b). Once ${\textit{Wi}}$ is reduced to $100$ while keeping ${Re}=1000$ , the magnitude of Reynolds stress becomes negligibly small compared with viscous and polymeric stresses. Although the $-3$ scaling of the energy spectrum and the turbulent signal are indicative of a turbulent flow field at this comparatively lower value of ${\textit{Wi}}=100$ , the Reynolds stress becomes negligibly small compared with the overall shear stress in the mid-plane. These not-so-smooth profiles of different shear stress components indicate a clear departure of the turbulent flow field from the inertial one. The distribution of the turbulent kinetic energy (TKE) $\mathcal K$ $= 0.5(\overline {u'^2} + \overline {v'^2} + \overline {w'^2})$ is shown in figures 4(c), 4(d) and 4(e) at the mid-plane where TKE is found to be maximum near the walls and decays to a minimum value at the centre. Interestingly, the average elongation of polymer molecules quantified by the trace of the conformation tensor $tr(\bar {B})$ shows a similar trend for different values of ${\textit{Wi}}$ , as shown in figures 4(f), 4(g) and 4(h). It is worth mentioning that the value of the applied pressure gradient decreases as ${\textit{Wi}}$ increases for a constant value of ${Re}$ , indicating a decrease in the drag force. As the chaotic flow is driven by elasticity, this unexpected drag reduction with an increase in ${\textit{Wi}}$ is attributed to the shear-thinning effect that dominates the drag increase caused by the pronounced chaotic flow.

Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) have shown that for the large Weissenberg numbers, the viscoelastic pipe flow can be linearly unstable to axisymmetric disturbances down to a critical Reynolds number of ${Re}_{cr} \approx 63$ . Wan, Sun & Zhang (Reference Wan, Sun and Zhang2021) and Buza, Page & Kerswell (Reference Buza, Page and Kerswell2022) later demonstrated that the polymer concentration ( $\beta$ ) is the main controlling factor. At a higher polymer concentration, the unstable region shrinks and the transition is mostly supercritical. A similar mechanism seems to be at play in the present single-phase duct flow, where transition occurs even in the absence of any finite perturbation at a high Weissenberg number. The subcritical transition does not occur in a single-phase flow due to the stabilising role of high polymer concentration ( $\beta =0.1$ ), as seen for the lower Weissenberg number cases for all the three Reynolds numbers.

Foggi Rota et al. (Reference Foggi Rota, Amor, Le Clainche and Rosti2024) have demonstrated that a single-phase viscoelastic flow can sustain a chaotic flow state in a straight channel at a Reynolds number as low as ${Re}=0.5$ ( ${\textit{Wi}}=50$ ) if started from a turbulent initial condition for a dilute polymer solution by using the FENE-P model with $\beta =0.9$ . The value of ${\textit{Wi}}=50$ was kept constant in their simulations while bringing down the Reynolds number from ${Re}=2800$ to ${Re}=0.5$ to keep the polymers in their stretched state. To observe whether the chaotic state is sustained in the present scenario of a duct flow with a concentrated polymer solution, extensive simulations are next performed by keeping the Reynolds number constant at ${Re}=1000$ and the Weissenberg number is reduced from ${\textit{Wi}}=1000$ to [ $500, 100, 50, 40, 10$ ]. Each of the lower ${\textit{Wi}}$ is started by using the initial conditions of the previous higher ${\textit{Wi}}$ case. It is found that the chaotic flow state is sustained until ${\textit{Wi}}=40$ . When ${\textit{Wi}}$ is further reduced to $10$ , the single-phase flow is laminarised. Next, the same simulations are repeated by fixing the Reynolds number at ${Re}=100$ . The Weissenberg number is reduced gradually from ${\textit{Wi}}=1000$ until $10$ and it is observed that the chaotic flow is sustained until $Wi\geqslant 40$ . The same limit of $Wi\geqslant 40$ is observed for the ${Re}=10$ cases to sustain a chaotic regime in this duct flow if perturbed initially. It can be concluded that if the flow is perturbed initially, the chaotic flow is sustained in this single-phase viscoelastic duct flow only if the value of $Wi\geqslant 40$ . Below this value, the single-phase flow is laminarised. This value is close to ${\textit{Wi}}=50$ used by Foggi Rota et al. (Reference Foggi Rota, Amor, Le Clainche and Rosti2024) in their channel flow simulations.

It is important to highlight that exceptionally high Weissenberg numbers required to trigger flow instability in the absence of any discrete external perturbations in this rectilinear channel may seem impractical. However, it should be noted that ${\textit{Wi}}$ is defined based on the average flow velocity in the channel. As shear rate decreases significantly towards the channel centre in this duct flow, the effective local Weissenberg number is actually much lower than the nominal value of ${\textit{Wi}}$ , which is based on the average flow velocity.

4.1.1. Turbulent kinetic energy

The turbulent kinetic energy ( ${\mathcal K}=({1}/{2}) \overline {u^{\prime}_i u^{\prime}_i}$ ) in a viscoelastic flow evolves by

(4.1) \begin{align} \frac {\partial {\mathcal K}}{\partial t} &= \overbrace {-\bar u_j \frac {\partial {\mathcal K}}{\partial x_{\!j}}}^{\mathcal A} - \overbrace {\frac {1}{2}\frac {\partial (\overline {u_i^{\prime} u_i^{\prime} u^{\prime}_j})}{\partial x_{\!j}}}^{\mathcal Q} - \overbrace {\frac {1}{\rho _o} \overline {u^{\prime}_i \frac {\partial p'}{\partial x_i}}}^{\mathcal R} - \overbrace {\left (\overline {u^{\prime}_i u^{\prime}_j} \frac {\partial \bar u_i}{\partial x_{\!j}}\right )}^{\mathcal P} + \overbrace {\frac {1}{2}\frac {\partial }{\partial x_{\!j}}\left(\nu _s \frac {\partial \overline {u^{\prime}_i u_i^{\prime}}}{\partial x_{\!j}}\right)}^{\mathcal D} \nonumber \\ &\quad - \underbrace {\nu _{s} \overline {\frac {\partial u^{\prime}_i}{\partial x_{\!j}}\frac {\partial u^{\prime}_i}{\partial x_{\!j}}}}_{{\mathcal \epsilon }} +\underbrace {\frac {1}{\rho _o}\left (\overline {u^{\prime}_i\frac {\partial \tau ^{\prime}_{ij}}{\partial x_{\!j}}}\right )}_{{\mathcal W}_p} + \underbrace {\overline {u^{\prime}_i f^{\prime}_i}}_{\mathcal B} , \end{align}

where $\mathcal A$ , $\mathcal Q$ , $\mathcal R$ , $\mathcal {P}$ , $\mathcal {D}$ , $\mathcal \epsilon$ , ${\mathcal W}_{p}$ and $\mathcal {B}$ represent advection by mean flow, transport by velocity fluctuations, transport by pressure, production by mean flow, viscous diffusion, viscous dissipation, polymer work and work done by body force, respectively. In a viscoelastic turbulent flow, ${\mathcal W}_{p}$ may change its sign and can serve as either dissipation or production depending upon the signs of the polymer stress fluctuations and fluctuating velocity gradients, as has been observed in both experimental (Ptasinski et al. Reference Ptasinski, Nieuwstadt, van den Brule and Hulsen2001) and numerical (Dallas, Vassilicos & Hewitt Reference Dallas, Vassilicos and Hewitt2010) works.

Figure 5. Contours of (a) advection by mean flow, (b) transport by velocity fluctuations, (c) transport by pressure, (d) production by mean flow, (e) viscous diffusion, (f) viscous dissipation, (g) polymer work and (h) time derivative of TKE are shown in a vertical cutting $xz$ -plane for the (left) ${Re}=1000, {\textit{Wi}}=1000$ and (right) ${Re}=100, {\textit{Wi}}=1000$ cases, respectively.

To get better insight of the EIT regime in this duct flow, these terms are averaged in time and in the streamwise ( $y$ ) direction once the flow reaches a statistically steady state. The constant contours of all the terms in (4.1) normalised by $v_\tau ^2$ are shown in figure 5 at two different values of ${Re}$ in a vertical cutting $xz$ -plane except for the body force term that is zero for this single-phase flow. For ${Re}=1000$ , the advection by mean flow (figure 5 a-left) shows small regions of positive peaks at the four corners and towards the centre of the duct. A pattern of four negative peaks is also observed slightly away from the walls at this high-Reynolds-number EIT flow. However, once the Reynolds number is lowered to ${Re}=100$ , the magnitude of the advection term becomes negligibly small. Similarly, a symmetric pattern of transport term is observed for the higher Reynolds number (figure 5 b-left) and the same term becomes negligibly small at lower inertia. At ${Re}=1000$ , the positive and negative peaks of the transport term remain closer to the channel walls except at the corners where it remains approximately zero. An opposite trend is observed for the pressure term (figure 5 c). Once the flow is dominated by the inertial effects ( ${Re}=1000$ ), the pressure term remains significant while, in an elastically dominated EIT regime ( ${Re}=100$ ), it shows a symmetric pattern of positive and negative zones, however, with a negligibly small magnitude. The positive peaks of the pressure term extend from four corners of the duct until the centre with negative peaks in the vicinity of the channel walls. Turbulence production by the mean flow is mainly produced near the walls except at the four corners (figure 5 d) and becomes zero in most of the central region of the duct with a symmetric pattern at ${Re}=1000$ . However, in the elastically dominated EIT flow ( ${Re}=100$ ), this pattern changes with four regions of positive peaks away from the walls. As the turbulence is mainly produced by the elastic effects at ${Re}=100$ , the magnitude of the production term also remains negligibly small. The diffusion term (figure 5 e) shows positive peaks adjacent to the walls except at the corners, immediately followed by the negative peaks before decaying to zero for ${Re}=1000$ whereas, for the ${Re}=100$ case, it is positive near the walls and remains chaotic in most of the remaining central portions of the duct. The dissipation term shows a similar pattern for both cases (figure 5 f), i.e. it is maximum at the walls except in the corners and becomes zero towards the duct centre. The most important term in this EIT regime is the work done by the elastic force (figure 5 g) as it can contribute towards both generation or dissipation of turbulence. At higher inertia ( ${Re}=1000$ ), the positive peak of the polymer work is observed to be closer to the four walls of the channel. It indicates that the viscoelastic stresses contribute to the turbulence production closer to the walls, but act as a dissipating sink away from it. The same term remains negligible at the centre. Interestingly, at lower inertia ( ${Re}=100$ ) where the elastic effects dominate, the polymer stresses contribute towards turbulence dissipation at the four corners and towards the central region of the duct. The turbulence production by polymer stresses is concentrated near the walls away from the corners. The viscous dissipation and the pressure terms follow the similar pattern, as seen in figures 5(f) and 5(c), respectively. Finally, the summation of all the terms in (4.1) approaches zero confirming that a statistically steady state is reached, as shown in figure 5(h).

Notably, the distribution of various components of TKE in the present EIT regime of a highly concentrated polymer solution at low ${Re}$ is markedly different than that in a dilute viscoelastic turbulent duct flow at high ${Re}$ as studied by Shahmardi et al. (Reference Shahmardi, Zade, Ardekani, Poole, Lundell, Rosti and Brandt2019), where the turbulent statistics of the flow were entirely dominated by the inertial effects.

4.2. EIT in a multiphase flow

For a given Reynolds number, a significantly high ${\textit{Wi}}$ is always required to achieve the EIT state in a single-phase flow in the absence of any external perturbations, as has been shown in the previous section. When ${Re}$ is low (e.g. ${Re}=10$ ), even a Weissenberg number as high as $1000$ cannot make the flow unstable (figure 2). Similarly, the single-phase flow is found to be essentially laminar at the lower values of ${\textit{Wi}}$ even when the Reynolds number is as high as ${Re}=1000$ .

We next examine the effects of bubble injection into the viscoelastic laminar channel flow. For this purpose, simulations are performed for three relatively low values of ${\textit{Wi}}=1$ , $5$ and $10$ at each value of ${Re}=10$ , $100$ and $1000$ . Note that the single-phase flow remains laminar for these combinations of ${\textit{Wi}}$ and ${Re}$ . Calculations are first carried out for the single-phase flow until a statistically steady state is reached. Then, spherical bubbles are injected instantaneously into the flow with a volume fraction of $3\,\%$ and a relative bubble size of $d_b/2h=0.2$ , where $d_b$ is the bubble diameter. The bubbles are initially injected randomly and uniformly in the entire duct. The capillary number is fixed at $Ca=0.01$ to keep the bubbles nearly spherical. The state of this multiphase flow field is continuously monitored at each time step using the two numerical probes as for the single-phase case. By introducing the bubbles into the flow, the flow is found to achieve a chaotic state even for Reynolds and Weissenberg numbers as low as ${Re}=10$ and ${\textit{Wi}}=5$ . Time histories of instantaneous velocity components are plotted in figure 6(a) for the ${Re}=1000$ and ${\textit{Wi}}=5$ case to show the transition from a laminar to an EIT state after the injection of bubbles. Bubbles are injected at approximately $t^*=400$ . As seen, the signals are very smooth until the injection of bubbles indicating a laminar flow. Once the bubbles are injected into the flow, all three components of the velocity exhibit random fluctuations indicating a transition to the EIT state. These fluctuations propagate throughout the domain over time, ultimately leading to a fully chaotic flow field. This chaotic flow state is visualised by plotting constant contours of vorticity ( $\boldsymbol{\omega }=\boldsymbol{\nabla }\times \boldsymbol{u}$ ) magnitude in figure 6(b) for the single-phase and multiphase regimes in a vertical cutting $xz$ -plane at the centre of the duct. Once the flow reaches a statistically steady state after the injection of bubbles, the simulations are continued at least for another $10$ flow-through time units ( $10\times 12h/u_o$ ) to collect data for further statistical analysis.

Figure 6. (a) Evolution of flow velocity and its transition to turbulent state once bubbles are injected into the flow. (b) Contours of vorticity magnitude in an $xz$ -plane at the centre of the domain for the single-phase and multiphase regimes ( ${Re}=1000, {\textit{Wi}}=5, \beta =0.1, Ca=0.01$ ).

Although not shown here due to space considerations, further simulations are performed to demonstrate the effect of bubble-induced streamline curvature on the destabilisation of the flow. For this purpose, simulations start from a two-phase Newtonian laminar flow at ${Re} = 100$ . Subsequently, the flow is made weakly viscoelastic by setting a small value of $Wi = 0.5$ and the simulation is continued until the viscoelastic stresses are fully developed and no flow destabilisation is observed. Then, the Weissenberg number is gradually increased first to ${\textit{Wi}}=1$ , then to ${\textit{Wi}}=2$ and eventually to ${\textit{Wi}}=5$ , while the flow state is continuously monitored. In this process, the simulations are carried out until a statistically steady state is reached before increasing ${\textit{Wi}}$ in each step. No transition occurred at low ${\textit{Wi}}$ , i.e. $Wi \leqslant 2$ . The flow remained laminar despite the presence of bubbles, confirming that the initial discontinuity (bubble injection) alone is insufficient to trigger transition to turbulence. Transition is finally observed at ${\textit{Wi}}=5$ . The transition is governed by the interplay of elasticity and streamline curvature but not the transient discontinuity. The destabilisation coincided with the development of strong streamline curvature around bubbles, consistent with prior studies showing that viscoelastic stresses amplify curvature-driven instabilities, as first observed and quantified by McKinley, Pakdel & Öztekin (Reference McKinley, Pakdel and Öztekin1996) and has also been confirmed in our earlier work as well (Naseer et al. Reference Naseer, Izbassarov, Ahmed and Muradoglu2024).

Figure 7(a) shows the distribution of first normal stress difference ( $N_1$ ) in the mid-plane for three different values of ${Re}=10$ , ${Re}=100$ and ${Re}=1000$ at ${\textit{Wi}}=10$ . A ‘chaotic arrowhead’ regime is observed for all three values of ${Re}$ at this relatively low ${\textit{Wi}}$ for which the single-phase flow remains fully laminar. Note that all bubbles are of the same size in figure 7(a) and the seemingly different appearance of the bubble sizes is simply caused by the out-of-plane bubbles at this particular time instant. The constant contours of $Q$ -criterion at the walls for the same three cases are depicted in figure 7(b), showing a chaotic pattern similar to that observed in the single-phase cases at a low value of ${\textit{Wi}}$ (figure 3 d). The contours of $Q$ -criterion at the wall are shown in figure 7(c) for ${Re}=10$ , $100$ and $1000$ at ${\textit{Wi}}=1$ . It is observed that, for this low ${\textit{Wi}}=1$ case, the flow remains laminar at the walls for the ${Re}=10$ and ${Re}=100$ cases. Small perturbations remain confined to the immediate surroundings of bubbles. The positive and negative structures of $Q$ -criterion are only visible at the walls once plotted at negligibly smaller values (i.e. $\pm 1\times 10^{-5}$ ). For ${Re}=1000$ , however, a chaotic flow regime is observed even at ${\textit{Wi}}=1$ as seen by the contours of $Q$ -criterion (figure 7 c).

Figure 7. (a) Contours of $N_1$ in the mid-plane for different values of ${Re}$ at ${\textit{Wi}}=10$ . Iso-surfaces of $Q$ - $\textrm{criterion}$ at the bottom walls of the channel for the multiphase EIT regime at (b) ${\textit{Wi}}=10$ and at (c) ${\textit{Wi}}=1$ for different values of ${Re}$ . For panel (b), the contours are plotted at $\pm {0.00001}$ , $\pm {0.001}$ and $\pm {0.01}$ for the ${Re}=10, 100$ and $1000$ cases with red colour (positive) and grey (negative), respectively. For panel (c), the contours are plotted at $\pm {0.00001}$ for ${Re}=10, 100$ and at $\pm {0.1}$ for the ${Re}=1000$ case. (d) Distribution of bubbles coloured by the magnitude of streamwise flow velocity in the channel for different values of ${Re}$ at ${\textit{Wi}}=1$ .

The bubble distributions in the duct are depicted in figure 7(d). It is interesting to see that, although the bubble distribution is essentially uniform in the other cases, the majority of the bubbles are aligned in the centre of the duct forming a string-shaped pattern for the lower values of ${Re}=10$ and ${\textit{Wi}}=1$ , which also explains the reason why the disturbances remain confined to the core region of the duct away from the walls for this case. In the case of non-colloidal spherical solid particles, Won & Kim (Reference Won and Kim2004) have experimentally shown that although the main driving force for the lateral migration of particles is the first normal stress difference, particle alignment is actually promoted by the shear-thinning effect for a particular rheological setting. In the present case of bubbles, the role of shear-thinning effect in the alignment of bubbles is investigated by performing another simulation at ${Re}=10$ and ${\textit{Wi}}=1$ with a higher shear thinning and the results are shown in figure 8. The change in the effective viscosity ( $\mu _e$ ) of the flow due to the shear-thinning effect at different shear rates ( $\dot {\gamma }$ ) is governed by the concentration of polymers ( $\beta$ ) and the mobility factor ( $\alpha$ ) in the Giesekus model. The slope of viscosity versus shear rate plot is controlled by the mobility factor $\alpha$ and the final value of viscosity at a high enough shear rate is determined by the concentration of polymers $\beta$ (Yoo & Choi Reference Yoo and Choi1989). The vertical line in figure 8(a) indicates the shear rate at the wall for the mid-plane of the duct for the ${Re}=10$ and ${\textit{Wi}}=1$ case. Once the mobility factor is increased to $\alpha =0.1$ , a random distribution of bubbles is observed instead of the string-shaped aggregation at the centre of the duct (figure 8 b). We thus conclude that, unlike the solid particles, a higher shear-thinning effect of the ambient fluid breaks up the alignment of bubbles at least for this particular rheological setting. Further investigation is required to reveal the exact mechanism behind the role of the shear-thinning effect on different forces acting on the bubbles and causing them to align in the central region of the duct.

Figure 8. (a) Change in the effective viscosity of a Giesekus fluid in a concentrated polymer solution. The vertical line shows the value of shear rate at the wall for the central $yz$ -plane of the present duct flow. (b) Distribution of bubbles in the duct for different values of shear-thinning parameter $\alpha$ ( ${Re}=10, {\textit{Wi}}=1, Ca=0.01, \beta =0.1$ ).

Figure 9(a) shows that, for the case of $({\textit{Re}},{\textit{Wi}})=(10,1)$ , the majority of bubbles accumulate in the core region, and the perturbations remain confined to immediate surroundings of the bubbles indicating a laminar flow, as shown by the contours of $Q$ -criterion. When the inertia is increased ( ${Re}=100$ ) at the same low value of ${\textit{Wi}}=1$ , the bubbles become more uniformly distributed, but the flow still remains laminar. Finally, once the inertia is sufficiently high ( ${Re}=1000$ ), the bubbles become randomly distributed in the entire duct and a fully chaotic flow regime is observed in the entire domain (figure 9 a) with a significant increase in the friction drag. Figure 9(b) shows a 1-D energy spectrum of the streamwise component of flow velocity for ${\textit{Wi}}=1$ , ${\textit{Wi}}=5$ and ${\textit{Wi}}=10$ at ${Re}=10$ . Once the flow reaches a statistically steady state, the velocity signal is collected from the point away from the centre ( $x^*=0.5, y^*=6, z^*=1$ ) to avoid the noise due to the bubble cluster. The flow velocity inside the bubbles is filtered out from the signal. A slope of $-2$ is observed for the ${\textit{Wi}}=5$ and $10$ cases when the flow is chaotic in the entire domain. This slope is the same as the temporal spectrum scaling of centreline velocity observed by Lellep, Linkmann & Morozov (Reference Lellep, Linkmann and Morozov2024) indicating a similar chaotic state. The energy spectra of the ${\textit{Wi}}=5$ and ${\textit{Wi}}=10$ cases overlap each other, showing a similar flow state for both the cases in the presence of bubbles. As the flow is not fully chaotic at ${\textit{Wi}}=1$ and its signal shows a very high intermittency, the same is manifested in its spectrum where the slope does not match with the slopes of the ${\textit{Wi}}=5$ and $10$ cases. For the ${Re}=100$ cases, the slope for ${\textit{Wi}}=1$ gets closer to the slopes in the ${\textit{Wi}}=5$ and $10$ cases, and finally at ${Re}=1000$ , all the three cases show a similar slope of $-2$ (figure 9 d). This slope of $-2$ is greater than the classical slope of $-5/3$ in the inertial turbulence, but still less than $-4$ reported for a purely elastic turbulence in the absence of inertia (Datta et al. Reference Datta2022) or at a very high value of ${\textit{Wi}}=1000$ observed in the single-phase cases (figure 4 a, c). As this multiphase EIT regime is achieved at the lower values of ${\textit{Wi}}$ , it is governed by both inertia and viscoelasticity, and hence the slope remains in between these two limits.

Figure 9. (a) Iso-surfaces of $Q$ -criterion are shown at different values of ${Re}$ for ${\textit{Wi}}=1$ . The contours are plotted at $\pm 0.002$ , $\pm 0.2$ and $\pm 2$ for the ${Re}=10, 100$ and $1000$ cases, respectively. The bubbles are marked with green colour, while the positive and negative values of $Q$ -criterion are marked with red and grey colours, respectively. (b, c, d) One-dimensional energy spectrum, (e, f, g) the distribution of various components of shear stress and (h, i, j) turbulent kinetic energy ( $\mathcal {K}$ ) are plotted in the mid-plane for the ${Re}=10$ (left), ${Re}=100$ (middle) and ${Re}=1000$ (right) cases, respectively ( $\beta =0.1$ and $Ca=0.01$ ).

The distribution of different components of shear stress from the channel wall towards the centre are shown in figures 9(e), 9(f) and 9(g) in the mid-plane of the duct for ${Re}=10, 100$ and $1000$ cases, respectively. As the flow is multiphase, an additional contribution from the surface tension force of the bubbles ( $\bar \tau _{{S}}$ ) is added towards the total shear stress balance. However, this contribution remains much smaller than the remaining components as the bubble volume fraction is a mere $3\,\%$ . Similar to the single-phase flow cases, these stress components are normalised by the local shear stress at the wall of the mid-plane. As seen in the figure, the chaotic flow regime is dominated by the viscoelastic stress ( $\bar \tau _{{P}}$ ) followed by the viscous stress due to the mean flow ( $\bar \tau _{{N}}$ ). These two stress components are highest at the wall and decay to zero towards the channel centre where the shear rate is minimum. The contribution of the Reynolds stress ( $\bar \tau _{{R}}$ ) remains negligible for the ${Re}=10$ and ${Re}=100$ cases. The Reynolds stress component $\bar \tau _{{R}}$ starts to show its effect once the inertia becomes high enough in the ${Re}=1000$ cases. As the source of flow instability is the curved streamlines around the bubble in this multiphase flow regime (McKinley et al. Reference McKinley, Pakdel and Öztekin1996), the transition to the EIT state greatly depends upon the distribution of bubbles across the channel. The role of viscoelasticity in promoting the bubble migration towards the channel centre or towards the wall is governed by the relative change in the convective time scale of the bubbles and the polymer relaxation time, as has been explained by Bothe et al. (Reference Bothe, Niethammer, Pilz and Brenn2022). A similar observation was also made in our earlier work (Naseer et al. Reference Naseer, Izbassarov, Ahmed and Muradoglu2024). The distribution of turbulent kinetic energy (TKE) for all these three values of ${Re}$ shows a peak value near the wall and minimum at the channel centre, as shown in figures 9(h), 9(i) and 9(j). By increasing the value of ${\textit{Wi}}$ , the magnitude of TKE also increases following the same qualitative trend in agreement with our previous study of viscoelastic turbulent single phase case (Izbassarov et al. Reference Izbassarov, Rosti, Brandt and Tammisola2021b ).

Next, further simulations are performed to examine the effects of size and number of bubbles on the transition and resulting chaotic flow for the $({\textit{Re}},{\textit{Wi}}) =(1000,1)$ case. First, the number of bubbles are increased from the previous volume fraction of $\varPhi =3\,\%$ to $\varPhi =9\,\%$ while keeping the bubble size the same, and then the bubble volume is decreased by $50\,\%$ while keeping the volume fraction the same at $\varPhi =3\,\%$ . As shown in figure 10, neither the number of bubbles nor their size has any significant effect on the scaling of energy spectrum except for a slight increase in the energy content at the smaller scales for $\varPhi =9\,\%$ . Moreover, no significant change in transition to a chaotic state is observed by varying the number of bubbles or their size.

Figure 10. Effects of number of bubbles and bubble size. The bubble distribution (left) and the evolution of flow velocity (right) for the void fraction of (a) $\varPhi =3\,\%$ , (b) $\varPhi =9\,\%$ for the same bubble size of $d_b/2h=0.2$ and (c) for the bubble size of $d_b/2h=0.2/\sqrt [3]{2}=0.159$ with $\varPhi =3\,\%$ . (d) One-dimensional energy spectra for the same cases ( ${Re} = 1000, Wi = 1, Ca = 0.01, \alpha = 0.001$ ).

Once the flow reaches a steady state, the applied pressure gradient is balanced by the total shear stress at the walls and the average wall shear stress can be simply computed as $\overline \tau _w = -({h}/{2})({{\rm d}p_o}/{{\rm d}y})$ . The change in the drag is quantified by $\Delta D=(\overline \tau _{w} - \tau _{w_o})/\tau _{w_o}$ , where $\tau _{w_o}$ is the wall shear stress of the corresponding Newtonian laminar flow for the same Reynolds number. Figure 11 summarises the percentage change in drag for all the multiphase simulations. The single-phase viscoelastic laminar flow has lower drag as compared with the Newtonian laminar flow for the same value of ${Re}$ . Understandably, as the inertia increases in this shear thinning viscoelastic fluid once the flow is made viscoelastic, the drag is reduced. For a particular value of ${\textit{Wi}}$ , the percentage reduction in drag is found to be the same for every value of ${Re}$ . However, the drag increases significantly once a viscoelastic flow transitions to a chaotic turbulent state after injecting the bubbles (figure 11). In particular, a rapid increase in the drag is observed at ${Re}=1000$ once the flow becomes turbulent for all three values of ${\textit{Wi}}$ . It is interesting to observe that, although the drag increases for the ${Re}=10$ and $100$ cases at ${\textit{Wi}}=5$ and $10$ upon transition to the EIT state, it still remains less than that of the respective Newtonian laminar state. Another interesting feature is also observed at ${\textit{Wi}}=1$ for the ${Re}=10$ and $100$ cases. Once bubbles are injected into the flow, minor fluctuations remain confined to the core region away from the walls where the bubbles collect and form a string-shaped pattern ( ${Re}=10$ ) or small clusters ( ${Re}=100$ ). The flow remains essentially laminar closer to walls and in the major portion of the duct. In this regime, the drag is reduced as compared with the single-phase laminar state. For this case, the bubbles forming a string-shaped pattern create a low viscosity corridor at the core region where they move faster and thus the streamwise component of flow velocity is slightly increased. As a result, the applied pressure gradient is reduced effectively to maintain the same flow rate. Furthermore, it is important to highlight that the magnitude of drag is small for this laminar flow, so the change is also small. At a higher inertia ( ${Re}=1000$ ), the flow becomes fully chaotic in the entire channel, and thus the drag is increased significantly even for the ${\textit{Wi}}=1$ case (figure 11).

Figure 11. Percentage change in drag measured by the change in total shear stress at the wall ( $\overline \tau _w$ ) as compared with the Newtonian laminar flow once the flow is made viscoelastic and is subsequently transitioned to turbulent state by injecting the bubbles.

It is interesting to note that for the $({\textit{Re}},{\textit{Wi}})=(10,10)$ and $(100,10)$ cases, the flow has become completely chaotic and the EIT state is achieved by injecting the bubbles into the flow; however, the drag is still less than the corresponding Newtonian laminar state, which is primarily attributed to the shear-thinning effect of viscoelastic fluid (figure 11).

4.2.1. TKE in a multiphase EIT regime

To get further insight into the multiphase EIT regime, all the terms in (4.1) are integrated in the streamwise direction and averaged in time after the flow reaches a statistically steady state. The results are plotted in figure 12 for ${Re}=100$ and ${Re}=1000$ at ${\textit{Wi}}=10$ . Compared with the single-phase case shown in figure 5, a prominent effect of bubbles at low Weissenberg number is observed in smoothing sharp gradients near the wall in all the terms that contribute to TKE. We note, however, that comparison between single-phase and multiphase cases should be interpreted only qualitatively since ${\textit{Wi}}$ is two orders of magnitude higher in the single-phase case. Figure 12 shows that the magnitudes of all terms are amplified as ${Re}$ is increased from 100 to 1000, but the amplification is more pronounced in advection ( $\mathcal A$ ), turbulent transport ( $\mathcal Q$ ), production by mean flow ( $\mathcal P$ ) and viscous diffusion ( $\mathcal D$ ) terms. There are also some qualitative differences in these terms in the ${Re}=100$ and ${Re}=1000$ cases. The peak values of $\mathcal A$ , $\mathcal Q$ , $P$ and $\mathcal \epsilon$ occur in the corners and in the middle of the side walls for ${Re}=100$ and ${Re}=1000$ , respectively. In contrast to the single-phase case, the maximum value of $\mathcal R$ occurs in the central portion of the duct due to the presence of bubbles (figure 12 c). Significant viscous dissipation observed in the channel centre is also attributed to the presence of bubbles there. The positive peaks of polymer work for the multiphase case shows that it still contributes towards the production instead of dissipation in some portions of the domain (figure 12 g). The body force term contributes significantly to the TKE budget now due to the work done by the surface tension (figure 12 h).

Figure 12. Contours of (a) advection by mean flow, (b) transport by velocity fluctuations, (c) transport by pressure, (d) production by mean flow, (e) viscous diffusion, (f) viscous dissipation, (g) polymer work and (h) body force terms are shown in a vertical cutting $xz$ -plane for (left) ${Re}=100, {\textit{Wi}}=10$ and (right) ${Re}=1000, {\textit{Wi}}=10$ multiphase cases, respectively.

Notably, the chaotic flow state at these lower values of Weissenberg numbers is only sustained in the presence of bubbles. If the bubbles are removed, the single-phase flow would laminarise, as observed in the simulations of single-phase flow where the chaotic flow state could not be sustained below $Wi\leqslant 40$ .

In the ET regime (very low ${Re}$ ), it is well established that a sufficiently strong extensional perturbation is required to destabilise the flow and trigger chaos. As inertia increases, the threshold for instability decreases and the flow becomes more susceptible to a wider range of perturbations – including those introduced by bubbles. In the present simulations, the presence of bubbles provides persistent, finite-amplitude disturbances that can maintain a chaotic state even as inertia becomes more significant. Thus, the bubbly chaotic state observed in the present study at moderate ${Re}$ can be interpreted as a form of ET that is transitioning towards EIT as inertia increases. This is further supported by the fact that the energy spectra do not yet reach the inertial $-5/3$ scaling, and the flow structures retain characteristic features of ET.

In summary, ET and EIT can be seen as two ends of a spectrum, with present results occupying an intermediate and transitional regime. The bubbly chaotic state in the present scenario likely represents a version of ET influenced by finite inertia and persistent multiphase perturbations, rather than a fully developed EIT state. The spectra and flow localisation support this interpretation.