Journal of Engineering
Volume 2016 (2016), Article ID 9718786, 10 pages
http://dx.doi.org/10.1155/2016/9718786

Research Article

Similarity Solution for High Weissenberg Number Flow of Upper-Convected Maxwell Fluid on a Linearly Stretching Sheet

Meysam Mohamadali and Nariman Ashrafi

Young Researchers and Elites Club, Science and Research Branch, Islamic Azad University, Tehran, Iran

Received 19 November 2015; Revised 9 March 2016; Accepted 24 April 2016

Academic Editor: Oronzio Manca

Copyright © 2016 Meysam Mohamadali and Nariman Ashrafi. 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.

High Weissenberg boundary layer flow of viscoelastic fluids on a stretching surface has been studied. The flow is considered to be steady, low inertial, and two-dimensional. Upon proper scaling and by means of an exact similarity transformation, the nonlinear momentum and constitutive equations of each layer transform into the respective system of highly nonlinear and coupled ordinary differential equations. Numerical solutions to the resulting boundary value problem are obtained using an efficient shooting technique in conjunction with a variable stepping method for different values of pressure gradients. It is observed that, unlike the Newtonian flows, in order to maintain a potential flow, normal stresses must inevitably develop. The velocity field and stresses distributions over plate are presented for difference values of pressure gradient and Weissenberg numbers.

1. Introduction

The flow of a liquid within a thin film over stretching plate is often encountered in most manufacturing processes. Examples include extrusion of plastic sheets, fabrication of adhesive tapes, and application of coating layers onto rigid substrates. Coating processes demand a smooth glossy surface to meet the requirements for best appearance and optimum service properties such as low friction, transparency, and strength. Due to the moving surface, the main flow is closed to the extruded material while the far field stays almost stagnant.

Most structured liquids, like polymers, show strong viscoelastic effects at small deformations, and their measurement is very useful as a physical probe of the microstructure. The subject of the viscoelastic flow was applied to developing pipeline designs such that the fluidity of high wax content crude oils is maintained [1]. More applications specific to the viscoelastic model include the modeling of plastics such as PET resins, which are used in, for example, film casting [2]. It has also been used to simulate the flow of polymer which is used for the process of wire-coating [3].

The analytical study of boundary layer flow due to a stretching sheet was initiated by Crane [4]. He assumed the velocity of the sheet to vary linearly as the distance from the slit and obtained an analytical solution. The work of Crane was subsequently extended mostly on both Newtonian and non-Newtonian (inelastic) boundary layer flows (see, e.g., [5–8]) and only a few works on viscoelastic and elastic boundary layer flows [9]. There is no report of such work on the stretching sheet. In this connection, however, Hassanien [10] studied the second-grade fluid boundary layer over a linearly stretching sheet. The study was on boundary layer approximations of Newtonian flows [11] in order to simplify the governing equations. Here, boundary layer equations were solved by a similarity method for elastic flows of Deborah numbers of up to 0.2.

The upper-convected Maxwell fluid is a class of viscoelastic fluid that can explain characteristics of fluid relaxation time. It excludes complicated effects of shear-dependent viscosity and thus allows one to emphasize the influence of fluid’s elasticity on characteristics of its boundary layer. So far, the exact solution corresponding to the unsteady flow of a Maxwell fluid induced by the impulsive motion of a plate between two side walls perpendicular to the plate is developed employing the Fourier sine transforms [12]. It is reported, in the mentioned research, that the velocity decreases by increasing the relaxation time while the magnitude of shear stress increases. Furthermore, Shateyi [13] studied the MHD flow of UCM past a vertical stretching sheet in a Darcian porous medium under the influence of thermophoresis, thermal radiation, and a uniform chemical reaction for Deborah as high as unity.

Moreover, the unsteady flow of Maxwell fluid induced over oscillating accelerated sheet was investigated in [14]. In another work, Ashraf et al. [15] used Homotopy method to simulate the flow and heat transfer of UCM fluid over a moving surface in a parallel free stream with the convective surface boundary condition. Recently time-dependent three-dimensional boundary layer flow of a Maxwell fluid over a stretched sheet has been investigated by Homotopy method [16].

In the above viscoelastic flows the governing equations are scaled by Reynolds number only (similar to Newtonian fluids). Effect of viscoelasticity is, therefore, not realistically presented. It is already reported that even at low Reynolds and high Weissenberg numbers a boundary layer develops in the flow of viscoelastic fluid [17]. This fact is also backed by experiments [18]. High Weissenberg flows mean long relaxation time in which the velocity of fluid vanishes at the wall and particles away from the wall travel long distances within one relaxation time so that particles close to the wall travel only a short distance. This leads to boundary layer in the shear stress [19]. The viscoelastic boundary layer is formed in a thin region closer to the wall in which the relaxation terms are recovered.

Up to now the boundary layer equations for the UCM fluids in two-dimensional creeping flow along a flat boundary for high Weissenberg numbers are derived [19, 20]. It was shown that scaling parameters in view of the high Weissenberg condition and taking the leading terms of the upper-convected Maxwell fluid governing equations result in the viscoelastic boundary layer development of order Wi ⁻¹. Similar studies on the Phan-Thien-Tanner (PTT) and the Giesekus fluids result in the boundary layer development of order Wi −1/3 and order Wi −1/2, respectively [21]. This phenomenon can also be physically interpreted that “elastic” boundary layers for the Phan-Thien-Tanner and Giesekus fluid are similar to those for the upper-convected Maxwell model and arise when the dimensionless parameter measuring the size of the quadratic term is small. In fact, if the quadratic term is not small, the PTT model will have “viscometric” boundary layers in which it behaves like a generalized Newtonian fluid [22]. For the Giesekus model, the viscometric behavior is different in that the shear stress remains bounded at infinite shear rate. Furthermore, Lie group theory is used in [9] by Atalık to obtain point symmetries of the boundary layer equations derived for the high Weissenberg number flow of UCM and Phan-Thien-Tanner (PTT) fluids.

Using an implicit function, the existence of solutions for viscoelastic boundary layer which arises from spatially periodic perturbations of uniform shear flow was addressed [23]. Also, the well-posedness of boundary layer equations for time-dependent flow of UCM fluid in the limit of high Weissenberg and Reynolds numbers was analyzed [24]. Furthermore, a systematic perturbation procedure to solve the initial value problem for creeping flow of the UCM fluid at high Weissenberg number is formulated [25].

For instance, citing an analogy between a viscoelastic medium and an electrically conducting fluid containing a magnetic field, Ogilvie and Proctor [26] showed that the dynamics of the Oldroyd-B fluid in the limit of large Deborah number correspond to that of a magnetohydrodynamic (MHD) fluid in the limit of large magnetic Reynolds number. In some aspects, the problem of high Weissenberg number asymptotic for viscoelastic flows is similar to high Reynolds number asymptotic for Newtonian flows. The boundary layer also arises in high Weissenberg number flows since the convected derivative terms become essential at a short distance from the wall, leading to the formation of the aforementioned sharp boundary layer in the stresses [19].

The aim of this work is study of the boundary layer formation in high Weissenberg creeping flow of UCM fluids past a stretching plate using similarity transformation. The stretching rate is assumed to be proportional to the ratio of horizontal distance on the direction. Using similarity transformation the partial differential governing equations are transformed into a set of ordinary differential equations. The ordinary differential equations are then integrated numerically using a Runge-Kutta subroutine and shooting technique. Typical results for the velocity and stress profiles are presented.

2. Governing Equations

The steady flow of a viscoelastic fluid over a (linearly) stretching sheet () is brought to attention here. Consider two-dimensional steady flow of an upper-convected Maxwell fluid occupying the half-plane . The fluid is flown by the movement of a thin elastic sheet emerging from a narrow slit at the origin of a Cartesian coordinate system under investigation shown schematically in Figure 1. The continuity and momentum equations for creeping flow are written as

where is the extrastress tensor, is the velocity field, and and are the pressure and density, respectively.

Figure 1: Schematic diagram showing the physical configuration and coordinate system.

For UCM fluid, the stress tensor, , can be related to the deformation-rate tensor as [27]

where and are the relaxation time and viscosity, respectively, and the symbol stands for the upper-convected derivative. Here, the Cartesian axes are assigned to the flow with the -axis being along the plate and the -axis normal to it (Figure 1). The governing equations may be rewritten in dimensionless form by introducing typical scales for length, velocity, stress, and pressure as follows:

where the capitals and primes represent reference and dimensionless values, respectively. The reference velocity of stretching sheet is assumed to be .

The dimensionless form of the governing equations is obtained by substituting the dimensionless parameter equation (3) into governing equations (1)–(2b):

where is the Weissenberg number characterizing the elastic effects [27]. The no-slip condition on the sheet and the far field condition boundary in dimensionless form are

The basic reason for the formation of viscoelastic boundary layer is quite simple to elaborate. The convected derivatives in the constitutive relation vanish at the wall, forcing the stresses to be viscometric. However, at high Weissenberg number, the convected derivative terms become important at a short distance from the wall, leading to the formation of a boundary layer in the stresses. To maintain the balance between inertial force, viscous and elastic stresses to the leading order, a self-consistent set of scalings of variables is proposed [19]:

where , , , , , and are to be determined. Implementing the scalings, the continuity equation becomes

Balancing the two velocity gradients above requires that . In light of the relation (7) and using the proposed scalings (6), the constitutive relations (4c) take the following form:

In the boundary layer approximation, taking order of velocity to be leads all terms of left hand sides of constitutive equations to remain, so that . This means that fewer terms are ignored and the approximate forms of boundary layer are much closer to the original form of boundary equations. It is observed, from the constitutive equations (8a), (8b), and (8c), that if is small of order , then the term can be comparable to the term and cannot be ignored. This observation suggests the following scalings for the viscoelastic boundary layer:

Substituting the above scalings into (4a), (4b), and (4c) and keeping only the leading order terms in the rescaled equations, the continuity and momentum equations in the boundary layer are obtained:

Similarly the constitutive equations changed to

3. Similarity Solution

In order to solve the governing equation subject to the viscoelastic boundary layer, several similarity transformations were tried among which the following set appears to be computationally useful:

where is similarity variable and , and are unknown similarity functions. Next using these translated stream functions in relation (13) the velocity components are given by

where the prime is derivation with respect to the similarity variable .

Clearly, velocity components and given in (14) satisfy the continuity equation (10). In terms of the new variables, the transformed momentum and stress equations, that is, (11a), (12a), (12b), and (12c), can be rewritten as

where the integration constant can be interpreted as a pressure gradient and the prime indicates derivation with respect to the similarity variable . The no-slip boundary conditions on the surface and far field velocity require that

The no-slip boundary conditions are imposed at the wall; however, due to the singularity of transformed equation on the wall at , it is not possible to integrate right from the wall [29, 30]. Expanding the stream function at a minute distance above the plate (meaning ) returns the following approximation for :

In the same manner, the other similarity functions are expanded as follows:

Substituting the expanded functions (17) and (18) for a very close distance above the wall of nearly into (15a), (15b), (15c), and (15d) the constants are evaluated as

It appears, in the above expressions, that the coefficients of the similarity functions contain independent parameters, and . In what follows, however, a relation between these two parameters can be established. The parameter is pertinent to the wall shear stress as it determines on the wall, while the parameter may be given a physical interpretation by substituting (17) and (18) in transformed momentum equation (15a) for in near-wall region, which gives

This shows that (0) is related to the coefficient of the pressure gradient C [29].

4. Numerical Method

The nonlinear differential equations (15a), (15b), (15c), and (15d) together with the boundary conditions (17)–(20) constitute a boundary value problem (BVP) which is solved numerically by the shooting technique (see, e.g., [31]). In doing so, inevitably the following function changes are introduced.

To numerically solve the above equations the infinite value of must be replaced by a finite value of (say ) which is sufficiently large to satisfy the asymptotic condition. Here the value of 1000 turned out to be appropriate. Consequently equations of (15a), (15b), (15c), and (15d) are integrated numerically by fourth-order Runge-Kutta scheme from to with (17)–(20) and guessed trail values (0) which should satisfy the right-end boundary condition (16). The Newton-Raphson scheme is employed to correct the arbitrary guess value such that the numerical solution will eventually satisfy the required boundary conditions to a precision of . For further details on the numerical procedure, the reader is referred to [31]. As the first set of results, the values of for various pressure gradients are tabulated in Table 1.

Table 1: Values of (0) for various pressure gradients.

(0)

0
−2.0215

−0.5
−2.110

−1
−2.140

−1.5
−2.156

−2
−2.162

5. Results and Discussion

In this section the boundary layer formed by the flow of viscoelastic fluid over stretching sheet is presented. The flow is at high Weissenberg number and low Reynolds number. Within the boundary layer, therefore, the flow variables such as velocity and shear and normal stresses are evaluated for several values of pressure gradients.

Figure 2(a) shows variation of velocity profiles with for different pressure gradients, . It can be seen that the velocity profiles decrease continuously to zero with the increase of the similarity parameter meaning that particles of fluid come to rest far from stretching surface. Additionally, it is clear that increasing the values of results in increasing the magnitude of velocity function causing velocity boundary layer to thicken. Physically it can be said that the increasing in pressure gradient leads to effect of surface velocity that is more sharply sensed far from the sheet.

(a)

(b)

(c)

(d)

Figure 2: Velocity and stress profiles for several values of pressure gradient constant, C.

Figure 2(b) reveals that, regardless of the pressure gradient, the normal stress function in the flow direction, , monotonically reaches a constant value. In another fashion, the normal stress component perpendicular to the flow, , grows by increasing as shown in Figure 2(c). Additionally, the values of both normal stress functions increase with increase of pressure gradient. Also, the value of shear stress, , decreases monotonically to zero as represented in Figure 2(d), implying that the stresses is nonexistent shear force far from the plate. Finally, it is observed that the wall shear stress increases in line with increase of pressure gradient.

To assess the accuracy of the present numerical subroutine, velocity and profiles for zero-pressure gradient () are compared to those obtained by Atalık [9]. It is found from Figure 3 that the velocity profile of current study decreases continuously to zero with the increase of the similarity parameter while Atalık’s result goes to a nonzero constant value which means nonzero velocity in far field. It seems that the current study could aim the outer far field condition better than the other one [9]. Moreover, Table 2 presents results of this comparison for wall shear stress and the value of . It can be observed from Table 2 that there are some discrepancies between the present results and those of Atalık which might lead to the above-mentioned phenomena in the far field flow.

Table 2: Comparison of (0) and for C = 0.

Present
Atalık [9]

(0)
−2.0215
−2

10⁻⁶
10⁻⁴

Figure 3: Comparison of similarity functions of velocity for of current study and Atalık [9].

In continuation, Figure 4 shows scaled stress components, , near the stretching wall in plane. The figures are arranged in a way that the pressure is increased incrementally in each row. Generally, the stress contours of normal stress and shear stress are denser near the wall. On the contrary, an opposite behavior is observed for second normal stress, . This phenomenon is maybe explained that, in high Weissenberg flows, far from the solid boundary, the elastic property of the flow is dominant which leads to formation of normal stress in the absence of velocity gradient and elimination of the shear stress. In this region the normal stress in the flow direction, , goes to a constant value to form a constant stress region far from the sheet. On the other hand the other normal stress component increases by distancing far from the plate. It is also observed that, far from the field, the shear stress becomes negligible. This fact could be attributed to the decreasing velocity gradient at a considerable distance from the surface. A point to make here is that, adjacent to the wall as the pressures gradient increases, the values of all stress components rise.

Figure 4: Scaled stress components in plane for values of .

In continuation, Figures 5 and 6 show dimensionless stress components, , near the solid wall in plane for and 100 and pressure gradients equal to 0 and −0.5, respectively. Generally, the stress contours are denser near the wall forming the viscoelastic boundary layer. This viscoelastic boundary layer becomes thinner with increase of number. This result shows that increasing elastic effects in the flow leads to the thinner boundary layer and the more variation of stresses in near the wall.

Figure 5: Dimensionless stress components on a stretching plate in the plane: for for 10 and 100.

Figure 6: Dimensionless stress components on a stretching plate in the plane: for for 10 and 100.

Furthermore to provide a better physical sense of the overall problem, the streamlines and velocity vector for a Newtonian flow derived from [28] are compared to those of viscoelastic flow equal to 10 in the absence of pressure gradient taken from current study. This comparison is shown in Figures 7 and 8. Conspicuously, the overall trend of the two flows is similar providing a sound basis for agreement between the current simulation and the physical evidence. As seen in Figure 7, far from the plate the velocity is essentially vertical. This may be pertinent to the normal stresses developed within the viscoelastic due to high elastic nature of high Weissenberg flow.

(a)

(b)

Figure 7: Development of (a) stream function and (b) velocity vector field for viscoelastic boundary layer ().

(a)

(b)

Figure 8: Development of (a) stream function and (b) velocity vector field for Newtonian fluid boundary layer () [28].

6. Conclusions

Stress boundary layer of highly elastic “creeping” flow of the upper-convected Maxwell fluid over a linear stretching sheet is studied numerically. Using similarity transformations, the scaled dimensionless momentum and constitutive equations are converted into a system of ordinary differential equations for which a numerical solution is sought. A custom-made shooting technique is devised to ensure arriving at credible solution to the derived highly coupled and highly nonlinear equations.

It can be concluded that, in general, the velocity boundary layer thickness increases with pressure gradient. As well, the values of both normal stress functions increase with increase of pressure gradient. However, the normal stress function in the flow direction, , monotonically reaches a constant value. In another fashion, the normal stress component perpendicular to the flow, , grows by increasing the parameter . Also, the value of shear stress, , decreases monotonically to zero, implying that the shear stress is nonexistent force far from the plate.

Additionally, in high Weissenberg flows, far from the solid boundary, the elastic property of the flow is dominant. This fact results in formation of normal stress in the absence of velocity gradient and elimination of the shear stress so that a “potential” flow is observed while contrary to the Newtonian flow the fluid velocity is essentially vertical to the plate. In this region the normal stress in the flow direction goes to a constant value and the other normal stress component increases by distancing from the plate. Finally, it is observed that the value of stress components increases by Weissenberg number.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

L. T. Wardhaugh and D. V. Boger, “Flow characteristics of waxy crude oils: application to pipeline design,” AIChE Journal, vol. 37, no. 6, pp. 871–885, 1991.

S. G. Hatzikiriakos, G. Heffner, D. Vlassopoulos, and K. Christodoulou, “Rheological characterization of polyethylene terephthalate resins using a multimode Phan-Tien-Tanner constitutive relation,” Rheologica Acta, vol. 36, no. 5, pp. 568–578, 1997.

V. Ngamaramvaranggul and M. F. Webster, “Simulation of pressure-tooling wire-coating flow with Phan-Thien/Tanner models,” International Journal for Numerical Methods in Fluids, vol. 38, no. 7, pp. 677–710, 2002.

L. J. Crane, “Flow past a stretching plate,” Zeitschrift für Angewandte Mathematik und Physik ZAMP, vol. 21, no. 4, pp. 645–647, 1970.

K. V. Prasad, S. R. Santhi, and P. S. Datti, “Non-Newtonian power-law fluid flow and heat transfer over a non-linearly stretching surface,” Applied Mathematics, vol. 3, no. 5, pp. 425–435, 2012.

H. Xu and S.-J. Liao, “Laminar flow and heat transfer in the boundary-layer of non-Newtonian fluids over a stretching flat sheet,” Computers & Mathematics with Applications, vol. 57, no. 9, pp. 1425–1431, 2009.

M. S. Abel, P. S. Datti, and N. Mahesha, “Flow and heat transfer in a power-law fluid over a stretching sheet with variable thermal conductivity and non-uniform heat source,” International Journal of Heat and Mass Transfer, vol. 52, no. 11-12, pp. 2902–2913, 2009.

C. Wang, “Analytic solutions for a liquid film on an unsteady stretching surface,” Heat and Mass Transfer, vol. 42, no. 8, pp. 759–766, 2006.

K. Atalık, “Group theoretical analysis and similarity solutions for stress boundary layers in viscoelastic flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 153, no. 1, pp. 62–71, 2008.

I. A. Hassanien, “Flow and heat transfer from a continuous surface in a parallel free stream of viscoelastic second-order fluid,” Applied Scientific Research, vol. 49, no. 4, pp. 335–344, 1992.

H. Schichting, Boundary Layer Theory, chapter 97, McGraw Hill, New York, NY, USA, 6th edition, 1964.

T. Hayat, C. Fetecau, Z. Abbas, and N. Ali, “Flow of a Maxwell fluid between two side walls due to a suddenly moved plate,” Nonlinear Analysis: Real World Applications, vol. 9, no. 5, pp. 2288–2295, 2008.

S. Shateyi, “A new numerical approach to MHD flow of a Maxwell fluid past a vertical stretching sheet in the presence of thermophoresis and chemical reaction,” Boundary Value Problems, vol. 2013, article 96, 2013.

C. Fetecau, M. Jamil, C. Fetecau, and I. Siddique, “A note on the second problem of Stokes for Maxwell fluids,” International Journal of Non-Linear Mechanics, vol. 44, no. 10, pp. 1085–1090, 2009.

M. B. Ashraf, T. Hayat, S. A. Shehzad, and A. Alsaedi, “Mixed convection radiative flow of three dimensional Maxwell fluid over an inclined stretching sheet in presence of thermophoresis and convective condition,” AIP Advances, vol. 5, no. 2, Article ID 027134, pp. 761–768, 2015.

M. Awais, T. Hayat, A. Alsaedi, and S. Asghar, “Time-dependent three-dimensional boundary layer flow of a Maxwell fluid,” Computers and Fluids, vol. 91, pp. 21–27, 2014.

K. R. Rajagopal, “Boundary layers in non-linear fluids,” in Trends in Applications of Mathematics to Mechanics, M. D. P. Monteivo Marques and J. F. Rodriques, Eds., vol. 77 of Pittman Monographs and Surveys in Pure and Applied Mathematics, Longman, 1995.

D. V. Boger and K. Walters, Rheological Phenomena in Focus, Elsevier Science, Amsterdam, Netherlands, 1993.

M. Renardy and X. Wang, “Boundary layers for the upper convected Maxwell fluid,” Journal of Non-Newtonian Fluid Mechanics, vol. 189-190, pp. 14–18, 2012.

M. Renardy, “High Weissenberg number boundary layers for the upper convected Maxwell fluid,” Journal of Non-Newtonian Fluid Mechanics, vol. 68, no. 1, pp. 125–132, 1997.

T. Hagen and M. Renardy, “Boundary layer analysis of the Phan-Thien-Tanner and Giesekus model in high Weissenberg number flow,” Journal of Non-Newtonian Fluid Mechanics, vol. 73, no. 1-2, pp. 181–189, 1997.

M. Renardy, “Prandtl boundary layers for the Phan-Thien TANner and Giesekus fluid,” Zeitschrift für angewandte Mathematik und Physik, vol. 66, no. 3, pp. 1061–1070, 2015.

M. Renardy, “Wall boundary layers for Maxwell liquids,” Archive for Rational Mechanics and Analysis, vol. 152, no. 2, pp. 93–102, 2000.

M. Renardy and X. Wang, “Well-posedness of boundary layer equations for time-dependent flow of non-Newtonian fluids,” Journal of Mathematical Fluid Mechanics, vol. 16, no. 1, pp. 179–191, 2014.

M. Renardy, “The initial value problem for creeping flow of the upper convected Maxwell fluid at high Weissenberg number,” Mathematical Methods in the Applied Sciences, vol. 38, no. 5, pp. 959–965, 2015.

G. I. Ogilvie and M. R. Proctor, “On the relation between viscoelastic and magnetohydrodynamic flows and their instabilities,” Journal of Fluid Mechanics, vol. 476, pp. 389–409, 2003.

R. B. Bird and R. C. Armstrong, Dynamics of Polymeric Liquids, vol. 1, John Wiley & Sons, New York, NY, USA, 2nd edition, 1987.

H. I. Andersson and V. Kumaran, “On sheet-driven motion of power-law fluids,” International Journal of Non-Linear Mechanics, vol. 41, no. 10, pp. 1228–1234, 2006.

J. D. Evans, “Re-entrant corner flows of the upper convected Maxwell fluid,” Proceedings of The Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, vol. 461, no. 2053, pp. 117–142, 2005.

M. Renardy, “A matched solution for corner flow of the upper convected Maxwell fluid,” Journal of Non-Newtonian Fluid Mechanics, vol. 58, no. 1, pp. 83–89, 1995.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes, vol. 77 of Fortran, Cambridge University Press, New York, NY, USA, 2nd edition, 2007.