A numerical study on nonlinear dynamics of oscillatory time-depended viscoelastic flow between infinite parallel plates: utilization of symmetric and antisymmetric Chandrasekhar functions

The one-dimensional viscoelastic fluid flow between two infinite parallel plates with oscillatory inlet condition is examined using the Johnson–Segalman model. The symmetric and antisymmetric Chandrasekhar functions in space are utilized to represent the velocity and stress fields. The non-dimensional form of the conservation laws in addition to the constitutive equations are solved numerically based on the Galerkin projection method. Two critical Weissenberg numbers (We) for various Reynolds numbers (Re) and viscosity ratios (ε) are obtained to determine the stable range of nonlinear system behavior. Moreover, for the unsteady case, the effects of Re, viscosity ratio of solvent to solution as well as We are investigated. According to the obtained results, 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. 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 intensely occurred in a short time period in contrary to the high Re case.


Introduction
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 [1,2]. In addition to the basic Maxwell model, several models based on the rheology theories are developed and implemented up to now (Bingham model [3] and Phan-Thien and Tanner (PTT) model [4]).
Moreover, Oldroyd-B model was used in [5] 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 [6,7] was investigated for low and high oscillation frequencies. The unsteady flow of blood through stenosed artery [8] 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 [9]. 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 [10]. 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) [11]. 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 [12].
In 2017, the effect of oscillatory vibration of a thin blade on the non-Newtonian fluid flow inside a channel is numerically examined [13]. 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 [14], 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. [15] 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 [16] tackled the problem of squeezing flow between parallel plates using the Homotopy analysis method. In another investigation, Hoshyar et al. [17] 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 [18]. 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.

The fluid flow laws
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 (η) is determined based on the superposition of Newtonian solvent viscosity (η 1 ) and the non-Newtonian polymer with viscosity of η 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.
Generally, the flow behavior is governed by the conservation laws of mass and linear momentum for an incompressible fluid.
Conservation of mass: 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: where λ is a dimensionless parameter (relaxation time), T is stress tensor, rU denotes the gradient of velocity vector and rU t is the transpose of the rU. Additionally, DT Dt can be calculated as: where ζ 2 ½0 À2 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. Taking U to be the maximum velocity in the flow direction, the nondimensional variables can be defined as: By substitution of these transformations into Eqs. (5), (6), and (7) the following equations are formed: The boundary condition is defined as: 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, ε ,and the oscillatory inlet velocity are defined as: Re ¼

The base flow governing equations
The steady state condition is determined via neglecting of time derivative terms in Eqs. (9), (10), (11), and (12). The steady state stress components are obtained as follows: where the A function is defined as: Eqs. (5), (6), and (7) can be rewritten by the definition of normal stresses difference (N) and a combination of normal stresses (Z) as: Moreover, to investigate the created disturbances in the flow, the following modifications are applied: For simplicity, the first term on the right-hand side of Eq. (9) is redefined as: Considering the relation between A and B parameters and the implementation of Eq. (23) modifications into Eqs. (9), (10), (11), and (12), the following governing equations are obtained: All of the unknown parameters in Eqs. (25), (26), and (27) are defined as:

Solution method
The solution of Eqs. (25), (26), and (27) 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), (34), and (35) as: where ϕ i ðyÞ is the even and odd Chandrasekhar functions, for i even and odd, respectively [12], . 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. u k and N k where k 2 ½1 À M): The coefficients in Eqs.  (41)). 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.

The bifurcation graph
In the present section the numerical solution of the governing equations (Eqs. (36), (37), (38), (39), (40), and (41)) 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 2 ½0 À10 and δ 2 ½0:2 À1:0 for limiting Re values of 0.2 and 10. It should be added that Δp and ε parameters are set to -1 and 0.04, respectively. According to Fig. 2 (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) (i.e. u', p' and τ ' ij ). As the value of these peaks is reduced the fluid resistance to the imposed forces and initial disturbances is decreased. Based on Fig. 2, three main zones can be defined as subcritical (We < We C1 ), critical (We C1 < We < We C2 ) and supercritical (We > We C2 Þ.
First and second critical We numbers are located in the range of ½0:2 À0:6 and ½0:4 À 0:8, respectively. Unstable base flow with nonlinear velocity profile is observed between these two critical We numbers.

Validation
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. [18], 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.  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 [18]. 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.

Transient simulation results
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). 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. The velocity profiles approximately resemble the Newtonian velocity profile, although a slight declination to the right is observed in the profiles.
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. 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.
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. 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).
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) and (39)) is the reason for such a behavior.
The effect of increasing of the oscillation frequency is presented in Fig. 7. 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(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).
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, that the velocity has a nonoscillatory 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), 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.
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(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.
The normal stress variable is depicted in Fig. 9(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-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.

The effect of Re on the transient behavior of flow parameters
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. 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.

The effect of ε on the transient behavior of flow parameters
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).
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).

Conclusion
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.

Author contribution statement
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.