Purely elastic linear instabilities in parallel shear flows with free-slip boundary conditions

We perform a linear stability analysis of viscoelastic plane Couette and plane Poiseuille flows with free-slip boundary conditions. The fluid is described by the Oldroyd-B constitutive model, and the flows are driven by a suitable body force. We find that both types of flow become linearly unstable, and we characterise the spatial structure of the unstable modes. By performing a boundary condition homotopy from the free-slip to no-slip boundaries, we demonstrate that the unstable modes are directly related to the least stable modes of the no-slip problem, destabilised under the free-slip situation. We discuss how our observations can be used to study recently discovered purely elastic turbulence in parallel shear flows.


Introduction
Stability of parallel shear flows of dilute polymer solutions is an outstanding problem of mechanics of complex fluids. While linear stability analyses performed with model viscoelastic fluids suggest that these flows are linearly stable (Gorodtsov & Leonov 1967;Wilson et al. 1999), it has been proposed that they can lose their stability through a finite-amplitude, sub-critical bifurcation (Morozov & van Saarloos 2005, similar to their Newtonian counterparts (Grossmann 2000). Recent experiments performed in pressure-driven straight channels support this scenario: For sufficiently high flow rates, there exists a critical strength of flow perturbations required to destabilise the flow (Pan et al. 2013), thus demonstrating the sub-critical nature of the transition. Beyond the transition, the flow exhibits features of purely elastic turbulence (Bonn et al. 2011;Pan et al. 2013;Qin & Arratia 2017), previously only reported in geometries with curved streamlines (Groisman & Steinberg 2000;Steinberg 2021). † Email address for correspondence: martin.lellep@ed.ac.uk ‡ Deceased on the 7 th of August 2019. Professor Eckhardt took active part in conceiving this project and interpreting the free-slip pCF results. He has not seen the pPF part of this work that we hope we continued in line with his high scientific standards. All shortcomings should be attributed to A.M.
Despite recent progress, the exact origin of purely elastic turbulence remains unknown. One of the major obstacles in understanding its mechanism is the absence of direct numerical simulations of confined flows at sufficiently high flow rates due to the high Weissenberg number problem (Owens & Phillips 2002), which is of purely numerical origin. Although significant progress has been made in developing numerical methods capable of controlling this problem (Alves et al. 2021), direct numerical simulations relevant to experiments in parallel shear flows (Bonn et al. 2011;Pan et al. 2013;Qin & Arratia 2017) are still lacking. Here, we propose to circumvent the high Weissenberg number problem by studying parallel shear flows with free slip boundary conditions that remove steep near-wall velocity gradients and significantly improve numerical stability.
The influence of wall slip on hydrodynamic stability of various flows has been extensively studied. The motivation behind these studies can be separated into two broad classes. The studies from the first class seek to approximate microscopically sound boundary conditions and thus address the question of linear stability of experimentally realisable flows. Such an approach has indicated, for instance, that the classical twodimensional Tollmien-Schlichting transition in plane channel flow of Newtonian fluids is significantly suppressed by the presence of wall slip (Gersting 1974;Lauga & Cossu 2005;Min & Kim 2005), although the non-normal energy growth mechanism was shown to be less affected (Lauga & Cossu 2005). When adopted to shear flows of polymeric fluids, this approach revealed the existence of a short wave-length instability (Black & Graham 1996, 2001, absent in the corresponding no-slip flows. The second class of studies employs slip boundary conditions exclusively as the means to simplify the underlying no-slip boundary condition problem. Free-slip boundary conditions usually adopted in those works are rather artificial, yet they yield significantly more tractable mathematical problems. This approach was pioneered by Lord Rayleigh (1916) in the context of Rayleigh-Bénard instability. More recently, Waleffe (1997) used free-slip boundary conditions to construct exact coherent states in parallel shear flows of Newtonian fluids that revolutionised our understanding of the transition to turbulence in those flows (Eckhardt et al. 2007;Tuckerman et al. 2020;Graham & Floryan 2021). The most important outcome of these studies is the observation that free-slip systems preserve the main phenomenology of their no-slip counterparts, albeit at increased values of the control parameter. Our work draws its motivation from this class of studies.
Here, we perform a linear stability analysis of plane Couette flow (pCF) and pressuredriven plane Poiseuille flow (pPF) with free-slip boundary conditions. We employ the Oldroyd-B model (Bird et al. 1987) in the absence of inertia, and demonstrate numerically that these flows develop linear instabilities absent in their no-slip counterparts. Importantly, we observe that the main features of no-slip flows are preserved under the free-slip conditions, and we discuss how these instabilities can be employed to gain insight into the nature of elastic turbulence in parallel geometries.

Problem setup
We study incompressible flows of a model viscoelastic fluid confined between two infinite parallel plates. We introduce a Cartesian coordinate system {x 1 , x 2 , x 3 } oriented along the streamwise, gradient, and spanwise directions, respectively; the plates are located at x 2 = ±L. The dynamics of the fluid are assumed to be described by the Oldroyd-B model (Bird et al. 1987), where v, τ , and p are the fluid velocity, polymeric Cauchy stress tensor, and the pressure, respectively, while (∇v) T denotes the transpose of the velocity gradient tensor; f is a body force to be specified below. These equations are rendered dimensionless using the channel's half-width L as the unit of length, the maximum velocity of the laminar profile V 0 as the unit velocity, and L/V 0 as the unit of time. The governing equations (2.1) feature three dimensionless groups, the Weissenberg number Wi = λV 0 /L, Reynolds number Re = ρLV 0 /η, and the viscosity ratio β = η s /η, where λ is the Maxwell relaxation time of the fluid, ρ is its density, while η s and η are the solvent viscosity and the total viscosity of the fluid, respectively. Unless stated otherwise, below we focus on the purely elastic limit with Re = 0.
The fluid is assumed to satisfy the following free-slip boundary conditions at the walls where ∂ 2 denotes the spatial derivative in the gradient direction. In this work, we study unidirectional laminar profiles v = {U (x 2 ), 0, 0}, where U is chosen to resemble the shape and symmetry of either plane Poiseuille flow (pPF) or plane Couette flow (pCF), while satisfying the boundary conditions, Eq.(2.2), The corresponding laminar components of the polymeric stress tensor are given by τ 11 = 2Wi (∂ 2 U (x 2 )) 2 and τ 12 = ∂ 2 U (x 2 ), with all other components being zero. To drive such flows in the presence of free-slip boundary conditions, we add the body force f = {−∂ 2 2 U (x 2 ), 0, 0} to the r.h.s. of Eq.(2.1b). To perform a linear stability analysis of these flows, we linearise Eqs.(2.1) around the laminar profile given above. In view of the viscoelastic analogue of the Squire's theorem (Bistagnino et al. 2007), that also holds for the free-slip boundary conditions studied here, we consider two-dimensional modal perturbations {δv(x 2 ), δτ (x 2 ), δp(x 2 )}e ik1x1 e σt of the velocity, stress, and pressure, respectively. Here, k 1 sets the spatial scale of the perturbation, and σ = σ r + iσ i denotes a complex temporal eigenvalue, with σ r and σ i being its real and imaginary part, respectively. The x 2 dependence of the perturbations are expanded in Chebyshev polynomials, and the linearised equations of motion are solved with a fully spectral Chebyshev-Tau method (Canuto et al. 1987) by employing the generalised eigenvalue solver eig from the scientific library SciPy for Python (Virtanen et al. 2020).

Results
The no-slip pCF and pPF of Oldroyd-B fluids are linearly stable for a wide range of β and Wi (Gorodtsov & Leonov 1967;Wilson et al. 1999 reported a purely elastic linear instability in no-slip pPF, however, its region of existence is confined to a small and, most likely, experimentally inaccessible part of parameter space with β > 0.9 and Wi > O(10 3 ). Surprisingly, here we find that free-slip pCF and pPF are linearly unstable for a wide range of β and Wi.
In Fig.1, we present examples of such instabilities by plotting the eigenvalue spectra for pCF and pPF. Similar to the no-slip problem (Wilson et al. 1999), the spectra consist of discrete physical eigenvalues and two lines of eigenvalues corresponding to continuous spectra with σ r = −1/Wi and σ r = −1/(βWi); the finite-resolution approximations to the latter take the balloon-like shapes visible in Fig. 1. The most profound feature of both spectra in Fig. 1 is that all leading eigenvalues have positive real parts, indicating the presence of a linear instability. Surprisingly, two of the three leading eigenvalues in pPF have numerically almost identical real parts for all values of the parameters we considered (see also Fig.6 below).
To demonstrate numerical convergence, in Fig.1, we present the eigenvalue spectra for two spectral resolutions, with M = 50 and M = 100 Chebyshev polynomials, respectively, and observe that the physical eigenvalues are well-converged. All calculations presented below are therefore carried out with M = 50. While this spectral resolution is quite low to yield converged results in the no-slip setting (Wilson et al. 1999;Khalid et al. 2021b), it is clearly sufficient for the free-slip boundary conditions. This observation supports our original rationale that the free-slip boundary conditions simplify numerical simulations of such flows by removing steep near-wall velocity gradients.
For no-slip boundary conditions, the least stable eigenmode of pCF and one of the three leading eigenmodes of pPF are mostly localised in the vicinity of the confining boundaries (Gorodtsov & Leonov 1967;Wilson et al. 1999), and are thus referred to as wall modes, while the other two pPF modes are mostly localised in the bulk (Wilson et al. 1999;Khalid et al. 2021b), i.e. they are centre modes. In Fig.2 we present the spatial profiles corresponding to the unstable free-slip eigenvalues discussed above, and observe that the wall versus centre mode distinction is significantly weakened. Although the symmetries of the stress distribution and the streamfunction are still different between these modes (compare the vortex structure in Fig.2(a)-(b) with (c)-(d), for instance), all modes now have significant presence throughout the domain. Notably, the fact that all these modes can be simultaneously unstable indicates the presence of several, rather different instability mechanisms.
In Fig.3  (1) Despite these similarities, we were unable to use the re-scaling procedure proposed in those works to collapse all our neutral stability data on a single master curve. Further, in Fig.4(a) we observe that Wi c increases rapidly as β → 1, indicating that a sufficiently strong velocity-stress coupling is required to drive the instability. We also note that Wi c is significantly smaller for pPF than for pCF, and that, surprisingly, the critical wave number k 1,c depends in different ways on β for these flows. As can be seen in Fig.4(b), in all cases, k 1,c ∼ O(1), indicating an instability lengthscale comparable with the wall-to-wall distance.
To assess the influence of weak inertia, in Fig.5 we plot the leading eigenvalue of pCF and pPF as a function of the Reynolds number. We see that while small amounts of inertia have stabilising effect on pCF, the low-β instabilities in pPF are only mildly suppressed. For pPF at sufficiently high β, the trend reverses, and we observe a stable purely elastic eigenmode becoming unstable for Re > 0. While a full Wi − β − Re critical surface is required to draw a definite conclusion regarding the influence of inertia on each mode, we note here that our observations are in line with the recent work on elasto-inertial instabilities (Chaudhary et al. 2019;Khalid et al. 2021a).
The spectra presented in Fig.1 bear a striking resemblance to the spectra of the corresponding flows of Oldroyd-B fluids with no-slip boundary conditions (Gorodtsov & Leonov 1967;Wilson et al. 1999) but for the position of the leading eigenvalues that appear to move to positive values of σ r in the free-slip case. To establish a connection between these eigenvalues and their no-slip counterparts, we perform a homotopy between the two types of boundary conditions. To this end, we introduce the following laminar profile that depends on a parameter α, (α = 0). A similar interpolation is imposed on the boundary conditions In Fig. 6 we show the superimposed spectra for this laminar profile as a function of α for pCF and pPF. The results of the homotopy procedure clearly indicate that the unstable modes do not appear either from infinity or the continuous spectra, but instead are continuously connected to the well-known least stable eigenvalues of no-slip parallel shear flows (Gorodtsov & Leonov 1967;Wilson et al. 1999).

Conclusions
In this work, we demonstrate that free-slip pCF and pPF of Oldroyd-B fluids are linearly unstable. The instabilities we observe are caused by neither the curvature of the base flow streamlines, where hoop stresses are known to lead to purely elastic instabilities (McKinley et al. 1996), nor by steep gradients of the material properties, as in the co-extrusion problems (Hinch et al. 1992). Suggestively, the base velocity profiles considered here change the sign of their second derivatives inside the domain, which might indicate a connection to Rayleigh's inflection point instability criterion for inviscid fluids (Drazin & Reid 2004). A similar situation occurs in purely elastic Kolmogorov flows (Boffetta et al. 2005), and it would be interesting to examine whether the instabilities observed there are related to those found in this work. By performing the homotopy between free-slip and no-slip boundary conditions, we demonstrate that the unstable eigenmodes are directly related to the least stable modes of the no-slip system. The exact way by which the introduction of an inflection point into the base profile would lead to destabilisation of the no-slip modes is currently unknown.
While the body forcing and the boundary conditions considered here are probably difficult, if not impossible, to realise in practice, the main importance of our work is not related to its direct experimental relevance. Instead, our contribution here is related to the observation that the free-slip problem considered here is significantly simpler to address numerically than its no-slip counterpart. Even at the linear level, a significantly lower number of Chebyshev modes is needed to resolve the eigenspectrum than in the no-slip case. A similar simplification is expected to also occur in fully non-linear direct numerical simulations and to alleviate the High Weissenberg number problem plaguing such studies. By rendering the leading eigenmodes unstable, free-slip flows offer an easier route to study strongly sub-critical structures that propagate from those eigenmodes in the no-slip case (Morozov & van Saarloos 2005. Once coherent states are resolved in direct numerical simulations of free-slip flows, a boundary condition homotopy, similar to the one employed here, can be used to track those solutions back to the no-slip conditions (Waleffe 2003). Recent work on elasto-inertial parallel shear flows (Garg et al. 2018;Chaudhary et al. 2019;Khalid et al. 2021a;Chaudhary et al. 2021) associates linear instabilities at high Wi and Re with the very leading eigenvalues studied here. As already mentioned above, the same instability can be traced to the purely elastic case, Re = 0, for sufficiently large values of β and Wi (Khalid et al. 2021b). It has also been shown that these instabilities are sub-critical in nature (Page et al. 2020). Taken together with the original proposal of Morozov & van Saarloos (2005, and the experimental and numerical studies of non-linear states in elasto-inertial flows (Samanta et al. 2013;Dubief et al. 2013;Sid et al. 2018;Choueiri et al. 2018;Lopez et al. 2019), this points towards the existence of exact coherent states related to the least stable eigenvalues, even though they are linearly stable for most values of β and Wi. Our work suggests using a homotopy from free-slip flows to discover these structures.
Finally, we note that our results suggest that all three leading eigenvalues of freeslip pPF can become unstable, see Fig.1, for example. In elasto-inertial flows, there is an ongoing discussion whether the most dynamically relevant coherent structures are related to wall or centre-line modes (Shekar et al. 2019;Datta et al. 2021). In Fig.7, we present the neutral stability curve for the first three eigenvalues and observe that by selecting the wavenumber and the Weissenberg number, one can control the type of the instability thus produced. As Fig.2 suggests, the different spatial symmetries of those modes should result in structures that connect to either wall or centre-line modes, when traced back to the no-slip conditions, thus allowing one to assess their relative importance independently.
Declaration of Interests. The authors report no conflict of interest.