- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Advances in Mechanical Engineering
Volume 2013 (2013), Article ID 612747, 9 pages
An Approach to Define the Heat Flow in Drilling with Different Cooling Systems Using Finite Element Analysis
Department of Mechanical Engineering, Federal University of São João del Rei, Praça Frei Orlando, 170 Centro, 36.307-352 São João del Rei, MG, Brazil
Received 5 May 2013; Revised 17 September 2013; Accepted 20 September 2013
Academic Editor: Shan-Tung Tu
Copyright © 2013 Carlos Henrique Lauro et al. 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.
The heat generated in the cutting zone with high-speed drilling causes damage in the machined part. The heat can affect the dimensions of the hole considering its diameter. Moreover, the heat reduces tool life of uncoated and coated tools. This paper shows experimental tests with high-speed drilling in hardened steel. Drilling was performed on AISI H13 steel with dimensions of 100 × 40 × 14 mm and 52 HRC. The work pieces were drilled with coated drills (TiAlN). A flooded lubricant system and the minimal quantity of lubricant (MQL) were applied to investigate the ability to remove heat from the cutting zone and to compare with dry tests. FEM was applied to define the heat flow and the coefficient of convection for the cooling systems. A steepest descent method was employed to minimize the difference between empirical and simulation data. The results showed that the simulation technique used to find values for heat flow and the coefficient of convection were close to the literature reference. In addition, the adjustment errors of the simulated temperature curves were less than 10% when compared with trial curves. Furthermore, the MQL showed a capability of cooling 3.5 times higher than that of the flooded system.
Drilling is among the most common machining processes. A standard drill has several parameters in its design such as the tip angle, chisel edge angle, chisel edge length, cutting lip length, and helix angle. Each parameter affects the heat generated, the cutting force, and the hole qualities in various ways . Cutting fluids are used in the machining process, mainly in the drilling process, in order to improve the tribological characteristics of the machining considering the increase of cooling and decrease of friction .
The high temperatures generated in the cutting zone decrease the tool life providing the decrease of the cutting speed. The investigation of thermal characteristics during the machining process is significant due to the increase of the material removal rate using high-speed cutting . The studies of the temperature during the drilling are important, because the effects of the heat are generally more severe in drilling than in the other cutting processes. The drilling process achieves high temperatures because it produces a hole in a work piece generating heat in a small region. The increase of temperatures reduces the tool life and causes serious damage to the surface quality of the products .
Park et al.  studied the distribution of the droplets in MQL system. According to the authors, the higher nozzle pressure provided more droplets, but the smaller droplets were obtained in all experimental tests. However, high air pressure does not provide the condition to wet the cutting zone, especially for great nozzle distances. Machining processes as drilling have several restrictions on the use of MQL system due to the difficulty of the air/oil mixture to reach the cutting area.
Bhowmick et al.  studied the MQL application in the tapping of Al-6.5% Si (319 Al) alloy. Although tapping is also a complex process to be lubricated, the authors supported that the use of MQL at the rate of 80 mL/h was similar to that of flooded system considering torque values and the temperature decreased to 5°C. Fromentin et al.  tested nine specific oils using MQL systems and two water-based emulsions and concluded that the efficiency of the oil depends mostly on its additives and less on its viscosity. Based on this, the authors supported that the performance of lubricant tested in form tapping is determined by the molecular structure of the sulfur additive and its thermal behavior. Thus, according to  the decrease of temperature is proportional to the amount of lubricant applied, and the MQL systems can be used to minimize the attrition and the cooling into the cutting zone with the same efficiency than the flooded condition.
The lead time for transformation of research in technology and practical applications is relatively long . This, along with other factors, is directly related to the methodologies used in research centers. Exploratory tests and their parameters have a direct influence on machining. The empirical approaches in the machining studies are usually expensive and require a lot of time. The alternatives used by modern researchers are the mathematical and computational models, in particular the finite elements analysis .
Lazoglu and Altintas  carried out experiments in interrupted orthogonal machining using a computer algorithm simulating a milling process. Wu and Di Han  predicted the maximum temperature in dry drilling using finite element models. The comparison between experimental tests and simulated was ones carried out by using an infrared camera. The values found for both experiments showed that the finite element method found good agreement between the measured and predicted temperatures.
Bono and Ni  analyzed the temperature profile along the cutting of a drill. The authors used a finite element method and demonstrated that the maximum temperature can occur near the chisel edge. Ulutan et al.  determined the three-dimensional temperature fields on the chip, tool, and work piece. The study was critical and previews the residual stress due to the high values of temperature.
According to Mackerle , the main object of the finite element analysis is to get a computational model for the forecasting of strain, stress, and other loads in the work piece, as well as the loads on tool cutting parameters. Several studies using FEA were conducted in machining, with the objective of measuring the internal strain [13–15] and the effects of the thermal distortions in the holes machined in dry condition and defining the maximum temperature value in the chisel and the lips of the drills [16–19].
Brandão et al. [20, 21] used the finite element method to define the heat flow in the milling of AISI H13 and AISI D2 and the heat flow in tapping of AISI H13 of hardened steels. The model proposed and implemented was capable of estimating the energy transferred from the chip formation to the work piece and the average coefficient of convection, based on the curve temperature versus time. Thus, in this paper, a finite element method was employed to define the heat “” and coefficient of convection “” using three cutting conditions; dry drilling, minimal quantity of lubricant, and flooded system.
2. Materials and Methods
2.1. Experimental Tests
Drilling tests were carried out in the AISI H13 hardened steel with 100 × 40 × 14 mm and 52 HRC. The average chemical composition of AISI H13 was 0.40% C, 0.95% Si, 0.31% Mn, 0.011% P, 0.006% S, 5% Cr, 0.12% Ni, 1.25% Mo, and 0.13% Cu. Solid carbide drills had single TiAlN coating with a total length of 89 mm and a helix length of 24 mm, according to DIN 6537K. A milling machine center with 7,500 rpm and 7.5 kW of main power was used. The experimental conditions were chosen according to the tool’s supplier and are shown in Table 1.
Figure 1 shows the positions of the thermocouples embedded at the distance of 0.8 mm from the diameter of the hole, denominated as Dist_0.8. The same configuration was used for the distances of 3.3 and 5.7 mm (Dist_3.3 and Dist_5.7). The insertion of the thermocouples in three different distances had the aim of determining the average temperature in three axial points of the work pieces. Thermocouple type “T” of copper/constantan was used due to the small diameter of the wire and the simplicity of insertion in the holes of the work pieces.
The configuration of the holes in the work piece was designed to maximize the use of them, as can be observed in Figure 1. The experiments were repeated three times for each condition of cooling. The three central holes were used to measure the thrust force and torque in another work not considered in this experiment. The aim was to measure the temperature throughout the work piece while disregarding the values of temperature on the tool and the chip.
The thermocouples were inserted in the holes with a diameter of 2.3 mm using a nylon tube as isolator according to Figure 2. The thermocouples were short and connected to an electronic circuit, which served as a signal conditioner, while a data acquisition cardboard captured the signals. The acquisition rate was 10 Hz, and the cutting fluid for both systems had a viscosity of 18 cSt at 40°C.
The MQL system was assembled with two nozzles 15 mm of the cutting zone. The supplier recommends an input pressure between 0.55 and 1.05 MPa, and it was regulated in the tests to 0.59 MPa. The systems were set for 20 mL/h, both with integral mineral oil, and assembled at 20 mm from the hole entrance.
The input pressure recommended by the supplier of the oil (MQL) should be a range between 0.55 and 1.05 MPa. Thus, the setup chosen was 0.59 MPa with an outflow of 20 mL/h. The cutting speed and feed rate for the dry tests were the same as those of the experiments using coolant. Figure 2 shows the layout for the experiments and Figure 3 shows a diagram of the MQL system.
The MQL system consists of a small amount of cutting fluid that is vaporized in the cutting process due to the high pressure of air. The device has a mixing chamber that receives oil and air with high pressure and injects the spray in the cutting process in a percentage of 20% of oil in the form of droplets. In Figure 3 the letter “A” in the circle represents the line that came from the machine oil tank and “B” represents the compressed air line that came from the accumulator located in the laboratory.
2.2. Finite Element Methodology
The objective was to provide theoretical temperature curves in order to minimize the error between the experimental and simulated curves using the finite element method. A computational routine was developed allowing the introduction of the values of the heat “” and coefficient of convection “” during the simulation. The same strategy was applied in the tapping process and allowed to define the heat flux and adjust the experimental and theoretical data Bono and Ni  and Brandão et al. . Heat source theory was adopted to calculate the temperature of a one-dimensional moving point. The heat source moved as a moving point source in a rod that loses heat by convection. The equations used in the simulation were Fourier’s law for heating considering the thermal contact conductance and Newton’s law for cooling, which states that the rate of heat loss of a body is proportional to the difference in temperatures between the body and its surroundings.
Initially, the values of heat were previously introduced in the finite element method in joules and then divided by the time of the simulation in seconds to set the value of heat flow in watts. The initial values were considered based on the studies carried out by Brandao et al.  on drilling processes. In this work the authors used a mathematical model with a good fit to define the heat flux in drilling. The values found in the previously mentioned work were used as an initial heat source. The initial boundary conditions were the environmental conditions at a room temperature of 24°C, which was the same recorded by the thermocouples. The simulation of the drilling process was developed using a mesh with 10 layers in the direction of the tool’s feed speed containing 96 elements in each layer. Each layer was divided into a central ring with 32 elements, an intermediate ring with 32 elements, and an external ring with 32 elements. Figure 4 shows a top view of the work piece demonstrating the three rings.
The number of 10 layers with 96 elements in each layer was based on a specific strategy because of the limitations of the software. This strategy was adopted because of the long time required by the simulation and because the FEM had a limited number of load steps available.
The simulation was performed in such a way that the heat “” was applied to the elements of the central ring in the first layer during the time which corresponded to the time of drilling divided by the number of layers. This heat was applied initially considering a specific value for “” in the program for all 32 elements simultaneously and based on the previous information . Afterwards, the heat was removed and the coefficient of convection “” was applied during the equivalent time. In the same way that was done with the heat “”, a value for “” was originally defined and applied during simulation.
Figure 5(a) shows the strategy used for the first stage. In the second stage of the simulation, which corresponds to a back layer, the heat “” came back to be applied to the elements of the central ring of the second layer and to the elements of the first layer of the intermediate ring. Simultaneously, the heat “” was removed and substituted for the coefficient of convection “,” as shown in Figure 5(b). Finally, the heat “” was applied simultaneously to the elements of the central ring of the third layer, the elements of the intermediate ring of the second layer, and the external ring of the first layer, whereby the heat “” was removed and the coefficient of convection “” was applied again, as shown in Figure 5(c).
The simulation used this procedure in the same order of applying heat “” and coefficient of convection “” in the direction of the feed rate at the end of the process when the tool passed along by the last layer of the external ring. This strategy was conducted with the objective of representing as closely as possible the drill geometry and the feed rate of the tool in the drilling.
All elements used in the simulation were considered with the same geometry in order to give accuracy in the simulation process. The use of identical geometry allowed the introduction of heat with a constant time. The time corresponds to the drilling process divided by the number of elements.
The temperature was measured in a ring with 0.8 mm more than the outer diameter of the holes. The positions of the measurement temperature point during the simulation were the nodes between each mesh. Considering the feed rate, temperatures were recorded between second/third layer for the first thermocouple, fifth/sixth for the second thermocouple, and seventh/eighth for the last thermocouple, considering the axial direction.
The values of “” and “” were initially estimated based on machining references. Those values corresponded to an average of 10% of the energy generated in the machining process as shown in (1). According to Shaw , the quantity of heat that goes to the tool corresponds to 12%, another parcel of heat that goes to the work piece is 10%, and the remaining energy that goes to the chips corresponds to 78% in average values, approximately. Thus, the increased heat in the work piece is where “” is the mass of the work piece, “” is the specific heat, and “” is the increase of the temperature throughout the work piece. The values of the coefficient of convection were previously established according to some references in the heat conduction Ozisik , Bejan , and Dewitt and Incropera . The simulations have generated nine curves where the index “” is the position of the thermocouple, , for each distance from the maximum diameter, that is, Pos_1, Pos_2, and Pos_3, at 0.8, 3.3, and 5.7 mm. The index “” is the distance of the thermocouple from the drill entrance, that is, T0, T1, and T2, at 3.0, 7.0, and 11.0 mm. After obtaining the simulated curves a comparison was performed to define the error with the experimental curves . Thus, the error for each comparison can be calculated according to where “” is the number of points during the simulation. In this stage, there is a mean quadratic error for each curve and an overall average error for each cooling system according to (3):
In order to define the best estimated values for “” and “” for minimizing the average error function Eravg, the steepest descent method was applied. This method estimates the first point of the function and moves to point in the direction of the regional downhill gradient . This method was adopted due to its simplicity, and it leads to a local minimum, which seems good enough since the energy, for example, has its limit, roughly between zero and the cutting power. The coefficient of convection has also some known limits [23, 24, 26].
To define the gradients in both directions, simulation steps in four new points were established around the first estimated one. Each error function was calculated, and the steepest gradient that followed down the error valley was found. At this point, the best combination of “” and “” had been found for each experimental condition. The initial energy was selected as 50 joules, and the coefficient of convection depends on the cooling system. The whole minimum search procedure was implemented using a program developed using a mathematical software.
The end criterion for stopping the fit tests was less than 0.1 for temperature values considering the or 5 steps in any direction. The machined surface represented by the face 6 receives energy from the moving source (drill) that represents the region where the material was removed, as seen in Figure 6. Figure 6 shows that in the same face 6 the red line represents the loss of heat according to the cooling system being simulated; this situation occurs just for the first element.
Considering the second layer, all cooling systems do not act directly. There is an exchange of energy in the surfaces 1, 2, 3, 4, and 5 with the other elements of the work piece, which are not affected by the cutting process.
3. Results and Discussion
The number of points recorded by the board acquisition was filtered proportionally to the simulation points. Thus, the experimental results were compared with the simulation data. This strategy was adopted due to the long time required by the simulation and because the finite element method was available for a limited number of load steps. Figures 7, 8, and 9 show the experimental temperature graphs in dry, MQL, and flooded tests for the position at 0.8 mm. It can be observed that all graphs of temperature have a similar behavior with differences only in the temperature peaks.
Although the drilling operation lasted about 30 seconds, the temperature was recorded for 300 seconds allowing complete stabilization. The peaks of temperature are different for each condition tested, being higher for the dry test and lower for the flooded test. However, the maximum peaks for the three systems, flooded, MQL, and dry, were very close to each other. Figures 7, 8, and 9 show that the temperature peaks occurred at around 30 seconds for the dry, MQL, and flooded condition, respectively.
Considering the theory, the highest temperature values for each thermocouple were registered when the drill passed nearest to its position. However, during the experiments, this situation did not occurr because the temperature peaks were recorded before the passing of the tool. That can be explained because, as the drill goes deep, more heat is produced and accumulated like successive heat waves being created. The thermocouples located ahead of the tool will receive and register the waves, since feed speed is relatively low.
Considering the three thermocouples, Figure 10 shows a comparison between the experimental temperature curves and simulation temperature curves. This temperature graph is in the MQL system, and it can be observed that the error is minimal for this cooling condition. The graphs for the experimental tests under the flooded and dry conditions showed a similar behavior, as can be seen in Figures 11 and 12. Thus, one can observe variations just in the maximum temperature peaks.
Figures 10, 11, and 12 show the variation of the temperature in the simulation for the all thermocouples (TMP1, TMP2, and TMP3). Considering the behavior of the thermocouples, the fitted values of temperature were not too similar to the others. This situation occurred because the first thermocouple is located closer to the upper surface, which provides a higher percentage of the cooling rate due to sprayed air oil particles.
According to the experimental results, the simulations show an increase of temperature from the first to the third thermocouple, in the early stages of the drilling process, due to a heat wave that arrives in a small interval of time. However, for the MQL system, the temperature values at the end of drilling are the highest due to the tip of the tool receiving heat for a long time and perfect cooling not occurring in the cutting zone. This situation can explain the difference in the fit between the experimental and simulated temperature curves. The software adjusted the “” and “” values simultaneously. If the adjustment of the convection coefficient has a great difference, the adjusted value of heat may be affected. However, according to Figures 10, 11, and 12 the difference is small, and its influence on heat is negligible.
Tables 2, 3, and 4 show the results of the simulation obtained with minimal errors for dry, flooded, and MQL, respectively. It is assumed that the starting point of the simulation began with 100 W/(m2·C) for the coefficient of convection. The value of heat was defined in function of the strain energy during a drilling process that is converted to thermal energy, according to Shaw . Considering the references, it is safe to assume as a first approximation that all the energy associated with chip formation is converted to heat. The energy retained in the work piece associated with this process is almost 10%. The strategy used in the simulated tests allows a variation of the heat around the coefficient of convection until a small error for the pair of variables is found. This strategy of simulation was repeated for each pair of variables shown in Tables 2, 3, and 4.
Table 5 shows a summary of the heat values and the coefficient of convection for the three cooling systems. The results of heat reproduced the same tendency of the temperature curves. It was observed that the heat increases more for the dry system than in the flooded system. The MQL system was located between the two extreme values. The difference between the dry and flooded tests, considering the minimum variation, was a range of 30 to 40 joules. Considering the drilling time as 30 seconds, only during the performance of the cooling systems, the heat flux was 1.3 watts, supporting that the data are very consistent. In the work carried out by Brandao et al.  values ranging from 0.26 to 1.03 W were found for dry and cooling conditions, respectively.
The value of the coefficient of convection was expected to have an inverse tendency of the heat value. However, for the MQL system, the coefficient of convection was bigger than that for the flooded. This situation can be explained due to the heat for the MQL system being greater than for the flooded. Therefore, as the two variables “” and “” were simulated together, the thermal balance of energy compensates the increase of one with the reduction of the other. However, the values for the coefficient of convection found were very different from the results of previous works. The highest value found for the convection coefficient is three times smaller than the smallest value found in previous work . This situation demonstrates that the prediction of coefficients of convection in the machining process is complex and depends on several parameters. However, the errors of minimization between experiment and simulation stayed below 10%. Based on this, the applied model can be considered to have simulated well the drilling process, but the use of a larger number of load steps can improve the methodology and provide a better fit.
In this work, an approach was presented whereby both heat and the coefficient of convection over three cooling systems could be obtained from simulation tests.(i)The simulation was capable of representing additionally the tendency of the increase of the temperature and different values for each cooling condition. These values were less for the flooded and higher for the dry condition.(ii)The values for heat between 30 and 40 joules provide an order of magnitude of 1 to 1.3 watts, considering the machining time. The coefficient of convection for the MQL system does not show a specific tendency of increasing according to the simulation, because the MQL system showed values 3.5 times higher than those of the flooded system.(iii)The coefficient of convection for the flooded and dry systems was somewhat different to that of the reference works. Considering the MQL system, there is no data sheet due to this system being relatively new. The errors of minimization between experiment and simulation stayed below 10% and can be considered to be acceptable, based on the theoretical model which only used 960 load steps.
The authors thank Gerdau Piratini and SANDVIK for providing the steel AISI H13 and tooling, respectively.
- M. Pirtini and I. Lazoglu, “Forces and hole quality in drilling,” International Journal of Machine Tools and Manufacture, vol. 45, no. 11, pp. 1271–1281, 2005.
- M. Soković and K. Mijanović, “Ecological aspects of the cutting fluids and its influence on quantifiable parameters of the cutting processes,” Journal of Materials Processing Technology, vol. 109, no. 1-2, pp. 181–189, 2001.
- E. Bagci and B. Ozcelik, “Finite element and experimental investigation of temperature changes on a twist drill in sequential dry drilling,” International Journal of Advanced Manufacturing Technology, vol. 28, no. 7-8, pp. 680–687, 2006.
- M. Bono and J. Ni, “The effects of thermal distortions on the diameter and cylindricity of dry drilled holes,” International Journal of Machine Tools and Manufacture, vol. 41, no. 15, pp. 2261–2270, 2001.
- K.-H. Park, J. Olortegui-Yume, M.-C. Yoon, and P. Kwon, “A study on droplets and their distribution for minimum quantity lubrication (MQL),” International Journal of Machine Tools and Manufacture, vol. 50, no. 9, pp. 824–833, 2010.
- S. Bhowmick, M. J. Lukitsch, and A. T. Alpas, “Tapping of Al-Si alloys with diamond-like carbon coated tools and minimum quantity lubrication,” Journal of Materials Processing Technology, vol. 210, no. 15, pp. 2142–2153, 2010.
- G. Fromentin, A. Bierla, C. Minfray, and G. Poulachon, “An experimental study on the effects of lubrication in form tapping,” Tribology International, vol. 43, no. 9, pp. 1726–1734, 2010.
- L. C. Brandao, R. T. Coelho, and C. H. Lauro, “Contribution to dynamic characteristics of the cutting temperature in the drilling process considering one dimension heat flow,” Applied Thermal Engineering, vol. 31, no. 17-18, pp. 3806–3813, 2011.
- H. Schulz, “The history of high-speed machining,” in Proceedings of 5th International Scientific Conference on Production Engineering, Opatija, Croatia, 1999.
- J. Mackerle, “Finite-element analysis and simulation of machining: a bibliography (1976–1996),” Journal of Materials Processing Technology, vol. 86, no. 1–3, pp. 17–44, 1998.
- I. Lazoglu and Y. Altintas, “Prediction of tool and chip temperature in continuous and interrupted machining,” International Journal of Machine Tools and Manufacture, vol. 42, no. 9, pp. 1011–1022, 2002.
- J. Wu and R. Di Han, “Technical paper: a new approach to predicting the maximum temperature in dry drilling based on a finite element model,” Journal of Manufacturing Processes, vol. 11, no. 1, pp. 19–30, 2009.
- M. Bono and J. Ni, “The location of the maximum temperature on the cutting edges of a drill,” International Journal of Machine Tools and Manufacture, vol. 46, no. 7-8, pp. 901–907, 2006.
- D. Ulutan, I. Lazoglu, and C. Dinc, “Three-dimensional temperature predictions in machining processes using finite difference method,” Journal of Materials Processing Technology, vol. 209, no. 2, pp. 1111–1121, 2009.
- M. H. Dirikolu, T. H. C. Childs, and K. Maekawa, “Finite element simulation of chip flow in metal machining,” International Journal of Mechanical Sciences, vol. 43, no. 11, pp. 2699–2713, 2001.
- K. C. Ee, O. W. Dillon Jr., and I. S. Jawahir, “Finite element modeling of residual stresses in machining induced by cutting using a tool with finite edge radius,” International Journal of Mechanical Sciences, vol. 47, no. 10, pp. 1611–1628, 2005.
- N. A. Abukhshim, P. T. Mativenga, and M. A. Sheikh, “Heat generation and temperature prediction in metal cutting: a review and implications for high speed machining,” International Journal of Machine Tools and Manufacture, vol. 46, no. 7-8, pp. 782–800, 2006.
- S. Y. Hong and Y. Ding, “Cooling approaches and cutting temperatures in cryogenic machining of Ti-6Al-4V,” International Journal of Machine Tools and Manufacture, vol. 41, no. 10, pp. 1417–1437, 2001.
- A. Niku-Lari, J. Lu, and J. F. Flavenot, “Measurement of residual-stress distribution by the incremental hole-drilling method,” Journal of Mechanical Working Technology, vol. 11, no. 2, pp. 167–188, 1985.
- L. C. Brandão, R. T. Coelho, and A. R. Rodrigues, “Experimental and theoretical study of workpiece temperature when end milling hardened steels using (TiAl)N-coated and PcBN-tipped tools,” Journal of Materials Processing Technology, vol. 199, no. 1, pp. 234–244, 2008.
- L. C. Brandão, R. T. Coelho, and A. T. Malavolta, “Experimental and theoretical study on workpiece temperature when tapping hardened AISI H13 using different cooling systems,” Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol. 32, no. 2, pp. 154–159, 2010.
- M. C. Shaw, Metal Cutting Principles, Oxford University, England, UK, 2005.
- N. M. Özisik, Heat Conduction, John Wiley & Sons, 2nd edition, 1980.
- A. Bejan, Heat Transfer, John Wiley & Sons, 2nd edition, 1996.
- W. H. Flannery, B. P. Teukolsky, and W. T. Vetterling, Numerical Recipes in Pascal—The Art of Scientific Computing, Cambridge University Press, Cambridge, UK, 1990.
- D. P. Dewitt and F. P. Incropera, Fundamentals of Heat and Mass Transfer, 3 edition, 1990.