Linear stability of the flow of a second order fluid past a wedge

The linear stability analysis of Rivlin–Ericksen fluids of second order is investigated for boundary layer flows, where a semi-infinite wedge is placed symmetrically with respect to the flow direction. Second order fluids belong to a larger family of fluids called order fluids, which is one of the first classes proposed to model departures from Newtonian behavior. Second order fluids can model non-zero normal stress differences, which is an essential feature of viscoelastic fluids. The linear stability properties are studied for both signs of the elasticity number K, which characterizes the non-Newtonian response of the fluid. Stabilization is observed for the temporal and spatial evolution of two-dimensional disturbances when K > 0 in terms of increase of critical Reynolds numbers and reduction of growth rates, whereas the flow is less stable when K < 0. By extending the analysis to three-dimensional disturbances, we show that a positive elasticity number K destabilizes streamwise independent waves, while the opposite happens for K < 0. We show that, as for Newtonian fluids, the non-modal amplification of streamwise independent disturbances is the most dangerous mechanism for transient energy growth, which is enhanced when K > 0 and diminished when K < 0. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0018300., s


I. INTRODUCTION
The mechanical behavior of many real fluids is described by Navier-Stokes theory. This theory is based on the assumption of a Newtonian constitutive equation. More specifically, the extra-stress tensor is expressed as a linear, isotropic function of the components of the velocity gradient. Many common fluids, such as water and air, can be assumed to be Newtonian. However, many rheologically complex fluids, such as polymer solutions, soaps, blood, paints, shampoo, and ketchup, exhibit a variety of phenomena that cannot be adequately described by a Newtonian constitutive equation.
Viscoelastic fluids are examples of non-Newtonian fluids, and they exhibit both viscous and elastic properties when undergoing deformation. Elasticity is the tendency of the material to return to its original shape once the external force is removed. Viscoelastic fluids undergo a gradual deformation and recovery when they are subjected to loading and unloading. The stress is directly proportional to neither the strain nor the rate of strain, and the relationship is more complex.
In this paper, we consider a subclass of differential type fluids known as the Rivlin-Ericksen fluids of second order. 1 In these models, only an infinitesimal part of the history of the deformation gradient has an influence on the stress. The extra stress is a function of the velocity gradient and its higher time derivatives. These materials lack a gradually fading memory, and they cannot represent the phenomenon of stress relaxation. However, they can predict non-zero normal stress differences, which is an important feature of viscoelastic fluids. As models to describe viscoelastic fluids, order fluids are suitable to describe slightly elastic fluids, where the fluid is only a small departure from the Newtonian fluid and flows for which the Rivlin-Ericksen tensors vary slowly. 2 In theoretical work, Rajagopal et al. 3 showed that it is possible to apply Prandtl's boundary layer theory to the case of a non-Newtonian fluid of second order. In particular, they showed that the equations of motion of a second order fluid can be satisfied by an irrotational flow, and they identified suitable assumptions to obtain a consistent theory. In the case of fluids of a differential type, the equations of motion are an order higher than the Navier-Stokes equations, and thus, the no-slip and no-penetration boundary conditions are insufficient to determine the solution completely. 4,5 The same is also true for the boundary layer approximation. In order to overcome this difficulty, in their study of an incompressible fluid of liquid B ′ near a stagnation point, Beard and Walters 6 suggested a perturbation method. This method was also adopted by Rajagopal et al. 7 in their analysis of the flow past a wedge of an incompressible fluid of second grade. The perturbation method reduces the order of the problem but is only valid for small values of the non-Newtonian parameter. This parameter multiplies the higher order spatial derivatives in the equation.
While studying flow near a stagnation point and flow past a wedge, Garg and Rajagopal 8,9 suggested that it would be preferable to use an augmented boundary condition justified by physically reasonable assumptions, also adopted by Vajravelu and Roper, 10 and by Vajravelu and Rollins. 11 Another difficulty that arises is the impossibility of finding a similarity solution to the boundary layer equations as in the Newtonian case, with the exception of stagnation flow. In this paper, we apply a pseudo-similarity transformation 9 to obtain the mean flow, which will be the starting point of the linear stability analysis.
The aim of this paper is to understand the stability properties of second order fluids in boundary layers. Specifically, a configuration of flow over a semi-infinite wedge is investigated (Fig. 1). One important motivation for studying the stability behavior of viscoelastic fluids, and in particular, polymer suspensions can be found in drag reduction in turbulent regime. [12][13][14] This phenomenon was first observed over 70 years ago. In turbulent boundary layers, dissolving a small quantity of long-chain flexible polymers in solution can reduce turbulent friction by a significant amount. For a summary of the historical findings and recent developments, we refer to the review by Xi. 15 First, we perform a classical linear stability analysis by solving the linearized system of equations around the steady mean flow, which is assumed to be parallel. The mean flow and the stability equations are solved numerically. For two-dimensional disturbances, the results are represented in terms of temporal growth rates, neutral stability curves, and critical Reynolds numbers. Later, the linear stability analysis is extended to three-dimensional disturbances.
Classical linear stability analysis is based on eigenvalues. However, in hydrodynamic stability and many other physical situations dominated by non-normal systems, eigenvalues prove to be misleading, and they do not describe correctly the whole dynamics. 16 In non-normal systems, such as Poiseuille, Couette, and Blasius flows, there can be short-time growth of energy even if all the eigenvalues decay exponentially. 17 This phenomenon is known as transient growth.
For Newtonian fluids, the possibility of transient growth has been known since the 1980s. 18 Some previous work has investigated transient growth of viscoelastic fluids in channel flows. 19 In this paper, we extend the analysis and consider the transient growth of second order fluids in boundary layers.

II. SECOND ORDER FLUIDS
The Cauchy stress tensor σ in a fluid of second order has the form 1 where p is the pressure, μ is the dynamic viscosity, and α 1 and α 2 (SI: kg/m) are the material moduli usually referred to as normal stress moduli. The spherical stress −pI is due to the constraint of incompressibility, while A 1 and A 2 are the Rivlin-Ericksen tensors of order 1 and order 2, respectively, defined by where v denotes the velocity field and D/Dt denotes the material time derivative. The sign of the material parameters in this model has been a source of some controversy. 20 In this paper, we consider both cases α 1 > 0 and α 1 < 0. The second order model with α 1 > 0 is studied because of its compatibility with thermodynamics. Since the form (1) is frame-indifferent, it can be used as an exact model. In this view, Dunn and Fosdick 21 and Fosdick and Rajagopal 22 justified some assumptions on the coefficients of the second order constitutive equation. In order for the fluid model to be compatible with thermodynamics, in the sense that all motions of the fluid satisfy the Clausius-Duhem inequality and the assumption that the specific Helmholtz free energy be a minimum in equilibrium, it is necessary that μ ≥ 0, α 1 ≥ 0, and α 1 + α 2 = 0. ( A detailed discussion of these assumptions can be found in the critical review of Dunn and Rajagopal. 20 The second order model with α 1 < 0 is studied because it gives the right sign for the first normal stress difference. 2 Moreover, in terms of linear stability, it is a consistent approximation to a stressrelaxing fluid, such as the Maxwell fluid, at small elasticity numbers and when the disturbance time scale is large compared to the characteristic time scale of the fluid. 23 To the best of our knowledge, little work has been performed on the stability of second order fluids in boundary layers unlike the situation for channel flows. Chun and Schwarz 24 studied the stability of the plane Poiseuille flow of a second order fluid (α 1 < 0). Their analysis yields an Orr-Sommerfeld equation modified by adding a non-Newtonian term. They showed that the critical Reynolds number decreases as the magnitude of the non-Newtonian Regarding the stability of other viscoelastic fluids, some results were obtained for channel flows but, to the best of our knowledge, not much has been done for boundary layer flows. Porteous and Denn 23 studied the linear stability analysis of plane Poiseuille flow for the second order (α 1 < 0) and Maxwell fluids. They showed a destabilization process due to elasticity. At high values of the elasticity number, the stability is qualitatively different from that for Newtonian fluids because it results from the second mode of the Orr-Sommerfeld equation.

Physics of Fluids
Sureshkumar and Beris 26 used an Arnoldi-based orthogonalization algorithm to investigate the linear stability of Poiseuille flow. The models investigated are Upper Convected Maxwell (UCM), Oldroyd-B, and Chilcott-Rallison fluids. The results show that the destabilization caused by elasticity for the UCM fluid is reduced when the effects of solvent viscosity and finite extensibility are taken into account. Zhang et al. 27 showed that when the polymer relaxation time is shorter than the instability time scale, the Poiseuille flow of FENE-P fluids appears to be less stable. However, in the opposite case, the strong elastic effect stabilizes the flow.

III. GOVERNING EQUATIONS
The field equations for an incompressible second order fluid can be derived by substituting expression (1) for the Cauchy stress into the balance of linear momentum, where ρ is the density of the fluid. Since the fluid is incompressible, we require all possible motions be isochoric, and hence, the continuity equation reduces to The geometric configuration considered consists of a wedge of angle γπ, which is placed symmetrically with respect to the direction of the uniform velocity field, as shown in Fig. 1 for γ > 0. The x-axis is chosen to be in the streamwise direction, the z-axis in the spanwise direction, and the y-axis in the wall-normal direction. Due to the symmetric nature of the problem, we can restrict our analysis to the case y ≥ 0. Note that if γ = 0, we recover the case of flow over a semiinfinite flat plate, while γ = 1 corresponds to the case of a stagnation point flow. When γ > 1, we have the flow into an acute corner, γ < 0 gives a flow past a corner, and 0 < γ < 1 is the flow past an acute wedge.

IV. MEAN FLOW
Applying the boundary layer approximation to the field equations (4) and (5) in the same way it is done for Newtonian fluids, we obtain where starred dependent and independent variables indicate dimensional variables. If the plate forms an angle γπ/2 with respect to the uniform velocity field (Fig. 1), the free-stream velocity varies with distance to the leading edge according to potential flow theory 28 as a power law, where a is a positive constant and the exponent m is related to the angle parameter: γ = 2m/(m + 1). After the boundary layer transformation, is a measure for the displacement thickness and ψ * is the stream function, the boundary layer equations (6) are transformed into the following local ordinary differential equation (ODE): where ′ indicates the derivative with respect to the boundary layer variable, η. The key idea is to solve Eq. (9) numerically for fixed values of x * in order to obtain a local solution. It can be easily seen that a similarity solution is possible only for stagnation point flow, 8,9 where m = 1. For the stability analysis, Eq. (9) will be transformed and the dependency on the streamwise position x * will be included in the elasticity parameter, which will be defined later in this section. The stability analysis is traditionally performed, for a Newtonian fluid, by choosing a fixed streamwise position x * = x 0 , as first proposed by Tollmien. 29 The approach consists of finding the longitudinal velocity at that station, ignoring the relatively small transverse velocity, and then solving the Orr-Sommerfeld equation for the resulting base profile.
Following the example of Schmid and Henningson, 30 we apply the same procedure to the second order fluid and we define a displacement thickness, δ 0 , at position x 0 , as follows: where C is a constant given by calculated in the Newtonian case. This choice was made in order to easily compare non-Newtonian solutions with Newtonian solutions. The Reynolds number Re 0 = Ue(x 0 )δ 0 /ν based on the displacement thickness satisfies the following relation: Using relation (12), Eq. (9) at the fixed position x 0 can be rewritten as where K 0 = α 1 /(ρδ 2 0 ) is a non-dimensional parameter known as the elasticity number that can be interpreted as representing the ratio of non-Newtonian normal stress forces to inertial forces. In fact, In order to have a valid boundary layer theory for non-Newtonian fluids of second order, it is necessary that not only the ratio of the inertial forces to the forces due to the tangential stresses be large (high Reynolds number), as in the Newtonian case, but also the ratio of the inertial forces to the forces due to the normal stresses should be large. 3 This implies the following assumptions: The mean flow for the stability analysis is non-dimensionalized using the free-stream velocity Ue(x 0 ), defined by (7), at the streamwise location x 0 . Hence, the velocity in the x-direction is The wall-normal velocity VB is It is clear that this flow is nearly parallel because the transverse velocity VB is smaller than UB by a factor of Re −1 0 , so it will be neglected in order to perform the stability analysis. This is a valid approximation when the Reynolds number Re 0 based on the displacement thickness is large.
Equation (13) is solved by applying the usual boundary conditions that ensure no-slip and no-penetration at the wall and matching with the free-stream velocity at infinity, augmented by the condition The additional boundary condition (15) is required since the differential equation (13) is of fourth order. Condition (15) is derived by imposing ∂u * ∂y * → 0 at infinity and is equivalent to requiring that the solution approaches the free-stream velocity smoothly far from the wall. 8,9 The effect of elasticity on the velocity profile changes with the geometrical configuration. For K 0 > 0, we can see from Figs. 2(a) and 2(b) that the velocity at all points in the boundary layer is larger in the non-Newtonian case for the flow over a flat plate (γ = 0) and the greater variation appears at the wall. Instead, for the second order model with K 0 < 0, the velocity at all points in the boundary layer is smaller in the non-Newtonian case for flow over a flat plate. Figures 2(c) and 2(d) show that, for a wedge angle of π/2, there is a smaller relative variation than for the flat plate observed in Figs. 2(a) and 2(b). When K 0 > 0, the non-Newtonian velocity is slightly smaller inside the boundary layer, while when K 0 < 0, the non-Newtonian velocity is larger. In both cases, the greater deviation from the Newtonian profile happens at a distance η ≈ 2 from the wall. In Figs. 2(e) and 2(f), we see that the effect of increasing |K 0 | for stagnation point flow (γ = 1) is the opposite of the flat plate case. Further investigations on the effects of elasticity on the mean flow are given by Cracco. 31 Note that the non-Newtonian parameter K 0 in Fig. 2 has been chosen to be large enough to be able to distinguish clearly the non-Newtonian effects on the mean flow. However, as already mentioned, we need |K 0 | ≪ 1 for the boundary layer theory to be valid.

V. MODAL ANALYSIS
For the purpose of the stability analysis, we scale the velocities with the constant free-stream velocity Ue(x 0 ) and the lengths with the displacement thickness δ 0 relative to the fixed location x 0 , defined by Eq. (10). The new dimensionless variables are In order to perform a local linear stability analysis, we assume the undisturbed flow to be steady and parallel, neglecting the transverse component of the velocity. The velocity of the base flow in the streamwise direction is taken to be UB(y) given by (14), i.e., the solution of the ODE (13) resulting from the boundary layer approximation at the fixed location x 0 , as shown in Sec. IV.
The governing equations (4) and (5) are linearized around the mean flow UB = (UB, 0, 0) T and written in terms of the disturbance wall-normal velocity v and the disturbance wall-normal vorticity, η = ∂u ∂z − ∂w ∂x . The linear system governing three-dimensional disturbances is obtained after the application of the normal mode form, as follows: where α and β are, respectively, the streamwise (x-direction) and spanwise (z-direction) wavenumbers and ω represents the frequency.
Defining q = (v,η) T , the problem to be solved is a linear system of the form where M and L are linear operators defined as follows: where D denotes the derivative with respect to y and k 2 = α 2 + β 2 .
We can see that in the Newtonian case, when K 0 = 0, the equation forv does not involve the wall-normal vorticityη and reduces to the well-known Orr-Sommerfeld equation when β = 0. Instead, the equation forη, also known as Squire's equation, is driven by solutions to the Orr-Sommerfeld equation through the forcing term βU ′v . In the Newtonian case, this term is responsible for an algebraic growth of energy and is referred to as the vortex tilting term. 32 We observe that, for a non-zero non-Newtonian parameter K 0 and a non-zero spanwise wavenumber β, system (17) is fully coupled due to the presence of a purely non-Newtonian coupling operator, LCN.
The system of equations (17) is subject to the boundary conditionsv = Dv =η = 0 at y = 0 and y → ∞.
The conditions at y = 0 are due to the no-slip and no-penetration at the rigid wall. The conditions at infinity derive from assuming that the disturbances tend to zero far from the surface of the plate. The linear stability equations (17) are solved using a Chebyshev collocation method. 33 The semi-infinite domain y ∈ [0, ∞) is mapped into the finite interval ξ ∈ [−1, 1] by means of the algebraic transformation ξ = y − 2 y + 2 .
All the numerical results are validated in the Newtonian limiting case by comparing with results in the literature. 30,34 A. Growth rates At first, we focus on two-dimensional disturbances where β = 0 and we solve the equation in the first row of system (17), which consists of a modified Orr-Sommerfeld equation. 24 In Fig. 3, the eigenvalues resulting from the linear temporal analysis of the flow over a flat plate (γ = 0) are displayed. In Fig. 3(a), we compare the eigenspectrum for the second order model with K = K 0 C 2 = 0.03 with eigenvalues obtained in the Newtonian case. The choice of Reynolds number Re = Re 0 /C = 580 and wavenumber α * = α/C = 0.179 [C defined by (11)] generates an unstable mode (i.e., ci > 0) in the Newtonian case, known as a Tollmien-Schlichting wave. We can see the stabilizing effect of elasticity that moves the unstable mode into the lower half-plane. Thus, the flow is temporally stable for the second order model for this choice of wavenumber, elasticity, and Reynolds numbers. In Fig. 3(b), we compare the eigenvalues for the model with K = −0.03 with the Newtonian eigenvalues for the same values of Reynolds number and wavenumber. We observe that in this case, elasticity is destabilizing since it pushes the unstable eigenvalues forward into the positive half-plane. We also note that the structure of the rest of the spectrum is different for the non-Newtonian models depending on the sign of K 0 .
Considering the flat plate configuration (γ = 0), Fig. 4(a) shows the temporal growth rate ω * i = ωi/C as a function of α * . We note that when K = 0.01, the maximum growth rate reduces dramatically from ω * i ≈ 1.8 × 10 −3 to about 10 −3 . Instead, when K = −0.01, the maximum growth rate increases to almost 3 × 10 −3 . In general, decreasing K extends the range of positive rates to shorter waves. Figure 4(b) shows the spatial growth rate −α * i as a function of frequency ω * . Again, we observe the marked stabilizing effect of elasticity in terms of growth rate reduction for K = 0.01. We observe that, for K = −0.05, the maximum growth rate increases, but not so dramatically. We also note that for some wavenumbers α * , the growth rate is actually smaller in the non-Newtonian case. The non-Newtonian effects in both models move the maximum to longer waves.
Since the spatial stability properties seem to reflect the temporal characteristics we choose to focus on the temporal problem. 31 Figure 5 shows the temporal growth rates in the Newtonian and non-Newtonian cases for different values of γ. In each case, we observe a reduction in temporal growth rate of the Tollmien-Schlichting waves due to elasticity for K > 0 and an increase in growth rate for K < 0.

B. Neutral stability curves
Temporal neutral stability curves define the region in the Re 0 -α plane where exponentially growing modes exist and where they do not. The region inside the curves represents instability, while the region outside corresponds to stability.
In order to plot neutral stability curves, we need to take into account that both Re 0 and K 0 depend on the location x 0 . If we decide to perform the stability analysis in which the Reynolds number varies depending on the distance x 0 from the leading edge to the location where the local stability analysis is performed, then we need Physics of Fluids ARTICLE scitation.org/journal/phf to write K 0 in terms of the Reynolds number and the base profile needs to be computed for each value of Re 0 . In the flat plate case (γ = 0), the non-Newtonian parameter based on the displacement thickness can be rewritten as Thus, we define the fixed quantity,K = α 1 a 2 ρν 2 , which is independent of x 0 so that K 0 (Re 0 ) =K/Re 2 0 . Figure 6(a) shows a comparison between the neutral stability curves in the Newtonian case and forK = ±10 3 for flow over a flat plate. This clearly shows the stabilizing effect of elasticity in the second order model withK > 0 in terms of increase of the critical Reynolds number. The non-Newtonian effects in the model with K < 0 promote the onset of instabilities. For high Reynolds numbers, the non-Newtonian neutral curves approach the Newtonian neutral curve. This behavior is expected, since when Re 0 → ∞, we have K 0 → 0.
In the case of a non-zero pressure gradient (γ ≠ 0), it is not possible to isolate Re 0 to vary the position x 0 only through the Reynolds number since we have For this reason, we decided to plot the neutral curves in Figs. 6(b)-6(e) by fixing the streamwise position at x 0 = 1. In this case, the interpretation must be different, the Reynolds number varies through a variation of the free-stream velocity Ue. Once again, whenK > 0, elasticity has the effect of reducing the region of two-dimensional instability as shown in Fig. 6 for different angle parameters. Wheñ K < 0, the instability happens at lower Reynolds numbers. Moreover, the neutral curves in the non-Newtonian case approach the Newtonian curves when the Reynolds number increases. It is worth noticing in Fig. 6(e) that, for the flow past a corner (γ = −0.14), as the Reynolds number increases, the non-Newtonian curves overlap the Newtonian curve. This means that the inviscid instability, which arises in the presence of an inflectional velocity profile, does not seem to be affected by non-Newtonian effects. Note that, for different values of γ, different values ofK are chosen in order to ensure that the Weissenberg number, Wi 0 , is of order 1 when the Reynolds number is close to critical for the onset of instability. This is to ensure that the boundary layer theory is valid, while the elasticity effects remain significant. 3

C. Critical Reynolds numbers
The critical Reynolds number is defined as the smallest Reynolds number for which there exists an exponentially unstable mode. We calculated the critical wavenumbers, αcr, and Reynolds numbers, Recr, for different values of γ, and the results are displayed in Table I. In order to be able to compare the non-Newtonian effect of elasticity for different values of γ we choose, as a measure of elasticity, the critical Weissenberg number, Wi 0,cr = K 0,cr Re 0,cr , defined with reference to the Newtonian critical Reynolds number Re 0,cr and the critical elasticity number K 0,cr =K/Re 2 0,cr . From Table I, we observe that the critical Reynolds numbers for Wi 0,cr > 0 increase for all values of γ considered, including the slightly negative value of γ that represents a profile with an inflection point. The effect is the opposite for Wi 0,cr < 0, where the critical Reynolds numbers decrease and the instability is anticipated for each value of the angle parameter γ.
Note that the magnitude of the critical Reynolds number Re 0,cr for the Newtonian case is strongly dependent upon the configuration characterized by γ. This strong dependence is maintained for the variation found in Re 0,cr when the non-Newtonian effects are introduced in the manner that we have described. For example, with a critical Weissenberg number Wi 0,cr = 0.5, for a flat plate (γ = 0), the increase or decrease in critical Reynolds number is of order 10, while for the stagnation point flow (γ = 1), it is of order 10 2 .

D. 3D results
We solved the three-dimensional eigenvalue problem (17). The results obtained are summarized by displaying the neutral stability curves in the α-β plane. A study of three-dimensional disturbances for fluids of second order is required. For parallel Newtonian flow, Squire's theorem justifies the study of two-dimensional instead of three-dimensional disturbances. Squire's theorem states that each three-dimensional mode corresponds to some two-dimensional mode at a lower Reynolds number. Therefore, to determine the critical Reynolds number, it is sufficient to study two-dimensional disturbances for Newtonian fluids. A result similar to Squire's theorem for a fluid of second order cannot be proven. Therefore, an extension to the study of three-dimensional disturbances is necessary. Figure 7 shows the contour plot of the temporal growth rate ωi in the Newtonian case for the flat plate (γ = 0). Figure 7(a) shows that the choice of a subcritical Reynolds number (Re 0 = 500) gives a stable flow. In Fig. 7(b), we increase the Reynolds number to Re 0 = 1000, and we can see an exponential instability for which ωi > 0, appearing at small spanwise wavenumbers. The red asterisk ( * ) represents the maximum growth rate reached in the α-β plane, i.e., maxα, β ωi. We can see that, in both cases, the maximum is reached for spanwise independent waves. This confirms Squire's theorem for Newtonian fluids. Figures 8(a) and 8(b) show the contour plots of the temporal growth rates ωi for the model with K < 0 and for the model with K > 0, respectively. We can see that when K > 0, there is a region of exponential instability for small streamwise wavenumbers and for a value of the Reynolds number (Re 0 = 500) that corresponds to a stable flow in the Newtonian case. In Fig. 8(c), we displayed the growth rates for a fixed and small α and for a fixed value of β in Fig. 8(d).
We observe how a positive elasticity number K destabilizes spanwise disturbances, while it stabilizes the Tollmien-Schlichting waves. The opposite happens for a negative value of K, which decreases the growth rates of mainly streamwise independent waves (α ≈ 0) and increases the growth rates of mainly spanwise independent waves (β ≈ 0). Figure 9 shows the results for a Reynolds number of Re 0 = 1000. The conclusions are the same, for K > 0, the Tollmien-Schlichting wave is slightly stabilized while growth rates near the α = 0 axis become larger. The opposite happens for K < 0. Figure 10 shows growth rates for the flow past a corner with γ = −0.14. The results are very similar to those for the flat plate. We do not report results for other values of the angle parameter γ since they are in line with the results we have discussed so far.

Physics of Fluids
To the best of our knowledge, the differing effects of elasticity on two-dimensional and three-dimensional disturbances have not been observed in the past for other viscoelastic models.

VI. NON-MODAL ANALYSIS
In this section, the initial-value problem that drives the development of disturbances is derived for the second order fluid. A formulation based on the initial-value problem enables us to study the behavior of general solutions not only of single eigenmodes. 30 Brandt 19 reviewed the main results in bypass transition for non-Newtonian fluids. Zhang et al. 27 performed the modal and non-modal linear analysis of the inertia-dominated channel flow of viscoelastic fluids modeled by Oldroyd-B and FENE-P closures. The authors observed that both modal and non-modal instabilities are enhanced when the polymer relaxation time is shorter than the instability timescale (i.e., for Weissenberg numbers, Wi ≲ 1), whereas the flow is more stable in the opposite case. In the subcritical regime, the non-modal amplification of streamwise elongated structures is still the most dangerous energy growth mechanism and is slightly enhanced by the presence of polymers. The lift-up effect is still the dominant instability mechanism also for viscoelastic fluids. 27 Hoda et al. 35 studied energy amplification in channel flows of Oldroyd-B fluids from an input-output point of view and found that increasing fluid elasticity through polymer contribution to the viscosity or the elasticity number enhances energy amplification. Once again, the disturbances that are most amplified are streamwiseelongated, with elasticity acting to reduce spanwise length scale.
Some other authors focused on weakly inertial or inertialess flows (e.g., Lieu et al. 36 and Page and Zaki 37 ), while more recently, others studied the secondary instability of streaks for viscoelastic fluids (e.g., Burshtein et al. 38 ).
After linearization around the mean flow UB = (UB, 0, 0) T , we take the normal mode form for the perturbations. However, unlike in Sec. V, we do not assume an exponential time-dependence. The initial value problem can be written as follows: where q(t, y) = (v(t, y),η(t, y)) T . The linear operators M, L are defined by Eq. (18). General solutions of the initial-value problem (19) are assumed to belong to the space S N spanned by a sufficient number N of eigenfunctions, which is defined as follows: where {q j }j are solutions of the eigenvalue problem (17). In other words, q ∈ S N can be expressed as where {kj}j are the coefficients of the expansion. This allows us to express the initial value problem (19) as N separated ordinary differential equations for the expansion coefficients, as follows: or in a more compact form, i.e., where k = (k 1 , . . . , kN) T and Ω = diag{ω 1 , . . ., ωN}. The simplified formulation (21) of the initial-value problem (19) is possible provided that the eigenspectrum is a complete set composed of discrete eigenmodes. For Newtonian fluids, it is known that if the domain is bounded, then the eigenspectrum is discrete, but for unbounded boundary layers, the spectrum is composed of a discrete and a continuous part. Although the discrete approximation differs from the exact representation, the sum of these eigenmodes correctly describes the solutions to the initial-value problem. 17 For Newtonian fluids, the completeness of the spectrum is proven by Gustavsson. 39 To the best of our knowledge, the completeness of the spectrum has not been proven yet for second order fluids or non-Newtonian fluids in general. We will not research this further in this paper, and we discretize the continuous spectrum for the second order models, as done by Butler and Farrell. 17 Therefore, particular attention is paid to ensure that the results are independent of the discretization parameter.
In order to determine the perturbation that grows the most in some sense, we need a way to quantify the growth. The energy norm is taken to be where M is defined by (18a). In order to quantify the transient growth, we define the maximum possible amplification of initial energy density as follows: where L 1 = −iM −1 L and L, M are the linear operators given by (18). Fixing the wavenumber vector (α, β), the function G represents the envelope of the energy evolution of all the initial perturbations, q 0 , with unit energy norm. At each moment in time, we maximize over all possible initial conditions. Note that traditional stability analysis focuses attention only on the eigenvalues of e −iΩt . These do not capture the whole behavior of G, which is also determined by the eigenvector matrix F and its inverse. Deducing the behavior of G from the eigenvalue matrix Ω alone is only valid when the similarity transformation given by F does not alter the norm, that is, when V is unitary and composed of orthogonal eigenvectors. This is the case when B is normal. If this is not the case, B is non-normal and short-time growth of perturbation energy is possible even though the matrix has stable eigenvalues. In order to compute the exponential norm (23), we use the decomposition (20). Thus, G can be calculated easily as follows: where σ 1 is the principal singular value of the matrix B = Fe −iΩt F −1 .
Employing the decomposition (20) provides an easy way to compute the maximum possible amplification G, which can be obtained by calculating the singular value decomposition (SVD) of the matrix B.

A. Global optima
We define the global optimal disturbance as the initial condition, q 0 , that maximizes the growth over time, i.e., Note that Gmax can only be defined when all the eigenvalues are stable. If an unstable mode exists, then G(t) → ∞ as t → ∞. We can also define the largest global growth obtained for any wavenumber vector as follows: The latter depends only on the base flow conditions and Re.
The results obtained have been validated by comparing with those found in the literature for Newtonian fluids. For this purpose, we refer to the book by Schmid and Henningson 30 and the paper by Corbett and Bottaro. 40 Figure 11 shows the contour plot of Gmax defined by (24) for the flat plate (γ = 0). The black line represents the neutral stability curve inside which an exponentially growing mode exists and where the maximum possible amplification is not defined or can be thought of as infinite. The Newtonian results in Fig. 11(a) are in agreement with the literature. 30,41 The largest global optimal growth defined by (25) is G Γ = 1515.6 reached at time t = 782 for α Γ = 0, β Γ = 0.65, as calculated by Corbett and Bottaro. 40 Figures 11(b) and 11(c) show the contour plot for K = 10 −4 and K = −10 −4 , respectively. These non-Newtonian parameters have been chosen as an example to show the non-Newtonian effects. We can see that the largest amplification of energy is still reached for streamwise independent disturbances, as in the Newtonian case. However, when K > 0, the amplification of energy is generally larger, and when K < 0, the amplification of energy is smaller than in the Newtonian case. Figure 12 shows the contour plot of Gmax for the flow past a wedge (γ = 0.5). The non-Newtonian effects on the transient growth are qualitatively similar to the flat plate case. Figure 13(a) displays the ratio of non-Newtonian Gmax to Newtonian Gmax for a fixed spanwise wavenumber β = 0.6 and varying Weissenberg number Wi 0 . We can observe that the non-Newtonian terms mostly affect streamwise independent disturbances, i.e., for α = 0. In Fig. 13(b), we can see that for K > 0, the global optima happen at larger times than in the Newtonian case, while for K < 0, the global optima happen at shorter times.
In Table II, we report the largest global optima G Γ defined in (25). For these calculations, we choose the momentum thickness scaling, following Corbett and Bottaro. 40 The reason is that when scaled using the momentum thickness, the spanwise wavenumber at which the largest global optimum is reached is independent of the mean flow conditions. Moreover, momentum thickness scaling accounts for the variation in t Γ (the time in which the optimal disturbance reaches its maximum) resulting from differences in the base flow. We choose to scale the lengths with the momentum thickness θ 0 relative to the fixed streamwise location x 0 , which is defined as follows: where δ is defined by Eq. (8) and θ Newt is the constant, calculated in the Newtonian case. We introduce Reynolds and Weissenberg numbers based on θ 0 as follows: Note that the following relations hold: where H = C/θ Newt is the shape factor defined as the ratio between displacement and momentum thickness, calculated in the Newtonian case. In Table II, we present the results obtained for Reynolds numbers Re θ = 166 and Re θ = 385. These Reynolds numbers have been chosen to compare the results with the ones obtained by Corbett and Bottaro. 40 Specifically, Re θ = 385 corresponds to the Reynolds number based on the displacement thickness Re 0 ≈ 1000 for the flat plate case. For all the flows considered, the largest global optimum is reached for streamwise-independent waves, i.e., α Γ = 0. We can see that, in the Newtonian case, when scaled with θ 0 , the spanwise wavenumber for G Γ appears to be independent of the mean flow condition characterized by γ and β θ ≈ 1/4. Note that, in the Newtonian case, the moment in time at which the largest global optimum is reached is about the same for all the positive angle parameters considered, t θ ≈ 880.
For flow past a corner (γ = −0.14), the maximum is reached at a larger time t θ ≈ 927. We observe how, for all the angle parameters

ARTICLE
scitation.org/journal/phf considered, the spanwise wavenumber β θ , the time t θ , and the largest possible amplification G Γ decrease when the model with Wi θ < 0 is selected and increase when the model with Wi θ > 0 is selected. Moreover, β θ appears to change approximately linearly with the Weissenberg number based on the momentum thickness. A Weissenberg number Wi 0 = ±0.05 produces a change in β θ of about 1% and Wi 0 = ±0.1 produces a change of about 2%. This linear dependence on the Weissenberg number manifests also on the time t θ and on the largest transient growth G Γ .

B. Optimal perturbations
We can determine the initial condition that reaches the maximum possible amplification at a given time t 0 by using the singular value decomposition (SVD) of the matrix B = Fe −it 0 Ω F −1 . The initial condition that reaches the global optimum Gmax at t = tmax defined by (24) is referred to as optimal disturbance. Figure 14 shows a comparison between optimal disturbances in the Newtonian and non-Newtonian cases for the stagnation point flow (γ = 1) and a Reynolds number Re 0 = 500. We choose a wavenumber vector (α, β) = (0, 0.6), which is close to the global optima. In Figs. 14(a)-14(c), |u| has been scaled such that max(|v 0,Newt |) = 1, and in Figs. 14(b)-14(d), |u| has been scaled such that max(|v max,Newt |) = 1.
We see that the optimal disturbances, in the non-Newtonian cases, have the same structure of streamwise-oriented vortices as in the Newtonian case. From Figs. 14(a)-14(c), we observe that the initial streamwise velocity |u 0 | is two orders of magnitude smaller than the cross-flow components, |v 0 |, |w 0 |. Figures 14(b)-14(d) show the evolved state of the optimal disturbances at t = tmax. The shape of the initial vortex is still present, although it has diffused outward away from the wall.
At t = tmax, the streamwise velocity |umax| is one order of magnitude larger than the cross-flow velocities, which indicates the presence of streaks. From Figs. 14(a) and 14(b), we see that when K > 0, the vortices are more diffused away from the wall, whereas when K < 0, the vortices are closer to the wall. Figures 14(c) and 14(d) show that, for the non-Newtonian fluid with K = −0.0001, the initial optimal streamwise velocity is larger than in the Newtonian case and at tmax, it grows more than in the Newtonian case. The behavior is the opposite when K = 0.0001. This is in agreement with the results obtained in Sec. VI A.
In Fig. 15, we plotted the streamwise vortices for K = −0.0001 and Re 0 = 1000. The solutions plotted are such that ∥q 0 ∥E = 1 and

VII. CONCLUSIONS
The linear stability analysis of the boundary layer flow of a viscoelastic fluid has been investigated. We applied a boundary layer theory to second order fluids in order to determine the mean flow. As for Newtonian fluids, this approach allowed us to simplify the governing equations. There is no similarity solution to the boundary layer equations with the exception of stagnation point flow. Therefore, a local similarity transformation was applied that yielded a local ODE depending on the streamwise position. 9 We reformulated this equation in order to represent the dependency on the location only through the non-Newtonian parameter K based on the displacement thickness at a fixed streamwise position. The equation obtained is then solved numerically and employed in the linear stability analysis. This approach is consistent with the fact that traditional linear stability analysis is a local analysis.
First, we considered two-dimensional disturbances, namely, disturbances that vary only in the streamwise and wall-normal directions. We solved numerically, using a Chebyshev collocation method, the modified Orr-Sommerfeld equation that governs the evolution of two-dimensional disturbances. The results were presented in terms of growth rates and neutral curves. For all values of the angle parameter γ, the non-Newtonian terms in the second order model with K > 0 stabilize the flow with respect to the Newtonian case, while they have the opposite effect for K < 0. Moreover, we determined the critical Reynolds number, which is the smallest Reynolds number for which there exists an exponentially unstable mode. For K > 0, there is a stabilizing effect in terms of an increase of the critical Reynolds numbers. The effect is the opposite for K < 0, where the instability is enhanced. The linear stability results for K < 0, which is the one that predicts the correct sign of the non-zero normal stress differences, are in qualitative agreement with those obtained by Sureshkumar and Beris 26 and Zhang 27 for the Poiseuille flow of other viscoelastic fluids.
When extending the analysis to three-dimensional disturbances, which can also vary in the spanwise direction, the non-Newtonian effects prove to be different. This differing behavior has not been observed before for other viscoelastic models. We showed that a positive elasticity number K destabilizes spanwise disturbances while it stabilizes the two-dimensional Tollmien-Schlichting waves. The opposite happens for a negative K, which decreases the growth rates of mainly streamwise independent waves and increases the growth rates of mainly spanwise independent waves.

ARTICLE scitation.org/journal/phf
In order to give a complete idea of the linear stability characteristics, the potential transient growth of energy cannot be ignored. To the best of our knowledge, the transient growth of viscoelastic fluids in boundary layers has not been investigated in the past. In this paper, the initial-value problem that drives the development of disturbances is derived for second order fluids. In the Newtonian case, our results are validated against those obtained by Schmid 41 for Blasius flow and by Corbett and Bottaro 40 for Falkner-Skan flows. We showed that, for the second order model with K > 0, an increase in the non-Newtonian parameter K provokes an increase in the maximum transient growth, G, while the second order model with K < 0 has the opposite behavior. The results are qualitatively similar for all values of the angle parameter, γ. The largest amplification of energy is still reached for streamwise independent disturbances (zero streamwise wavenumber), as in the Newtonian case.
Non-Newtonian terms mostly affect streamwise independent disturbances. For K > 0, the global optimum, Gmax, is reached for larger times and shorter waves (larger spanwise wavenumber) than in the Newtonian case. On the contrary, for K < 0, the global optimum is reached for shorter times and longer waves (smaller spanwise wavenumber).
The second order models are not used in industrial applications, where other constitutive equations are preferred. However, they can represent non-zero normal stress differences, and they give an idea of the stability characteristics of viscoelastic fluids in boundary layers. The natural progression of this work is the investigation of the linear stability properties of rheologically more complex models in boundary layers. Part of this work was performed using the computational facilities of the Advanced Research Computing@Cardiff (ARCCA) Division, Cardiff University.

DATA AVAILABILITY
The data that support the findings of this study are openly available via the Cardiff University data repository at https://doi.org/10.17035/d.2020.0112209117.