# Spectral Local Linearisation Approach for Natural Convection Boundary Layer Flow

Authors: S. S. Motsa; Z. G. Makukula; S. Shateyi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The present work introduces a spectral local linearisation method (SLLM) to solve a natural convection boundary layer flow problem with domain transformation. It is customary to find solutions of semi-infinite interval problems by first truncating the interval and subsequently applying a suitable numerical method. However, this gives rise to increased error terms in the numerical solution. Carrying out a transformation of the semi-infinite interval problems into singular problems posed on a finite interval can avoid the domain truncation error and enables the efficient application of collocation methods. The SLLM is based on linearising and decoupling nonlinear systems of equations into a sequence or subsystems of differential equations which are then solved using spectral collocation methods. A comparative study between the SLLM and existing results in the literature was carried out to validate the results. The method has shown to be a promising efficient tool for nonlinear boundary value problems as it gives converging results after very few iterations.


---

## Body

## 1. Introduction

Boundary layer problems involve a rapid change in the value of a physical variable over a limited region of time or space and therefore are governed by highly nonlinear equations. They are very common in physics and have numerous applications in extrusion, melt spinning, food processing, mechanical forming processes, and many other fields.In recent years, a host of different methods have been developed with the aim of finding the most efficient method for solving the often complex nonlinear equation systems governing the fluid flows. Boundary layer problems are very difficult to solve because besides having very thin regions where there is rapid change of the fluid properties they are defined on unbounded domains. Analytical methods are being used to solve boundary layer problems but converge very slowly for some boundary layer problems, particularly those with very large parameters. Traditional numerical techniques such as finite differences, finite element methods, Keller-Box method [1], Runge-Kutta with shooting methods [2–4], and spectral collocation methods [5–7] are still being used in solving nonlinear BVPs arising in fluid mechanics and other engineering applications. However, numerical approaches also encounter problems in resolving the solution of the governing equations in the very thin regions and in cases where there are singularities or multiple solutions hence the need for the search for ways of improving numerical methods.In this paper we present a new numerical technique based on linearising and decoupling nonlinear systems of equations into a sequence of or subsystems of differential equations which are then solved using spectral collocation methods. The spectral local linearisation method (SLLM) is used here to solve the problem of natural convection boundary layer flow that arises in glass-fibre production process. Free convection boundary layer problems arise in many applications including buoyancy-induced motions in the atmosphere, oceans, and quasi-solid bodies such as earth. The problem under study was earlier considered in the works of Kuiken [8, 9], Liao [10], and Motsa et al. [11]. In this work, however, the problem is first transformed from the semi-infinite domain [0,∞) into the finite interval [0,1] to avoid using a scaling parameter to invoke the boundary condition at infinity. The skin friction and Nusselt number are generated by the SLLM and validated by comparing them with those of Kuiken [8].


## 2. Problem Statement

Consider the laminar free-convection flow of an incompressible viscous fluid about a flat plate parallel to the direction of the generating body force. Figure1 presents the physical model of the flow. The basic governing boundary layer equations, following Kuiken [9], Ghotbi et al. [12], and Ali et al. [13], are as follows:
(1)∂u∂x+∂u∂y=0,(2)u∂u∂x+v∂u∂y=ν∂2u∂y2+gβ(T-T∞),(3)u∂T∂x+v∂T∂y=α∂2T∂y2,
with boundary conditions
(4)aty=0;u=v=0,T=T0,asy⟶∞;u⟶0,T⟶T∞,
where u and v are the velocity components along the x and y-axes, respectively, ν is the kinematic viscosity of the fluid, g is the acceleration due to gravity, β is the coefficient of thermal expansion, T is the fluid temperature, T∞ is the temperature of the surrounding fluid, and α is the thermal diffusivity.Figure 1
Free convection flow over a vertical plate.Applying the transformations below following [9, 12, 13]
(5)ψ=[gβ(Tw-T∞)ν2x02]1/4f(η),T=T∞+(Tw-T∞)[xx0-x]3θ(η),η=[gβ(Tw-T∞)x03ν2]y(x0-x),
where ψ is the stream function defined as
(6)u=∂ψ∂y,v=-∂ψ∂x,f(η) and θ(η) are similarity functions, (2) and (3) become
(7)f′′′(η)-f′2(η)+θ(η)=0,θ′′(η)-3Prf′(η)θ(η)=0,
along with the boundary conditions
(8)f(0)=f′(0)=0,θ(0)=1,f′(∞)=θ(∞)=0.Prime denotes differentiation with respect to the similarity variableη and Pr is the Prandtl number. It was reported in [9] that both f(η) and θ(η) decay algebraically as η→+∞. Equations (7) can be simplified by using the following transformations [10]:
(9)ξ=1+λη,F(ξ)=f′(η),G(ξ)=θ(η),
where λ is a spatial scale parameter. Substituting  (9) in (7) and (8) results in
(10)λ2F′′(ξ)-F2(ξ)+G(ξ)=0,λ2G′′(ξ)-3PrF(ξ)G(ξ)=0,
with boundary conditions
(11)F(1)=0,G(1)=1,F(∞)=G(∞)=0,
where the prime now denotes differentiation with respect to ξ.


## 3. Domain Transformation

In this section we transform the problem from the semi-infinite domain into the finite domain[0,1]. The transformation procedure proposed by Hammerling et al. [14] transforms a general BVP of the form
(12)f(x,u,u′,u′′,…,u(n))=0,x∈[1,∞),
into the domain [0,1] by change of the independent variable. A transformation of the form x→1/x transforms the semi-infinite interval into [0,1]. For the general case, the function and its derivatives are then expressed as
(13)u(1x)=U(x)u′(1x)=-x2U′(x)u′′(1x)=-x2(-x2U′(x))′u′′′(1x)=-x2(-x2(-x2U′(x))′)′⋮u(n)(1x)=-x2[-x2(-x2(…(-x2U′(x))′)′)′]′.︸(n-1)-productsThe boundary conditions are transformed in a similar manner with an exception of a problem with Neumann boundary conditions at the free end. In that case the problem has to be reduced to have admissible boundary conditions [14]. Previous studies carried out by Budd et al. [15] and Kitzhofer et al. [16] showed that transforming a problem into a finite domain before it is being solved facilitates the efficient use of collocation methods, and classical convergence orders can be observed.Using (13), (10) together with boundary conditions (11) is transformed into
(14)λ2(ξ4F1′′+2ξ3F1′)-F12+G1=0,λ2(ξ4G1′′+2ξ3G1′)-3PrF1G1=0,(15)F1(1)=0,G1(1)=1,F1(0)=0,G1(0)=0.


## 4. Spectral Local Linearisation Method (SLLM)

The transformed system of (14) and (15) is solved in this section using the spectral local linearisation method. The SLLM is based on developing a decoupled iterative scheme that is then chronologically solved using spectral methods. Local linearisation is applied to the equations before generating the iterative scheme in a manner similar to the Gauss-Seidel approach of decoupling linear algebraic systems. Consider a system of n nonlinear ordinary differential equations in n unknown functions Fi(η), i=1,2,…,n written as a sum of its linear and nonlinear terms as follows:
(16)Li[F1,F2,…,Fn]+Ni[F1,F2,…,Fn]=0,i=1,2,…,n.Let the subscriptsr and r+1 represent the previous and current iteration, respectively. Local linearisation of Ni at the previous iteration about Fi is then carried out as follows:
(17)Ni[F1,F2,…,Fn]=Ni[F1,r,F2,r,…,Fn,r]+∂Ni∂Fi[F1,r,F2,r,…,Fn,r](Fi-Fi,r),
so that, at the current iteration, (16) becomes
(18)Li[F1,r+1,F2,r+1,…,Fn,r+1]+∂Ni∂Fi[⋯]Fi,r+1=∂Ni∂Fi[⋯]Fi,r-Ni[F1,r,F2,r,…,Fn,r],
where [⋯] represents [F1,r,F2,r,…,Fn,r]. Using the Gauss-Seidel approach, whereby the updated solutions Fj(j<i), obtained at previous equations, are used to find a solution Fi at the current iteration level Fi,r+1, the local linearisation iteration scheme can be generated as follows:
(19)L1[F1,r+1,F2,r,…,Fn,r]+∂N1∂F1[⋯]F1,r+1=∂N1∂F1[⋯]F1,r-N1[F1,r,…,Fn,r],L1[F1,r+1,F2,r+1,F3,r,…,Fn,r]+∂N2∂F2[⋯]F2,r+1=∂N2∂F2[⋯]F2,r-N2[F1,r+1,F2,r…,Fn,r],⋮Ln[F1,r+1,F2,r+1,…,Fn,r+1]+∂Nn∂Fn[⋯]Fn,r+1=∂Nn∂Fn[⋯]Fn,r-Nn[F1,r+1,…,Fn-1,r+1,Fn,r],
where [⋯] at the ith iteration represents [F1,r+1,F2,r+1,…,Fi-1,r+1,Fi,r,…,Fn,r]. Hence, starting from an initial approximation F1,0,F2,0,…,Fn,0, the iterative scheme (19) is solved until convergence is reached for all unknowns.Following the guidelines above, we generate the iterative scheme for (14) with the corresponding boundary conditions as follows:
(20)λ2(ξ4F1,r+1′′+2ξ3F1,r+1′)-2F1,rF1,r+1=-F1,r2-G1,r,F1,r+1(0)=0,F1,r+1(1)=0,λ2(ξ4G1,r+1′′+2ξ3G1,r+1′)-3PrF1,r+1G1,r+1=0,G1,r+1(0)=0,G1,r+1(1)=1.Since the coefficient parameters and the right hand side of (20) are known (from previous iterations), the equation system (20) can easily be solved using any numerical method such as finite differences, finite elements, Runge-Kutta-based shooting methods, or collocation methods. In this work, the equations 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
(21)Xj=cosπjN,j=0,1,…,N,
where N is the number of collocation points used (see, e.g., [17, 18]). In order to implement the method, the physical region [0,1] is transformed into the region [-1,1] using the following mapping:
(22)ξ=X+12,-1≤X≤1.
The unknown functions F1,r+1 and G1,r+1 are approximated at the collocation points by
(23)F1,r+1(X)≈∑k=0NF1,r+1(Xk)Tk(Xj),G1,r+1(X)≈∑k=0NG1,r+1(Xk)Tk(Xj),00000000000000.00j=0,1,…,N,
where N is the number of collocation points and Tk is the kth Chebyshev polynomial defined as
(24)Tk(X)=cos[kcos-1(X)].The derivatives of the variables at the collocation points are represented as(25)d2F1,r+1dξ2=∑k=0NDkj2F1,r+1(Xk),d2G1,r+1dξ2=∑k=0NDkj2G1,r+1(Xk),0000000000.0j=0,1,…,N,
where D=2𝒟 with 𝒟 being the Chebyshev spectral differentiation matrix (see, e.g., [17, 18]). Substituting (22)–(25) in (20) we obtain the following system:
(26)A1F1,r+1=B1,F1,r+1(XN)=0,F1,r+1(X0)=0,A2G1,r+1=B2,G1,r+1(XN)=0,G1,r+1(X0)=1,
where
(27)A1=λ2diag(ξ4)D2+2λ2diag(ξ3)D-2diag(F1,r),B1=-F1,r2-G1,r,B2=λ2diag(ξ4)D2+2λ2diag(ξ3)D-3Prdiag(F1,r+1),B2=0.In the above equationsF1,r+1 and G1,r+1 are the values of functions F1,r+1 and G1,r+1, respectively, and 0 is a (N+1)×1 zero vector.Equations (26) are then solved chronologically starting from the initial approximate solutions
(28)F1,0(ξ)=γ(ξ2-ξ3),G1,0(ξ)=ξ4,
which are chosen to satisfy the boundary conditions (15). These initial guesses were suggested in [10] and are presented here after being transformed into the finite interval. The parameter γ appearing in (28) is the so-called auxiliary parameter which is used to control convergence of the solution.The convergence rate of the SLLM algorithm can be significantly improved by applying the successive overrelaxation (SOR) technique on (20). Under the SOR framework, a convergence controlling relaxation parameter ω is introduced and the SLLM scheme for finding, say Y, is modified to
(29)A1Yr+1=(1-ω)A1Yr+ωB1.The results in the following section show that, forω<1, applying the SOR method improves the efficiency and accuracy of the SLLM.


## 5. Results and Discussion

Numerical computations for the solution of the governing equations (14) and (15) were carried out using the spectral local linearization method. Values of the skin friction and the Nusselt number were generated at different Prandtl numbers. Results for the typical velocity and temperature profiles were also presented. In order to measure the accuracy and convergence of the SLLM, the present results were verified by comparing them with those obtained by Kuiken [8]. The plots were generated using 50 collocation points unless specified.As pointed out in Section4 that γ is a convergence controlling parameter, a proper selection of γ is necessary. To choose proper values of γ we use the norm of residual function of F1(ξ) (ResF) and that of G1(ξ) (ResG) where ResF and ResG are plotted as functions of γ:
(30)ResF=∥λ2(ξ4F1,r+1′′+2ξ3F1,r+1′)-F1,r+12+G1,r+1∥∞,ResG=∥λ2(ξ4G1,r+1′′+2ξ3G1,r+1′)-3PrF1,r+1G1,r+1∥∞.The valid region ofγ where the solution converges becomes the horizontal part of the residual curve. Numerical iterations are carried out until a representative solution residual falls below a set tolerance. In this work the tolerance level was set to be 10-10. Illustrative residual curves are shown in Figures 2 and 3. For example, in Figure 2, the valid values of γ that will give converging results for λ=2 are λ≥30 and about λ≥12 for λ=1.Figure 2
The residual ofF1(ξ) at different values of λ for Pr=0.1 when N=50.Figure 3
The residual ofG1(ξ) at different values of λ for Pr=0.1 when N=50.Figures4 and 5 show the velocity and temperature profiles, respectively, when the Prandtl number Pr is varied. The dimensionless velocity profiles are observed to decrease with increase in Pr in Figure 4. Because of the nature of the boundary conditions (slip boundary condition not considered), the velocity profile is seen to increase to reach a particular maximum point and then reduces to zero as η→∞. Similarly, as Pr increases, the temperature profiles also decrease as can be seen from Figure 5. For larger values of the Prandtl number, the thermal boundary layer thickness tends to be shorter because it is when the effective thermal diffusivity is low for a fixed kinematic viscosity. The results are consistent with results in the literature (see, e.g., [10, 12]).Figure 4
The effect of the Prandtl number on the velocity profile.Figure 5
The effect of the Prandtl number on the temperature profile.Tables1 and 2 present a comparison of the SLLM solutions of the skin friction and the Nusselt number for the flow at different values of the Prandtl number against those of Kuiken [8]. The effect of the Prandtl number on the skin friction coefficient and the Nusselt number is further shown in Figures 6 and 7. Convergence to the numerical solutions obtained by Kuiken [8] is achieved for up to eight decimal places. Values of the skin friction are accurate for all values of Pr while the Nusselt number is accurate for up to six decimal places in some cases. The skin friction decreased with increase in the Prandtl numbers while the heat transfer rates increased with increased Prandtl numbers. The tables also reveal the advantage of relaxation on the method. Convergence rates are improved remarkably for each case. Worth noting also is the advantage of transforming the domain into the finite interval. There was no need to search for a suitable scaling parameter to invoke the boundary conditions at infinity while a minimal number of collocation points were used to generate the results.Table 1
Comparison of SLLM solutions forf′′(0) against those of Kuiken [8] for different values of Pr when N=50.
Pr
γ
λ
Basic SLLM
ω
SLLM with SOR (ω)
Reference [8]
Iter.
f
′
′
(
0
)
Iter.
f
′
′
(
0
)
0.001
40
0.1
39
1.12313813
0.85
19
1.12313813
1.12313813
0.01
50
0.1
43
1.06338086
0.90
19
1.06338086
1.06338086
0.1
30
0.1
37
0.92408304
0.85
19
0.92408304
0.92408304
1
10
0.1
25
0.69321163
0.90
13
0.69321163
0.69321163
10
5
1
36
0.44711652
0.85
13
0.44711652
0.44711652
100
5
1
46
0.26452354
0.85
14
0.26452354
0.26452354
1000
5
1
49
0.15129020
0.85
15
0.15129020
0.15129020
10000
5
1
50
0.08554085
0.85
16
0.08554085
0.08554085Table 2
Comparison of the values of the SLLM solutions for-θ′(0) with those of Kuiken [8] for different values of Pr generated using N=50.
Pr
γ
λ
Basic SLLM
ω
SLLM with SOR (ω)
Reference [8]
Iter.
-
θ
′
(
0
)
Iter.
-
θ
′
(
0
)
0.001
30
0.1
36
0.04680746
0.85
24
0.04680746
0.04680746
0.01
30
0.1
41
0.13576074
0.85
19
0.13576074
0.13576074
0.1
30
0.5
37
0.35005967
0.90
18
0.35005967
0.35005967
1
5
1
27
0.76986120
0.90
13
0.76986119
0.76986120
10
5
1
36
1.49709921
0.85
13
1.49709921
1.49709921
100
5
1
46
2.74688550
0.85
14
2.74688549
2.74688550
1000
5
1
49
4.93494763
0.85
15
4.93494756
4.93494763
10000
5
1
50
8.80444927
0.85
16
8.80444960
8.80444927Figure 6
Variation of the skin friction withPr when λ=1 and γ=5.Figure 7
Variation of the heat transfer rate withPr when λ=1 and γ=5.


## 6. Conclusion

The problem of natural convection boundary layer flow with heat transfer has been considered in the current study. The coupled nonlinear system of differential equations was transformed into a finite domain and solved using the spectral local linearisation method (SLLM). The convergence rate of the SLLM algorithm was improved by applying the successive overrelaxation (SOR) technique. Solving the problem in a finite domain avoided the errors that come with approximating the scaling parameter for the boundary conditions at infinity leading to more accurate results. In this work the efficiency of domain transformation on spectral methods was revealed as accurate results were achieved while using relatively small numbers of collocation points. The flow velocity, temperature profiles, skin friction, and rate of surface heat transfer were successfully computed. Both the velocity and temperature of the flow were seen to be reduced by increased Prandtl numbers. The effect of the Prandtl number was seen to reduce the skin friction while increasing the rate of heat transfer. The SLLM results were also found to be consistent with results available in the literature. The method gave converging results for very large Prandtl numbers. While quite simple and easy to implement, the method has produced accurate results at quite remarkable convergence rates, thus displaying its potential as an improved numerical method for solving nonlinear boundary value problems.



---

*Source: 765013-2013-09-04.xml*