# On the Successive Linearisation Approach to the Flow of Reactive Third-Grade Liquid in a Channel with Isothermal Walls

Authors: S. S. Motsa; O. D. Makinde; S. Shateyi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The nonlinear differential equations modeling flow of a reactive third-grade
liquid between two parallel isothermal plates is investigated using a novel hybrid of numerical-analytical
scheme known as the successive linearization method (SLM). Numerical and graphical
results obtained show excellence in agreement with the earlier results reported in the literature. A
comparison with numerical results generated using the inbuilt MATLAB boundary value solverbvp4c demonstrates that the new SLM approach is a very efficient technique for tackling highly
nonlinear differential equations of the type discussed in this paper.


---

## Body

## 1. Introduction

The rheological properties of many fluids used in industrial and engineering processes do exhibit non-Newtonian behaviour [1, 2]. Meanwhile, the study of heat transfer plays an important role during the handling and processing of non-Newtonian fluids [3–5]. A complete thermodynamics analysis of the constitutive function for fluid of the differential type with the third-grade fluid being a special case has been performed by Fosdick and Rajagopal [6]. Similar studies with respect to non-Newtonian fluid are also reported by Makinde [7–9]. Moreover, the thermal boundary layer equations for non-Newtonian third-grade fluid constitute a nonlinear problem, and their solutions in space provide an insight into an inherently complex physical process in the system. In most cases, the nonlinear nature of the model equations precludes its exact solution. In recent time, several approximation techniques have been developed to tackle this problem [10–13], for example, the Adomian decomposition method, the variation iteration method, the improved finite differences method, the spectral method, and so forth. The ideas of developing new hybrids of numerical-analytical scheme to tackle nonlinear differential equations have experienced a revival. One such trend is the spectral homotopy analysis method that has recently been reported in [14, 15] which is a hybrid between the standard homotopy analysis method [16] and the Chebyshev spectral collocation method [17–19]. Other novel strategies involve using the Padé technique to improve the radius of convergence of the analytical methods of solution. Recent studies that make use of the Padé technique include the Hermite-Padé [7], the Homotopy-Padé [16, 20, 21], and the Hankel-Padé approaches.The purpose of the present work is to present a new method, called the successive linearisation method (SLM), of solving nonlinear boundary value problems. We demonstrate the applicability of SLM in tackling nonlinear differentiation equations modeling the flow of a reactive third-grade liquid between two parallel isothermal plates. The mathematical formulation of the problem is established in Section2. In Section 3 we introduce and apply some rudiments of SLM technique. Both numerical and graphical results are presented and discussed quantitatively with respect to various parameters embedded in the system in Section 4. A limited parametric study comparing numerical results generated using MATLAB’s bvp4c boundary value solver is compared with the SLM results, and good agreement is observed. Using the SLM approach, multiple solutions which were theoretically proved to exist for such problems in [7] are also generated.


## 2. Mathematical Formulation

Figure1 depicts the problem geometry. We consider the steady flow of an incompressible third-grade reactive fluid placed between two parallel isothermal plates. It is assumed that the flow is hydrodynamically and thermally fully developed under the action of a constant axial pressure gradient.Figure 1
Geometry of the problem.Following [1, 4, 5, 7–9], the dimensionless governing equations for the momentum and energy balance can be written as
(1)d2Wdy2+6γd2Wdy2(dWdy)2=-1,(2)d2θdy2+λ[eθ/(1+εθ)+m(dWdy)2×(1+2γ(dWdy)2)]=0,
where W is the dimensionless velocity component, y is the dimensionless normal coordinate, θ is the dimensionless temperature, and λ, ε, γ, and m represent the Frank-Kamenetskii parameter, activation energy parameter, the dimensionless non-Newtonian parameter, and the viscous heating parameter, respectively. The additional Arrhenius kinetics term in energy balance equation (2) is due to [3]. The appropriate boundary conditions in dimensionless form are given as follows: the surface of the channel is fixed, impermeable, and maintained at a given temperature:
(3)W(1)=0,θ(1)=0,aty=1,
and the symmetry condition along the centerline, that is,
(4)dWdy=dθdy=0,aty=0.We have employed the following nondimensional quantities in (1)–(4):
(5)θ=E(T-T0)RT02,y-=ya,λ=QEAa2C0e-E/RT0T02Rk,W=uUG,m=μG2U2eE/RT0QAa2C0,ε=RT0E,G=-a2μUdPdx,γ=β3U2G2a2μ,
where T is the absolute temperature, U is the fluid characteristic velocity, T0 is the plate temperature, k is the thermal conductivity of the material, Q is the heat of reaction, A is the rate constant, E is the activation energy, R is the universal gas constant, C0 is the initial concentration of the reactant species, a is the channel half width, β3 is the material coefficient, P is the modified pressure, and μ is the fluid dynamic viscosity coefficient. In the following sections, (1)–(4) are solved numerically using the new successive linearization technique.


## 3. Successive Linearisation Method (SLM)

In this section, we apply the proposed linearisation method of solution, hereinafter referred to as the successive linearisation method (SLM), to solve the governing equations (1) and (2). Before applying the SLM, we first note that the problem can be significantly simplified by finding the explicit analytical solution for the derivative dW/dy. First, (1) is rewritten as
(6)ddy[dWdy+2γ(dWdy)3]=-1.
The above equation is then integrated on both sides, and the symmetry boundary condition dW/dy(0)=0 is used to evaluate the resulting integrating constant to give
(7)2γ(dWdy)3+(dWdy)+y=0.
We note that (7) is a cubic equation which can have either one or three real solutions. If only positive values of γ are considered, (7) will have a unique real solution which can be computed using Maple and is given by
(8)dWdy=16γ[ψ(y)]1/3-1[ψ(y)]1/3,
where
(9)ψ(y)=(-54y+632+27γy2γ)γ2.The analytical result is very important because, when evaluated aty=1, it gives an explicit analytical expression for the skin friction coefficient (Cf). A close inspection of (8) indicates that if Cf has a critical point, then (γc,Cf)=(-3/2,-2/27). This critical point was reported as a bifurcation point in [7]. It is worth noting that dW/dy, and hence the solution W(y), is only valid when γ≥-2/27.Since the momentum equation (1) is decoupled from the energy equation (2), we solve for the velocity W(y) first then substitute the result in the energy equation to obtain θ(y). To solve W(y), we write (1) as
(10)d2Wdy2=f(y),
where f(y) is a known explicit function of y given by
(11)f(y)=-11+6γ(dW/dy)2,
with dW/dy given by (8). Equation (10) can easily be integrated using any numerical method. The SLM is based on the assumption that the unknown function θ(y) can be expanded as
(12)θ(y)=Ti(y)+∑n=0i-1‍θn(y),i=1,2,3,…,
where Ti are unknown functions and θn are successive approximations whose solutions are obtained recursively from solving the linear part of the equation that results from substituting (12) in the governing equations (2) using θ0(y) as an initial approximation. The linearisation technique is based on the assumption that Ti becomes increasingly smaller as i becomes larger, that is,
(13)limi→∞Ti=0.
Substituting (12) in (2) gives
(14)Ti′′+λexp[Ti+∑n=0i-1θn1+ε(Ti+∑n=0i-1θn)]=g(y)-∑n=0i-1θn′′,
where g(y) is a known function (from (8) and (9)) given by
(15)g(y)=-λm(dWdy)2[1+2γ(dWdy)2].
We choose θ0=0 as an initial approximation which is chosen to satisfy the boundary conditions. The subsequent solutions for θn  (n≥1) are obtained by successively solving the linearised form of (14) which are given as
(16)θi′′+ai-1θi=ri-1,
subject to the boundary conditions
(17)θi(-1)=θi(1)=0,
where
(18)ri-1=g(y)-[∑n=0i-1‍θn′′+λexp(∑n=0i-1θn1+ε∑n=0i-1θn)∑n=0i-1‍],ai-1=λ(1+ε∑n=0i-1θn)2exp[∑n=0i-1θn1+ε∑n=0i-1θn].Once each solution forθi  (i≥1) has been found from successively solving (11) and for each i, the approximate solutions for θ(y) are obtained as
(19)θ(y)≈∑n=0M‍θn(y),
where M is the order of SLM approximation.We remark that the coefficient parameterai-1 and the right-hand side ri-1 of (16) for i=1,2,3,… are known (from previous iterations). Thus, (16) can easily be solved using analytical means (whenever possible) or any numerical methods such as finite differences, finite elements, Runge-Kutta-based shooting methods, or collocation methods. In this work, (10) and (16) are solved using the Chebyshev spectral collocation method. This method is based on approximating the unknown functions by the Chebyshev interpolating polynomials in such a way that they are collocated at the Gauss-Lobatto points defined as
(20)yj=cosπjN,j=0,1,…,N,
where N is the number of collocation points used (see, e.g., [17, 19]). The derivatives are approximated at the collocation points by
(21)d2Wdy2=∑k=0N‍Dkj2W(yk),d2θidy2=∑k=0N‍Dkj2θi(yk),j=0,1,…,N,
where D is the Chebyshev spectral differentiation matrix (see, e.g., [17, 19]). Substituting (21) in (10) and (16) leads to matrix equations given by
(22)D2W=F,W(y0)=W(yN)=0,(23)ATi=Ri-1,θi(y0)=θi(yN)=0,
in which A is a (N+1)×(N+1) square matrix and T, W, F, and R are (N+1)×1 column vectors defined by
(24)A=D2+ai-1,(25)F=[f(y0),f(y1),…,f(yN-1),f(yN)]T,(26)Ti=[θi(y0),θi(y1),…,θi(yN-1),θi(yN)]T,(27)W=[W(y0),W(y1),…,W(yN-1),W(yN)]T,(28)Ri-1=[ri-1(y0),ri-1(y1),…,ri-1(yN-1),ri-1(yN)]T.
In the above definitions, ai-1 is a diagonal matrix of size (N+1)×(N+1) and the superscript T denotes transpose. After modifying the matrix system (22) and (23) to incorporate boundary conditions, the solutions for W(y) and θi(y) are obtained as
(29)W=(D2)-1F,Ti=A-1Si-1.


## 4. Results

In this section, we present the results showing the velocity distribution and temperature distribution, for different values of the governing parameters. To check the accuracy of the proposed successive linearisation method (SLM), comparison is made with numerical solutions obtained using the MATLAB routinebvp4c, which is an adaptive Lobatto quadrature scheme. All the SLM results presented in this work were generated using N=60 collocation points.Figure2 depicts the effect of non-Newtonian parameter (γ) on both the velocity and temperature profiles. Generally, both fluid velocity and temperature profiles attained their maximum values along the channel centerline and minimum at the walls satisfying the boundary conditions. Moreover, a gradual decrease in the magnitude of fluid velocity and temperature profiles is noticed with an increase in the value of γ. This can be attributed to the fact that as γ increases, the fluid viscosity increases leading to a decrease in the flow rate. In Figure 3, we observed that the fluid temperature generally increases with an increase in the value of the Frank-Kamenetskii parameter (λ) due to the Arrhenius kinetics. The effect of viscous dissipation parameter (m) on the fluid temperature is displayed in Figure 4. The internal heat generation due to viscous heating increases as the parameter value of m increases leading to a general increase in the fluid temperature. Meanwhile, the possibility of a lower and upper solution branches is also highlighted in Figures 3 and 4. This can be attributed to the nonlinear nature of the Arrhenius kinetics in the governing thermal boundary layer equation (2). It is noteworthy that the fluid temperature decreases with an increase in the activation energy parameter (ε) as illustrated in Figure 5. As ε increases, the fluid becomes less volatile, and its activation energy decreases.Velocity and temperature profiles whenγ=0.1,1,5 when m=1, λ=0.3, and ε=0.1. SLM results (circles), bvp4c (solid line).
(a)
(b)Temperature profiles whenλ=0.4,0.5,0.6 when m=1, γ=0.1, and ε=0.1. SLM results (circles), bvp4c (solid line).
(a)
(b)Upper branch and lower branch temperature profiles whenm=0,5,10,15,20 when λ=0.5, γ=1, and ε=0.1. SLM results (circles), bvp4c (solid line).
(a)
(b)Figure 5
Temperature profiles whenε=0,1,5,20 when λ=0.5, γ=0.1, and m=1. SLM results (circles), bvp4c (solid line).A slice of the bifurcation diagram forγ>0 in the (λ,-θ′(1)) plane is shown in Figure 6. It represents the variation of wall heat flux with the Frank-Kamenetskii parameter (λ). In particular, for every 0≤ε≤0.1 there is a critical value λc (a turning point) such that, for 0≤λ<λc there are two solutions. This result is in perfect agreement with the one reported in Makinde [7].Figure 6
Bifurcation diagram whenγ=0.1, ε=0.1, and m=1.In Table1, we show the computations illustrating the comparison between the SLM results at different orders and the bvp4c numerical results for wall heat flux Nu for various values of m, γ, ε, and λ. It can be seen from the table that the SLM results are in very good agreement with the bvp4c numerical results.Table 1
Comparison between the SLM results at different orders and thebvp4c numerical results for wall heat flux   Nu for various values of  m, γ, ε, and λ.
m
γ
ε
λ
1st order
2nd order
3rd order
bvp4c
2
1
0.1
0.1
0.14948880
0.14958883
0.14958883
0.14958883
4
1
0.1
0.1
0.19550530
0.19565401
0.19565401
0.19565401
6
1
0.1
0.1
0.24152181
0.24172906
0.24172906
0.24172906
8
1
0.1
0.1
0.28753831
0.28781405
0.28781405
0.28781405
1
1
0.1
0.1
0.12648054
0.12655991
0.12655991
0.12655991
1
1
0.1
0.2
0.26159825
0.26233446
0.26233447
0.26233447
1
1
0.1
0.3
0.40652882
0.40944698
0.40944728
0.40944728
1
1
0.1
0.4
0.56267511
0.57092737
0.57093089
0.57093089
1
1
0.1
0.1
0.12648054
0.12655991
0.12655991
0.12655991
1
3
0.1
0.1
0.12168410
0.12175980
0.12175980
0.12175980
1
5
0.1
0.1
0.11956272
0.11963676
0.11963676
0.11963676
1
10
0.1
0.1
0.11691095
0.11698289
0.11698289
0.11698289
1
1
1
0.5
0.73172481
0.71360852
0.71359336
0.71359336
1
1
5
0.5
0.73172481
0.66332507
0.66312547
0.66312547
1
1
10
0.5
0.73172481
0.64501234
0.64484620
0.64484620
1
1
20
0.5
0.73172481
0.63173135
0.63164065
0.63164065


## 5. Conclusion

In this work, we employed a very powerful new linearisation technique, known as the successive linearisation method (SLM), to investigate the flow of reactive third-grade liquid in a channel with isothermal walls. The SLM results for the governing flow parameters were compared with results obtained using MATLAB’sbvp4c function, and excellent agreement was observed. From the results obtained in the study, the following was observed. (i)
An increase in the non-Newtonian parameter (γ) leads to a gradual decrease in the magnitude of fluid velocity and temperature profiles.(ii)
The fluid temperature generally increases with an increase in the value of the Frank-Kamenetskii parameter (λ).(iii)
The internal heat generation due to viscous heating increases as the parameter valuem increases leading to a general increase in the fluid temperature.It was also shown that the governing nonlinear equations admit multiple solutions. Using the SLM approach, lower and upper branch solutions were obtained and discussed.



---

*Source: 635392-2013-05-22.xml*