2 Mathematical formulation

Consider an incompressible two-dimensional viscoelastic liquid film flow with density $\unicode[STIX]{x1D70C}$ , limiting viscosity $\unicode[STIX]{x1D707}$ and surface tension $\unicode[STIX]{x1D70E}$ on a horizontal plane, where the motion of the viscoelastic liquid starts due to an oscillation of the plane in the horizontal direction with a velocity $\hat{\boldsymbol{i}}U_{0}\cos \unicode[STIX]{x1D714}t$ , where $\hat{\boldsymbol{i}}$ is a unit vector in the streamwise direction. Figure 1 shows a schematic diagram of the present model, where the origin is located at the free surface, and the $x$ and $y$ axes are directed along the streamwise flow and the cross-streamwise flow respectively. Walter’s liquid $B^{\prime \prime }$ , which is an approximation to the first order in elasticity from the Newtonian behaviour and possesses a short or rapidly fading memory, is considered in the current study. The liquid $B^{\prime \prime }$ satisfies the following constitutive equation of state (Beard & Walters Reference Beard and Walters1964; Andersson & Dahl Reference Andersson and Dahl1999):

where $\unicode[STIX]{x1D70F}_{ij}$ is the stress tensor, $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta, $p$ is the pressure and $\unicode[STIX]{x1D70F}_{ij}^{\prime }$ is defined by

where $2\unicode[STIX]{x1D707}\unicode[STIX]{x1D626}_{ij}$ is the Newtonian viscous stress contribution, $2M_{0}(\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}t)\unicode[STIX]{x1D626}_{ij}$ is the elastic stress contribution, $M_{0}$ is the viscoelastic coefficient, $\unicode[STIX]{x1D626}_{ij}=(\unicode[STIX]{x2202}_{i}u_{j}+\unicode[STIX]{x2202}_{j}u_{i})/2$ is the strain rate tensor and the co-rotational derivative of the strain rate tensor $\unicode[STIX]{x1D626}_{ij}$ is defined by

Obviously, the constitutive equation (2.1) recovers the Newtonian behaviour when the viscoelastic coefficient $M_{0}$ is set to zero. The general form of the constitutive equation can be found in Oldroyd (Reference Oldroyd1950), Beard & Walters (Reference Beard and Walters1964), Bird et al. (Reference Bird, Curtiss, Armstrong and Hassager1987) and Barnes, Hutton & Walters (Reference Barnes, Hutton and Walters1989). The rheological model given in the constitutive equation (2.1) is obtained by removing the relaxation time parameter, and the reduced model is referred to as Walter’s liquid $B^{\prime \prime }$ (Beard & Walters Reference Beard and Walters1964; Andersson & Dahl Reference Andersson and Dahl1999). A mixture of polymethyl methacrylate in pyridine with density $\unicode[STIX]{x1D70C}=0.98\times 10^{3}~\text{kg}~\text{m}^{-3}$ , limiting viscosity $\unicode[STIX]{x1D707}=0.79~\text{N}~\text{s}~\text{m}^{-2}$ and viscoelastic coefficient $M_{0}=0.04~\text{N}~\text{s}^{2}~\text{m}^{-2}$ is an example of Walter’s liquid $B^{\prime \prime }$ (Walters Reference Walters1960; Andersson & Dahl Reference Andersson and Dahl1999).

Figure 1. Schematic diagram of a viscoelastic liquid on an oscillating plane.

The mass and momentum conservation equations that govern the flow are

where $g_{x}=0$ and $g_{y}=-g$ , $g$ is the gravitational acceleration. The associated boundary conditions are as follows. (i) At the plane $y=-d$ ,

(ii) At the free surface $y=h(x,t)$ , the balances of tangential and normal stresses,

where $\boldsymbol{n}=(n_{x}\hat{\boldsymbol{i}}+n_{y}\hat{\boldsymbol{j}})$ is the unit normal vector, $\boldsymbol{t}=(t_{x}\hat{\boldsymbol{i}}+t_{y}\hat{\boldsymbol{j}})$ is the unit tangent vector and $\hat{\boldsymbol{j}}$ is a unit vector in the cross-streamwise direction. (iii) The kinematic boundary condition at the free surface $y=h(x,t)$ is given by

where $f(x,t)=h(x,t)-y$ . Since the present study is an unsteady flow problem and the primary flow occurs due to the horizontal oscillation of the plane, an appropriate scale for time is the inverse of imposed frequency. Next, the velocity scale is obtained from the balance of the inertia term and the viscous friction term. Accordingly, the governing equations (2.4)–(2.8) are normalized by choosing $d$ as the length scale, $\unicode[STIX]{x1D708}/d$ as the velocity scale, $1/\unicode[STIX]{x1D714}$ as the time scale and $\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}^{2}/d^{2}$ as the pressure scale. The proposed characteristic scales lead to a set of dimensionless parameters: the Reynolds number, $Re=U_{0}d/\unicode[STIX]{x1D708}$ , shows the effect of the oscillation amplitude, the Galileo number, $\unicode[STIX]{x1D712}=gd^{3}/(2\unicode[STIX]{x1D708}^{2})$ , shows the effect of gravity, the capillary number, $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}^{2}/(\unicode[STIX]{x1D70E}d)=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D708}/d)/\unicode[STIX]{x1D70E}$ , shows the effect of surface tension, the viscoelastic parameter, $M=M_{0}/(\unicode[STIX]{x1D70C}d^{2})$ , shows the effect of viscoelasticity and the Womersley number, $\unicode[STIX]{x1D6FD}=\sqrt{\unicode[STIX]{x1D714}d^{2}/(2\unicode[STIX]{x1D708})}$ , shows the effect of the imposed frequency. The dimensionless viscoelastic parameter $M$ compares the square of the elastic length $M_{0}/\unicode[STIX]{x1D70C}$ with the square of the mean film thickness $d^{2}$ , while the dimensionless imposed frequency $2\unicode[STIX]{x1D6FD}^{2}$ compares the square of the mean film thickness $d^{2}$ with the square of the Stokes layer thickness $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}$ , or, equivalently, it compares the viscous time scale $d^{2}/\unicode[STIX]{x1D708}$ with the frequency time scale $1/\unicode[STIX]{x1D714}$ .

Consider a unidirectional time periodic parallel flow with constant film thickness. Consequently, the governing equations (2.4)–(2.8) are simplified to the following dimensionless form:

The base flow solution can be expressed as

where $S=[(1+4M^{2}\unicode[STIX]{x1D6FD}^{4})^{1/2}+2M\unicode[STIX]{x1D6FD}^{2}]$ and $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6FD}\![\{(1+4M^{2}\unicode[STIX]{x1D6FD}^{4})^{1/2}-2M\unicode[STIX]{x1D6FD}^{2}\}/(1+4M^{2}\unicode[STIX]{x1D6FD}^{4})\!]^{1/2}$ . Here, $\text{Re}[\cdots \,]$ represents the real part of that complex function. It should be noted that the base flow is strongly dependent on the viscoelastic parameter $M$ . In the absence of the viscoelastic parameter ( $M=0$ ), the base flow velocity is identical to that of Or (Reference Or1997). Moreover, the base velocity coincides with that of Dandapat & Gupta (Reference Dandapat and Gupta1975) when the characteristic scales are the same.

3 Orr–Sommerfeld boundary value problem

The primary instability is studied by considering a two-dimensional infinitesimal disturbance $q=Q+\tilde{q}$ on the base flow, where $q$ represents the flow variables, $Q$ represents the basic quantities and $\tilde{q}$ represents the perturbation quantities. By substituting $q$ in the dimensionless form of the governing equations (2.4)–(2.8) and linearizing with respect to the base state, we obtain the following linearized perturbation equations:

Now, the stream function $\tilde{\unicode[STIX]{x1D713}}$ is introduced by using the relations $\tilde{u} =\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D713}}$ and $\tilde{v}=-\unicode[STIX]{x2202}_{x}\tilde{\unicode[STIX]{x1D713}}$ . Next, the solution of the perturbation equations (3.1)–(3.7) is assumed in the normal mode form $\tilde{\unicode[STIX]{x1D713}}(x,y,t)=\unicode[STIX]{x1D719}(y,t)\exp (\text{i}kx)$ and $\tilde{h}(x,t)=\unicode[STIX]{x1D702}(t)\exp (\text{i}kx)$ , where $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D702}$ are respectively the amplitude of the perturbation stream function and the deformed free surface. Here, $k$ is the wavenumber. By inserting $\tilde{\unicode[STIX]{x1D713}}$ and $\tilde{h}$ into the perturbation equations (3.1)–(3.7) and eliminating the pressure term from the momentum equations, we obtain the following time-dependent Orr–Sommerfeld boundary value problem (OS BVP) for a viscoelastic liquid:

where ${\mathcal{D}}=\unicode[STIX]{x2202}_{y}$ and ${\mathcal{L}}=\unicode[STIX]{x2202}_{yy}-k^{2}$ are differential operators. When the viscoelastic parameter $M$ vanishes, the present OS BVP (3.8)–(3.12) reduces to the time-dependent OS BVP for a Newtonian liquid (Or Reference Or1997). The difference in the coefficients from the OS BVP obtained by Dandapat & Gupta (Reference Dandapat and Gupta1975) is attributed to the choice of various characteristic scales.

4 Long-wavelength expansion

The time-dependent OS BVP (3.8)–(3.12) with periodic coefficients is resolved based on the long-wavelength expansion ( $k\rightarrow 0$ ) along with the Floquet theory proposed by Yih (Reference Yih1968). Consequently, the solution is assumed in the following form:

where the Floquet exponent $\unicode[STIX]{x1D6FF}$ is of the form

To take into account the effect of surface tension in the first-order formulation, we consider the capillary number $\unicode[STIX]{x1D701}\sim O(k^{2})$ . By inserting (4.1) and (4.2) into the OS BVP (3.8)–(3.12) and collecting the leading-order $O(k^{0})$ terms, we obtain the following set of equations:

Since $\unicode[STIX]{x1D702}_{0}(t)$ is a periodic function of $t$ , the kinematic boundary condition (4.7) leads to an admissible solution $\unicode[STIX]{x1D6FF}_{0}=0$ . Without loss of generality, we choose $\unicode[STIX]{x1D702}_{0}=1$ . Otherwise, $\unicode[STIX]{x1D6FF}_{0}\neq 0$ supplies a damped Floquet mode demonstrated by Yih (Reference Yih1968), and that damped mode is not of interest here. By substituting the values of $\unicode[STIX]{x1D6FF}_{0}$ and $\unicode[STIX]{x1D702}_{0}$ and solving (4.3)–(4.6), the leading-order solution can be expressed as

Obviously, $\unicode[STIX]{x1D719}_{0}$ relies on the viscoelastic parameter $M$ . In the limit $M\rightarrow 0$ , the reduced solution (4.8) is identical to that of Yih (Reference Yih1968) when the characteristic scales are the same. For the first-order $O(k)$ approximation, we have the following set of equations:

Since $\unicode[STIX]{x1D719}_{0}(y,t)$ , $U(y,t)$ and $\unicode[STIX]{x1D702}_{0}(t)$ are all periodic functions of time, to obtain a periodic solution $\unicode[STIX]{x1D702}_{1}(t)$ from the kinematic boundary condition (4.13), we must have $\unicode[STIX]{x1D6FF}_{1}=0$ . As a result, (4.13) reduces to the following form:

which yields a periodic solution

where $\text{Im}[\cdots \,]$ represents the imaginary part of that complex function. Indeed, the first-order equations (4.9)–(4.12) possess steady and unsteady terms. As the Floquet exponent is independent of time, the solution of only first-order steady equations is sufficient to compute the next-order Floquet exponent $\unicode[STIX]{x1D6FF}_{2}$ . Therefore, we shall focus solely on the steady equations in the subsequent calculations. After some mathematical manipulation, the solution of the first-order steady equations can be expressed as

where superscript ‘ ${\mathcal{S}}$ ’ denotes the steady solution and

where $R_{1}$ and $R_{2}$ are of the form

Now, in order to compute the Floquet exponent $\unicode[STIX]{x1D6FF}_{2}$ , we consider the second-order $O(k^{2})$ approximation of the kinematic boundary condition at $y=0$ ,

Again, $\unicode[STIX]{x1D6FF}_{2}$ is independent of time, so we must have

Inserting the values of $\unicode[STIX]{x1D6FF}_{0}$ , $\unicode[STIX]{x1D6FF}_{1}$ and $\unicode[STIX]{x1D6FF}_{2}$ , the Floquet exponent $\unicode[STIX]{x1D6FF}$ can be written as

where VECT, GT and ST represent the viscoelastic coupling term, gravity term and surface tension term respectively. Equation (4.25) demonstrates that the Floquet exponent $\unicode[STIX]{x1D6FF}$ is simply a combination of three terms $VECT,GT$ and $ST$ , and its magnitude will reduce in the presence of gravity and surface tension. In other words, gravity and surface tension exhibit a stabilizing effect on the long-wave mode.

We have noticed that the leading- and first-order solutions of the kinematic boundary condition (3.12) yield $\unicode[STIX]{x1D6FF}_{0}=0$ and $\unicode[STIX]{x1D6FF}_{1}=0$ . Hence, in the long-wavelength limit $k\rightarrow 0$ , the amplitude [ $\propto \text{e}^{\unicode[STIX]{x1D6FF}t}=\text{e}^{k^{2}\unicode[STIX]{x1D6FF}_{2}+O(k^{3})}$ ] of an infinitesimal disturbance will grow or decay exponentially with time if $\unicode[STIX]{x1D6FF}_{2}>\text{or}<0$ , or, equivalently, if the following criterion is satisfied:

Therefore, one can conclude that the long-wavelength disturbance will be stable if $\unicode[STIX]{x1D6FF}_{2}<0$ ; otherwise, it will be unstable if $\unicode[STIX]{x1D6FF}_{2}>0$ . In the absence of elasticity and for a low magnitude of surface tension, equation (4.26) coincides with that of Yih (Reference Yih1968) and Or (Reference Or1997), where $2\unicode[STIX]{x1D712}/Re^{2}=Fr^{-2}=gd/U_{0}^{2}$ , and $L(\unicode[STIX]{x1D6FD})=[3A^{2}(1+\unicode[STIX]{x1D6FE}^{2})/\unicode[STIX]{x1D6FD}^{2}](R_{2}-R_{1})$ in Yih’s notation. In addition, the expression (4.26) is identical to that of Dandapat & Gupta (Reference Dandapat and Gupta1975) when the characteristic scales are the same.

Figure 2. (a) The neutral curves in the $(\unicode[STIX]{x1D6FD},~Re)$ plane for different values of the viscoelastic parameter when $\unicode[STIX]{x1D712}=1$ and $\unicode[STIX]{x1D701}^{\ast }=0$ . Here, the thick solid, dashed and thin solid lines stand for $M=0.01$ , $M=0.005$ and $M=0.0$ respectively. (b) The variation of $L(\unicode[STIX]{x1D6FD})$ as a function of $\unicode[STIX]{x1D6FD}$ when $\unicode[STIX]{x1D712}=1$ and $M=0.01$ . Here, the thick solid, dashed and thin solid lines stand for $\unicode[STIX]{x1D701}^{\ast }=0$ , $\unicode[STIX]{x1D701}^{\ast }=0.5$ and $\unicode[STIX]{x1D701}^{\ast }=1$ respectively. The points stand for the result of Dandapat & Gupta (Reference Dandapat and Gupta1975). (c) The variation of $(VECT+GT+ST)$ as a function of $\unicode[STIX]{x1D6FD}$ when $Re=10$ , $\unicode[STIX]{x1D712}=1$ and $\unicode[STIX]{x1D701}^{\ast }=1$ . Here, the thick solid, dashed and thin solid lines stand for $M=0.01$ , $M=0.005$ and $M=0$ respectively.

In order to analyse the effect of surface tension on the long-wave mode independently, a new parameter is introduced: $\unicode[STIX]{x1D701}^{\ast }=k^{2}/\unicode[STIX]{x1D701}$ represents the effect of surface tension. The neutral stability criterion $\unicode[STIX]{x1D6FF}_{2}=0$ leads to an analytical expression of the Reynolds number as a function of the dimensionless imposed frequency $2\unicode[STIX]{x1D6FD}^{2}$ ,

Figure 2(a) shows that $\text{U}$ -shaped neutral curves appear in the separated bandwidths of the imposed frequency when $k\rightarrow 0$ . Therefore, the long-wavelength instability prevails only in these frequency bandwidths. Outside of these bandwidths, all infinitesimal disturbances are stable. This result is completely consistent with that of Dandapat & Gupta (Reference Dandapat and Gupta1975). For this reason, the stabilizing effect of the viscoelastic parameter was found by Dandapat & Gupta (Reference Dandapat and Gupta1975) in these stable bandwidths of the imposed frequency. Moreover, in each unstable frequency bandwidth, there exists a critical Reynolds number $(Re=U_{0}d/\unicode[STIX]{x1D708})$ , or, equivalently, a critical amplitude of the horizontal oscillation above which the long-wave modes are unstable. The interesting result is that the number of unstable bandwidths decreases with increasing value of the viscoelastic parameter. However, in the different unstable frequency bandwidths, the dominant mode of instability intensifies with increasing value of the viscoelastic parameter. Figure 2(b) demonstrates the variation of $10\times L(\unicode[STIX]{x1D6FD})$ with $\unicode[STIX]{x1D6FD}$ when $M=0.01$ is fixed. If the surface tension is incorporated into the $O(k)$ approximation of the long-wavelength analysis, $L(\unicode[STIX]{x1D6FD})$ , or, equivalently, the growth rate of the long-wave mode, diminishes with increasing value of $\unicode[STIX]{x1D701}^{\ast }$ . Therefore, one can conclude that the viscoelastic parameter has a destabilizing influence and the surface tension has a stabilizing influence on the long-wave mode in the unstable frequency bandwidths for an oscillatory viscoelastic film flow. The stabilizing effect of surface tension also supports the analytical expression given in (4.25). The variation of $10\times L(\unicode[STIX]{x1D6FD})$ with $\unicode[STIX]{x1D701}^{\ast }$ in the higher-frequency zone ( $\unicode[STIX]{x1D6FD}\geqslant 2$ ) can be found as an inset in figure 2(b). Further, the present result captures the result of Dandapat & Gupta (Reference Dandapat and Gupta1975) very well when the surface tension parameter $\unicode[STIX]{x1D701}^{\ast }$ is set to zero. In order to provide further evidence of the destabilizing effect of viscoelasticity on the long-wave mode, the sum ( $VECT+GT+ST$ ) of the different terms in the expression of the Floquet exponent (4.25) is depicted in figure 2(c) for different values of $M$ when the Galileo number $\unicode[STIX]{x1D712}$ and the modified capillary number $\unicode[STIX]{x1D701}^{\ast }$ are fixed. Indeed, the sum is enhanced with increasing value of the viscoelastic parameter, or, equivalently, the magnitude of the viscoelastic coupling term increases with $M$ because both $GT$ and $ST$ are constants. This fact results in growth of the Floquet exponent and supports the destabilizing effect of $M$ . In other words, the destabilizing effect of the viscoelastic parameter on the long-wave mode can also be described somewhat in the spirit of Smith (Reference Smith1990), Huang & Khomami (Reference Huang and Khomami2001) and Wei (Reference Wei2005). If an infinitesimal disturbance is considered on the unsteady base flow, a perturbation shear stress, which is a combination of the Newtonian viscous stress and the elastic stress (see (3.5)), evolves at the free surface to maintain the stress-free surface. As a result, the energy, which is greater in comparison with the Newtonian case, is transferred from the base flow to the disturbed surface through the tangential stress balance equation (3.5). In addition, this perturbation shear stress drives an unsteady advective flow under the disturbance with non-zero perturbation streamwise and cross-streamwise velocity components because the leading-order stream function is non-zero ( $\unicode[STIX]{x1D719}_{0}\neq 0$ ). The energy of this advective flow is supplied by the work done by the perturbation shear stress. The perturbation cross-streamwise velocity component pushes the fluid towards the disturbance crest through the elastic and viscous stresses, while the hydrostatic pressure and the surface tension pull the fluid away from the disturbance crest. The net force intensifies the growth of an infinitesimal disturbance, and this yields a destabilizing effect of the viscoelastic parameter.

5 Stability analysis for an arbitrary wavenumber

Numerical solution of the OS BVP (3.8)–(3.12) is carried out for disturbances of arbitrary wavenumbers. The Chebyshev spectral collocation method (Schmid & Henningson Reference Schmid and Henningson2001) is used to recast the OS BVP into a matrix equation (Or Reference Or1997; Or & Kelly Reference Or and Kelly1998),

where $\unicode[STIX]{x1D731}=[\unicode[STIX]{x1D719}_{0},\unicode[STIX]{x1D719}_{1},\ldots ,\unicode[STIX]{x1D719}_{m},\unicode[STIX]{x1D702}]^{\text{T}}$ is a column matrix and $\mathbb{A}$ , $\mathbb{B}$ , $\mathbb{F}_{c}$ and $\mathbb{F}_{s}$ are $(m+1)\times (m+1)$ square matrices. As the Chebyshev polynomials are defined over the domain $[-1,1]$ , the liquid layer domain $y\in [-1,0]$ is shifted to $z\in [-1,1]$ by applying a transformation $y=(z-1)/2$ . Consequently, the corresponding derivatives are replaced by ${\mathcal{D}}\rightarrow 2{\mathcal{D}}$ , ${\mathcal{D}}^{2}\rightarrow 4{\mathcal{D}}^{2}$ , and so on. Next, the matrix equation (5.1) is resolved based on the Floquet theory along with the Newton–Raphson iterative scheme (Or Reference Or1997). In this method, any disturbance can be expressed in the form of a truncated complex Fourier series,

where $\unicode[STIX]{x1D6F7}_{n}$ represents the constant coefficient vectors and $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FF}_{r}+\text{i}\unicode[STIX]{x1D6FF}_{i}$ is the complex Floquet exponent. By substituting (5.2) into (5.1) and collecting the coefficients of Fourier components $\exp [(\text{i}n+\unicode[STIX]{x1D6FF})t]$ , we obtain a matrix difference equation,

where $\mathbb{F}=(\mathbb{F}_{c}+\text{i}\mathbb{F}_{s})/2$ and $\mathbb{F}^{\ast }$ is the complex conjugate of $\mathbb{F}$ . At $n=k$ , the difference equation (5.3) is reduced to the following form:

where $\mathbb{R}_{k}=-[\mathbb{A}-2\unicode[STIX]{x1D6FD}^{2}(\unicode[STIX]{x1D6FF}+\text{i}k)\mathbb{B}]^{-1}\mathbb{F}^{\ast }$ . By successive backward substitution, we can show

At $n=-k$ , the difference equation (5.3) is reduced to the following form:

where $\mathbb{R}_{-k}=-[\mathbb{A}-2\unicode[STIX]{x1D6FD}^{2}(\unicode[STIX]{x1D6FF}-\text{i}k)\mathbb{B}]^{-1}\mathbb{F}$ . Similarly, by successive forward substitution, we can show

Using (5.5) and (5.7), at $n=0$ , (5.3) can be recast into a homogeneous algebraic system,

Hence, for a non-trivial solution of $\unicode[STIX]{x1D6F7}_{0}$ , the determinant must vanish, which implies

Here, we shall focus solely on the synchronous solutions $(\unicode[STIX]{x1D6FF}_{i}=0)$ . In particular, the matrix equation (5.1) belongs to a class of equations that cannot describe the subharmonic solutions (Or Reference Or1997).

Figure 3. (a) The neutral curve in the ( $\unicode[STIX]{x1D6FD}$ , $Re$ ) plane when $\unicode[STIX]{x1D712}=1$ and $\unicode[STIX]{x1D701}=0.05$ . Here, the thick and thin lines stand for $M=0.005$ and $M=0$ respectively. The branch points are indicated by solid circles ( $M=0.005$ ) and stars ( $M=0$ ). (b) The variation of the critical wavenumber as a function of $\unicode[STIX]{x1D6FD}$ . Here, the solid and dashed lines stand for $M=0.005$ and $M=0$ respectively. The points stand for the result of Or (Reference Or1997).

All of the numerical computations are performed for 21 Chebyshev modes ( $m=21$ ) and 28 Fourier modes ( $k=14$ ), which is sufficient to achieve accurate numerical results. Further, if the number of Chebyshev modes is increased by 2, the results change by less than $0.1\,\%$ only in the higher-frequency regime ( $\unicode[STIX]{x1D6FD}>6$ ). The validation of the numerical result and the long-wavelength analytical result can be found in table 1. In order to explore the finite wavenumber results and to compare with those for a Newtonian liquid (Or Reference Or1997), the numerical test is carried out for two different values of the viscoelastic parameter $M$ when $\unicode[STIX]{x1D712}=1$ and $\unicode[STIX]{x1D701}=0.05$ are fixed. Here, we choose $M=0$ and $M=0.005$ . The results are illustrated in figure 3. Unlike the long-wave analytical result, the U-shaped neutral curves appear in different frequency ranges at $k\rightarrow 0$ . If the wavenumber is increased gradually, branch points are detected on the U-shaped neutral curves. This fact indicates that one branch may continue with the Reynolds number at long wavelength ( $k\rightarrow 0$ ), while another branch may continue with the frequency at finite wavenumber. Indeed, the new finite- $k$ neutral curves, indicated by an oblique thick line ( $M=0.005$ ) and an oblique thin line ( $M=0$ ), emerge from the branch points and grow monotonically with $\unicode[STIX]{x1D6FD}$ until the curves intersect the next U-shaped neutral curve. In fact, these new finite- $k$ oblique neutral curves are no longer confined to the regime of long-wave unstable frequency bandwidths. Instead, these curves lie in the regime of long-wave stable frequency bandwidths. Consequently, the long-wave stabilizing effect of viscoelasticity in these stable frequency bandwidths is no longer valid when the wavenumber is finite. In particular, in these stable frequency bandwidths, the instability is dominated by the finite wavenumber modes. In addition, in the long-wave frequency bandwidths, the most unstable mode is intensified in the presence of the viscoelastic parameter, which fully agrees with the analytical result discussed in § 4. Obviously, the viscoelastic parameter has a significant impact on the instability in the higher-frequency zone $\unicode[STIX]{x1D6FD}>3$ ; i.e. the presence of viscoelasticity makes the intensities of the long-wave and finite wavenumber modes stronger. The separation distance between branch points corresponding to Newtonian liquid and non-Newtonian liquid grows successively with increasing value of $M$ . Moreover, the finite- $k$ oblique neutral curves emerging from the branch points have finite critical wavenumbers $k_{c}$ , which alter with the dimensionless oscillation frequency $\unicode[STIX]{x1D6FD}$ (see figure 3 b). It should be noted that the viscoelastic parameter has a significant influence also on the finite critical wavenumbers $k_{c}$ . It is observed that the critical wavenumbers corresponding to all of the branches of the finite- $k$ neutral curves grow monotonically with $\unicode[STIX]{x1D6FD}$ , and this fact is in contrast to the result for Newtonian liquid (Or Reference Or1997), where the critical wavenumbers $k_{c}$ associated with the higher branches of the finite- $k$ neutral boundaries approach a constant value asymptotically at sufficiently large frequency $\unicode[STIX]{x1D6FD}$ .

Now, the neutral curves are analysed thoroughly in the local neighbourhood of the branch points. The first branch point is detected at $\unicode[STIX]{x1D6FD}=2.4925$ , which is different from the value of $\unicode[STIX]{x1D6FD}=2.563$ obtained for a Newtonian liquid (Or Reference Or1997). Figure 4(a) shows the neutral curves in the $(k,Re)$ plane for two different values of $\unicode[STIX]{x1D6FD}$ . The lower solid neutral curve shows a monotonic shape at a slightly lower value, $\unicode[STIX]{x1D6FD}=2.4922$ , than the branch point. Therefore, this neutral curve has a single minimum at $k\rightarrow 0$ , which signals that the long-wave mode is more unstable than the finite wavenumber mode. On the contrary, the upper solid neutral curve has a non-monotonic shape at a slightly higher value, $\unicode[STIX]{x1D6FD}=2.4935$ , than the branch point. The local minimum is detected at $k=0.22$ , which signals that the finite wavenumber mode is more dangerous than the long-wave mode. Obviously, the long-wave instability switches to the finite wavenumber instability exactly at the branch point. In fact, this switching process proceeds in a continuous fashion, which is illustrated by a dash-dotted line in figure 4(a). The variation of the neutral curves in the local neighbourhood of the second branch point, $\unicode[STIX]{x1D6FD}=5.1590$ , is shown in figure 4(b). It should be noted that both curves show a non-monotonic shape, which is in contrast to the previous result. When $\unicode[STIX]{x1D6FD}=5.1586$ , the neutral curve has a local minimum at $k\rightarrow 0$ , which indicates the dominant effect of the long-wave mode. On the other hand, the local minimum shifts to $k=0.46$ at $\unicode[STIX]{x1D6FD}=5.1610$ , which indicates the dominant effect of the finite wavenumber mode. Again, the long-wave instability switches to the finite wavenumber instability through the second branch point, but in a discontinuous fashion. It seems that a competition occurs between the long-wave mode and the finite wavenumber mode in the vicinity of the branch points. Figure 4(c) illustrates the resulting neutral curves in the local neighbourhood of the third branch point at $\unicode[STIX]{x1D6FD}=7.72$ . In this case, the lower solid neutral curve shows a monotonic shape at a slightly lower value, $\unicode[STIX]{x1D6FD}=7.7165$ , than the branch point. Therefore, this neutral curve has a single minimum at $k\rightarrow 0$ , which signals that the long-wave mode is more unstable than the finite wavenumber mode. On the contrary, the upper solid neutral curve has a non-monotonic shape at a slightly higher value, $\unicode[STIX]{x1D6FD}=7.725$ , than the branch point. The local minimum appears at $k=0.073$ , which is a small value in comparison with the previous two cases. However, the dash-dotted curve representing the stability limit grows monotonically with wavenumber $k$ . Therefore, in this case also, the long-wave instability switches to the finite wavenumber instability through the third branch point. However, this switching process proceeds in a continuous fashion as in the first case.

Figure 4. (a) The neutral curves in the local neighbourhood of the first branch point when $\unicode[STIX]{x1D712}=1$ , $M=0.005$ and $\unicode[STIX]{x1D701}=0.05$ . (b) The neutral curves in the local neighbourhood of the second branch point. (c) The neutral curves in the local neighbourhood of the third branch point.

Now, the numerical test is performed for a comparatively large value of the viscoelastic parameter, $M=0.01$ , as proposed by Dandapat & Gupta (Reference Dandapat and Gupta1975). The associated results are depicted in figure 5(a). Here, we shall focus only on the third U-shaped neutral curve because the finite wavenumber neutral curve emerging from the branch point of the third U-shaped neutral curve is no longer an oblique line. Therefore, the finite wavenumber neutral curve requires thorough investigation close to the third branch point in the $(k,Re)$ plane. To do this, the numerical test is performed with $\unicode[STIX]{x1D6FD}=7.94$ , a value on the long-wave neutral curve. Figure 5(b) shows the associated neutral curve whose local minimum appears at $k\rightarrow 0$ , and this fact indicates the dominant effect of the long-wave mode on the instability. If $\unicode[STIX]{x1D6FD}$ continues to increase gradually along the finite wavenumber neutral curve through the branch point, the local minimum shifts to the finite wavenumber region, and this fact indicates the dominant effect of the finite wavenumber mode on the instability. Again, the long-wave instability switches to the finite wavenumber instability. The associated neutral curves exhibit a tongue-like shape. In this case, there exists a stable bandwidth because the finite wavenumber neutral curve emerging from the branch point of the third U-shaped long-wave neutral curve never intersects the next U-shaped long-wave neutral curve.

Figure 5. (a) The neutral curves in the ( $\unicode[STIX]{x1D6FD}$ , $Re$ ) plane when $M=0.01$ , $\unicode[STIX]{x1D712}=1$ and $\unicode[STIX]{x1D701}=0.05$ . Here, the thick and thin lines represent the results corresponding to the long-wavelength and finite wavelength analyses. The branch points are indicated by solid circles. (b) The scenario of neutral curves in the $(k,Re)$ plane for different values of $\unicode[STIX]{x1D6FD}$ .