FLUIDS 3 , 113301 ( 2018 ) Characterization of superimposed instabilities in the planar extensional flow of viscoelastic fluids

We simulate the isothermal, incompressible, time-dependent flow of Boger fluids, using the FENE-CR constitutive model, in the optimized shape cross-slot extensional rheometer [Haward et al., Phys. Rev. Lett. 109, 128301 (2012)]. We uncover a family of predominantly elastic instabilities resulting from the interaction of stationary asymmetric and time-dependent flow transitions. The superposition of these instabilities, with varying degrees of relative amplitude, produces five distinct flow regimes which are classified by their elasticity, suggesting that the dynamic system is situated in the vicinity of a triple point in the inertia-elasticity-viscosity parameter space. A detailed characterization of the various associated modes of instability is provided. Spectral analysis of the first component of the velocity vector suggests that the flow regimes featuring loss of birefringence strand integrity become chaotic and locally turbulent near the center of the cross-slot, above a certain flow rate threshold.


I. INTRODUCTION
Viscoelastic fluids have been known to display a wide array of flow instabilities in several geometries, as discussed in early literature surveys by Larson [1] and Shaqfeh [2].One such geometry is the microscale planar cross-slot, which became the subject of intensified research after the discovery by Arratia et al. [3] of a stationary elastic instability, wherein at sufficiently high Deborah number De, the flow field becomes asymmetric as a result of accumulated extensional stress along the outlet centerline of the slot, forming a thin layer of large normal stresses, the so-called birefringence strand.Further increases in flow rate cause a transition to time-dependent flow.These instabilities were shortly after confirmed to be purely elastic in nature, by a set of inertialess simulations using the upper-convected Maxwell model [4].Earlier computational work in the cross-slot device (e.g., Harris and Rallison [5] and Remmelgas et al. [6]) was performed on a truncated mesh, such that only a quarter or half of the full cross-slot was simulated with imposed symmetry boundary conditions, thus precluding the emergence of flow asymmetries.Similarly, analytical studies such as the work of Öztekin et al. [7] focused only on a portion of the extensional flow equivalent to one of the cross-slot outlets, in order to assess the characteristics of the time-dependent flow transition in the absence of asymmetries, which had not been reported at the time.Experimentally, the first indirect evidence of flow asymmetries can be found in the work of Gardner et al. [8], who reported asymmetric velocity profiles along cross-slot outlet channels (cf.Fig. 6 in Ref. [8]).
Since the work by Arratia et al. [3], renewed interest has emerged in the cross-slot geometry.Further computational work has confirmed that the stationary asymmetric instability may occur for shear-thinning Phan-Thien-Tanner fluids [9] or other bounded extensional viscosity model fluids, such as the Finitely Extensible Nonlinear Elastic (FENE) model with either Peterlin (-P) or Chilcott-Rallison (-CR) closures [10].Rocha et al. [10] have also studied the effect of solvent viscosity ratio, β = η s /(η s + η p ), where η s and η p are, respectively, the solvent and polymer viscosities, and extensibility, L 2 , and concluded that the onset of the stationary asymmetric transition occurs at lower De for lower β and/or higher L 2 .Xi and Graham [11] also simulated the flow of a FENE-P fluid in a cross-slot device and have not found a stationary instability, instead reporting a direct transition to time-dependent flow.Their choice of a high solvent viscosity ratio, β = 0.95, together with the previously mentioned effects of β and L 2 on the onset of stationary instabilities, suggests that a minimum threshold in the relative strength of extensional flow relative to shearing effects is necessary to trigger the steady asymmetric transition.Experimental work performed on cross-slot flows of wormlike micellar solutions [12,13] seems to confirm this hypothesis: highly viscoelastic fluids sequentially undergo both transitions, whereas weakly viscoelastic fluids directly evolve from steady symmetric to unsteady flow as the flow rate increases.This is further supported by the experimental work of Sousa et al. [14] who studied the flow stability of various viscoelastic fluids with different β.Experimental [14] and computational [15] investigations of cross-slot flows with different channel aspect ratios have shown that the stationary asymmetric transition occurs only in channels with sufficient depth, i.e., high aspect ratio, whereas channels which are relatively shallow usually evolve directly from steady symmetric to unsteady flow as De is increased.The relative proximity of the top and bottom bounding walls in shallow channels causes an increase in the relative magnitude of shear effects [14], therefore suppressing the stationary asymmetric transition that is linked to the extensional nature of the flow.However, as the aspect ratio decreases, the second, time-dependent transition occurs at progressively lower Deborah numbers [15], indicating that the additional shear-induced normal stresses near the corners are themselves destabilizing, although leading to a different type of instability.
Whether it be by manipulation of fluid rheology or the geometric aspect ratio of the cross-slot device, evidence suggests extension and shear compete in order to determine the nature of flow instabilities.However, studying the two effects separately has proven to be difficult.According to the Pakdel-McKinley criterion [16], inspired by the micromechanical model proposed by Larson et al. [17], elastic instabilities are a consequence of curved streamlines acting on polymer chains not fully aligned with the streamwise direction, causing them to stretch and produce tensile stress (cf.Fig. 12 in Ref. [17]).Disturbances in the alignment of polymers are thus amplified, eventually leading to unstable flow.Alternatively, as discussed in Ref. [15], a flow which is rotational in nature, despite having curved streamlines, cannot generate the aforementioned differential stretching of polymer chains owing to the uniform angular velocity of the flow field or, for that matter, the lack of extension thought to be the cause of the stationary asymmetric transition.Ideally, isolating the influence of shear and extension in elastic instabilities would require a flow geometry capable of producing a nearly purely extensional flow field while retaining streamline curvature.Attempts made by Cruz et al. [9] and Rocha et al. [10] to reduce the influence of shear by rounding the corners of the cross-slot geometry have also caused the steady asymmetric transition to be postponed to higher flow rates, slightly when the curvature of rounded corners is 5% of the channel width D [9] and pronouncedly for rounding curvature equivalent to 50% of D [10].Thus the extensional character of standard planar cross-slot flows is also predicated on the sharpness of its corners.
Introduced by Alves [18] and experimentally validated by Haward et al. [19,20], the optimized shape cross-slot extensional rheometer (OSCER device), depicted in Fig. 1, is a variant of the cross-slot geometry, numerically optimized to produce an homogeneous extensional field (u x , u y ) = (εx, −εy) with deformation rate ε = 1.5U 0 /(15H ) along a wide region in the center of the geometry [cf.Fig. 1(b) in Ref. [19]] and, consequently, a well-defined homogeneous birefringence strand along the outlet axis.Although initially proposed as an extensional rheometer, the OSCER device has also been used as a platform for the investigation of elastic instabilities.(a) Geometry of the optimized shape cross-slot extensional rheometer (OSCER device).The halfwidth H is used as the characteristic length scale.(b) Detail of the computational mesh.The deformation of computational cells is the result of the application of an optimization algorithm to the mesh of a planar cross-slot, described in Ref. [18].(c) Zoomed-out view of the mesh, showing the inlet and outlet channels.Cells in these channels are progressively larger in the streamwise direction towards the inlets and outlets.measured the birefringence strand formed by the flow of various poly(ethylene oxide) solutions with varying proportions of glycerol and of hyaluronic acid aqueous solutions.Depending on the magnitude of the Elasticity number El, representative of the ratio between elastic and inertial forces, which due to the shear-thinning nature of the test fluids was defined at a flow rate corresponding to the respective critical Weissenberg number Wi crit of each flow transition, flows may either become steady asymmetric for El crit > 1, or undergo a transition into a time-dependent, inertio-elastic regime for El crit < 1, with both Reynolds number Re crit > 10 and Wi crit > 1.While the high elasticity transition is analogous to the stationary asymmetric transition observed in standard planar cross-slots, the low elasticity transition is characterized by what appears to be rapid oscillation of the birefringence strand, in a manner which resembles the behavior of a vibrating taut string [21].Furthermore, Haward et al. [22] studied the flow in the OSCER device using monodisperse, dilute atactic-polystyrene dissolved in dioctyl phthalate.Their test fluids were weakly shear-thinning and despite having a relatively high solvent viscosity ratio, β > 0.5, generate strong yet thin birefringence strands due to the large extensibility of polymer chains.The authors reported a new type of instability [22], wherein the stagnation point becomes laterally displaced and oscillates between the central position of the slot and one of the outlets.Further increases in flow rate cause a symmetry breaking transition, analogous to the steady asymmetric transition in standard cross-slots.However, the resulting flow is unsteady in what appears to be a superposition of the lateral displacement oscillations with a base asymmetric flow field, mixed with aperiodic modes of fluctuation (cf.Supplemental movie S4 in Ref. [22]).
The observation of symmetry breaking transitions in the OSCER device [21,22] indicates that the shearing produced by the corners in a standard cross-slot is not a necessary condition for the underlying stationary asymmetric transition.Furthermore, in the near absence of inertia or shearing effects, of which the latter is an integral part of the micromechanical model of elastic instabilities proposed by Larson et al. [17], hitherto unknown time-dependent instabilities [22] may occur and eventually become superimposed on an asymmetric flow field.In this work, we simulate the twodimensional (2D) flow of dilute, highly extensible, FENE-CR fluids in the OSCER device, under conditions of low to negligible inertia, and attempt to characterize the various emerging transitions and their superimposed states.Considering the high aspect ratio of the OSCER device, specifically AR = depth/width = 10.5 [19], we expect that 2D flow simulations are a reasonable approximation to experimental flows, since the influence of the top and bottom bounding walls should be negligible near the center plane of the device.

II. NUMERICAL METHODS
We simulate the isothermal, incompressible, time-dependent flow of monodisperse, dilute, longchained polymer solutions, modeled by a finitely extensible nonlinear elastic (FENE) dumbbell model [23] with the Chilcott-Rallison (-CR) closure [24].The relevant rheological properties of the model fluids are given in Table I and are based on real Boger fluids of the kind developed by Odell and Carrington [25].The governing equations are respectively the continuity and momentum equations, and the constitutive equations of the FENE-CR model, In the above equations, η s is the solvent viscosity, η p is the polymer shear viscosity contribution, τ is the extra-stress tensor, L 2 is the extensibility of polymer chains, λ 0 is the zero shear-rate relaxation time, and g is an auxiliary expression related to the decrease in effective relaxation time λ eff = λ 0 /g.Since the FENE-CR model is a constant shear viscosity approximation, it is applicable only to the simulation of at most weakly shear-thinning fluids.As mentioned in Table I, some of the rheological data used in this work was kindly provided by the authors of Ref. [22] for other fluids tested but not reported in Ref. [22].A cursory power law fit of the slopped portion of their shear viscosity curves (cf.Fig. 4 in Ref. [22]) indicates that the flow behavior index n varies between 0.95 and 1, confirming a very weak shear-thinning character.The FENE-CR model is more commonly formulated using the conformation tensor A as follows: with the following relation between τ and A: The governing equations are solved using an implicit, second-order, finite volume numerical method using a collocated mesh arrangement, which is described in detail elsewhere [26], adapted from the SIMPLEC algorithm [27].The method includes the CUBISTA high-resolution scheme for the treatment of advective terms [28] in the momentum and constitutive equations, and the log-conformation technique for the transformation of the constitutive equation [29], following the methodology originally developed by Fattal and Kupferman [30], in order to increase the numerical stability.Furthermore, to increase the accuracy, the second-order backward differencing scheme [31] is used for time integration.At the walls, no-slip boundary conditions are assumed for the velocity components, and the extra-stress components are linearly extrapolated from the first two adjacent cells along the wall-normal direction, as described in Ref. [32].Uniform velocity profiles and null extra-stress profiles are assumed at the inlets.Neumann boundary conditions of the form n • ∇φ = 0 are assumed at the outlets, except for pressure, which is linearly extrapolated from the first two upstream cells.Accumulation of round-off error at machine level precision suffices to generate flow instabilities, so no disturbances are artificially introduced.The numerical method is implemented in double-precision FORTRAN.
The block structured mesh, depicted in Fig. 1, comprises a total of 94 233 cells, split over five blocks, one block with 21 008 cells for each of the inlet and outlet channels and a fifth central block with 10 201 cells.For each of the inlet and outlet blocks, the cell length along the streamwise direction is progressively larger towards the inlets and outlets.The center of the device is occupied by a square cell, with side length of 0.0616H.
We define three nominal dimensionless numbers, respectively, the Weissenberg, Reynolds, and Elasticity numbers, where U 0 is both the uniform inlet velocity and the characteristic velocity of the flow, η 0 = η s + η p is the solution shear viscosity, and d h is the hydraulic diameter.Since we simulate flow on a 2D geometry with bottom-top symmetry, the experimental equivalent would be a very deep channel, therefore the hydraulic diameter is twice the channel width, or in this case four times the half-width H, which we define as the characteristic length scale.These dimensionless quantities are designated as nominal since, for time-dependent flow, the effective relaxation time and strain rate at the center of the device are variable.
For each model fluid, we simulate the flow at progressively higher flow rates until the flow is clearly time-dependent.Due to the emergence of very large spatial and temporal gradients, in both velocity and stress fields, error propagation and eventual numerical divergence of simulations are common.We are thus forced to compute using very small time steps, usually t/λ 0 = 0.0001−0.0005,and to iterate five to fifteen times per time step, specifically to allow our pressurecorrection algorithm a higher degree of convergence and to reduce the explicitness of the method.Simulations are performed on an Intel Xeon E5-2643 CPU running at 3.30 GHz.Using a sequential implementation of the numerical method and depending on the time step and number of iterations per time step, the CPU computes 5λ 0 −75λ 0 of real time per simulation per month.Time elapsed between flow start-up and the development of periodic or quasiperiodic flow varies substantially, between 10λ 0 and 50λ 0 .

A. Classification of instability regimes
Based on the qualitative behavior of flow instabilities, the nominal Elasticity number El nom is used to divide the parameter space into five flow regimes.

Very low elasticity, El nom = {1.6, 2.3}
Figure 2 depicts an illustrative example of flow dynamics near the lower end of the elasticity spectrum, at Wi nom = 10.4.The stagnation point oscillates periodically along the horizontal centerline, with average position at the geometric center (x, y) = (0, 0).The amplitude of these fluctuations is small, with peak velocities measured at the center not larger than u x /U 0 ≈ ±0.05.The values of the normalized first normal stress difference N 1 /(η 0 U 0 /2H ) also oscillate in a similar manner, with no loss of birefringence strand integrity.Experimentally, such a strand would appear to brighten and dim periodically, albeit slightly.The periodic character of this flow is further reinforced by the converging phase-plane plot in Fig. 2(d), suggesting that the dynamics of the flow are governed by a limit cycle, as reported elsewhere for standard planar cross-slots [11].Inertia is low yet not negligible, with 2.2 < Re nom < 6.4 for all reported time-dependent flows in this regime.In Fig. 2, and in the remaining figures, t = 0 corresponds to the onset of periodic or quasiperiodic behavior and not the initial simulation time, which can correspond to either the start-up of flow or the previous result of a simulation at a lower flow rate, used as the initial state of the flow field for the new simulation, in order to shorten start-up times.

Low elasticity, El nom = {3.2, 4.1}
An example of flow behavior in the low elasticity regime is depicted in Fig. 3 for Wi nom = 3.4.An initially horizontal birefringence strand undergoes a loss of symmetry and begins rotating, either clockwise or counterclockwise (cf.video S4 [33]).The flow field becomes progressively more asymmetric as time evolves until eventually the degree of asymmetry reaches a critical level and the birefringence strand breaks, followed by the downstream advection of what remains of the formerly contiguous strand.The breaking process is accompanied by the brief appearance of a secondary flow near the center of the geometry.The flow field then returns to a symmetric configuration, and a new birefringence strand forms, resetting the cycle.The amplitude of fluctuations is significant, with normalized velocity at the geometric center peaking at u x /U 0 ≈ ±2 and the normalized first normal stress difference N 1 /(η 0 U 0 /2H ) periodically approaching values near zero.The sequence of events culminating in this loss of strand integrity strongly resembles the start-up of stationary  Oscillations in velocity and normal stresses at the center of the OSCER device are quasiperiodic, albeit with very low amplitude, such that the normalized velocity is mostly contained within the range u x /U 0 ≈ ±0.005.Therefore the position of the stagnation point is nearly constant, and consequently the aforementioned oscillations may be thought of as residual.However, as shown in Fig. 5(b), the flow assumes a locally asymmetric configuration near the center of the OSCER device.Although small, this asymmetry is distinguishable from a similar snapshot of the flow field for the very low elasticity regime, shown in Fig. 5(a).Thus the intermediate elasticity regime is characterized by a quasistationary, locally asymmetric flow field, which again indicates the presence of a secondary instability interfering with the expected behavior of the stationary asymmetric transition.Figure 4(d) suggests a converging limit cycle (cf.video S5 [33]).Attempts to study higher flow rates, namely, Wi nom = {8.2,8.8, 9.4, 10.0}, resulted in the systematic numerical divergence of simulations, both with the method described in this work and with an alternative method using, in lieu of the CUBISTA scheme [28], a WENO implementation [34] for the treatment of advective terms in the momentum and constitutive equations.The characteristic behavior in the high elasticity regime is illustrated in Fig. 6 for Wi nom = 3.6.The periodic oscillations in this regime appear to be a combination of dynamics from the very low elasticity and low elasticity regimes.The time series of normalized velocity is a periodic wave, punctuated by narrow peaks at its maxima and minima, which indicate the breaking of the birefringence strand.However, whereas in previous examples the position of the stagnation point remained near or at the center of the geometry, except during loss of strand integrity, in this regime the horizontal location of the stagnation point oscillates significantly, in the range x/H = ±5.Qualitatively, the motion of the birefringence strand is perhaps best described as analogous to the propagation of a wave in a finite elastic string, or rather a finite elastic sheet since the simulated geometry is equivalent to a deep channel, with the important distinction that the strand periodically breaks and reforms whenever the wave would be reflected by the fixed ends of the sheet.This string-sheet analogy has been previously used to describe the experimental behavior of birefringence strands in the OSCER device [21], albeit referring to the flow of a shear-thinning, 3000 ppm poly(ethylene oxide) solution in water, with rheological properties which are significantly different from those typical of the high elasticity regime, namely, L 2 = 3360, β = 0.26, λ 0 = 0.004s, η 0 = 0.0038 Pa s, tested at a flow rate corresponding to Re = 64.4,whereas the example in Fig. 6 was simulated at Re nom = 0.45.contamination by aperiodic modes of oscillation.At regular intervals, jets of fluid pierce the birefringence strand (cf.video S10 [33]).Consequently, as seen in Fig. 5(c), the strand is fragmented into segments, which upon downstream advection generate an apparent rippling pattern in Fig. 7(a).While the instabilities in lower elasticity regimes are to some extent similar to previous observations in planar cross-slot flows, either individually or as a superposition of known phenomena, the very high elasticity regime is to the best of our knowledge a new observation.Although Haward et al. [22] studied the experimental counterpart of the simulation depicted in Fig. 7, they were unable to describe the instantaneous dynamics of flow in detail.Based on the typical time step of our simulations, t/λ 0 = 0.0001, we estimate that a data acquisition rate in the range of 10-100 kHz may be required for detailed experimental characterization of a fluid with the same properties as aPS16-1400 ppm, specifically the relaxation time.

Stationary asymmetric flow
In addition to the various time-dependent instabilities, a few stationary asymmetric flows were also observed, of which the most prominent example is shown in Fig. 5(d).This was unexpected, since all model fluids have a high solvent viscosity ratio, β > 0.5, and as shown elsewhere [14] for various solutions of poly(acrylamide) in glycerol and water, fluids with β > 0.05 undergo a direct transition to unsteady flow without an intermediary steady asymmetric state.Since the dynamics of cross-slot flows are the result of the relative contributions of extension and shear effects, a geometry such as the OSCER device, optimized to maximize the former, is able to generate the strong compressive/extensional normal stresses which are thought to be the cause of the stationary asymmetry instability in planar cross-slots [4], despite the dilute character of our model fluids.
A stability map including all the fluids and flow rates discussed in this work is shown in Fig. 8.As discussed by Haward et al. [21,22], the nominal Elasticity number represents the trajectory through the Wi-Re parameter space for a given fluid.Representation of the Weissenberg number axis as (1 − β )Wi nom , where (1 − β ) is the polymer viscosity ratio, highlights an underlying pattern to the various time-dependent instability transitions, all of which occur in the range 0.33 < (1 − β )Wi nom < 1.06, which is substantially narrower than the corresponding Wi nom range, 1.2 < Wi nom < 6.4.

B. Distribution of normal stresses at the center
Figure 9 shows the probability distribution of the normalized first normal stress difference, ), at the center of the OSCER device, constructed from data cropped to an integer number of periods, for the flows used to describe the various elasticity regimes in the previous section.Each of these flows is compared to a steady-state case obtained at a lower flow rate.The evolution of the distribution of N * 1 with Wi nom provides information concerning the overall state of the flow field.For instance, low amplitude instabilities, such as those seen for the very low elasticity regime, lead to a narrow distribution of N * 1 , depicted in Fig. 9(a).Furthermore, whenever the birefringence strand assumes an asymmetric configuration, the values of N * 1 decrease concomitantly.For instance, as shown in Fig. 5(d to be smaller than the value of N * 1 for the formerly stationary flow, and the extent of this reduction varies considerably.For instance, in Fig. 9(d) there is an eightfold decrease in the average N * 1 at the onset of time-dependent instabilities.On the other hand, for example in Fig. 9(e), the onset of time-dependent flow is accompanied by a twofold decrease in the average value of N * 1 .Individually, both the degree of asymmetry and the amplitude of time-dependent fluctuations contribute to the decrease in the average of N * 1 .However, the superimposed effect of both types of instability is not trivial, and the interpretation of N * 1 distributions may require complementary information.Supplemental Material is provided in the form of space-time diagrams of N * 1 along the horizontal centerline of the OSCER geometry [33].Some of the weakly time-dependent flows may appear stationary upon inspection of these space-time plots, owing to the low amplitude of the normalized first normal stress difference fluctuations.
Concerning the very low elasticity regime, depicted in Fig. 9(a), besides the aforementioned narrow distribution of N * 1 , there is a slight decrease in the time-averaged value of N * 1 .The intermediate elasticity regime is depicted in Fig. 9(c).We group this regime with the very low elasticity flows due to various similarities, namely, the persistent integrity of the birefringence strand, which correlates with the low amplitude of velocity fluctuations.However, these velocity fluctuations are one order of magnitude lower than in the flows underlying Fig. 9(a), therefore the decrease of the modal N * 1 bin in Fig. 9(c) is an indicator of increased asymmetry, which nevertheless remains mild, as shown in Fig. 5(b).
The low elasticity regime is depicted in Fig. 9(b).Two primary local maxima are identifiable, the N * 1 bins [0, 40[ and [200, 240[, separated by a valley of low probability.The latter bin is part of a larger cluster which itself has two additional secondary maxima, [280, 320[ and [360, 400[.These four N * 1 bins correlate to different stages in the formation, rotation and breakup of the birefringence strand (cf.video S3 [33]).The [360, 400[ bin is the most probable range for N * 1 when the strand reforms after a breaking event, and as may be inferred from Fig. 3(a), there is a slight delay before the strand starts rotating.The [280, 320[ bin of N * 1 corresponds to a momentary arrest in the rotation of the strand, also visible in Fig. 3(a) when the width of the brighter, central region starts decreasing at a lower rate.The [200, 240[ bin is the most probable range for N * 1 when the strand reaches maximum tilt, and finally the [0, 40[ bin corresponds to the temporary absence of the strand immediately following a breaking event.The cycle deducible from the N * 1 distribution, namely the absence, de novo formation, rotation and maximum tilt preceding the momentary loss of the birefringence strand, shows the superposition of a limit cycle with an asymmetry.
The high elasticity regime is depicted in Fig. 9(d).This regime is notable since, in addition to the rotation of the birefringence strand, the position of the stagnation point also changes substantially in a continuous fashion, as shown in Fig. 6(a).Therefore, most of the time, the tilted strand does not intersect the center of the OSCER device, leading to a distribution of N * 1 that is dominated by a probability maximum in the lower range [0, 40[ and a second, significantly less probable local maximum in the bin [320, 360[, corresponding to the translation of the strand across the geometric center.The intermediate and high elasticity regimes illustrate how the superposition of two modes of instability can lead to different results depending on the corresponding amplitudes.Table I shows that the high amplitude limit cycle underlying Fig. 9(d) emerges with only small changes in the relaxation time and in the solvent viscosity ratio relative to Fig. 9(c).The high elasticity regime is also, by a wide margin, the regime wherein the accumulation of N * 1 is largest during stationary symmetric flow, followed by the intermediate regime, occupied by the single fluid aPS10-700 ppm, and closely by the fluid aPS16-350 ppm, shown in Fig. 9(f).Not coincidentally, these two fluids border the high elasticity regime in the El nom parameter space.The emergence of the birefringence strand fragmentation mechanism in the very high elasticity regime precludes the rotation of the strand, characteristic of various other preceding regimes, further illustrating the competition between different types of instability seen throughout this work.
The very high elasticity regime is depicted in Fig. 9(e).The N * 1 histogram resembles the characteristic shape of a Poisson distribution, with a maximum in the bin [40, 80[, followed by a series of trailing bins with progressively lower probability.After the fragmentation of the birefringence strand, the accumulation of normal stresses in the newly formed strand proceeds at an accelerating pace, as may be inferred from Fig. 7(c), leading to a decrease of the probability of the various progressively higher N * 1 bins (cf.video S9 [33]).Furthermore, the N * 1 bin [0, 40[ has lower probability than [40, 80[, suggesting that the strand does not completely disintegrate.

C. Spectral analysis of velocity fluctuations
Power density spectra of the normalized first component of the velocity vector, u x /U 0 , sampled at coordinates (x, y) = (H, 0), are depicted in Figs.10-14.Sampling at these coordinates instead of at the center, (x, y) = (0, 0), has the advantage of facilitating comparisons with experiments, since it is at the center of cross-slots that the noise-to-signal ratio in experimental measurements of the velocity is highest.The mean of each time series was subtracted from the signal prior to transformation, and no windowing of the signal was used other than cropping the sample to an integer number of periods, since most time series have a discernible dominant frequency.The power density is calculated as 2 , where FFT(f ) is the Fourier transform of the signal, f is the frequency, and M is the number of samples in the discrete time domain signal.In order to assess the possibility of aliasing caused by the finite cell size of computational meshes, we estimated the residence time in the sampled cell as t res = x cell /u x,rms .The characteristic local velocity is calculated as the root mean square of u x instead of the arithmetic mean because u x may often change sign, leading to unreasonably high values of t res .From the Sampling Theorem, a discrete sequence of samples can perfectly reconstruct a continuous signal if the sampling frequency is at least twice the highest frequency component of the signal.In other words, the cutoff frequency f c above which frequencies are poorly resolved due to the finite cell size of the mesh is estimated as 2f c ≈ 1/t res and is shown in dimensionless form in Table II for relevant flows.The estimated frequency cutoff, λ 0 f c , was borne in mind in the evaluation of the spectral decay of velocity fluctuations, as an upper limit on the fitted data range.
The very low elasticity regime, depicted in Fig. 10, is characterized by a principal frequency component, which is often the fundamental unit of a harmonic series and, specifically in Fig. 10(b), may be contaminated by a variable amount of weakly nonperiodic fluctuations.These fluctuations appear to increase in relative magnitude with the flow rate, up to a certain limit, above which a regularization of flow occurs.This regularization is most visible in Fig. 10(a) when the flow rate increases from Wi nom = 8.4 to Wi nom = 9.4 and in Fig. 10(b) a similar transition is seen between Wi nom = 7.6 and Wi nom = 8.4.A similar regularization is also seen for the intermediate elasticity regime, depicted in Fig. 11, when the flow rate increases from Wi nom = 5.8 to Wi nom = 6.4,although the latter spectra are primarily composed of broadband peaks instead of a well-defined harmonic series.The two regimes share other features, such as the trivial decay of the power spectra (data not shown) and more broadly the low amplitude of velocity oscillations and the persistent integrity of the birefringence strand.Qualitatively, emerging aperiodic components of motion appear to be suppressed by their periodic counterparts, thus the very low and intermediate regimes do not display the sudden instabilities seen in other regions of the inertia-elasticity-viscosity parameter space.
The low elasticity regime is depicted in Fig. 12.The overall periodicity of the flow is visible in the initial portion of the spectra, such as the peak at Wi nom = 3.4.However, the relative power density of the aperiodic modes of oscillation is now much larger, and the power law decay range has an exponent, α, estimated as the slope in log-log coordinates, of approximately α ≈ −4.9, on average.The fluid aPS10-350 ppm is the only instance in this work where the length of the time domain signal was deemed insufficient for the calculation of power spectra, essentially due to the inconsistent amplitude and, to a lesser extent, periodicity of u x .Nevertheless, the inconsistent amplitude of the cycles underlying the aPS10-350 ppm flows effectively serves as a secondary indicator of greater instability.
The high elasticity regime is depicted in Fig. 13.The low frequency range is dominated by a harmonic series comprising the odd multiples of the fundamental frequency.This frequency composition is characteristic of certain regular waveforms, such as the square wave or triangle wave, and indeed the limit cycle in Fig. 6(d) confirms the well-defined periodicity of this regime.The decay range has a power law exponent of about α ≈ −4.3.In the very high elasticity regime, depicted in Fig. 14, as the flow rate increases, the formerly predominantly harmonic spectra become increasingly irregular in the low frequency range, without significant changes in the power law decay rate, α ≈ −5.9.Note that there are no harmonic series in Fig. 14(a) possibly because we increased the flow rate in steps corresponding to Wi nom = 0.6 and may have skipped this regime.Concerning the specific spectral decay estimates for each flow, given in Table II, we underline that, even though the values shown are accurate for the specific fit range and method, some uncertainty exists on the correct boundaries of the decay range.
For elastic turbulence, the velocity power spectrum decay rate is expected to be faster than α = −3 for unbounded, isotropic, and homogeneous flows [35], and early experimental observations point to a decay rate between α = −3 and α = −3.5 [36,37].Given the lack of a universal definition of turbulence, we resort to one parameter to quantify the phenomenon, the turbulence intensity, T, calculated along the downstream centerline of the OSCER device, here defined as where the operators • and • denote, respectively, the magnitude and the average over time.Note that periodic flows may also generate a high value of T, therefore a separate method is necessary to evaluate if a given flow qualifies as chaotic.From the preceding spectral analysis, the very low and intermediate elasticity regimes are characterized by regular, low amplitude fluctuations and lack a broad range of excited frequencies, one of the commonly cited characteristics of turbulence.On the other hand, the high elasticity regime is unequivocally periodic, and in particular is defined by a well-resolved harmonic series.For the low and very high elasticity regimes, only the highest flow rates seem to qualify as chaotic, mildly for the former regime and strongly for the latter.The corresponding turbulence intensity plots are shown in Fig. 15.As expected, the results for each FIG. 13.Power spectra obtained from time series of the normalized first velocity component u x /U 0 sampled at coordinates (x, y) = (H, 0), for the high elasticity regime.The degree of vertical shifting of the PSD is indicated in the legend.In panel (a.1) the frequency axis is linear and the power density axis is logarithmic, whereas in panel (a.2) both axes are logarithmic.TABLE II.Summary of findings concerning the power spectra of velocity fluctuations in the OSCER device, for the first normalized velocity component u x /U 0 sampled at coordinates (x, y) = (H, 0).The power law decay α was estimated via a linear fit, log 10 (PSD) = αlog 10 (λ 0 f ) + c, of the spectral data within the range [λ 0 f low , λ 0 f high ], and ±ε α is the corresponding 95% confidence interval, rounded up to the nearest decimal place.λ 0 f c is the estimated cutoff above which frequency components of the signal cannot be accurately resolved due to the finite cell size of computational meshes.The values of α ± ε α are accurate for the specific fit range and method; however, some uncertainty exists concerning the exact boundaries of the decay range.regime differ, with the low elasticity regime having mild values of T near the center of the OSCER device, and the very high elasticity regime exhibiting 2.4-fold higher values of T over the same region.However, for all examples, T decays to a residual value near the outlet of the geometry, indicating an ephemeral character of the unstable flow near the center, induced by the high normal stresses therein.An additional observation concerns the relation between the principal frequency and the characteristic flow rate, depicted in dimensionless form in Fig. 16.It is unsurprising that fluids with higher elasticity develop instabilities at lower flow rates and consequently with lower thresholds for excitable frequencies, as shown in Fig. 16(a).We define the nominal effective relaxation time, λ eff,nom , which is the theoretical effective relaxation time, λ eff = λ 0 /g, with g given by Eq. ( 6), calculated at the characteristic deformation rate ε = 1.5U 0 /(15H ) of the corresponding simulation.Multiplication of the principal frequency f max by λ eff,nom yields a narrow distribution of values, shown in Fig. 16(b), with λ eff,nom f max = 0.075.This value is a reasonable approximation of the principal frequency component nearly up to the onset of chaotic flow, above which the identification of a characteristic frequency is no longer possible and a power law spectral decay over a wide range of frequencies is observed.

IV. DISCUSSION
Central to the phenomenology described in this work is the idea that different types of instability can coexist in a superimposed state.Moreover, the two primary types of transition we discuss, periodic time-dependent and stationary asymmetric, often compete in order to determine the dynamics of the flow.This behavior is common in many other nonequilibrium systems [38,39].In the field of fluid dynamics, two recent studies are of particular interest to this discussion.Wen et al. [40] studied the flow of a shear-thinning, 0.15% (w/w) aqueous solution of xanthan gum in a cylindrical pipe.The authors measured streamwise velocity profiles for Reynolds numbers corresponding to laminar, transition, and turbulent flows.In the transition regime a loss of symmetry occurs and the location of peak velocity shifts away from the center of the pipe.Symmetry is then restored in the turbulent regime.This flow rate dependent restoration of symmetry is also seen at least partially in several of the flows discussed throughout this work, most notably concerning fluid aPS16-350 ppm in the very high elasticity regime.  of symmetry occurs.In Ref. [40], the restoration of symmetry is shown to be a consequence of the onset of turbulence.By analogy, we posit that the restoration of symmetry induced by the onset of time-dependent instabilities for the fluid aPS16-350 ppm constitutes indirect evidence of the onset of the chaotic flow regime.
Varshney and Steinberg [41] studied the flow of a 100 ppm, 18 MDa polyacrylamide aqueous solution in 62% (w/w) sucrose and 1% (w/w) NaCl, past two widely spaced cylindrical obstacles, with normalized spacing E = 3.3.The creeping viscoelastic flow generates a sequence of two elastic instabilities, a loss of time-reversal symmetry and at higher flow rates a superimposed loss of mirror symmetry.The authors then speculate on why the second instability was not reported in previous studies, such as in Refs.[42,43].Below a critical normalized cylinder spacing E = 1.25 only the loss of time-reversal symmetry is present, leading to the conclusion that a codimension-2 point probably occurs in the range 1.25 < E < 3.3.Previous experimental [14] and computational [15] work on creeping viscoelastic flows in standard planar cross-slots with variable aspect ratio indicates a reverse sequence of transitions.An initial loss of symmetry is followed by the onset of time-dependent flow at higher flow rates.However, for sufficiently shallow channels, the flow transitions directly to an unsteady state.Based on our results, it is possible that a superposition of instabilities also occurs in standard planar cross-slots.

V. CONCLUSIONS
on simulations of the flow of Boger fluids in the OSCER device, we identified a family of predominantly elastic, time-dependent instabilities.The different types of flow dynamics may be grouped by their nominal Elasticity number.For El nom = {1.6,2.3}, the flow is characterized by low amplitude, periodic oscillations, with a high degree of flow field symmetry and no loss of birefringence strand integrity.For El nom = {3.2,4.1}, the flow field gradually loses symmetry over time as the birefringence strand rotates either clockwise or counterclockwise, eventually resulting in the breaking and subsequent reconstitution of the strand.For El nom = 5.6, the flow FIG. 1. (a)Geometry of the optimized shape cross-slot extensional rheometer (OSCER device).The halfwidth H is used as the characteristic length scale.(b) Detail of the computational mesh.The deformation of computational cells is the result of the application of an optimization algorithm to the mesh of a planar cross-slot, described in Ref.[18].(c) Zoomed-out view of the mesh, showing the inlet and outlet channels.Cells in these channels are progressively larger in the streamwise direction towards the inlets and outlets.

FIG. 2 .
FIG. 2. Illustrative results for the very low elasticity regime.The example given is the fluid aPS7-350 ppm at Wi nom = 10.4 and El nom = 1.6.The initial, transient segment of the simulation was discarded, and time was reset to zero.(a), (b) Space-time diagrams of the normalized first normal stress difference, respectively, along the horizontal and vertical centerlines, with the evolving position of the stagnation point overlaid on the former.(c) Continuous line: time series of the first velocity component u x at the center of the OSCER device, normalized by the characteristic velocity U 0 ; dotted line: time series of the normalized first normal stress difference sampled at the same location.(d) Phase-plane plot constructed with the time series shown in the previous panel.

3 . 6 Figure 4
Figure 4 illustrates typical flow behavior in the intermediate elasticity regime for Wi nom = 7.6.Oscillations in velocity and normal stresses at the center of the OSCER device are quasiperiodic, albeit with very low amplitude, such that the normalized velocity is mostly contained within the range u x /U 0 ≈ ±0.005.Therefore the position of the stagnation point is nearly constant, and consequently the aforementioned oscillations may be thought of as residual.However, as shown

FIG. 3 .
FIG. 3. Illustrative results for the low elasticity regime.The example given is the fluid aPS7-1400 ppm at Wi nom = 3.4 and El nom = 3.2.The initial, transient segment of the simulation was discarded, and time was reset to zero.(a), (b) Space-time diagrams of the normalized first normal stress difference, respectively, along the horizontal and vertical centerlines, with the evolving position of the stagnation point overlaid on the former.(c) Continuous line: time series of the first velocity component u x at the center of OSCER device, normalized by the characteristic velocity U 0 ; dotted line: time series of the normalized first normal stress difference sampled at the same location.(d) Phase-plane plot constructed with the time series shown in the previous panel.

FIG. 4 .
FIG. 4. Illustrative results for the intermediate elasticity regime.The example given is the fluid aPS10-700 ppm at Wi nom = 7.6 and El nom = 5.6.The initial, transient segment of the simulation was discarded, and time was reset to zero.(a), (b) Space-time diagrams of the normalized first normal stress difference, respectively, along the horizontal and vertical centerlines, with the evolving position of the stagnation point overlaid on the former.(c) Continuous line: time series of the first velocity component u x at the center of OSCER device, normalized by the characteristic velocity U 0 ; dotted line: time series of the normalized first normal stress difference sampled at the same location.(d) Phase-plane plot constructed with the time series shown in the previous panel.

Figure 7
Figure7illustrates the typical behavior near the higher end of the elasticity spectrum assessed in this work, at Wi nom = 2.0.The flow may still be described as periodic, albeit with significant

FIG. 7 .
FIG. 7. Illustrative results for the very high elasticity regime.The example given is the fluid aPS16-1400 ppm at Wi nom = 2.0 and El nom = 32.2.The initial, transient segment of the simulation was discarded, and time was reset to zero.(a), (b) Space-time diagrams of the normalized first normal stress difference, respectively, along the horizontal and vertical centerlines, with the evolving position of the stagnation point overlaid on the former.(c) Continuous line: time series of the first velocity component u x at the center of OSCER device, normalized by the characteristic velocity U 0 ; dotted line: time series of the normalized first normal stress difference sampled at the same location.(d) Phase-plane plot constructed with the time series shown in the previous panel.

FIG. 8 .
FIG. 8. Stability diagram showing the evolution of each fluid in the Wi-Re parameter space.Empty circles represent stationary symmetric flow, whereas filled squares indicate a stationary asymmetric field configuration.Each of the five time-dependent regimes identified in this is represented as follows: filled circles, very low elasticity regime; filled diamonds, low elasticity regime; filled triangles, intermediate elasticity regime; empty triangles, high elasticity regime; empty diamonds, very high elasticity regime.

13 FIG. 9 .
Figure9shows the probability distribution of the normalized first normal stress difference, N * 1 = N 1 /(η 0 U 0 /2H ), at the center of the OSCER device, constructed from data cropped to an integer number of periods, for the flows used to describe the various elasticity regimes in the previous section.Each of these flows is compared to a steady-state case obtained at a lower flow rate.The evolution of the distribution of N * 1 with Wi nom provides information concerning the overall state of the flow field.For instance, low amplitude instabilities, such as those seen for the very low elasticity regime, lead to a narrow distribution of N * 1 , depicted in Fig.9(a).Furthermore, whenever the birefringence strand assumes an asymmetric configuration, the values of N * 1 decrease concomitantly.For instance, as shown in Fig.5(d), the fluid aPS16-350 ppm develops a stationary asymmetry at Wi nom = 3.4, leading to a threefold decrease in the N * 1 values due to this asymmetry, as shown in Fig.9(f).Another important observation concerns the relation between N * 1 and time-dependent instabilities.At the onset of time-dependent flow, the average of the newly formed distribution of N * 1 tends 113301-13

FIG. 11 .
FIG. 11.Power spectra obtained from time series of the normalized first velocity component u x /U 0 sampled at coordinates (x, y) = (H, 0), for the intermediate elasticity regime.The degree of vertical shifting of the PSD is indicated in the legend.The frequency axis is linear, and the power density axis is logarithmic.

FIG. 12 .
FIG. 12. Power spectra obtained from time series of the normalized first velocity component u x /U 0 sampled at coordinates (x, y) = (H, 0), for the low elasticity regime.The degree of vertical shifting of the PSD is indicated in the legend.Both the frequency and power density axes are logarithmic.

FIG. 15 .
FIG. 15.Turbulence intensity along the downstream centerline of the OSCER device for the highest simulated flow rate for each of the test fluids.(a) Low elasticity regime and (b) very high elasticity regime.A combined linear/log scale is used in the x/H axis.
Figure 5(d)  and the accompanying changes in the distribution of normalized N 1 values at the center of the OSCER device shown in Fig.9(f) demonstrate a loss of symmetry.With the onset of time-dependent instabilities characteristic of the very high elasticity regime, the birefringence strand is no longer able to rotate and thus a restoration

FIG. 16 .
FIG.16.Principal frequency f max (a) multiplied by the zero shear-rate relaxation time λ 0 or (b) multiplied by the nominal effective relaxation time λ eff,nom , plotted against the dimensionless flow rate expressed as the nominal Weissenberg number Wi nom .λ eff,nom is calculated at the center of the OSCER device for a hypothetical stationary extensional flow with characteristic strain rate ε = 1.5U 0 /15H .Only the flows with a clearly defined maximum power density peak are represented, usually in the form of the principal component of a harmonic series.

TABLE I .
[22]erties of the model fluids.Solvent shear viscosity for all fluids is η s = 0.059 Pa s, and density is ρ = 981 kg m −3 .The nominal Elasticity number El nom is calculated using zero shear-rate properties.Additional data were kindly provided by the authors of Ref.[22].