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

Authors: Seyed Ali Madani Tonekaboni; Ramin Abkar; Reza Khoeilar
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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 Section4 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 are(2.1)Continuity:∂u∂x+∂v∂y=0,x-momentum:ρ(u∂u∂x+v∂u∂y)=-∂p∂x+∂τxx∂x+∂τxy∂y,y-momentum:ρ(u∂v∂x+v∂v∂y)=-∂p∂y+∂τxy∂x+∂τyy∂y,
where p is pressure, ρ is density of fluid, and u and v are velocity components in x and y directions, respectively. Also,  τxx, τxy, τyx, and τyy are components of stress matrix. These variables can be obtained for Walters-B viscoelastic fluid by following equation:(2.2)τ=2μ0d-2k0d∇
which is the constitutive equation of Walters’ B fluid and μ0 is the dynamic viscosity, k0 is the short-memory coefficient, and d is the deformation rate tensor. Also d∇ is upper-convected derivative of a tensor and can be obtained from this equation:(2.3)d∇=∂d∂t+(v⃗⋅∇)d-d⋅∇VT-∇V⋅d.
So, the stress components in steady conditions are(2.4)τij=2μ0dij-2k0(DdijDt-∂vj∂xmdim-∂vi∂xmdmj),
where dij can be obtained as follows:(2.5)dij=12[∂Vi∂xj+∂Vj∂xi].
On the basis of the above relation, the stress components can be written as follows:(2.6)τxx=2μ0∂u∂x-2k0[u∂2u∂x2+v∂2u∂x∂y-2{(∂u∂x)2+12∂u∂y(∂u∂y+∂v∂x)}],τyx=τxy=μ0(∂u∂y+∂v∂x)-2k0[12u(∂2u∂x∂y+∂2v∂x2)+12v(∂2u∂y2+∂2v∂x∂y)-{∂u∂x∂v∂x+∂u∂y∂v∂y}],τyy=2μ0∂v∂y-2k0[∂2v∂t∂x+u∂2v∂x∂y+v∂2v∂y2-2{12∂v∂x(∂u∂y+∂v∂x)+(∂v∂y)2}].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 becomex-momentum:(2.7)u∂u∂x+v∂u∂y=-1ρ∂p∂x+μ0ρ(∂2u∂x2+∂2u∂y2)-k0ρ[u(∂3u∂x3+∂3u∂x∂y2)+v(∂3u∂x2∂y+∂3u∂y3)-∂u∂x(3∂2u∂x2-∂2u∂y2)-∂u∂y(∂2u∂x∂y+∂2v∂x2)-2∂v∂x∂2u∂x∂y],y-momentum:(2.8)u∂v∂x+v∂v∂y=-1ρ∂p∂y+μ0ρ(∂2v∂x2+∂2v∂y2)-k0ρ[u(∂3v∂x3+∂3v∂x∂y2)+v(∂3v∂x2∂y+∂3v∂y3)-∂v∂y(3∂2v∂y2-∂2v∂x2)-∂v∂x(∂2v∂x∂y+∂2u∂y2)-2∂u∂y∂2v∂x∂y].
The x- and y-momentum equations can be simplified as follows after employing boundary layer approximation:(2.9)u∂u∂x+v∂u∂y=-1ρdpdx+μ0ρ(∂2u∂y2)-k0ρ[u∂3u∂x∂y2+v∂3u∂y3+∂u∂x∂2u∂y2-∂u∂y∂2u∂x∂y]∂p∂y=0.On the basis of viscous flow theory, the change in pressure across boundary layer isO(δ2). 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:(2.10)-1ρ∂p∂x=U∞dU∞dx,
where U∞ is the mainstream velocity. So, boundary layer equation for this fluid flow is obtained as follows:(2.11)u∂u∂x+v∂u∂y=U∞dU∞dx+υ(∂2u∂y2)-k0ρ[u∂3u∂x∂y2+v∂3u∂y3+∂u∂x∂2u∂y2-∂u∂y∂2u∂x∂y].
Because of the similarity of the velocity profiles in different cross-sections, it is concluded that the only possible form of U∞ under these conditions for obtaining the solution of (2.11) is(2.12)U∞=ε̇x,
where ε̇ is a constant parameter [37]. This procedure has been described in [37] completely.Thex-momentum equation can be nondimensionalized using stream function Ψ. So, the velocity components can be written as(2.13)u=∂Ψ∂y,v=-∂Ψ∂x.
The stream function can be made dimensionless as follows:(2.14)Ψ(x,y)=υε̇xf(η),
where f(η)f is the dimensionless stream function, υ is the kinematic viscosity, and η is the similarity variable, defined by(2.15)η(x,y)=(ε̇υ)y.
By substituting (2.11)–(2.15) into (2.9), dimensionless equation can be obtained:(2.16)f′′′+ff′′+1-f′2+K(ffiv-2f′f′′′+f′′2)=0,
where K is the local Weissenberg number defined by(2.17)K=k0ε̇μ0.
On the basis of no slip condition, no permeability, and u(y→∞)=ε̇x, three boundary conditions of (2.16) are(2.18)f(0)=0,f′(0)=0,f′(∞)=1.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.19)τxy=μ0∂u∂y-k0(u∂2u∂x∂y+v∂2u∂y2+2∂u∂x∂u∂y).
By substituting u and v 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:(2.20)τxy∣y=0=ε̇3/2xμ0υf′′(0).


### 2.2. Blasius Flow

Equation of motion for Blasius flow can be obtained using (2.1) that leads to (2.9).Thex-momentum equation of (2.9) can be nondimensionalized using stream function Ψ. So, the velocity components can be written as(2.21)u=∂Ψ∂y,v=-∂Ψ∂x.
The stream function can be made dimensionless as follows:(2.22)Ψ(x,y)=υxU∞f(η),
where f is the dimensionless stream function, υ is the kinematic viscosity, and η is the similarity variable, defined by(2.23)η(x,y)=U∞υxy.
So, nondimensional equation of motion for Blasius flow can be obtained on the basis of dp/dx=0 as follows:(2.24)f′′′+12ff′′+K(ffiv+2f′f′′′-f′′2)=0,
where K is the local Weissenberg number defined by(2.25)K(x)=k0U∞2μ0x.
Boundary conditions of (2.24) are(2.26)f(0)=0,f′(0)=0,f′(∞)=1.
By substituting u and v 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:(2.27)τxy∣y=0=μ0U∞U∞υxf′′(0).


### 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:(2.28)f(0)=0,f′(0)=1,f′(∞)=0.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:(2.29)τxy∣y=0=(1-K)μ0U∞U∞υxf′′(0).


## 2.1. Stagnation-Point Flow

Cauchy equations are employed for obtaining boundary layer equation of two-dimensional stagnation-point flow. Steady Cauchy equations are(2.1)Continuity:∂u∂x+∂v∂y=0,x-momentum:ρ(u∂u∂x+v∂u∂y)=-∂p∂x+∂τxx∂x+∂τxy∂y,y-momentum:ρ(u∂v∂x+v∂v∂y)=-∂p∂y+∂τxy∂x+∂τyy∂y,
where p is pressure, ρ is density of fluid, and u and v are velocity components in x and y directions, respectively. Also,  τxx, τxy, τyx, and τyy are components of stress matrix. These variables can be obtained for Walters-B viscoelastic fluid by following equation:(2.2)τ=2μ0d-2k0d∇
which is the constitutive equation of Walters’ B fluid and μ0 is the dynamic viscosity, k0 is the short-memory coefficient, and d is the deformation rate tensor. Also d∇ is upper-convected derivative of a tensor and can be obtained from this equation:(2.3)d∇=∂d∂t+(v⃗⋅∇)d-d⋅∇VT-∇V⋅d.
So, the stress components in steady conditions are(2.4)τij=2μ0dij-2k0(DdijDt-∂vj∂xmdim-∂vi∂xmdmj),
where dij can be obtained as follows:(2.5)dij=12[∂Vi∂xj+∂Vj∂xi].
On the basis of the above relation, the stress components can be written as follows:(2.6)τxx=2μ0∂u∂x-2k0[u∂2u∂x2+v∂2u∂x∂y-2{(∂u∂x)2+12∂u∂y(∂u∂y+∂v∂x)}],τyx=τxy=μ0(∂u∂y+∂v∂x)-2k0[12u(∂2u∂x∂y+∂2v∂x2)+12v(∂2u∂y2+∂2v∂x∂y)-{∂u∂x∂v∂x+∂u∂y∂v∂y}],τyy=2μ0∂v∂y-2k0[∂2v∂t∂x+u∂2v∂x∂y+v∂2v∂y2-2{12∂v∂x(∂u∂y+∂v∂x)+(∂v∂y)2}].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 becomex-momentum:(2.7)u∂u∂x+v∂u∂y=-1ρ∂p∂x+μ0ρ(∂2u∂x2+∂2u∂y2)-k0ρ[u(∂3u∂x3+∂3u∂x∂y2)+v(∂3u∂x2∂y+∂3u∂y3)-∂u∂x(3∂2u∂x2-∂2u∂y2)-∂u∂y(∂2u∂x∂y+∂2v∂x2)-2∂v∂x∂2u∂x∂y],y-momentum:(2.8)u∂v∂x+v∂v∂y=-1ρ∂p∂y+μ0ρ(∂2v∂x2+∂2v∂y2)-k0ρ[u(∂3v∂x3+∂3v∂x∂y2)+v(∂3v∂x2∂y+∂3v∂y3)-∂v∂y(3∂2v∂y2-∂2v∂x2)-∂v∂x(∂2v∂x∂y+∂2u∂y2)-2∂u∂y∂2v∂x∂y].
The x- and y-momentum equations can be simplified as follows after employing boundary layer approximation:(2.9)u∂u∂x+v∂u∂y=-1ρdpdx+μ0ρ(∂2u∂y2)-k0ρ[u∂3u∂x∂y2+v∂3u∂y3+∂u∂x∂2u∂y2-∂u∂y∂2u∂x∂y]∂p∂y=0.On the basis of viscous flow theory, the change in pressure across boundary layer isO(δ2). 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:(2.10)-1ρ∂p∂x=U∞dU∞dx,
where U∞ is the mainstream velocity. So, boundary layer equation for this fluid flow is obtained as follows:(2.11)u∂u∂x+v∂u∂y=U∞dU∞dx+υ(∂2u∂y2)-k0ρ[u∂3u∂x∂y2+v∂3u∂y3+∂u∂x∂2u∂y2-∂u∂y∂2u∂x∂y].
Because of the similarity of the velocity profiles in different cross-sections, it is concluded that the only possible form of U∞ under these conditions for obtaining the solution of (2.11) is(2.12)U∞=ε̇x,
where ε̇ is a constant parameter [37]. This procedure has been described in [37] completely.Thex-momentum equation can be nondimensionalized using stream function Ψ. So, the velocity components can be written as(2.13)u=∂Ψ∂y,v=-∂Ψ∂x.
The stream function can be made dimensionless as follows:(2.14)Ψ(x,y)=υε̇xf(η),
where f(η)f is the dimensionless stream function, υ is the kinematic viscosity, and η is the similarity variable, defined by(2.15)η(x,y)=(ε̇υ)y.
By substituting (2.11)–(2.15) into (2.9), dimensionless equation can be obtained:(2.16)f′′′+ff′′+1-f′2+K(ffiv-2f′f′′′+f′′2)=0,
where K is the local Weissenberg number defined by(2.17)K=k0ε̇μ0.
On the basis of no slip condition, no permeability, and u(y→∞)=ε̇x, three boundary conditions of (2.16) are(2.18)f(0)=0,f′(0)=0,f′(∞)=1.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.19)τxy=μ0∂u∂y-k0(u∂2u∂x∂y+v∂2u∂y2+2∂u∂x∂u∂y).
By substituting u and v 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:(2.20)τxy∣y=0=ε̇3/2xμ0υf′′(0).


## 2.2. Blasius Flow

Equation of motion for Blasius flow can be obtained using (2.1) that leads to (2.9).Thex-momentum equation of (2.9) can be nondimensionalized using stream function Ψ. So, the velocity components can be written as(2.21)u=∂Ψ∂y,v=-∂Ψ∂x.
The stream function can be made dimensionless as follows:(2.22)Ψ(x,y)=υxU∞f(η),
where f is the dimensionless stream function, υ is the kinematic viscosity, and η is the similarity variable, defined by(2.23)η(x,y)=U∞υxy.
So, nondimensional equation of motion for Blasius flow can be obtained on the basis of dp/dx=0 as follows:(2.24)f′′′+12ff′′+K(ffiv+2f′f′′′-f′′2)=0,
where K is the local Weissenberg number defined by(2.25)K(x)=k0U∞2μ0x.
Boundary conditions of (2.24) are(2.26)f(0)=0,f′(0)=0,f′(∞)=1.
By substituting u and v 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:(2.27)τxy∣y=0=μ0U∞U∞υxf′′(0).


## 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:(2.28)f(0)=0,f′(0)=1,f′(∞)=0.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:(2.29)τxy∣y=0=(1-K)μ0U∞U∞υxf′′(0).


## 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)f′=p,p′=Q
should be used to rewritte (2.16) as follows:(3.2)kfQ′′+(1-2kp)Q′+fQ+1-p2+kQ2=0,
and boundary conditions (2.18) become(3.3)f(0)=0,p(0)=0,p(∞)=1.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 Q and Q′.The following equations should be used for approximating the first- and second-order derivatives:(3.4)Dyi=1h(∇-12∇2-16∇3-⋯)yi+1,D2yi=1h2(∇2-12∇4-16∇5-⋯)yi+1,
where D 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:(3.5)1h(14fi+1+56fi-32fi-1+12fi-2-112fi-3)=pi,1h(14pi+1+56pi-32pi-1+12pi-2-112pi-3)=Qi,kfih(56Qi+1+54Qi-13Qi-1+76Qi-2-12Qi-3+112Qi-4)+1-2kpih(14Qi+1+56Qi-32Qi-1+12Qi-2-112Qi-3)+fiQi+1-pi2+kQi2=0.
Equation (3.5) can be explicitly solved for f,p,andQ at the (i+1)th point in a supposing mesh generation. Using the subscript “p” for the predicted values, we can write(3.6)fi+1(p)=13(-10fi+18fi-1-6fi-2+fi-3+12hpi),pi+1(p)=13(-10pi+18pi-1-6pi-2+pi-3+12hQi),Qi+1(p)=[10kfi+3h(1-2kpi)]-1[12h2(fiQi+1-pi2+kQi2)kfi(15Qi+4Qi-1-14Qi-2+6Qi-3-Qi-4)+h(1-2kpi)(-10Qi+18Qi-1-6Qi-2+Qi-3)-12h2(fiQi+1-pi2+kQi2)].
These steps were about the predicting the values of fi+1,pi+1,andQi+1. So, the following operators are necessary to obtain the correct values of these variables:(3.7)Dyi=1h(∇+12∇2+13∇3+14∇4⋯)yi,D2yi=1h2(∇2+∇3+1112∇4+56∇5+⋯)yi.
On the basis of (3.7), (3.1) and (3.2) become(3.8)fi+1=125(48fi-36fi-1+16fi-2-3fi-3+12hpi+1(p)),pi+1=125(48pi-36pi-1+16pi-2-3pi-3+12hQi+1(p)),Qi+1=[45kfi+1(p)+25h(1-2kpi+1(p))]-1[12h2(fi+1(p)Qi+1(p)+1-pi+1(p)2+kQi+1(p)2)kfi+1(p)(154Qi-214Qi-1+156Qi-2-61Qi-3+10Qi-4)+h(1-2kpi+1(p))(48Q-36Qi-1+16Q-3Qi-3)-12h2(fi+1(p)Qi+1(p)+1-pi+1(p)2+kQi+1(p)2)].Therefore, we have a fourth-order predictor-corrector method for obtaining the values off,f′ and f′′ at a particular point based with respect to their values at the previous five mesh points. The values of f and f′are given at η=0. Since f′ 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 f′′ so that the boundary condition at infinity is satisfied.Assume that the amount of the initial guess off′′ is(3.9)f′′(0)=s.The amounts of inf,p,andQ 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 η=±hη=±2h. Taylor series for expansion of of  f around η=0 results in the following equation:(3.10)f(η)=∑i=05ciηi,
where ci are constant and should be obtained using the boundary values of f,p,andQ at η=0. Henceforth, some of these constants are(3.11)c0=0,c1=0,c2=12s,c3=-Ks2+16,c4=0,c5=s2120,
where c0,c1,andc2 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:(3.12)fi+1(p)=13(-10fi+18fi-1-6fi-2+fi-3+12hpi),pi+1(p)=13(-10pi+18pi-1-6pi-2+pi-3+12hQi),Qi+1(p)=[10kfi+3h(1+2kpi)]-1[12h2(12fiQi-kQi2)kfi(15Qi+4Qi-1-14Qi-2+6Qi-3-Qi-4)+h(1+2kpi)(-10Qi+18Qi-1-6Qi-2+Qi-3)-12h2(12fiQi-kQi2)].
Corrector:(3.13)fi+1=125(48fi-36fi-1+16fi-2-3fi-3+12hpi+1(p)),pi+1=125(48pi-36pi-1+16pi-2-3pi-3+12hQi+1(p)),Qi+1=[45kfi+1(p)+25h(1+2kpi+1(p))]-1[12h2(12fi+1(p)Qi+1(p)-kQi+1(p)2)kfi+1(p)(154Qi-214Qi-1+156Qi-2-61Qi-3+10Qi-4)+h(1+2kpi+1(p))(48Q-36Qi-1+16Q-3Qi-3)-12h2(12fi+1(p)Qi+1(p)-kQi+1(p)2)].
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(3.14)f(η)=∑i=05ciηi,
where(3.15)c0=0,c1=0,c2=12s,c3=-Ks26,c4=0,c5=-s2240,
where c0,c1,andc2 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 c0,c1, and so forth, are different from Blasius flow that using (2.19) and (2.21) are obtained as(3.16)c0=0,c1=1,c2=12s,c3=Ks26(2K+1),c4=-s48(3K+1),c5=6K2+4K+1240(2K+1)(3K+1)(4K+1).
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.


## 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)f′=p,p′=Q
should be used to rewritte (2.16) as follows:(3.2)kfQ′′+(1-2kp)Q′+fQ+1-p2+kQ2=0,
and boundary conditions (2.18) become(3.3)f(0)=0,p(0)=0,p(∞)=1.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 Q and Q′.The following equations should be used for approximating the first- and second-order derivatives:(3.4)Dyi=1h(∇-12∇2-16∇3-⋯)yi+1,D2yi=1h2(∇2-12∇4-16∇5-⋯)yi+1,
where D 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:(3.5)1h(14fi+1+56fi-32fi-1+12fi-2-112fi-3)=pi,1h(14pi+1+56pi-32pi-1+12pi-2-112pi-3)=Qi,kfih(56Qi+1+54Qi-13Qi-1+76Qi-2-12Qi-3+112Qi-4)+1-2kpih(14Qi+1+56Qi-32Qi-1+12Qi-2-112Qi-3)+fiQi+1-pi2+kQi2=0.
Equation (3.5) can be explicitly solved for f,p,andQ at the (i+1)th point in a supposing mesh generation. Using the subscript “p” for the predicted values, we can write(3.6)fi+1(p)=13(-10fi+18fi-1-6fi-2+fi-3+12hpi),pi+1(p)=13(-10pi+18pi-1-6pi-2+pi-3+12hQi),Qi+1(p)=[10kfi+3h(1-2kpi)]-1[12h2(fiQi+1-pi2+kQi2)kfi(15Qi+4Qi-1-14Qi-2+6Qi-3-Qi-4)+h(1-2kpi)(-10Qi+18Qi-1-6Qi-2+Qi-3)-12h2(fiQi+1-pi2+kQi2)].
These steps were about the predicting the values of fi+1,pi+1,andQi+1. So, the following operators are necessary to obtain the correct values of these variables:(3.7)Dyi=1h(∇+12∇2+13∇3+14∇4⋯)yi,D2yi=1h2(∇2+∇3+1112∇4+56∇5+⋯)yi.
On the basis of (3.7), (3.1) and (3.2) become(3.8)fi+1=125(48fi-36fi-1+16fi-2-3fi-3+12hpi+1(p)),pi+1=125(48pi-36pi-1+16pi-2-3pi-3+12hQi+1(p)),Qi+1=[45kfi+1(p)+25h(1-2kpi+1(p))]-1[12h2(fi+1(p)Qi+1(p)+1-pi+1(p)2+kQi+1(p)2)kfi+1(p)(154Qi-214Qi-1+156Qi-2-61Qi-3+10Qi-4)+h(1-2kpi+1(p))(48Q-36Qi-1+16Q-3Qi-3)-12h2(fi+1(p)Qi+1(p)+1-pi+1(p)2+kQi+1(p)2)].Therefore, we have a fourth-order predictor-corrector method for obtaining the values off,f′ and f′′ at a particular point based with respect to their values at the previous five mesh points. The values of f and f′are given at η=0. Since f′ 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 f′′ so that the boundary condition at infinity is satisfied.Assume that the amount of the initial guess off′′ is(3.9)f′′(0)=s.The amounts of inf,p,andQ 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 η=±hη=±2h. Taylor series for expansion of of  f around η=0 results in the following equation:(3.10)f(η)=∑i=05ciηi,
where ci are constant and should be obtained using the boundary values of f,p,andQ at η=0. Henceforth, some of these constants are(3.11)c0=0,c1=0,c2=12s,c3=-Ks2+16,c4=0,c5=s2120,
where c0,c1,andc2 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:(3.12)fi+1(p)=13(-10fi+18fi-1-6fi-2+fi-3+12hpi),pi+1(p)=13(-10pi+18pi-1-6pi-2+pi-3+12hQi),Qi+1(p)=[10kfi+3h(1+2kpi)]-1[12h2(12fiQi-kQi2)kfi(15Qi+4Qi-1-14Qi-2+6Qi-3-Qi-4)+h(1+2kpi)(-10Qi+18Qi-1-6Qi-2+Qi-3)-12h2(12fiQi-kQi2)].
Corrector:(3.13)fi+1=125(48fi-36fi-1+16fi-2-3fi-3+12hpi+1(p)),pi+1=125(48pi-36pi-1+16pi-2-3pi-3+12hQi+1(p)),Qi+1=[45kfi+1(p)+25h(1+2kpi+1(p))]-1[12h2(12fi+1(p)Qi+1(p)-kQi+1(p)2)kfi+1(p)(154Qi-214Qi-1+156Qi-2-61Qi-3+10Qi-4)+h(1+2kpi+1(p))(48Q-36Qi-1+16Q-3Qi-3)-12h2(12fi+1(p)Qi+1(p)-kQi+1(p)2)].
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(3.14)f(η)=∑i=05ciηi,
where(3.15)c0=0,c1=0,c2=12s,c3=-Ks26,c4=0,c5=-s2240,
where c0,c1,andc2 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 c0,c1, and so forth, are different from Blasius flow that using (2.19) and (2.21) are obtained as(3.16)c0=0,c1=1,c2=12s,c3=Ks26(2K+1),c4=-s48(3K+1),c5=6K2+4K+1240(2K+1)(3K+1)(4K+1).
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 Table1.Table 1
Variation off′′(0) in two-dimensional stagnation-point flow (ηmax=20).
K (Weissenberg number)h=0.1h=0.05h=0.02h=0.01Serth [38]01.2013888261.213497571.213566436—1.2325870.051.2946180791.294642261.2946466061.2946461.2946460.11.3695009921.369534491.3695410151.3695391.3695380.21.5603145291.587324931.5873281251.5873281.5873320.31.0888995321.138504452.1108215332.110818—On the basis of Table1, it is obvious that step size (h) must be smaller than 0.05 in this problem to obtain the reliable values of  f′′(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 K≤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 f′′(0). So, it is obvious that, for larger values of  f′′(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 f′′(0) are based on step sizes smaller than 0.05.The boundary layers in stagnation-point flow for different Weissenberg numbers are compared in Figure1. 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 K≥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:k=0, dashed line: k=0.05, long-dashed line: k=0.1).


### 4.2. Blasius Flow

Quantities off′′(0) for different Weissenberg numbers are obtained, for h=0.02, as shown in Table 2. It is shown that increasing the amount of K causes to reduce the amount of  f′′(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 K is opted 20, that is, similar to the solution of stagnation-point flow.Table 2
Variation off′′(0) in Blasius flow (ηmax=20,h=0.02).
K (Weissenberg number)f′′(0)00.33205734250.10.29707534790.20.26832115170.30.24440628050.40.22446540830.50.20688659660.60.19237102500.70.18043045040.80.16982551570.90.159889488210.1505795288Figure3 shows the growth of boundary layers in Blasius flow for different Weissenberg numbers. As K increases, the boundary layer thickness decreases. It should be noted that a nonphysical phenomenon occurred for larger values of Weissenberg number,  K≥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 off′′(0) in two-dimensional Sakiadis flow is clarified in Table 3. So, it is obvious that values of f′′(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 off′′(0) in Sakiadis flow (ηmax=20,h=0.02).
K (Weissenberg number)f′′(0)0−0.443491970.2−0.456582030.4−0.475299450.6−0.505254210.7−0.530166503The most important result of this study is shown in Figure4. 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 off′′(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 off′′(0) for 2nd and 4th order of approximation of boundary conditions in stagnation-point flow (h=0.01).
K (Weissenberg number)f′′(0)2nd-order boundary4th-order boundary01.23213891.23258760.051.29412571.29464670.11.36891911.36953890.21.58633051.58732760.32.10780142.1108185Table 5
Comparison off′′(0) for 2nd and 4th order of approximation of boundary conditions in Blasius flow (h=0.02).
K (Weissenberg number)f′′(0)2nd-order boundary4th-order boundary00.33191497800.33205734250.10.29696151730.29707534790.20.26822814940.26832115170.30.24432891840.24440628050.40.22439956660.22446540830.50.20682987210.20688659660.60.19232215880.19237102500.70.18038703910.18043045040.80.16978637690.16982551570.90.15985383980.159889488210.15054714200.1505795288Table 6
Comparison off′′(0) for 2nd and 4th order of approximation of boundary conditions in Sakiadis flow (h=0.02).
K (Weissenberg number)f′′(0)2nd-order boundary4th-order boundary0−0.44364105−0.443491970.2−0.45671180−0.456582030.4−0.47542854−0.475299450.6−0.50540725−0.505254210.7−0.53035644−0.530166503It is obvious that the amounts off′′(0) for all three problems are more accurate by the 4th order of approximation in boundary conditions.


## 4.1. Stagnation-Point Flow

The results of the solution of stagnation-point flow with the proposed fourth-order method are presented in Table1.Table 1
Variation off′′(0) in two-dimensional stagnation-point flow (ηmax=20).
K (Weissenberg number)h=0.1h=0.05h=0.02h=0.01Serth [38]01.2013888261.213497571.213566436—1.2325870.051.2946180791.294642261.2946466061.2946461.2946460.11.3695009921.369534491.3695410151.3695391.3695380.21.5603145291.587324931.5873281251.5873281.5873320.31.0888995321.138504452.1108215332.110818—On the basis of Table1, it is obvious that step size (h) must be smaller than 0.05 in this problem to obtain the reliable values of  f′′(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 K≤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 f′′(0). So, it is obvious that, for larger values of  f′′(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 f′′(0) are based on step sizes smaller than 0.05.The boundary layers in stagnation-point flow for different Weissenberg numbers are compared in Figure1. 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 K≥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:k=0, dashed line: k=0.05, long-dashed line: k=0.1).


## 4.2. Blasius Flow

Quantities off′′(0) for different Weissenberg numbers are obtained, for h=0.02, as shown in Table 2. It is shown that increasing the amount of K causes to reduce the amount of  f′′(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 K is opted 20, that is, similar to the solution of stagnation-point flow.Table 2
Variation off′′(0) in Blasius flow (ηmax=20,h=0.02).
K (Weissenberg number)f′′(0)00.33205734250.10.29707534790.20.26832115170.30.24440628050.40.22446540830.50.20688659660.60.19237102500.70.18043045040.80.16982551570.90.159889488210.1505795288Figure3 shows the growth of boundary layers in Blasius flow for different Weissenberg numbers. As K increases, the boundary layer thickness decreases. It should be noted that a nonphysical phenomenon occurred for larger values of Weissenberg number,  K≥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 off′′(0) in two-dimensional Sakiadis flow is clarified in Table 3. So, it is obvious that values of f′′(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 off′′(0) in Sakiadis flow (ηmax=20,h=0.02).
K (Weissenberg number)f′′(0)0−0.443491970.2−0.456582030.4−0.475299450.6−0.505254210.7−0.530166503The most important result of this study is shown in Figure4. 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 off′′(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 off′′(0) for 2nd and 4th order of approximation of boundary conditions in stagnation-point flow (h=0.01).
K (Weissenberg number)f′′(0)2nd-order boundary4th-order boundary01.23213891.23258760.051.29412571.29464670.11.36891911.36953890.21.58633051.58732760.32.10780142.1108185Table 5
Comparison off′′(0) for 2nd and 4th order of approximation of boundary conditions in Blasius flow (h=0.02).
K (Weissenberg number)f′′(0)2nd-order boundary4th-order boundary00.33191497800.33205734250.10.29696151730.29707534790.20.26822814940.26832115170.30.24432891840.24440628050.40.22439956660.22446540830.50.20682987210.20688659660.60.19232215880.19237102500.70.18038703910.18043045040.80.16978637690.16982551570.90.15985383980.159889488210.15054714200.1505795288Table 6
Comparison off′′(0) for 2nd and 4th order of approximation of boundary conditions in Sakiadis flow (h=0.02).
K (Weissenberg number)f′′(0)2nd-order boundary4th-order boundary0−0.44364105−0.443491970.2−0.45671180−0.456582030.4−0.47542854−0.475299450.6−0.50540725−0.505254210.7−0.53035644−0.530166503It is obvious that the amounts off′′(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.



---

*Source: 861508-2012-01-26.xml*