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

Authors: J. L. González-Santander
Journal: Mathematical Problems in Engineering
 (2016)

---

## 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).


---

## Body

## 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 of2
l width (m in SI units), which moves at a speed v
→
f
=
v
f
i
→ (m s−1) over a semi-infinite solid surface (see Figure 1). Setting the Cartesian coordinate system X
Y
Z fixed to wheel, as it is shown in Figure 1, the temperature field of the workpiece T
t
,
x
,
z with respect to the room temperature T
0 (K) must satisfy the convective heat equation [3, §1.7(2)]:(1)
∂
T
∂
t
=
k
∂
2
T
∂
x
2
+
∂
2
T
∂
z
2
-
v
f
∂
T
∂
x
,where k is the thermal diffusivity (m2 s−1). Since initially the workpiece is at room temperature, (1) is subjected to an homogeneous initial condition:(2)
T
0
,
x
,
z
=
0
.Figure 1
Setup in surface down-grinding.Wet grinding can be modeled assuming a constant heat transfer coefficienth (W m−2 K−1) on the workpiece surface [4]. Assuming as well a dimensionless heat flux profile f
x within the contact area between wheel and workpiece, we have the following boundary condition:(3)
k
0
∂
T
∂
z
t
,
x
,
0
=
h
T
t
,
x
,
0
-
q
f
x
θ
l
-
x
θ
l
+
x
,where θ
x denotes the Heaviside function, k
0 (W m−1 K−1) is the thermal conductivity, and q (W m−2) is the average heat flux entering into the workpiece along the contact width 2
l; thus,(4)
1
2
l
∫
-
l
l
f
x
d
x
=
1
.In the case of dry grinding, that is,h
=
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, l
i
m
t
→
∞
T
t
,
x
,
0, can be found [6].For the wet case, that is,h
>
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 T
t
,
x
,
z satisfying (1)–(3) for any constant h
≥
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–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. Section2 derives the solution of the value boundary problem stated in (1)–(3) for any analytic heat flux profile f
x. 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 f
x 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 X
′
Y
′
Z
′ (see Figure 2), where, at t
=
0, 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 functionG
t
,
r
→
′
;
t
0
,
r
→
0
′ is interpreted in this case as the temperature field rise evaluated at t
,
r
→
′
=
t
,
x
′
,
y
′
,
z
′ in a semi-infinite body (z
>
0) in which, at position r
→
0
′
=
x
0
′
,
y
0
′
,
z
0
′, an instantaneous point heat source of energy E (J) appears at the instant t
0
<
t and where the surface (z
=
0) is subject to a constant heat transfer coefficient h [3, §14.9 (4)]:(5)
G
t
,
r
→
′
;
t
0
,
r
→
0
′
=
E
4
π
ρ
c
k
t
-
t
0
exp
⁡
-
x
′
-
x
0
′
2
+
y
′
-
y
0
′
2
4
k
t
-
t
0
exp
⁡
-
z
′
-
z
0
′
2
/
4
k
t
-
t
0
+
exp
⁡
-
z
′
+
z
0
′
2
/
4
k
t
-
t
0
2
π
k
t
-
t
0
-
h
~
exp
⁡
h
~
z
′
+
z
0
′
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
′
+
z
0
′
2
k
t
-
t
0
,where ρ is the density of the body (kg m−3) and c its specific heat (J kg−1 K−1) and where we have defined theradiation coefficient (m−1) as(6)
h
~
=
h
k
0
.A formal derivation of (5) is given in [18, Appendix].According to Figure2, a point P of the heat source has x
0
,
y
0
,
0 coordinates in X
Y
Z and x
0
′
,
y
0
′
,
0 in X
′
Y
′
Z
′. The relationship between both coordinates systems is (7)
x
′
=
x
+
v
f
t
-
t
0
,
y
′
=
y
,
z
′
=
z
.Therefore, the temperature rise in the workpiece at coordinatest
,
x
,
y
,
z due to a point P of the heat source with respect to the coordinate system fixed to the wheel X
Y
Z is given by(8)
G
t
,
x
,
y
,
z
=
E
4
π
k
0
t
-
t
0
exp
⁡
-
x
+
v
f
t
-
t
0
-
x
0
2
+
y
-
y
0
2
4
k
t
-
t
0
e
x
p
-
z
2
/
4
k
t
-
t
0
2
π
k
t
-
t
0
-
h
~
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
.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 surfacez
=
0 parallel to the y-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 they-axis. Let us assume in (8) that the energy per unit length of this heat source q
L (J m−1) is constant:(9)
d
E
=
q
L
d
y
0
;thus, integrating in y
0
∈
-
∞
,
∞, we obtain the field temperature due to an instantaneous line heat source as(10)
T
l
i
n
e
t
,
x
,
z
=
q
L
2
π
ρ
c
k
t
-
t
0
exp
⁡
-
x
+
v
f
t
-
t
0
-
x
0
2
4
k
t
-
t
0
exp
⁡
-
z
2
4
k
t
-
t
0
-
h
~
π
k
t
-
t
0
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
.


### 2.2. Instantaneous Strip Source

Let us consider now a strip source of variable heat strength along thex-axis, according to the following profile:(11)
q
L
=
q
B
f
x
0
d
x
0
,where q
B (J m−2) is the energy that this heat source leaves at the instant t
0 per surface unit of the strip source. Notice that the f function in (11) is dimensionless, as in (3). Also, f 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, x
0
∈
-
l
,
l, we have the following temperature field in the workpiece:(12)
T
s
t
r
i
p
t
,
x
,
z
=
q
B
2
π
ρ
c
k
t
-
t
0
exp
⁡
-
z
2
4
k
t
-
t
0
-
h
~
π
k
t
-
t
0
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
∫
-
l
l
f
x
0
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
.Assuming that thef function is analytic, we can expand it in its Taylor series, so that, taking into account also that f is dimensionless, we have(13)
f
x
0
=
g
x
0
l
=
∑
m
=
0
∞
g
m
0
m
!
x
0
l
m
;thus, the integral given in (12) can be written as (14)
∫
-
l
l
f
x
0
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
=
∑
m
=
0
∞
g
m
0
l
m
m
!
∫
-
l
l
x
0
m
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
.According to the result given in Appendix (A.9), we have (15)
∫
-
l
l
x
0
m
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
=
-
2
k
t
-
t
0
m
+
1
2
∑
n
=
0
m
m
n
-
x
+
v
f
t
-
t
0
2
k
t
-
t
0
m
-
n
F
n
u
,
t
-
t
0
u
=
x
-
l
x
+
l
, where(16)
F
n
u
,
t
=
sign
n
+
1
⁡
u
+
v
f
t
2
k
t
γ
n
+
1
2
,
u
+
v
f
t
2
k
t
2and where γ
a
,
z denotes thelower incomplete gamma function [19, Eqn. 45:1:2]. Substituting back these results, we arrive at(17)
T
s
t
r
i
p
t
,
x
,
z
=
q
B
2
π
ρ
c
k
t
-
t
0
exp
⁡
-
z
2
4
k
t
-
t
0
-
h
~
π
k
t
-
t
0
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
∑
m
=
0
∞
g
m
0
x
+
v
f
t
-
t
0
l
m
∑
n
=
0
m
1
n
!
m
-
n
!
-
2
k
t
-
t
0
x
+
v
f
t
-
t
0
n
F
n
u
,
t
-
t
0
u
=
x
-
l
x
+
l
.


### 2.3. Continuous Strip Source

Letq (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 d
t
0, we have(18)
q
B
=
q
d
t
0
.Substitution of (18) in (17) and integration over time t
0
∈
0
,
t (remember that t
0
<
t) result in (19)
T
t
,
x
,
z
=
q
2
π
ρ
c
∫
0
t
exp
⁡
-
z
2
/
4
k
t
′
k
t
′
-
h
~
π
k
t
′
e
h
~
z
+
h
~
2
k
t
′
erfc
⁡
h
~
k
t
′
+
z
2
k
t
′
∑
m
=
0
∞
g
m
0
x
+
v
f
t
′
l
m
∑
n
=
0
m
1
n
!
m
-
n
!
-
2
k
t
′
x
+
v
f
t
′
n
F
n
u
,
t
′
u
=
x
-
l
x
+
l
d
t
′
,where we have performed the change of variables t
′
=
t
-
t
0. Considering now the dimensionless variables(20)
X
=
x
s
,
Z
=
z
s
,
L
=
l
s
,
H
=
h
~
s
,
τ
=
v
f
t
2
k
,
T
=
π
ρ
c
v
f
2
q
T
,where s is a characteristic length(21)
s
=
2
k
v
f
,we finally arrive at(22)
T
τ
,
X
,
Z
=
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
∑
m
=
0
∞
g
m
0
L
m
∑
n
=
0
m
-
2
w
n
X
+
2
w
2
m
-
n
n
!
m
-
n
!
F
n
u
,
w
u
=
X
-
L
X
+
L
d
w
, where(23)
F
n
u
,
w
=
sign
n
+
1
⁡
u
2
w
+
w
γ
n
+
1
2
,
u
2
w
+
w
2
.


## 2.1. Instantaneous Line Source

Let us consider now an instantaneous line heat source parallel to they-axis. Let us assume in (8) that the energy per unit length of this heat source q
L (J m−1) is constant:(9)
d
E
=
q
L
d
y
0
;thus, integrating in y
0
∈
-
∞
,
∞, we obtain the field temperature due to an instantaneous line heat source as(10)
T
l
i
n
e
t
,
x
,
z
=
q
L
2
π
ρ
c
k
t
-
t
0
exp
⁡
-
x
+
v
f
t
-
t
0
-
x
0
2
4
k
t
-
t
0
exp
⁡
-
z
2
4
k
t
-
t
0
-
h
~
π
k
t
-
t
0
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
.


## 2.2. Instantaneous Strip Source

Let us consider now a strip source of variable heat strength along thex-axis, according to the following profile:(11)
q
L
=
q
B
f
x
0
d
x
0
,where q
B (J m−2) is the energy that this heat source leaves at the instant t
0 per surface unit of the strip source. Notice that the f function in (11) is dimensionless, as in (3). Also, f 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, x
0
∈
-
l
,
l, we have the following temperature field in the workpiece:(12)
T
s
t
r
i
p
t
,
x
,
z
=
q
B
2
π
ρ
c
k
t
-
t
0
exp
⁡
-
z
2
4
k
t
-
t
0
-
h
~
π
k
t
-
t
0
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
∫
-
l
l
f
x
0
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
.Assuming that thef function is analytic, we can expand it in its Taylor series, so that, taking into account also that f is dimensionless, we have(13)
f
x
0
=
g
x
0
l
=
∑
m
=
0
∞
g
m
0
m
!
x
0
l
m
;thus, the integral given in (12) can be written as (14)
∫
-
l
l
f
x
0
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
=
∑
m
=
0
∞
g
m
0
l
m
m
!
∫
-
l
l
x
0
m
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
.According to the result given in Appendix (A.9), we have (15)
∫
-
l
l
x
0
m
exp
⁡
-
x
-
x
0
+
v
f
t
-
t
0
2
k
t
-
t
0
2
d
x
0
=
-
2
k
t
-
t
0
m
+
1
2
∑
n
=
0
m
m
n
-
x
+
v
f
t
-
t
0
2
k
t
-
t
0
m
-
n
F
n
u
,
t
-
t
0
u
=
x
-
l
x
+
l
, where(16)
F
n
u
,
t
=
sign
n
+
1
⁡
u
+
v
f
t
2
k
t
γ
n
+
1
2
,
u
+
v
f
t
2
k
t
2and where γ
a
,
z denotes thelower incomplete gamma function [19, Eqn. 45:1:2]. Substituting back these results, we arrive at(17)
T
s
t
r
i
p
t
,
x
,
z
=
q
B
2
π
ρ
c
k
t
-
t
0
exp
⁡
-
z
2
4
k
t
-
t
0
-
h
~
π
k
t
-
t
0
e
h
~
z
+
h
~
2
k
t
-
t
0
erfc
⁡
h
~
k
t
-
t
0
+
z
2
k
t
-
t
0
∑
m
=
0
∞
g
m
0
x
+
v
f
t
-
t
0
l
m
∑
n
=
0
m
1
n
!
m
-
n
!
-
2
k
t
-
t
0
x
+
v
f
t
-
t
0
n
F
n
u
,
t
-
t
0
u
=
x
-
l
x
+
l
.


## 2.3. Continuous Strip Source

Letq (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 d
t
0, we have(18)
q
B
=
q
d
t
0
.Substitution of (18) in (17) and integration over time t
0
∈
0
,
t (remember that t
0
<
t) result in (19)
T
t
,
x
,
z
=
q
2
π
ρ
c
∫
0
t
exp
⁡
-
z
2
/
4
k
t
′
k
t
′
-
h
~
π
k
t
′
e
h
~
z
+
h
~
2
k
t
′
erfc
⁡
h
~
k
t
′
+
z
2
k
t
′
∑
m
=
0
∞
g
m
0
x
+
v
f
t
′
l
m
∑
n
=
0
m
1
n
!
m
-
n
!
-
2
k
t
′
x
+
v
f
t
′
n
F
n
u
,
t
′
u
=
x
-
l
x
+
l
d
t
′
,where we have performed the change of variables t
′
=
t
-
t
0. Considering now the dimensionless variables(20)
X
=
x
s
,
Z
=
z
s
,
L
=
l
s
,
H
=
h
~
s
,
τ
=
v
f
t
2
k
,
T
=
π
ρ
c
v
f
2
q
T
,where s is a characteristic length(21)
s
=
2
k
v
f
,we finally arrive at(22)
T
τ
,
X
,
Z
=
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
∑
m
=
0
∞
g
m
0
L
m
∑
n
=
0
m
-
2
w
n
X
+
2
w
2
m
-
n
n
!
m
-
n
!
F
n
u
,
w
u
=
X
-
L
X
+
L
d
w
, where(23)
F
n
u
,
w
=
sign
n
+
1
⁡
u
2
w
+
w
γ
n
+
1
2
,
u
2
w
+
w
2
.


## 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 (24)
f
0
x
=
g
0
x
l
=
1
. Thus, in (22), we have only the summand which corresponds to m
=
0, so, from (23) and using the property given in Appendix (A.12), we have(25)
F
0
u
,
w
=
sign
⁡
u
2
w
+
w
γ
1
2
,
u
2
w
+
w
2
=
π
erf
⁡
u
2
w
+
w
.Therefore, (22) reduces to(26)
T
0
τ
,
X
,
Z
=
π
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
erf
⁡
X
+
L
2
w
+
w
-
erf
⁡
X
-
L
2
w
+
w
d
w
,where the superscript (
0
) denotes that we are considering a constant heat flux profile. Particularizing (26) to H
=
0 and performing the change of variables s
=
2
w
2, we obtain the following expression for the case of dry grinding:(27)
T
d
r
y
0
τ
,
X
,
Z
=
1
2
π
2
∫
0
2
τ
2
exp
⁡
-
Z
2
4
s
erf
⁡
X
+
L
+
s
2
s
-
erf
⁡
X
-
L
+
s
2
s
d
s
s
,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]:(28)
lim
τ
→
∞
⁡
T
d
r
y
0
τ
,
X
,
Z
=
π
∫
0
∞
exp
⁡
-
Z
2
4
w
2
erf
⁡
X
+
L
2
w
+
w
-
erf
⁡
X
-
L
2
w
+
w
d
w
(29)
=
∫
X
-
L
X
+
L
e
-
w
K
0
w
2
+
Z
2
d
w
,where K
0
x is the zeroth-ordermodified 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(30)
lim
τ
→
∞
⁡
T
0
τ
,
X
,
Z
=
∫
X
-
L
X
+
L
e
-
w
K
0
w
2
+
Z
2
d
w
-
π
H
e
H
Z
∫
0
∞
w
e
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
erf
⁡
X
+
L
2
w
+
w
-
erf
⁡
X
-
L
2
w
+
w
d
w
, which is the expression given by [4]. Recently, in [21], the expression given in (30) particularized on the surface (Z
=
0) can be expressed in closed form as(31)
lim
τ
→
∞
⁡
T
0
τ
,
X
,
0
=
π
∑
n
=
0
∞
-
H
/
2
n
Γ
n
+
1
/
2
Yu
n
/
2
⁡
u
X
-
L
X
+
L
,where(32)
Yu
ν
⁡
x
=
x
x
ν
e
x
2
ν
+
1
K
ν
x
+
sign
⁡
x
K
ν
+
1
x
,
x
≠
0
,
0
,
x
=
0
.It is worth noting that (31) does not converge for high Biot numbers, that is, H
>
1.


### 3.2. Linear Heat Flux Profile

For a linear heat flux profile satisfying (4), we have (33)
f
1
x
=
g
x
l
=
1
+
x
l
; thus, (22) is reduced to(34)
T
1
τ
,
X
,
Z
=
1
L
∫
0
τ
d
w
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
π
X
+
L
+
2
w
2
e
r
f
u
2
w
+
w
+
2
w
exp
⁡
-
u
2
w
+
w
2
u
=
X
-
L
X
+
L
, 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: (35)
F
1
u
,
w
u
=
X
-
L
X
+
L
=
γ
1
,
u
2
w
+
w
2
u
=
X
-
L
X
+
L
=
-
exp
⁡
-
u
2
w
+
w
2
u
=
X
-
L
X
+
L
.Performing in (34) the change of variables r
=
2
w
2 and the limit τ
→
∞, we obtain the following expression for the stationary regime, which is given in [7, Eqns.3.46-47]:(36)
lim
τ
→
∞
⁡
T
1
τ
,
X
,
Z
=
1
2
L
∫
0
∞
exp
⁡
-
Z
2
2
r
π
2
u
X
+
L
+
r
erfc
⁡
X
+
L
+
r
2
r
-
erfc
⁡
X
-
L
+
r
2
r
+
exp
⁡
-
X
+
L
+
r
2
2
u
-
exp
⁡
-
X
-
L
+
r
2
2
r
d
r
-
π
H
4
L
∫
0
∞
exp
⁡
-
Z
2
2
r
experf
⁡
H
r
2
r
X
+
L
+
r
erfc
⁡
X
-
L
+
r
2
r
-
erfc
⁡
X
+
L
+
r
2
r
+
2
r
π
exp
⁡
-
X
+
L
+
r
2
2
r
-
exp
⁡
-
X
-
L
+
r
2
2
r
d
r
, where the following function has been defined: (37)
experf
⁡
x
=
exp
⁡
x
2
erfc
⁡
x
.


### 3.3. Triangular Heat Flux Profile

For a triangular heat flux profile satisfying (4) (see Figure 3), we have(38)
f
Δ
x
=
g
Δ
x
l
=
2
λ
+
1
x
l
+
1
,
x
∈
-
l
,
l
m
a
x
2
λ
-
1
x
l
-
1
,
x
∈
l
m
a
x
,
l
,where we have defined the following dimensionless parameter: (39)
λ
=
l
m
a
x
l
∈
-
1,1
.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: (40)
∫
-
l
l
f
x
0
exp
⁡
-
A
+
B
x
0
2
d
x
0
=
2
λ
+
1
∫
-
l
l
m
a
x
x
0
l
+
1
exp
⁡
-
A
+
B
x
0
2
d
x
0
+
2
λ
-
1
∫
l
m
a
x
l
x
0
l
-
1
exp
⁡
-
A
+
B
x
0
2
d
x
0
, where we have set (41)
A
=
x
+
v
f
t
-
t
0
2
k
t
-
t
0
,
B
=
-
1
2
k
t
-
t
0
.By using formulas (A.13) and (A.15) given in Appendix, some algebra leads to the following result:(42)
∫
-
l
l
f
x
0
exp
⁡
-
A
+
B
x
0
2
d
x
0
=
1
l
B
2
φ
A
-
B
l
λ
+
1
+
2
φ
A
+
B
l
m
a
x
λ
2
-
1
-
φ
A
+
B
l
λ
-
1
, where we have defined(43)
φ
x
=
π
x
erf
⁡
x
+
exp
⁡
-
x
2
.Substituting (42) into (12), integrating over time t
0
∈
0
,
t, and performing the change of variables t
′
=
t
-
t
0, we get the time-dependent temperature field due to a triangular heat flux profile as (44)
T
Δ
t
,
x
,
z
=
2
q
π
ρ
c
∫
0
t
exp
⁡
-
z
2
/
4
k
t
′
k
t
′
-
π
h
~
e
h
~
z
+
h
~
2
k
t
′
erfc
⁡
h
~
k
t
′
+
z
2
k
t
′
k
t
′
l
φ
x
+
l
+
v
f
t
′
/
2
k
t
′
λ
+
1
+
2
φ
x
-
λ
l
+
v
f
t
′
/
2
k
t
′
λ
2
-
1
-
φ
x
-
l
+
v
f
t
′
/
2
k
t
′
λ
-
1
d
t
′
, which, in dimensionless variables (20), is rewritten as(45)
T
Δ
τ
,
X
,
Z
=
4
L
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
w
φ
X
+
L
/
2
w
+
w
λ
+
1
+
2
φ
X
-
λ
L
/
2
w
+
w
λ
2
-
1
-
φ
X
-
L
/
2
w
+
w
λ
-
1
d
w
.


### 3.4. Parabolic Heat Flux Profile

For a parabolic heat flux profile satisfying (4) (see Figure 3), we have(46)
f
2
x
=
g
2
x
l
=
3
4
1
+
x
l
2
;thus, (22) is reduced to (47)
T
2
τ
,
X
,
Z
=
3
2
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
1
2
+
X
+
2
w
2
L
+
X
+
2
w
2
2
2
L
2
F
0
u
,
w
-
2
w
L
1
+
X
+
2
w
2
L
F
1
u
,
w
+
2
w
2
L
2
F
2
u
,
w
u
=
X
-
L
X
+
L
d
w
, 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(48)
T
2
τ
,
X
,
Z
=
3
4
L
2
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
π
2
w
2
+
L
+
2
w
2
+
X
2
erf
⁡
u
2
w
+
w
+
2
w
2
L
+
2
X
+
2
w
2
-
u
exp
⁡
-
u
2
w
+
w
2
u
=
X
-
L
X
+
L
d
w
.


## 3.1. Constant Heat Flux Profile

For a constant heat flux profile satisfying (4), we have (24)
f
0
x
=
g
0
x
l
=
1
. Thus, in (22), we have only the summand which corresponds to m
=
0, so, from (23) and using the property given in Appendix (A.12), we have(25)
F
0
u
,
w
=
sign
⁡
u
2
w
+
w
γ
1
2
,
u
2
w
+
w
2
=
π
erf
⁡
u
2
w
+
w
.Therefore, (22) reduces to(26)
T
0
τ
,
X
,
Z
=
π
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
erf
⁡
X
+
L
2
w
+
w
-
erf
⁡
X
-
L
2
w
+
w
d
w
,where the superscript (
0
) denotes that we are considering a constant heat flux profile. Particularizing (26) to H
=
0 and performing the change of variables s
=
2
w
2, we obtain the following expression for the case of dry grinding:(27)
T
d
r
y
0
τ
,
X
,
Z
=
1
2
π
2
∫
0
2
τ
2
exp
⁡
-
Z
2
4
s
erf
⁡
X
+
L
+
s
2
s
-
erf
⁡
X
-
L
+
s
2
s
d
s
s
,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]:(28)
lim
τ
→
∞
⁡
T
d
r
y
0
τ
,
X
,
Z
=
π
∫
0
∞
exp
⁡
-
Z
2
4
w
2
erf
⁡
X
+
L
2
w
+
w
-
erf
⁡
X
-
L
2
w
+
w
d
w
(29)
=
∫
X
-
L
X
+
L
e
-
w
K
0
w
2
+
Z
2
d
w
,where K
0
x is the zeroth-ordermodified 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(30)
lim
τ
→
∞
⁡
T
0
τ
,
X
,
Z
=
∫
X
-
L
X
+
L
e
-
w
K
0
w
2
+
Z
2
d
w
-
π
H
e
H
Z
∫
0
∞
w
e
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
erf
⁡
X
+
L
2
w
+
w
-
erf
⁡
X
-
L
2
w
+
w
d
w
, which is the expression given by [4]. Recently, in [21], the expression given in (30) particularized on the surface (Z
=
0) can be expressed in closed form as(31)
lim
τ
→
∞
⁡
T
0
τ
,
X
,
0
=
π
∑
n
=
0
∞
-
H
/
2
n
Γ
n
+
1
/
2
Yu
n
/
2
⁡
u
X
-
L
X
+
L
,where(32)
Yu
ν
⁡
x
=
x
x
ν
e
x
2
ν
+
1
K
ν
x
+
sign
⁡
x
K
ν
+
1
x
,
x
≠
0
,
0
,
x
=
0
.It is worth noting that (31) does not converge for high Biot numbers, that is, H
>
1.


## 3.2. Linear Heat Flux Profile

For a linear heat flux profile satisfying (4), we have (33)
f
1
x
=
g
x
l
=
1
+
x
l
; thus, (22) is reduced to(34)
T
1
τ
,
X
,
Z
=
1
L
∫
0
τ
d
w
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
π
X
+
L
+
2
w
2
e
r
f
u
2
w
+
w
+
2
w
exp
⁡
-
u
2
w
+
w
2
u
=
X
-
L
X
+
L
, 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: (35)
F
1
u
,
w
u
=
X
-
L
X
+
L
=
γ
1
,
u
2
w
+
w
2
u
=
X
-
L
X
+
L
=
-
exp
⁡
-
u
2
w
+
w
2
u
=
X
-
L
X
+
L
.Performing in (34) the change of variables r
=
2
w
2 and the limit τ
→
∞, we obtain the following expression for the stationary regime, which is given in [7, Eqns.3.46-47]:(36)
lim
τ
→
∞
⁡
T
1
τ
,
X
,
Z
=
1
2
L
∫
0
∞
exp
⁡
-
Z
2
2
r
π
2
u
X
+
L
+
r
erfc
⁡
X
+
L
+
r
2
r
-
erfc
⁡
X
-
L
+
r
2
r
+
exp
⁡
-
X
+
L
+
r
2
2
u
-
exp
⁡
-
X
-
L
+
r
2
2
r
d
r
-
π
H
4
L
∫
0
∞
exp
⁡
-
Z
2
2
r
experf
⁡
H
r
2
r
X
+
L
+
r
erfc
⁡
X
-
L
+
r
2
r
-
erfc
⁡
X
+
L
+
r
2
r
+
2
r
π
exp
⁡
-
X
+
L
+
r
2
2
r
-
exp
⁡
-
X
-
L
+
r
2
2
r
d
r
, where the following function has been defined: (37)
experf
⁡
x
=
exp
⁡
x
2
erfc
⁡
x
.


## 3.3. Triangular Heat Flux Profile

For a triangular heat flux profile satisfying (4) (see Figure 3), we have(38)
f
Δ
x
=
g
Δ
x
l
=
2
λ
+
1
x
l
+
1
,
x
∈
-
l
,
l
m
a
x
2
λ
-
1
x
l
-
1
,
x
∈
l
m
a
x
,
l
,where we have defined the following dimensionless parameter: (39)
λ
=
l
m
a
x
l
∈
-
1,1
.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: (40)
∫
-
l
l
f
x
0
exp
⁡
-
A
+
B
x
0
2
d
x
0
=
2
λ
+
1
∫
-
l
l
m
a
x
x
0
l
+
1
exp
⁡
-
A
+
B
x
0
2
d
x
0
+
2
λ
-
1
∫
l
m
a
x
l
x
0
l
-
1
exp
⁡
-
A
+
B
x
0
2
d
x
0
, where we have set (41)
A
=
x
+
v
f
t
-
t
0
2
k
t
-
t
0
,
B
=
-
1
2
k
t
-
t
0
.By using formulas (A.13) and (A.15) given in Appendix, some algebra leads to the following result:(42)
∫
-
l
l
f
x
0
exp
⁡
-
A
+
B
x
0
2
d
x
0
=
1
l
B
2
φ
A
-
B
l
λ
+
1
+
2
φ
A
+
B
l
m
a
x
λ
2
-
1
-
φ
A
+
B
l
λ
-
1
, where we have defined(43)
φ
x
=
π
x
erf
⁡
x
+
exp
⁡
-
x
2
.Substituting (42) into (12), integrating over time t
0
∈
0
,
t, and performing the change of variables t
′
=
t
-
t
0, we get the time-dependent temperature field due to a triangular heat flux profile as (44)
T
Δ
t
,
x
,
z
=
2
q
π
ρ
c
∫
0
t
exp
⁡
-
z
2
/
4
k
t
′
k
t
′
-
π
h
~
e
h
~
z
+
h
~
2
k
t
′
erfc
⁡
h
~
k
t
′
+
z
2
k
t
′
k
t
′
l
φ
x
+
l
+
v
f
t
′
/
2
k
t
′
λ
+
1
+
2
φ
x
-
λ
l
+
v
f
t
′
/
2
k
t
′
λ
2
-
1
-
φ
x
-
l
+
v
f
t
′
/
2
k
t
′
λ
-
1
d
t
′
, which, in dimensionless variables (20), is rewritten as(45)
T
Δ
τ
,
X
,
Z
=
4
L
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
w
φ
X
+
L
/
2
w
+
w
λ
+
1
+
2
φ
X
-
λ
L
/
2
w
+
w
λ
2
-
1
-
φ
X
-
L
/
2
w
+
w
λ
-
1
d
w
.


## 3.4. Parabolic Heat Flux Profile

For a parabolic heat flux profile satisfying (4) (see Figure 3), we have(46)
f
2
x
=
g
2
x
l
=
3
4
1
+
x
l
2
;thus, (22) is reduced to (47)
T
2
τ
,
X
,
Z
=
3
2
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
1
2
+
X
+
2
w
2
L
+
X
+
2
w
2
2
2
L
2
F
0
u
,
w
-
2
w
L
1
+
X
+
2
w
2
L
F
1
u
,
w
+
2
w
2
L
2
F
2
u
,
w
u
=
X
-
L
X
+
L
d
w
, 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(48)
T
2
τ
,
X
,
Z
=
3
4
L
2
∫
0
τ
exp
⁡
-
Z
2
4
w
2
-
π
H
w
e
H
Z
+
H
2
w
2
erfc
⁡
H
w
+
Z
2
w
π
2
w
2
+
L
+
2
w
2
+
X
2
erf
⁡
u
2
w
+
w
+
2
w
2
L
+
2
X
+
2
w
2
-
u
exp
⁡
-
u
2
w
+
w
2
u
=
X
-
L
X
+
L
d
w
.


## 4. Maximum Temperature

From a physical point of view, the maximum temperatureT
m
a
x 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 X
m
a
x must be on the surface Z
=
0, within the contact zone between wheel and workpiece, X
m
a
x
∈
-
L
,
L. It is worth noting that recently the latter has been mathematically proved in [22] for constant and linear heat flux profiles. Therefore, since T
τ
,
X
,
Z is a differentiable function in X, a convenient way to evaluate maximum temperature is to solve numerically all points X
∗ that satisfies(49)
lim
τ
→
∞
⁡
∂
T
τ
,
X
∗
,
0
∂
X
=
0
,
X
∗
∈
-
L
,
L
,and then the maximum is given by (50)
T
m
a
x
=
max
X
∗
⁡
lim
τ
→
∞
⁡
T
τ
,
X
∗
,
0
.In Section6, we will check numerically that for the most common heat flux profiles proposed in the literature (constant, linear, triangular, and parabolic) the root X
∗ is unique, so X
∗
=
X
m
a
x.For a constant heat flux profile, (49) reads(51)
0
=
∫
0
∞
1
w
-
π
H
e
H
2
w
2
erfc
⁡
H
w
exp
⁡
-
u
2
w
+
w
2
u
=
X
∗
-
L
X
∗
+
L
d
w
.Similarly, for a linear heat flux profile, we have(52)
0
=
∫
0
∞
1
-
π
H
w
e
H
2
w
2
erfc
⁡
H
w
π
X
∗
+
L
+
2
w
2
e
r
f
u
2
w
+
w
-
X
∗
+
L
-
u
w
exp
⁡
-
u
2
w
+
w
2
u
=
X
∗
-
L
X
∗
+
L
d
w
.For a triangular profile, we obtain(53)
0
=
∫
0
∞
1
-
π
H
w
e
H
2
w
2
erfc
⁡
H
w
erf
⁡
X
∗
+
L
/
2
w
+
w
λ
+
1
+
2
erf
⁡
X
∗
-
λ
L
/
2
w
+
w
λ
2
-
1
-
erf
⁡
X
∗
-
L
/
2
w
+
w
λ
-
1
d
w
,and for a parabolic one(54)
0
=
∫
0
∞
1
-
π
H
w
e
H
2
w
2
erfc
⁡
H
w
π
X
∗
+
L
+
2
w
2
e
r
f
u
2
w
+
w
+
X
∗
+
L
-
u
2
2
w
+
2
w
exp
⁡
-
u
2
w
+
w
2
u
=
X
∗
-
L
X
∗
+
L
d
w
.


## 5. Relaxation Time

The stationary regime is asymptotically reached atτ
→
∞; that is,(55)
lim
τ
→
∞
⁡
∂
T
τ
,
X
,
Z
∂
τ
=
0
.In order to avoid thermal damage, the most important point is the location of the maximum temperatureX
m
a
x. 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
,
∞:(56)
∂
T
τ
¯
,
X
m
a
x
,
0
∂
τ
=
η
≈
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).


### 5.1. Constant Heat Flux Profile

For a constant heat flux profile, according to (26), (56) is written as(57)
η
=
π
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
.Let us solve (57) approximately. On the one hand, we find in the literature [20, Eqn.7.1.23] the following asymptotic expansion: (58)
π
z
e
z
2
erfc
⁡
z
≈
1
+
∑
m
=
1
∞
-
1
m
2
m
-
1
!
!
2
z
2
m
,
z
⟶
∞
; thus, expanding up to m
=
1 and taking z
=
H
τ
¯, we have(59)
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
≈
1
2
H
2
τ
¯
2
,
τ
¯
⟶
∞
.On the other hand, expanding in Taylor series the following function up to the first order, we have(60)
erf
⁡
τ
¯
+
X
m
a
x
±
L
2
τ
¯
≈
erf
⁡
τ
¯
+
X
m
a
x
±
L
π
τ
¯
e
-
τ
¯
2
,
τ
¯
⟶
∞
; thus,(61)
π
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
≈
2
L
τ
¯
e
-
τ
¯
2
,
τ
¯
⟶
∞
.Substituting (59) and (61) in (57), we get the following approximated equation:(62)
τ
¯
3
e
τ
¯
2
≈
L
H
2
η
.By using the LambertW function [23], (62) can be solved explicitly, obtaining the following approximated expression for the relaxation time:(63)
τ
¯
≈
3
2
W
2
3
L
H
2
η
2
/
3
.Notice that (63) does not depend on X
m
a
x, 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(64)
η
L
=
π
X
m
a
x
+
L
+
2
τ
¯
2
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
-
2
τ
¯
exp
⁡
-
X
m
a
x
+
L
2
τ
¯
+
τ
¯
2
-
exp
⁡
-
X
m
a
x
-
L
2
τ
¯
+
τ
¯
2
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
.In order to solve (64) approximately, let us expand in Taylor series the following function up to first order:(65)
exp
⁡
-
τ
¯
+
X
m
a
x
±
L
2
τ
¯
2
≈
e
-
τ
¯
2
-
X
m
a
x
±
L
e
-
τ
¯
2
,
τ
¯
⟶
∞
; thus,(66)
exp
⁡
-
X
m
a
x
+
L
2
τ
¯
+
τ
¯
2
-
exp
⁡
-
X
m
a
x
-
L
2
τ
¯
+
τ
¯
2
≈
2
L
e
-
τ
¯
2
,
τ
¯
⟶
∞
.Now, substituting (59), (61), and (66) in (64), we arrive at (67)
τ
¯
3
e
τ
¯
2
≈
X
m
a
x
+
L
H
2
η that can be solved explicitly as(68)
τ
¯
≈
3
2
W
2
3
X
m
a
x
+
L
H
2
η
2
/
3
.


### 5.3. Triangular Heat Flux Profile

Let us consider now a triangular heat flux profile, so, according to (45), (56) reads as(69)
η
L
4
=
φ
X
m
a
x
+
L
/
2
τ
¯
+
τ
¯
λ
+
1
+
2
φ
X
m
a
x
-
λ
L
/
2
τ
¯
+
τ
¯
λ
2
-
1
-
φ
X
m
a
x
-
L
/
2
τ
¯
+
τ
¯
λ
-
1
τ
¯
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
.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):(70)
φ
χ
2
τ
¯
+
τ
¯
=
π
χ
2
τ
¯
+
τ
¯
erf
⁡
χ
2
τ
¯
+
τ
¯
+
exp
⁡
-
χ
2
τ
¯
+
τ
¯
2
≈
π
χ
2
τ
¯
+
τ
¯
erf
⁡
τ
¯
+
1
+
χ
2
2
τ
¯
2
e
-
τ
¯
2
,
τ
¯
⟶
∞
.Taking into account (59) and (70) in (69), after some algebra, we arrive at the following approximated equation: (71)
τ
¯
3
e
τ
¯
2
≈
2
L
H
2
η that can be solved explicitly as(72)
τ
¯
≈
3
2
W
2
3
2
L
H
2
η
2
/
3
.As in the constant case, (72) does not depend on X
m
a
x; 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(73)
4
L
2
3
η
=
2
τ
¯
X
m
a
x
+
L
+
2
τ
¯
2
exp
⁡
-
X
m
a
x
+
L
2
τ
¯
+
τ
¯
2
-
X
m
a
x
+
3
L
+
2
τ
¯
2
exp
⁡
-
X
m
a
x
-
L
2
τ
¯
+
τ
¯
2
+
π
2
τ
¯
2
+
X
m
a
x
+
L
+
2
τ
¯
2
2
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
.Taking into account asymptotic formulas (59), (61), and (65), some algebra leads (73) to the following approximate equation: (74)
4
L
3
η
≈
X
m
a
x
+
L
2
+
4
X
m
a
x
τ
¯
2
H
2
τ
¯
3
e
-
τ
¯
2
≈
4
X
m
a
x
H
2
τ
¯
e
-
τ
¯
2
,
τ
¯
⟶
∞
, which can be solved as(75)
τ
¯
≈
1
2
W
18
X
m
a
x
H
4
L
2
η
2
.


## 5.1. Constant Heat Flux Profile

For a constant heat flux profile, according to (26), (56) is written as(57)
η
=
π
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
.Let us solve (57) approximately. On the one hand, we find in the literature [20, Eqn.7.1.23] the following asymptotic expansion: (58)
π
z
e
z
2
erfc
⁡
z
≈
1
+
∑
m
=
1
∞
-
1
m
2
m
-
1
!
!
2
z
2
m
,
z
⟶
∞
; thus, expanding up to m
=
1 and taking z
=
H
τ
¯, we have(59)
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
≈
1
2
H
2
τ
¯
2
,
τ
¯
⟶
∞
.On the other hand, expanding in Taylor series the following function up to the first order, we have(60)
erf
⁡
τ
¯
+
X
m
a
x
±
L
2
τ
¯
≈
erf
⁡
τ
¯
+
X
m
a
x
±
L
π
τ
¯
e
-
τ
¯
2
,
τ
¯
⟶
∞
; thus,(61)
π
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
≈
2
L
τ
¯
e
-
τ
¯
2
,
τ
¯
⟶
∞
.Substituting (59) and (61) in (57), we get the following approximated equation:(62)
τ
¯
3
e
τ
¯
2
≈
L
H
2
η
.By using the LambertW function [23], (62) can be solved explicitly, obtaining the following approximated expression for the relaxation time:(63)
τ
¯
≈
3
2
W
2
3
L
H
2
η
2
/
3
.Notice that (63) does not depend on X
m
a
x, 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(64)
η
L
=
π
X
m
a
x
+
L
+
2
τ
¯
2
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
-
2
τ
¯
exp
⁡
-
X
m
a
x
+
L
2
τ
¯
+
τ
¯
2
-
exp
⁡
-
X
m
a
x
-
L
2
τ
¯
+
τ
¯
2
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
.In order to solve (64) approximately, let us expand in Taylor series the following function up to first order:(65)
exp
⁡
-
τ
¯
+
X
m
a
x
±
L
2
τ
¯
2
≈
e
-
τ
¯
2
-
X
m
a
x
±
L
e
-
τ
¯
2
,
τ
¯
⟶
∞
; thus,(66)
exp
⁡
-
X
m
a
x
+
L
2
τ
¯
+
τ
¯
2
-
exp
⁡
-
X
m
a
x
-
L
2
τ
¯
+
τ
¯
2
≈
2
L
e
-
τ
¯
2
,
τ
¯
⟶
∞
.Now, substituting (59), (61), and (66) in (64), we arrive at (67)
τ
¯
3
e
τ
¯
2
≈
X
m
a
x
+
L
H
2
η that can be solved explicitly as(68)
τ
¯
≈
3
2
W
2
3
X
m
a
x
+
L
H
2
η
2
/
3
.


## 5.3. Triangular Heat Flux Profile

Let us consider now a triangular heat flux profile, so, according to (45), (56) reads as(69)
η
L
4
=
φ
X
m
a
x
+
L
/
2
τ
¯
+
τ
¯
λ
+
1
+
2
φ
X
m
a
x
-
λ
L
/
2
τ
¯
+
τ
¯
λ
2
-
1
-
φ
X
m
a
x
-
L
/
2
τ
¯
+
τ
¯
λ
-
1
τ
¯
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
.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):(70)
φ
χ
2
τ
¯
+
τ
¯
=
π
χ
2
τ
¯
+
τ
¯
erf
⁡
χ
2
τ
¯
+
τ
¯
+
exp
⁡
-
χ
2
τ
¯
+
τ
¯
2
≈
π
χ
2
τ
¯
+
τ
¯
erf
⁡
τ
¯
+
1
+
χ
2
2
τ
¯
2
e
-
τ
¯
2
,
τ
¯
⟶
∞
.Taking into account (59) and (70) in (69), after some algebra, we arrive at the following approximated equation: (71)
τ
¯
3
e
τ
¯
2
≈
2
L
H
2
η that can be solved explicitly as(72)
τ
¯
≈
3
2
W
2
3
2
L
H
2
η
2
/
3
.As in the constant case, (72) does not depend on X
m
a
x; 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(73)
4
L
2
3
η
=
2
τ
¯
X
m
a
x
+
L
+
2
τ
¯
2
exp
⁡
-
X
m
a
x
+
L
2
τ
¯
+
τ
¯
2
-
X
m
a
x
+
3
L
+
2
τ
¯
2
exp
⁡
-
X
m
a
x
-
L
2
τ
¯
+
τ
¯
2
+
π
2
τ
¯
2
+
X
m
a
x
+
L
+
2
τ
¯
2
2
erf
⁡
X
m
a
x
+
L
2
τ
¯
+
τ
¯
-
erf
⁡
X
m
a
x
-
L
2
τ
¯
+
τ
¯
1
-
π
H
τ
¯
e
H
2
τ
¯
2
erfc
⁡
H
τ
¯
.Taking into account asymptotic formulas (59), (61), and (65), some algebra leads (73) to the following approximate equation: (74)
4
L
3
η
≈
X
m
a
x
+
L
2
+
4
X
m
a
x
τ
¯
2
H
2
τ
¯
3
e
-
τ
¯
2
≈
4
X
m
a
x
H
2
τ
¯
e
-
τ
¯
2
,
τ
¯
⟶
∞
, which can be solved as(75)
τ
¯
≈
1
2
W
18
X
m
a
x
H
4
L
2
η
2
.


## 6. Numerical Results

Table1 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, A
l
2
O
3 (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
k
0
13
60.5
46
k
4.23
×
1
0
-
6
1.77
×
1
0
-
5
1.51
×
1
0
-
5
Grinding  regime
h
27.29
×
1
0
4
4.1
×
1
0
5
1.3
×
1
0
5
2
l
2.66
×
1
0
-
3
1.4
×
1
0
-
3
2.5
×
1
0
-
3
q
5.89
×
1
0
7
1.4
×
1
0
7
1.8
×
1
0
7
v
f
0.53
3.3
×
1
0
-
2
3.3
×
1
0
-
2
T
0
300
300
300
Dimensionless parameters
H
0.335
7.29
2.59
L
83.41
0.65
1.36
T
23.02
79.23
114.3
η
1
0
-
3
1
0
-
3
1
0
-
3Figure3 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).Table2 shows the value of maximum temperature T
m
a
x and its location x
m
a
x
/
l 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 -
L
,
L. 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’sNMaximize command to the surface temperature in the stationary regime, that is, l
i
m
τ
→
∞
T
τ
,
X
,
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.Table 2
Location and value of maximum temperature.
Data 1
Data 2
Data 3
Constant  profile
x
m
a
x
/
l
-
0.973
-
0.301
-
0.501
T
m
a
x  [K]
477.6
329.1
407.6
Linear  profile
x
m
a
x
/
l
0.609
0.687
0.605
T
m
a
x  [K]
537.3
341.3
442.4
Triangular  profile
x
m
a
x
/
l
0.459
0.530
0.459
T
m
a
x  [K]
551.6
343.6
450.4
Parabolic  profile
x
m
a
x
/
l
0.748
0.799
0.751
T
m
a
x  [K]
612.8
353.4
481.4Figures4–6 show the surface temperature (Z
=
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 (x
/
l
=
±
1). This is because the heat source velocity v
f 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).Table3 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.
t
-  [s]
Data 1
Data 2
Data 3
Constant  profile
Exact
5.90
×
1
0
-
3
9.89
×
1
0
-
2
2.00
×
1
0
-
1
Approx.
6.06
×
1
0
-
4
1.11
×
1
0
-
1
1.92
×
1
0
-
1
Linear  profile
Exact
1.66
×
1
0
-
3
1.04
×
1
0
-
1
1.82
×
1
0
-
1
Approx.
1.08
×
1
0
-
4
5.60
×
1
0
-
3
1.77
×
1
0
-
2
Triangular  profile
Exact
1.90
×
1
0
-
3
1.05
×
1
0
-
1
1.83
×
1
0
-
1
Approx.
6.42
×
1
0
-
4
1.36
×
1
0
-
1
2.19
×
1
0
-
1
Parabolic  profile
Exact
1.25
×
1
0
-
3
1.05
×
1
0
-
1
1.83
×
1
0
-
1
Approx.
5.31
×
1
0
-
4
2.10
×
1
0
-
1
2.79
×
1
0
-
1Figures7–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, t
¯
>
2
l
/
v
f. 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.



---

*Source: 5387612-2016-04-24.xml*