Preface

Chapter 8 **Disintegration Kinetics of Microbial Cells 169** Marek Solecki and Monika Solecka

Chapter 9 **Mass Transfer in Multiphase Systems 189** Badie I. Morsi and Omar M. Basha

da Silva

**VI** Contents

**Engineering 245**

Chapter 10 **Ion Exchange Fundamentals and New Challenges 219**

Chapter 11 **Application of Mass Transfer Models in Environmental**

Pen-Chi Chiang and Shu-Yuan Pan

Maria Angélica Simões Dornellas de Barros, Marcelino Luiz

Gimenes, Melissa Gurgel Adeodato Vieira and Meuris Gurgel Carlos

Mass transfer is the basis for the occurrence of many phenomena: physical, chemical, biolog‐ ical, as well as many others from related disciplines and areas of human activity. The trans‐ fer of different particle matter sets based on specific transport principles is described by analogous laws, congruent models, and similar relations. Research results obtained from mass transfer have a wide interdisciplinary application and, as a consequence, they have a great impact on the progress. The effects of intensive research carried out in selected direc‐ tions are often applied spectacularly to other areas. Therefore, the creation of a platform for exchange of observation, experience, and achievements in various fields in the area of mass transfer contributes to the development and dissemination of knowledge and improves processes and technical means to carry them out.

This book is aimed at presenting the most recent achievements of mass transfer theory and methods by stages of modelling different processes. Development of knowledge occurs cy‐ clically, generally reaching successive stages: observation, hypothesis, experiment, and for‐ mulation of theory. Searching for appropriate mathematical description must be preceded by a critical assessment of the present state and definition of demand for a new solution. Built models require validation. Such research stages from various academic disciplines are described in further chapters of this book. I hope that they will trigger further interactions between different fields of science. The book contains a wide variety of issues related to ad‐ vancements in modelling mass transfer. I would like to thank all the authors for their contri‐ bution to facilitate the understanding of complex problems and issues. The book is dedicated to those in the chemical, biotechnological, pharmaceutical, and nanotechnology industries as well as to students of natural sciences, technical, environmental and to employ‐ ees in companies which manufacture machines for the above-mentioned industries. This work can also be a useful source for researchers and engineers dealing with mass transfer and related issues.

> **Dr. Marek Solecki** Department of Process Equipment Lodz University of Technology, Lodz, Poland

**Chapter 1**

**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**

### Sandile Motsa

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/60444

#### **Abstract**

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.

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

**Keywords:** Unsteady boundary layer flow, heat and mass transfer, bivariate spec‐ tral collocation, impulsive flow

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **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 quasilinearisation (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 differ‐ ences for discretizing in one of the independent variables and, instead, apply spectral collo‐ cation 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 colloca‐ tion 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,

$$\begin{split} \frac{\hat{\sigma}^{3}f}{\hat{\sigma}\eta^{3}} + \frac{\eta}{2} (1-\xi) \frac{\hat{\sigma}^{2}f}{\hat{\sigma}\eta^{2}} + f \frac{\hat{\sigma}^{2}f}{\hat{\sigma}\eta^{2}} \Big[ \xi \Big(\frac{m+3}{6}\Big) - \Big(1-\xi\Big) \Big(\frac{m-3}{6}\Big) \log\left(1-\xi\right) \Big] + \frac{m}{3}\xi \Big[1-\left(\frac{\hat{\sigma}f}{\hat{\sigma}\eta}\right)^{2} \Big] \\ \lambda\xi\big(\theta+N\phi\big) = \xi\Big(1-\xi\Big) \frac{\hat{\sigma}^{2}f}{\hat{\sigma}\eta\hat{\sigma}\xi} + \frac{\xi}{3} (m-3)(1-\xi)\log\left(1-\xi\right) \Big[ \frac{\hat{\sigma}^{2}f}{\hat{\sigma}\eta^{2}} \frac{\hat{\sigma}f}{\hat{\sigma}\xi} - \frac{\hat{\sigma}f}{\hat{\sigma}\eta} \frac{\hat{\sigma}^{2}f}{\hat{\sigma}\eta\hat{\sigma}\xi} \Big], \end{split} \tag{1}$$

$$\begin{split} &\frac{\hat{\mathcal{O}}^{2}\theta}{\hat{\mathcal{O}}\eta^{2}} + \frac{\eta}{2}Pr(1-\xi)\frac{\hat{\mathcal{O}}\theta}{\hat{\mathcal{O}}\eta} + Pr\hat{f}\frac{\hat{\mathcal{O}}\theta}{\hat{\mathcal{O}}\eta} \Big[\xi\bigg(\frac{m+3}{6}\big) - (1-\xi)\bigg(\frac{m-3}{6}\big)\log(1-\xi)\bigg] - Pr\left(\frac{2m-3}{3}\right)\xi\theta\frac{\hat{\mathcal{O}}\theta}{\hat{\mathcal{O}}\eta} \\ &= \xi(1-\xi)Pr\frac{\hat{\mathcal{O}}\theta}{\hat{\mathcal{O}}\xi} + Pr(1-\xi)\log(1-\xi)\xi\Big(\frac{m-3}{3}\big) \Big[\frac{\partial\theta}{\partial\eta}\frac{\partial\hat{f}}{\partial\xi} - \frac{\partial\hat{f}}{\partial\eta}\frac{\partial\theta}{\partial\xi}\Big], \end{split} \tag{2}$$

$$\begin{split} &\frac{\hat{\sigma}^{2}\phi}{\hat{c}\eta^{2}} + \frac{\eta}{2} \text{Sc}(1-\xi) \frac{\hat{c}\phi}{\hat{c}\eta} + \text{S}\xi\dagger \frac{\hat{c}\phi}{\hat{c}\eta} \Big[\xi\bigg(\frac{m+3}{6}\big) - (1-\xi)\bigg(\frac{m-3}{6}\big) \log(1-\xi)\bigg] - \text{Sc}\left(\frac{2m-3}{3}\right)\xi\phi \frac{\hat{c}\xi}{\hat{c}\eta} \\ &= \xi(1-\xi)\text{Sc}\frac{\hat{c}\phi}{\hat{c}\xi} + \text{Sc}(1-\xi)\log(1-\xi)\xi\bigg(\frac{m-3}{3}\bigg) \left[\frac{\partial\phi}{\partial\eta}\frac{\partial\xi}{\partial\xi} - \frac{\partial\xi}{\partial\eta}\frac{\partial\phi}{\partial\xi}\right], \end{split} \tag{3}$$

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

**1. Introduction**

2 Mass Transfer - Advancement in Process Modelling

converge quadratically.

with existing results from literature.

**2. Governing equations**

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 quasilinearisation (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

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 differ‐ ences for discretizing in one of the independent variables and, instead, apply spectral collo‐ cation 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 colloca‐ tion 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

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

$$\begin{aligned} \partial\_{\tau} f\left(0, \xi\right) = 0, \quad \frac{\partial f\left(0, \xi\right)}{\partial \eta} = 0, \quad \theta\left(0, \xi\right) = 1, \quad \phi\left(0, \xi\right) = 1, \\\frac{\partial f\left(\infty, \xi\right)}{\partial \eta} = 1, \quad \theta\left(\infty, \xi\right) = 0, \quad \phi\left(\infty, \xi\right) = 0. \end{aligned} \tag{4}$$

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

$$\text{Re}\_{\perp}^{1/2}\mathbb{C}\_{f} = \mathbb{Z}\xi^{-1/2}f''(0,\xi),\tag{5}$$

$$\text{Re}\_{\perp}^{-1/2} \text{Nu} = -\xi^{-1/2} \theta'(0, \xi), \tag{6}$$

$$\text{Re}\_L^{-1/2} S \text{lt} = -\xi^{-1/2} \phi'(0, \xi). \tag{7}$$

#### **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 *ϕ*(*η*,*ξ*)

$$\mathbf{1}\Omega\_k \left[ \mathbf{F}, T, H \right] = \mathbf{0}, \quad \text{for} \quad k = \mathbf{1}, \mathbf{2}, \mathbf{3}, \tag{8}$$

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

$$\begin{split} F &= \left\{ f, \frac{\widehat{\partial}f}{\widehat{\partial}\eta}, \frac{\widehat{\partial}^2 f}{\widehat{\partial}\eta^2}, \frac{\widehat{\partial}^3 f}{\widehat{\partial}\eta^3}, \frac{\widehat{\partial}^2 f}{\widehat{\partial}\xi\widehat{\partial}\eta} \right\}, \\ T &= \left\{ \theta, \frac{\widehat{\partial}\theta}{\widehat{\partial}\eta}, \frac{\widehat{\partial}^2 \theta}{\widehat{\partial}\eta^2}, \frac{\widehat{\partial}\theta}{\widehat{\partial}\xi} \right\}, \\ H &= \left\{ \phi, \frac{\widehat{\partial}\phi}{\widehat{\partial}\eta}, \frac{\widehat{\partial}^2 \phi}{\widehat{\partial}\eta^2}, \frac{\widehat{\partial}\phi}{\widehat{\partial}\xi} \right\}. \end{split}$$

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

#### **LLM Algorithm**


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

$$\left[F \cdot \nabla\_{f} \mathfrak{Q}\_{1} \Big\lbrack \mathbb{F}\_{r'} T\_{r'} H\_{r} \Big\rbrack = \mathbb{F}\_{r} \cdot \nabla\_{f} \mathfrak{Q}\_{1} \Big\lbrack \mathbb{F}\_{r'} T\_{r'} H\_{r} \Big\rbrack - \mathfrak{Q}\_{1} \Big\lbrack \mathbb{F}\_{r'} T\_{r'} H\_{r} \Big\rbrack \Big\rbrack \tag{9}$$

On the New Bivariate Local Linearisation Method for Solving Coupled Partial Differential Equations in … http://dx.doi.org/10.5772/60444 5

$$T \cdot \nabla\_{\vartheta} \Omega\_2 \Big[ F\_{r+1}, T\_{r}, H\_r \Big] = T\_r \cdot \nabla\_{\vartheta} \Omega\_2 \Big[ F\_{r+1}, T\_{r}, H\_r \Big] - \Omega\_2 \Big[ F\_{r+1}, T\_{r}, H\_r \Big],\tag{10}$$

$$\left[H \cdot \nabla\_{\boldsymbol{\phi}} \Omega\_3 \Big[\boldsymbol{F}\_{r+1}, \boldsymbol{T}\_{r+1}, \boldsymbol{H}\_r\right] = \boldsymbol{H}\_r \cdot \nabla\_{\boldsymbol{\phi}} \Omega\_3 \Big[\boldsymbol{F}\_{r+1}, \boldsymbol{T}\_{r+1}, \boldsymbol{H}\_r\Big] - \Omega\_3 \Big[\boldsymbol{F}\_{r+1}, \boldsymbol{T}\_{r+1}, \boldsymbol{H}\_r\Big].\tag{11}$$

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

$$\begin{split} \nabla &= \left\{ \nabla\_{f'} \nabla\_{\boldsymbol{\theta}'} \nabla\_{\boldsymbol{\phi}} \right\}\_{'} \\ \nabla\_{f} &= \left\{ \frac{\partial}{\partial f'}, \frac{\partial}{\partial f'}, \frac{\partial}{\partial f''}, \frac{\partial}{\partial f'''}, \frac{\partial}{\partial f\_{\boldsymbol{\xi}}'} \right\}\_{'} \\ \nabla\_{\boldsymbol{\phi}} &= \left\{ \frac{\partial}{\partial \boldsymbol{\theta}'}, \frac{\partial}{\partial \boldsymbol{\theta}'}, \frac{\partial}{\partial \boldsymbol{\theta}'''}, \frac{\partial}{\partial \boldsymbol{\theta}\_{\boldsymbol{\xi}}} \right\}\_{'} \\ \nabla\_{\boldsymbol{\phi}} &= \left\{ \frac{\partial}{\partial \boldsymbol{\phi}'}, \frac{\partial}{\partial \boldsymbol{\phi}'}, \frac{\partial}{\partial \boldsymbol{\phi}'''}, \frac{\partial}{\partial \boldsymbol{\phi}\_{\boldsymbol{\xi}}} \right\}. \end{split}$$

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

$$a\_{1,r}f\_{r+1}^{\prime\prime\prime} + a\_{1,r}\left(\eta\_r\xi\right)f\_{r+1}^{\prime\prime} + a\_{2,r}\left(\eta\_r\xi\right)f\_{r+1}^{\prime} + a\_{3,r}\left(\eta\_r\xi\right)f\_{r+1} + a\_{4,r}\left(\eta\_r\xi\right)\frac{\hat{\sigma}\_{r+1}^{f}}{\hat{\sigma}\_{\xi}^{\xi}} = \mathcal{R}\_{1,r}\left(\eta\_r\xi\right),\tag{12}$$

$$\theta\_{r+1}{}^{\prime\prime} + b\_{1,r} \left(\eta\_{r} \xi \right) \theta\_{r+1}{}^{\prime} + b\_{2,r} \left(\eta\_{r} \xi \right) \theta\_{r+1} + b\_{3,r} \left(\eta\_{r} \xi \right) \frac{\partial \theta\_{r+1}}{\partial \xi} = \mathcal{R}\_{2,r} \left(\eta\_{r} \xi \right),\tag{13}$$

$$c\_{r+1}"+c\_{1,r}\left(\eta\_{r}\xi\right)\phi\_{r+1}"+c\_{2,r}\left(\eta\_{r}\xi\right)\phi\_{r+1}+c\_{3,r}\left(\eta\_{r}\xi\right)\frac{\partial\phi\_{r+1}}{\partial\xi}=\mathcal{R}\_{3,r}\left(\eta\_{r}\xi\right),\tag{14}$$

where the primes denote partial differentiation with respect to *η* and the coefficients *ai*,*r*,*bi*,*<sup>r</sup>* and *ci*,*<sup>r</sup>* (*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

1/2 1/2 = (0, ). *Re Sh <sup>L</sup>* xfx

The derivation of the solution method is described for a general system of non-linear PDEs in

where Ω1, Ω2 and Ω<sup>3</sup> are non-linear operators that represent the non-linear PDEs (1), (2) and

*fff f F f*

=, , , ,

ì ü ï ï ¶¶ ¶ í ý ï ï î þ ¶ ¶ ¶ ì ü ï ï ¶¶ ¶ í ý ï ï î þ ¶ ¶ ¶

h h h

q qq

h

f ff

q

f h

*T*

*H*

=, , , .

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

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

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

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

ë û ë ûë û ×Ñ W - W (9)

1 11 ,, = ,, ,, , *f rr r r f rr r rr r F FT H F FT H FT H* ×Ñ W é ù é ùé ù

=,, , , ,

ì ü ï ï ¶¶ ¶ ¶ í ý ï ï ¶ ¶ ¶ ¶ ¶ î þ

23 2 2 3

 x h

 x h

x h

this section. Consider a system of 3 coupled PDEs in *f* (*η*,*ξ*), *θ*(*η*,*ξ*) and *ϕ*(*η*,*ξ*)

**3. Derivation of solution method**

4 Mass Transfer - Advancement in Process Modelling

(3), respectively and *F* ,*T* ,*H* are defined as

linearisation method (LLM):

previous iteration. ↣ *Tr*+1.

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

**LLM Algorithm**

*Fr*+1.

equations (8) gives


, , = 0, = 1,2,3, *<sup>k</sup>* W é ù *F T H for k* ë û (8)

$$\begin{pmatrix} f\_{r+1}\begin{pmatrix} 0, \xi \end{pmatrix} = 0, & f\_{r+1}'\begin{pmatrix} 0, \xi \end{pmatrix} = 0, & \theta\_{r+1}\begin{pmatrix} 0, \xi \end{pmatrix} = 1, & \phi\_{r+1}\begin{pmatrix} 0, \xi \end{pmatrix} = 1,\tag{15}$$

$$\theta\_r f\_{r+1}'\left(\infty,\xi\_{\mathfrak{p}}^r\right) = 1,\quad \theta\_{r+1}\left(\infty,\xi\_{\mathfrak{p}}^r\right) = 0,\quad \phi\_{r+1}\left(\infty,\xi\_{\mathfrak{p}}^r\right) = 0.\tag{16}$$

Starting from given initial approximations, denoted by *f* <sup>0</sup>(*η*,*ξ*),*θ*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

$$f\left(\eta,\xi\right) \approx \sum\_{m=0}^{N\_x} \sum\_{j=0}^{N\_l} f\left(\tau\_{m'},\zeta\_j'\right) L\_m\left(\tau\right) L\_j\left(\zeta'\right),\tag{17}$$

where the functions *L <sup>m</sup>*(*τ*) and *L <sup>j</sup>* (*ζ*) are the well-known characteristic Lagrange cardinal polynomials. The function (17) interpolates *f* (*η*,*ξ*) at the collocation points (known as Cheby‐ shev-Gauss-Lobatto) defined by

$$\text{for } \tau\_i = \cos\left(\frac{\pi\,i}{N\_\pm}\right), \quad \mathcal{L}\_j = \cos\left(\frac{\pi\,j}{N\_\pm}\right), \quad i = 0, 1, \dots, N\_\pm; \ j = 0, 1, \dots, N\_\pm. \tag{18}$$

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 transforma‐ tions have been used to transform *η* ∈ 0,*η∞* and *ξ* ∈ 0,*L <sup>t</sup>* to *τ* ∈ −1,1 and *ζ* ∈ −1,1 , respec‐ tively. Here *η∞* is a finite value that is introduced to facilitate the application of the numerical method at infinity and *L <sup>t</sup>* 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

$$\mathbf{a}\_{i}\mathbf{F}\_{r+1,i} + \mathbf{a}\_{4,i}\sum\_{j=0}^{N\_{t}} \mathbf{d}\_{i,j} \mathbf{D}\mathbf{F}\_{r+1,j} + \mathbf{a}\_{5,i}\sum\_{j=0}^{N\_{t}} \mathbf{d}\_{i,j}\mathbf{F}\_{r+1,j} = \mathbf{R}\_{1,i} \tag{19}$$

$$\mathbf{B}\_i \mathbb{Q}\_{r+1,i} + b\_{3,i} \sum\_{j=0}^{N\_t} d\_{i,j} \mathbb{Q}\_{r+1,j} = \mathbf{R}\_{2,i} \,\prime \tag{20}$$

$$\mathbf{C}\_{i}\mathbf{F}\_{r+1,i} + c\_{3,i}\sum\_{j=0}^{N\_{\rm t}}\mathbf{d}\_{i,j}\mathbf{F}\_{r+1,j} = \mathbf{R}\_{3,i'} \tag{21}$$

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

$$\begin{split} &F\_{r+1,i} = \left[f\_{r+1,i}(\tau\_{0})\_{\prime }f\_{r+1,i}(\tau\_{1})\_{\prime }, \dots, f\_{r+1,i}(\tau\_{N\_{\underline{x}}})\right]^{\mathrm{T}}, \\ &Q\_{r+1,i} = \left[\Theta\_{r+1,i}(\tau\_{0})\_{\prime }\Theta\_{r+1,i}(\tau\_{1})\_{\prime }, \dots, \Theta\_{r+1,i}(\tau\_{N\_{\underline{x}}})\right]^{\mathrm{T}}, \\ &F\_{r+1,i} = \left[\phi\_{r+1,i}(\tau\_{0})\_{\prime }\phi\_{r+1,i}(\tau\_{1})\_{\prime }, \dots, \phi\_{r+1,i}(\tau\_{N\_{\underline{x}}})\right]^{\mathrm{T}}, \\ &A\_{i} = D^{3} + a\_{1,i}D^{2} + a\_{2,i}D + a\_{3,i}, \\ &B\_{i} = D^{2} + b\_{1,i}D + b\_{2,i}, \\ &C\_{i} = D^{2} + c\_{1,i}D + c\_{2,i}, \\ \end{split}$$

*am*,*r*(*ξ<sup>i</sup>* )(*m*=1,2,3,4,5) is the diagonal matrix of the vector *am*,*r*(*τ*0),*am*,*r*(*τ*1),…,*am*,*r*(*τNx* ) *<sup>T</sup>* . The matrices *b*1,*<sup>i</sup>* ,*b*2,*<sup>i</sup>* ,*c*1,*<sup>i</sup>* ,*c*2,*<sup>i</sup>* are also diagonal matrices corresponding to *b*1,*r*,*b*2,*r*,*c*1,*r*,*c*2,*<sup>r</sup>* evaluated at collocation points. Additionally, *R*1,*<sup>i</sup>* , *R*2,*<sup>i</sup>* and *R*3,*<sup>i</sup>* are (*Nx* + 1)×1 vectors corresponding to *R*1,*<sup>r</sup>*, *R*2,*<sup>r</sup>* and *R*3,*<sup>r</sup>*, 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:

$$
\begin{bmatrix} A\_{0,0} & A\_{0,1} & \cdots & A\_{0,N\_t} \\ A\_{1,0} & A\_{1,1} & \cdots & A\_{1,N\_t} \\ \vdots & \vdots & \ddots & \vdots \\ A\_{N\_t,0} & A\_{N\_t,1} & \cdots & A\_{N\_t,N\_t} \end{bmatrix} \begin{bmatrix} F\_{r+1,0} \\ F\_{r+1,1} \\ \vdots \\ F\_{r+1,N\_t} \end{bmatrix} = \begin{bmatrix} R\_{1,0} \\ R\_{1,1} \\ \vdots \\ R\_{1,N\_t} \end{bmatrix},\tag{22}
$$

where

Starting from given initial approximations, denoted by *f* <sup>0</sup>(*η*,*ξ*),*θ*0(*η*,*ξ*) and *ϕ*0(*η*,*ξ*), equations (12) – (14) are solved iteratively for *r* =1,2,…, until approximate solutions that are consistent

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

( ) ( ) ( ) ( )

= cos , = cos , = 0,1, , ; = 0,1, , . *i j x t*

p  t z

,, ,

polynomials. The function (17) interpolates *f* (*η*,*ξ*) at the collocation points (known as Cheby‐

*<sup>i</sup> <sup>j</sup> i Nj N N N*

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 transforma‐ tions have been used to transform *η* ∈ 0,*η∞* and *ξ* ∈ 0,*L <sup>t</sup>* to *τ* ∈ −1,1 and *ζ* ∈ −1,1 , respec‐ tively. Here *η∞* is a finite value that is introduced to facilitate the application of the numerical

Substituting equation (17) in equation (12)-(14) and making use of the derivatives formulas for

1, 4, , 1, 5, , 1, 1, =0 =0

> 1, 3, , 1, 2, =0

1, 3, , 1, 3, =0

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

*N N t t i r i i ij r j i ij r j i j j*

> *Nt i r i i ij r j i j*

*Nt i r i i ij r j i j*

= ,

*B bd R* + + Q Q <sup>+</sup> å (20)

*C cd R* + + F F <sup>+</sup> å (21)

*A F a d DF a d F R* + ++ + + å å (19)

= ,

= ,

*mjm j*

 t z

» åå (17)

(*ζ*) are the well-known characteristic Lagrange cardinal

K K (18)

=0 =0

*m j f f LL*

*x t*

Lagrange functions at Gauss-Lobatto points given in [2, 8], results in

 z æ ö æö ç ÷ ç÷ è ø èø

p

method at infinity and *L <sup>t</sup>* is the largest value of *ξ*.

hx

*N N x t*

to within a certain tolerance level are obtained.

solution for *f* (*η*,*ξ*) takes the form

6 Mass Transfer - Advancement in Process Modelling

where the functions *L <sup>m</sup>*(*τ*) and *L <sup>j</sup>*

shev-Gauss-Lobatto) defined by

t

$$\begin{aligned} A\_{i,j} &= A\_i + a\_{4,i}d\_{i,i}D + a\_{5,i}d\_{i,i}I\_{,i} \\ A\_{i,j} &= a\_{4,i}d\_{i,i}D + a\_{5,i}d\_{i,i}I\_{,i} \quad \text{where } i \neq j\_{,i} \end{aligned}$$

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

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

$$
\begin{bmatrix}
\mathcal{B}\_{0,0} & \mathcal{B}\_{0,1} & \cdots & \mathcal{B}\_{0,N\_t} \\
\mathcal{B}\_{1,0} & \mathcal{B}\_{1,1} & \cdots & \mathcal{B}\_{1,N\_t} \\
\vdots & \vdots & \ddots & \vdots \\
\mathcal{B}\_{N\_t,0} & \mathcal{B}\_{N\_t,1} & \cdots & \mathcal{B}\_{N\_t,N\_t}
\end{bmatrix}
\begin{bmatrix}
\mathcal{Q}\_{r+1,0} \\
\mathcal{Q}\_{r+1,1} \\
\vdots \\
\mathcal{Q}\_{r+1,N\_t}
\end{bmatrix} = 
\begin{bmatrix}
\mathcal{R}\_{2,0} \\
\mathcal{R}\_{2,1} \\
\vdots \\
\mathcal{R}\_{2,N\_t}
\end{bmatrix},\tag{23}
$$

$$
\begin{bmatrix}
\mathbb{C}\_{0,0} & \mathbb{C}\_{0,1} & \cdots & \mathbb{C}\_{0,N\_{\ell}} \\
\mathbb{C}\_{1,0} & \mathbb{C}\_{1,1} & \cdots & \mathbb{C}\_{1,N\_{\ell}} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbb{C}\_{N\_{\ell},0} & \mathbb{C}\_{N\_{\ell},1} & \cdots & \mathbb{C}\_{N\_{\ell},N\_{\ell}}
\end{bmatrix}
\begin{bmatrix}
\mathbb{F}\_{r+1,0} \\
\mathbb{F}\_{r+1,1} \\
\vdots \\
\mathbb{F}\_{r+1,N\_{\ell}-1}
\end{bmatrix} = \begin{bmatrix}
\mathbb{R}\_{3,0} \\
\mathbb{R}\_{3,1} \\
\vdots \\
\mathbb{R}\_{3,N\_{\ell}}
\end{bmatrix} \tag{24}
$$

where

$$\begin{aligned} B\_{i,j} &= B\_i + b\_{3,i} d\_{i,i} I\_{\prime} \\ B\_{i,j} &= b\_{3,i} d\_{i,j} I\_{\prime} \quad \text{when } i \neq j, \\ C\_{i,j} &= C\_i + c\_{3,i} d\_{i,i} I\_{\prime} \\ C\_{i,j} &= c\_{3,i} d\_{i,j} I\_{\prime} \quad \text{when } i \neq j, \end{aligned}$$

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,

$$E\_f = \max\_{0 \le i \le N\_{\tilde{t}}} \left\| F\_{r+1,i} - F\_{r,i} \right\|\_{\alpha'} \quad E\_\theta = \max\_{0 \le i \le N\_{\tilde{t}}} \left\| \mathbb{Q}\_{r+1,i} - \mathbb{Q}\_{r,i} \right\|\_{\alpha'} \quad E\_\phi = \max\_{0 \le i \le N\_{\tilde{t}}} \left\| F\_{r+1,i} - F\_{r,i} \right\|\_{\alpha} \tag{25}$$

#### **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 com‐ putations 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 10<sup>−</sup><sup>7</sup> were determined using equation (25).

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.

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

On the New Bivariate Local Linearisation Method for Solving Coupled Partial Differential Equations in … http://dx.doi.org/10.5772/60444 9

**Figure 1.** Effect of *λ* and *Pr* on the velocity profile *f* ′ (*η*,*ξ*)

where

between values of successive iterations,

8 Mass Transfer - Advancement in Process Modelling

**4. Results and discussion**

results to within a tolerance level of 10<sup>−</sup><sup>7</sup>

results generated using the proposed BSLLM.

, 3, , , 3, , , 3, , , 3, ,

+

*ij i i ii ij i ij ij i i ii ij i ij*

*C C cdI*

+

1, , 1, , 1, , 00 0 <sup>=</sup> max , = max , = max . *<sup>f</sup> r i ri r i ri r i ri i N i N i N tt t*

¥ ¥¥ £ £ £ £ £ £

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 com‐ putations 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

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.

Figs. 3, 4 and 5 show the effect of buoyancy ratio *N* on the velocity, temperature and concen‐ tration profiles, respectively, when *Pr* =0.7, *Sc* =0.94 and *λ* =5 at *ξ* =0.5. The graphs are qualita‐ tively similar to those reported in [7]. These figures are a qualitative validation of the numerical

*E FF E E* +++ q

*B B bdI*

= ,

= ,

=, ,

¹

¹

 f


were determined using equation (25).

*B b d I when i j*

*C c d I when i j*

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

=, ,

**Figure 2.** Effect of *λ* and *Pr* on the temperature profile *θ*(*η*,*ξ*)

**Figure 3.** Effect of *N* velocity profile *f* ′ (*η*,*ξ*)

**Figure 4.** Effect of *N* on the temperature profile *θ*(*η*,*ξ*)

On the New Bivariate Local Linearisation Method for Solving Coupled Partial Differential Equations in … http://dx.doi.org/10.5772/60444 11

**Figure 5.** Effect of *N* concentration profile *ϕ*(*η*,*ξ*)

**Figure 3.** Effect of *N* velocity profile *f* ′

10 Mass Transfer - Advancement in Process Modelling

**Figure 4.** Effect of *N* on the temperature profile *θ*(*η*,*ξ*)

(*η*,*ξ*)

**Figure 6.** Effect of *λ* on the skin friction coefficient *Cf ReL* 1

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 lineari‐ sation 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.

#### **Appendix**

#### **Definition of coefficients**

$$a\_{1,r} = \frac{\eta}{2}(1-\xi) + \beta\_1(\xi)f\_r - \beta\_2(\xi)\frac{\partial f\_r}{\partial \xi}\_{\xi} \tag{26}$$

$$a\_{2,r} = -\frac{2m}{3} \xi f\_r' + \beta\_2(\xi) \frac{\partial f\_r'}{\partial \xi}\_{\mathfrak{T}}' \tag{27}$$

$$a\_{3,r} = \beta\_1(\xi) f\_{r'}^{\prime} \tag{28}$$

$$a\_{4,r} = -\xi(1-\xi) + \beta\_2(\xi)f'\_{r,r} \tag{29}$$

$$a\_{\mathfrak{s}\_{r'}} = -\mathcal{B}\_2(\xi) f\_{r'}^{''} \tag{30}$$

$$R\_{1,r} = \beta\_1(\xi) f\_r f\_r'' - \frac{m}{3} \xi \left[ 1 + (f\_r')^2 \right] - \lambda \xi (\theta\_r + N\phi\_r) - \beta\_2(\xi) \left[ f\_r' \frac{\partial f\_r}{\partial \xi} - f\_r' \frac{\partial f\_r'}{\partial \xi} \right] \tag{31}$$

$$\beta\_1(\xi) = \xi \left(\frac{m+3}{6}\right) - (1-\xi) \left(\frac{m-3}{6}\right) \log(1-\xi),\tag{32}$$

On the New Bivariate Local Linearisation Method for Solving Coupled Partial Differential Equations in … http://dx.doi.org/10.5772/60444 13

$$
\beta\_2 \beta\_\sharp (\xi) = \xi (1 - \xi) \log(1 - \xi) \frac{(m - 3)}{3} \,\mathrm{},\tag{33}
$$

$$b\_{1,r} = \frac{\eta}{2} Pr(1-\xi) + Pr\beta\_1(\xi)f\_r - Pr\beta\_2(\xi)\frac{\hat{\sigma}f\_r}{\hat{\sigma}\xi},\tag{34}$$

$$b\_{2,r} = -Pr\frac{(2m-3)}{3} \xi f\_{r,r}' \tag{35}$$

$$b\_{3,r} = -\xi (1 - \xi) Pr + Pr \beta\_2(\xi) f'\_{r'} \tag{36}$$

$$c\_{1,r} = \frac{\eta}{2} \text{Sc}(1-\xi) + \text{Sc}\mathcal{B}\_1(\xi)f\_r - \text{Sc}\mathcal{B}\_2(\xi)\frac{\hat{c}f\_r}{\hat{c}\xi} \tag{37}$$

$$\mathcal{L}\_{2,r} = -\mathrm{Sc}\frac{(2m-3)}{3}\xi f\_{r}^{\prime}{}\_{\prime} \tag{38}$$

$$\mathcal{L}c\_{3,r} = -\xi (1-\xi) \mathcal{S}c + \mathcal{S}c \mathcal{B}\_2(\xi) f'\_{r},\tag{39}$$

#### **Author details**

#### Sandile Motsa\*

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

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 lineari‐ sation 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

1, <sup>1</sup> <sup>2</sup> = (1 ) ( ) ( ) , <sup>2</sup>


3, 1 = () , *r r a f* b x

4, <sup>2</sup> = (1 ) ( ) , *r r a f* - -+ x x bx

> 5, 2 = () , *r r a f* b x

1, 1 <sup>2</sup> = () 1 ( ) ( ) () , <sup>3</sup>

3 3 ( )= (1 ) log(1 ), 6 6 *m m*

x

æö æö + - - - - ç÷ ç÷

 f b x*N ff*

lx q

*r r r r r r r r <sup>m</sup> f f R ff f*

2

 x

1

bx x  bx x

*r r <sup>f</sup> a f*

> 2, 2 <sup>2</sup> <sup>=</sup> () , <sup>3</sup>

*r r <sup>m</sup> <sup>f</sup> a f* x bx x

x bx

h

*r*

¶ ¢ - +¢ ¶ (27)

¢¢ (28)

¢¢ (30)

 x

*r r*

x

é ù ¶ ¶ ¢ ¢¢- + - +- - é ù ¢ ¢¢ ¢ ê ú ë û ¶ ¶ ë û (31)

x

èø èø (32)

¢ (29)

¶ (26)

¶

*r*

trend in the skin friction coefficient for near the stagnation region can be observed.

differential equations with fluid mechanics applications.

**5. Conclusion**

12 Mass Transfer - Advancement in Process Modelling

**Appendix**

**Definition of coefficients**

bx

Address all correspondence to: sandilemotsa@gmail.com

School of Mathematics, Statistics and Computer Science, University of KwaZulu-Natal, Private Bag X, Scottsville, Pietermaritzburg, South Africa

#### **References**

