Mathematical Problems in Engineering
Volume 2016 (2016), Article ID 5387612, 14 pages
http://dx.doi.org/10.1155/2016/5387612

Research Article

Maximum Temperature and Relaxation Time in Wet Surface Grinding for a General Heat Flux Profile

J. L. González-Santander

Universidad Católica de Valencia “San Vicente mártir”, C/ Guillem de Castro 94, 46001 Valencia, Spain

Received 26 November 2015; Accepted 24 February 2016

Academic Editor: John D. Clayton

Copyright © 2016 J. L. González-Santander. 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.

We solve the boundary-value problem of the heat transfer modeling in wet surface grinding, considering a constant heat transfer coefficient over the workpiece surface and a general heat flux profile within the friction zone between wheel and workpiece. We particularize this general solution to the most common heat flux profiles reported in the literature, that is, constant, linear, parabolic, and triangular. For these cases, we propose a fast method for the numerical computation of maximum temperature, in order to avoid the thermal damage of the workpiece. Also, we provide a very efficient method for the numerical evaluation of the transient regime duration (relaxation time).

1. Introduction

Grinding is a machining industrial process of metallic plates that consist in an abrasive wheel rotating at a high speed over the surface of a workpiece. The abrasion that produces the grinding wheel removes the surface material of the metallic plate being ground. During the grinding process, most of the energy is converted into heat, which accumulates within the contact zone between the wheel and the workpiece [1]. The high temperatures reached during this machining process may produce an unacceptable decrease in the quality of workpieces, such as burning or residual stresses [2]. In order to avoid this thermal damage of the workpiece, liquid coolant is usually injected into it, so power generation by friction is reduced and cooling by convection occurs. Therefore, a theoretical model that can predict the maximum temperature in wet grinding for the online monitoring process is demanded.

In surface grinding, the heat transfer inside the workpiece usually is modeled by a strip heat source infinitely long and of width (m in SI units), which moves at a speed (m s⁻¹) over a semi-infinite solid surface (see Figure 1). Setting the Cartesian coordinate system fixed to wheel, as it is shown in Figure 1, the temperature field of the workpiece with respect to the room temperature (K) must satisfy the convective heat equation [3, §1.7(2)]:

where is the thermal diffusivity (m² s⁻¹). Since initially the workpiece is at room temperature, (1) is subjected to an homogeneous initial condition:

Figure 1: Setup in surface down-grinding.

Wet grinding can be modeled assuming a constant heat transfer coefficient (W m⁻² K⁻¹) on the workpiece surface [4]. Assuming as well a dimensionless heat flux profile within the contact area between wheel and workpiece, we have the following boundary condition:

where denotes the Heaviside function, (W m⁻¹ K⁻¹) is the thermal conductivity, and (W m⁻²) is the average heat flux entering into the workpiece along the contact width ; thus,

In the case of dry grinding, that is, , boundary-value problem (1)–(3) can be solved in integral form for the case of a constant heat flux profile [5, Eqn.22]. Also, considering a general heat flux profile, a closed form expression for the surface temperature in the stationary regime, that is, , can be found [6].

For the wet case, that is, , (1)–(3) has been solved for constant [4] and linear [7] heat flux profiles, both in the stationary regime. The scope of this paper is to generalize these results and calculate, in terms of an integral expression, the time-dependent temperature field satisfying (1)–(3) for any constant , considering any heat flux profile. From this general expression, we will recover the expressions for the particular cases found in the literature (constant [8, 9] and linear [7] heat flux profiles), calculating also solutions for other heat flux profiles found in the literature, such as triangular [10–15] and parabolic [16]. These latter solutions do not seem to be reported in the literature. Further, we will provide a simple and rapid method for the calculation of the maximum temperature reached in the workpiece and some approximated expressions for the duration of the transient regime, for all the heat flux profiles considered (i.e., constant, linear, triangular, and parabolic). These approximate expressions seem to be novel as well. It is worth noting that the computation of the relaxation time is essential for analyzing the transient regime in which the wheel is engaged and disengaged from the workpiece (cut in and cut out) [17].

This paper is organized as follows. Section 2 derives the solution of the value boundary problem stated in (1)–(3) for any analytic heat flux profile . Section 3 particularizes the general solution of the previous section to the most common heat flux profiles reported in the literature, that is, constant, linear, and parabolic. We also give an expression for a triangular heat flux profile, although in this case is not an analytic function. Section 4 provides a very efficient method for searching the maximum temperature in the workpiece for the different heat flux profiles considered. In Section 5, for each heat flux profile considered, we propose the root searching of an equation in order to evaluate numerically the duration of the transient regime. For this purpose, we provide analytical approximations to these roots, to be used as starting iteration point in the root finding. Section 6 presents some numerical simulations of the workpiece surface temperature in the stationary regime as well as the maximum temperature computation by using the expressions derived in previous sections. We also plot the surface temperature evolution during the transient regime. Our conclusions are summarized in Section 7.

2. General Heat Flux Profile Solution

In order to solve the problem given in (1)–(3), we can build the solution by using the method of Green’s function. For this purpose, let us denote the Cartesian coordinates fixed to the workpiece as (see Figure 2), where, at , both coordinates systems are overlapped.

Figure 2: Relationship between the coordinates systems, one fixed to the wheel and the other one to the workpiece.

Green’s function is interpreted in this case as the temperature field rise evaluated at in a semi-infinite body () in which, at position , an instantaneous point heat source of energy (J) appears at the instant and where the surface () is subject to a constant heat transfer coefficient [3, §14.9 (4)]:

where is the density of the body (kg m⁻³) and its specific heat (J kg⁻¹ K⁻¹) and where we have defined the radiation coefficient (m⁻¹) as

A formal derivation of (5) is given in [18, Appendix].

According to Figure 2, a point of the heat source has coordinates in and in . The relationship between both coordinates systems is

Therefore, the temperature rise in the workpiece at coordinates due to a point of the heat source with respect to the coordinate system fixed to the wheel is given by

The derivation of the solution of (1)–(3) by using (8) follows three steps:

(i)Superposition in space to give temperature due to instantaneous line source, which acts on the surface parallel to the -axis.(ii)Superposition of line sources to give temperature due to an instantaneously acting strip source on the surface.(iii)Superposition in time to give temperature due to continuously acting strip source.

2.1. Instantaneous Line Source

Let us consider now an instantaneous line heat source parallel to the -axis. Let us assume in (8) that the energy per unit length of this heat source (J m⁻¹) is constant:

thus, integrating in , we obtain the field temperature due to an instantaneous line heat source as

2.2. Instantaneous Strip Source

Let us consider now a strip source of variable heat strength along the -axis, according to the following profile:

where (J m⁻²) is the energy that this heat source leaves at the instant per surface unit of the strip source. Notice that the function in (11) is dimensionless, as in (3). Also, takes into account the heat flux profile entering into the workpiece within the contact zone with the wheel. Substituting now (11) in (10) and integrating over the strip width, , we have the following temperature field in the workpiece:

Assuming that the function is analytic, we can expand it in its Taylor series, so that, taking into account also that is dimensionless, we have

thus, the integral given in (12) can be written as

According to the result given in Appendix (A.9), we have

where

and where denotes the lower incomplete gamma function [19, Eqn. 45:1:2]. Substituting back these results, we arrive at

2.3. Continuous Strip Source

Let (W m⁻²) be the energy per time unit and per surface unit that the strip heat source leaves over the workpiece surface. Therefore, during an infinitesimal interval of time , we have

Substitution of (18) in (17) and integration over time (remember that ) result in

where we have performed the change of variables . Considering now the dimensionless variables

where is a characteristic length

we finally arrive at

where

3. Particular Cases

In this section, we derive, from general expression (22), particular expressions for the most common heat flux profiles considered in the literature (i.e., constant, linear, triangular, and parabolic). For the constant and linear cases, we obtain expressions already given in the literature, as mentioned before.

3.1. Constant Heat Flux Profile

For a constant heat flux profile satisfying (4), we have

Thus, in (22), we have only the summand which corresponds to , so, from (23) and using the property given in Appendix (A.12), we have

Therefore, (22) reduces to

where the superscript denotes that we are considering a constant heat flux profile. Particularizing (26) to and performing the change of variables , we obtain the following expression for the case of dry grinding:

which is the result given in [5, Eqn.22]. Undoing the change of variables and performing the limit , we get the stationary regime for the dry case [5, Eqn.7]:

where is the zeroth-order modified Bessel function of the second kind [20, Sect. 9.6]. Therefore, taking into account (29) in (26), we obtain the stationary regime for the wet case as

which is the expression given by [4]. Recently, in [21], the expression given in (30) particularized on the surface () can be expressed in closed form as

where

It is worth noting that (31) does not converge for high Biot numbers, that is, .

3.2. Linear Heat Flux Profile

For a linear heat flux profile satisfying (4), we have

thus, (22) is reduced to

where the superscript indicates a linear heat flux profile and where we have applied (25). Also, according to (A.14), the following result has been applied:

Performing in (34) the change of variables and the limit , we obtain the following expression for the stationary regime, which is given in [7, Eqns.3.46-47]:

where the following function has been defined:

3.3. Triangular Heat Flux Profile

For a triangular heat flux profile satisfying (4) (see Figure 3), we have

where we have defined the following dimensionless parameter:

Figure 3: Heat flux profiles considered.

Although we cannot expand (38) in its Taylor series, we can redo the calculation of Section 2.2, performing the integral given in (12) as follows:

where we have set

By using formulas (A.13) and (A.15) given in Appendix, some algebra leads to the following result:

where we have defined

Substituting (42) into (12), integrating over time , and performing the change of variables , we get the time-dependent temperature field due to a triangular heat flux profile as

which, in dimensionless variables (20), is rewritten as

3.4. Parabolic Heat Flux Profile

For a parabolic heat flux profile satisfying (4) (see Figure 3), we have

thus, (22) is reduced to

where the superscript denotes that we are considering a parabolic heat flux profile. By using the properties given in Appendices (A.12), (A.14), and (A.17), we arrive at

4. Maximum Temperature

From a physical point of view, the maximum temperature must be reached at the stationary regime , because the longer the heat source is acting, the greater the temperature in the workpiece is. Moreover, the location of maximum temperature must be on the surface , within the contact zone between wheel and workpiece, . It is worth noting that recently the latter has been mathematically proved in [22] for constant and linear heat flux profiles. Therefore, since is a differentiable function in , a convenient way to evaluate maximum temperature is to solve numerically all points that satisfies

and then the maximum is given by

In Section 6, we will check numerically that for the most common heat flux profiles proposed in the literature (constant, linear, triangular, and parabolic) the root is unique, so .

For a constant heat flux profile, (49) reads

Similarly, for a linear heat flux profile, we have

For a triangular profile, we obtain

and for a parabolic one

5. Relaxation Time

The stationary regime is asymptotically reached at ; that is,

In order to avoid thermal damage, the most important point is the location of the maximum temperature . Therefore, for estimating how rapidly the stationary regime is reached (i.e., the relaxation time of the transient regime), we can solve the following equation for :

In this section, we are going to provide approximated explicit solutions of (56) for the different heat flux profiles considered. Each of these approximate solutions will be used as starting iteration point for the numerical root searching of (56).

5.1. Constant Heat Flux Profile

For a constant heat flux profile, according to (26), (56) is written as

Let us solve (57) approximately. On the one hand, we find in the literature [20, Eqn.7.1.23] the following asymptotic expansion:

thus, expanding up to and taking , we have

On the other hand, expanding in Taylor series the following function up to the first order, we have

thus,

Substituting (59) and (61) in (57), we get the following approximated equation:

By using the Lambert function [23], (62) can be solved explicitly, obtaining the following approximated expression for the relaxation time:

Notice that (63) does not depend on , so we do not have to compute (57) if we want an estimation of the relaxation time for a constant heat flux profile.

5.2. Linear Heat Flux Profile

Considering now a linear heat flux profile, according to (34), (56) is written as

In order to solve (64) approximately, let us expand in Taylor series the following function up to first order:

thus,

Now, substituting (59), (61), and (66) in (64), we arrive at

that can be solved explicitly as

5.3. Triangular Heat Flux Profile

Let us consider now a triangular heat flux profile, so, according to (45), (56) reads as

Equation (69) can be solved approximately expanding in the Taylor series the following function up to first order, taking into account the function definition (43) and asymptotic formulas (60) and (65):

Taking into account (59) and (70) in (69), after some algebra, we arrive at the following approximated equation:

that can be solved explicitly as

As in the constant case, (72) does not depend on ; moreover, it does not depend on either.

5.4. Parabolic Heat Flux Profile

Finally, for a parabolic heat flux profile, according to (48), (56) is expressed as

Taking into account asymptotic formulas (59), (61), and (65), some algebra leads (73) to the following approximate equation:

which can be solved as

6. Numerical Results

Table 1 shows three sets of parameters (SI units) for the numerical simulations. Data set 1 considers a titanium alloy VT20 workpiece, whose thermal properties are given in [24]. The grinding regime for this simulation can be found in [25]. Data sets 2 and 3 consider plain carbon steel and aluminum oxide, (sapphire), as workpiece material, respectively. The thermal properties and the grinding regimes for these data sets can be found in [26].

Table 1: Simulation parameters in SI units.

Data 1
Data 2
Data 3

Workpiece
Material
VT20
Steel
Sapphire

Grinding regime

Dimensionless parameters

Figure 3 presents the different heat flux profiles considered. All of them have been normalized according to (4); thus, they have got the same heat flux on average. Therefore, we can compare how the distribution of the heat flux within the contact zone affects the graph of the workpiece surface temperature. Also, for all the numerical simulations, we have taken a value of for the triangular profile (45).

Table 2 shows the value of maximum temperature and its location within the friction zone. For this purpose, we have computed the roots of (51)–(54) by using Brent’s method [27], taking as starting interval . Notice that, for every data set, the maximum temperature is ordered depending on the heat flux profile. This order, from lowest to highest, corresponds to a constant, linear, triangular, and parabolic heat flux profile, respectively. This qualitative behavior of maximum temperature agrees with the heat flux distribution of each profile (see Figure 3), in which maximum temperature is higher where input of heat is also higher. Therefore, the location of maximum temperature agrees with the location of maximum heat flux. Thus, considering the same input data, the ordering of these locations from left to right with respect to the type heat flux of profile is as follows: constant, triangular, linear, and parabolic; that is, the more the heat flux profile is skewed to the leading edge, the nearer the maximum temperature is located in that edge. It is worth noting that the numerical evaluation of maximum temperature by the root searching of (51)–(54) is much faster than applying directly Mathematica’s NMaximize command to the surface temperature in the stationary regime, that is, . For instance, taking data set 1 as input parameters and a constant heat flux profile, NMaximize command is ≈ 125 times slower than the numerical root finding of (51). Moreover, NMaximize command fails if we take a triangular heat flux profile for all the data sets of Table 1.

Table 2: Location and value of maximum temperature.

Data 1
Data 2
Data 3

Constant profile

[K]

Linear profile

[K]

Triangular profile

[K]

Parabolic profile

[K]

Figures 4–6 show the surface temperature () in the stationary regime () for the three data sets given in Table 1, respectively. For the numerical evaluation of these plots, we have used the formulas given in Section 3: (26) for constant heat flux profile, (34) for linear profile, (45) for triangular profile, and (48) for a parabolic one. It is worth noting that the numerical integration of the linear and the parabolic cases is quite efficient if we use the double exponential strategy [28]. The shapes of the graphs are quiet similar to the different data sets of input parameters, except at the leading and trailing edges (). This is because the heat source velocity in Figure 4 is much faster than in Figures 5 and 6. Also, Figures 4–6 show that the skewness of the temperature curves follows the corresponding heat flux profile.

Figure 4: Surface temperature for a VT20 titanium alloy workpiece (data 1).

Figure 5: Surface temperature for a carbon steel workpiece (data 2).

Figure 6: Surface temperature for a sapphire workpiece (data 3).

Table 3 shows the relaxation times for the different heat flux profiles considered. For this purpose, we have used the relaxation times approximations given in (63), (68), (72), and (75) for the different heat flux profiles, as starting iteration points for the root searching of (57), (64), (69), and (73), respectively. In general, the approximation formulas provide the order of magnitude of the exact root. Also, the convergence to the exact value is extremely rapid (≈0.01 s).

Table 3: Relaxation times for the different heat flux profiles considered.

[s]
Data 1
Data 2
Data 3

Constant profile
Exact

Approx.

Linear profile
Exact

Approx.

Triangular profile
Exact

Approx.

Parabolic profile
Exact

Approx.

Figures 7–9 show the temperature evolution on the workpiece surface, taking a triangular heat flux profile for the three data sets of Table 1. The black solid line in these figures corresponds to the stationary regime and the gray solid line to the relaxation time (see Table 3). In the case of Figures 8 and 9, the plot corresponding to the relaxation time and the stationary regime are overlapped. For these cases, this is due to the fact that the relaxation time is greater than the time that the heat source needs to cover the friction zone width; that is, . The plots are exponentially discretized over time because the transient regime occurs very rapidly. Notice that the first stages of the temperature evolution resemble the triangular heat flux profile, this being more clear for data set 1 in Figure 7.

Figure 7: Evolution of the surface temperature in a titanium alloy VT20 workpiece (data set 1).

Figure 8: Evolution of the surface temperature in a carbon steel workpiece (data set 2).

Figure 9: Evolution of the surface temperature in a sapphire workpiece (data set 3).

7. Conclusions

By using the Green function, we have solved the boundary-value problem that models the heat transfer in wet surface grinding, considering an arbitrary heat flux profile and a constant heat transfer coefficient on the workpiece surface. We have particularized this result to the most common heat flux profiles reported in the literature, that is, constant, linear, parabolic, and triangular. In constant case (26) and linear case (34), we recover expressions given in the literature. Nonetheless, the triangular case (45) and parabolic case (48) do not seem to be found in the literature.

From the analytical expressions found for the surface temperature in the stationary regime for the different heat flux profiles considered, we propose the root searching of some equations in integral form (see (51)–(53)), for the maximum temperature computation. It turns out that this method is very fast (≈1 s) and also much more rapid than the direct numerical maximization of the surface temperature in the stationary regime. Therefore, it provides a very useful aid to prevent thermal damage in the online monitoring grinding process.

Finally, we have provided a very rapid method (≈0.01 s) for the computation of the relaxation time for the different heat flux profiles considered. For the initial guess of this root searching, we have given analytical expressions that in general provide the order of magnitude of the relaxation time.

Appendix

Let us calculate the integral

performing the change of variables , we get

Taking into account Newton’s formula [19, Eqn. 6:14:1], we can rewrite (A.2) as

Now, consider the integral

When , perform the change of variables , so

where

denotes the lower incomplete gamma function [19, Eqn. 45:3:1]. Notice that, performing the change of variables , we have

Therefore, for , we have

Applying (A.8) to (A.3) results in

where we have defined

According to the formula [19, Sect. 45:4]

and knowing that [19, Eqn. 40:5:1] and [19, Eqn. 8:0:1], for , (A.10) is reduced to

so

For , applying directly [19, Sect. 45:4], we have

so

Similarly, for , taking into account (A.12) and (A.14) and applying the recursion formula [19, Eqn. 45:5:1]

we arrive at

so

Competing Interests

The author declares that there are no competing interests regarding the publication of this paper.

Acknowledgments

The author wishes to thank the financial support received from Generalitat Valenciana under Grant GVA/2015/007.

References

S. Malkin and C. Guo, Grinding Technology: Theory and Application of Machining with Abrasives, Industrial Press, New York, NY, USA, 2008.

J. B. J. W. Hegeman, J. T. M. De Hosson, and G. De With, “Grinding of WC—Co hardmetals,” Wear, vol. 248, no. 1-2, pp. 187–196, 2001.

H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Oxford Science Publications, Oxford, UK, 1986.

N. R. DesRuisseaux and R. D. Zerkle, “Temperature in semi-infinite and cylindrical bodies subjected to moving heat surfaces and surface cooling,” Journal of Heat Transfer, vol. 92, pp. 456–464, 1970.

J. C. Jaeger, “Moving sources of heat and the temperature at sliding contracts,” Proceedings of the Royal Society of New South Wales, vol. 76, pp. 204–224, 1942.

J. L. González-Santander, “Calculation of an integral arising in dry flat grinding for a general heat flux profile. Application to maximum temperature evaluation,” Journal of Engineering Mathematics, vol. 88, pp. 137–160, 2014.

W. J. Sauer, Thermal aspects of grinding [PhD dissertation], Carnegie Mellon University, Pittsburgh, Pa, USA, 1971.

M. C. Shaw, “Temperatures in cutting and grinding,” American Society of Mechanical Engineering, vol. 146, pp. 17–24, 1990.

A. S. Lavine, S. Malkin, and T. C. Jen, “Thermal aspects of grinding with CBN wheels,” CIRP Annals—Manufacturing Technology, vol. 38, no. 1, pp. 557–560, 1989.

C. Guo and S. Malkin, “Analysis of energy partition in grinding,” Journal of Engineering for Industry, vol. 117, no. 1, pp. 55–61, 1995.

L. C. Zhang, T. Suto, T. H. Noguchi, and T. Waida, “An overview of applied mechanics in grinding,” Manufacturing Review, vol. 4, pp. 261–273, 1992.

M. Mahdi and L. Zhang, “Applied mechanics in grinding—VI. Residual stresses and surface hardening by coupled thermo-plasticity and phase transformation,” International Journal of Machine Tools and Manufacture, vol. 38, no. 10-11, pp. 1289–1304, 1998.

L. C. Zhang, T. Suto, H. Noguchi, and T. Waida, “Applied mechanics in grinding. Part II: modelling of elastic modulus of wheels and interface forces,” International Journal of Machine Tools and Manufacture, vol. 33, no. 2, pp. 245–255, 1993.

C. Guo, Y. Wu, V. Varghese, and S. Malkin, “Temperatures and energy partition for grinding with vitrified CBN wheels,” CIRP Annals, vol. 48, no. 1, pp. 247–250, 1999.

I. Zarudi and L. C. Zhang, “A revisit to some wheel—workpiece interaction problems in surface grinding,” International Journal of Machine Tools and Manufacture, vol. 42, no. 8, pp. 905–913, 2002.

W. B. Rowe, S. C. E. Black, B. Mills, H. S. Qi, and M. N. Morgan, “Experimental investigation of heat transfer in grinding,” CIRP Annals—Manufacturing Technology, vol. 44, no. 1, pp. 329–332, 1995.

J. González-Santander, “Analytic solution for maximum temperature during cut in and cut out in surface dry grinding,” Applied Mathematical Modelling, vol. 40, no. 3, pp. 2356–2367, 2016.

J. L. González-Santander, J. M. Isidro, and G. Martín, “An analysis of the transient regime temperature field in wet grinding,” Journal of Engineering Mathematics, vol. 90, pp. 141–171, 2015.

K. Oldham, J. Myland, and J. Spanier, An Atlas of Functions, Springer, New York, NY, USA, 2nd edition, 2008.

M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, USA, 1972.

J. L. González-Santander and G. Martín González, “Closed form expression for the surface temperature in wet grinding. Application to maximum temperature evaluation,” Journal of Engineering Mathematics, vol. 90, pp. 173–193, 2015.

J. L. González-Santander and G. Martín, “A theorem for finding maximum temperature in wet grinding,” Mathematical Problems in Engineering, vol. 2015, Article ID 150493, 13 pages, 2015.

R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, “On the Lambert W function,” Advances in Computational Mathematics, vol. 5, no. 4, pp. 329–359, 1996.

V. I. Murav'ev, A. V. Yakimov, and A. V. Chernyshev, “Effect of deformation, welding, and electrocontact heating on the properties of titanium alloy VT20 in pressed and welded structures,” Metal Science and Heat Treatment, vol. 45, no. 11-12, pp. 419–422, 2003.

D. L. Skuratov, Y. L. Ratis, I. A. Selezneva, J. Pérez, P. Fernández de Córdoba, and J. F. Urchueguía, “Mathematical modelling and analytical solution for workpiece temperature in grinding,” Applied Mathematical Modelling, vol. 31, no. 6, pp. 1039–1047, 2007.

A. S. Lavine, “A simple model for convective cooling during the grinding process,” Journal of Engineering for Industry, vol. 110, no. 1, pp. 1–6, 1988.

R. P. Brent, “An algorithm with guaranteed convergence for finding a zero of a function,” in Algorithms for Minimization without Derivatives, chapter 4, Prentice-Hall, Englewood Cliffs, NJ, USA, 1973.

H. Takahasi and M. Mori, “Double exponential formulas for numerical integration,” Publications of the Research Institute for Mathematical Sciences, vol. 9, pp. 721–741, 1973.