3.1. Transition to the VMS regime with increasing $El$

We first investigate the gradual approach to the VMS regime as the elasticity of the fluid increases. For this initial simulation, an aspect ratio of $\varGamma =12.56$ was considered, whereas the maximum polymer extension was set to $L = 100$. A Newtonian simulation (i.e. $\beta = 1$) was initially run to calculate a TVF pattern with six pairs of counter-rotating vortices. Using this state as initial condition, the fluid's elasticity was slowly and steadily increased at a uniform rate, $El = 10^{-3}t/Re$, until a dynamical regime characterized by merging and splitting events was found. We would like to stress that the protocol followed in this simulation differs from that in the experiments, where the VMS regime is achieved by increasing $Re$ while keeping a constant elasticity (Lacassagne et al. Reference Lacassagne, Cagney, Gillissen and Balabani2020). Our study hence offers a different perspective into the pathway leading to this flow regime and allows to identify the gradual transformation the flow undergoes as the working fluid becomes more elastic.

Panel (a) in figure 1 provides an overview of the structural and dynamical changes induced by the polymers on the initial TVF as $El$ increases. It shows a space–time diagram of the radial velocity $u$ at mid-gap along the axial direction $z$, where time has been replaced by its corresponding $El$ value. Red areas represent fluid motion from the inner to the outer cylinder, i.e. outflows, whereas blue regions indicate fluid moving from the outer to the inner cylinder, i.e. inflows. Panel (b) places emphasis in the range of $El$ values for which a complex spatio-temporal dynamics happens. The stationary pattern of vortices becomes unstable at $El \approx 0.29$, leading to periodic oscillations of the vortex pairs along the $z$-axis. The onset of the VMS regime takes place soon after, at $El \approx 0.315$, as the dynamics of the distinct vortex pairs decouples and these begin to move independently in the axial direction. It will be shown later in § 3.5 that this threshold is sensitive to the number of vortices of the initial condition and the aspect ratio used in the simulations. Merging events, where two vortex pairs coalesce to form a single vortex pair, are indicated as dashed (green) rectangles in figure 1(b). These events fully dominate the dynamics in the initial phase of the VMS regime (for $0.32 < El < 0.34$) and since they occur simultaneously at different axial locations, the total number of vortex pairs in the system rapidly decreases. After this initial phase ($El > 0.34$), merging events coexist with events where a vortex pair branches into two, i.e. splitting events, shown as dash-dotted (purple) rectangles in the figure, as well as with regions where the dynamics becomes transiently chaotic (see for instance the flow region between $4 < z/d < 7$ for $0.34 < El < 0.36$ or $0 < z/d < 2.5$ for $0.38 < El < 0.40$). The number of vortex pairs fluctuates between two and four in this phase.

Figure 1. Space–time plot representing the magnitude of the radial velocity $u$ at mid-gap along the axial direction $z$. Panels (a) and (b) correspond to a simulation performed at $Re = 95$ where, starting from a Newtonian TVF, $El$ was slowly increased with time ($El = 10^{-3}t/Re$). Note that the time $t$ has been replaced by the corresponding $El$ values on the horizontal axis. Panel (a) shows the variation of $u$ from Newtonian flow (i.e. $El = 0$) up to the largest $El$ value simulated ($El = 0.50$), whereas panel (b) shows in more detail the range of $El$ values for which a complex spatio-temporal dynamics take place. Panel (c) illustrates the spatio-temporal dynamics when $El$ is kept constant after the VMS regime is achieved. The case exemplified corresponds to $El = 0.32$. Red and blue areas indicate outflows and inflows respectively. Note that periodic boundary conditions are used in $z$.

It is important to note that the occurrence of VMS events does not depend on the continuous increase of $El$ with time. If $El$ is held constant after the VMS regime is achieved, the simulations show the same dynamic events just described: vortex merging, vortex splitting and transient chaotic motion, reflecting that these are temporal characteristics of the flow that occur when $El$ exceeds a certain critical threshold. This is demonstrated in figure 1(c), which shows a space–time map for a simulation where $El$ has been fixed to $0.32$. Interestingly, the VMS events observed in simulations with constant $El$ are similar to those observed in simulations where $El$ varies with time. The reason (which will be discussed later in the paper) is that increasing $El$ has little influence on the vortices in this flow regime. As a result, space–time diagrams corresponding to simulations where $El$ changes with time not only feature the various flow regimes obtained when $El$ is varied but also provide an accurate representation of the VMS events.

Another important feature that is clearly illustrated in figure 1(a) is the strong impact that increasing $El$ has on the structure of the vortex pairs. We anticipate here that these structural changes are key to understanding the physics underlying merging and splitting events. Hence, before getting into detail about the dynamics in the VMS regime, it is convenient to present a comprehensive study about the influence of elasticity on the TVF.

3.2. Viscoelastic modification of the TVF

A well-known property of viscoelastic flows with curved streamlines is the appearance of a radially inward force which is caused by the elastic stresses arising from the stretching of the polymer molecules by the primary flow (Groisman & Steinberg Reference Groisman and Steinberg1998). This force has been identified as the driving source of a number of instabilities in curvilinear flows of highly elastic polymer solutions, which are usually known as purely elastic instabilities (Shaqfeh Reference Shaqfeh1996). The mechanism underlying these instabilities has been discussed in detail and verified in many studies, particularly in flow regimes where inertial effects are negligible ($Re \to 0$). It is, however, reasonable to expect that this elastic force will also have an influence in parameter regimes where Newtonian flows become unstable due to inertial forces. The stationary pattern of Newtonian Taylor vortices used as starting solution in our calculations is one such case: it arises from a centrifugal instability of the purely azimuthal primary flow (Taylor Reference Taylor1923). This instability mechanism is expected to persist in the viscoelastic case as $El$ is slowly increased starting from the Newtonian limit. However, the structure of the Taylor vortex pattern is likely to be modified by the competition between the centrifugal and elastically induced forces. Additionally, if the fluid's elasticity becomes sufficiently large, the elastic instability mechanism might replace the centrifugal mechanism, leading to a flow state that is elastic in nature but whose structure could be modified by the presence of inertial effects. In this section we show that this is indeed the case in our simulations.

Figure 2 shows colour maps of the radial velocity, $u$, illustrating the dependence of the flow structure as $El$ increases from the Newtonian limit ($El = 0$) up to the regime in which the flow exhibits spatio-temporal behaviour. Note that, to save space, in all figures illustrating flow patterns throughout the paper, the system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction, so that the inner (outer) cylinder is located at the bottom (top) of each panel and the positive $z$-direction goes from right to left (see the coordinate system in the bottom panel of figure 2). The structure of the TVF pattern in the Newtonian case (top panel) shows a small asymmetry between outflows and inflows, as the axial extent of the inflows is slightly greater than that of the outflows. This characteristic fully reverses as elasticity comes into play. The axial extent of the inflows decreases with increasing $El$ and these become eventually confined to strong jets that extend over narrow regions in the axial direction. Conversely, the axial extent of the outflows increases notably (note that they become nearly four times larger than the inflows for $El \geq 0.12$) and the magnitude of $u$ in these regions decreases substantially as $El$ increases (the range of values of $u$ corresponding to each panel is specified in the caption).

Figure 2. Structural variation of the TVF pattern as $El$ increases from the Newtonian case (i.e. $El = 0.00$) up to values near the onset of spatio-temporal dynamics. Each panel shows a colour map of the radial velocity $u$ in a meridional plane $(z/d,r/d) \in [0,4{\rm \pi} ] \times [3.35,4.35]$, where red regions indicate outflows (i.e. positive velocity), blue regions represent inflows (i.e. negative velocity) and the zero velocity has been set as grey. The colour scale is based on the maximum and minimum values of each case and hence differs among the different panels. There are four positive and four negative contours evenly distributed across the full range of values in each case: 1. $El = 0.00$, $u$ in $[-0.037,0.059]$; 2. $El = 0.03$, $u$ in $[-0.043,0.047]$; 3. $El = 0.05$, $u$ in $[-0.048,0.029]$; 4. $El = 0.10$, $u$ in $[-0.048,0.014]$; 5. $El = 0.15$, $u$ in $[-0.082,0.015]$; 6. $El = 0.25$, $u$ in $[-0.072,0.011]$. The system is shown rotated by $90^{\circ}$ in the counterclockwise direction with respect to its original position. The location of the inner and outer cylinders, $r_i$ and $r_o$, respectively, as well as the locations of the top ($z/d = 4{\rm \pi}$) and bottom ($z/d =0$) of the system are indicated in the bottom panel.

It was postulated in a previous experimental study that this strong asymmetry between inflows and outflows might be caused by the work done by the elastically induced force (Groisman & Steinberg Reference Groisman and Steinberg1998). To verify this hypothesis quantitatively, figure 3 shows axial profiles of the elastic force, hereafter denoted as $F_e$ (panels a,c,e), and its associated work $F_e u$ (panels b,d,f), obtained at the mid-gap for the last three cases shown in the figure 2. As seen, the profiles of $F_e$ are always negative, reflecting that $F_e$ is an inward force, and they exhibit strong peaks in the inflows whose magnitude increases with increasing $El$. Since $F_e$ acts in the same direction as $u$ in the inflows, it does positive work on the flow in these regions. This circumstance implies that the strong peaks of $F_e$ will result in large positive work (see the peaks of $F_e u$ in panels b,d,f), which enhances the fluid motion in the inflows and creates the strong localized jets that appear as $El$ increases. In the outflows, on the contrary, $F_e$ acts in opposition to the fluid motion and therefore does a negative work on the flow. This characteristic explains the decay in the magnitude of $u$ that is observed in the outflows as $El$ increases. The axial extent of the inward jets decreases as the magnitude of the peaks grows, which evidently entails an increase in the axial extent of the outflows and creates the asymmetry between inflows and outflows observed in figure 2.

Figure 3. Axial profiles of the elastic force $F_e$ (left panels) and its associated work $F_e u$ (right panels) obtained at the mid-gap for $El$ values matching those of the last three panels in figure 2. The elastic force is calculated as $F_e = ({(1-\beta )}/{Re}) (\partial _r T_{rr} + {(T_{rr}-T_{\theta \theta })}/{r} + \partial _z T_{rz})$. A dashed line has been added at $F_e u = 0$ to help identify inflows ($F_e u > 0$) and outflows ($F_e u < 0$).

In addition to the emergence of this asymmetry, a second transformation takes place inside the outflows. The region where the largest positive velocity occurs, which in the Newtonian case is located at the centre of the outflows, separates in the viscoelastic cases into two identical regions which are symmetric with respect to the central symmetry plane of the outflow. These new regions of maximum positive velocity move away from each other as $El$ increases and approach progressively the adjacent inflows. When the elasticity is sufficiently large ($El \geq 0.12$), the strong inflow jets are flanked by these regions of maximum positive velocity, whereas the flow in the central part of the outflows becomes nearly uniform in the axial direction.

A remarkable feature of this structural transition is the fact that the changes are most pronounced in the low $El$ regime ($El < 0.12$). This observation is quantitatively confirmed by the axial profiles of $u$ at the mid-gap shown in figure 4. Profiles corresponding to $El < 0.12$ (shown in panel (a)) differ markedly and clearly reflect strong changes in both the magnitude of $u$ (particularly in the outflows) and the axial extent of inflows and outflows. However, for $El \geq 0.12$ (see panel (b)), the differences among profiles are small and mainly occur in the magnitude of $u$, which keeps slightly decreasing (increasing) in the outflows (inflows) with increasing $El$. Further quantitative evidence of this behaviour is given in figure 5. Panel (a) in this figure shows the dependence of the axial extent of the inflow and outflow regions at the midgap with increasing $El$. It is apparent that the largest variation in the extent of these regions (which are naturally inversely proportional) occur within the low $El$ regime, for $0.05 < El < 0.1$. Likewise, the sharpest change in the ratio between the maximal velocity of outflows and inflows (shown in panel (b)) also happens at very low $El$ values ($El < 0.05$), where the strength of the outflows decays strongly. These observations clearly evidence that $F_e$ exerts a surprisingly strong influence on the flow structure even in weakly elastic fluids. Another interesting feature revealed by the two panels of figure 5 is that the flow characteristics appear to exhibit asymptotic behaviour. As $El$ increases, the sizes of the outflows and inflows approach values close to $L_{out}/d \approx 1.7$ and $L_{in}/d \approx 0.38$, respectively, whereas the ratio between the maximal velocities in outflows and inflows appears to level off at approximately $0.13$. This observation is also corroborated by the velocity profiles, which seem to be gradually converging with increasing $El$ (see panel (b) in figure 4).

Figure 4. Profiles of the radial velocity $u$ at mid-gap along the $z$-axis. Panel (a) shows profiles for $El < 0.12$, where profound changes in the structure of the flow pattern are observed in figure 2, whereas panel (b) focuses on the range of $El$ values ($El \geq 0.12$) where the variation in the profiles is small. Note that in both panels $u$ is normalized with the largest absolute value among all cases which is obtained for $El = 0.12$ and corresponds to $|u_{max}| = 0.0824$.

Figure 5. Quantification of the changes in the structure and strength of the vortices as $El$ increases. Panel (a) shows the variation in the axial extent of inflows and outflows at the midgap in units of the gap-width, $d$, as the fluid becomes more elastic. Panel (b) displays the ratio between the maximal values of $u$ in the outflows and inflows as a function of $El$. The dashed line indicates the approximate value this ratio seems to approach asymptotically as $El$ increases.

An important distinction between the low and high $El$ regimes (i.e. $El < 0.12$ and $El \geq 0.12$) is the magnitude of $u$ in the inflows. As seen in figure 4(a), the maximum velocity of the inflows in the low $El$ regime increases initially with increasing $El$, but eventually converges to a value close to $u/|u_{max}| = -0.55$. This value is substantially lower than those shown by the profiles in the high $El$ regime (see figure 4b), where $u/|u_{max}|$ ranges from $-1$ at $El = 0.125$ (when it is maximal) to $\sim -0.78$ at $El > 0.25$ (when the profiles seem to converge). The transition between both regimes can be clearly identified in the space–time plot of figure 1(a) as a sudden change in the colour intensity that takes place at $El \sim 0.12$. The abrupt nature of this transition strongly suggests that it may be caused by a change in the physical mechanism associated with the instability of the primary flow. To test this hypothesis we examine the integral energy budgets. For viscoelastic flows, the energy balance reads (Dallas, Vassilicos & Hewitt Reference Dallas, Vassilicos and Hewitt2010; Dubief, Terrapon & Soria Reference Dubief, Terrapon and Soria2013),

(3.1)\begin{equation} \int_V \mathcal{P} \, {\rm d} V - \int_V \epsilon \, {\rm d} V - \int_V \varPi_e \, {\rm d} V = 0, \end{equation}

where $\mathcal {P}$ is the kinetic energy production, $\epsilon$ is the viscous dissipation rate and $\varPi _e$ denotes the work done by the elastic stresses. These quantities were calculated using the following expressions:

(3.2)$$\begin{gather} \mathcal{P} ={-}\overline{u'v'} \frac{\partial \bar{v}}{\partial r} + \overline{u'v'} \frac{\bar{v}}{r}, \end{gather}$$

(3.3)$$\begin{gather}\epsilon = \frac{2\beta}{Re} \overline{\boldsymbol{\mathsf{S}}':\boldsymbol{\mathsf{S}}'}, \end{gather}$$

(3.4)$$\begin{gather}\varPi_e = \frac{1-\beta}{Re} \overline{\boldsymbol{\mathsf{S}}':\boldsymbol{\mathsf{T}}'}. \end{gather}$$

Here, the overline denotes axially averaged quantities, $\boldsymbol{\mathsf{S}}' = (\boldsymbol {\nabla } \boldsymbol {v}' + \boldsymbol {\nabla } \boldsymbol {v}'^{\rm T})/2$ is the rate of strain tensor and the prime symbol indicates deviations of the velocity or polymer stress tensor from their axially averaged values ($\bar {\boldsymbol {v}}=(0,\bar {v},0)$ and $\bar {\boldsymbol{\mathsf{T}}}$). It must be clarified that, although this equation was derived in the context of turbulent flow, it also applies to steady and axisymmetric vortex flow (the derivation of the equation for this particular case is given in Appendix A). The first and second integrals in (3.1) are always positive, meaning that they act as source and sink terms of the energy balance, respectively (note that there is minus sign in front of the second integral). The sign of the third integral can be positive or negative. If it is positive, this term has a negative contribution to the balance and thus polymers act to dissipate the fluid's kinetic energy. By contrast, if it is negative, polymers act as an energy source. The variation of the values yielded by these integrals with increasing $El$ is shown in figure 6. As expected, the behaviour of the polymers changes drastically at $El \sim 0.12$, consistent with the transition between the low and high elasticity regimes. At low $El$ values, polymers play a dissipative role, helping the viscous forces to damp the centrifugally induced vortices. However, at $El \sim 0.1$, the work done by the elastic stresses changes sign and the net contribution of the polymers to the energy balance becomes positive, indicating that they inject energy into the flow through the elastic stresses. The amount of energy that the polymers supply to the system is initially very small (for $0.1 \leq El < 0.12$) but increases suddenly when $El \sim 0.12$. After this transition occurs, the energetic contribution of the polymers becomes the dominant energy source and its magnitude continues increasing with increasing $El$. In contrast, the energy production due to inertial mechanisms remains small and decreases very gradually as $El$ increases. From this analysis, it is clear that the nature of the mechanism driving the instability indeed changes from being centrifugal ($El < 0.12$) to being elastic ($El \geq 0.12$), a characteristic that sets a clear distinction between the two regimes investigated so far.

Figure 6. Variation of the integral energy budgets with increasing $El$ before the spatio-temporal dynamics sets in. Here, $\mathcal {P}$ denotes the kinetic energy production by the inertial forces, $\epsilon$ stands for the viscous dissipation and $\varPi _e$ is the work done by the elastic stresses. The dotted (orange) vertical line indicates the $El$ value at which the transition between the low and high $El$ regimes takes place. Hereafter, these regimes are denoted as centrifugally and elasticity dominated regimes, respectively.

Finally, to facilitate comparison between the flow structure in the high $El$ regime and other elastically induced stationary patterns previously reported, the streamlines of the flow $\psi$ at $El = 0.25$ are shown in the bottom panel of figure 7. These are naturally very different from those in the Newtonian case (also shown for comparison in the top panel) and reflect again the structural changes just discussed. As seen, unlike the Newtonian case, where the vortices are nearly equidistant, in the elastically dominated regime they appear to be arranged in pairs, with their cores being located very close to one another. We note that this type of structure has been previously reported in the literature and it is usually known as DW (Groisman & Steinberg Reference Groisman and Steinberg1997; Lange & Eckhardt Reference Lange and Eckhardt2001; Thomas et al. Reference Thomas, Sureshkumar and Khomami2006, Reference Thomas, Khomami and Sureshkumar2009), due to its similarity to the shape of a magnetic dipole. However, there are a couple of important differences between the structures described in previous works and the one presented here. A characteristic shared by all previous studies is that DW appear after a hysteretic transition, when $Re$ is decreased starting from a flow state driven by an elastic instability. In fact, it is often stated in the literature that flow deceleration is a necessary condition to observe these structures (Groisman & Steinberg Reference Groisman and Steinberg1997; Lange & Eckhardt Reference Lange and Eckhardt2001; Thomas et al. Reference Thomas, Sureshkumar and Khomami2006; Lacassagne et al. Reference Lacassagne, Cagney, Gillissen and Balabani2020). The present study shows that this is not the case and that, at least in the regime investigated here, these structures may also appear for a fixed $Re$ if the elasticity of the working fluid is sufficiently large so that the elastic instability mechanism replaces the centrifugal mechanism. In the low $Re$ regimes where most previous studies were conducted, DW appear to be localized in the axial direction, i.e. there are regions where the flow is laminar interspersed between distinct DW. When the distance between DW becomes less than $5d$, they approach each other and coalesce. This characteristic is so far absent in the present simulation. Despite the distance between DW being substantially lower than $5d$, they remain stationary and form a pattern of equally spaced structures along the axial direction. A possible reason for such difference will be discussed later in § 3.6. We would like to finally note that similar arrangements of DW, yet not stationary, have also been reported in the literature, which were dubbed OS (Groisman & Steinberg Reference Groisman and Steinberg1996; Thomas et al. Reference Thomas, Sureshkumar and Khomami2006, Reference Thomas, Khomami and Sureshkumar2009).

Figure 7. Comparison between the flow streamlines $\psi$ in the centrifugally (top) and elastically (bottom) dominated regimes. Positive (negative) values are represented as red (blue), whereas the zero value is shown as grey. The colour scale is based on the maximum and minimum values of each case and hence differs between the panels. There are four positive and four negative contours evenly distributed across the full range of values in each case: 1. $El = 0.00$, $\psi$ in $[-0.660,0.660]$; 2. $El = 0.25$, $\psi$ in $[-0.295,0.295]$. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction.

3.3. Onset of spatio-temporal dynamics

The stationary pattern of DW loses its stability at a Hopf bifurcation which takes place at $El \sim 0.29$ leading to an axial oscillation of the vortices. This is illustrated in figure 8 through colour maps of $u$ taken at five equally spaced time instants within a period (shown as circles in figure 9a). The displacement of the vortices is relatively small, yet clearly discernible in the figure by looking at the axial position of the inflows. It should be recalled that the system is shown rotated by $90\,^{\circ}$ counterclockwise, and so the upwards (downwards) motion of the inflows corresponds to leftwards (rightwards) motion in these figures. Starting from the state where the inflows are at their lowest axial positions (panel ${\rm A}$), it is observed that the inflows move first axially upwards (panel ${\rm B}$), reach the position of maximum displacement (panel ${\rm C}$) and subsequently move downwards (panel ${\rm D}$), returning eventually to the initial state (panel ${\rm E}$, which is identical to ${\rm A}$). A notable difference with respect to the stationary case is the breakdown of the symmetry in the outflows. When the vortices move axially upwards, the lower half of the outflow (see region enclosed by the (green) dashed rectangle in the panel ${\rm B}$) remains similar to that in the stationary case (see bottom panel in figure 2), however, the upper half (marked by the purple dash-dotted rectangle in panel ${\rm B}$) notably changes due to an increase in the maximal velocity next to the inflow. The opposite is observed when the vortices move downwards. The upper half (shown as a green dashed rectangle in panel ${\rm D}$) of the outflow remains as in the stationary case, whereas the lower half (purple dash-dotted rectangle in panel ${\rm D}$) takes a similar form to that of the upper half during the upward motion. As a consequence, flow states moving axially upwards and downwards, where the vortices are located at the same axial positions, exhibit antisymmetric outflows (that is the case, for example, for states ${\rm B}$ and ${\rm D}$ shown in the figure).

Figure 8. Colour maps of $u$ illustrating the axial oscillation of the vortex pairs that emerges from the Hopf bifurcation. Letters from A to E are used to show the correspondence between each panel and the time series shown in figure 9. Positive (negative) velocity is represented as red (blue), whereas the zero velocity is shown as grey. There are four positive and four negative contours evenly distributed in $u \in [-0.0780,0.0095]$. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction. Note that the flow patterns shown in ${\rm A}$ and ${\rm E}$ are identical, whereas those shown in B and D only differ in that their outflows are antisymmetric.

Figure 9. Periodic motion arising from the instability of the stationary pattern of DW. (a) Time series of the axial velocity $w$ for $El = 0.3$ calculated at $r/d = 4.15$ and $z/d = {\rm \pi}$ . The circles indicate time instants for which colour maps of $u$ are presented in figure 8. (b) Power spectral density (PSD). The frequency is normalized with the elastic frequency, $f_e$.

The frequency of the oscillation was determined by applying fast Fourier transform to a time series of the axial velocity $w$ obtained at a radial location close to the outer cylinder (figure 9a). The power spectral density is shown in figure 9(b), where the frequency is normalized with the elastic frequency, $f_e$. Following Lacassagne et al. (Reference Lacassagne, Cagney, Gillissen and Balabani2020), $f_e$ was calculated as $f_e = 2c_e/k_{avg}$, where $c_e$ denotes the wave celerity, $c_e = \sqrt {\nu /\lambda }$, and $k_{avg}$ is the average spatial wavelength of the vortex flow pattern. The spectrum shows a pronounced peak at $f_e/3$, which clearly indicates that this is the dominant frequency of the oscillation. Other peaks with smaller amplitudes are also observed. However, they correspond in all cases to other sub-harmonics of $f_e$ and are therefore commensurate with the dominant frequency. It should be noted that $f_e/3$ is the dominant frequency for the particular case where the flow pattern has ${\rm six}$ vortex pairs. For other flow patterns with different numbers of vortex pairs, the oscillation is characterized by other subharmonics of $f_e$. The fact that frequencies appear in the spectrum as subharmonics of $f_e$ is in full agreement with the experimental observations (Lacassagne et al. Reference Lacassagne, Cagney, Gillissen and Balabani2020) and reflects once again the elastic nature of the instabilities taking place at these $El$ values.

From a mechanistic perspective, the periodic up and down motion of the vortices is just a consequence of the physics described in the previous section. The distance between the centres of the vortices on either side of the inflows keeps decreasing (albeit very gradually) as $El$ increases, leaving a gap between DW where the radial velocity is increasingly weak. The axial velocity, whose role before the instability onset is simply to transport the fluid vertically near the cylinders (see top panel in figure 10), eventually penetrates into these intermediate regions, connecting adjacent vortices to each other in the outflow region and giving rise to an axial wave. This wave propagates first axially upwards (middle panel in figure 10) and subsequently reflects back and travels axially downwards (bottom panel in figure 10), thereby creating a standing wave. The interaction between standing wave and vortices lead to the axial oscillation of the latter illustrated above.

Figure 10. Colour maps of $w$ illustrating the standing wave causing the axial oscillation of the flow pattern. The top panel shows $w$ for a stationary state at $El = 0.25$, whereas the middle and bottom panels correspond to the states B and D indicated in figure 9. Positive (negative) values of $w$ are shown as red (blue), whereas that for $w=0$ is shown as grey. From top to bottom, there are ${\rm four}$ positive and ${\rm four}$ negative contours evenly distributed in $w \in [-0.030,0.030]$, $[-0.026,0.030]$ and $[-0.030,0.026]$, respectively. Two additional contours at $w = \pm 0.003$ are also added to help visualize the axial waves. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction.

3.4. VMS dynamics

With a further increase in $El$, the dynamics of the different DW decouples and these begin to travel independently along the axial direction. A merging event happens when two DW move towards each other and eventually coalesce into a single entity. This phenomenon is illustrated in figure 11, which shows the variation of the radial velocity over the initial stage of the VMS regime in the simulation performed at constant $El=0.32$. Two merging events occur simultaneously in the time window shown (corresponding to the green dashed rectangle in figure 1c). In the initial phase of the merging process (up to $t \approx 500$), the second and fifth DW starting from the top (indicated by black leftwards arrows) leave their axial locations and begin to move upwards, i.e. leftwards in the rotated figure. This motion becomes apparent by comparing the positions of the inflows between the first ($t=0$) and second ($t=475$) panels. The inflows of the second and fifth DW have notably moved towards the top of the system by $t = 475$ and differ from the others in that they are tilted slightly upwards (a distinctive feature of the DW which are moving upwards). The inflows of the other DW on the other hand remain at their axial positions and only exhibit small oscillations caused by the instability discussed in the section above. As the cores of the DW travelling upwards approach the cores of the adjacent DW, they experience an attractive interaction which results in the first and fourth DW (indicated by red rightwards arrows) travelling axially downwards, i.e. rightwards in the rotated figure. Note that, as opposed to the DW moving upwards, the inflows of the DW moving downwards are tilted downwards. The attractive interaction between DW becomes stronger as their cores get close to each other, leading to a rapid increase of their travelling speeds. This characteristic results in the typical parabolic shape exhibited by the merging events in the space–time plot shown in figure 1. Another feature that is clearly illustrated in figure 11 is the increase in the tilting angle of the inflows as the DW approach one another. This angle is initially very small (see e.g. panel for $t = 730$) but increases rapidly as the mutual interaction between DW becomes stronger, reaching a value of approximately $30\,^{\circ}$ with respect to the radial axis by the time when the merging of the DW takes place (see panel for $t = 836$).

Figure 11. Colour maps of $u$ illustrating the simultaneous occurrence of two merging events at the beginning of the VMS regime. They correspond to a simulation performed at constant $El = 0.32$, whose space–time diagram was shown in figure 1(c). The dashed (green) rectangle in the latter figure indicates the time window that is illustrated in the current figure. Positive (negative) values of $u$ are shown as red (blue), whereas $u=0$ is shown as grey. In each panel there are ${\rm four}$ positive and ${\rm four}$ negative contours evenly distributed in $u \in [-0.071,0.013]$. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction.

After the merging events are completed, a flow pattern characterized by four DW emerges (see panel for $t = 865$), which retains a discrete translational symmetry along the $z$-direction (i.e. the flow pattern remains invariant if it is shifted by $2{\rm \pi}$ in the axial direction). As can be seen from figure 1(c), this flow pattern undergoes subsequent merging events, which occur again simultaneously, when $t$ is between $865$ and $2320$. After this second pair of merging events, the axial symmetry of the flow is fully broken (not shown) and the dynamics of the distinct DW is fully decoupled. It is only after the latter happens that the sequence of merging and splitting events begins and the number of DW in the system may either grow or decay. This initial cascade of merging events leading to a complete symmetry breaking is always observed in the simulations at the beginning of the VMS regime and can thus be interpreted as a transitional stage where the coupling between DW is fully broken. It must be finally noted that the time scale associated with merging events is highly variable and depends crucially on the distance between the cores of the two DW undergoing the merging event. They may occur very slowly, extending over nearly $1000$ advective time units, such as the merging process exemplified in figure 11, or very quickly, within $50$ to $100$ advective time units, like the ones discussed in the next paragraph.

A central aspect of the dynamics in the VMS regime is the emergence of transient chaotic motion interspersed between merging and splitting events. As seen in figure 1, this chaotic dynamics appears localized in the vicinity of the region where a merging event has taken place and it is characterized by the repeated emergence of closely spaced, weak vortices which merge shortly after they form, creating a quick succession of merging and splitting events. The characteristic cycle of this transient dynamics is exemplified in figure 12, which shows the evolution of the flow pattern in a narrow temporal window of the simulation performed at $El = 0.32$ (see the solid line brown rectangle in figure 1c). In this example, the chaotic dynamics occurs in the central part of the system, framed by a dashed (orange) rectangle in the figure, and has little effect on the DW located outside this region. When a merging event is accomplished (see upper panel), the energy released by the DW which is eliminated is transferred to an irregular wavy motion. This wave transports the energy axially upwards and downwards from the location where the merging took place and results in the formation of new vortices (the inflows of the newly created vortex pairs are indicated with arrows in the intermediate panel). The amount of time it takes from the end of the merging event until the new vortices are fully formed normally ranges between $20$ and $30$ advective time units. The new vortices are, however, just a small distance apart from each other, so that they undergo a strong attractive interaction and quickly merge (see lower panel). The energy released after the new merging events is again redistributed and the process just described starts over again. The time span between consecutive merging events is of the order of $100$ advective time units, whereas the entire time period over which the chaotic dynamics typically extends ranges from $1000$ to $2000$ advective time units.

Figure 12. Colour maps of $u$ illustrating the flow patterns associated with the chaotic dynamics that appears interspersed between slow merging events and splitting events. The dynamics displayed in this figure corresponds to the simulation with constant $El=0.32$, and more specifically, to the time window shown as a solid line (brown) rectangle in figure 1(c). Positive (negative) values of $u$ are shown as red (blue), whereas $u=0$ is shown as grey. There are ${\rm four}$ positive and ${\rm four}$ negative contours evenly distributed in $u \in [-0.073,0.013]$. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction.

The main splitting events (i.e. events where a new, strong and persistent DW is created) are normally preceded by chaotic motion and take place when the distance between newly created DW, as well as the distance between these and the other DW in the system, become sufficiently large so that their mutual interaction is weak. An example of a splitting event taking place between $t = 5300$ and $t = 5470$ in the simulation performed at $El=0.32$ (see the dash-dotted purple rectangle in figure 1c) is shown in figure 13. The splitting event occurs in the region marked by the (red) dotted rectangle. The DW that appears within this region in the upper panel of the figure (indicated by a leftwards arrow) is the result of a merging event which has just completed. The other two DW in the figure (which are enclosed in a green dashed rectangle) are moving towards each other as a part of a merging process that will be completed after the splitting event takes place. Hence, the distance between the core of the topmost DW and those of the other two DW becomes increasingly large with time. This enables that when a new DW appears in the space left between them (see intermediate panel, where the new DW is indicated by a rightwards arrow), it can be sufficiently far apart from its neighbours to avoid a strong attractive interaction. As a result, the new DW does not undergo any merging events shortly (in contrast to what happens during the transient chaotic motion) and its strength increases until it becomes comparable to that of the other DW in the system (see lower panel). The new DW is connected to the topmost DW through the outflows, which is an indication that these DW will eventually undergo a merging event. However, such event happens at $t \sim 5930$ (see figure 1c), nearly $500$ advective time units after the new DW first appeared. In general, the DW created after primary splitting events are found to persist for several hundred advective time units, as opposed to the vortices created during the chaotic dynamics which never persist longer than a few tens of advective time units.

Figure 13. Colour maps of $u$ illustrating the occurrence of a splitting event. The example shown corresponds to the simulation with constant $El=0.32$ and it is indicated as a dash-dotted (purple) rectangle in figure 1(c). Positive (negative) values of $u$ are shown as red (blue), whereas $u=0$ is shown as grey. There are ${\rm four}$ positive and ${\rm four}$ negative contours evenly distributed in $u \in [-0.077,0.015]$. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction.

The dynamics illustrated by these examples repeats successively with time, creating a chaotic regime characterized by continuous changes in the number of DW (the regime that has been dubbed VMS regime). Power spectral characterization of this flow regime is shown in figure 14. The spectra were computed by applying a fast Fourier transform to time series of the radial velocity at the mid-gap obtained from simulations where $El$ was held constant. The PSD increases at the lowest frequencies until it reaches a maximum (indicated as a dashed orange line in the figure). For the range of $El$ values investigated, the frequency at which the maximum takes place is close to $f_e$ and increases with increasing $El$, from $0.87 f_e$ at $El = 0.32$ to $1.36 f_e$ at $El = 0.5$. After the peak, the PSD decreases monotonically and exhibits to a good approximation a power law decay range with a decay rate of $-3$ (the best fit to the data yields exponents ranging between $-2.82$ and $-3.33$). This exponent is in agreement with the universal spectral decay rate that was theoretically predicted for EIT (Fouxon & Lebedev Reference Fouxon and Lebedev2003), which has been recently verified in experiments (Yamani et al. Reference Yamani, Keshavarz, Raj, Zaki, McKinley and Bischofberger2021), and thereby identifies the VMS regime as a class of elasto-inertial turbulent states. It must also be remarked that the power spectra obtained at different $El$ values are similar. This observation suggests that flows in the VMS regime are not significantly affected by changes in $El$.

Figure 14. PSD obtained from simulations in the VMS regime at three distinct values of $El$. The left panels show the spectra in a linear–log plot, with the frequency normalized with the elastic frequency, $f_e$, whereas the right panels show the spectra in log–log scale and the frequency is non-dimensionalized with the inverse of the advective time.

3.5. Influence of the domain size and other computational parameters

Since VMS events arise from interactions among DW, it is natural to wonder about the impact that changes in the length-to-gap aspect ratio $\varGamma$ (and consequently in the number of DW the system may contain) may have on the findings reported in the previous subsections. To investigate this aspect, a set of simulations where $\varGamma$ was varied between $9$ and $16$ has been conducted. The protocol followed in these simulations is the same as that described at the beginning of § 3.1: they are started from a Newtonian TVF and $El$ is slowly increased with time at a uniform rate, $El = 10^{-3}t/Re$, where $Re$ was again fixed to $95$. The influence of the initial condition has also been examined. To this end, two or three simulations have been performed for each value of $\varGamma$ considered and these have been started from states having a different number of vortex pairs. A complex spatio-temporal dynamics consistent with the VMS regime has been observed in all cases. However, the $El$ threshold at which merging events first appear changes significantly depending on the domain size and the initial condition. This is illustrated in figure 15, which shows space–time diagrams of $u$ at the mid-gap along the $z$ direction for simulations where (a) $\varGamma = 10$ and the initial condition has ${\rm five}$ vortex pairs, (b) $\varGamma = 14$ and the initial condition has ${\rm eight}$ vortex pairs and (c) $\varGamma = 14$ and the initial condition has ${\rm six}$ vortex pairs. Note that, analogously to figure 1(a), time has been replaced by its corresponding $El$ value. As seen, the first merging events are accomplished at notably different values of $El$ in each case: $El \approx 0.16$ in (a), $El \approx 0.21$ in (b) and $El \approx 0.28$ in (c), which in turn differ from the onset of merging events reported in § 3.1 ($El \approx 0.32$ for $\varGamma = 12.56$, using a state with ${\rm six}$ vortex pairs as initial condition). It is therefore evident that this feature is very sensitive to the domain size and the initial condition used in the simulations. The same is true for the onset of chaotic motion. For the $\varGamma$ values and initial conditions investigated, the first occurrence of a merging event has been found to range between $El \approx 0.15$ and $El \approx 0.32$, whereas chaotic motion has been first observed at $El$ values ranging from $0.25$ to $0.42$.

Figure 15. Space–time diagrams showing the magnitude of $u$ at the mid-gap along the $z$ direction in simulations where the axial length of the system and/or the number of vortex pairs of the initial condition were varied with respect to the simulations shown in figure 1. The cases illustrated correspond to (a) $\varGamma = 10$ using an initial condition with ${\rm five}$ vortex pairs, (b) $\varGamma = 14$ using an initial condition with ${\rm eight}$ vortex pairs and (c) $\varGamma = 14$ using an initial condition with ${\rm six}$ vortex pairs. As in the simulations presented in figure 1, $El$ has been slowly increased with time ($El = 10^{-3}t/Re$) and the latter has been replaced by its corresponding $El$ value in the horizontal axis of the space–time plots. Red and blue areas indicate outflows and inflows respectively. Note that periodic boundary conditions are used in $z$.

The characteristics of the transition towards the VMS regime are not significantly altered by changes in $\varGamma$ or the initial condition. As $El$ increases initially from the Newtonian limit, a centrifugally dominated regime is identified in all cases, where the axial extent of the inflows (outflows) gradually decreases (increases) with increasing $El$. This behaviour continues until the elastic instability sets in abruptly and the intensity of the inflows becomes much higher. This can be seen in figure 15 as a sudden change in the colour intensity that happens at $El \approx 0.125$. It is interesting to note that, despite the significant variation in the $El$ values at which the VMS events occur, the threshold of the elastic instability remains nearly unchanged with varying $\varGamma$ or initial condition (it is always found to occur at $El$ values ranging from $0.12$ to $0.13$). Another feature that is shared by all simulations regardless of the domain size and initial condition is the existence of an initial stage of merging events, with some of them taking place simultaneously, that precede the appearance of chaotic motion. A previously unreported event, where three DW merge simultaneously, has been observed in some simulations within this initial stage (see figure 15c). This particular type of merging event occurs when there is a group of ${\rm three}$ equal DW in which the upper and lower DW move towards the central DW. The forces exerted by the upper and lower DW on the central DW are equal and act in opposite senses, so that the central DW does not move and eventually the three DW merge simultaneously. The main difference between the transition scenarios illustrated in figure 15 and that described in the previous subsections is the absence of the regime characterized by the small axial oscillations of the flow pattern. This regime preceded the onset of the VMS regime in the simulation for $\varGamma = 12.56$. The same axial oscillations are observed in figure 15 at large $El$ values. However, merging and splitting events occur prior to the emergence of the oscillations in these cases. This clearly indicates that the occurrence of these oscillations is not a necessary condition for the VMS dynamics to exist, but an additional dynamics that occurs at large $El$ values and may or not coexist with the VMS events.

The strong dependence of the onset of the VMS events on the initial condition implies that these are highly nonlinear phenomena. It is thus rational to expect that hysteretic behaviour will be observed in the simulations if the control parameter is varied in the reverse direction, i.e. if $El$ is decreased with time. To examine this possibility, we have conducted a simulation where $El$ has been decreased with time at the same rate as it was increased in the previous simulations. The simulation was initiated from the flow state obtained at $El = 0.392$ in the simulation presented in § 3.1 (indicated by a dashed line in figure 1a) and the same parameter values as in § 3.1 were used. To facilitate the description, we denote the simulation in which $El$ increases (decreases) with time as forwards (backwards) simulation. Figure 16 shows the space–time diagram corresponding to the backwards simulation. It becomes apparent from the comparison between this figure and the space–time diagram of the forwards simulation (figure 1a,b) that there is a strong hysteresis. The flow in the backwards simulation eventually returns to the initial state of the forwards simulation (a TVF with six pairs of vortices), but it follows a completely different path, characterized by VMS events and flow states that differ from those observed in the forwards simulation. It is worth noting that the initial cascade of merging events that precede the onset of chaotic motion in forwards simulations is absent in the backwards simulation. This reflects the irreversible character of the symmetry-breaking processes that take place over this initial stage of the VMS regime. Another notable difference is observed in the transition between the centrifugally and elastically dominated regimes. This occurs at a lower $El$ value ($El \approx 0.09$) than in the forwards simulations ($El \approx 0.125$) and not only entails a sudden change in the strength of the vortices, but also a change in their number (from three vortex pairs in the elastically dominated regime to six vortex pairs in the centrifugally dominated regime). Once the centrifugally dominated regime is achieved, flow states obtained in the forwards and backwards simulations for the same $El$ values become identical. This observation reflects the linear nature of the centrifugal instability and evidence that the hysteretic effects observed in the simulations are solely due to the subcritical character of the DW.

Figure 16. Space–time diagram showing the magnitude of $u$ at the mid-gap along the $z$ direction in a simulation where $El$ was slowly decreased with time at the same rate as it was increased in the simulations previously presented ($El = 10^{-3}t/Re$). The parameters used in this simulation are the same as those used in the simulation shown in § 3.1. The initial condition corresponds to the state obtained in the simulation where $El$ was slowly increased with time for $El = 0.392$ (indicated by an orange dashed line in figure 1a). The flow states observed when $El$ decreases differ from those obtained when $El$ increases, evidencing the existence of hysteresis. Red and blue areas indicate outflows and inflows, respectively. Note that periodic boundary conditions are used in $z$.

We next examine whether the occurrence of VMS events depends on the extensibility of the polymers used in the simulations. To this end, we have conducted a set of simulations where the maximum polymer extension $L$ was varied between $30$ and $250$ while keeping the other parameters as in § 3.1. It has been found that the choice of $L$ has an important impact on the characteristics of the elastic instability. This is illustrated in figure 17, which shows the contribution of the polymers to the integral energy budget ($\varPi _e$) in the range of $El$ values where the elastic instability takes place. As shown earlier in figure 6, the onset of the elastic instability leads to a marked increase in the value of $\varPi _e$, which replaces the energy production due to inertia ($\mathcal {P}$) as the leading term that balances the viscous dissipation ($\epsilon$). This characteristic is absent in simulations where the extensibility of the polymers is low (see $L = 30$ case in the figure). In these simulations, $\varPi _e$ increases very gradually with increasing $El$ and remains negligible with respect to $\mathcal {P}$ and $\epsilon$ for the entire range of $El$ investigated (up to $El = 0.5$). This implies that at these elasticity levels the elastic instability does not occur in these cases. As a result, the VMS dynamics does not take place and the flow at high $El$ values is characterized by elastically modified Taylor vortices (not shown). A clear increase in $\varPi _e$ consistent with the emergence of an elastic instability is observed in the simulations when these are performed using $L \gtrapprox 70$. The $El$ threshold at which the instability sets in decreases slightly with increasing $L$ (from $El \approx 0.135$ for $L=70$ to $El \approx 0.115$ for $L=250$). The transition between the centrifugally and elasticity dominated regimes is initially rather smooth (see $L = 70$ case) but it becomes increasingly sharp as $L$ increases. As the transition gets sharper the magnitude of $\varPi _e$ increases substantially, leading to increasingly strong vortices and causing the emergence of spatio-temporal dynamics right after the transition in cases where very large extensibility is considered (see $L=250$ case, where oscillations in the value of $\varPi _e$ arise simultaneously with the transition).

Figure 17. Contribution of the elastic stresses to the integral kinetic energy budget $\varPi _e$ as a function of the maximum polymer extension $L$ in the range of $El$ values where the transition between the centrifugally and elastically dominated regimes takes place.

The onset of the VMS dynamics (which has been observed in the simulations where $L \geq 90$) also takes place at smaller values of $El$ as $L$ increases. It is interesting that, although the elastic instability is observed for $L = 70$, the space–time diagram of this simulation (shown in figure 18a) does not show any sign of spatio-temporal complexity. This might reflect that the emergence of VMS events requires elasticity levels higher than those simulated here when $L$ is close to the critical value for which the elastic instability emerges. The most striking difference in the characteristics of the VMS regime with respect to those observed in the previous simulations occurs when highly extensible polymers are used. Here, after the initial cascade of merging events is accomplished, the dynamics is characterized by a sequence of VMS events that exhibits a clear periodicity (see figure 18(b), which shows the space–time diagram for the simulation using $L = 250$). These structures are closely reminiscent of the flame-like patterns observed in previous studies (Thomas et al. Reference Thomas, Sureshkumar and Khomami2006, Reference Thomas, Khomami and Sureshkumar2009; Liu & Khomami Reference Liu and Khomami2013), with the difference that, in these studies, the flow was non-axisymmetric and the flame-like dynamics was superposed with a rotating wave, whereas in the present study the flow is axisymmetric and therefore the rotating wave is absent.

Figure 18. Space–time diagram showing the magnitude of $u$ at the mid-gap along the $z$ direction in simulations where the maximum polymer extension was set to $L=70$ (a) and $L=250$ (b). The rest of the parameters were kept as in § 3.1. Red and blue areas indicate outflows and inflows, respectively.

We have finally investigated the effect of varying the rate at which $El$ is increased in the simulations. The increase of the $El$ value with time ($El = \alpha t/Re$) can be understood as a continuous perturbation that is imposed on the system, where the parameter $\alpha$ (the rate of increase) regulates the amplitude of the perturbation. We have conducted simulations varying $\alpha$ between $7.5 \times 10^{-4}$ and $5 \times 10^{-3}$, while keeping the other parameters as in § 3.1. The VMS regime (with characteristics similar to those reported in the previous subsections) was found in the simulations where $7.5 \times 10^{-4} \lessapprox \alpha \lessapprox 1.25 \times 10^{-3}$. However, when $\alpha$ was set to higher values, the VMS dynamics did not occur. Panel (a) of figure 19 illustrates the space–time plot as a function of $El$ for a simulation where $\alpha = 1.6 \times 10^{-3}$. As seen, the variation of the vortex pattern with increasing $El$ is initially analogous to that observed when $\alpha = 1 \times 10^{-3}$ (figure 1a). The transition between the centrifugally and elasticity dominated regimes taking place at $El \approx 0.12$ is clearly identified by the sudden change of the colour intensity of the vortices. Also as in the simulation for $\alpha = 1 \times 10^{-3}$, the vortex pattern becomes unstable at $El \approx 0.29$, resulting in small axial oscillations of the DW. These oscillations persist until several merging events take place simultaneously and break the axial symmetry of the flow pattern (which occurs at $El \approx 0.36)$. However, the flow regime that emerges after the symmetry-breaking process differs from the VMS regime. The DW remain at approximately the same axial positions and their number does not change with increasing $El$, nor with time if $El$ is held constant, as shown in panel (b) of the same figure. Similar to the VMS regime, localized transient chaotic dynamics is often observed in this flow regime (see regions enclosed by the dashed red rectangles). As shown in figure 20, a particular feature of the flow structure in this flow regime is the emergence of small scale vortices near the inner cylinder (see regions demarcated by the dashed green rectangles). These vortices are consistent with the elastic Görtler vortices identified by Song et al. (Reference Song, Teng, Liu, Ding, Lu and Khomami2019) at higher Reynolds numbers. The finding of a flow regime distinct from the VMS regime in these simulations reflects a well-known feature of viscoelastic TCF: the coexistence of various flow regimes for the same values of the control parameters. The simulations may converge to one regime or the other depending on the amplitude and type of the perturbations imposed. Our study shows that, to capture the VMS in the simulations, the amplitude of the perturbation cannot be too large, otherwise the VMS regime is replaced by the flow regime just described.

Figure 19. Example of the dynamics observed in the simulations when $\alpha$ (i.e. the rate of increase in $El$) is larger than $1.25 \times 10^{-3}$. (a) Space–time diagram showing the magnitude of $u$ at the mid-gap along the $z$ direction in a simulation where $\alpha = 1.6 \times 10^{-3}$. The rest of the parameters are as in § 3.1. (b) Space–time plot obtained when $El$ is held constant at $0.5$ in the simulation shown in (a). The dashed (red) rectangles demarcate regions of transient chaotic dynamics. Red and blue areas indicate outflows and inflows, respectively.

Figure 20. Colour map of $u$ illustrating the instantaneous flow structure at $El = 0.5$ in a simulation where $\alpha = 5\times 10^{-3}$. Positive (negative) values of $u$ are shown as red (blue), whereas $u=0$ is shown as grey. There are $10$ contours evenly distributed in $u \in [-0.065,0.015]$. The system is shown rotated by $90\,^{\circ}$ in the counterclockwise direction. Dashed (green) rectangles are used to highlight the small scale vortices that emerge near the inner cylinder.

3.6. Variation of the angular momentum transport with increasing fluid elasticity level

An important question raised by the observations above is why the dynamics of the distinct vortex pairs decouples when the polymer elasticity exceeds a certain threshold. It is well known that in two-dimensional vortex systems the conservation of angular momentum imposes strong restrictions on the motion of the vortices and the mean separation among them remains generally nearly constant (Aref Reference Aref1983; Batchelor Reference Batchelor2000). It may therefore appear surprising that vortex pairs in the VMS regime can freely move through the system, either toward each other or away from each other, without drastic changes in the system's energy. To explain this seeming inconsistency, it is instructive to examine the impact of increasing $El$ in the different contributions to the flux of angular momentum ($J^{\omega }$) across the annular gap. In a viscoelastic fluid flow, $J^{\omega }$ can be split into three terms

(3.5)\begin{equation} J^{\omega} = r^3 \left[\underbrace{\overline{u\omega}}_{J^{\omega}_c} - \underbrace{\frac{\beta}{Re} \partial_r\bar{\omega}}_{J^{\omega}_d} - \underbrace{\frac{(1-\beta)}{Re}\frac{\overline{T_{r \theta}}}{r}}_{J^{\omega}_p}\right], \end{equation}

where $\omega = v/r$ is the angular velocity and the bar symbol indicates averaging over the axial direction (for $El$ values where the flow is steady) or over the axial direction and time (when the flow is non-steady or chaotic). Here, $J^{\omega }_c$ denotes the convective transport of angular momentum, which is associated with the vortices, $J^{\omega }_d$ is the diffusive transport due to viscosity and $J^{\omega }_p$ is the angular momentum transport caused by polymer stresses.

Although the above equation was originally derived for turbulent flow (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007; Song et al. Reference Song, Teng, Liu, Ding, Lu and Khomami2019), it is straightforward to show that it also applies to steady vortex flow (see Appendix B for a step by step derivation under steady and axisymmetric conditions). Hence, it can be used to study the variation in the contributions of the different transport mechanisms as $El$ increases from the Newtonian limit up to the largest value simulated within the VMS regime. This is shown in figure 21 for a radial location near the mid-gap using the data obtained from the simulation presented in § 3.1 (similar analyses for other simulations are shown in Appendix C). Since the dynamics in the VMS regime is chaotic and $J^{\omega }$ is a fluctuating quantity, several additional simulations at constant $El$ had to be performed in order to obtain meaningful values of $J^{\omega }$ and its contributing terms in this flow regime. The initial conditions for these simulations were flow states obtained in the simulation with slowly varying $El$. Starting from these solutions, the $El$ values were fixed and the chaotic flow was simulated in each case for $20000$ advective time units. The momentum fluxes were then calculated by averaging over this long time period.

Figure 21. Variation of the axially and time averaged angular momentum current ($J^{\omega }$) and its contributing terms, $J^{\omega }_c$, $J^{\omega }_d$ and $J^{\omega }_p$, obtained near mid-gap, as $El$ increases. The (orange) vertical dotted lines demarcate the different flow regimes observed in the present study: I. the centrifugally dominated regime, II. the elastically dominated regime characterized by steady patterns of equally spaced DW, III. the elastically dominated regime characterized by spatio-temporal dynamics.

As seen in the figure, three clear stages can be distinguished in the behaviour of the angular momentum fluxes as $El$ increases, which are consistent with the different flow regimes identified throughout our study. In the first stage, coinciding with the centrifugally dominated regime, $J^{\omega }$ remains nearly constant. While its main contribution stems from the diffusive transport ($J^{\omega }_d$), the convective flux ($J^{\omega }_c$) also plays an important role close to the Newtonian limit ($El \to 0$). However, due to the dissipative nature of the polymer activity in this flow regime, the contribution of the vortices to $J^{\omega }$ decays with increasing $El$ and it is gradually replaced by the angular momentum flux due to the polymer stresses ($J^{\omega }_p$). Note that, although it may seem surprising that the contribution of $J^{\omega }_d$ exceeds that of $J^{\omega }_c$ in a Newtonian (or low $El$) vortex flow, this happens because the simulation is conducted at a Reynolds number ($Re=95$) which is very close to the onset of the Taylor vortices ($Re = 89$). Here, the vortices are still weak and the momentum transport is dominated by molecular diffusion (as in the laminar regime). As $Re$ increases, the contribution of $J^{\omega }_c$ near the mid-gap becomes increasingly large and the contribution of $J^{\omega }_d$ decreases, so that the former eventually dominates the momentum transport in this regime (not shown). The onset of the elasticity dominated regime ($El \sim 0.12$) is accompanied by an abrupt increase in $J^{\omega }_p$, which becomes the leading contribution to the angular momentum flux. Moreover, the magnitude of $J^{\omega }_p$ keeps increasing as $El$ increases, leading to a substantial increase in $J^{\omega }$ with respect to that in the centrifugally dominated regime. The convective and diffusive fluxes, on the other hand, exhibit a slight increase and decrease, respectively, when the elastic instability sets in, and subsequently decrease very gradually with increasing $El$. The onset of the spatio-temporal dynamics ($El \sim 0.29$) brings a significant drop in the total angular momentum flux. This is mainly associated with an initial decay in $J^{\omega }_p$ that takes place during the initial cascade of merging events where the different vortex pairs become fully decoupled. Despite its initial decrease, $J^{\omega }_p$ is still the main contribution to $J^{\omega }$, and after the initial phase of merging events, its magnitude increases with increasing $El$ over the entire VMS regime. The convective contribution $J^{\omega }_c$ also decays during the initial phase of merging events. However, it appears to fluctuate around a constant value, $J^{\omega }_c \approx 0.013$, with increasing $El$. A similar behaviour is observed in $J^{\omega }_d$, which increases initially with the emergence of the VMS dynamics but remains subsequently nearly constant as $El$ increases.

The analysis of the angular momentum fluxes just presented allows us to propose an answer to the question posed at the beginning of this section. As shown, polymer stresses are very efficient in transporting angular momentum (provided that the level of elasticity in the working fluid is sufficiently large) and so the contribution of the vortices in this regard, which is essential in the Newtonian case, is only marginal in the viscoelastic case. The amount of momentum carried by the vortices becomes increasingly small as $El$ increases until it eventually reaches a nearly constant level, $J^{\omega }_c \approx 0.013$, in all simulations where the VMS is found. On the basis of this observation, we suggest that, when the angular momentum carried by the vortices drops to this level, the constraints imposed by the angular momentum conservation on the vortices break and these may decouple from the rest of the system. This limit would mark the beginning of the VMS dynamics and could also be interpreted as the minimum amount of angular momentum that the DW must carry to form a pattern of steady vortices. The latter interpretation offers an explanation for the question of why the DW do not merge at lower $El$ values. As noted in § 3.1, arrangements of equally spaced, steady DW have not been so far experimentally reported. In fact, it has been inferred from the experiments that two DW coalesce when the distance between them is lower than $5d$, a characteristic that would render the formation of these arrangements of DW unfeasible. The reason for this apparent contradiction may lie in the fact that these experiments were conducted at low values of $Re$, where the flow in the Newtonian case is centrifugally stable (i.e. molecular diffusion suffices to transport angular momentum). The amount of angular momentum carried by the DW at these low $Re$ is expected to be very small and hence it is reasonable to assume that it might be below the threshold required to observe these arrangements. Another important remark concerning the angular momentum fluxes in the VMS regime is that the two Newtonian contributions ($J^{\omega }_c$ and $J^{\omega }_d$) remain nearly constant over the entire regime. The increase in the angular momentum flux taking place in this flow regime is thus only due to the contribution of the polymer stresses, which continue increasing with increasing $El$. This circumstance strongly suggests that the dynamics in the VMS regime is fully dominated by elastic effects (as already noted in the experimental study by Lacassagne et al. Reference Lacassagne, Cagney, Gillissen and Balabani2020) and raises the question about a possible relationship between flow states in the VMS regime and those driven by pure elastic instabilities in the inertialess limit.