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

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 themost commonheat flux profiles reported in the literature, that is, constant, linear, parabolic, and triangular. For these cases, we propose a fastmethod for the numerical computation ofmaximum 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).


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 2ℓ width (m in SI units), which moves at a speed ⃗ V  = V  ⃗  (m s −1 ) 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  0 (K) must satisfy the convective heat equation [3, §1.7(2)]: where  is the thermal diffusivity (m 2 s −1 ).Since initially the workpiece is at room temperature, (1) is subjected to an homogeneous initial condition: Wet grinding can be modeled assuming a constant heat transfer coefficient ℎ (W m −2 K −1 ) 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,  0 (W m −1 K −1 ) is the thermal conductivity, and  (W m −2 ) is the average heat In the case of dry grinding, that is, ℎ = 0, 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, lim →∞ (, , 0), can be found [6].
For the wet case, that is, ℎ > 0, (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 ℎ ≥ 0, 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][11][12][13][14][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.

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  = 0, both coordinates systems are overlapped.
Green's function (, ⃗   ;  0 , ⃗   0 ) is interpreted in this case as the temperature field rise evaluated at (, ⃗   ) = (,   ,   ,   ) in a semi-infinite body ( > 0) in which, at position ⃗   0 = (  0 ,   0 ,   0 ), an instantaneous point heat source of energy  (J) appears at the instant  0 <  and where the surface ( = 0) is subject to a constant heat transfer coefficient ℎ [3, §14.9 (4)]: where  is the density of the body (kg m −3 ) and  its specific heat (J kg −1 K −1 ) and where we have defined the radiation coefficient (m −1 ) as A formal derivation of ( 5) is given in [18,Appendix].According to Figure 2, a point  of the heat source has ( 0 ,  0 , 0) coordinates in  and (  0 ,   0 , 0) 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  = 0 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.

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 −1 ) is constant: thus, integrating in  0 ∈ (−∞, ∞), we obtain the field temperature due to an instantaneous line heat source as

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 −2 ) is the energy that this heat source leaves at the instant  0 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,  0 ∈ (−ℓ, ℓ), 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

Continuous Strip Source.
Let  (W m −2 ) 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  0 , we have Substitution of ( 18) in ( 17) and integration over time where we have performed the change of variables   =  −  0 .Considering now the dimensionless variables where  is a characteristic length we finally arrive at where ) . (23)

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.

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  = 0, so, from (23) and using the property given in Appendix (A.12), we have Therefore, ( 22) reduces to where the superscript (0) denotes that we are considering a constant heat flux profile.Particularizing (26) to  = 0 and performing the change of variables  = 2 2 , 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  0 () 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 ( = 0) can be expressed in closed form as where It is worth noting that (31) does not converge for high Biot numbers, that is,  > 1.

Linear Heat Flux Profile.
For a linear heat flux profile satisfying (4), we have thus, ( 22) is reduced to where the superscript (1) 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  = 2 2 and the limit  → ∞, we obtain the following expression for the stationary regime, which is given in [7,: where the following function has been defined:

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: 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 Mathematical Problems in Engineering 7 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  0 ∈ (0, ), and performing the change of variables   =  −  0 , we get the time-dependent temperature field due to a triangular heat flux profile as which, in dimensionless variables (20), is rewritten as (45)

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 (2) 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

Maximum Temperature
From a physical point of view, the maximum temperature T max 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  max must be on the surface  = 0, 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 T(, , ) 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  * =  max .
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 (54)

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  max .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  ∈ (0, ∞): 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).

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  = 1 and taking  = , we have On the other hand, expanding in Taylor series the following function up to the first order, we have Substituting ( 59) and ( 61) in (57), we get the following approximated equation: By using the Lambert W function [23], (62) can be solved explicitly, obtaining the following approximated expression for the relaxation time: ).
Notice that (63) does not depend on  max , so we do not have to compute (57) if we want an estimation of the relaxation time for a constant heat flux profile.

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: Now, substituting (59), (61), and (66) in (64), we arrive at that can be solved explicitly as ). (68)

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  max ; moreover, it does not depend on  either.

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

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, Al 2 O 3 (sapphire), as workpiece material, respectively.The thermal properties and the grinding regimes for these data sets can be found in [26].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  = 0.7 for the triangular profile (45).
Table 2 shows the value of maximum temperature  max and its location  max /ℓ 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, lim →∞ T(, , 0).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.Figures 4-6 show the surface temperature ( = 0) 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 (/ℓ = ±1).This is because the heat source velocity V  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.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).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,  > 2ℓ/V  .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.

Conclusions
By using the Green function, we have solved the boundaryvalue 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.

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

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

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

Table 1 :
Simulation parameters in SI units.

Table 2 :
Location and value of maximum temperature.