1. Introduction {#sec1}
===============

Over the past few decades, the focus on viscoelastic fluids behavior especially in oscillatory conditions has been increased substantially due to their wide use in many industries. The significant complexities in mathematical modeling of viscoelastic fluid flow at practical scenarios is the main challenge in determination of flow characteristics.

The response of viscoelastic fluid to the oscillating force was first examined using basic models such as Maxwell model \[[@bib1], [@bib2]\]. In addition to the basic Maxwell model, several models based on the rheology theories are developed and implemented up to now (Bingham model \[[@bib3]\] and Phan-Thien and Tanner (PTT) model \[[@bib4]\]).

Moreover, Oldroyd- B model was used in \[[@bib5]\] to study the flow of polymer dilute solvents on a flexible surface in the range of low Reynolds number. Fluid parameters including solution to solvent viscosity ratio and relaxation time was investigated using the Weissenberg number. According to their results, the flow was stable prior to reaching the maximum *We*.

One of the most familiar examples of oscillating non-Newtonian flood is blood. In the category of physiological applications, the viscoelastic nature of blood \[[@bib6], [@bib7]\] was investigated for low and high oscillation frequencies. The unsteady flow of blood through stenosed artery \[[@bib8]\] was also simulated, based on an oscillatory pressure gradient as the driving force. Blood has also been treated as a non-Newtonian fluid characterized by the Oldroyd-B and Cross models \[[@bib9]\]. According to the evaluated results, the magnitude of normal stress is reported to be much higher than shear stress and the increasing of relaxation time results in higher values of normal stress, as well.

The non-Newtonian viscoelastic fluid between parallel oscillating plates occurs in several applications namely, lubrication of mechanical devices, biomedical engineering and oil industry \[[@bib10]\]. Fluid near the plates would oscillate harmonically with the plate in similar time period, while the amplitude of fluid motion rises slowly across the direction of the guiding surface.

Incompressible fluid flow with porosity effects are examined due to sine and cosine oscillations of infinite plates and the exact solutions were obtained for velocity field corresponding to shear stress using integral transforms techniques (Laplace and Fourier sine transformations) \[[@bib11]\]. Based on the results, the velocity field and the shear stress are decreasing functions of oscillations amplitude. In addition, the dependence of shear stress on the affecting parameters such as relaxation and retardation time, for hydro-magnetic, viscoelastic fluid passing over an oscillating surface is investigated using Oldroyd model for free convection flow \[[@bib12]\].

In 2017, the effect of oscillatory vibration of a thin blade on the non-Newtonian fluid flow inside a channel is numerically examined \[[@bib13]\]. The power-law model for shear thinning and shear thicking fluid is implemented and the heat transfer enhancement due to the presence of vibrating blade is investigated.

Despite the fact that up to now, several investigations focus on the non-Newtonian fluid flow between parallel plates, the lack of a comprehensive examination regarding the flow stability analysis and parameter study is still existed. Some of the recent researches focusing on the Newtonian and non-Newtonian flow between parallel plates will be critically reviewed hereafter. Feng et al. in 2017 \[[@bib14]\], numerically simulate the generalized Oldroyd-B fluid flow between parallel plates. The main novelty of the mentioned study is the utilization of time and space fractional calculus to expand the generality of the proposed method. However, the simulated problem is limited to moving boundary condition and the effect of fluid oscillations at the inlet is not accounted. Moreover, the investigation is merely focused on the numerical scheme proposition and the physical parameters study is the significant missing data. Keimanesh et al. \[[@bib15]\] also conduct similar investigation on the non-Newtonian fluid between two parallel plates using the multi-step differential transform method for non-oscillatory flow condition.

The effect of parallel plates' movement normal to the fluid flow is also examined by several researchers. Rashidi et al. in 2008 \[[@bib16]\] tackled the problem of squeezing flow between parallel plates using the Homotopy analysis method. In another investigation, Hoshyar et al. \[[@bib17]\] focused on unsteady incompressible Newtonian fluid flow between two parallel plates using the same method. Despite the fact these studies are dealing with the flow between parallel plates subjected to external movement effects, the non-Newtonian as well as oscillatory fluid motion is still remained as an intact expansion of the pointed out studies.

The oscillating Couette--Poiseuille flow for the Oldroyd B fluid is examined by Ma et al. in 2019 \[[@bib18]\]. The proposed analytical method is utilized to investigate the inertial, viscous, and elastic effects on the flow behavior and stress responses. However, the range of flow stability and the bifurcation curves, which can be implemented to significantly increase the knowledge of flow behavior, is not examined in the mentioned research. Therefore, a comprehensive study of oscillatory and non-Newtonian fluid flow using Johnson--Segalman model and the flow stability analysis are still intact problems, which should be examined for specific applications.

The purpose of the present study is to investigate the steady state and transient behavior of the viscoelastic fluid with time varying inlet condition. The effect of non-linear system is considered through evaluation of time dependent governing equation with Galerkin projection method. The effect of inlet boundary condition fluctuations is investigated via introduction of oscillatory Weissenberg function. Regarding the critical Weissenberg numbers captured form the steady state solution, first the system unsteady behavior in constant Weissenberg number is obtained, then two higher and lower Weissenberg number values with respect to the critical magnitudes are chosen as the amplitude of the sinusoidal Weissenberg number for following examined cases. Moreover, the velocity and stresses plots are analyzed in various intervals of Weissenberg number. The effect of Reynolds number and polymer solute to solution viscosity ratio on velocity, normal and shear stresses are the other examined parameters.

2. Theory {#sec2}
=========

2.1. The fluid flow laws {#sec2.1}
------------------------

The present study investigates the one-dimensional flow of a viscoelastic fluid between two infinite parallel plates. The plates are stationary and separated by a distance of *d*. The isothermal fluid with single relaxation time, constant viscosity and time dependent inlet velocity is examined. Moreover, the viscosity of the polymer solution ($\text{η}$) is determined based on the superposition of Newtonian solvent viscosity ($\text{η}$~1~) and the non-Newtonian polymer with viscosity of $\text{η}$~2~ as the solute fluid. The flow is considered as fully developed with fluctuations occur merely in one direction and no-slip condition at the walls. The schematic of the described model is illustrated in [Fig. 1](#fig1){ref-type="fig"}.Fig. 1Schematic of the oscillatory flow passage through the channel.Fig. 1

Generally, the flow behavior is governed by the conservation laws of mass and linear momentum for an incompressible fluid.

Conservation of mass:$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0$$

Conservation of momentum in x-direction$$\rho\left( \frac{\partial u}{\partial t} \right) = \frac{- \partial P}{\partial x} + \frac{\partial\tau_{xy}}{\partial y} + \eta_{2}\frac{\partial^{2}u}{\partial y^{2}}$$where *u* and *v* are velocity components in tangential and normal directions relative to the base flow, *P* is the pressure, *t* denotes the time, *τ* is the shear stress and *η*~*2*~ is the viscosity of the solvent.

The constitutive equation is selected according to the Johnson-Segalman model as:$$\lambda\frac{DT}{Dt} + T = \eta_{1}\left\lbrack {\nabla U + \left( \nabla U \right)}^{t} \right\rbrack$$where *λ* is a dimensionless parameter (relaxation time), *T* is stress tensor, *∇U* denotes the gradient of velocity vector and $\nabla U^{t}$ is the transpose of the *∇U*. Additionally, $\frac{DT}{Dt}$ can be calculated as:$$\frac{DT}{Dt} = \frac{\partial T}{\partial t} + U.\nabla T - \left( {1 - \frac{\zeta}{2}} \right)\left( {\left( {\nabla U} \right)^{t}.T + T.\nabla U} \right) + \frac{\zeta}{2}\left( {\nabla U.T + T.\left( {\nabla U} \right)^{t}} \right)$$where $\text{ζ} \in \left\lbrack 0 - 2 \right\rbrack$ is a non-dimensional fluid parameter which is a measure of non-uniform motion distribution in the stress tensor.

For the present specific case the Johnson-Segalman equation after some simplifications is converted to the following equations.$$T_{xx} = - \lambda\frac{\partial T_{xx}}{\partial t} - \lambda\left( \zeta - 2 \right)\frac{\partial u}{\partial y}T_{xy}$$$$T_{xy} = \eta_{1}\frac{\partial u}{\partial y} - \lambda\left\{ {\frac{\partial T_{xy}}{\partial t} + \frac{\zeta}{2}T_{xx}\frac{\partial u}{\partial y} + \left( {\frac{\zeta}{2} - 1} \right)\frac{\partial u}{\partial y}T_{yy}} \right\}$$$$T_{yy} = - \ \lambda\left\{ {\frac{\partial T_{yy}}{\partial t} + \zeta\frac{\partial u}{\partial y}T_{xy}} \right\}$$

2.2. Non-dimensional form of the governing equations {#sec2.2}
----------------------------------------------------

Taking *U* to be the maximum velocity in the flow direction, the non-dimensional variables can be defined as:$$u = u^{\ast}\frac{d}{\lambda},\text{~\!}p = \ p^{\ast}\frac{\eta_{1}U}{d},\text{~\!}T = \frac{T^{\ast}\eta_{1}}{\lambda},\text{~\!}t = t^{\ast}\lambda,\text{~\!}y = y^{\ast}d,\text{~\!}x = x^{\ast}d$$

By substitution of these transformations into Eqs. [(5)](#fd5){ref-type="disp-formula"}, [(6)](#fd6){ref-type="disp-formula"}, and [(7)](#fd7){ref-type="disp-formula"} the following equations are formed:$$\overline{Re}\frac{\partial u}{\partial t} = \  - We\left( t \right)\ \frac{dP}{dx} + \ \frac{\partial\tau_{xy}}{\partial y} + \ \varepsilon\ \frac{\partial^{2}u}{\partial y^{2}}$$$$\frac{\partial\tau_{xx}}{\partial t} = \  - \tau_{xx} + \ \left( {2 - \zeta} \right)\frac{\partial u}{\partial y}\tau_{xy}$$$$\frac{\partial\tau_{yy}}{\partial t} = - \tau_{yy} - \zeta\frac{\partial u}{\partial y}\tau_{xy}$$$$\frac{\partial\tau_{xy}}{\partial t} = \  - \tau_{xy} + \frac{\partial u}{\partial y} - \frac{\zeta}{2}\frac{\partial u}{\partial y}\tau_{xx} + \left( {1 - \frac{\zeta}{2}} \right)\frac{\partial u}{\partial y}\tau_{yy}$$

The boundary condition is defined as:$$u\left( {y = \mp \frac{1}{2},t} \right) = 0$$

Additionally, the three important similarity groups in the present problem, namely, the Reynolds number, *Re*, the Weissenberg number, *We*, the solvent to solute viscosity ratio, $\text{ε}$ ,and the oscillatory inlet velocity are defined as:$$U = U_{{^\circ}\ }sin\left( \omega t \right)$$$$We\left( t \right) = \frac{U_{{^\circ}\ }\lambda\ }{d}sin\left( \omega t \right)$$$$\overline{Re} = \frac{d^{2}\rho}{\eta_{1}\lambda}$$$$\varepsilon = \frac{\eta_{2}}{\eta_{1}}$$

2.3. The base flow governing equations {#sec2.3}
--------------------------------------

The steady state condition is determined via neglecting of time derivative terms in Eqs. [(9)](#fd9){ref-type="disp-formula"}, [(10)](#fd10){ref-type="disp-formula"}, [(11)](#fd11){ref-type="disp-formula"}, and [(12)](#fd12){ref-type="disp-formula"}. The steady state stress components are obtained as follows:$$\tau_{xy}^{b} = \frac{u_{'y}}{1 + 2\zeta u_{'y}^{2} - \zeta^{2}u_{'y}^{2}} = \frac{u_{'y}}{1 + \zeta\left( 2 - \zeta \right)u_{'y}^{2}} = \frac{A^{2}y^{2}}{1 + \zeta\left( 2 - \zeta \right)A^{2}y^{2}}$$$$\tau_{xx}^{b} = \frac{\left( 2 - \zeta \right)u_{'y}^{2}}{1 + \zeta\left( 2 - \zeta \right)u_{'y}^{2}} = \frac{\left( 2 - \zeta \right)A^{2}y^{2}}{1 + \zeta\left( 2 - \zeta \right)A^{2}y^{2}}$$$$\tau_{yy}^{b} = \frac{- \zeta u_{'y}^{2}}{1 + \zeta\left( 2 - \zeta \right)u_{'y}^{2}} = \frac{- \zeta A^{2}y^{2}}{1 + \zeta\left( 2 - \zeta \right)A^{2}y^{2}}$$where the *A* function is defined as:$$A = \frac{We\left( t \right)}{\varepsilon}\frac{dP}{dx}$$

Eqs. [(5)](#fd5){ref-type="disp-formula"}, [(6)](#fd6){ref-type="disp-formula"}, and [(7)](#fd7){ref-type="disp-formula"} can be rewritten by the definition of normal stresses difference (*N*) and a combination of normal stresses (*Z*) as:$$\begin{matrix}
{N = \tau_{xx} - \tau_{yy}} & , & {Z = \frac{\zeta}{2}\tau_{xx} + \left( 1 - \frac{\zeta}{2} \right)\tau_{yy}} \\
\end{matrix}$$

Moreover, to investigate the created disturbances in the flow, the following modifications are applied:$$\begin{matrix}
{u = u^{b} + u^{'},} & {\tau_{ij} = \tau_{ij}^{b} + \tau_{ij}^{'},} & {p = p^{b} + p^{'},} \\
\end{matrix}$$

For simplicity, the first term on the right-hand side of [Eq. (9)](#fd9){ref-type="disp-formula"} is redefined as:$$B\left( t \right) = - We\left( t \right)\frac{dP}{dx}$$

Considering the relation between *A* and *B* parameters and the implementation of [Eq. (23)](#fd23){ref-type="disp-formula"} modifications into Eqs. [(9)](#fd9){ref-type="disp-formula"}, [(10)](#fd10){ref-type="disp-formula"}, [(11)](#fd11){ref-type="disp-formula"}, and [(12)](#fd12){ref-type="disp-formula"}, the following governing equations are obtained:$$\overline{Re}u_{'t} = S_{'y} + S_{'y}^{b} + \varepsilon u_{'yy}$$$$N_{'t} = - N + 2\left( {A\left( t \right)yS + u_{'y}S^{b} + u_{'y}S} \right)$$$$S_{'t} = - S + u_{'y} + \alpha\left( {A\left( t \right)yN + u_{'y}N^{b} + u_{'y}N} \right)$$

All of the unknown parameters in Eqs. [(25)](#fd25){ref-type="disp-formula"}, [(26)](#fd26){ref-type="disp-formula"}, and [(27)](#fd27){ref-type="disp-formula"} are defined as:$$N^{b} = \frac{2{u_{'y}^{b}}^{2}}{1 + \zeta\left( 2 - \zeta \right){u_{'y}^{b}}^{2}}$$$$S^{b} = \frac{u_{'y}^{b}}{1 + \zeta\left( 2 - \zeta \right){u_{'y}^{b}}^{2}}$$$$u_{'y}^{b} = Ay$$$$A = \frac{We}{\varepsilon}\frac{dP}{dx}$$$$\alpha = \ \zeta\left( {\frac{\zeta}{2} - 1} \right)$$

2.4. Solution method {#sec2.4}
--------------------

The solution of Eqs. [(25)](#fd25){ref-type="disp-formula"}, [(26)](#fd26){ref-type="disp-formula"}, and [(27)](#fd27){ref-type="disp-formula"} is determined using the Galerkin projection method. To do so, *U(y,t)*, *S(y,t)* and *N(y,t)* are introduced by Chandrasekhar Function which satisfies the no-slip boundary condition.

Therefore, the general series representation of the velocity and stress differences is given by Eqs. [(33)](#fd33){ref-type="disp-formula"}, [(34)](#fd34){ref-type="disp-formula"}, and [(35)](#fd35){ref-type="disp-formula"} as:$$u\left( {y,t} \right) = \sum\limits_{i = 1}^{M}{U_{i}\left( t \right)\phi_{i}\left( y \right)}$$$$N\left( {y,t} \right) = \sum\limits_{i = 1}^{M}{N_{i}\left( t \right)\phi_{i}\left( y \right)}$$$$S\left( {y,t} \right) = \sum\limits_{i = 1}^{M}{S_{i}\left( t \right)\phi_{i}^{'}\left( y \right)}$$where $\phi_{i}\left( y \right)$ is the even and odd Chandrasekhar functions, for $i$ even and odd, respectively \[[@bib12]\], $\phi_{i}^{'}\left( y \right)$ are defined as $\phi_{i}^{'}\left( y \right) = \left( {1/\alpha_{i}} \right)\left( {d\phi_{i}/dy} \right)$, where $\alpha_{i}$ are constants and *M* is the number of modes. The suitable level of truncation that is imposed hereafter in the calculations is *M* = 2 which leads to a six-dimensional model.

The first step in Galerkin projection method is the substitution of Eqs. [(33)](#fd33){ref-type="disp-formula"}, [(34)](#fd34){ref-type="disp-formula"}, and [(35)](#fd35){ref-type="disp-formula"} into Eqs. [(25)](#fd25){ref-type="disp-formula"}, [(26)](#fd26){ref-type="disp-formula"}, and [(27)](#fd27){ref-type="disp-formula"}. Then each equation is multiplied by the suitable mode and the integration is performed in the range related to the problem geometry (et. $\left. y \in \left\lbrack \frac{- 1}{2} - \frac{1}{2} \right\rbrack \right)$. Hence, a set of ordinary, coupled and non-linear differential equations are derived with time dependent coefficients.

The projection method leads to explicit expressions for time derivatives (et. $\text{u}_{\text{k}}$ and $\text{N}_{\text{k}}$ where $\text{k} \in \left\lbrack 1 - \text{M} \right\rbrack$):$$\overline{Re}\frac{\partial U_{1}}{\partial t} = S_{1}\alpha_{1}D_{1} + S_{2}\alpha_{2}D_{2} + \varepsilon U_{1}\alpha_{1}^{2}D_{1} + \varepsilon U_{2}\alpha_{2}^{2}D_{2} + E_{1}$$$$\overline{Re}\frac{\partial U_{2}}{\partial t} = S_{1}\alpha_{1}D_{3} + S_{2}\alpha_{2}D_{4} + \varepsilon U_{1}\alpha_{1}^{2}D_{3} + \varepsilon U_{2}\alpha_{2}^{2}D_{4} + E_{2}$$$$\frac{\partial N_{1}}{\partial t} = - N_{1} + 2A\left( t \right)S_{1}D_{5} + 2A\left( t \right)S_{2}D_{6} + 2\alpha_{1}U_{1}D_{10} + 2\alpha_{2}U_{2}D_{11} + 2\alpha_{1}U_{1}S_{1}D_{7} + 2\alpha_{1}U_{1}S_{2}D_{8} + 2\alpha_{2}U_{2}S_{1}D_{8} + 2\alpha_{2}U_{2}S_{2}D_{9}$$$$\frac{\partial N_{2}}{\partial t} = - N_{2} + 2A\left( t \right)S_{1}D_{12} + 2A\left( t \right)S_{2}D_{13} + 2\alpha_{1}U_{1}D_{14} + 2\alpha_{2}U_{2}D_{15} + 2\alpha_{1}U_{1}S_{1}D_{16} + 2\alpha_{1}U_{1}S_{2}D_{17} + 2\alpha_{2}U_{2}S_{1}D_{17} + 2\alpha_{2}U_{2}S_{2}D_{18}$$$$\frac{\partial S_{1}}{\partial t} = - S_{1} + \alpha_{1}U_{1}\left( {1 + \alpha E_{3}} \right) + \alpha\alpha_{2}U_{2}E_{4} + \alpha A\left( t \right)N_{2}D_{19} + \alpha\alpha_{1}U_{1}N_{1}D_{20} + \alpha\alpha_{2}U_{2}N_{2}D_{21}$$$$\frac{\partial S_{2}}{\partial t} = - S_{2} + \alpha_{2}U_{2}\left( {1 + \alpha E_{5}} \right) + \alpha A\left( t \right)N_{1}D_{22} + \alpha\alpha_{1}U_{1}E_{6} + \alpha\alpha_{1}U_{1}N_{2}D_{23} + \alpha\alpha_{2}U_{2}N_{1}D_{24}$$

The coefficients in Eqs. [(36)](#fd36){ref-type="disp-formula"}, [(37)](#fd37){ref-type="disp-formula"}, [(38)](#fd38){ref-type="disp-formula"}, [(39)](#fd39){ref-type="disp-formula"}, [(40)](#fd40){ref-type="disp-formula"}, and [(41)](#fd41){ref-type="disp-formula"} are defined as:$$\begin{matrix}
{\text{D}_{1} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{pp}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{2} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{2}^{\text{pp}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{3} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{pp}}\varnothing_{2}\text{dy}}}} \\
{{\mspace{9mu}\text{D}}_{4} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{2}^{\text{pp}}\varnothing_{2}\text{dy}}}} & {{\mspace{9mu}\text{E}}_{1} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\text{S}_{\text{'y}}^{\text{b}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{E}}_{2} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\text{S}_{\text{'y}}^{\text{b}}\varnothing_{2}\text{dy}}}} \\
{{\mspace{9mu}\text{D}}_{5} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{6} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\text{y}\varnothing_{2}^{\text{p}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{7} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{1}^{\text{p}}\varnothing_{1}\text{dy}}}} \\
{{\mspace{9mu}\text{D}}_{8} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{2}^{\text{p}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{9} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{2}^{\text{p}}\varnothing_{2}^{\text{p}}\varnothing_{1}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{10} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{1}\text{S}^{\text{b}}\text{dy}}}} \\
{{\mspace{9mu}\text{D}}_{11} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{2}^{\text{p}}\varnothing_{1}\text{S}^{\text{b}}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{12} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\text{y}\varnothing_{1}^{\text{p}}\varnothing_{2}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{13} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\text{y}\varnothing_{2}^{\text{p}}\varnothing_{2}\text{dy}}}} \\
{{\mspace{9mu}\text{D}}_{14} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{2}\text{S}^{\text{b}}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{16} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{1}^{\text{p}}\varnothing_{2}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{15} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{2}^{\text{p}}\varnothing_{2}\text{S}^{\text{b}}\text{dy}}}} \\
{{\mspace{9mu}\text{D}}_{17} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{1}^{\text{p}}\varnothing_{2}^{\text{p}}\varnothing_{2}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{18} = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\varnothing_{2}^{\text{p}}\varnothing_{2}^{\text{p}}\varnothing_{2}\text{dy}}}} & {{\mspace{9mu}\text{D}}_{19} = \text{y}\frac{\left( \frac{\varnothing_{2}}{\varnothing_{1}^{\text{p}}} \right)}{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}} \\
{{\mspace{9mu}\text{D}}_{20} = \frac{\left( \frac{\varnothing_{1}^{\text{p}}\varnothing_{1}}{\varnothing_{1}^{\text{p}}} \right)}{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}} & {{\mspace{9mu}\text{D}}_{21} = \frac{\left( \frac{\varnothing_{2}^{\text{p}}\varnothing_{2}}{\varnothing_{1}^{\text{p}}} \right)}{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}} & {{\mspace{9mu}\text{D}}_{22} = \text{y}\frac{\left( \frac{\varnothing_{1}}{\varnothing_{2}^{\text{p}}} \right)}{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}} \\
{{\mspace{9mu}\text{D}}_{23} = \frac{\left( \frac{\varnothing_{1}^{\text{p}}\varnothing_{2}}{\varnothing_{2}^{\text{p}}} \right)}{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}} & {{\mspace{9mu}\text{D}}_{24} = \frac{\left( \frac{\varnothing_{2}^{\text{p}}\varnothing_{1}}{\varnothing_{2}^{\text{p}}} \right)}{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}} & {{\mspace{9mu}\text{E}}_{3} = \text{N}^{\text{b}}\frac{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}} \\
{{\mspace{9mu}\text{E}}_{4} = \text{N}^{\text{b}}\frac{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{1}^{\text{p}}} \right)}} & {{\mspace{9mu}\text{E}}_{5} = \text{N}^{\text{b}}\frac{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}} & {{\mspace{9mu}\text{E}}_{6} = \text{N}^{\text{b}}\frac{\left( \frac{\varnothing_{1}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}{\left( \frac{\varnothing_{2}^{\text{p}}}{\varnothing_{2}^{\text{p}}} \right)}} \\
\end{matrix}$$where $\varnothing_{\mathbf{i}}^{\mathbf{p}\mathbf{p}}$ denotes second derivative of $\varnothing_{\mathbf{i}}$. The () operator in Eq. [(42)](#fd42){ref-type="disp-formula"} is defined as:$$\left( \frac{\text{a}}{\text{b}} \right) = {\int\limits_{\frac{- 1}{2}}^{\frac{1}{2}}{\text{ab~\!}\mathbb{d}\text{y}}}$$

Briefly, the nonlinear governing equations are transformed into a set of simultaneous ordinary differential equations (ODE) using the Galerkin method (Eqs. [(36)](#fd36){ref-type="disp-formula"}, [(37)](#fd37){ref-type="disp-formula"}, [(38)](#fd38){ref-type="disp-formula"}, [(39)](#fd39){ref-type="disp-formula"}, [(40)](#fd40){ref-type="disp-formula"}, and [(41)](#fd41){ref-type="disp-formula"}). The coefficients of the ODE set can be obtained using the integration of the Chandrasekhar function modes. The numerical simulation of the Mathematica software is performed with the aid of a computer system with following characteristics: Intel Core i7 <CPU@2.8GHz>, 16 GB of DDR4 RAM. Each simulation is converged in about 3 min to the final results.

3. Results & discussion {#sec3}
=======================

3.1. The bifurcation graph {#sec3.1}
--------------------------

In the present section the numerical solution of the governing equations (Eqs. [(36)](#fd36){ref-type="disp-formula"}, [(37)](#fd37){ref-type="disp-formula"}, [(38)](#fd38){ref-type="disp-formula"}, [(39)](#fd39){ref-type="disp-formula"}, [(40)](#fd40){ref-type="disp-formula"}, and [(41)](#fd41){ref-type="disp-formula"}) is presented for various *ξ* and *We* numbers. The variations of velocity, normal and shear stresses (namely *U*~*1*~, *U*~*2*~, *N*~*1*~, *N*~*2*~, *S*~*1*~ and *S*~*2*~) are calculated for $We \in \left\lbrack 0 - 10 \right\rbrack$ and $\delta \in \left\lbrack {0.2 - 1.0} \right\rbrack$ for limiting *Re* values of 0.2 and 10. It should be added that $\text{Δp}$ and *ε* parameters are set to -1 and 0.04, respectively.

According to [Fig. 2](#fig2){ref-type="fig"} (the bifurcation graph), the *U*~*1*~ variable has two sequential maximum and minimum values for relatively low *We* numbers which is due to the disturbance terms in [Eq. (23)](#fd23){ref-type="disp-formula"} (*i.e.* $u\text{'}$, $p\text{'}$ and $\text{τ}_{\text{ij}}^{\text{'}}$). As the value of these peaks is reduced the fluid resistance to the imposed forces and initial disturbances is decreased. Based on [Fig. 2](#fig2){ref-type="fig"}, three main zones can be defined as subcritical ($We < We_{C1}$), critical ($We_{C1} < We < We_{C2}$) and supercritical ($\left. We > We_{C2} \right)$.Fig. 2The variations of a) *U*~*1*~, and b) *U*~*2*~ as a function of *We* for $\text{ε} = 0.04$,$\text{~\!Re} = 0.2\text{~\!}$ and $\Delta\text{p} = - 1$.Fig. 2

First and second critical *We* numbers are located in the range of $\left\lbrack 0.2 - 0.6 \right\rbrack$ and $\left\lbrack 0.4 - 0.8 \right\rbrack$, respectively. Unstable base flow with nonlinear velocity profile is observed between these two critical *We* numbers.

3.2. Validation {#sec3.2}
---------------

It should be noted that, the oscillatory problem of current study using the Johnson--Segalman model is not tackled in the available literature yet. Thus, the validation is performed by the analytical data of oscillating Couette--Poiseuille flow of upper-convected Maxwell fluid of Ma et al. \[[@bib18]\], which has approximately similar problem definition in comparison to the present study. The velocity profile for Re = 10 at four different time instances during one time period of oscillations are calculated and illustrated in [Fig. 3](#fig3){ref-type="fig"}.Fig. 3Velocity profile for oscillating Couette--Poiseuille flow of upper-convected Maxwell fluid at various time instances.Fig. 3

It should be noted that the ξ and π parameters are dimensionless depth and time variables, respectively. According to the results, it can be concluded that the velocity profile of the current study is in close agreement with the analytical data of \[[@bib18]\]. The difference between the numerical results of the current study with the analytical data is below 4% at all of the time instances through the channel depth.

3.3. Transient simulation results {#sec3.3}
---------------------------------

To examine the nonlinear behavior of the system in the transient period, the time dependent outputs of the governing equation from the initial condition is studied. The effect of boundary fluctuations is considered by dedication of a sinusoidal function to the *We* ([Eq. 15](#fd15){ref-type="disp-formula"}). To do so, two *We* numbers in the stable sub/supercritical zones are selected based on the steady solution (et. 0.15 and 15) and the variations of velocity and stress components are examined for both constant and oscillatory *We* numbers. Moreover, the *Re* number is assumed to be 0.2 while the case of *Re* = 10 is also investigated to determine the effect of *Re* number on the fluid flow characteristics.

The velocity profile variations from the initial condition to the steady configuration is presented in [Fig. 4](#fig4){ref-type="fig"}. The velocity profiles approximately resemble the Newtonian velocity profile, although a slight declination to the right is observed in the profiles.Fig. 4The transient behavior of flow in the channel in subcritical zone for $\text{Re} = 0.2$, $\text{We} = 0.15$, $\Delta\text{p} = - 1$, $\text{ε} = 0.04$ and $\text{ζ} = 0.2$.Fig. 4

Considering the final configuration of the velocity profiles, it can be stated that from t = 0 to t = 3 the velocity magnitude is overgrowing, while it is decreased slightly to t = 4 when the final profile is almost shaped.

For better visualization, the y-t-u diagram is presented in [Fig. 5](#fig5){ref-type="fig"}. Considering the results, the fluid has a tendency to fluctuate immediately after the initial condition. However, the oscillations are damped to zero after limited waves near t = 4. In addition, it should be noted that the velocity is set to zero at the walls due to the no-slip condition.Fig. 5The variations of a) velocity, b) normal and c) shear stress as a function of y and t for subcritical *We* with parameters of $\text{ε} = 0.04$, $\text{~\!Re} = 0.2$, $\nabla\text{p} = - 1$, We = 0.15 and, ζ = 0.2.Fig. 5

The normal stress is equal to zero at channel walls and has an increasing trend marching toward the channel midsection, where the maximum is occurred at y = 0.12. Moreover, the early fluctuations as well as the existence of steady state after t = 4.5 is noticeable.

The shear stress as it is expected, has its maximum near the wall (y = 0.5), which is intensified until t = 2.0. Subsequently, the shear stress maximum is reduced slightly toward its stable value (t \> 4.0), which is similar to the Newtonian flow.

The effect of temporal oscillations of the inlet velocity magnitude is illustrated in [Fig. 6](#fig6){ref-type="fig"}. The oscillatory behavior of the inlet boundary is simulated for low and high frequencies (*i.e.* sin (0.5t) and sin (20t)). It can be observed that the low frequency variations of the inlet velocity magnitude have the same effect as the non-fluctuating boundary condition. The same trend can be seen for the shear stress, which leads to the conclusion that the Weissenberg fluctuations especially at low frequencies has negligible effect on both velocity and its derivatives ([Fig. 6](#fig6){ref-type="fig"}).Fig. 6The variations of a) velocity, b) normal and c) shear stress as a function of y and t for subcritical *We* with parameters of $\text{ε} = 0.04$, $\text{Re} = 0.2$, $\nabla\text{p} = - 1$, *We=0.15\*Sin(0.5t)* and ζ = 0.2.Fig. 6

In contrary to the velocity and shear stress, the normal stress is completely affected by the low frequency oscillations of *We*. The normal stress fluctuations have variable time periods from 50 to 200. The effect of inlet sinusoidal *We* on the damping and exciting terms in the *N* equation (Eqs. [(38)](#fd38){ref-type="disp-formula"} and [(39)](#fd39){ref-type="disp-formula"}) is the reason for such a behavior.

The effect of increasing of the oscillation frequency is presented in [Fig. 7](#fig7){ref-type="fig"}. The amplitude of *We* is chosen to be 0.15 and its frequency is set to be 20/2π. Velocity and shear stress behaviors are similar to the previous case with inlet velocity frequency of 0.5/2π ([Fig. 6](#fig6){ref-type="fig"}(a) and (c)). However, the amplitude of normal stress fluctuations significantly reduced (to about 0.005) with the period of 100. After t = 30 the oscillations turn into stable periodic waves with a peak value (0.022) slightly lower than the previous case (0.032) with lower oscillation frequency ([Fig. 6](#fig6){ref-type="fig"}).Fig. 7The variations of a) velocity, b) normal and c) shear stress as a function of y and t for subcritical zone with parameters of $\text{ε} = 0.04$, $\text{Re} = 0.2$, $\nabla\text{p} = - 1$, *We=0.15\*Sin(20t)* and ζ = 0.2.Fig. 7

In the current section the characteristics of the viscoelastic flow for *We* variations (constant and low frequency oscillations) is examined in the supercritical range. With the operational conditions of *We* = 15, *ε* = 0.04 and *Re* = 0.2, it is depicted in [Fig. 8](#fig8){ref-type="fig"}, that the velocity has a non-oscillatory behavior and will be uniform after about t = 30, similar to the normal and shear stress diagrams. Comparison to related figures for subcritical condition ([Fig. 5](#fig5){ref-type="fig"}), one can deduce that the location of maximum velocity relocated from channel center to near walls and the required time to stabilization of the flow is substantially increased.Fig. 8The variations of a) velocity, b) normal and c) shear stress as a function of y and t for subcritical zone with parameters of $\text{ε} = 0.04$, $\text{Re} = 0.2$, $\nabla\text{p} = - 1$, *We* = 15 and ζ = 0.2.Fig. 8

The normal stress increases till t = 30 after which a nearly uniform trend is observed with its highest value occurred at y = 0.2. It is worth noting that in contrary to the subcritical case, the normal stress is positive for *We* \> *We*~*C1*~.

Finally, the effect of boundary fluctuations on the viscoelastic flow characteristics in supercritical zone is investigated. According to [Fig. 9](#fig9){ref-type="fig"}(a) the velocity magnitude initially increased and after t = 10 a slight reduction to its steady profile at t = 20 occurred with velocity maximum located at y = 0.2.Fig. 9The variations of a) velocity, b) normal and c) shear stress as a function of y and t for subcritical zone with parameters of $\text{ε} = 0.04$, $\text{Re} = 0.2$, $\nabla\text{p} = - 1$, *We=15\*Sin(0.1t)* and, ζ = 0.2.Fig. 9

The normal stress variable is depicted in [Fig. 9](#fig9){ref-type="fig"}(b). Except of two instabilities between t = 5 and t = 10, the flow is almost stable with zero value of normal stress at the channel walls and maximum at y = 0.2. Regarding the shear stress in the supercritical range of *We*, after some pulsation for t = 5 to 10, its magnitude reaches the final value ([Fig. 9](#fig9){ref-type="fig"}-c).

For the supercritical case especially at high oscillations frequencies, the flow field becomes unstable. The reason of such behavior is due to the relative increase of the *We* in comparison to *ε* parameter, which in turn enhances the magnitude of non-linear terms and leads to instability of flow field.

3.4. The effect of Re on the transient behavior of flow parameters {#sec3.4}
------------------------------------------------------------------

The transient response of the flow parameters (*U*~*1*~, *N*~*1*~ and *S*~*1*~) for *We* = 0.15, *Δp* = -1, *ζ* = 0.2 and *ε* = 0.04 to the variations of *Re* number is presented in [Fig. 10](#fig10){ref-type="fig"}. For small *Re* the convergence to steady condition due to viscosity effect is uniform and rapid with several high amplitude oscillations. In the case with Re = 0.2 after three consecutive vibrations at t = 8 the variables become steady, however the mentioned steady state condition doesn\'t occur before t = 10 for Re = 4. The numerical data show the same trend for stress components. Therefore, the *Re* number determines the number of oscillations as well as convergence rate prior to the stable condition.Fig. 10Transient behavior of a) *U1*, b) *N1* and c) *S1*, for various *Re* with *We* = 0.15, *Δp* = -1, *ζ* = 0.2 and *ε* = 0.04.Fig. 10

3.5. The effect of *ε* on the transient behavior of flow parameters {#sec3.5}
-------------------------------------------------------------------

The effect of *ε* (the viscosity ratio of solvent to soluble) on the velocity, shear and normal stress is examined in the present section.

According to *U*~*1*~ and *U*~*2*~ curves which are obtained for *We* = 0.15 and *Re* = 0.2, by increasing the *ε* parameter from 0.04 to 0.1 the stable velocity magnitude as well as initial fluctuations and the time needed for the system to reach steady state is increased ([Fig. 11](#fig11){ref-type="fig"}).Fig. 11Transient behavior of a) *U*~*1*~ and b) *U*~*2*~, for various *ε* with *We* = 0.15, *Δp* = -1, *ζ* = 0.2 and *Re* = 0.2.Fig. 11

Besides, increasing of *ε* by the assumption of constant viscosity for the Newtonian solute (water), means a reduction in the viscosity of the Non-Newtonian soluble. Hence, the pointed out growths are meaningful as a results of reduction of solution viscosity and its lower capability in damping of the flow fluctuations.

Increasing of *ε* results in generation of higher normal stresses. Moreover, the maximum and minimum peaks at the initial time are intensified. The same observation can be detected for the shear stress ([Fig. 12](#fig12){ref-type="fig"}).Fig. 12Transient behavior of a) *N*~*1*~, b) *N*~*2*~, c) *S*~*1*~ and d) *S*~*2*~, for various *ε* with *We* = 0.15, *Δp* = -1, *ζ* = 0.2 and *Re* = 0.2.Fig. 12

4. Conclusion {#sec4}
=============

1)The examination of the governing equations of the viscoelastic flow with Johnson-Segalman model denotes that the flow behavior between two critical *We* numbers becomes unstable and the nonlinear velocity profile emerges.2)The inspection of the *Re* number effect implies that by enhancement of *Re*, the flow becomes stable for all of the simulated *We* and the early velocity and stress peaks diminish.3)In the subcritical zone, the velocity profile for non-oscillatory viscoelastic fluid is similar to the Newtonian fluid. However, the velocity and stress profiles overgrowing the stable profile at some early instances and damp to the final condition afterwards.4)Increasing of oscillations frequency in subcritical zone, the same as low frequency case, has almost no effect on the velocity and its gradient. Nevertheless, the normal stress amplitude of oscillations is reduced.5)For non-oscillatory, supercritical case in contrary to the subcritical situation in which the velocity maximum occurs at the midline, its location shifted toward the walls and the stabilization period increased.6)The *Re* number determines the number of oscillations and the needed time prior to the steady condition. For lower *Re*, due to higher effect of viscosity the initial fluctuations are intense and occurred in a short time period in contrary to the high *Re* case.

Declarations {#sec5}
============

Author contribution statement {#sec5.1}
-----------------------------

Reza Roohi, Nariman Ashrafi, Sepideh Samghani, Mohammad Najafi: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.

Funding statement {#sec5.2}
-----------------

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Competing interest statement {#sec5.3}
----------------------------

The authors declare no conflict of interest.

Additional information {#sec5.4}
----------------------

No additional information is available for this paper.
