**American Journal of Computational Mathematics**

Vol.05 No.02(2015), Article ID:57239,16 pages

10.4236/ajcm.2015.52013

Numerical Solution of MHD Flow of Micropolar Fluid with Heat and Mass Transfer towards a Stagnation Point on a Vertical Plate

N. T. El-Dabe^{1}, A. Y. Ghaly^{1}, R. R. Rizkallah^{1}, K. M. Ewis^{2}, A. S. Al-Bareda^{1*}

^{1}Department of Mathematics, Faculty of Education, Ain Shams University, Cairo, Egypt

^{2}Department of Engineering Mathematics and Physics, Faculty of Engineering, El-Fayoum University, El-Fayoum, Egypt

Email: ^{*}ameen_azeez@hotmail.com

Copyright © 2015 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

Received 21 April 2015; accepted 14 June 2015; published 18 June 2015

ABSTRACT

The paper investigates the numerical solution of problem of magnetohydrodynamic (MHD) micropolar fluid flow with heat and mass transfer towards a stagnation point on a vertical plate. In this study, we consider both strong concentrations (n = 0) and weak concentrations (n = 1/2). The governing equations have been transformed into nonlinear ordinary differential equations by applying the similarity transformation and have been solved numerically by using the finite difference method (FDM) and analytically by using (DTM). The effects of various governing parameters, namely, material parameter, radiation parameter, magnetic parameter, Prandtl number, Schmidt number, chemical reaction parameter and Soret number on the velocity, microrotation, temperature and concentration have been computed and discussed in detail through some figures and tables. In order to verify the accuracy of the present results, we have compared these results with the analytical solutions by using the differential transform method (DTM) and the multi-step differential transform method (MDTM). It is observed that this approximate numerical solution is in good agreement with the analytical solution.

**Keywords:**

Finite Difference Method (FDM), Differential Transform Method (DTM), Micropolar Fluid, MHD, Heat and Mass Transfer, Stagnation Flow, Chemical Reaction, Radiation

1. Introduction

The hydromagnetic stagnation point flows past heated or cooled bodies have attracted many researchers due to their applications in most of the engineering and natural processes. Examples include plasma studies, blood flow problems, MHD generators, and the cooling of an infinite metallic plate in a cooling bath. The classical two dimensional stagnation point flow on a flat plate first studied by Hiemenz [1] was extended to axisymmetric case by Homann [2] . The problem of two dynamic impinging stagnation flows of two fluids of different densities was investigated by Wang [3] and it was shown that the flow field depends heavily on the ratios of viscosity and density. The influence of an external magnetic field on Hiemenz flow was studied by Ariel [4] . Two dimensional MHD steady stagnation point flow towards a stretching surface was analyzed by Mahapatra and Gupta [5] . Ishak et al. [6] theoretically studied the similarity solutions for the steady magnetohydrodynamic flow towards a stagnation point on a vertical surface immersed in an incompressible micropolar fluid. Lokendra Kumar [7] has analyzed the effect of MHD flow of micropolar fluid towards a stagnation point on a vertical stretching sheet. Ishak et al. [8] investigated the steady two-dimensional stagnation point flow of an incompressible micropolar fluid towards a vertical stretching sheet.

The theory of micropolar fluids which is originally formulated by Eringen [9] can be used to explain the flow of crystals, animal blood, paints, polymers, etc. The theory introduces new material parameters, an additional in dependent vector field the microrotation―and new constitutive equations which must be solved simultaneously with the usual equations for Newtonian flow. Ramachandran et al. [10] studied laminar mixed convection in two-dimensional stag nation flows around surfaces. He considered both cases of an arbitrary wall temperature and arbitrary surface heat flux variations and found that they are versed flow developed in the buoyancy opposing flow region, and dual solutions are found to certain range of the buoyancy parameter. Hassanien et al. [11] extended Ramachandran’s work to micropolar fluid. They considered both assisting and opposing flows, but the existence of dual solutions was not reported [12] . Devi et al. extended the problem posed by Ramachandran et al. [10] to the unsteady case, and they found that dual solution exists for a certain range of the buoyancy parameter when the flow is opposing. Similar problem for steady and unsteady cases, for a vertical surface immersed in a micropolar fluid was investigated by Lok et al. [13] [14] . Existence of dual solutions was reported in [13] only for the opposing flow regime. Mahapatra et al. [15] examined heat transfer in stagnation-point flow on stretching sheet with viscous dissipation effect. Attia [16] studied the hydromagnetic stagnation point flow on porous stretching sheet with suction and injection. El-dabe et al. [17] studied the peristaltic motion of incompressible micropolar fluid through a porous medium in a two-dimensional channel under the effects of heat absorption and chemical reaction in the presence of magnetic field. Effects of thermal radiation on magnetohydrodynamic (MHD) flow of a micropolar fluid towards a stagnation point on a vertical Plate were studied by Olanrewaju et al. [18] . Pop et al. [19] examined the flow over a stretching sheet near a stagnation point taking radiation effect. Olanrewaju and Adesanya [20] studied the effects of radiation and viscous dissipation on stagnation flow of a micropolar fluid towards a vertical permeable surface. Manjoolatha et al. [21] studied the radiation and mass transfer effects on MHD flow of a micropolar fluid towards a stagnation point on a vertical stretching sheet.

The main objective of the present work is to study the numerical solutions of problem of magnetohydrodynamic (MHD) micropolar fluid flow with heat and mass transfer towards a stagnation point on a vertical plate. In this study, we consider both strong concentrations (n = 0) and weak concentrations (n = 1/2). The governing equations have been transformed into nonlinear ordinary differential equations by applying the similarity transformation and have been solved numerically by using the finite difference method (FDM) and analytically by using (DTM). The effects of various parameters of the problem on fluid velocity, microrotation, temperature and concentration have been computed and discussed in detail through some figures and tables. Comparisons with analytical solution by using (DTM) of present work are performed and show that the present results apply with results of analytical solution, as well as the comparisons with previous published results and are found to be in good agreement.

2. Mathematical Formulation

We consider the steady, two-dimensional flow of an incompressible electrically conducting micropolar fluid near the stagnation point on a vertical heated stretching sheet placed in the (y = 0) of a cartesian coordinate system with the x-axis along the sheet as shown in Figure 1. The fluid occupies the half plane (y > 0). It is assumed that the velocity of the flow external to the boundary layer U, the temperature of the sheet and concentration of sheet are proportional to the distance x from the stagnation point, i.e. where a, b and c are positive constants, with being the uniform temperature of the fluid and with being the uniform concentration of the fluid. A uniform magnetic field of strength is assumed to be applied in the positive y-direction normal to the plate. The magnetic Reynolds number of the flow is taken to be small enough so that the induced magnetic field is negligible.

Under the usual boundary layer approximation, the governing equations are given as follows:

(1)

(2)

(3)

(4)

(5)

The boundary conditions are as follows:

(6)

where u and v are the velocity components along the x and y direction, g is the acceleration due to gravity, T is the fluid temperature, is the plate temperature, is the ambient temperature, C is the fluid concentration in the boundary layer, is the plate concentration, is the ambient fluid concentration, is the kine-

Figure 1. Physical model and coordinate system.

matic viscosity of the fluid, is the thermal expansion coefficient, is the coefficient of expansion with concentration, is the magnetic field of constant strength, is the electric conductivity of the fluid, is the radiative heat flux, is the specific heat at constant pressure of the fluid, is the mass diffusion coefficient, is the mean fluid temperature, is the thermal diffusivity ratio, is concentration susceptibility, and is the chemical reaction coefficient. Furthermore, are respectively the dynamic viscosity, vortex viscosity (or the microrotationviscosity), fluid density, microinertia density, microrotation vector (or angular velocity), spingradient viscosity and thermal diffusivity.

We follow the work of many authors by assuming that where is the

material parameter. This assumption is involved to allow the field of equations predicts the correct behavior in the limiting case when the microstructure effects become negligible and the total spin N reduces to the angular velocity (see Ishak et al. [6] ).

By using the Rosseland approximation (Brewster [22] ), the radiative heat flux is given by given by

(7)

where is the Stefan-Boltzmann constant and is the mean absorption coefficient. It should be noted that by using the Rosseland approximation, the present analysis is limited to optically thick fluids. If temperature differences within the flow are sufficiently small, then Equation (7) can be linearized by expanding into the Taylor series about, which after neglecting higher order terms takes the form

(8)

In view of Equations (7) and (8), Equation (4) reduces to

(9)

where k is the constant thermal conductivity.

3. Method of Solutions

To solve system of Equations (1)-(5), we will consider the following similarity transformations:

(10)

where is the independent similarity variable, the dimensionless stream function, the dimensionless micorotation, the dimensionless temperature and the dimensionless concentration. Further, is the stream function which is defined in the usual way

(11)

where prime denotes differentiation with respect to.

By using Equation (10) in the Equations (1)-(5) and (9), we get the following ordinary differential equations in dimensionless form:

(12)

(13)

(14)

(15)

where is the magnetic parameter, is the thermal buoyancy parameter, sign has the

same meaning as in Equation (2). For ג > 0, buoyancy forces act in the direction of the mainstream and fluid is accelerated in the manner of a favorable pressure gradient (assisting flow). When ג < 0, buoyancy forces oppose

the motion, retarding the fluid in the boundary layer, acting as an adverse pressure gradient (opposing flow), is the solutal buoyancy parameter, is the local thermal Grashof number, is the local soutal Groshof number, is the local Reynolds number, is the Eckert number, is the Prandtl number, , is the radiation parameter, is the Schmidt number, is the Chemical reaction parameter, is the Dufour number and is the Soret number.

Also, the subjected boundary conditions Equation (6) will now take the form:

(16)

The skin-friction coefficient, local Nusselt number and Sherwood number are defined as:

(17)

In which are given by

(18)

Substituting by similarity transformations and Equation (18) into Equation (17), we get

(19)

3.1. The Differential Transform Method

The differential transformation of an analytical function for one variable is defined as [23] .

(20)

where f(t) is the original function and F(k) is the transformed function.

The differential inverse transformation of F(k) is defined as:

(21)

Combining Equations (20) and (21), we obtain

(22)

From Equations (20)-(22), it can be seen that the differential transformation method is derived from Taylor's series expansion, but the method does not calculate the derivatives representatively. However, the relative derivatives are calculated by an iterative way which is described by the transformed equations of the original function. For implementation purposes, the function f(t) is expressed by a finite series and Equation (21) can be written as

(23)

By Equation (20), the following theorems can be deduced:

Theorem 1. If then

Theorem 2. If then

Theorem 3. If then

Theorem 4. If then

Theorem 5. If then

Theorem 6. If then

Theorem 7. If then

3.2. Basic Concepts of the Multi-Step Differential Transform Method (MDTM)

When the DTM is used for solving differential equations with the boundary condition at infinity or problems that have highly non-linear behavior, the obtained results were found to be incorrect (when the boundary-layer variable go to infinity, the obtained series solutions are divergent). Besides that, power series are not useful for large values of the independent variable.

To overcome this shortcoming, the MDTM that has been developed for the analytical solution of the differential equations is presented in this section. For this purpose, the following non-linear initial-value problem is considered:

(24)

subject to the initial conditions, for.

Let [0, T] be the interval over which we want to find the solution of the initial-value problem (24). In actual applications of the DTM, the approximate solution of the initial value problem (24) can be expressed by the following finite series:

(25)

The multi-step approach introduces a new idea for constructing the approximate solution. Assume that the interval [0, T] is divided into M subintervals,
of equal step size h = (T/M) by using the nodes t_{m} = mh. The main ideas of the MDTM are as follows. First, we apply the DTM to Equation (24) over the interval [0, t_{1}], we will obtain the following approximate solution:

(26)

using the initial conditions. For and at each subinterval we will use the initial conditions and apply the DTM to Equation (23) over the interval, where

in Equation (24) is replaced by. The process is repeated and generates a sequence of approximate solutions, for the solution f(t):

(27)

where. In fact, the MDTM assumes the following solution:

(28)

The new algorithm, MDTM, is simple for computational performance for all values of h. It is easily observed that if the step size h = T, then the MDTM reduces to the classical DTM. As we will see in the next section, the main advantage of the new algorithm is that the obtained series solution converges for wide time regions and can approximate non-chaotic or chaotic solutions.

3.3. Analytical Solution Be the MDTM

By applying the MDTM to Equations (12)-(15), the following recursive relations in each sub-domain, is given.

(29)

(30)

(31)

(32)

where are the differential transform of

The differential transformed boundary conditions in Equation (16) to:

(33)

where are constants. These constants are computed from the boundary condition.

Moreover, substituting Equation (33) into Equations (29)-(32) and by using the recursive method, we can calculate other values of.

Hence, substituting all, into Equation (23), we obtain series solutions.

The velocity, microrotation, temperature and concentration distributes achieved with the aid of MATHAMATICA application software.

3.4. Numerical Solution

The governing ordinary differential Equations (10)-(12) are being highly nonlinear and difficult to solve analytically. To solve these coupled equations, we use a finite difference based numerical algorithm. Following Ashraf et al. [24] , we reduce the order of Equation (12) by one with the help of the substitution:

(34)

Equations (12)-(15) in view of Equation (34) can be written as:

(35)

(36)

(37)

(38)

The boundary conditions:

(39)

Now, we can write Equations (35)-(38) in the following finite difference method form:

(40)

(41)

(42)

(43)

The velocity, microrotation, temperature and concentration distributes at all interior nodal points computed by successive applications of the above finite difference equations and these are achieved with the aid of MATLAB application software.

4. Results and Discussion

In this section, we present our results in tabular and graphical form is shown in Figures 2-16 to illustrate the influence of physical parameters such as, material parameter K, magnetic parameter M, thermal buoyancy ג, solutal buoyancy parameter, radiation parameter, Prandtl number, Dufour number, Schmidtnumber, chemical reaction and Soret number, on the velocity, microrotation, temperature and concentration when n = 0 and n = 1/2 for both the cases of assisting and opposing flows.

Figure 2 and Figure 3 illustrate the effect of the material parameter K on the normal velocity profiles. It is seen that the normal velocity profile decreases with an increasing in the material parameter K for both assisting and opposing flows. It is clear from the figure that the normal velocity profile is higher for assisting flow when compared with the normal velocity profile for opposing flow. Figures 4-7 represent the velocity profiles for the flow parameters. It is observed that the velocity profile increases with increasing the magnetic parameter M, while it decreases with increasing the material parameter K for both assisting and opposing flows. It is clear from the figures that the velocity is higher for assisting flow when compared with the velocity for opposing flow. Figures 8-10 represent the microrotation distribution for the flow parameters. Figure 8 illustrate the effect of the material parameter K on the microrotation when n = 1/2. It is clear that the microrotation increases with an increase K for assisting flow. For opposing flow the microrotation increases with an increase K in the vicinity of the surface but the reverse happens away from surface when n = 1/2 as shown in Figure 9, while the microrotation decreases with an increase K for both assisting and opposing flows when n = 0 as shown in Figure 10.

The effect of the flow parameters on the temperature field is represented by Figures 11-13. It is noted that temperature increases with an increase R_{a} and D_{u}, while it decrease with increasing P_{r} for both assisting and op-

Figure 2. Normal velocity profiles for some values of K when n = 0.5.

Figure 3. Normal velocity profiles for some values of K when n = 0.

Figure 4. Velocity profiles for some values of M when n = 0.5.

Figure 5. Velocity profiles for some values of M when n = 0.

Figure 6. Velocity profiles for some values of K when n = 0.5.

posing flows. It is clear from the figures that the temperature is higher for opposing flow when compared with the temperature for assisting flow. Figures 14-16 represent the concentration profiles for the flow parameters. It is noted that concentration increases with an increases S_{0}, while it decrease with the increases S_{c} and K_{r} for both assisting and opposing flows. It is clear from the figures that the concentration is higher for opposing flow when

Figure 7. Velocity profiles for some values of K when n = 0.

Figure 8. Angular velocity profiles for some values of K when n = 0.5.

Figure 9. Angular velocity profiles for some values of K when n = 0.5.

compared with the concentration for assisting flow.

In order to verify the accuracy of the numerical solution of present work by using (FDM) we have compared

Figure 10. Angular velocity profiles for some values of K when n = 0.

Figure 11. Temperature profiles
for some values of P_{r}.

Figure 12. Temperature profiles
for some values of R_{a}.

these results for with Ramachandran et al. [10] , Lok et al. [13] , Ishak et al. [6] and Olanrewaju et al. [18] . It is observed that this numerical solution in excellent agreement with the results previous obtained values as shown in Table 1 and Table 2.

Finally, Tables 3-5 show comparison between numerical solution by using (FDM) and the analytical solu-

Figure 13. Temperature profiles
for some values of D_{u}.

Figure 14. Concentration profiles
for some values of S_{c}.

Figure 15. Concentration profiles
for some values of K_{r}.

tion using (DTM). It is observed that this approximate numerical solution is in excellent agreement with the results of the analytical solution of by using (DTM).

Figure 16. Concentration profiles for some values of S0.

Table 1. Comparison between the previous obtained values of
by Ramachandran et al. [10] , Lok et al. [13] , Ishak et al. [6] Olanrewaju et al. [18] and our results for K = 0, M = 0, ג = 1 and R_{a} = 0.

Table 2. Comparison between the previous obtained values of
by Ramachandran et al. [10] , Lok et al. [13] , Ishak et al. [6] Olanrewaju et al. [18] and our results for K = 0, M = 0, ג = 1 and R_{a} = 0.

Table 3. Comparison between the previous obtained values of where ג = 1,

Table 4. Comparison between the previous obtained values of where ג = 1,

Table 5. Comparison between the previous obtained values of where ג = 1,

5. Conclusions

In this work, we have obtained the numerical solution of problem of magnetohydrodynamic (MHD) micropolar fluid flow with heat and mass transfer towards a stagnation point on a vertical plate. In this study, we consider both strong concentrations (n = 0) and weak concentrations (n = 1/2). The resulting partial differential equations which describe the problem are transformed into ordinary differential equations by using a similarity transformation and then solved numerically by using the finite difference method (FDM) and analytically by using DTM. A representative set of numerical results for velocity, microrotation, temperature and concentration profiles is presented graphically and discussed. The figures and tables clearly show that the results by using (FDM) are in excellent agreement with previously published works and with the results of the analytical solution by using DTM. The important results for this study are summarized as follows:

・ The normal velocity profile decreases with the increase of the material parameter K for both assisting and opposing flows.

・ The velocity distribution increases with the increase of M, while it decreases with the increase of K for both assisting and opposing flows.

・ The microrotation distribution increases with increase of K for assisting flow, while it increases with an increase of K in the vicinity of the surface but the reverse happens away from surface for opposing flows, when n = 1/2.

・ The microortation decreases with an increase K for both assisting and opposing flows when n = 0.

・ The temperature distribution increases with the increase of R_{a} and D_{u}, while it decreases with the increases of P_{r} for both assisting and opposing flows.

・ The concentration distribution increases with increase of S_{0}, while it decreases with the increases of S_{c} and K_{r} for both assisting and opposing flow.

References

- Hiemenz, K. (1911) Die grenzschicht an einem in dem gleichformingen ussigkeitsstrom eingetauchten gerade kreiszylinder. Dingler Polytechnic Journal, 326, 321-340.
- Homann, F. (1936) Der einuss grosser zahigkeit bei der stromung um den zylinder und umdie kugel. Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 16, 153-164. http://dx.doi.org/10.1002/zamm.19360160304
- Wang, C.Y. (1987) Impinging Stagnation Flows. Physics of Fluids, 30, 915-917. http://dx.doi.org/10.1063/1.866345
- Ariel, P.D. (1994) Hiemenz ow in hydromagnetics. Acta Mechanica, 103, 31-43. http://dx.doi.org/10.1007/BF01180216
- Mahapatra, T.R. and Gupta, A.S. (2001) Magnetohydrodynamic Stagnation Point Flow towards a Stretching Sheet. Acta Mechanica, 152, 191-196. http://dx.doi.org/10.1007/BF01176953
- Ishak, A., Nazar, R. and Pop, I. (2008) Magnetohydrodynamics (MHD) Flow of a Micropolar Fluid towards a Stagnation Point on a Vertical Surface. Computers and Mathematics with Applications, 56, 3188-3194. http://dx.doi.org/10.1016/j.camwa.2008.09.013
- Kumar, L., Singh, B., Kumar, L. and Bhargava, R. (2011) Finite Element Solution of MHD Flow of Micropolar Fluid towards a Stagnation Point on a Vertical Stretching Sheet. International Journal of Applied Mathematics and Mechanics, 7, 14-30.
- Ishak A., Nazar R. and Pop, I. (2006) Mixed Convection Boundary Layers in the Stagnation Point Flow toward a Stretching Vertical Sheet. Meccanica, 41, 509-518. http://dx.doi.org/10.1007/s11012-006-0009-4
- Eringen, A.C. (2001) Microcontinuum Field Theories. II: Fluent Media. Springer, New York.
- Ramachandran, N., Chem, T.S. and Armaly, B.F. (1988) Mixed Convection in a Stagnation Flows Adjacent to a Vertical Surfaces. Journal of Heat Transfer, 110, 373-377. http://dx.doi.org/10.1115/1.3250494
- Hassanien, I. and Gorla, R.S.R. (1990) Combined Forced and Free Convection in Stagnation Flows of Micropolar Fluids over Vertical Non-Isothermal Surfaces. International Journal of Engineering Science, 28, 783-792. http://dx.doi.org/10.1016/0020-7225(90)90023-C
- Devi, C.D.S., Takhar, H.S. and Nath, G. (1991) Unsteady Mixed Convection Flow in Stagnation Region Adjacent to a Vertical Surface. Heat and Mass Transfer, 26, 71-79. http://dx.doi.org/10.1007/bf01590239
- Lok, Y.Y., Amin, N., Campean, D. and Pop, I. (2005) Steady Mixed Convection Flow of a Micropolar Fluid Near the Stagnation Point on a Vertical Surface. International Journal of Numerical Methods for Heat & Fluid Flow, 15, 654- 670. http://dx.doi.org/10.1108/09615530510613861
- Lok, Y.Y., Amin, N. and Pop, I. (2006) Unsteady Mixed Convection Flow of a Micropolar Fluid Near the Stagnation Point on a Vertical Surface. International Journal of Thermal Science, 11, 49-57. http://dx.doi.org/10.1016/j.ijthermalsci.2006.01.015
- Mahapatra, T.R. and Gupta, A.S. (2002) Heat Transfer in Stagnation-Point Flow towards a Stretching Sheet. Heat Mass Transfer, 38, 517-523. http://dx.doi.org/10.1007/s002310100215
- Attia, H.A. (2003) Hydrodynamic Stagnation Point Flow with Heat Transfer over a Permeable Surface. Arabian Journal for Science and Engineering, 28IB, 107-112.
- El-Dabe, N.T., Kamel, K.A., Galila abd-Allah, M. and Ramadan, S.F. (2013) Heat Absorption and Chemical Reaction Effects on Peristaltic Motion of Micropolar Uid through a Porous Medium in the Presence of Magnetic Field. The African Journal of Mathematics and Computer Science Research, 6, 94-101.
- Olanrewaju, P.O., Okedayo, G.T. and Gbadeyan, J.A. (2011) Effects of Thermal Radiation on Magnetohydrodynamic (MHD) Flow of a Micropolar Fluid towards a Stagnation Point on a Vertical Plate. International Journal of Applied Science and Technology, 6, 219-230.
- Pop, S.R., Grosan, T. and Pop, I. (2005) Radiation Effect on the Flow near the Stagnation Point of a Stretching Sheet. Technische Mechanik, 25, 100-106.
- Olanrewaju, P.O. and Adesanya, A.O. (2011) Effects of Radiation and Viscous Dissipation on Stagnation Flow of a Micropolar Fluid towards a Vertical Permeable Surface. Australian Journal of Basic and Applied Sciences, 5, 2279- 2289.
- Manjoolatha, E., Bhaskar Reddy, N. and Poornima, T. (2013) Radiation and Mass on MHD Flow of a Micropolar Uid towards a Stagnation Point on a Vertical Stretching Sheet. International Journal of Engineering Research and Applications (IJERA), 3, 195-204.
- Brewster, M.Q. (1992) Thermal Radiative Transfer and Properties. John Wiley & Sons Ltd., New York.
- Zhou, J.K. (1986) Differential Transformation and Its Application for Electrical Circuits. Huazhong University Press, Wuhan.
- Ashraf, M., Kamal, M.A. and Syed, K.S. (2009) Numerical Simulation of a Micropolar Fluid between a Porous Disk and Non-Porous Disk. Applied Mathematical Modelling, 33, 1933-1943. http://dx.doi.org/10.1016/j.apm.2008.05.002

NOTES

^{*}Corresponding author.