Spatio temporal linear stability of viscoelastic free shear flows: dilute regime

TL;DR

This paper conducts a detailed linear stability analysis of viscoelastic free shear flows using the Oldroyd-B model, revealing how viscoelasticity influences flow stability and transition to instability.

Contribution

It provides the first quantitative linear stability analysis of dilute viscoelastic free shear flows, highlighting the effects of Weissenberg number on flow stability and instability mechanisms.

Findings

Viscoelasticity reduces the maximum growth rate of instabilities.

The unstable spectrum shifts towards longer waves with increasing Weissenberg number.

Flow becomes either absolutely or convectively unstable at very low Reynolds numbers.

Abstract

We report the temporal and spatio-temporal stability analyses of anti-symmetric, free shear, viscoelastic flows obeying the Oldroyd-B constitutive equation in the limit of low to moderate Reynolds number and Weissenberg number. The resulting fourth order Orr-Sommerfeld equation is reduced to a set of six auxiliary equations which are numerically integrated starting from the rescaled far-field conditions, i.~e., via the Compound Matrix Method. Numerical results indicate that with increasing Weissenberg number: (a) the peak of the maximal growth rate (i. e., the maximum value of the imaginary component of the growth rate) is reduced, (b) the entire range of unstable spectrum is shifted towards longer waves (i.~e., the entire region of temporal instability is gradually concentrated near zero wavenumber), (c) the vorticity structure contours are dilated and (d) the residual Reynolds stresses are diminished. All these observations suggest that viscoelasticity reduces the temporal instability but does not completely suppress it. The Briggs idea of analytic continuation is deployed to classify regions of temporal stability, absolute and convective instabilties as well as evanescent (false) modes, in a finite range of Reynolds number, Weissenberg number and the viscosity coefficient. The main result is that, unlike Newtonian fluids, the free shear flow of dilute polymeric liquids are either (absolutely/convectively) unstable for all Reynolds number or the transition to instability occurs at very low Reynolds number, a finding attributed to the fact that viscoelasticity aggravates free surface flow instabilities. Although this transitional pathway connecting the temporally stable state to the elastoinertially unstable state have been identified by some in-vitro experiments, but until now, has not been quantified theoretically via linear stability analysis.

Equations

$$\nabla \cdot \mathbf{v}=0,\rho \frac{D\mathbf{v}}{Dt}=-\nabla p+\mathbf{v}\cdot \tau ,$$

$$\tau =\eta [\nu \mathbf{D}+(1-\nu )\mathbf{A}],$$

$$\frac{\partial A}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{A}-\nabla {\mathbf{v}}^T\cdot \mathbf{A}-\mathbf{A}\cdot \nabla \mathbf{v}=\frac{\mathbf{D}-\mathbf{A}}{\lambda },$$

$$U(y)=\tanh (y),\Omega (y)={\tanh }^2(y)-1,\Psi (y)=\log (\cosh (y)),$$

$$\psi (x,y,t)=\Psi (y)+\varphi (y)e^{i(\alpha x-\omega t)},\tau (x,y,t)=\mathbf{T}(y)+\phi (y)e^{i(\alpha x-\omega t)},$$

$$\displaystyle i\big{\{}(\alpha U-\omega)(\phi^{\prime\prime}-\alpha^{2}\phi)-\alpha U^{\prime\prime}\phi\big{\}}-\frac{\nu}{Re}\left((\frac{d}{dy})^{2}-\alpha^{2})\right)^{2}\phi=\frac{1-\nu}{\mathcal{F}Re}\sum_{n=0}^{4}c_{n}\phi^{(n)},$$

$$\displaystyle\begin{split}c_{0}(y)&=\alpha^{4}+\mathcal{F}^{(4)}-\alpha^{2}\bigg{(}\frac{\mathcal{F}^{\prime\prime}}{\mathcal{F}}\bigg{)}(\mathcal{F}-1)-2\alpha^{2}\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}^{2}(\mathcal{F}^{2}+1)-4\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}\mathcal{F}^{\prime\prime\prime}-3\frac{(\mathcal{F}^{\prime\prime})^{2}}{\mathcal{F}}\\ &+~{}4\frac{(\mathcal{F}^{\prime})^{4}}{\mathcal{F}^{2}}-6\mathcal{F}^{\prime\prime}\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}^{2}(\mathcal{F}-1),\\ c_{1}(y)&=-2\alpha^{2}\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}(\mathcal{F}-1)+4\bigg{(}\frac{\mathcal{F}^{\prime\prime}}{\mathcal{F}}\bigg{)}\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}(\mathcal{F}-1)^{2}-4\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}^{2}(\text{D}_{y}\mathcal{F})(\mathcal{F}-1)\\ &+~{}2\bigg{(}\frac{\mathcal{F}^{\prime\prime\prime}}{\mathcal{F}}\bigg{)}(\mathcal{F}-1),\\ c_{2}(y)&=-2\alpha^{2}+3\bigg{(}\frac{\mathcal{F}^{\prime\prime}}{\mathcal{F}}\bigg{)}(\mathcal{F}-1)+2\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}^{2}(\mathcal{F}-1)^{2},\\ c_{3}(y)&=2\bigg{(}\frac{\mathcal{F}^{\prime}}{\mathcal{F}}\bigg{)}(\mathcal{F}-1),\\ c_{4}(y)&=1,\end{split}$$

$$\varphi ={\varphi }^{\prime \prime }=0$$

$$\varphi ={\varphi }^{\prime }=0$$

$$D\left(-i\frac{\partial }{\partial x},i\frac{\partial }{\partial t},\mathbf{M}\right)u(x,t)=\delta (x)\delta (t),$$

$$G(x,t)=\frac{1}{4{\pi }^2}{\int }_L{\int }_F\frac{e^{i(\alpha x-\omega t)}}{D(\alpha ,\omega ,\mathbf{M})}d\alpha d\omega$$

$$G(x,t)\sim -\frac{1}{\sqrt{2\pi }}\frac{e^{i[\pi /4+{\alpha }_{\ast }x-{\omega }_{\ast }t]}}{\frac{\partial D}{\partial \omega }{\left[\frac{d^2\omega }{d{\alpha }^2}\right]}^{1/2}},$$

$$G(x,t)\stackrel{t\to \infty }{\to }\infty ,\text{along the ray}\text{nicefrac}xt=0,$$

$$G(x,t)\stackrel{t\to \infty }{\to }0,\text{along the ray}\text{nicefrac}xt=0.$$

$${\varphi }^{(4)}-2{\alpha }^2{\varphi }^{\prime \prime }+{\alpha }^4\varphi =\frac{iRe{\mathcal{F}}_{\infty }}{1-\nu +\nu {\mathcal{F}}_{\infty }}(\alpha -\omega )({\varphi }^{\prime \prime }-{\alpha }^2\varphi ),$$

$$\varphi =a_1{\varphi }_1+a_3{\varphi }_3.$$

$$\left({\varphi }_1{\varphi }_3^{\prime \prime }-{\varphi }_1^{\prime \prime }{\varphi }_3\right){|}_{y=0}=0,$$

$$y_1$$

$$y_2$$

$$y_3$$

$$y_4$$

$$y_5$$

$$y_6$$

$$y_1^{\prime }$$

$$y_2^{\prime }$$

$$y_3^{\prime }$$

$$y_4^{\prime }$$

$$y_5^{\prime }$$

$$+\frac{\mathcal{F}(i\text{Re}[(\alpha U-\omega ){\alpha }^2+\alpha U^{\prime \prime }]+{\alpha }^4\nu )y_1}{(1-\nu +\mathcal{F}\nu )},$$

$$y_6^{\prime }$$

$$y_1$$

$$y_2$$

$$y_3$$

$$y_4$$

$$y_5$$

$$y_6$$

$$y_1$$

$$y_2$$

$$y_3$$

$$y_4$$

$$y_5$$

Full text

Keywords: Orr-Sommerfeld equation, Compound Matrix method, Oldroyd B flows, Viscoelastic stabilization, Absolute and convective instabilities

1 Introduction

Understanding the hydrodynamic stability and transition of free shear flows of dilute polymer solutions is paramount to both the fundamental theory of viscoelastic liquids [1, 2] and their industrial applications, especially those arising in microfluidic mixing [3], viscoelastic thin films [4] and viscoelastic stabilization via polymer addition [5]. Likewise, understanding the mechanism governing the transition to turbulence in shearing flows of viscoelastic biofluids is decisive in resolving the swelling instability transition in mucus [6], cartilage [7], adhesion-fragmentation transition in cells [8, 9, 10] and bending instability transition in soft tissues [11]. Classical approaches to the viscoelastic instability studies involve identifying an equilibrium state, whose stability is studied through eigenvalue analysis by linearizing the governing equations. The eigenvalue analysis seek the least stable eigenmode of the linearized mass and momentum conservation equations suitably transformed to the Orr-Sommerfeld Equation (OSE) [12, 13] via parallel flow approximation. One of the features of traditional eigenvalue analysis is that the disturbance field is assumed to grow either in space or in time. Huerre and Monkewitz [14] have attempted to apply the combined spatio-temporal theory to a restricted class of Newtonian mixing layers with the goal of determining a general criterion whereby a class of flows can be analyzed by either the spatial or temporal theory by inspecting the dispersion relation (DRP) in the complex wave number-frequency plane. Some other applications as well as in-vivo and in-silico stability studies of shear flows of Newtonian and non-Newtonian liquids are listed in [15, 16, 17]. However, the models and experiments listed in these references fail to present a comprehensive, spatio-temporal analyses of plane, free shear flows in the form of an instability-phase diagram. Therefore a motivation of this article is to illustrate this phase diagram in the viscoelastic parameter space.

The characterization of fluid flows as being absolutely or convectively unstable and the method of spatio-temporal analysis by progressive moving of the isocontours in the complex frequency and wavenumber plane was first proposed by Briggs [18] and later by Bers [19] in the context of plasma physics. Rallison et al. detailed separate studies of temporal as well as spatial instability in channel flows of dilute suspensions, including stratified flows [20] and multilayer Couette and Poiseuille flows [21], obeying either the Oldroyd-B, Upper Convected Maxwell or the FENE fluid constitutive equations. The literature on theoretical studies of the spatio-temporal linear analysis of viscoelastic flows are more recent. Govindarajan highlighted the role of the spatio-temporal hierarchical structure of complex fluids, in determining the alterations of flow stability [22]. A series of early [23] and later experiments in viscoelastic, wall bounded flows [24], extensional flows [25], torsional flows [26] and in Taylor-Couette flows [27], outline the transition to turbulence at Reynolds numbers ($Re$) much lower than the Newtonian threshold, thereby dubbing this phenomena as ‘early turbulence’. Consequently, this ‘new’ transitional pathway, connecting the laminar to a novel elastoinertial turbulent state is speculated to be either captured only by full Direct Numerical Simulations (DNS) [28] or by purely elastic subcritical (or non-linear) instability analysis [29]. However, this article offers a crucial conceptual explanation of the above mentioned transitional pathway to early turbulence via a linear, spatio-temporal stability analysis.

Through linear stability analysis, we aim to address the following intriguing questions: What is the critical flow / polymer relaxation conditions for the onset of instability ? and more crucially What is the linear spatio-temporal, time asymptotic response of the flow at the critical value of the material parameters ? In the next section, we delineate the model of the free shear layer flow coupled with the Oldroyd-B constitutive relation for the extra elastic stress tensor (§2.1) and the details of the linear stability analysis via the fourth order OSE (§2.2). A thorough description of the Briggs contour integral method to determine the existence of absolute instability, convective instability and evanescent modes in spatio-temporal analysis in included in this section. In §3, we introduce the Compound Matrix Method to numerically integrate the resultant system of stiff differential equations for the auxiliary variables emerging from the Orr-Sommerfeld equation. Section 4 highlights the simulation results of the temporal (§4.1) and the spatio-temporal (§4.2) stability analyses of anti-symmetric, free shear flows of polymeric liquids. This article concludes with a brief discussion of the implication of these results and the focus of our future direction (§5).

2 Mathematical model and linear stability analysis

Unlike Newtonian solvents, the stability predictions of polymeric liquids depend upon the details of the equations relating stress to the shear rate. The Oldroyd-B constitutive equations has its range of applicability limited to dilute solutions and moderate shear rates [2]. It predicts no shear thinning, a constant first normal stress coefficient and a zero second normal stress coefficient. The numerical solution obtained by linearizing the Navier Stokes along with the extra stress constitutive equation (or the OSE stability equation) is reliable in describing the characteristics of the initial stages of the mixing layer transition [1]. The linear stability results can be utilized as the initial-boundary conditions for DNS [13] or the large eddy simulation of polymeric liquids [28] in order to further study the flow evolving process downstream.

2.1 Problem definition

For the mixing layer flow configuration it is customary to assume $\mathcal{U}_{1}$ (respectively $\mathcal{U}_{2}$) as the free-stream velocity of the upper (lower) flow. We denote $\mathcal{U}_{0}=\frac{1}{2}(\mathcal{U}_{1}-\mathcal{U}_{2})$ as the free-stream velocity in a reference frame moving with the average velocity of the flow (i. e., $\frac{1}{2}(\mathcal{U}_{1}+\mathcal{U}_{2})$) and $\delta$, the momentum thickness of the mixing layer. The continuity and the momentum equations for an incompressible flow are,

[TABLE]

where ${\bf v}$ is the velocity vector, $\rho$ the density, $p$ the isotropic pressure and $\tau$ the extra stress tensor. Equations (1) are closed through the evolution equation for the extra stress tensor $\tau$, written as the sum of the viscous, Newtonian stress, $\tau^{s}$ (= $\eta_{s}{\bf D}$, ${\bf D}$ is the shear rate tensor) and the elastic stress, $\tau^{p}$ (= $\eta_{p}{\bf A}$), where $\eta_{s},\eta_{p}$ are the solvent viscosity and the polymeric contribution to the shear viscosity, respectively. We let $\eta(=\eta_{s}+\eta_{p})$ and $\nu(=\eta_{s}/\eta)$ denote the total viscosity and the viscous contribution to the total viscosity of the fluid, respectively. The constitutive equation for the extra stress tensor, $\tau$, obeying the Oldroyd-B model then becomes,

[TABLE]

where the tensor ${\bf A}$ satisfies the Upper Convected Maxwell equation,

[TABLE]

$\lambda$ is the polymer relaxation time.

2.2 Linear Stability Analysis

Using $\mathcal{U}_{0}$ and $\delta$ as the reference velocity and length scale, respectively, we characterize the non-dimensional form of equations (1), by the dimensionless numbers, $Re=\nicefrac{{\rho\delta\mathcal{U}_{0}}}{{\eta}}$, and the Weissenberg number, $We=\nicefrac{{\lambda\mathcal{U}_{0}}}{{\delta}}$. We assume that the mean flow is two-dimensional and quasi-parallel with its variation entirely in the direction normal to the flow, i. e.,

[TABLE]

where $U(y),\Omega(y),\Psi(y)$ are the dimensionless streamwise mean velocity (with zero spanwise mean velocity), spanwise mean vorticity and the associated streamfunction, respectively. Further, assume that the mean flow supports a two-dimensional disturbance field. The streamfunction and the extra stress tensor are represented by the base state profile $(\Psi(y),\mathbf{T}(y))$ plus a small perturbation, Fourier transformed in $x$ and $t$ as follows,

[TABLE]

where $\phi(y),\varphi(y)$ are the spanwise perturbations in the streamfunction and the extra stress tensor and $\alpha,\omega$ are the complex wavenumber and angular frequency, respectively. Rewriting equation (1) in the streamfunction-vorticity formulation and utilizing equations (2), (3), (4), (5), we arrive at the equation governing the perturbation of the streamfunction, given by the familiar fourth order OSE [30],

[TABLE]

and $\mathcal{F}=1+iWe(\alpha U-\omega)$, solved together with the vanishing boundary conditions for the free shear flow disturbance ($\phi\rightarrow 0,\phi^{\prime}\rightarrow 0$ as $y\rightarrow\pm\infty$). Limiting the domain of integration to the upper half of the flow and analyzing only the anti-symmetric disturbance, the boundary conditions are revised as follows,

[TABLE]

In § 3, we delineate a numerical solution procedure for solving the eigenvalue equation (6) together with the boundary conditions (7a)(7b). In the rest of this section we describe a mechanism to classify the spatio-temporal flow instabilities.

The spatio-temporal evolution of a localized disturbance (located at the origin of the $x$-$t$ plane) is illustrated by considering the response of a given base velocity profile, $u(x,t)$, to an impulse excitation [14],

[TABLE]

where $D(\alpha,\omega,{\bf M})=0$ is the DRP relation and ${\bf M}$ is the vector of material and fluid parameters. The solution to equation (8) is dictated by the Green’s function,

[TABLE]

where $F$ and $L$ are the Fourier contour in $\alpha$ and the Laplace contour in the $\omega$ plane. The Fourier contour integral is placed parallel to the Real($\alpha$) axis (in the ensuing description we denote real/imaginary components with subscript r/i, respectively), while the Laplace contour integral is placed above all singularities of the DRP in the $\omega$-plane so as to satisfy the causality condition (i. e., $G\equiv 0$ if $t<0$). An important criteria for the understanding of instability entails the study of the flow behavior in the ‘long term’ (i. e., $t\rightarrow\infty$) where an analytical solution of equation (8) is possible. In the asymptotic limit of long time, the integration of equation (9) is analytically accomplished by using the method of stationary phase [31], i. e.,

[TABLE]

where $\alpha_{*}$ is the saddle point in the $\alpha$ plane (i. e., the root of $\frac{\partial D}{\partial\alpha}=0$) and $\omega_{*}$ is the corresponding branch point in the $\omega$ plane satisfying the DRP relation. From equation (10), we can surmise the condition for which the flow will be absolutely unstable (signifying the growth of disturbance in both upstream and downstream direction from the origin, otherwise known as the ‘resonance mode’ [32]), i. e.,

[TABLE]

versus when the flow will be convectively unstable (where disturbances are swept downstream from the source and given sufficient time these disturbances decay at any fixed position in space, also known as the ‘driven mode’), i. e.,

[TABLE]

Our analysis also reveal the presence of evanescent modes (or false modes) in the flow field, which is the non-propagating mode or the locally concentrated mode [33], described later in this section.

The necessary (but not sufficient) condition for the presence of absolute instability is the vanishing characteristic of the group velocity, ${\it v}_{g}$, at the saddle point in the $\alpha$-plane or the branch point in the $\omega$-plane (${\it v}_{g}=\frac{\partial\omega}{\partial\alpha}=\nicefrac{{(\frac{\partial D}{\partial\alpha})}}{{(\frac{\partial D}{\partial\omega})}}=0$ such that $\omega=D(\alpha)$). But the group velocity is zero at every saddle point, especially where the two $\alpha$-branches meet, independent of whether the branches originate from the same half of the $\alpha$-plane (i. e., when evanescent modes are detected) or not. To overcome this inadequacy, Briggs [18] devised the idea of analytic continuation in which the Laplace contour $L$, in equation (9), is deformed towards the $\omega_{r}$ axis of the complex $\omega$-plane, with the simultaneous adjustment of the Fourier contour $F$ in the $\alpha$-plane to maintain the separation of the $\alpha$-branches; those which originate from the top half (the upstream modes with $\alpha_{i}>0$) from those which originate from the bottom half of the $\alpha$-plane (or the downstream modes). The deformation of the $F$ contour (while preserving causality) is inhibited, however, when the paths of the two $\alpha$-branches originating from the opposite halves of the $\alpha$-plane intersect each other, leading to the appearance of saddle points which are the pinch point, $\alpha^{\text{p}inch}$. The concurrent branch point appearance in the $\omega$-plane is the cusp point, $\omega^{\text{c}usp}$ (i. e., $D(\alpha^{\text{p}inch},\omega^{\text{c}usp})\!=\!\frac{\partial D(\alpha^{\text{p}inch},\omega^{\text{c}usp})}{\partial\alpha}\!=\!0$ but $\frac{\partial^{2}D(\alpha^{\text{p}inch},\omega^{\text{c}usp})}{\partial\alpha^{2}}\!\neq\!0$). Kupfer [34] employed this local mapping procedure to conceptualize the stability characteristics of this branch point. Near a ‘reasonably close’ neighborhood of the pinch point, a local Taylor expansion yields a DRP relation which has a second-order algebraic form in the $\omega$-plane (and which is a first order saddle point in the $\alpha$-plane), i. e., $(\omega-\omega^{\text{c}usp})\sim(\alpha-\alpha^{\text{p}inch})^{2}$. This period-doubling characteristics of the map causes the $\alpha_{i}$-contours to ‘rotate’ around $\omega^{\text{c}usp}$, forming a cusp. In the $\omega$-plane, we draw a ray parallel to the $\omega_{i}$-axis from the cusp point such that it intersects the image of the F-contour (or $\alpha_{i}=0$ curve) and count the number of intersections (consequently, count the number of times both $\alpha$-branches cross the $\alpha_{r}$-axis before forming a pinch point in the $\alpha$-plane, as shown in Figure 5). If the ray drawn from the cusp point intersects the image of the F contour in the $\omega$-plane (or if either one or both the $\alpha$-branches cross the $\alpha_{r}$-axis) even number of times, then the flow dynamics correspond to an evanescent mode. Otherwise, in the case of odd intersections the observed cusp point is genuine, leading to either absolutely unstable system (in the upper half of the $\omega$-plane) or convectively unstable system (in the lower half of the $\omega$-plane); provided the system is temporally unstable.

3 Solution to the eigenvalue problem

The eigenvalues (i. e., the values of ($\alpha,\omega$) satisfying equation (6)) are found by examining the consequence of the far stream boundary condition, given by equation (7b), on the solution structure of the OSE. In the limit $y\rightarrow\infty$, we have $U(y)=1$ and $U^{\prime\prime}(y)=0$, which reduces the OSE to the following constant coefficient ODE [12],

[TABLE]

where $\mathcal{F}_{\infty}=1+We(\alpha-\omega)$. The solution to equation (13) can be obtained in the form $\phi=e^{\lambda y}$, with the characteristic roots given by $\lambda_{1,2}=\mp\alpha$ and $\lambda_{3,4}=\mp q$, where $q={\displaystyle\left[\alpha^{2}+\frac{{\it i}Re\mathcal{F}_{\infty}}{1-\nu+\nu\mathcal{F}_{\infty}}(\alpha-\omega)\right]^{1/2}}$. The fourth order OSE (6) will have four fundamental solutions, $\{\phi_{i}\}^{4}_{i=1}$, i. e., $\phi={\displaystyle\sum^{4}_{i=1}a_{i}\phi_{i}}$. To satisfy the boundary condition (7b), one must have $a_{2}=a_{4}=0$ for real ($\alpha,q$) ¿ 0 which warrants a general solution of the form

[TABLE]

Equation (14) is a non-trivial, admissible solution of the OSE, satisfying the zero boundary conditions at the centerline (i. e., at $y=0$, equation (7a)) if and only if the determinant of the associated matrix of the linear algebraic system vanishes at $y=0$, i. e.,

[TABLE]

which is the dispersion relation of the problem. In the asymptotic limit of $Re\rightarrow\infty$, the eigenmodes of the free shear layer instability problems are such that $|q|\gg|\alpha|$. This enormous contrast between the two sets of roots of the characteristic equation (13) is the source of stiffness causing the fundamental solutions of the OSE to vary by several orders of magnitude in the entire physical domain, and thus necessitate the use of Compound Matrix Method (CMM) [35]. In CMM, one works with a set of auxiliary variables which are combinations of the fundamental solutions $\phi_{1}$ and $\phi_{3}$, namely

[TABLE]

satisfying the corresponding first order ODEs,

[TABLE]

Since, at $y\rightarrow\infty$: $\phi_{1}\sim e^{-\alpha y}$ and $\phi_{3}\sim e^{-qy}$, we can estimate the free stream values of the unknowns,

[TABLE]

Since all the variables in equation (18) attenuate at the same rate, the stiffness of the problem is removed by rescaling the variables, say with respect to $y_{1}$, so that the (rescaled) free-stream conditions for solving equations (17) are

[TABLE]

The solution for equations (17) is obtained by marching backward from the free stream to the centerline, with the initial conditions provided by equations (19). A suitable value of ($\alpha,\omega$), or the eigenvalue for a given combination of material parameters M is obtained by enforcing the DRP relation (equation (15)) in auxiliary variables, i. e.,

[TABLE]

The eigenfunction, $\phi$, are found using a suitable linear combination of $\phi_{1}$ and $\phi_{3}$, via equations (14), (16), (17), such that the characteristic roots at free stream (i. e., the solution to equation (13)) have a negative real part and the coefficients in the ODE do not vanish at the centerline. An equation satisfying these constraints was suggested by Davey [36] as follows,

[TABLE]

4 Results

In our in-silico experiments, the zeros of the DRP relation (equation (15), (20)) were found in the $\alpha/\omega$ plane inside the region $|\omega_{r}|\leq 1.5,|\omega_{i}|\leq 1.5,|\alpha_{r}|\leq 3.0$ and $|\alpha_{i}|\leq 0.1$, using a discrete step-size of $\Delta\alpha_{\text{r}/i}=\Delta\omega_{\text{r}/i}=5\times 10^{-3}$. The free stream boundary conditions were applied at $\eta=12$, which leads to the numerical value of the momentum thickness, $\delta=0.30685$. Numerical integration from the free stream to the shear layer centerline was accomplished via the fourth order Runge Kutta time integration with discrete step-size of $\Delta y=6.67\times 10^{-4}$ in the tranverse direction. The eigenfunctions are found by numerically integrating equation (21) from the centerline to the free stream starting from the centerline conditions, equation (7a). The presence of only two centerline conditions imply that the eigenfunctions are scaled with respect to the centerline value of $\phi^{\prime}$ (i .e., $\phi^{\prime}(0)$). Hence, the results discussed in §4.1, 4.2, involving these eigenfunctions, are qualitative.

4.1 Temporal stability analysis

In this section we explore the stability of the solution of the OSE (equation (6)) by exclusively assigning $\omega$ to be a complex number. The temporal stability analysis in the inviscid limit was earlier studied by Azaiez [30], but owing to the limitations in their analysis some of the conclusions obtained in their work need to be revisited. In Figure 1 we present the maximum temporal growth rate of the disturbance (i. e., $\omega^{\text{Temp}}_{i}$, the largest positive imaginary component of any root of the DRP relation (15)) versus the wavenumber, $\alpha_{r}$. Note that, as $We$ is increased, the peak of the maximum growth rate, $\omega^{\text{max}}_{i}$ as well as the range of the unstable wavenumbers (i. e., the range of $\alpha_{r}$ in which $\omega^{\text{Temp}}_{i}>0$) is reduced and the entire temporally unstable spectrum is shifted towards longer waves (i. e., the region of instability is gradually concentrated around $\alpha_{r}=0$). All these observations suggest a mechanism of elastic stabilization in non-Newtonian fluids which we have discussed below. Later, in §4.2, we outline how the nature of this instability changes from absolute instability to convective instability.

Further characterization of the presence of viscoelasticity on the flow is highlighted via the ratio of the transverse to the streamwise velocity fluctuations (Figure 2(a)) and the variations of the magnitude of the vorticity disturbance (Figure 2(b)). These figures are presented for six different combinations of ($Re,We$) evaluated at the peak of the maximum growth rate curves, $\alpha^{\text{max}},\omega^{\text{max}}$ (refer Figure 1 caption for the respective numerical values). Figure 2(a) depicts the relative strength of the transverse velocity fluctuations versus tranverse spatial coordinates. We remark that at a lower value of $We$ (e. g., $We=1.0$) these relative fluctuations are significant within the viscous mixing region ($y\leq 1$) whereas at higher $We$ values (e. g., $We=40.0$) these fluctuations become important outside the mixing layer region and in a region where the fluid is essentially moving at the free-stream velocity (i. e., at $We=40.0$ the velocity fluctuation ratio attains its maximum at $y\approx 3.8$ where the mean velocity U${}_{0}(y)$ is 99% of its free-stream value). This study suggests that in traditionally non-Newtonian shear flows, the dynamical effects of the disturbance on the mean flow are negligible when compared to an equivalent Newtonian case [13]. This observation is further corroborated by observing the magnitude of the vorticity disturbance (Figure 2(b)). While at lower $We$ ($We\leq 10.0$) the region of maximum magnitude of vorticity disturbance is inside the range $y\in[0,1]$, at higher $We$ this maximum occurs in the range, $y\geq 1$.

More insights about this mechanism of stabilization can be gained by comparing the contours of the rate of vorticity production (drawn at identically equal levels) for an arbitrary time in the streamwise versus transverse direction, at $Re=400,We=1.0$ (Figure 3(a)) and $Re=400,We=40.0$ (Figure 3(b)) evaluated at the eigenvalues corresponding to the peak of the growth rate curves (refer Figure 1 caption). Note that as $We$ increases, the vorticity structures become larger and the contour spacing becomes wider, thereby suggesting an enhanced drag reduction.

An equivalent outcome can be drawn by examining the root mean square of streamwise velocity fluctuations, $u_{\text{rms}}(0,y)$ (Figure 4(a)) and Residual Reynolds shear stress (Figure 4(b)) at a fixed streamwise location, $x_{0}=0$. A larger drag reduction rate (associated with higher $We$) begets a larger $u_{\text{rms}}$ (e. g., $u_{\text{rms}}\approx 1.5,6.0$ at $We=1.0,40$ respectively) and a lower residual Reynolds stresses in an absolute sense. To summarize the viscoelastic free shear flows are less unstable in the presence of viscoelastic additives.

4.2 Spatio-temporal stability analysis

The spatio-temporal analysis, where both $\omega$ and $\alpha$ are designated as complex numbers, is a step which is necessary to reveal the stability phase diagram in the flow-material parameter space (i. e., in the ($Re,We,\nu$)-space, Figure 8). In the onset, we demonstrate an example of a cusp point in the $\omega$-plane (Figure 5(a)) and the saddle point in the $\alpha$-plane (Figure 6) which correspond to each other via the local angle-doubling map as indicated by Kupfer [34] (details in § 2.2). To this end, a particular set of parameters are chosen for which the flow is absolutely unstable: $Re=100,We=1,\nu=0.9$. First, the image of the dispersion relation at $\alpha_{i}=0$ is obtained (Figure 5(a)). This curve is the root of the DRP relation (equation (15)) with the maximum positive imaginary component of the frequency, $\omega_{i}^{\text{max}}$, as outlined in §4.1. As we decrease the imaginary component of the wavenumber (for downstream mode in this example), the corresponding image in the $\omega$-plane is observed. When $\alpha_{i}$ is sufficiently negative (e. g., $\alpha_{i}=-0.0301$), we detect the appearance of the cusp point at $\omega^{\text{cusp}}=0.6268+0.2362{\it i}$. The pinch point in the $\alpha$-plane (which is also the saddle point) corresponding to the cusp point in the $\omega$-plane arises at $\alpha^{\text{pinch}}=-0.041-0.0301{\it i}$ (Figure 6). The pinch point is obtained by drawing the isocontours of $\omega_{r}$ (Figure 6(a)) and $\omega_{i}$ (Figure 6(b)) in the $\alpha$-plane. Consequently, by drawing a ray parallel to the $\omega_{i}$ axis (Figure 5(a)), we note that this ray intersects $\alpha_{i}=0$ curve only once (i. e., odd number of times), implying that the cusp point is genuine. Since the imaginary component of this cusp point is positive ($\omega_{i}^{\text{cusp}}=0.2362$), the free shear viscoelastic fluid flow at $Re=100,We=1,\nu=0.9$ is absolutely unstable.

Evanescent modes are also encountered in our analysis and to illustrate their existence we consider another example of the cusp point formation at $Re=100,\,\,We=1\,\,\,\nu=0.5$ (Figure 5(b)). In this case, a ray emerging from the cusp point and parallel to the $\omega_{i}$ axis, intersects the $\alpha_{i}=0$ curve twice (i. e., even number of times), thereby establishing that this cusp point corresponds to an evanescent mode. Evanescent modes (also known as the direct resonance mode [37]) arises if the two coalescing modes (given by the solution to equation (8)) originate from waves propagating in the same direction (detailed discussion in §2.2). Linear stability theory predicts that these disturbances decay in the asymptotic limit of long time, but in the short term an algebraic growth associated with such a mode may be decisive to cause the transition to turbulence in receptivity studies [32].

In order to determine the range of $Re$ (for fixed $\nu$, and different $We$) for which the flow regimes are absolutely unstable, convectively unstable or temporally stable, we plot the absolute growth rate in Figure 7 (i. e., $\omega_{i}^{\text{c}usp}$, which is the growth rate at the cusp point obtained by deforming the Fourier contour in the $\alpha$-plane starting from the peak of the temporal growth rates, as explained in the above two examples). We define $Re_{s}$ ($We_{s}$) as the maximum critical Reynolds number (minimum critical Weissenberg number) for the flow to be temporally stable. Notice that the magnitude of the absolute growth rate (wherever it exists) predominantly decreases with increasing $Re$ and with increasing $We$. Although the combined effect of these parameters on the flow stability are often reported in terms of the elasticity number ($E_{\it l}=We/Re$, which corresponds to the ratio of the time scale for the elastic stress evolution to the time scale for diffusion of momentum [30]), we elaborate the impact of these flow-material parameters via the stability phase-diagram, discussed next.

We depict the boundaries of the temporally stable regions (S), convective instabilities (C), evanescent modes (E) and absolute instabilities (A) within a selected range of flow-material parameter space, i. e., $Re\in[0.1,500]$, $We\in[1,1000]$ and $\nu=0.35$ (Figure 8(a)), $\nu=0.20$ (Figure 8(b)). The boundaries of temporally stable regions are estimated to reside within the range $We\in[We_{s}=505,1000]$, $Re\in[0,Re_{s}=24.9]$ at $\nu=0.35$ (refer Figure 8(a)) and $We\in[We_{s}=455,1000]$, $Re\in[0,Re_{s}=20.1]$ at $\nu=0.20$ (refer Figure 8(b)), respectively. The linear stability phase diagram for the free shear flow systems dominated by viscous stresses (i. e., the case in which $\nu>0.35$), have revealed only the absolutely unstable and the evanescent modes within an identical range of the flow parameters. Hence, the discussion of those phase diagrams are omitted in this article. Two conclusions can be deduced from these phase diagrams. First, since the boundary of the absolutely unstable region, (A, E), is not parallel to the y-axis, the authors surmise that the flow will be absolutely unstable for sufficiently high $Re$ and irrespective of $We$. This assertion is supported by previous numerical studies of inviscid viscoelastic mixing layer in the limit of high $Re,We$ such that the elasticity number is held fixed and finite [30]. Second, unlike Newtonian fluids, viscoelastic liquids are either (absolutely/convectively) unstable for all $Re$ (e. g., consider the region $We<455$ in Figure 8(a) and 8(b)) or the transition to instability occurs at very low $Re$ (e. g., $Re_{s}=24.9,20.1$ at $\nu=0.35,0.20$ and $We=1000$, respectively). Further in the latter case, the transition is direct from temporally stable region to absolutely unstable region for elastic stress dominated fluids (i. e., the case of $\nu=0.35$, Figure 8(a)) while this transition is more gradual, ${\bf S}\rightarrow{\bf C}\rightarrow{\bf A}$ for viscous stress dominated fluids (i. e., the case of $\nu=0.20$, Figure 8(b)). Although absolute instability is far more vicious than convective instability as the occurrence of absolute instability pervades the entire flow domain, convective instability can also cause transition to turbulence if sufficient spatial distance is given for the disturbance to grow [14]. The second conclusion can be attributed to the fact that viscoelasticity dramatically exacerbates free surface flow instabilities [28]. This destabilization arises from a combination of large tangential tension and the curvature along the free surface leading to an unstable stress gradient.

We recapitulate our discussion by indicating some experimental findings which corroborate our numerical results, starting from the polymer drag reduction predictions in the early 1980s by Hibbert who used Laser-Doppler anemometry to notice the tendency for the enhancement of the vorticity structures and the diminution of small scale turbulence in free shear flow of polymeric liquids [38]. The flow visualization experiments of Riediger on surfactant blended fluids showed a longer lifetime of the typical structures with the mixing layer, i. e. roll-up and pairing, as well as a tendency for the suppression of the mixing layer growth [39]. More recently conducted experiments by Rahgozar et al. [40] have reported an order of magnitude reduction in the natural (unforced) transition $Re$ of boundary-free, uniformly sheared, dilute polymer solutions; a phenomena rebaptized as ‘early turbulence’ [28]. Subsequently, this new transitional pathway to the elastoinertial turbulent state have been experimentally demonstrated to exist over a much wider range of flow-material parameters [41].

5 Conclusions

This investigation addresses the linear, temporal and spatio-temporal analyses of free shear flows of dilute polymeric liquids with anti-symmetric centerline conditions for low to moderate Reynolds and Weissenberg numbers. §2 presented the free shear viscoelastic flow as well as the elements of the linear stability analysis via the solution of the Orr Sommerfeld equation, including the illustration of the Briggs contour integral method to determine the existence of absolutely unstable, convectively unstable and evanescent modes. §3 demonstrated the steps of the Compounded Matrix Method, utilized to numerical solve the resultant system of stiff differential equations. The numerical outcome, in a selected range of Reynolds and Weissenberg number, highlight that with increasing Weissenberg number the peak of the maximal growth rate is reduced, the region of temporal instability is gradually concentrated near zero wavenumber, the vorticity structures of the disturbance become larger and more widely spaced and the residual Reynolds stresses are lowered in absolute sense, all these observations indicating a mechanism of elastic stabilization in non-Newtonian free shear flows of dilute solutions. The spatio-temporal analysis reveal that either the flow is unstable for all Reynolds number or the transition to instability transpires at very low Reynolds number stipulating the significance of the unraveling of the polymer chains at the initiation of the free surface instabilities. Although the Oldroyd-B model describes well the behavior of polymeric liquids composed of a low concentration of high molecular weight polymer in a very viscous Newtonian solvent at moderate shear rates, more realistic description of polymeric solution and melts, including the effect of shear thinning, non-zero normal stress coefficient and stress overshoot in transient flows is missing and hence deserves a full numerical exploration in the near future, using a constitutive relation which incorporates these above mentioned features.

Acknowledgements:

Authors thank Prof. Tapan Sengupta (Dept. Aerospace Engg, IIT Kanpur) for introducing the Compound Matrix Method. SS acknowledges the financial support of the grant DST ECR/2017/000632. DB acknowledges the financial support of her fellowship CSIR 09/1117(0004)/2017-EMR-I which helped her run the numerical simulations.