On the Study of Viscoelastic Walters' B Fluid in Boundary Layer Flows

Abstract

Viscoelastic Walters' B fluid flows for three problems, stagnation-point flow, Blasius flow, and Sakiadis flow, have been investigated. In each problem, Cauchy equations are changed to a nondimensional differential equations using stream functions and with assumption of boundary layer flow. The fourth-order predictor-corrector finite-difference method for solving these nonlinear differential equations has been employed. The results that have been obtained using this method are compared with the results of the last studies, and it is clarified that this method is more accurate. It is also shown that the results of last study about Sakiadis flow of Walter's B fluid are not true. In addition, the effects of order of discretization in the boundaries are investigated. Moreover, it has been discussed about the valid region of Weissenberg numbers for the second-order approximation of viscoelastic fluids in each case of study.

1. Introduction

Boundary layer flows of non-Newtonian fluids can be a beneficial approach for modeling several manufacturing processes of industry such as the aerodynamic extrusion of plastic sheets, the cooling of metallic plates, fabrication of adhesive tapes, application of coating and layers onto rigid substrates. Some of articles that have investigated these kinds of fluid flows are given in [1–24]. These types of phenomenon can be considered in two branches, one about the geometry of flow and it’s boundary conditions and the other about the viscoelastic model. If we take a glance at the literature, we find that there are many models for considering fluids that have both elastic and viscous properties. One of these models is Walters’ B that has been used in several studies and for modeling many phenomena [25–33]. Study of the flow problems of this class of fluids not only is important technologically, but is also challenging to engineers and applied mathematicians who are interested in obtaining accurate solutions.

Constitutive equations of viscoelastic fluids usually generate higher-order derivative terms in the momentum equations that make them more difficult to solve in comparison with Newtonian fluids. There are some analytical and numerical methods for investigating these kinds of fluids such as homotopy analysis method [29] and spectral method [34]. Because of the apparent nonavailability of extra boundary conditions, researchers tend to develop a regular perturbation solution of the problem, taking the solution for the Newtonian fluid as the primary solution and the first-order perturbed solution as the secondary solution. A classical example of this technique is the analysis of the two-dimensional flow of Walters’ B fluid about a stagnation point given by Beard and Walters [11]. Perturbation approach can lead to erroneous results not only quantitatively but also qualitatively which is shown by some researchers such as Ariel [35]. He proposed a predictor-corrector method for solving linear and nonlinear differential equations such as boundary layer stagnation-point flow. This method can solve a differential equation by an explicit method with fourth order of approximation and without any restriction of convergence.

In this paper, similarity solutions of three viscoelastic Walters’ B boundary layer flow problems, stagnation-point flow, Blasius flow and Sakiadis flow, are obtained and all of the equations are changed to nondimensional forms in the second section. The mentioned fourth-order predictor-corrector method for solving these nonlinear differential equations is presented in the third section, and the results are compared with last studies in the fourth section. It is obtained that this method is more accurate and straightforward. It is also shown that the results of last article about Sakiadis flow by Sadeghy and Sharifi [36] are not true. Moreover, valid regions of Weissenberg numbers for the second-order model in each problem are obtained, and the effects of order of discretization in the boundaries are investigated.

2. Mathematical Formulation

In this section, nondimensional equations of motions for three different boundary layer flows, stagnation-point flow, Blasius flow, and Sakiadis flow, are obtained. These equations are solved in Section 4 using a predictor-corrector method.

2.1. Stagnation-Point Flow

Cauchy equations are employed for obtaining boundary layer equation of two-dimensional stagnation-point flow. Steady Cauchy equations areContinuity:𝜕𝑢+𝜕𝑥𝜕𝑣𝑢𝜕𝑦=0,𝑥-momentum:𝜌𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦=−𝜕𝑝+𝜕𝑥𝜕𝜏_𝑥𝑥+𝜕𝑥𝜕𝜏_𝑥𝑦,𝑢𝜕𝑦𝑦-momentum:𝜌𝜕𝑣𝜕𝑥+𝑣𝜕𝑣𝜕𝑦=−𝜕𝑝+𝜕𝑦𝜕𝜏_𝑥𝑦+𝜕𝑥𝜕𝜏_𝑦𝑦,𝜕𝑦 (2.1) where 𝑝 is pressure, 𝜌 is density of fluid, and 𝑢 and 𝑣 are velocity components in 𝑥 and 𝑦 directions, respectively. Also, 𝜏_𝑥𝑥, 𝜏_𝑥𝑦, 𝜏_𝑦𝑥, and 𝜏_𝑦𝑦 are components of stress matrix. These variables can be obtained for Walters-B viscoelastic fluid by following equation:𝜏=2𝜇₀𝑑−2𝑘𝑑₀^∇ (2.2) which is the constitutive equation of Walters’ B fluid and 𝜇₀ is the dynamic viscosity, 𝑘₀ is the short-memory coefficient, and 𝑑 is the deformation rate tensor. Also 𝑑^∇ is upper-convected derivative of a tensor and can be obtained from this equation:𝑑^∇=𝜕𝑑+()𝜕𝑡𝑣⃗⋅∇𝑑−𝑑⋅∇𝑉^𝑇−∇𝑉⋅𝑑. (2.3) So, the stress components in steady conditions are𝜏_𝑖𝑗=2𝜇₀𝑑_𝑖𝑗−2𝑘₀𝐷𝑑_𝑖𝑗−𝐷𝑡𝜕𝑣_𝑗𝜕𝑥_𝑚𝑑_𝑖𝑚−𝜕𝑣_𝑖𝜕𝑥_𝑚𝑑_𝑚𝑗, (2.4) where 𝑑_𝑖𝑗 can be obtained as follows:𝑑_𝑖𝑗=12𝜕𝑉_𝑖𝜕𝑥_𝑗+𝜕𝑉_𝑗𝜕𝑥_𝑖. (2.5) On the basis of the above relation, the stress components can be written as follows:𝜏_𝑥𝑥=2𝜇₀𝜕𝑢𝜕𝑥−2𝑘₀𝑢𝜕²𝑢𝜕𝑥²𝜕+𝑣²𝑢𝜕𝑥𝜕𝑦−2𝜕𝑢𝜕𝑥²+12𝜕𝑢𝜕𝑦𝜕𝑢+𝜕𝑦𝜕𝑣,𝜏𝜕𝑥_𝑦𝑥=𝜏_𝑥𝑦=𝜇₀𝜕𝑢+𝜕𝑦𝜕𝑣𝜕𝑥−2𝑘₀12𝑢𝜕²𝑢+𝜕𝜕𝑥𝜕𝑦²𝑣𝜕𝑥²+12𝑣𝜕²𝑢𝜕𝑦²+𝜕²𝑣−𝜕𝑥𝜕𝑦𝜕𝑢𝜕𝑥𝜕𝑣+𝜕𝑥𝜕𝑢𝜕𝑦𝜕𝑣,𝜏𝜕𝑦_𝑦𝑦=2𝜇₀𝜕𝑣𝜕𝑦−2𝑘₀𝜕²𝑣𝜕𝜕𝑡𝜕𝑥+𝑢²𝑣𝜕𝜕𝑥𝜕𝑦+𝑣²𝑣𝜕𝑦²1−22𝜕𝑣𝜕𝑥𝜕𝑢+𝜕𝑦𝜕𝑣+𝜕𝑥𝜕𝑣𝜕𝑦². (2.6)

These equations are the relation of components of stress tensor and velocity components. Because of requirement for derivatives of these stress components, these derivatives are presented.

By substituting these relations in (2.1), Cauchy equations become

𝑥-momentum:𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢1𝜕𝑦=−𝜌𝜕𝑝+𝜇𝜕𝑥𝜌⁰𝜕²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²−𝑘₀𝜌𝑢𝜕³𝑢𝜕𝑥³+𝜕³𝑢𝜕𝑥𝜕𝑦²𝜕+𝑣³𝑢𝜕𝑥²+𝜕𝜕𝑦³𝑢𝜕𝑦³−𝜕𝑢3𝜕𝜕𝑥²𝑢𝜕𝑥²−𝜕²𝑢𝜕𝑦²−𝜕𝑢𝜕𝜕𝑦²𝑢+𝜕𝜕𝑥𝜕𝑦²𝑣𝜕𝑥²−2𝜕𝑣𝜕𝜕𝑥²𝑢,𝜕𝑥𝜕𝑦 (2.7)𝑦-momentum:𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣1𝜕𝑦=−𝜌𝜕𝑝+𝜇𝜕𝑦𝜌⁰𝜕²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²−𝑘₀𝜌𝑢𝜕³𝑣𝜕𝑥³+𝜕³𝑣𝜕𝑥𝜕𝑦²𝜕+𝑣³𝑣𝜕𝑥²+𝜕𝜕𝑦³𝑣𝜕𝑦³−𝜕𝑣3𝜕𝜕𝑦²𝑣𝜕𝑦²−𝜕²𝑣𝜕𝑥²−𝜕𝑣𝜕𝜕𝑥²𝑣+𝜕𝜕𝑥𝜕𝑦²𝑢𝜕𝑦²−2𝜕𝑢𝜕𝜕𝑦²𝑣.𝜕𝑥𝜕𝑦 (2.8) The 𝑥- and 𝑦-momentum equations can be simplified as follows after employing boundary layer approximation:𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢1𝜕𝑦=−𝜌𝑑𝑝+𝜇𝑑𝑥𝜌⁰𝜕²𝑢𝜕𝑦²−𝑘₀𝜌𝑢𝜕³𝑢𝜕𝑥𝜕𝑦²𝜕+𝑣³𝑢𝜕𝑦³+𝜕𝑢𝜕𝜕𝑥²𝑢𝜕𝑦²−𝜕𝑢𝜕𝜕𝑦²𝑢𝜕𝑥𝜕𝑦𝜕𝑝𝜕𝑦=0. (2.9)

On the basis of viscous flow theory, the change in pressure across boundary layer is 𝑂(𝛿²). Therefore, the pressure gradient term in (2.9) can be obtained from the flow just outside the boundary layer. So, following equation is obtained in this region and under the steady flow condition:−1𝜌𝜕𝑝𝜕𝑥=𝑈_∞𝑑𝑈_∞,𝑑𝑥 (2.10) where 𝑈_∞ is the mainstream velocity. So, boundary layer equation for this fluid flow is obtained as follows:𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦=𝑈_∞𝑑𝑈_∞𝜕𝑑𝑥+𝜐²𝑢𝜕𝑦²−𝑘₀𝜌𝑢𝜕³𝑢𝜕𝑥𝜕𝑦²𝜕+𝑣³𝑢𝜕𝑦³+𝜕𝑢𝜕𝜕𝑥²𝑢𝜕𝑦²−𝜕𝑢𝜕𝜕𝑦²𝑢.𝜕𝑥𝜕𝑦 (2.11) Because of the similarity of the velocity profiles in different cross-sections, it is concluded that the only possible form of 𝑈_∞ under these conditions for obtaining the solution of (2.11) is𝑈_∞=𝜀𝑥̇, (2.12) where 𝜀̇ is a constant parameter [37]. This procedure has been described in [37] completely.

The 𝑥-momentum equation can be nondimensionalized using stream function Ψ. So, the velocity components can be written as𝑢=𝜕Ψ𝜕𝑦,𝑣=−𝜕Ψ.𝜕𝑥 (2.13) The stream function can be made dimensionless as follows:Ψ𝑥,𝑦=√𝜐𝜀𝑥̇𝑓𝜂, (2.14) where 𝑓(𝜂)𝑓 is the dimensionless stream function, 𝜐 is the kinematic viscosity, and 𝜂 is the similarity variable, defined by𝜂𝑥,𝑦=𝜀̇𝜐𝑦. (2.15) By substituting (2.11)–(2.15) into (2.9), dimensionless equation can be obtained:𝑓^′′′+𝑓𝑓^′′+1−𝑓^′2+𝐾𝑓𝑓^𝑖𝑣−2𝑓^′𝑓^′′′+𝑓^2′′=0, (2.16) where 𝐾 is the local Weissenberg number defined by𝑘𝐾=⁰𝜀̇𝜇₀. (2.17) On the basis of no slip condition, no permeability, and 𝑢(𝑦→∞)=𝜀𝑥̇, three boundary conditions of (2.16) are𝑓(0)=0,𝑓^′(0)=0,𝑓^′(∞)=1. (2.18)

Order of dimensionless equation of motion, (2.16), is the fourth while the momentum equations are the third order. This is because of derivative operations that are applied to momentum equations for obtaining dimensionless equation (2.16).

In addition to velocities and stream functions that are common parameters for investigating boundary layer problems, shear stress as another important parameter should obtained in the wall. Therefore, the following equation for shear is presented:𝜏_𝑥𝑦=𝜇₀𝜕𝑢𝜕𝑦−𝑘₀𝑢𝜕²𝑢𝜕𝜕𝑥𝜕𝑦+𝑣²𝑢𝜕𝑦²+2𝜕𝑢𝜕𝑥𝜕𝑢.𝜕𝑦 (2.19) By substituting 𝑢 and 𝑣 on the basis of (2.13) and (2.14) on (2.19), shear stress on the solid boundary can be obtained using (2.18) as follows:𝜏_𝑥𝑦∣_𝑦=0=𝜀̇^3/2𝑥𝜇₀√𝜐𝑓^′′(0). (2.20)

2.2. Blasius Flow

Equation of motion for Blasius flow can be obtained using (2.1) that leads to (2.9).

The 𝑥-momentum equation of (2.9) can be nondimensionalized using stream function Ψ. So, the velocity components can be written as𝑢=𝜕Ψ𝜕𝑦,𝑣=−𝜕Ψ.𝜕𝑥 (2.21) The stream function can be made dimensionless as follows:Ψ𝑥,𝑦=𝜐𝑥𝑈_∞𝑓𝜂, (2.22) where 𝑓 is the dimensionless stream function, 𝜐 is the kinematic viscosity, and 𝜂 is the similarity variable, defined by𝜂𝑥,𝑦=𝑈_∞𝜐𝑥𝑦. (2.23) So, nondimensional equation of motion for Blasius flow can be obtained on the basis of 𝑑𝑝/𝑑𝑥=0 as follows:𝑓^′′′+12𝑓𝑓^′′+𝐾𝑓𝑓^𝑖𝑣+2𝑓^′𝑓^′′′−𝑓^2′′=0, (2.24) where 𝐾 is the local Weissenberg number defined by𝐾(𝑥)=𝑘₀𝑈_∞2𝜇₀𝑥. (2.25) Boundary conditions of (2.24) are𝑓(0)=0,𝑓^′(0)=0,𝑓^′(∞)=1. (2.26) By substituting 𝑢 and 𝑣 from (2.21) and (2.22) on (2.19) and employing boundary conditions that are presented in (2.26), shear stress on the solid boundary can be obtained as follows:𝜏_𝑥𝑦∣_𝑦=0=𝜇₀𝑈_∞𝑈_∞𝑓𝜐𝑥^′′(0). (2.27)

2.3. Sakiadis Flow

Sakiadis flow for viscoelastic Walters’ B fluid has been investigated by Sadeghy and Sharifi [36] at 2004. The results that are presented in [36] are not true and reliable which are based on false nondimensional equation of motion that is derived in [36]. In this paper, the true equation of motion is presented and reliable results are obtained.

Sakiadis flow is similar to Blasius flow and the only difference of these fluid flows is their boundary conditions that are presented in following relations:𝑓(0)=0,𝑓^′(0)=1,𝑓^′(∞)=0. (2.28)

In the next section, the predictor-corrector method for analysis of nondimensional equations of these flows is described.

In addition, shear stress on the solid boundary can be obtained similar to (2.27) and using boundary conditions that are presented in (2.28) as follows:𝜏_𝑥𝑦∣_𝑦=0=(1−𝐾)𝜇₀𝑈_∞𝑈_∞𝑓𝜐𝑥^′′(0). (2.29)

3. Method of Solution

In this section, fourth-order predictor-corrector method that was proposed by Ariel [35] has been described. It should be noted that stagnation-point flow has been investigated by Ariel [35] using this method. So, solution of this flow is just obtained to show the true using of this approach.

3.1. Two-Dimensional Stagnation-Point Flow

It should be noted that because of the singularity of (2.16) and (2.19) at 𝜂=0, some usual methods, such as shooting method, cannot be employed for solving these equations. Henceforth, the predictor-corrector method, which not only is an explicit method without requiring for solving any matrix but also converges without any convergence condition, is used in this study. For clarity, this method is described by applying to (2.16) that is the nondimensional equation of motion for stagnation-point flow.

The auxiliary parameters𝑓^′=𝑝,𝑝^′=𝑄 (3.1) should be used to rewritte (2.16) as follows:𝑘𝑓𝑄^′′+1−2𝑘𝑝𝑄^′+𝑓𝑄+1−𝑝²+𝑘𝑄²=0, (3.2) and boundary conditions (2.18) become𝑓(0)=0,𝑝(0)=0,𝑝(∞)=1. (3.3)

The key of the algorithm is to retrain the second-order derivative in (3.2) and not to further decompose it into a pair of differential equations in 𝑄 and 𝑄′.

The following equations should be used for approximating the first- and second-order derivatives:𝐷𝑦_𝑖=1ℎ1∇−2∇²−16∇³−⋯𝑦_𝑖+1,𝐷²𝑦_𝑖=1ℎ²∇²−12∇⁴−16∇⁵−⋯𝑦_𝑖+1, (3.4) where 𝐷 is the derivative operator and ∇ is the backward operator.

These operators are employed for discretizing (3.1) and (3.2) which lead to the following equations:1ℎ14𝑓_𝑖+1+56𝑓_𝑖−32𝑓_𝑖−1+12𝑓_𝑖−2−1𝑓12^𝑖−3=𝑝_𝑖,1ℎ14𝑝_𝑖+1+56𝑝_𝑖−32𝑝_𝑖−1+12𝑝_𝑖−2−1𝑝12^𝑖−3=𝑄_𝑖,𝑘𝑓_𝑖ℎ56𝑄_𝑖+1+54𝑄_𝑖−13𝑄_𝑖−1+76𝑄_𝑖−2−12𝑄_𝑖−3+1𝑄12^𝑖−4+1−2𝑘𝑝_𝑖ℎ14𝑄_𝑖+1+56𝑄_𝑖−32𝑄_𝑖−1+12𝑄_𝑖−2−1𝑄12^𝑖−3+𝑓_𝑖𝑄_𝑖+1−𝑝²_𝑖+𝑘𝑄²_𝑖=0. (3.5) Equation (3.5) can be explicitly solved for 𝑓,𝑝,and𝑄 at the (𝑖+1)th point in a supposing mesh generation. Using the subscript “𝑝” for the predicted values, we can write𝑓^(𝑝)_𝑖+1=13−10𝑓_𝑖+18𝑓_𝑖−1−6𝑓_𝑖−2+𝑓_𝑖−3+12ℎ𝑝_𝑖,𝑝^(𝑝)_𝑖+1=13−10𝑝_𝑖+18𝑝_𝑖−1−6𝑝_𝑖−2+𝑝_𝑖−3+12ℎ𝑄_𝑖,𝑄^(𝑝)_𝑖+1=10𝑘𝑓_𝑖+3ℎ1−2𝑘𝑝_𝑖⁻¹𝑘𝑓_𝑖15𝑄_𝑖+4𝑄_𝑖−1−14𝑄_𝑖−2+6𝑄_𝑖−3−𝑄_𝑖−4+ℎ1−2𝑘𝑝_𝑖−10𝑄_𝑖+18𝑄_𝑖−1−6𝑄_𝑖−2+𝑄_𝑖−3−12ℎ²𝑓_𝑖𝑄_𝑖+1−𝑝²_𝑖+𝑘𝑄²_𝑖. (3.6) These steps were about the predicting the values of 𝑓_𝑖+1,𝑝_𝑖+1,and𝑄_𝑖+1. So, the following operators are necessary to obtain the correct values of these variables:𝐷𝑦_𝑖=1ℎ1∇+2∇²+13∇³+14∇⁴⋯𝑦_𝑖,𝐷²𝑦_𝑖=1ℎ²∇²+∇³+11∇12⁴+56∇⁵+⋯𝑦_𝑖. (3.7) On the basis of (3.7), (3.1) and (3.2) become𝑓_𝑖+1=12548𝑓_𝑖−36𝑓_𝑖−1+16𝑓_𝑖−2−3𝑓_𝑖−3+12ℎ𝑝^(𝑝)_𝑖+1,𝑝_𝑖+1=12548𝑝_𝑖−36𝑝_𝑖−1+16𝑝_𝑖−2−3𝑝_𝑖−3+12ℎ𝑄^(𝑝)_𝑖+1,𝑄_𝑖+1=45𝑘𝑓^(𝑝)_𝑖+1+25ℎ1−2𝑘𝑝^(𝑝)_𝑖+1⁻¹𝑘𝑓^(𝑝)_𝑖+1154𝑄_𝑖−214𝑄_𝑖−1+156𝑄_𝑖−2−61𝑄_𝑖−3+10𝑄_𝑖−4+ℎ1−2𝑘𝑝^(𝑝)_𝑖+148𝑄−36𝑄_𝑖−1+16𝑄−3𝑄_𝑖−3−12ℎ²𝑓^(𝑝)_𝑖+1𝑄^(𝑝)_𝑖+1+1−𝑝^(𝑝)2_𝑖+1+𝑘𝑄^(𝑝)2_𝑖+1. (3.8)

Therefore, we have a fourth-order predictor-corrector method for obtaining the values of 𝑓,𝑓^′ and 𝑓^′′ at a particular point based with respect to their values at the previous five mesh points. The values of 𝑓 and 𝑓^′are given at 𝜂=0. Since 𝑓^′ is given at infinity, (2.16) can be thought as a boundary value problem (BVP). Therefore, we can solve these equations by changing the initial guess of 𝑓^′′ so that the boundary condition at infinity is satisfied.

Assume that the amount of the initial guess of 𝑓^′′ is𝑓^′′(0)=𝑠. (3.9)

The amounts of in 𝑓,𝑝,and𝑄 at five extra nodes are needed to start the solution of the discretized equation. So, by applying the Taylor series around zero, we can find the amount of these parameters in 𝜂=±ℎ𝜂=±2ℎ. Taylor series for expansion of of 𝑓 around 𝜂=0 results in the following equation:𝑓𝜂=⁵_𝑖=0𝑐_𝑖𝜂^𝑖, (3.10) where 𝑐_𝑖 are constant and should be obtained using the boundary values of 𝑓,𝑝,and𝑄 at 𝜂=0. Henceforth, some of these constants are𝑐₀=0,𝑐₁=0,𝑐₂=12𝑠,𝑐₃=−𝐾𝑠²+16,𝑐₄=0,𝑐₅=𝑠²,120 (3.11) where 𝑐₀,𝑐₁,and𝑐₂ are obtained on the basis of boundary values of parameters at 𝜂=0 and others are obtained using (2.16).

3.2. Blasius Flow

Nondimensional equation of Blasius flow (2.19) can be solved by the proposed fourth-order predictor-corrector method. So, the predictor and corrector relations are obtained as follows.

Predictor:𝑓^(𝑝)_𝑖+1=13−10𝑓_𝑖+18𝑓_𝑖−1−6𝑓_𝑖−2+𝑓_𝑖−3+12ℎ𝑝_𝑖,𝑝^(𝑝)_𝑖+1=13−10𝑝_𝑖+18𝑝_𝑖−1−6𝑝_𝑖−2+𝑝_𝑖−3+12ℎ𝑄_𝑖,𝑄^(𝑝)_𝑖+1=10𝑘𝑓_𝑖+3ℎ1+2𝑘𝑝_𝑖⁻¹𝑘𝑓_𝑖15𝑄_𝑖+4𝑄_𝑖−1−14𝑄_𝑖−2+6𝑄_𝑖−3−𝑄_𝑖−4+ℎ1+2𝑘𝑝_𝑖−10𝑄_𝑖+18𝑄_𝑖−1−6𝑄_𝑖−2+𝑄_𝑖−3−12ℎ²12𝑓_𝑖𝑄_𝑖−𝑘𝑄²_𝑖. (3.12) Corrector:𝑓_𝑖+1=12548𝑓_𝑖−36𝑓_𝑖−1+16𝑓_𝑖−2−3𝑓_𝑖−3+12ℎ𝑝^(𝑝)_𝑖+1,𝑝_𝑖+1=12548𝑝_𝑖−36𝑝_𝑖−1+16𝑝_𝑖−2−3𝑝_𝑖−3+12ℎ𝑄^(𝑝)_𝑖+1,𝑄_𝑖+1=45𝑘𝑓^(𝑝)_𝑖+1+25ℎ1+2𝑘𝑝^(𝑝)_𝑖+1⁻¹𝑘𝑓^(𝑝)_𝑖+1154𝑄_𝑖−214𝑄_𝑖−1+156𝑄_𝑖−2−61𝑄_𝑖−3+10𝑄_𝑖−4+ℎ1+2𝑘𝑝^(𝑝)_𝑖+148𝑄−36𝑄_𝑖−1+16𝑄−3𝑄_𝑖−3−12ℎ²12𝑓^(𝑝)_𝑖+1𝑄^(𝑝)_𝑖+1−𝑘𝑄^(𝑝)2_𝑖+1. (3.13) It is obvious that these equations can be solved explicitly.

Quantities of extra nodes that are necessary for solving above equations are obtained similar to those obtained for stagnation-point flow using Taylor series. So, we have𝑓𝜂=⁵_𝑖=0𝑐_𝑖𝜂^𝑖, (3.14) where𝑐₀=0,𝑐₁=0,𝑐₂=12𝑠,𝑐₃=−𝐾𝑠²6,𝑐₄=0,𝑐₅𝑠=−²,240 (3.15) where 𝑐₀,𝑐₁,and𝑐₂ are obtained by employing the boundary values of parameters at 𝜂=0 and others are obtained using (2.19).

3.3. Sakiadis Flow

It is mentioned that the only difference between Sakiadis flow and Blasius flow is just in their boundary conditions. So, their predictor and corrector relations are same and we can use those of Blasius flow for solving (2.19) which is the nondimensional equation of motion of Sakiadis flow. On the other hand, the constants 𝑐₀,𝑐₁, and so forth, are different from Blasius flow that using (2.19) and (2.21) are obtained as𝑐₀=0,𝑐₁=1,𝑐₂=12𝑠,𝑐₃=𝐾𝑠²6(2𝐾+1),𝑐₄𝑠=−(),𝑐483𝐾+1₅=6𝐾²+4𝐾+1().2402𝐾+1)(3𝐾+1)(4𝐾+1 (3.16) In the next section, the predictor-corrector method that is described in this section is applied to the mentioned problems, stagnation-point flow, Blasius flow, and Sakiadis flow. Therefore, the results are obtained and compared with the results of the other papers that are obtained by different methods. In addition, the results presented by Sadeghy and Sharifi [36] in 2004 are reformed.

4. Result and Discussion

4.1. Stagnation-Point Flow

The results of the solution of stagnation-point flow with the proposed fourth-order method are presented in Table 1.

Table 1 

Variation of 𝑓^′′(0) in two-dimensional stagnation-point flow (𝜂_max=20).

On the basis of Table 1, it is obvious that step size (ℎ) must be smaller than 0.05 in this problem to obtain the reliable values of 𝑓^′′(0) which is the key point of the solution of this problem. In addition, larger values 𝜂_max cause some errors in the results for smaller value of Weissenberg number which can be omitted using smaller values of 𝜂_max which are chosen as 10 in these situations. Because of these errors, far from the wall, the numerical error becomes larger and makes an unstable numerical procedure. Thus, for 𝐾≤0.05 we used 𝜂_max=10, and for greater values of this viscoelastic parameter, 𝜂_max has been set equal to 20 to obtain the suitable results.

Moreover, the comparison of values of shear stress at solid boundary can be beneficial in this problem. On the basis of (2.20), this parameter in a particular location of the domain depends linearly on 𝑓^′′(0). So, it is obvious that, for larger values of 𝑓^′′(0), shear stress on the solid boundary is larger than for smaller ones. Therefore, it can be concluded from Table 1 that this parameter increases as Weissenberg number increases. It should be noted that valid values of 𝑓^′′(0) are based on step sizes smaller than 0.05.

The boundary layers in stagnation-point flow for different Weissenberg numbers are compared in Figure 1. As Weissenberg number increases, the boundary layer thickness decreases. In addition, Figure 1 clarifies the differences between the quantities of nondimensional horizontal velocity at any particular nondimensional parameter 𝜂. On the basis of this figure, values of nondimensional horizontal velocity increases as Weissenberg number increases.

Figure 1 

Nondimensional horizontal velocity of stagnation-point flow with respect to nondimensional parameter 𝜂.

As we see in the literature [38, 39], this ODE equation has been solved with other numerical scheme, perturbation method, keller-box method, and so forth, nonphysical phenomena, which are overshoot in the boundary layer, and occurred for 𝐾≥0.2. It is clarified that Walters’ B model and its constitutive equation are not suitable for this range of Weissenberg number. Figure 2 shows stream lines for different values of Weissenberg number. The differences of boundary layer thicknesses for different values of Weissenberg numbers are clarified in this figure.

Figure 2 

Effect of Weissenberg number of the stream line near stagnation point (solid line: 𝑘=0, dashed line: 𝑘=0.05, long-dashed line: 𝑘=0.1).

4.2. Blasius Flow

Quantities of 𝑓^′′(0) for different Weissenberg numbers are obtained, for ℎ=0.02, as shown in Table 2. It is shown that increasing the amount of 𝐾 causes to reduce the amount of 𝑓^′′(0). Therefore, on the basis of (2.27), shear stress at the solid boundary is smaller for larger values of Weissenberg number. The amount of 𝜂_max for all values of 𝐾 is opted 20, that is, similar to the solution of stagnation-point flow.

Table 2 

Variation of 𝑓^′′(0) in Blasius flow (𝜂_max=20,ℎ=0.02).

Figure 3 shows the growth of boundary layers in Blasius flow for different Weissenberg numbers. As 𝐾 increases, the boundary layer thickness decreases. It should be noted that a nonphysical phenomenon occurred for larger values of Weissenberg number, 𝐾≥0.3, which clarifies that this model is useless for investigating this problem in this this range of Weissenberg number.

Figure 3 

Nondimensional horizontal velocity of Blasius flow with respect to nondimensional parameter 𝜂.

The differences between the values of nondimensional horizontal velocity at any particular nondimensional parameter 𝜂 is shown in Figure 3. So, it is obvious that this variable decreases as Weissenberg number increases.

4.3. Sakiadis Flow

Variation of 𝑓^′′(0) in two-dimensional Sakiadis flow is clarified in Table 3. So, it is obvious that values of 𝑓^′′(0) become larger as Weissenberg number increases. Therefore, similar to Blasius flow and on the basis of (2.29), shear stress at the solid boundary is smaller for larger values of Weissenberg number with respect to smaller values.

Table 3 

Variation of 𝑓^′′(0) in Sakiadis flow (𝜂_max=20,ℎ=0.02).

The most important result of this study is shown in Figure 4. This figure that presents the growth of boundary layer of Sakiadis flow clarifies that results of Sadeghy and Sharifi [36] are unsatisfactory and clearly false. It is obvious that, for larger values of Weissenberg number, more than 0.3, Walters’ B model is not a reliable model for investigating the viscoelastic Sakiadis flow. But, Sadeghy and Sharifi [36] had concluded that this model is satisfactory for larger values even more than 1.

Figure 4 

Nondimensional horizontal velocity of Sakiadis flow with respect to nondimensional parameter 𝜂.

4.4. Order of Approximation of Boundary Conditions

One of the fundamental principles of numerical methods to obtain an arbitrary order of approximation is the sameness of order of discretization both in the entire domain and in the boundaries. The predictor-corrector method that is used in this study is a the fourth order method. So, boundary condition also should be approximated in the fourth-order. Amounts of 𝑓^′′(0) for 2nd and 4th order of approximation for calculating boundary conditions in stagnation point, Blasius, and Sakiadis flows are clarified in Tables 4, 5, and 6. On the basis of these tables, the effect of the approximation of boundary conditions is clarified.

Table 4 

Comparison of 𝑓^′′(0) for 2nd and 4th order of approximation of boundary conditions in stagnation-point flow (ℎ=0.01).

Table 5 

Comparison of 𝑓^′′(0) for 2nd and 4th order of approximation of boundary conditions in Blasius flow (ℎ=0.02).

Table 6 

Comparison of 𝑓^′′(0) for 2nd and 4th order of approximation of boundary conditions in Sakiadis flow (ℎ=0.02).

It is obvious that the amounts of 𝑓^′′(0) for all three problems are more accurate by the 4th order of approximation in boundary conditions.

5. Conclusion

In this study, the fourth-order predictor-corrector method is employed for solving the three viscoelastic fluid flow problems. Equations of motion for these problems, stagnation-point flow, Blasius flow, and Sakiadis flow, are obtained using Cauchy equation of motion, with assumption of boundary layer flow and turned into nondimensional forms. The discretized equations on the approach of the predictor-corrector method are solved, and the results are compared with the results of last studies, that is, clarified the high accuracy of the method. Also it has been shown that as elasticity increases, the stress on the solid boundary increases in stagnation-flow problem and decreases in Blasius and Sakiadis flows. In addition, it is shown that the results of the last study about Sakiadis flow by Sadeghy and Sharifi [36] are not true. Moreover, it has been shown that the effect of order of discretizations in the boundaries is an important issue. Finally, it is shown that there are restrictions for valid region of Weissenberg numbers and the second-order method of viscoelastic fluids failed out of these ranges.