Self-sustained elastoinertial Tollmien-Schlichting waves

Direct simulations of two-dimensional plane channel flow of a viscoelastic fluid at Reynolds number Re = 3000 reveal the existence of a family of attractors whose structure closely resembles the linear Tollmien-Schlichting (TS) mode, and in particular exhibits strongly localized stress fluctuations at the critical layer position of the TS mode. At the parameter values chosen, this solution branch is not connected to the nonlinear TS solution branch found for Newtonian flow, and thus represents a new solution family that is nonlinearly self-sustained by viscoelasticity. The ratio between stress and velocity fluctuations is in quantitative agreement for the attractor and the linear TS mode, and increases strongly with Weissenberg number, Wi. For the latter, there is a transition in the scaling of this ratio as Wi increases, and the Wi at which the nonlinear solution family comes into existence is just above this transition. Finally, evidence indicates that this branch is connected through an unstable solution branch to two-dimensional elastoinertial turbulence (EIT). These results suggest that, in the parameter range considered here, the bypass transition leading to EIT is mediated by nonlinear amplification and self-sustenance of perturbations that excite the Tollmien-Schlichting mode.


Introduction
Adding minute quantities (parts per million) of long-chain polymer additives can dramatically change the turbulent flow of Newtonian fluids, the most significant consequence being the considerable drop in friction factor, which is commonly referred to as the Toms effect (Virk 1975;White & Mungal 2008;Graham 2014). Accompanying this macroscopic change is a structural change to the flow. At high levels of viscoelasticity, Samanta et al. (2013) and Dubief et al. (2013) have shown that trains of weak spanwise-oriented flow structures with inclined sheets of polymer stretch dominate the flow, denoting this regime as elastoinertial turbulence (EIT). In further contrast to the 3D structures that sustain Newtonian turbulence, Sid et al. (2018) demonstrated that EIT is fundamentally 2D in nature by showing that the structures sustaining 2D EIT in channel flow simulations are similar to those in 3D. Choueiri et al. (2018) experimentally studied the path to EIT in pipe flow by varying polymer concentration at fixed Reynolds number, Re. For sufficiently low Re, they observed an initial drop in friction factor (i.e. modification of Newtonian turbulence) as concentration increased, followed by re-laminarization and eventually by a reentrant transition to EIT, where the flow had very different structure from Newtonian turbulence. These observations point at two distinct types of turbulence in dilute polymer solutionsone that is suppressed by viscoelasticity (Newtonian turbulence) and one that is promoted (EIT). Shekar et al. (2019) corroborated this observation of a reentrant transition to EIT in simulations of channel flow with increasing Weissenberg number, Wi, the ratio between the polymer relaxation time scale and the shear time scale. They further showed that close to its inception, EIT exhibits localized polymer stress fluctuations that bear strong resemblance to critical layer structures predicted by linear analyses, i.e., sheetlike fluctuations localized at wall-normal locations where the disturbance wavespeed equals the base flow velocity. In particular, they demonstrated that the fluctuation structure corresponding to the dominant spectral content strongly resembles the viscoelastic extension of the linear Tollmien-Schlichting wave (TSW). This is perhaps a surprising result, as the flow in the parameter regime considered is linearly stable, and in Newtonian turbulence, the TS mode plays a very limited role. Some light is shed on this issue through resolvent analysis, i.e., determination of the response of the linearized dynamics to harmonic-intime disturbances, which shows that the linear TS mode becomes highly amplified in the presence of viscoelasticity. This strong amplification implies that even very weak disturbances may be sufficient to trigger the nonlinear effects necessary to sustain EIT.
We note that similar structures have been observed by other researchers in different contexts. Page & Zaki (2015) analyzed the evolution of vortical perturbations in 2D viscoelastic simple shear flow. Their analysis reveals a viscoelastic analogue of the Newtonian Orr mechanism. This "reverse-Orr" mechanism generates tilted sheets of polymer stress fluctuations resembling those seen at EIT and thus may play some role in this phenomenon.
Because prior work on elastoinertial turbulence reveals structures similar to those seen in the linear Tollmien-Schliching mode, the present work focuses on Tollmien-Schlichting waves, but in the fully nonlinear context of self-sustained solutions in the channel flow geometry. (In the parameter regime here, the laminar flow is always linearly stable.) After introducing the formulation and computational methods, we show how the Newtonian nonlinear TS wave branch is modified by viscoelasticity, resulting in its disappearance as Wi increases. At still higher Wi, however, we demonstrate the onset of a new, viscoelasticity-driven, nonlinear solution branch that strongly resembles the linear Tollmien-Schlichting mode, and illustrate how it is related to the TS mode of linear theory and to elastoinertial turbulence.

Formulation
This study focuses on two-dimensional pressure-driven channel flow with constant mass flux. The x and y axes are aligned with the streamwise and wall-normal directions, respectively. Lengths are scaled by the half channel height l, so the dimensionless channel height L y = 2. The domain is periodic in x with length L x . Velocity v is scaled with the Newtonian laminar centerline velocity U ; time t with l/U , and pressure p with ρU 2 , where ρ is the fluid density. The polymer stress tensor τ p is related to the polymer conformation tensor α through the FENE-P constitutive relation, which models each polymer molecule as a pair of beads connected by a nonlinear spring with maximum extensibility b.
We solve the momentum, continuity and FENE-P equations: Here Re = ρU l/(η s + η p ), where η s and η p are the solvent and polymer contributions to the zero-shear rate viscosity. The viscosity ratio β = η s /(η s + η p ). We fix β = 0.97 and b = 6400. Since 1 − β is proportional to polymer concentration and b to the number of monomer units, this parameter set corresponds to a dilute solution of a high molecular weight polymer. The Weissenberg number Wi = λU/l, where λ is the polymer relaxation time.
For the nonlinear direct numerical simulations (DNS) described below, a finite difference scheme and a fractional time step method are adopted for integrating the Navier-Stokes equation. Second-order Adams-Bashforth and Crank-Nicolson methods are used for convection and diffusion terms, respectively. The FENE-P equation is discretized using a high resolution central difference scheme (Kurganov & Tadmor 2000;Vaithianathan et al. 2006;Dallas et al. 2010). No artificial diffusion is applied. Resolution tests were performed to ensure convergence of statistics. A typical resolution for the following results is (N x , N y ) = (79, 402).
We also consider the linearized evolution of infinitesimal perturbations to the laminar state with given streamwise wavenumber k. Two approaches are used. The first is classical linear stability analysis, in which solutions of the formφ(y) exp [ik(x − ct)] are sought, resulting in an eigenvalue problem for the complex wavespeed c at a given k. In the present study,ˆalways indicates deviation from the laminar state. If all c i < 0, which is the case for all conditions considered in the present study, the flow is linearly stable. A linearized version of the DNS code was also developed using the numerical schemes described above. Results were validated against linear stability analysis and agreement to three significant digits was obtained for the value of c for the viscoelastic TS mode at the parameters of interest.
The second linear approach used here determines the linear response of the laminar flow to external forcing with given wavenumber k and frequency ω using the resolvent operator of the linearized equations (Schmid 2007;McKeon & Sharma 2010). The norm used in the resolvent calculations is where A is the conformation tensor in the laminar state. The second term provides a measure of the conformation tensor perturbation magnitude that is motivated by the non-Euclidean geometry of positive-definite tensors (Hameduddin et al. 2019). For both the linear stability and resolvent analyses, the equations are discretized with a Chebyshev pseudospectral method using 401 Chebyshev polynomials.

Origin of Newtonian and viscoelastic nonlinear Tollmien-Schlichting attractors
In Newtonian flow, a family of nonlinear Tollmien-Schlichting waves bifurcates subcritically from the laminar branch at Re ≈ 5772 with L x ≈ 2π/1.02 ≈ 6.15. The lower limit of the parameter regime for which this solution family exists is Re ≈ 2800, L x ≈ 2π/1.3 ≈ 4.83 (Jiménez 1990). Furthermore, in prior work on elastoinertial turbulence (Shekar et al. 2019), as well more recent simulations in long two-dimensional domains, a strong peak in the spatial spectrum is found at L x ≈ 5. Based on these observations, all of the results presented in this study will be at Re = 3000, L x = 5. In the Newtonian limit at these parameters, there are upper and lower branch solutions (which merge in a saddle-node bifurcation as Re is lowered); the upper branch traveling wave solution is linearly stable with respect to two-dimensional perturbations and is thus easily computed via DNS. We call this solution branch, including its viscoelastic extension, the Newtonian Nonlinear Tollmien-Schlichting attractor (NNTSA). (The word "attractor" is chosen rather than "wave" because, depending on parameters one can observe a pure traveling wave state or one with periodic or nonperiodic modulations.) On increasing W i, the self-sustained nonlinear viscoelastic TS wave at Re = 3000 develops sheets of high polymer stretch resembling near wall structures seen at EIT. This structure originates in the nonlinear Kelvin cat's eye kinematics of TS waves at finite amplitude, as detailed in Shekar et al. (2019). Some of the observations made in that paper are included here, as they form the background for the new main results of the present study. Figure 1(a) illustrates this point with a snapshot ofα xx on the NNTSA branch at W i = 3.
At the parameters chosen, the solution branch originating in the self-sustained Newtonian TS wave bifurcates to a periodic orbit at Wi ≈ 3.5 (cf. Lee & Zaki (2017)) before turning back into a traveling wave and losing existence beyond Wi = 3.875, evidently in a saddle-node bifurcation yielding a lower branch TS wave solution that becomes the Newtonian solution as Wi → 0. Consistent with a saddle-node bifurcation, if the solution at Wi = 3.875 is used as an initial condition for a simulation at slightly higher Wi, the flow laminarizes. This bifurcation scenario is shown on Fig. 2 in terms of average wall shear rate vs. Wi. Because we have not directly computed the lower branch, it is indicated as a dashed curve.
As shown in Shekar et al. (2019), if Wi is large, sufficiently energetic initial conditions lead to 2D EIT. For comparison with the structure on the NNTSA branch, Figure 1(b) shows a snapshot ofα xx on the EIT branch at Wi = 15. Fig. 2(a) also shows the mean wall shear rate for the EIT solution branch, which loses existence at finite amplitude when

Wi
13. The bifurcation underlying this transition is presumably also of saddle-node form.
The central observation of the present paper arises from considering what happens just below the onset of the EIT regime at Wi ≈ 13. We do this by using a velocity and stress field from EIT at Wi = 13 as an initial condition for a run at Wi = 12. This initial condition persists as a slowly decaying form of EIT for hundreds of time units, consistent with behavior just beyond a saddle-node bifurcation. As time increases further, the structure continues to decay, but does not ultimately reach the laminar state. Instead, it evolves to a nontrivial attractor state that is very nearly a traveling wave, and in particular strongly resembles the linear TS mode at these parameters. We call this new state the viscoelastic nonlinear Tollmien-Schlichting attractor (VNTSA).

Viscoelastic nonlinear Tollmien-Schlichting attractor
To elaborate on the relationship between the VNTSA and the linear TS mode, we now describe 2D simulation results at Re = 3000, Wi = 13, L x = 5, i.e. close to the point where the 2D EIT branch first comes into existence as shown in the bifurcation diagram (Figure 2a). EIT and the VNTSA are coexisting attractors at these parameter values. Figure 3 shows the evolution of the L 2 norm ofα xx starting from an initial condition consisting of the laminar state plus some amplitude of the linear TS mode for this parameter set. This mode, with = 1, is shown in Figure 4a. The structure of the velocity field is virtually unchanged from the Newtonian case and the polymer conformations are strongly localized to the critical layer positions y = ±0.82. Sufficiently small perturbations, e.g. = 1, decay to the laminar state, as they must since that state is linearly stable. However, larger perturbations = 10 and = 100, where nonlinear mechanisms play a role,α xx settle to a finite value corresponding to the VNTSA.
For comparison, the dashed lines on Figure 3 show the linearized evolution starting from the same initial conditions; these all decay to laminar, illustrating the role of nonlinearity in the transition to the VNTSA. This state is robust: initial perturbation amplitudes over a wide range will evolve to it. However, initial conditions with very large Figure 3: Time evolution of the L 2 -norm ofα xx for Re = 3000, Wi = 13, L x = 5, starting from an initial condition of laminar state + × TS-mode. Dashed lines correspond to linearized runs starting from the same initial conditions for = 10 and 100.
magnitudes (e.g. = 6000) evolve to EIT: as noted above, both EIT and VNTSA are attractors at the chosen parameters (as is the laminar state). Figure 4b is a snapshot showing the typical fluctuation structure of the VNTSA at Wi = 13. The streamwise conformationα xx has tilted sheets highly localized near y = ±0.82 and contours of wall normal velocityv span the entire channel. This structure bears strong resemblance to the TS mode shown in Figure 4a. The VNTSA is thus a weakly nonlinear self-sustaining state whose primary structure is the viscoelastic TS mode. We elaborate in the following section on the linear TS mode and its connections to the VNTSA.
In the VNTSA state, the velocity fluctuations are very weak, and the mean wall shear rate displays a very small change from laminar. This can be understood on the grounds that changes of the mean wall shear rate correspond to fluctuations with k x = 0, which arise only due to nonlinear interactions. Since the primary velocity structure is very weak, the nonlinear effects will be even weaker. To illustrate nonlinear effects, Figures 4c and  4d, respectively, show the kL x /2π = 1, 2, spatial Fourier components of the snapshot shown in 4b. Figure 4c closely resembles the TS mode, with a slight symmetry-breaking across the centerline y = 0. The structure at kL x /2π = 2 also displays polymer stress fluctuations localized around the critical layer position, an observation that also holds for higher wavenumbers.
Having established the structure of the flow on the VNTSA branch, we now illustrate the bifurcation scenario of this solution branch by continuing in Wi. The VNTSA branch loses existence at finite amplitude (i.e. in a saddle-node bifurcation) for Wi 6, as we have confirmed both by using the Wi = 6 solution as an initial condition for simulations at lower Wi and by running simulations starting from the laminar state perturbed by the TS mode with small . For W i < 6 all these initial conditions decay to laminar. On increasing Wi, the VNTSA branch loses existence beyond Wi ≈ 49, and initial conditions that land on the VNTSA for Wi = 49 evolve to EIT at Wi = 50. These observations suggest that the VNTSA turns around and forms an unstable branch that joins up with the unstable lower branch of EIT. Due to the small amplitude of the VNTSA branch, the bifurcation scenario associated with it is shown separately in Figure 2b, using the L 2 norm ofα xx as the amplitude measure. The hypothesized unstable branch connecting VNTSA and EIT is shown schematically with dashed lines on the bifurcation diagrams and the asterisks indicate how they are presumably connected in moving from Figure 2a to Figure 2b. Figure 5a shows the fluctuation structure of the VNTSA at Wi = 8, close to the point where it first comes into existence. The structure closely resembles the TS mode and does not change appreciably with time. The flow is almost a pure nonlinear traveling wave with some weak non-periodicity, as indicated in the power density plot of the wall normal velocity at position (0, 0.825) shown in figure 5b. The spectrum is mainly composed of the dominant TS mode frequency and its higher harmonics. The dynamics and structures get more complicated as Wi increases. Figure 5c shows a typical snapshot at Wi = 20, which clearly is more complex than a TS mode. However, at this Wi, the VNTSA still intermittently displays clear TS-like structures such as the snapshot in Figure 5d.

Linear analyses
In this section we elaborate on the linearized problem and its connection to the attractors described above using linear stability and resolvent analyses. The spectrum corresponding to disturbances with wavelengths equal to the DNS box size, i.e. k = 2π/5, has a least-stable eigenvalue at c ≈ 0.32 − 0.010i, and the associated eigenfunction is the viscoelastic extension of the TS mode. For low values of Wi, the mode is less stable than its Newtonian counterpart, while for Wi 2, it becomes more stable with increasing elasticity; this non-monotonic behavior has been reported by Zhang et al. (2013), who attribute it to viscoelastic modification of the phase difference between u and v. However, over the range of Wi considered here, the eigenvalue varies by less than 1% of the Newtonian value. Importantly, the mode never becomes unstable in the range of parameters considered here, confirming the observations in Section 3.2 that finite amplitude disturbances are required to trigger transition to EIT or the VNTSA. However, linear instabilities not related to the TS mode have been found in other regions of parameter space (Garg et al. 2018), implying the possibility of different attractor families in those regions.
A measure of the relative importance of the conformation tensor and velocity disturbances is the ratio of the peak amplitudes ofα xx , (the largest component of the conformation tensor), andv. This ratio is shown in Figure 6a. Two distinct regimes are apparent, with the transition between the two occurring at Wi ≈ 3.1. The low Wi regime scales as Wi 2 , which is the same scaling as in linear shear flow. The amplitude ratio above the change in slope does not exhibit power law scaling. The change in slope at Wi ≈ 3.1, can be understood by examining the α xx mode shapes, the magnitudes of which are plotted in Figure 6b for several values of Wi in the range shown in Figure 6a. For small Wi, the disturbance is largest at the wall and decays rapidly away from it. Therefore, the Wi 2 scaling in this regime can be explained by the fact that the leadingorder approximation of the base flow very near the wall is simple shear. As Wi increases, this value decreases, while a new local maximum emerges and grows, becoming the global maximum just above Wi = 3; the arrow in the figure indicates the profile where this occurs. Upon further increase in Wi, the maximum gradually shifts away from the wall, and the modes become increasingly localized around the location of the critical layer y c , at which the real part of the wavespeed equals the base flow velocity. The critical layer for Wi = 13 is indicated by the vertical dashed line. This suggests that a critical layer mechanism is responsible for the change in scaling at large Wi, though at present we do not understand the specific origin of this result. Interestingly, the Wi at which the VNTSA comes into existence is only slightly larger than that at which the transition to critical layer scaling occurs.
Also shown in Figure 6a is the amplitude ratio computed from the VNTSA for several values of Wi. Excellent agreement between the linear and nonlinear results quantitatively reinforces the TS-mode-like nature of the VNTSA. Additionally, the profile of |α xx |, averaged in the streamwise direction and over many snapshots, for the VNTSA at Wi = 13 is shown by the thick red line in Figure 6b, and the blue line highlights the linear mode for the same Wi. The VNTSA profile exhibits the same localization, and the location of the peak value is in close agreement with the critical layer location. Figure 6c shows the first two singular values of the resolvent operator for k and c corresponding to the linear TS mode. Shekar et al. (2019) showed that such modes are the most-amplified 2D disturbances in this parameter regime and that the leading response mode is nearly identical to the TS eigenmode; for this reason the resolvent modes are not plotted separately. The substantial increase in the leading singular value with Wi indicates that this amplification becomes much stronger with increasing elasticity, and consequently that considerably smaller disturbances may be sufficient to trigger selfsustaining nonlinear mechanisms. Further, the symmetry of the flow geometry about y = 0 means that resolvent modes typically come in pairs having similar amplification, with one mode having a symmetricv response and the other having an antisymmetricv response. However, the growing separation between the first and second singular values with increasing Wi indicates that this pairing is broken by elasticity, and that the symmetry exhibited by the TS mode is preferred in terms of linear amplification.

Conclusion
This study focuses on two-dimensional plane channel flow of a very dilute polymer solution at Re = 3000. At sufficiently high Wi, elastoinertial turbulence is observed in this parameter regime, and the focus of the present work is to make progress toward understanding the structures and mechanisms underlying the dynamics in this regime. We report here the existence of a new attractor that is based on the viscoelastic linear Tollmien-Schlichting mode and is nonlinearly sustained by viscoelastic stresses. We denote this as the viscoelastic nonlinear Tollmien-Schlichting attractor (VNTSA). At the parameters considered here, this solution branch is not connected to the Newtonian branch of nonlinear self-sustained Tollmien-Schlichting waves; it would be interesting to learn whether they become connected at higher Re. In a domain of dimensionless length 5, this solution comes into existence at finite but very small amplitude when Wi 6, increasing in amplitude until Wi ≈ 49 where it loses existence again. At higher Wi, initial conditions corresponding to this solution branch at lower Wi evolve into elastoinertial turbulence. In general, we find not pure nonlinear traveling waves, but, until Wi is large, the nonperiodic fluctuations are very small. The connection of the VNTSA to the linear TS mode is established via their strong structural similarities, including a quantitive agreement between the relative magnitudes of the velocity and stress fluctuations. The value of Wi at which the VNTSA comes into existence is close to where the relative amplitude of the stress and velocity fluctuations for the linear TS mode undergoes a change in scaling. Above this transition the stress fluctuations become highly localized at the position of the critical layer.
Taken together, these results suggest that, at least in the parameter range considered here, the bypass transition leading to EIT is mediated by nonlinear amplification and self-sustenance of perturbations that excite the Tollmien-Schlichting mode. Gaining an understanding of the mechanism underlying this phenomenon will shed light on the origin of elastoinertial turbulence.