Open access peer-reviewed chapter

On the New Bivariate Local Linearisation Method for Solving Coupled Partial Differential Equations in Some Applications of Unsteady Fluid Flows with Heat and Mass Transfer

Written By

Sandile Motsa

Submitted: 28 October 2014 Reviewed: 23 February 2015 Published: 22 October 2015

DOI: 10.5772/60444

From the Edited Volume

Mass Transfer - Advancement in Process Modelling

Edited by Marek Solecki

Chapter metrics overview

1,759 Chapter Downloads

View Full Metrics


This work presents a new numerical approach for solving unsteady two-dimensional boundary layer flow with heat and mass transfer. The flow model is described in terms of a highly coupled and nonlinear system of partial differential equations that models the problem of unsteady mixed convection flow over a vertical cone due to impulsive motion. The proposed method of solution uses a local linearisation approach to decouple the original system of PDEs to form a sequence of equations that can be solved in a computationally efficient manner. Approximate functions defined in terms of bivariate Lagrange interpolation polynomials are used with spectral collocation to approximate the solutions of the decoupled linearised equations. To test the accuracy and to validate the results of the proposed method, numerical error analysis and convergence tests are conducted. The present results are also compared with results from published literature for some special cases. The proposed algorithm is shown to be very accurate, convergent and very effective in generating numerical results in a computationally efficient manner.


  • Unsteady boundary layer flow
  • heat and mass transfer
  • bivariate spectral collocation
  • impulsive flow

1. Introduction

In this investigation we revisit the problem of unsteady mixed convection flow over a vertical cone due to impulsive motion that was previously discussed in Singh and Roy [7]. In the flow model, the external stream is impulsively set into motion and the surface temperature is suddenly changed from its ambient temperature. The resulting constitutive equations describing the flow properties are a system of highly coupled non-linear partial differential equations (PDEs) which cannot be solved exactly. In [7], the solution of the PDEs system was approximated numerically using implicit finite differences after linearising using the quasi-linearisation (QLM) technique of Bellman and Kalaba [1]. The QLM simplifies the solution process by simply reducing the original PDEs to a linearised form which, nevertheless, is coupled. In view of the coupled nature of the QLM iteration scheme, matrices of very large dimensions are obtained when using any numerical method is used for discretization. For systems of coupled equations, the method proposed by [3] was used in [7] to decouple the governing equations. There is generally a large overhead in computing resources when performing computational tasks on large matrix systems. To avoid dealing with large matrix systems, relaxation methods can be used. However, the disadvantage of relaxation methods is that they convergence much slower than quasi-linearisation schemes which are known to converge quadratically.

A linearisation method, termed the local linearisation method (LLM), was recently introduced in [4] as an efficient method for solving coupled systems of non-linear ordinary differential equations that model boundary layer equations. This method was extended to a PDE system in [5] where the Chebyshev spectral collocation method [2, 8] was used for discretization in one independent variable and finite differences was used for discretization in the second independent variable of the PDE. In this study, we present a new approach that seeks to improve on the original LLM approach of [5] by eliminating the need for using finite differences for discretizing in one of the independent variables and, instead, apply spectral collocation independently in all independent variables of the PDE system. Finite differences require very fine grids (with a large number of grid points) to give more accurate solutions. In contrast, spectral methods are computationally less expensive, converge faster and are more accurate than finite difference methods particularly for problems with smooth solutions. The collocation method applied in this work uses bivariate Lagrange interpolation polynomials as basis functions. The proposed method converges fast and gives very accurate results which are obtained in a computationally efficient manner. The accuracy is established through residual and solution error analysis. Further validation of present results in established by comparing with existing results from literature.


2. Governing equations

The model under consideration is that of an unsteady mixed convection flow over a vertical cone that is impulsively set into motion to cause unsteadiness in the flow (see [6, 7]). It is assumed that buoyancy forces are present due to temperature and concentration variation of the fluid flow. The governing boundary layer momentum, energy and concentration equations were reduced to dimensionless form in [7] using transformations that were initially proposed in Williams and Rhyne [9] to give,


where Pr is the Prandtl number, Sc is the Schmidt number, f is the dimensionless stream function, θ and ϕ are the dimensionless temperature and concentration, N=λ/λ* is the ratio of the thermal (λ) and concentration λ* buoyancy parameter.

The boundary conditions are


Other quantities of physical interest include the local skin friction coefficient, Nusselt number and Sherwood numbers which are given, in dimensionless forms as,


3. Derivation of solution method

The derivation of the solution method is described for a general system of non-linear PDEs in this section. Consider a system of 3 coupled PDEs in f(η,ξ), θ(η,ξ) and ϕ(η,ξ)


where Ω1, Ω2 and Ω3 are non-linear operators that represent the non-linear PDEs (1), (2) and (3), respectively and F,T,H are defined as


Equation (8) can be simplified and decoupled by using the following algorithm called local linearisation method (LLM):

LLM Algorithm

  1. Solve for F in 1st equation assuming that T and H are known from previous iteration. Fr+1.

  2. With Fr+1 now known, solve for T in 2nd equation, assuming that H is known from previous iteration. Tr+1.

  3. Solve for H in 3rd equation Hr+1.

The LLM algorithm was recently reported in [4] as an efficient method for solving coupled system of non-linear ordinary differential equations that model boundary layer equations. The method was later extended to partial differential equations in [5]. Applying the algorithm in equations (8) gives


The above equation forms a system of three decoupled linear PDEs that are to be solved iteratively for f(η,ξ), θ(η,ξ) and ϕ(η,ξ), where is a vector of partial derivatives defined as


Applying the LLM iteration scheme (9) - (11) on the governing non-linear PDE system (1) - (3) gives


where the primes denote partial differentiation with respect to η and the coefficients ai,r,bi,r and ci,r (i=1,2,3,4) are known coefficients whose quantities are known from the previous iteration level r. The definition of the coefficients is given in the Appendix.

The boundary conditions are given by


Starting from given initial approximations, denoted by f0(η,ξ),θ0(η,ξ) and ϕ0(η,ξ), equations (12) – (14) are solved iteratively for r=1,2,, until approximate solutions that are consistent to within a certain tolerance level are obtained.

To solve the linearised equations (12) – (16), an approximate solution defined in terms of bivariate Lagrange interpolation polynomials in sought. For example, the approximate solution for f(η,ξ) takes the form


where the functions Lm(τ) and Lj(ζ) are the well-known characteristic Lagrange cardinal polynomials. The function (17) interpolates f(η,ξ) at the collocation points (known as Chebyshev-Gauss-Lobatto) defined by


The choice of collocation points (18) makes it possible for the Chebyshev spectral collocation method to be used as a solution procedure. The Chebyshev collocation method requires that the domain of the problem be transformed to [1,1]×[1,1]. Accordingly, linear transformations have been used to transform η[0,η] and ξ[0,Lt] to τ[1,1] and ζ[1,1], respectively. Here η is a finite value that is introduced to facilitate the application of the numerical method at infinity and Lt is the largest value of ξ.

Substituting equation (17) in equation (12)-(14) and making use of the derivatives formulas for Lagrange functions at Gauss-Lobatto points given in [2, 8], results in


where di,j (i,j=0,1,,Nt) are entries of the standard Chebyshev differentiation matrix d=[di,j] of size (Nt+1)×(Nt+1) (see, for example [2, 8]), D=(2/ηe)[Dr,s] (r,s=0,1,2,,Nx) with [Dr,s] being an (Nx+1)×(Nx+1) Chebyshev derivative matrix. The vectors and matrices are defined as


am,r(ξi) (m=1,2,3,4,5) is the diagonal matrix of the vector [am,r(τ0),am,r(τ1),,am,r(τNx)]T. The matrices b1,i,b2,i,c1,i,c2,i are also diagonal matrices corresponding to b1,r,b2,r,c1,r,c2,r evaluated at collocation points. Additionally, R1,i, R2,i and R3,i are (Nx+1)×1 vectors corresponding to R1,r, R2,r and R3,r, respectively, evaluated at the collocation points.

After imposing the boundary conditions for i=0,1,,Nt, equation (20) can be written as the following matrix equation:




where I is an (Nx+1)×(Nx+1) identity matrix.

Applying the same approach to (20) and (21) gives




Starting from a given initial guess, the approximate solutions for f(η,ξ), θ(η,ξ) and ϕ(η,ξ) are obtained by iteratively solving the matrix equations (22), (23) and (24), in turn, for r=0,1,2,. Convergence to the expected solution can be affirmed by considering the norm of the difference between successive iterations. If this quantity is less than a given tolerance level, the algorithm is assumed to have converged. The following error norms are defined for the difference between values of successive iterations,


4. Results and discussion

The local linearisation method (LLM) described by equations (12) - (14 was used to generate approximate numerical solutions for the governing systems of equations (1) - (3). The linearised equations were subsequently solved using the bivariate spectral collocation method as described in the previous section. The whole solution process is termed bivariate spectral local linearisation method (BSLLM). In this section, we present the results of the numerical computations for the various flow profiles. The BSLLM results are compared with published results from literature. In computing the numerical results presented in this paper, Nx=60 and Nt=15 collocation points in the η and ξ domain, respectively, were found to be sufficient to give accurate and consistent results. The minimum number of iterations required to give consistent results to within a tolerance level of 107 were determined using equation (25).

Figure 1.

Effect of λ and Pr on the velocity profile f(η,ξ)

Figure 2.

Effect of λ and Pr on the temperature profile θ(η,ξ)

The effects of buoyancy parameter (λ) and Prandtl number (Pr) on velocity and temperature profiles are shown in Figs. 1 and 2, respectively for when N=0.5, Sc=0.94 at ξ=0.5. It can be seen from the figures that the results from the present work match those reported in [7] exactly. We remark that in [7], the quasi-linearisation approach was used with implicit finite differences as a solution method. Details of the analysis of the effects of various governing parameters on the governing equations (1) – (3) can be found in [7] and have not been repeated here.

Figure 3.

Effect of N velocity profile f(η,ξ)

Figure 4.

Effect of N on the temperature profile θ(η,ξ)

Figure 5.

Effect of N concentration profile ϕ(η,ξ)

Figs. 3, 4 and 5 show the effect of buoyancy ratio N on the velocity, temperature and concentration profiles, respectively, when Pr=0.7, Sc=0.94 and λ=5 at ξ=0.5. The graphs are qualitatively similar to those reported in [7]. These figures are a qualitative validation of the numerical results generated using the proposed BSLLM.

Figure 6.

Effect of λ on the skin friction coefficient CfReL1/2

The effect of λ on the skin friction coefficient is shown in Fig. 6 for Pr=0.7, Sc=0.22 and N=0.5. Again, it is observed that the results match those of [7]. In particular, the oscillating trend in the skin friction coefficient for near the stagnation region can be observed.


5. Conclusion

The purpose of the current study was to develop a bivariate Lagrange polynomial based spectral collocation method for solving system of coupled non-linear partial differential equations. The applicability of the proposed method, termed bivariate spectral local linearisation method (BSLLM) was tested on the well-known problem of unsteady mixed convection flow over a vertical cone due to impulsive motion. The validity of the approximate numerical results was confirmed against known results from literature. The results of this study were found to be qualitatively congruent with those from published literature. The study revealed that the proposed method can be used as a viable approach for solving coupled partial differential equations with fluid mechanics applications.



Definition of coefficients



  1. 1. R.E. Bellman and R.E. Kalaba, Quasi-linearisation and nonlinear boundary-value problems, Elsevier, New York, 1965.
  2. 2. C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1988.
  3. 3. K. Inoyue, A. Tate, Finite differences version of quasilinearization applied to boundary layer equations, AIAA J. 12 (1974) 558-560.
  4. 4. S.S. Motsa, A new spectral local linearization method for nonlinear boundary layer flow problems, J. Applied Math. vol. 2013, Article ID 423628, 15 pages, 2013. doi:10.1155/2013/423628
  5. 5. S.S.Motsa, On the practical use of the spectral homotopy analysis method and local linearisation method for unsteady boundary-layer flows caused by an impulsively stretching plate, Numer. Algor. 66 (2014) 865-883
  6. 6. R. Seshadri, N. Sreeshylan, G. Nath, Unsteady mixed convection flow in the stagnation region of a heated vertical plate due to impulsive motion, Int. J. Heat Mass Transfer 45 (2002) 1345-1352.
  7. 7. P.R. Singh, S. Roy, Unsteady mixed convection flow over a vertical cone due to impulsive motion, International Journal of Heat and Mass Transfer 50 (2007) 949-959
  8. 8. L. N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000.
  9. 9. J.C. William, T.B. Rhyne, Boundary layer development on a wedge impulsively set into motion, SIAM J. Appl. Math. 38 (1980) 215-224.

Written By

Sandile Motsa

Submitted: 28 October 2014 Reviewed: 23 February 2015 Published: 22 October 2015