Journal of Fluids

Volume 2013 (2013), Article ID 604893, 9 pages

http://dx.doi.org/10.1155/2013/604893

Research Article

Slip-Flow and Heat Transfer in a Porous Microchannel Saturated with Power-Law Fluid

Yazan Taamneh¹ and Reyad Omari²

¹Department of Mechanical Engineering, Tafila Technical University, P.O. Box 179, Tafila 66110, Jordan

²Department of Mathematics, Al-Balqa Applied University, Irbid University College, P.O. Box 19117, Irbid 19110, Jordan

Received 18 May 2013; Accepted 5 October 2013

Academic Editor: Ciprian Iliescu

Copyright © 2013 Yazan Taamneh and Reyad Omari. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This study aims to numerically examine the fluid flow and heat transfer in a porous microchannel saturated with power-law fluid. The governing momentum and energy equations are solved by using the finite difference technique. The present study focuses on the slip flow regime, and the flow in porous media is modeled using the modified Darcy-Brinkman-Forchheimer model for power-law fluids. Parametric studies are conducted to examine the effects of Knudsen number, Darcy number, power law index, and inertia parameter. Results are given in terms of skin friction and Nusselt number. It is found that when the Knudsen number and the power law index decrease, the skin friction on the walls decreases. This effect is reduced slowly while the Darcy number decreases until it reaches the Darcy regime. Consequently, with a very low permeability the effect of power law index vanishes. The numerical results indicated also that when the power law index decreases the fully-developed Nusselt number increases considerably especially, in the limit of high permeability, that is, nonDarcy regime. As far as Darcy regime is concerned the effects of the Knudsen number and the power law index of the fully-developed Nusselt number is very little.

1. Introduction

Fluid flow and heat transfer in porous media has been a subject of continuous interest during past decades because of the wide range of engineering applications. In addition to conventional applications including solar receivers, building thermal insulation materials, packed bed heat exchangers, and energy storage units, new applications in the emerging field of microscale heat transfer have existed. However, microchannels are now used in several industries and equipment such as cooling of electronic package, microchannel heat sinks, microchannel heat exchanger, microchannel fabrication, and cooling, and heating of different devices [1–5].

One of the major difficulties in trying to predict the gaseous transport in micron sized devices can be attributed to the fact that the continuum flow assumption implemented in the Navier-Stokes equations breaks down when the mean free path of the molecules (λ) is comparable to the characteristic dimension of the flow domain. Under these conditions, the momentum and heat transfer start to be affected by the discrete molecular composition of the gas and a variety of noncontinuum or rarefaction effects are likely to be exhibited such as velocity slip and temperature jump at the gas-solid interface. Velocity profiles, fluid flow rate, boundary wall shear stresses, temperature profiles, heat transfer rates, and Nusselt number are all influenced by the noncontinuum regime.

However, there is a certain limit of the channel size with which one can still apply Navier-Stokes equations with some modifications on the boundary conditions [6]. This is the case when Knudsen number (Kn) is in the range , and the flow under such condition is called slip-flow. Knudsen number is defined as the ratio of the molecular mean free path to the characteristic length of the system. It is also used to measure of the degree of rarefaction of gases encountered in flows through narrow passages, and also to measure the degree of the validity of the continuum model.

The continuum model is valid for very small Knudsen number flows (). While the Knudsen number increases, the rarefaction effects become more pronounced, and eventually the continuum assumption breaks down. Therefore, several researchers have suggested that the well-accepted no-slip boundary condition may not be suitable for flows at the micro- and nanoscale [7–10]. Recently, lots of mechanisms have been proposed to explain this phenomenon. In fact, Arkilic et al. [6], Beskok and Karniadaksi [11], and Sparrow and Lin [12], have found that the Navier-Stokes equations, when combined with velocity-slip boundary conditions, yield results for pressure drop and friction factor that are in agreement with experimental data for some microchannel flows.

The appropriate flow and heat transfer models for a given gas flow problem depend on the range of Knudsen number. A classification of the different gas flow regimes is given as follows: for the continuum flow, for the slip flow, for the transition flow, and for the free molecular flow. In this study, the slip flow regime () is considered and the modified extended Darcy-Brinkman-Forchheimer model for power-law fluid is employed to describe the flow behavior in porous medium.

 Convection heat transfer in circular and noncircular microchannels has been solved over the years [13–15]. In these studies, the effects of velocity slip and temperature jump at the wall and viscous dissipation were considered. The main consequence was that the velocity slip and temperature jump have opposite effects on heat transfer. Although the velocity slip tends to increase the Nusselt number, the temperature jump tends to reduce it. The inclusion of the viscous heating increases the Nusselt number for the fluid being cooled and decreases it for the fluid being heated.

Laminar forced convection of Newtonian fluid flow in microchannels filled with a porous medium has been solved by numerical and analytical means over the years [16–19]. In these studies, the effect of Knudsen number, Darcy number, Forchheimer number, and Reynolds number on the velocity slip and temperature jump at the wall were considered. The main results were that the skin friction had been increased by (i) decreasing the Knudsen number, (ii) increasing the Darcy number, and (iii) decreasing the Forchheimer number. Heat transfer was found to be (i) decreased as the Knudsen and Forchheimer numbers increase and (ii) increased as the Reynolds and Darcy numbers increase.

A theoretical and numerical analysis of the fully-developed forced convection in a porous channel saturated with a power-law fluid in porous channel has been investigated recently [20–24]. A closed form boundary layer solutions using the integral method was obtained for velocity profiles, temperature fields, and fully-developed Nusselt number. The theoretical solutions can be used to predict primary characteristics of physical phenomena associated with forced convection of nonNewtonian fluids in porous media. These solutions were convenient to serve as a benchmark for more complicated numerical solutions. The results indicated that the nonDarcy regime, the effects of power law index on hydrodynamics and heat transfer behavior in the porous channel are significant, whereas in the Darcy regime the effects of Darcy number were predominant. Researchers have also found that, in the nonDarcy regime, the Nusselt number increases and the pressure drop decreases as the power law index decreases. Consequently, the combined use of a highly permeable porous matrix with shear thinning fluid appeared to be promising as a heat transfer augmentation technique.

However, little of information on the related literature regarding the flow and heat transfer of power law fluids through porous microchannels. That said, in this study, the forced convection of heat and fluid flow of power law fluids through parallel plate microchannels filled with porous media were considered. The aim of the present study is to investigate the effects of Knudsen number, Darcy number and the inertia parameter on the hydrodynamic and thermal behavior of a power law fluid flow between infinitely long parallel-plates microchannels filled with porous media.

2. Mathematical Formulation

The analysis is carried out for unsteady state, incompressible and laminar forced convection flow between parallel-plates microchannel filled with porous medium and heated with uniform wall temperature at the walls. The flow is assumed to be hydrodynamically fully-developed. The porous medium is saturated with a single phase nonNewtonian fluid described by the power law model and assumed to be in local thermodynamic equilibrium with the fluid.

As a result of the continuity equation, the flow is a unidirectional and is expressed in terms of the axial velocity alone. That is, the velocity component in the y-direction vanishes and the velocity component in the x-direction (denoted by u) becomes dependent on . In addition, the flow is assumed to be thermally developing under constant pressure gradient driving force; therefore, the temperature becomes the function of (x, y) only. The present study focuses on the slip flow regime (), and therefore the Navier-Stokes equations and energy equation combined with slip/jump boundary conditions have been applied. The physical properties of the solid matrix and of the fluid are assumed to be constant except for the viscosity of the power-law fluid which depends on the shear rate. In the present study, a fibrous or foam-metal material is considered such that the porosity and permeability are assumed to be constant even close to the walls. On the other hand, it is assumed that the porous medium is isotropic and homogeneous. Finally, viscous dissipation is neglected in the energy equation. The governing equations can be written as follows [22].

Continuity is

Momentum is

Energy is

In the momentum equation, a modified Darcy’s law for power-law fluids was used where represents the consistency of the power-law fluid and is the modified permeability which depends on the structure of the porous medium and on the power law index of the fluid [23]. The linear approximation of Darcy’s law can directly be derived from the fluid macroscopic, Navier-Stokes, and momentums balance equation [25].

In the energy equation, the thermal dispersion conductivity of the porous media is assumed to be constant and is incorporated into the effective thermal conductivity. The axial heat conduction effects are usually negligible for nearly parallel flows. Momentum and energy transfer between the liquid molecules and the surface requires specification of interactions between the impinging molecules and the surface. From the macroscopic point of view, it is sufficient to know some average parameters in terms of the so-called tangential momentum () and thermal () accommodation coefficients. These coefficients describe the gas-surface interaction and are functions of the composition and temperature of the gas, the gas velocity over the surface, and the solid surface temperature, chemical state, and roughness. The accommodation coefficients take any value between 0 and 1, where these values represent specular reflection and diffuse reflection, respectively. For most engineering applications, values of the accommodation coefficients are near unity [9].

Under the above assumptions and by using the nondimensional variables listed in the nomenclature, the equations of motion and energy equation are reduced to the following form:

where is the inertia parameter and is the modified inertia parameter. Also, the boundary conditions, in nondimensional form, are as follows:

The quantities of primary interest in this study are the friction factor and Nusselt number. These are defined as follows:

where is the dimensionless mean temperature, is the hydraulic diameter, is the local heat transfer coefficient and is thermal conductivity of the fluid.

3. Numerical Method

The governing equations are solved numerically using the finite difference technique. The governing momentum and energy equations are not coupled, consequently the numerical solution proceeds by first solving the velocity distribution from the momentum equation, and then solving the energy equation for the temperature distribution.

The momentum equation (2) is parabolic partial differential equation (unsteady, one-dimensional diffusion equation) if the time term is left in the equation. One way to solve this type of equations is by using the well-known Backward-Time Central-Space (BTCS) method [26]. This particular method is also called the fully implicit method. The finite difference equation (FDE) which approximates the partial differential equation is obtained by replacing the exact partial derivative by the backward-time approximation, while the exact partial derivative is replaced by the centered-space approximation. The steady state solution is obtained by marching in time until no further significant change in the solution is obtained with additional marching steps. The present numerical method is advantageous over other available numerical methods in that it is unconditionally stable.

In a similar manner, the energy equation is discretized using the same numerical scheme. In contrast to the momentum equation, it should be mentioned that the steady state form of the energy equation is discretized directly and therefore a time dependent solution was not obtained. When the momentum equation is discretized and applied at every point in the finite difference grid, a system of nonlinear finite difference (algebraic) equations was obtained. To overcome the difficulties in solving such a system, the nonlinear term (last term) in (2) should be linearized. To do so, the well-known time lagging technique is used [26].

Based on the above approach, the resulting systems of algebraic equations obtained by discretizing the momentum and energy equations are tri-diagonal, which are best solved by using Thomas algorithm. The adequacy of the grid is verified by comparing the results of different grid sizes. A mesh refinement study was carried out in order to ensure grid independent solutions. It was found that the obtained numerical solution for the momentum equation is invariant beyond a grid size of 75 points in the y-direction. Therefore, all velocity profiles are obtained using this grid size. Similar refinement study was carried out for the energy equation. It was found that a grid size of is adequate for the accuracy and any increase in the number of grid points would result in an insignificant effect on the results.

4. Result and Discussion

In order to verify the validity and accuracy of the numerical model, the present numerical results were compared with corresponding integral solution results for the case of fully developed forced convection in porous macrochannel saturated with a power-law fluid [23]. Figure 1 shows a comparison between the two solutions where the agreement is very good at the Darcy and nonDarcy regimes.

Figure 1: Comparison of results for fully developed Nusselt number with analytical results of Chen and Hadim. .

The effect of the Knudsen number and the power law index on the axial fully-developed velocity profiles is shown in Figure 2 for a microscopic inertial coefficient equal to 0. From this figure, it is clear that as the number increases, the velocity slip at the wall increases regardless of the power law index value. This is because the increase in Kn number can be due to increase in the mean free path of the molecules, which, in turn, decreases the retarding effect at the wall and thus yields larger flow rates near the channel walls. It is obvious from Figure 2 that as the Kn numberincreases or even when the power law index decreases, the flow velocity near the walls increases while the peak velocity at the centerline decreases to satisfy mass conservation. Figure 2 also shows that the effects of Kn number is more evident for shear thinning fluids () and higher values of Kn.

(a)

(b)

Figure 2: Effect of Nudsen number and power law indexon the axial velocity profile for .0, .0. (a) , (b) .

The effect of the Darcy number and the power law index on the axial fully developed velocity profiles is shown in Figure 3 for . As the Darcy number decreases, a flat velocity profile occurs in the core region due to bulk damping caused by the presence of the porous matrix and the viscous effects near the walls. In Figure 3, it is shown that the effects of Darcy number are more considerable for shear thickening fluid () than for shear thinning fluids (). Thoroughly inspecting Figure 3, it is obvious that as the power law index decreases, the velocity gradient near the wall increases while the peak velocity at the centerline decreases to satisfy mass conservation. It is noted that the effect of power law index become insignificant at low Darcy number ().

(a)

(b)

Figure 3: Effect of Darcy number and power law index on the axial velocity profile for .0, .001. (a) , (b) .

In contrast, Figure 4 shows the effect of the inertia parameter and the power law index on the axial fully developed velocity profile for . As the inertia parameter increases, the fully developed velocity profile in channel becomes flattened (velocity gradient near the walls increases), and the boundary layer behavior is more pronounced [22]. However, the flattening of the velocity profile is relatively weaker for shear thinning fluids compared with shear thickening fluids.

(a)

(b)

Figure 4: Effects of inertia parameter and power law index on the axial velocity profile for .0, .001. (a) , (b) .

The combined effects of the Knudsen number, power law index and Darcy number on the skin friction, are clearly presented in Figure 5 for .0. As the Kn number increases, the skin friction decreases. It is known that any increase in Kn number would increase in number due to the increase in the flow velocity, which is a result of the reduction in the retardation effect of the wall. On the other hand, the skin friction value decreases due to both the decrease in the velocity gradient at the wall and the increase in the flow velocity. The net result of the effect of Kn number is found to decrease value. This means that the reduction in skin friction value is more significant so that it overcomes the increase in Re number. The effect of the Knudsen number is more significant at larger values of Darcy number. In Figure 5 also the effects of the Kn number diminish while the Darcy number decreases until they become negligible in the Darcy regime (i.e., very low permeability) due to the obvious effects of the porous matrix. As the Darcy number increases, the skin friction decreases sharply, especially in the Darcy regime, thus approaching the asymptotic value for the flow inside the clear channel. As the power law index decreases, the skin friction also decreases significantly. In the nonDarcy regime, it is noted that the skin friction for shear thickening fluids is about five times the skin friction for shear thinning fluids. In fact, this conclusion supports the numerical results obtained by [22].

Figure 5:  Effects of Knudsen number and power law index on the skin friction with Darcy number for .0.

The combined effects of the Kn number, power law index, and the inertia parameter on the skin friction are clearly presented in Figure 6 for . At a fixed value of Kn number. The figure shows that for small values of the inertia parameter, the skin friction remains constant and then increases at a faster rate. When the inertia parameter increases, the velocity gradient near the wall increases therefore increasing the skin friction near the walls. Also, the effect of the Kn value on the skin friction becomes more definite at higher values of inertia parameter because at higher values of inertia parameter, the velocity gradient near the wall gets relatively higher. In the case of very low permeability, as the power law index increases, the skin friction increases.

(a)

(b)

Figure 6: Effect of Knudsen number and power law index on the skin friction with inertia parameter for .0. (a) (b) .

The effect of the Kn number and power law index on the variation of the fully developed Nusselt number with inertia parameter is shown in Figure 7. Figure 7 shows that at relatively lower Kn number (), the inertia effects are very definite while Kn number increases (), these effects become very little. Moreover, at relatively small values of inertia parameter, the fully developed Nusselt number remains constant and depends on the power law index and Knudsen number. As the inertia parameter increases, the fully developed Nusselt number increases and then approaches an asymptotic value. It is clear from Figure 7 that for lower Knudsen number case (), shear thinning fluids results in higher Nusselt number due to the large velocity gradients near the walls. Also, the effects of the power law index on the fully-developed Nusselt number becomes more evident at lower values of inertia parameter because, at lower values of inertia parameter, the velocity near the walls are relatively lower leading to lower convection heat transfer rate. This implies that any change in the power law index at lower inertia parameter would lead to a significant change in the fully-developed Nusselt number, especially at lower Kn number.

Figure 7:   Effects of Knudsen number and power law index on the variation of the fully developed Nusselt number with inertia parameter.

Figure 8 shows the combined effects of the Kn number, power law index, and Darcy number on the fully-developed Nusselt number when . It is noted that the effects of Kn number diminishes as the Darcy number decreases until they become negligible in the Darcy regime (i.e., very low permeability) because of the obvious effects of the porous matrix. When the Knvalue increases, the fully developed Nusselt number decrease regardless of the power law index. This is because when normally the Kn number increases, the jump in temperature at the wall increases and a lesser amount of energy is consequently transferred from the wall to the adjacent fluid. As the Darcy number increases, the fully developed Nusslet number for shear thickening () decreases more steeply than at shear thinning () due to the velocity profile attained at the different power law index. When Darcy number increases, the asymptotic regime is identified as the fully developed Nusselt number which approaches the asymptotic value for the flow in clear channel case.

(a)

(b)

Figure 8: Effects of Knudsen number and power law index on the variation of the fully developed Nusselt number with Darcy number (a) (b) .

5. Conclusion

The present numerical solutions are conducted for steady laminar forced convection flow between parallel-plate microchannels filled with porous medium and saturated with a power-law fluid. In this study, the slip flow regime () is considered, and the modified extended Darcy-Brinkman-Forchheimer model is adopted to describe the hydrodynamic and thermal behavior of the power-law fluid in porous microchannels. The present study reports the effects of power law index, Knudsen number, Darcy number, microscopic inertial coefficient on the flow, and heat transfer in the microchannel. The results indicate that in the case of high permeability regime (nonDarcy), the effects of the Knudsen number and the power law index on the flow and heat transfer in the porous microchannel are significant. However, in case of low permeability regime, the effects of Darcy number become more evident. In case of high permeability, the skin friction and the Nusselt number decreases while the Knudsen number and the power law index increase.

Nomenclature

: | Specific heat, J/kgK
: | Inertia coefficient
: | Modified inertia coefficient
: | Darcy number,
: | Hydraulic diameter,
: | Channel height, m
: | Local heat transfer coefficient
: | Intrinsic permeability of the porous medium,
: | Modified permeability of the porous medium,
Kn: | Knudsen number, ()
: | Thermal conductivity of the fluid, W/m-K
: | Effective thermal conductivity of the fluid saturated porous medium, W/m-K
: | Power law index, m
: | Local Nusselt number, ()
: | Fully developed Nusselt number
: | Pressure, Pa
: | Dimensionless pressure,
: | Prandtl number,
: | Reynolds number,
: | Temperature, K
: | Time, s
: | Reference time, (/)
: | Nondimensional axial velocity, ()
: | Axial velocity component, m/s
: | Dimensionless axial coordinate, (/)
: | Axial coordinate, m
: | Dimensionless transverse coordinate, (/)
: | Transverse coordinate, m.

Greek Symbols

: | Thermal diffusivity ()
: | Mean free path of the gas molecules
: | Tangential momentum accommodation coefficient
: | Thermal accommodation coefficient
: | Specific heat ratio, ()
: | Microinertia parameter, ()
: | Porosity of the porous media
: | Dynamic viscosity
: | Consistency index of the power law fluid
: | Fluid density
: | Nondimensional temperature,
: | Dimensionless mean temperature,
: | Dimensionless time,
: | Shear stress at the wall (.

Subscripts

: | Effective (i.e., fluid-saturated porous medium)
: | Fluid
: | Mean
: | Pressure
: | Solid
: | Volume
: | Wall.

References

G. Karniadakis, A. Beskok, and N. Aluru, Micro Flows and Nano Flows: Fundamentals and Simulation, Springer, New York, NY, USA, 2005.

G. Karniadakis and A. Beskok, Micro Flows, Springer, New York, NY, USA, 2002.

D. Nield and A. Bejan, Convection in Porous Media, Springer, New York, NY, USA, 1999.

K. Vafai and C. L. Tien, “Boundary and inertia effects on flow and heat transfer in porous media,” International Journal of Heat and Mass Transfer, vol. 24, no. 2, pp. 195–203, 1981.

K. Watanabe, Y. Yanuar, and H. Mizunuma, “Slip of Newtonian fluids at slid boundary,” JSME International Journal B, vol. 41, no. 3, pp. 525–529, 1998.

E. B. Arkilic, K. S. Breuer, and M. A. Schmidt, “Gaseous flow in microchannels,” in Proceedings of the Applications of Micro-Fabrication to Fluid Mechanics (ASME '94), vol. 197, pp. 57–66, November 1994.

F. Ezquerra Larrodé, C. Housiadas, and Y. Drossinos, “Slip-flow heat transfer in circular tubes,” International Journal of Heat and Mass Transfer, vol. 43, no. 15, pp. 2669–2680, 2000.

R. F. Barron, X. Wang, T. A. Ameel, and R. O. Warrington, “The Graetz problem extended to slip-flow,” International Journal of Heat and Mass Transfer, vol. 40, no. 8, pp. 1817–1823, 1997.

S. Yu and T. A. Ameel, “Slip-flow heat transfer in rectangular microchannels,” International Journal of Heat and Mass Transfer, vol. 44, no. 22, pp. 4225–4234, 2001.

J. Liu, Y.-C. Tai, and C.-M. Ho, “MEMS for pressure distribution studies of gaseous flows in microchannels,” in Proceedings of the IEEE Micro Electro Mechanical Systems Conference, pp. 209–215, February 1995.

A. Beskok and G. E. Karniadakis, “Simulation of heat and momentum transfer in complex microgeometries,” Journal of Thermophysics and Heat Transfer, vol. 8, no. 4, pp. 647–655, 1994.

E. M. Sparrow and S. H. Lin, “Laminar heat transfer in tubes under slip-flow conditions,” Journal of Heat Transfer, vol. 84, pp. 363–369, 1962.

T. A. Ameel, R. F. Barron, X. Wang, and R. O. Warrington Jr., “Laminar forced convection in a circular tube with constant heat flux and slip flow,” Microscale Thermophysical Engineering, vol. 1, no. 4, pp. 303–320, 1997.

W. Qu, G. M. Mala, and D. Li, “Heat transfer for water flow in trapezoidal silicon microchannels,” International Journal of Heat and Mass Transfer, vol. 43, no. 21, pp. 3925–3936, 2000.

Y. Zhu and S. Granick, “Rate-dependent slip of Newtonian liquid at smooth surfaces,” Physical Review Letters, vol. 87, no. 9, Article ID 096105, pp. 961051–961054, 2001.

O. M. Haddad, M. A. Al-Nimr, and Y. Taamneh, “Hydrodynamic and thermal behavior of gas flow in microchannels filled with porous media,” Journal of Porous Media, vol. 9, no. 5, pp. 403–414, 2006.

O. M. Haddad, M. A. Al-Nimr, and M. S. Sari, “Forced convection gaseous slip flow in circular porous micro-channels,” Transport in Porous Media, vol. 70, no. 2, pp. 167–179, 2007.

C.-H. Chen, “Thermal transport characteristics of mixed pressure and electro-osmotically driven flow in micro- and nanochannels with Joule heating,” Journal of Heat Transfer, vol. 131, no. 2, pp. 1–10, 2009.

M. Kaviany, “Laminar flow through a porous channel bounded by isothermal parallel plates,” International Journal of Heat and Mass Transfer, vol. 28, no. 4, pp. 851–858, 1985.

A. Nakayama and A. V. Shenoy, “Non-Darcy forced convective heat transfer in a channel embedded in a non-Newtonian inelastic fluid-saturated porous medium,” Canadian Journal of Chemical Engineering, vol. 71, no. 1, pp. 168–173, 1993.

A. V. Shenoy, “Non-Newtonian fluid heat transfer in porous media,” Advances in Heat Transfer, vol. 24, pp. 101–190, 1994.

G. Chen and H. A. Hadim, “Forced convection of a power-law fluid in a porous channel—numerical solutions,” Heat and Mass Transfer, vol. 34, no. 2-3, pp. 221–228, 1998.

G. Chen and A. Hadim, “Forced convection of power-law fluid in porous channel-Integral solution,” Journal of Porous Media, vol. 2, pp. 59–70, 1999.

A. Hadim, “Forced convection in a porous channel with localized heat sources,” ASME Journal of Heat Transfer, vol. 116, no. 2, pp. 465–472, 1994.

S. Sorek, D. Levi-Hevroni, A. Levy, and G. Ben-Dor, “Extensions to the macroscopic Navier-Stokes equation,” Transport in Porous Media, vol. 61, no. 2, pp. 215–233, 2005.

J. D. Hoffman, Numerical Methods for Engineers and Scientists, McGraw Hill, New York, NY, USA, 1992.