Vortex merging and splitting events in viscoelastic Taylor Couette flow

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 $\Gamma=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 Taylor vortex flow 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., 2020). 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 Taylor vortex flow 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$ values. 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 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 decouple and these begin to move independently in the axial direction. It will be shown later in section 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 the panel (b) of the figure 1. 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.

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 the panel (c) of the figure 1, 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 the panel (a) of the figure 1 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 Taylor vortex flow.

Figure 2 shows colormaps 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$ degrees 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 Taylor vortex flow 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).

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, 1998). To verify this hypothesis quantitatively, figure 3 shows axial profiles of the elastic force, hereafter denoted as $F_{e}$ (left panels), and its associated work $F_{e}u$ (right panels), 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 the right panels), which enhances the fluid motion in the inflows and create 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.

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 the 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. The 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 the 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 behavior. 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).

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 the panel $(a)$ of figure 4, 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 panel $(b)$ in figure 4), 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 color 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 et al., 2010; Dubief et al., 2013),

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

Here, the overline denotes axially averaged quantities, $S^{\prime}=(\nabla\mathbf{v}^{\prime}+\nabla\mathbf{v}^{\prime 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 ($\overline{\mathbf{v}}=(0,\overline{v},0)$ and $\overline{\mathbf{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 the appendix A). The first and second integrals in equation 4 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.

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 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 diwhirls (DW) (Groisman & Steinberg, 1997; Lange & Eckhardt, 2001; Thomas et al., 2006, 2009), due to its similarity with 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, 1997; Lange & Eckhardt, 2001; Thomas et al., 2006; Lacassagne et al., 2020). 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 elasticiy 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 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 is 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 section 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 oscillatory strips (Groisman & Steinberg, 1996; Thomas et al., 2006, 2009).

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 color maps of $u$ taken at five equally spaced time instants within a period (shown as circles in the panel $(a)$ of figure 9). 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 degrees 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 $A$), it is observed that the inflows move first axially upwards (panel $B$), reach the position of maximum displacement (panel $C$) and subsequently move downwards (panel $D$), returning eventually to the initial state (panel $E$, which is identical to $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 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 the panel 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 the panel D) of the outflow remains as in the stationary case, whereas the lower half ((purple) dash-dotted rectangle in the panel 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 the states B and D shown in the figure).

The frequency of the oscillation was determined by applying the Fast Fourier transform to a time series of the axial velocity $w$ obtained at a radial location close to the outer cylinder (figure 9 $(a)$). 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. (2020), $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 $6$ vortex pairs. For other flow patterns with different number 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., 2020) 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.

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\pi$ in the axial direction). As can be seen from the panel (c) of figure 1, 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 are 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 the 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 the 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 the panel (c) of the figure 1). 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.

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 the panel (c) of the figure 1) 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 neighbors 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 is an indication that these DW will eventually undergo a merging event. However, such event happens at $t\sim 5930$ (see panel (c) in the figure 1), 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.

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 the figure 14. The spectra were computed by applying the fast Fourier transform to time series of the radial velocity at the mid-gap obtained from simulations where $El$ was held constant. The power spectral density 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.87f_{e}$ at $El=0.32$ to $1.36f_{e}$ at $El=0.5$. After the peak, the power spectral density 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 elasto-inertial turbulence (Fouxon & Lebedev, 2003), which has been recently verified in experiments (Yamani et al., 2021), and thereby identifies the VMS regime as a class of elasto-inertial turbulent states. It must also be remarked the similarity of the power spectra obtained at different $El$ values. This observation suggests that flows in the VMS regime are not significantly affected by changes in $El$.

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 $\Gamma$ (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 $\Gamma$ 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 the section 3.1: they are started from a Newtonian Taylor vortex flow and the $El$ number 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 $\Gamma$ 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 the figure 15, which shows space time diagrams of $u$ at the mid-gap along the $z$ direction for simulations where (a) $\Gamma=10$ and the initial condition has $5$ vortex pairs, (b) $\Gamma=14$ and the initial condition has $8$ vortex pairs and (c) $\Gamma=14$ and the initial condition has $6$ vortex pairs. Note that analogously to the 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 the $El$ number 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 the section 3.1 ($El\approx 0.32$ for $\Gamma=12.56$, using a state with $6$ 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 $\Gamma$ 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$.

The characteristics of the transition towards the VMS regime are not significantly altered by changes in $\Gamma$ 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 the figure 15 as a sudden change in the color 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 $\Gamma$ 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 panel (c) of figure 15). This particular type of merging event occurs when there is a group of $3$ 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 the 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 $\Gamma=12.56$. The same axial oscillations are observed in the 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 reversed direction, i.e. if the $El$ number 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 the section 3.1 (indicated by a dashed line in the panel (a) of the figure 1) and the same parameter values as in the section 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 (panels (a) and (b) of the figure 1) that there is a strong hysteresis. The flow in the backwards simulation eventually returns to the initial state of the forwards simulation (a Taylor vortex flow with 6 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 3 vortex pairs in the elastically dominated regime to 6 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.

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 section 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 the figure 17, which shows the contribution of the polymers to the integral energy budget ($\Pi_{e}$) in the range of $El$ values where the elastic instability takes place. As shown earlier in the figure 6, the onset of the elastic instability leads to a marked increase in the value of $\Pi_{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, $\Pi_{e}$ increases very gradually with increasing $El$ and remains negligible with respect to $\mathcal{P}$ and $\epsilon$ for the entire range of $El$ numbers 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 $\Pi_{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 $\Pi_{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 $\Pi_{e}$ arise simultaneously with the transition).

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 the panel (a) of the figure 18) 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 the panel (b) of the figure 18, 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., 2006, 2009; Liu & Khomami, 2013), 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.

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\cdot 10^{-4}$ and $5\cdot 10^{-3}$, while keeping the other parameters as in section 3.1. The VMS regime (with characteristics similar to those reported in the previous subsections) was found in the simulations where $7.5\cdot 10^{-4}\lessapprox\alpha\lessapprox 1.25\cdot 10^{-3}$. However, when $\alpha$ was set to higher values, the VMS dynamics did not occur. The panel (a) of the figure 19 illustrates the space time plot as a function of $El$ for a simulation where $\alpha=1.6\cdot 10^{-3}$. As seen, the variation of the vortex pattern with increasing $El$ is initially analogous to that observed when $\alpha=1\cdot 10^{-3}$ (panel (a) of the figure 1). 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 color intensity of the vortices. Also as in the simulation for $\alpha=1\cdot 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 the 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 the 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. (2019) 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 it is replaced by the flow regime just described.

An important question raised by the observations above is why the dynamics of the distinct vortex pairs decouple when the polymers 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 (Batchelor, 2000; Aref, 1983). 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 in three terms,

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 et al., 2007; Song et al., 2019), 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 the largest value simulated within the VMS regime. This is shown in the figure 21 for a radial location near the mid-gap using the data obtained from the simulation presented in the section 3.1 (similar analyses for other simulations are shown in the 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.

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 polymers 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 to the question of why the DW do not merge at lower $El$ values. As noted in section 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 continues 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. (2020)) 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.