- Email: [email protected]
A nonlinear inverse problem in estimating the heat flux of the disc in a disc brake system
Applied Thermal Engineering 31 (2011) 2439e2448
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
A nonlinear inverse problem in estimating the heat flux of the disc in a disc brake system Yu-Ching Yang, Wen-Lih Chen* Clean Energy Center, Department of Mechanical Engineering, Kun Shan University, No. 949, Dawan Rd., Yongkang Dist., Tainan City 710, Taiwan ROC
a r t i c l e i n f o
a b s t r a c t
Article history: Received 26 January 2011 Accepted 7 April 2011 Available online 16 April 2011
In this study, an inverse algorithm based on the conjugate gradient method and the discrepancy principle is applied to estimate the unknown space- and time-dependent heat flux of the disc in a disc brake system from the knowledge of temperature measurements taken within the disc. In the direct problem, the specific heat and thermal conductivity are functions of temperature; hence this is a nonlinear inverse problem. The temperature data obtained from the direct problem are used to simulate the temperature measurements, and the effect of the errors in these measurements upon the precision of the estimated results is also considered. Results show that an excellent estimation on the space- and time-dependent heat flux can be obtained for the test cases considered in this study. The current methodology can be also applied to the prediction of heat flux in the pads, and then the heat partition coefficient between the disc and the pad in a brake system can be calculated. Ó 2011 Elsevier Ltd. All rights reserved.
Keywords: Inverse problem Disc brake Heat flux Conjugate gradient method
1. Introduction In recent decades, disc brakes have increasingly become more popular and gradually found their ways in many different types of vehicles, ranging from light motorcycles to heavy road trucks or even trains. They offer several advantages over drum brakes, including better stopping performance (disc cooled readily), easier to control (not self-applying), less prone to brake fade, and quicker recovery from immersion, which largely contributed to their popularity. However, after long repetitive braking, brake fade could still set in and compromise the performance of a disc brake due to the change of friction characteristics caused by temperature rise and overheating of brake components. Overheated components could further lead to more problems such as thermal cracks, plastic deformation, brake fluid vaporization, and premature wear, and the life span of a disc brake could be shortened as a result. Therefore, to accurately predict the temperature rise in a disc brake system is of eminent importance in the early design stage. Yet, owing to the interactions of several physical phenomena, the prediction of disc brake performance is highly complicated. First, disc braking is a very dynamic and highly transient phenomenon in which extremely large temperature gradients can be generated at the vicinity of disc/ pad interface in a fraction of a second during a high-g deceleration and result in excessive thermal stresses, a phenomenon termed as
* Corresponding author. Tel.: þ886 62050496; fax: þ886 62050509. E-mail address: [email protected] (W.-L. Chen). 1359-4311/$ e see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2011.04.008
“thermal shock”. Thermal shock can give rise to surface cracks and plastic deformation on the disc, making disc surfaces uneven. As a consequence, the heat dissipation area at the discepad interface is largely reduced and becomes non-uniform [1]. Frictional heating also produces thermal deformation and elastic contact which, in turn, largely affect the contact pressure and temperature on sliding interface. These phenomena are termed “thermomechanical coupling or thermomechanical effects”, and are one of the major causes of creating hot spots on the sliding interface. Hot spots could further induce frictional vibration known as “judder”, which affects driving comfort and compromises driving safety if inexperienced drivers react to it improperly [2]. Second, asperities and wear particles can accumulate and form a thin layer called “third body” at the discepad interface. The presence of the third body not only changes friction properties at the discepad interface but also acts as a thermal resistance to the heat dissipated into the disc. Third, as the application of disc brakes has been widened, the geometries of brake components are becoming more sophisticated. For example, the use of drilled-through holes and slots on disc is common in motorcycle’s disc brake systems, and pads made of individual elements have been widely employed on the disc brake systems of high-speed trains. This development toward complex geometry makes thermal analysis even more difficult. Since there exists complex coupling of tribological, mechanical, and thermal behaviors in a disc brake, a comprehensive analysis of its performance is extremely difficult, if not impossible. Then the problem is generally simplified by decoupling the thermal behavior from others through the introduction of some boundary conditions
2440
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
Nomenclature c h J J0 k p pmax q q0 q0 r r1 r2 r3 S T T0 TN t tf z
specific heat (J kg1 K1) convection heat transfer coefficient (W m2 K1) functional gradient of functional thermal conductivity (W m1 K1) direction of descent maximum pressure in the pad (N m2) heat flux density (Wm2) initial guess of heat flux density (W m2) heat flux density at time t ¼ 0 (W m2) radial coordinate (m) inner disc radius (m) inner pad radius (m) outer radius of the disc and the pad (m) frictional surface (m2) temperature (K) initial temperature (K) surrounding temperature (K) time coordinate (s) final time (s) axial coordinate (m)
at the discepad interface. In thermal analysis, accurate estimation of the total energy absorbed by the brake system and the distributions of this energy among the disc and the pads is very crucial. The former can be reasonably assumed to be the same as the reduction of kinetic energy of the vehicle during the braking process, while the determination of the latter is not straightforward. Two interface postulations, perfect contact and imperfect contact, are commonly adopted to resolve this issue. In the perfect contact postulation, the surface temperatures between the disc and the pad at the interface are equal, and with a very short braking time, both disc and pad can further be regarded as semi-infinite solids. Then the ratio between the heat fluxes dissipated into the disc and pad can be obtained analytically [3,4]. However, as the time scale of a brake problem is very short, it is highly transient in nature. Therefore thermal equilibrium at the interface, a state required for equal surface temperature to exist, cannot be maintained at each moment. In the imperfect postulation, on the other hand, the third body is present between the disc and the pad at the discepad interface. It acts as a thermal contact resistance, and the friction heat is assumed generated at the pad surface and then dissipated into the disc through the third body. As a consequence, the surface temperatures of the disc and pad are not equal, which complies with some experimental observations of the clear presence of temperature gradients between pad and disc surfaces. Yet in this practice, difficulty arises to accurately determine the value of the thermal resistance of the third body. It requires some experiments specifically designed to measure thermal resistance under various braking conditions [5]. Furthermore, there have been arguments on the deficiencies of this approach, mainly the overlook of thermal interface inertia and the assumption of a linear temperature gradient across the third body [6]. From the brief review above, it is clear that a more precise friction model for the disc/pad interface is yet to be formulated; and due to the complexity of the problem, developing such a model will be heavily relied on the calibrations with experimental data, especially to establish the relationship among the heat dissipated into the disc and the mechanical, tribological, and thermal parameters involved in a disc brake process.
z2
half disc thickness (m)
Greek symbols D small variation quality a thermal diffusivity (m2 s1) b step size g conjugate coefficient x thermal effusivity (J m2 K1 s0.5) h very small value q0 pad contact angle (rad) l variable used in adjoint problem m friction coefficient r density (kg m3) s standard deviation s transformed time coordinate f heat partition coefficient 6 random variable u0 initial angular velocity of the disc (s1) Superscripts/subscripts K iterative number d disc p pad
The studies of inverse heat conduction problem (IHCP) in recent years have offered convenient alternatives to experimental measurements to obtain accurate thermo-physical quantities such as heat sources, material’s thermal properties, and boundary temperature or heat flux distributions, in many heat conduction problems [7e10]. In these cases, the direct heat conduction problems are concerned with the determination of temperature at interior points of a region when the initial and boundary conditions, heat generation, and material properties are specified, whereas, the IHCP involves the determination of the surface conditions, energy generation, thermo-physical properties, etc., from the knowledge of temperature measurements taken within the body. To date, a variety of analytical and numerical techniques have been developed for the solution of the IHCP, for example, the conjugate gradient method (CGM) [10e13], the Tikhonov regularization method [14], and the genetic algorithm [15], etc. The estimation of heat flux has been the main theme of a number of studies. For example, Niliot and Lefevre [16] adopted a parameter estimation approach based on the boundary element method to solve the inverse problem of point heat source identification. Huang and Lo [17] estimated the internal heat flux of housing for high-speed motors by the CGM. Lee et al. [18] applied the conjugate gradient method to estimate the unknown space- and timedependent heat source in aluminum-coated tapered optical fibers for scanning near-field optical microscopy, by reading the transient temperature data at the measurement positions. Jin and Marin [19] presented the use of the method of fundamental solutions (MFS) for recovering the heat source in steady-state heat conduction problems from boundary temperature and heat flux measurements. Recently, Chen and Yang [20] estimated the heat flux at the roller/workpiece interface during a rolling process by the use of CGM. Recent developments on thermal couple technology and infrared techniques such as pyrometry and thermography have dramatically increased the spatial resolution to the order of millimeter and reduced the response time down to a few micro seconds [21]. Litos et al. [22] developed a measuring system for experimental research on thermomechanical coupling. This system
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
includes a fast temperature measurement subsystem consisting of infrared photon detectors, optical fibers, an inductive sensor, and a data logger. The system allows observation on temperature field revolution on disc, positions of hot spots, evolution of pressure and amount of vibrations. These developments signify that the temperature evolution at the discepad interface can be measured in great details not only spatially but also temporally. Combining such temperature measurement technique with an inverse method would allow the heat flux dissipated into the disc to be estimated very accurately, giving researchers an edge to develop more precise formulations to describe the relationship between heat flux and operational parameters such as sliding velocity, friction coefficient, and discepad interface pressure distribution; hence the performance of a disc brake can be predicted more accurately. Therefore, the focus of the present study is to develop a two-dimensional inverse analysis for estimating the heat dissipated into the disc during a braking process using temperature measurements taken within or on the surface of the disc. The heat flux is crucial to study and to predict the performance of a disc brake system. Due to the wide range of temperature variation in a disc brake process, the assumption of constant thermal properties is not adequate and could lead to significant errors in analytical or numerical solutions. Therefore, to increase the level of accuracy of the mathematical models for a disc brake system, the adoption of temperaturedependent thermal properties is necessary. In this study, we present the conjugate gradient method and the discrepancy principle [23] to estimate the heat flux into the disc by using the simulated temperature measurements. Here, temperaturedependent conductivity and specific heat are adopted, making the proposed inverse method more applicable to a disc brake process. Subsequently, the distributions of temperature in the disc can be determined as well. The conjugate gradient method with an adjoint equation, also called Alifanov’s iterative regularization method, belongs to a class of iterative regularization techniques, which mean the regularization procedure is performed during the iterative processes, thus the determination of optimal regularization conditions is not needed. The CGM is derived based on the perturbation principles and transforms the inverse problem into the solutions of three problems, namely, the direct, the sensitivity, and the adjoint problems, which will be discussed in detail in the following sections. On the other hand, the discrepancy principle is used to terminate the iteration process in the conjugate gradient method. 2. Analysis 2.1. Direct problem Fig. 1(a) shows a car disc brake system which is separated from wheel assembly to show the disc and the pads, while a schematic of the disc and a pad in sliding contact is given in Fig. 1(b). During the course of braking, the frictional heat is generated at the disc/pad interfaces and then dissipated into the disc and the pads. If the heat dissipated into the disc can be estimated correctly, the heat dissipated into the pads also can be calculated accurately because the heat generated at the disc/pad interface is the sum of the heat dissipated into the disc and pads. Therefore in this study, we further separate the pads from the disc and only consider half of the disc as the computational domain due to that the discepad configuration is symmetrical to the center line. The heat dissipated into the disc, either directly from the discepad interface (perfect contact) or through the third body (imperfect contact), can be treated as boundary heat flux into the disc. Since the disc rotates very fast during a braking process, the temperature gradient in q-direction can be neglected. Then the heat flux can be assumed to be applied
2441
Fig. 1. A disc brake assembly; (a) real assembly (picture taken from ref. [4]), (b) model configuration.
on the entire circumferential domain, and the problem is further reduced to an axial-symmetric problem where temperature is a function of r, z, and t only. This approach has been widely adopted in analytical studies [24] and has been proved by experimental investigations [25]. To illustrate the methodology for developing expressions for the use in determining the unknown space- and time-dependent frictional heat flux at the discepad interface in a disc brake system, we consider the following heat transfer problem of the disc as shown in Fig. 2. Under the condition of axisymmetry, the governing equations and the associated boundary and initial conditions for the disc can be written as [4,24]:
1 v vTðr; z; tÞ v vTðr; z; tÞ kðTÞr þ kðTÞ r vr vr vz vz vTðr; z; tÞ ; in r1 < r < r3 ; 0 < z < z2 ; t > 0; vt
(1)
vTðr;z;tÞ ¼ h½TN Tðr;z;tÞ; at r ¼ r1 ;0 < z < z2 ;t > 0; vr
(2)
¼ rcðTÞ kðTÞ
2442
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
The direct problem considered here is concerned with the determination of the medium temperature when the frictional heat flux q(r,t), thermo-physical properties of the disc, and initial and boundary conditions are known. 2.2. Inverse problem For the inverse problem, the function q(r,t) is regarded as being unknown, while everything else in Eqs. (1)e(7) is known. In addition, temperature readings taken within the disc are considered available. The objective of the inverse analysis is to predict the unknown space- and time-dependent function of frictional heat flux, q(r,t), merely from the knowledge of these temperature readings. Let the measured temperature at the measurement position ðr; zÞ ¼ ðri ; zm Þ and time t be denoted by Yi ðri ; zm ; tÞ, i ¼ 1eM, where M represents the number of thermocouples and zm is the z coordinate of all the thermocouples. Then this inverse problem can be stated as follows: by utilizing the above mentioned measured temperature data Yi ðri ; zm ; tÞ, the unknown q(r,t) is to be estimated over the specified space and temporal domain. The solution of the present inverse problem is to be obtained in such a way that the following functional is minimized:
J½qðr; tÞ ¼
Ztf X M t¼0
½Ti ðri ; zm ; tÞ Yi ðri ; zm ; tÞ2 dt;
(8)
i¼1
where Ti ðri ; zm ; tÞ is the estimated (or computed) temperature at the measurement location ðr; zÞ ¼ ðri ; zm Þ. In this study, Ti ðri ; zm ; tÞ are determined from the solution of the direct problem given previously by using an estimated qK ðr; tÞ for the exact qðr; tÞ, here qK ðr; tÞ denotes the estimated quantities at the Kth iteration. tf is the final time of measurement. In addition, in order to develop expressions for the determination of the unknown q(r,t), a “sensitivity problem” and an “adjoint problem” are constructed as described below. 2.3. Sensitivity problem and search step size Fig. 2. Schematic of the disc.
kðTÞ
vTðr;z;tÞ ¼ h½TN Tðr;z;tÞ; at r ¼ r3 ;0 < z < z2 ;t > 0; vr
(3)
kðTÞ
vTðr; z; tÞ ¼ 0 at z ¼ 0; r1 < r < r3 ; t > 0; vz
(4)
kðTÞ
vTðr;z;tÞ ¼ h½TN Tðr;z;tÞ; at z ¼ z2 ;r1 < r < r2 ;t > 0; vz
(5)
vTðr; z; tÞ ¼ qðr; tÞ at z ¼ z2 ; r2 < r < r3 ; t > 0; kðTÞ vz
(6)
Tðr; z; 0Þ ¼ T0 ; for r1 < r < r3 ; 0 < z < z2 ; t ¼ 0;
(7)
where TN is the surrounding temperature and T0 is the initial temperature. r is the density, k(T) and c(T)are the temperaturedependent thermal conductivity and specific heat, respectively, while h is the convection heat transfer coefficient at the outer surfaces. q(r,t) in Eq. (6) is the frictional heat flux into the disc. If the contact pressure between the pad and disc is assumed constant, the frictional heat flux is a function of time and space variable r. This is because work done by friction force grows as r increases. However, if uniform wear contact is adopted, the frictional heat flux is just a function of time and is independent of the space variable r.
The sensitivity problem is obtained from the original direct problem defined by Eqs. (1)e(7) in the following manner: It is assumed that when q(r,t) undergoes a variation q(r,t), T(r,z,t), k(T), and c(T) are perturbed by DT, Dk, and Dc, respectively. Then replacing in the direct problem q(r,t) by qðr; tÞ þ Dqðr; tÞ, T by T þ DT, k(T) by k(T) þ Dk, and c(T) by c(T) þ Dc, subtracting from the resulting expressions the direct problem, using the fact that Dk ¼ ½dkðTÞ=dTDT and Dc ¼ ½dcðTÞ=dTDT, and neglecting the second-order terms, the following sensitivity problem for the sensitivity function DT can be obtained.
1 v v v2 r ½kðTÞDTðr; z; tÞ þ 2 ½kðTÞDTðr; z; tÞ r vr vr vz v ¼ ½rcðTÞDTðr; z; tÞ; in r1 < r < r3 ; 0 < z < z2 ; t > 0; vt
(9)
v ½kðTÞDTðr;z;tÞ ¼ hDTðr;z;tÞ; at r ¼ r1 ;0 < z < z2 ;t > 0; vr
(10)
v ½kðTÞDTðr;z;tÞ ¼ hDTðr;z;tÞ; at r ¼ r3 ;0 < z < z2 ;t > 0; vr
(11)
v ½kðTÞDTðr; z; tÞ ¼ 0; at z ¼ 0; r1 < r < r3 ; t > 0; vz
(12)
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
2443
v ½kðTÞDTðr;z;tÞ ¼ hDTðr;z;tÞ; at z ¼ z2 ;r1 < r < r2 ;t > 0; vz
(13)
kðTÞ
vlðr; z; tÞ ¼ hlðr; z; tÞ; at r ¼ r1 ; 0 < z < z2 ; t > 0; vr
(19)
v ½kðTÞDTðr;z;tÞ ¼ Dqðr;tÞ; at z ¼ z2 ;r2 < r < r3 ;t > 0; vz
(14)
kðTÞ
vlðr; z; tÞ ¼ hlðr; z; tÞ at r ¼ r3 ; 0 < z < z2 ; t > 0; vr
(20)
DTðr; z; 0Þ ¼ 0; for r1 < r < r3 ; 0 < z < z2 ; t ¼ 0:
(15)
kðTÞ
vlðr; z; tÞ ¼ 0 at z ¼ 0; r1 < r < r3 ; t > 0; vz
(21)
The sensitivity problem of Eqs. (9)e(15) can be solved by the same method as the direct problem of Eqs. (1)e(7).
kðTÞ
vlðr; z; tÞ ¼ hlðr; z; tÞ at z ¼ z2 ; r1 < r < r2 ; t > 0; vz
(22)
2.4. Adjoint problem and gradient equation
kðTÞ
vlðr; z; tÞ ¼ 0 at z ¼ z2 ; r2 < r < r3 ; t > 0; vz
(23)
To formulate the adjoint problem, Eq. (1) is multiplied by the Lagrange multiplier (or adjoint function) l, and the resulting expressions are integrated over the time and correspondent space domains. Then the results are added to the right hand side of Eq. (8) to yield the following expression for the functional J½qðr; tÞ:
Zz2
Ztf
Zr3 X M
J½qðr; tÞ ¼
½Ti ðr; z; tÞ
i¼1
t¼0 z¼0 r¼r1
lðr; z; tÞ ¼ 0 for r1 < r < r3 ; 0 < z < z2 ; t ¼ tf :
The adjoint problem is different from the standard initial value problem in that the final time condition at time t ¼ tf is specified instead of the customary initial condition at time t ¼ 0. However, this problem can be transformed to an initial value problem by the transformation of the time variable as s ¼ tf t. Then the adjoint problem can be solved by the same method as the direct problem. Finally the following integral term is left:
Yi ðr; z; tÞ2 dðr ri Þdðz zm Þdrdzdt Ztf
Zz2
DJ ¼
1 v vT r$lðr; z; tÞ$ kðTÞr r vr vr
Zr3
þ t¼0 z¼0 r¼r1
Zz2
Zr3 X M
t¼0 z¼0 r¼r1
(16)
DJ ¼
i¼1
Zr3
Ztf
rJ 0 ðr; tÞDqðr; tÞdrdt;
(26)
t¼0 r¼r2
where J0 (t) is the gradient of the functional J(q). A comparison of Eqs. (25) and (26) leads to the following form:
J 0 ðr; tÞ ¼ lðr; z2 ; tÞ:
(27)
The following iteration process based on the conjugate gradient method is now used for the estimation of q(r,t) by minimizing the above functional J[q(r,t)]:
Yi ðr; z; tÞDT dðr ri Þdðz zm Þdrdzdt ( Ztf Zz2 Zr3 v v l$ þ $r$ ½kðTÞDT vr vr
K qKþ1 ðr; tÞ ¼ qK ðr; tÞ b pK r; t ;
)
v2 v þ r$ 2 ½kðTÞDT r$ ½rcðTÞDT drdzdt vt vz
(17)
We can integrate the second triple-integral term in Eq. (17) by parts, utilizing the initial and boundary conditions of the sensitivity problem. Then DJ is allowed to go to zero. The vanishing of the integrands containing DT leads to the following adjoint problem for the determination of lðr; z; tÞ:
pK ðr; tÞ ¼ J 0 Kðr; tÞ þ gK pK1 r; t ;
gK ¼ (18)
(28)
(29)
which is conjugation of the gradient direction J 0 Kðr; tÞ at iteration K and the direction of descent pK1 ðr; tÞ at iteration K 1. The conjugate coefficient gK is determined from:
Ztf X M h
M 1 X þ $ 2½Ti ðr; z; tÞ Yi ðr; z; tÞdðr ri Þdðz zm Þ r i¼1
K ¼ 0; 1; 2; .;
where bK is the search step size in going from iteration K to iteration Kþ1, and pK(r,t) is the direction of descent (i.e., search direction) given by:
1 v vlðr; z; tÞ v2 lðr; z; tÞ vlðr; z; tÞ þ rcðTÞ $kðTÞ$ r þ kðTÞ$ r vr vr vt vz2
¼ 0; in r1 < r < r3 ; 0 < z < z2 ; t > 0;
(25)
2.5. Conjugate gradient method for minimization
2½Ti ðr; z; tÞ
t¼0 z¼0 r¼r1
r lðr; z2 ; tÞDqðr; tÞdrdt:
t¼0 r¼r2
where dð,Þ is the Dirac function. The variation DJ is derived by perturbing q(r,t) by qðr; tÞ þ Dqðr; tÞ, T by T þ DT, k(T) by k(T) þ Dk, and c(T) by c(T) þ Dc, respectively, in Eq. (16). Subtracting from the resulting expressions the original Eq. (16), using the fact that Dk ¼ ½dkðTÞ=dTDT and Dc ¼ ½dcðTÞ=dTDT, and neglecting the secondorder terms, we thus find:
Ztf
Zr3
Ztf
From the definition used in Ref. [7], we have:
v vT vT kðTÞ rcðTÞ drdzdt þ vz vz vt
DJ½qðr; tÞ ¼
(24)
t¼0
i2 J 0 Kðri ; tÞ dt
i¼1
Ztf X M h i2 J 0 K 1ðri ; tÞ dt t¼0
i¼1
with g0 ¼ 0:
(30)
2444
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
The convergence of the above iterative procedure in minimizing the functional J is proved in Ref. [7]. To perform the iterations according to Eq. (28), we need to compute the step size bK and the gradient of the functional J 0 Kðr; tÞ. The functional J½qKþ1 ðr; tÞ for iteration K þ 1 is obtained by rewriting Eq. (8) as:
Ztf X M h h i i2 K Ti qK b pK Yi ðri ;zm ;tÞ dt; J qKþ1 ðr;tÞ ¼
Inner disc diameter Outer disc diameter Mean sliding radius Sliding length Disc thickness Pad thickness Initial speed Braking time Deceleration Total energy Convection coefficient
(31)
i¼1
t¼0
where we replace qKþ1 by the expression given by Eq. (28). If K temperature Ti ðqK b pK Þ is linearized by a Taylor expansion, Eq. (31) takes the form:
h i J qKþ1 ðr; tÞ ¼
Table 1 Brake dimensions and operating conditions.
Ztf t¼0
M h X
K Ti ðqK Þ b DTi ðpK Þ
i¼1
i2 Yi ðri ; zm ; tÞ dtdt;
ð32Þ
properties are taken as constant. Furthermore, under the assumption of constant deceleration and uniform wear contact between the pad and the disc, the boundary heat flux in Eq. (6) is a function of time only and can be written as:
Ztf X M
bK ¼
t¼0
DTi pK
h i Ti qK Yi dt
i¼1
Z tf X M h
i2 DTi ðp Þ dt
t¼0
:
(33)
2.6. Stopping criterion If the problem contains no measurement errors, the traditional check condition specified as
(34)
where h is a small specified number, can be used as the stopping criterion. However, the observed temperature data contains measurement errors; as a result, the inverse solution will tend to approach the perturbed input data, and the solution will exhibit oscillatory behavior as the number of iteration is increased [26]. Computational experience has shown that it is advisable to use the discrepancy principle [23] for terminating the iteration process in the conjugate gradient method. Assuming Ti ðri ; zm ; tÞ Yi ðri ; zm ; tÞys, the stopping criteria h by the discrepancy principle can be obtained from Eq. (8) as
h ¼ Ms2 tf ;
(36)
mfpmax r2 u0 ;
(37)
with
q0 ¼
q0 2p
where pmax is the maximum pressure distributed in the pad and f is the heat partition coefficient between the disc and the pad defined as:
f¼
x d Sd ; xd Sd þ xp Sp
(38)
K
i¼1
J qKþ1 < h;
! t 1 ; tf
qðtÞ ¼ q0
ðqK Þ
where Ti is the solution of the direct problem at ðr; zÞ ¼ ðri ; zm Þ by using estimated qK(r,t) for exact q(r,t) at time t. The sensitivity function DTi ðpK Þ is taken as the solution of Eqs. (9)e(15) at the measured position ðr; zÞ ¼ ðri ; zm Þ by letting Dq ¼ pK [7]. The search step size bK is determined by minimizing the functional given by Eq. (32) with respect to bK. The following expression can be obtained:
0.132 (m) 0.227 (m) 0.094.5 (m) 0.037 (m) 0.011 (m) 0.010 (m) 100 (km h1) 3.96 (s) 7 (ms2) 165 (kJ) 60 (W m2 K1)
(35)
where s is the standard deviation of the measurement error. Then the stopping criterion is given by Eq. (34) with h determined from Eq. (35). The computational procedure is the same as in [20], hence there is no need to repeat here. 3. Results and discussion A bench-mark test case in [24] is first considered to validate the code used in this study. The details of the dimensions and operating conditions are given in Tables 1 and 2. In this test case, the thermal
where xd and xP are the thermal effusivities of the disc and the pad, and and Sd and Sp are frictional contact surfaces of the disc pffiffiffiffiffiffiffi ffi the pad, respectively. Thermal effusivity x is defined as x ¼ krc. The numerical procedure in this paper is based on the unstructured-mesh, fully collocated, finite-volume code, ‘USTREAM’ developed by W.-L. Chen. This is the descendent of the structuredmesh, multi-block code of ‘STREAM’ [27]. The computational mesh consists of 860 cells, where there are 20 and 43 cells, respectively, in the z- and r-directions. This mesh has been proved in a gridindependent test to achieve grid-independent solutions. The comparison of mean sliding radius surface temperature variations is given in Fig. 3. To find an appropriate time-step interval to achieve time-step-interval independent solutions, the results are computed by using three different numbers of time steps, namely 100, 200, and 300. It can be seen that the results by the three time-step numbers are almost identical and are in very good agreement with Newcomb’s analytical solution. This indicates that the time step number of 100, corresponding to a time-step interval of 0.0396 s, is good enough to ensure time-step-interval independent solutions in this problem. In addition, the test also versifies the correctness of the code used in this study. To test the effect of adopting temperature-dependent thermophysical properties, the disc is assumed made of AISI-1045 steel, which is commonly used in friction elements. Since the density of AISI-1045 steel has relatively weak temperature dependence Table 2 Thermo-physical properties of disc and pad.
Conductivity (W m1 K1) Density (kg m3) Specific heat (J kg1 K1)
Disc
Pad
43.5 7850 445
12 2500 900
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
2445
a
Fig. 3. Comparison disc mean sliding radius surface temperature variations.
over a wide temperature range [28], a constant average value of 7844 kg m3 is adopted. However, the specific heat and thermal diffusivity from the experimental measurements of Touloukian [29] are used here, and Fig. 4 shows the variations of specific heat and thermal conductivity with respect to temperature from 280 K to 1300 K. The mean sliding radius temperature variations computed with temperature-dependent thermo-physical properties are also shown in Fig. 3. Although the results are not significantly different from those adopting constant thermo-physical properties, the deviation, about 5 K at most places, is clearly present. To ensure high level of accuracy, temperature-dependent thermo-physical properties should be used, and in the following, all results were obtained using temperature-dependent thermo-physical properties. The objective of this article is to validate the current approach when used in estimating the unknown heat flux into the disc during a braking process accurately with no prior information on the functional form of the unknown quantities, a procedure called function estimation. In this study, we did not have an experiment setup to measure the temperatures on the disc. Instead, we substitute an assumed interface heat flux such as Eq. (36) into the direct problem of Eqs. (1)e(7) to calculate the temperatures at the locations where the temperature measurements are taken. The results are termed the computed temperature Yexact ðri ; zm ; tÞ. However, in reality, temperature measurements always contain some degree of error, whose magnitude depends upon the particular measuring method employed. To take measurement errors into account, a random error noise is added to the above computed temperature Yexact ðri ; zm ; tÞ to obtain the measured temperature Yi ðri ; zm ; tÞ. Hence, the measured temperature Yi ðri ; zm ; tÞ is expressed as:
Yi ðri ; zm ; tÞ ¼ Yexact ðri ; zm ; tÞ þ 6s;
(39)
where 6 is a random variable within 2.576 to 2.576 for a 99% confidence bounds, and s is the standard deviation of the measurement. The measured temperature Yi ðri ; zm ; tÞ generated in such way is called the simulated measurement temperature. At first, the temperature measurement is assumed located at zm ¼ 0.0054 m, equivalent to 0.1 mm inside the disc, and along the friction interface (see Fig. 2). The comparison between the inverse and exact heat flux, using Eq. (36), at three different time steps is given in Fig. 5. The inverse results were obtained with an initial guess of q0(r,t) ¼ 1.0 105 W m2 and measurement errors s ¼ 0.0.
b
Fig. 4. Measured AISI-1045 steel specific heat and thermal conductivity variations against temperature. (a) Specific heat variations, and (b) thermal conductivity variations.
Fig. 5 shows that for this simple heat flux function, the inverse estimation is almost identical to the exact value. Next, we adopt the uniform pressure assumption in Talati and Jalalifar [4], where the heat flux function is written as:
!
r t 1 ; qðr; tÞ ¼ q0 r2 tf
(40)
Eq. (40) indicates that the heat flux is now a function of radius and time. The results at three different time steps are shown in Fig. 6, and yet again very good agreement between the inverse and the exact values can be seen. Here the inverse results were obtained with an initial guess of q0(r,t) ¼ 1.0 105 W m2 and measurement errors s ¼ 0.0. Once the heat flux has been estimated correctly, the temperature evolution anywhere on the disc can be calculated accurately too. This is shown in Fig. 7, where the comparison of temperature evolution at four different locations on the disc surface is given. It is noticeable that the variations of calculated and exact temperature are almost identical.
2446
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
Fig. 5. Estimated q(t) (Eq. (37)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
Fig. 7. Estimated disc surface temperature variations at 4 different radial positions with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
In order to investigate the influence of measurement error upon the estimated results, Fig. 8 shows the estimated values of the unknown function q(r,t), obtained with the initial guesses q0(r,t) ¼ 1.0 105 W m2, and measurement error of deviation s ¼ 0.005. The deviation s ¼ 0.005 corresponds to 1.28% of measurement error, which is roughly in the range of 4e6 K over the entire disc surface temperature variation during the braking process. This magnitude of temperature measurement error is huge according to the accuracy standard of modern thermal couples, whose precision is normally within 1 K. Despite this unrealistically large measurement error, the results indicate that this temperature measurement error in general creates only about 2% of error on heat flux magnitude, demonstrating that the accuracy of the current method is not very sensitive to temperature measurement error. Since the direct problem is a transient heat transfer problem, there have been studies reporting the sensitivity of inverse accuracy on the relative distance between the estimated and measurement quantities in inverse transient heat conduction problems [30]. Under some
circumstances, the accuracy deteriorates rapidly as the measurement location moves away from the estimated quantity location. Therefore in the current inverse problem, this issue needs to be investigated. To examine the sensitivity of inverse method’s accuracy on measurement location, the thermal couple is assumed to be located respectively at three other positions, zm ¼ 0.005 m, 0.0035 m, and 0.0001 m. The estimated heat flux distributions, together with those of zm ¼ 0.0054 m, are plotted in Fig. 9. These results were all obtained with an initial guess of q0(r,t) ¼ 1.0 105 W m2 and measurement errors s ¼ 0.0, and they obviously show that the inverse accuracy deteriorates progressively as the measurement position moves away from the friction interface, especially at the vicinity of r ¼ r2. In the case of zm ¼ 0.0001 m, the inverse heat flux at t ¼ 3.66 s is in general half the magnitude of the exact value, resulting in an discrepancy even larger than that produced by taking measurement error into account.
Fig. 6. Estimated q(r,t) (Eq. (40)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
Fig. 8. Estimated q(r,t) (Eq. (40)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.005.
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
Fig. 9. Estimated q(r,t) (Eq. (40)) with four different measurement locations, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
2447
Fig. 11. Estimated q(r,t) (Eq. (41)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
a In order to demonstrate the capability of the present method in obtaining an accurate estimation of arbitrary form of the unknown heat flux, we consider another case of q(r,t), based on the thermoelastic simulation in Choi and Lee [2], which is calculated by:
qðr; tÞ ¼ mr uðtÞpðr; tÞ;
b
(41)
where u(t) and p(r,t) are distributions taken from Figs. 4 and 6 in [2]. The heat generation during the braking stage is shown in Fig. 10(a). Please note that due to data digitizing and scaling factors, the absolute magnitude of heat flux shown here is not identical to that in [2] but the tendency of data variation is preserved. In their study, heat conduction and elastic equations are coupled through the heat generated at discepad interface calculated by Eq. (41). Therefore, the current inverse method can be applied to thermoelastic approach to estimate the heat flux into the disc and then evaluate the heat partition between the disc and pad at the friction interface. Fig. 10(b) shows the estimated results of q(r,t), obtained with zm ¼ 0.0054 m, the initial guesses q0(r,t) ¼ 1.0 105 W m2, and measurement error of deviation s ¼ 0.0. Excellent agreement between the exact and inverse heat flux can be seen by comparing Fig. 10(a) and (b). More detailed comparison on the heat flux distributions along the radial direction at three different time steps given in Fig. 11 further confirms the accuracy of the inverse estimation. The results indicate that the current inverse method performs very well even on this highly complicated heat flux function. 4. Conclusion
Fig. 10. The heat flux function taken from Choi and Lee [2]; (a) exact q(r,t) distributions, and (b) inverse q(r,t) distributions.
An inverse algorithm based on the conjugate gradient method and the discrepancy principle was successfully applied for the solution of the inverse problem to determine the unknown space- and time-dependent heat flux into the disc during a braking process, while knowing the temperature history at some measurement locations. Due to the wide range of temperature variation in a braking process, temperature-dependent specific heat and thermal conductivity are used here; hence the governing equation of the direct problem becomes nonlinear. Numerical results confirm that the method proposed herein can accurately estimate the space- and time-dependent heat flux and temperature distributions for the problem even involving the inevitable measurement errors. The
2448
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
proposed inverse algorithm does not require prior information for the functional form of the unknown quantities to perform the inverse calculation, and excellent estimated values can be obtained for the considered problem. However, the results also show that the inverse accuracy is very sensitive to the distance between the estimated and measurement quantities. The temperature measurement location should be put as close to the friction interface as possible because the accuracy of the inverse method deteriorates rapidly as the measurement location moves away from the friction interface. Acknowledgements This work was supported by the National Science Council, Taiwan, Republic of China, under the grant numbers NSC 98-2221E-168-035 and NSC 98-2221-E-168-034. References [1] P. Dufrenoy, Two-/three-dimensional hybrid model of the thermomechanical behavior of disc brakes, Proc. IMechE Part F: J. Rail Rapid Transit. 218 (2004) 17e30. [2] J.H. Choi, I. Lee, Finite element analysis of transient thermoelastic behaviors in disc brakes, Wear 257 (2004) 47e58. [3] C.H. Gao, X.Z. Lin, Transient temperature field analysis of a brake in a nonaxisymmetric three-dimensional model, J. Mater. Process. Technol. 129 (2002) 513e517. [4] F. Talati, S. Jalalifar, Analysis of heat conduction in a disk brake system, Heat Mass Transfer 45 (2009) 1047e1059. [5] N. Laraqi, Phenomene de constriction thermique dans les contacts glissants (thermal constriction phenomenon in sliding contacts), Int. J. Heat Mass Transfer 39 (1996) 3717e3724. [6] D. Majcherczak, P. Dufrenoy, M. Nait-Abdelaziz, Third body influence on thermal friction contact problems: application to braking, ASME J. Tribol. 127 (2005) 89e95. [7] O.M. Alifanov, Inverse Heat Transfer Problem. Springer-Verlag, New York, 1994. [8] C.H. Huang, C.C. Shih, A shape identification problem in estimating simultaneously two interfacial configurations in a multiple region domain, Appl. Therm. Eng. 26 (2006) 77e88. [9] Y.C. Yang, Simultaneously estimating the contact heat and mass transfer coefficients in a double-layer hollow cylinder with interface resistance, Appl. Therm. Eng. 27 (2007) 501e508. [10] H.T. Chen, H.C. Wang, Estimation of heat-transfer characteristics on a rectangular fin under wet conditions, Int. J. Heat Mass Transfer 51 (2008) 2123e2138. [11] W.L. Chen, Y.C. Yang, S.S. Chu, Estimation of heat generation at the interface of cylindrical bars during friction process, Appl. Therm. Eng. 29 (2009) 351e357.
[12] Y.C. Yang, W.L. Chen, An iterative regularization method in simultaneously estimating the inlet temperature and heat-transfer rate in a forced-convection pipe, Int. J. Heat Mass Transfer 52 (2009) 1928e1937. [13] W.L. Chen, Y.C. Yang, Simultaneous estimation of heat-transfer rate and coolant fluid velocity in a transpiration cooling process, Int. J. Therm. Sci. 49 (2010) 1407e1416. [14] C. Le Niliot, F. Lefevre, Multiple transient point heat sources identification in heat diffusion: application to numerical two- and three-dimensional problems, Numer. Heat Transfer Part B 39 (2001) 277e301. [15] M. Raudensky, K.A. Woodbury, J. Kral, Genetic algorithm in solution of inverse heat conduction problems, Numer. Heat Transfer Part B 28 (1995) 293e306. [16] C. Le Niliot, F. Lefevre, A parameter estimation approach to solve the inverse problem of point heat sources identification, Int. J. Heat Mass Transfer 47 (2004) 827e841. [17] C.H. Huang, H.C. Lo, A three-dimensional inverse problem in estimating the internal heat flux of housing for high speed motors, Appl. Therm. Eng. 26 (2006) 1515e1529. [18] H.L. Lee, W.J. Chang, W.L. Chen, Y.C. Yang, An inverse problem of estimating the heat source in tapered optical fibers for scanning near-field optical microscopy, Ultramicroscopy 107 (2007) 656e662. [19] B. Jin, L. Marin, The method of fundamental solutions for inverse source problems associated with the steady-state heat conduction, Int. J. Numer. Methods Eng. 69 (2007) 1570e1589. [20] W.L. Chen, Y.C. Yang, Inverse problem of estimating the heat flux at the roller/ workpiece interface during a rolling process, Appl. Therm. Eng. 30 (2010) 1247e1254. [21] J. Thevenet, M. Siroux, B. Desmet, Measurements of brake disc surface temperature and emissivity by two-color pyrometry, Appl. Therm. Eng. 30 (2010) 753e759. [22] P. Litos, M. Honner, V. Lang, J. Bartik, M. Hynek, A measuring system for experimental research on the thermomechanical coupling of disc brakes, Proc. IMechE, Part D: J. Automobile Eng. 222 (2008) 1247e1257. [23] O.M. Alifanov, E.A. Artyukhin, Regularized numerical solution of nonlinear inverse heat-conduction problem, J. Eng. Phys. 29 (1975) 934e938. [24] T.P. Newcomb, Transient temperatures attained in disk brakes, Br. J. Appl. Phys. 10 (1959) 339e340. [25] D. Majcherczak, P. Dufrenoy, Y. Berthier, Tribological, thermal and mechanical coupling aspects of the dry sliding contact, Tribol. Int. 40 (2007) 834e843. [26] O.M. Alifanov, Application of the regularization principle to the formulation of approximate solution of inverse heat conduction problem, J. Eng. Phys. 23 (1972) 1566e1571. [27] F.S. Lien, W.L. Chen, M.A. Leschziner, A multiblock implementation of a nonorthogonal, collocated finite volume algorithm for complex turbulent flows, Int. J. Numer. Methods Fluids 23 (1996) 567e588. [28] R.D. Pelhlke, A. Jeyarajan, H. Wada, Report No. NSF/MAE-82028, NSF Applied Research Division, University of Michigan (1982). [29] Y.S. Touloukian, Thermophysical Properties of High Temperature Solid Materials. Macmillan, New York, 1967. [30] W.L. Chen, Y.C. Yang, Inverse prediction of frictional heat flux and temperature in sliding contact with a protective strip by iterative regularization method, Appl. Math. Model. 35 (2011) 2874e2886.
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
A nonlinear inverse problem in estimating the heat flux of the disc in a disc brake system Yu-Ching Yang, Wen-Lih Chen* Clean Energy Center, Department of Mechanical Engineering, Kun Shan University, No. 949, Dawan Rd., Yongkang Dist., Tainan City 710, Taiwan ROC
a r t i c l e i n f o
a b s t r a c t
Article history: Received 26 January 2011 Accepted 7 April 2011 Available online 16 April 2011
In this study, an inverse algorithm based on the conjugate gradient method and the discrepancy principle is applied to estimate the unknown space- and time-dependent heat flux of the disc in a disc brake system from the knowledge of temperature measurements taken within the disc. In the direct problem, the specific heat and thermal conductivity are functions of temperature; hence this is a nonlinear inverse problem. The temperature data obtained from the direct problem are used to simulate the temperature measurements, and the effect of the errors in these measurements upon the precision of the estimated results is also considered. Results show that an excellent estimation on the space- and time-dependent heat flux can be obtained for the test cases considered in this study. The current methodology can be also applied to the prediction of heat flux in the pads, and then the heat partition coefficient between the disc and the pad in a brake system can be calculated. Ó 2011 Elsevier Ltd. All rights reserved.
Keywords: Inverse problem Disc brake Heat flux Conjugate gradient method
1. Introduction In recent decades, disc brakes have increasingly become more popular and gradually found their ways in many different types of vehicles, ranging from light motorcycles to heavy road trucks or even trains. They offer several advantages over drum brakes, including better stopping performance (disc cooled readily), easier to control (not self-applying), less prone to brake fade, and quicker recovery from immersion, which largely contributed to their popularity. However, after long repetitive braking, brake fade could still set in and compromise the performance of a disc brake due to the change of friction characteristics caused by temperature rise and overheating of brake components. Overheated components could further lead to more problems such as thermal cracks, plastic deformation, brake fluid vaporization, and premature wear, and the life span of a disc brake could be shortened as a result. Therefore, to accurately predict the temperature rise in a disc brake system is of eminent importance in the early design stage. Yet, owing to the interactions of several physical phenomena, the prediction of disc brake performance is highly complicated. First, disc braking is a very dynamic and highly transient phenomenon in which extremely large temperature gradients can be generated at the vicinity of disc/ pad interface in a fraction of a second during a high-g deceleration and result in excessive thermal stresses, a phenomenon termed as
* Corresponding author. Tel.: þ886 62050496; fax: þ886 62050509. E-mail address: [email protected] (W.-L. Chen). 1359-4311/$ e see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2011.04.008
“thermal shock”. Thermal shock can give rise to surface cracks and plastic deformation on the disc, making disc surfaces uneven. As a consequence, the heat dissipation area at the discepad interface is largely reduced and becomes non-uniform [1]. Frictional heating also produces thermal deformation and elastic contact which, in turn, largely affect the contact pressure and temperature on sliding interface. These phenomena are termed “thermomechanical coupling or thermomechanical effects”, and are one of the major causes of creating hot spots on the sliding interface. Hot spots could further induce frictional vibration known as “judder”, which affects driving comfort and compromises driving safety if inexperienced drivers react to it improperly [2]. Second, asperities and wear particles can accumulate and form a thin layer called “third body” at the discepad interface. The presence of the third body not only changes friction properties at the discepad interface but also acts as a thermal resistance to the heat dissipated into the disc. Third, as the application of disc brakes has been widened, the geometries of brake components are becoming more sophisticated. For example, the use of drilled-through holes and slots on disc is common in motorcycle’s disc brake systems, and pads made of individual elements have been widely employed on the disc brake systems of high-speed trains. This development toward complex geometry makes thermal analysis even more difficult. Since there exists complex coupling of tribological, mechanical, and thermal behaviors in a disc brake, a comprehensive analysis of its performance is extremely difficult, if not impossible. Then the problem is generally simplified by decoupling the thermal behavior from others through the introduction of some boundary conditions
2440
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
Nomenclature c h J J0 k p pmax q q0 q0 r r1 r2 r3 S T T0 TN t tf z
specific heat (J kg1 K1) convection heat transfer coefficient (W m2 K1) functional gradient of functional thermal conductivity (W m1 K1) direction of descent maximum pressure in the pad (N m2) heat flux density (Wm2) initial guess of heat flux density (W m2) heat flux density at time t ¼ 0 (W m2) radial coordinate (m) inner disc radius (m) inner pad radius (m) outer radius of the disc and the pad (m) frictional surface (m2) temperature (K) initial temperature (K) surrounding temperature (K) time coordinate (s) final time (s) axial coordinate (m)
at the discepad interface. In thermal analysis, accurate estimation of the total energy absorbed by the brake system and the distributions of this energy among the disc and the pads is very crucial. The former can be reasonably assumed to be the same as the reduction of kinetic energy of the vehicle during the braking process, while the determination of the latter is not straightforward. Two interface postulations, perfect contact and imperfect contact, are commonly adopted to resolve this issue. In the perfect contact postulation, the surface temperatures between the disc and the pad at the interface are equal, and with a very short braking time, both disc and pad can further be regarded as semi-infinite solids. Then the ratio between the heat fluxes dissipated into the disc and pad can be obtained analytically [3,4]. However, as the time scale of a brake problem is very short, it is highly transient in nature. Therefore thermal equilibrium at the interface, a state required for equal surface temperature to exist, cannot be maintained at each moment. In the imperfect postulation, on the other hand, the third body is present between the disc and the pad at the discepad interface. It acts as a thermal contact resistance, and the friction heat is assumed generated at the pad surface and then dissipated into the disc through the third body. As a consequence, the surface temperatures of the disc and pad are not equal, which complies with some experimental observations of the clear presence of temperature gradients between pad and disc surfaces. Yet in this practice, difficulty arises to accurately determine the value of the thermal resistance of the third body. It requires some experiments specifically designed to measure thermal resistance under various braking conditions [5]. Furthermore, there have been arguments on the deficiencies of this approach, mainly the overlook of thermal interface inertia and the assumption of a linear temperature gradient across the third body [6]. From the brief review above, it is clear that a more precise friction model for the disc/pad interface is yet to be formulated; and due to the complexity of the problem, developing such a model will be heavily relied on the calibrations with experimental data, especially to establish the relationship among the heat dissipated into the disc and the mechanical, tribological, and thermal parameters involved in a disc brake process.
z2
half disc thickness (m)
Greek symbols D small variation quality a thermal diffusivity (m2 s1) b step size g conjugate coefficient x thermal effusivity (J m2 K1 s0.5) h very small value q0 pad contact angle (rad) l variable used in adjoint problem m friction coefficient r density (kg m3) s standard deviation s transformed time coordinate f heat partition coefficient 6 random variable u0 initial angular velocity of the disc (s1) Superscripts/subscripts K iterative number d disc p pad
The studies of inverse heat conduction problem (IHCP) in recent years have offered convenient alternatives to experimental measurements to obtain accurate thermo-physical quantities such as heat sources, material’s thermal properties, and boundary temperature or heat flux distributions, in many heat conduction problems [7e10]. In these cases, the direct heat conduction problems are concerned with the determination of temperature at interior points of a region when the initial and boundary conditions, heat generation, and material properties are specified, whereas, the IHCP involves the determination of the surface conditions, energy generation, thermo-physical properties, etc., from the knowledge of temperature measurements taken within the body. To date, a variety of analytical and numerical techniques have been developed for the solution of the IHCP, for example, the conjugate gradient method (CGM) [10e13], the Tikhonov regularization method [14], and the genetic algorithm [15], etc. The estimation of heat flux has been the main theme of a number of studies. For example, Niliot and Lefevre [16] adopted a parameter estimation approach based on the boundary element method to solve the inverse problem of point heat source identification. Huang and Lo [17] estimated the internal heat flux of housing for high-speed motors by the CGM. Lee et al. [18] applied the conjugate gradient method to estimate the unknown space- and timedependent heat source in aluminum-coated tapered optical fibers for scanning near-field optical microscopy, by reading the transient temperature data at the measurement positions. Jin and Marin [19] presented the use of the method of fundamental solutions (MFS) for recovering the heat source in steady-state heat conduction problems from boundary temperature and heat flux measurements. Recently, Chen and Yang [20] estimated the heat flux at the roller/workpiece interface during a rolling process by the use of CGM. Recent developments on thermal couple technology and infrared techniques such as pyrometry and thermography have dramatically increased the spatial resolution to the order of millimeter and reduced the response time down to a few micro seconds [21]. Litos et al. [22] developed a measuring system for experimental research on thermomechanical coupling. This system
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
includes a fast temperature measurement subsystem consisting of infrared photon detectors, optical fibers, an inductive sensor, and a data logger. The system allows observation on temperature field revolution on disc, positions of hot spots, evolution of pressure and amount of vibrations. These developments signify that the temperature evolution at the discepad interface can be measured in great details not only spatially but also temporally. Combining such temperature measurement technique with an inverse method would allow the heat flux dissipated into the disc to be estimated very accurately, giving researchers an edge to develop more precise formulations to describe the relationship between heat flux and operational parameters such as sliding velocity, friction coefficient, and discepad interface pressure distribution; hence the performance of a disc brake can be predicted more accurately. Therefore, the focus of the present study is to develop a two-dimensional inverse analysis for estimating the heat dissipated into the disc during a braking process using temperature measurements taken within or on the surface of the disc. The heat flux is crucial to study and to predict the performance of a disc brake system. Due to the wide range of temperature variation in a disc brake process, the assumption of constant thermal properties is not adequate and could lead to significant errors in analytical or numerical solutions. Therefore, to increase the level of accuracy of the mathematical models for a disc brake system, the adoption of temperaturedependent thermal properties is necessary. In this study, we present the conjugate gradient method and the discrepancy principle [23] to estimate the heat flux into the disc by using the simulated temperature measurements. Here, temperaturedependent conductivity and specific heat are adopted, making the proposed inverse method more applicable to a disc brake process. Subsequently, the distributions of temperature in the disc can be determined as well. The conjugate gradient method with an adjoint equation, also called Alifanov’s iterative regularization method, belongs to a class of iterative regularization techniques, which mean the regularization procedure is performed during the iterative processes, thus the determination of optimal regularization conditions is not needed. The CGM is derived based on the perturbation principles and transforms the inverse problem into the solutions of three problems, namely, the direct, the sensitivity, and the adjoint problems, which will be discussed in detail in the following sections. On the other hand, the discrepancy principle is used to terminate the iteration process in the conjugate gradient method. 2. Analysis 2.1. Direct problem Fig. 1(a) shows a car disc brake system which is separated from wheel assembly to show the disc and the pads, while a schematic of the disc and a pad in sliding contact is given in Fig. 1(b). During the course of braking, the frictional heat is generated at the disc/pad interfaces and then dissipated into the disc and the pads. If the heat dissipated into the disc can be estimated correctly, the heat dissipated into the pads also can be calculated accurately because the heat generated at the disc/pad interface is the sum of the heat dissipated into the disc and pads. Therefore in this study, we further separate the pads from the disc and only consider half of the disc as the computational domain due to that the discepad configuration is symmetrical to the center line. The heat dissipated into the disc, either directly from the discepad interface (perfect contact) or through the third body (imperfect contact), can be treated as boundary heat flux into the disc. Since the disc rotates very fast during a braking process, the temperature gradient in q-direction can be neglected. Then the heat flux can be assumed to be applied
2441
Fig. 1. A disc brake assembly; (a) real assembly (picture taken from ref. [4]), (b) model configuration.
on the entire circumferential domain, and the problem is further reduced to an axial-symmetric problem where temperature is a function of r, z, and t only. This approach has been widely adopted in analytical studies [24] and has been proved by experimental investigations [25]. To illustrate the methodology for developing expressions for the use in determining the unknown space- and time-dependent frictional heat flux at the discepad interface in a disc brake system, we consider the following heat transfer problem of the disc as shown in Fig. 2. Under the condition of axisymmetry, the governing equations and the associated boundary and initial conditions for the disc can be written as [4,24]:
1 v vTðr; z; tÞ v vTðr; z; tÞ kðTÞr þ kðTÞ r vr vr vz vz vTðr; z; tÞ ; in r1 < r < r3 ; 0 < z < z2 ; t > 0; vt
(1)
vTðr;z;tÞ ¼ h½TN Tðr;z;tÞ; at r ¼ r1 ;0 < z < z2 ;t > 0; vr
(2)
¼ rcðTÞ kðTÞ
2442
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
The direct problem considered here is concerned with the determination of the medium temperature when the frictional heat flux q(r,t), thermo-physical properties of the disc, and initial and boundary conditions are known. 2.2. Inverse problem For the inverse problem, the function q(r,t) is regarded as being unknown, while everything else in Eqs. (1)e(7) is known. In addition, temperature readings taken within the disc are considered available. The objective of the inverse analysis is to predict the unknown space- and time-dependent function of frictional heat flux, q(r,t), merely from the knowledge of these temperature readings. Let the measured temperature at the measurement position ðr; zÞ ¼ ðri ; zm Þ and time t be denoted by Yi ðri ; zm ; tÞ, i ¼ 1eM, where M represents the number of thermocouples and zm is the z coordinate of all the thermocouples. Then this inverse problem can be stated as follows: by utilizing the above mentioned measured temperature data Yi ðri ; zm ; tÞ, the unknown q(r,t) is to be estimated over the specified space and temporal domain. The solution of the present inverse problem is to be obtained in such a way that the following functional is minimized:
J½qðr; tÞ ¼
Ztf X M t¼0
½Ti ðri ; zm ; tÞ Yi ðri ; zm ; tÞ2 dt;
(8)
i¼1
where Ti ðri ; zm ; tÞ is the estimated (or computed) temperature at the measurement location ðr; zÞ ¼ ðri ; zm Þ. In this study, Ti ðri ; zm ; tÞ are determined from the solution of the direct problem given previously by using an estimated qK ðr; tÞ for the exact qðr; tÞ, here qK ðr; tÞ denotes the estimated quantities at the Kth iteration. tf is the final time of measurement. In addition, in order to develop expressions for the determination of the unknown q(r,t), a “sensitivity problem” and an “adjoint problem” are constructed as described below. 2.3. Sensitivity problem and search step size Fig. 2. Schematic of the disc.
kðTÞ
vTðr;z;tÞ ¼ h½TN Tðr;z;tÞ; at r ¼ r3 ;0 < z < z2 ;t > 0; vr
(3)
kðTÞ
vTðr; z; tÞ ¼ 0 at z ¼ 0; r1 < r < r3 ; t > 0; vz
(4)
kðTÞ
vTðr;z;tÞ ¼ h½TN Tðr;z;tÞ; at z ¼ z2 ;r1 < r < r2 ;t > 0; vz
(5)
vTðr; z; tÞ ¼ qðr; tÞ at z ¼ z2 ; r2 < r < r3 ; t > 0; kðTÞ vz
(6)
Tðr; z; 0Þ ¼ T0 ; for r1 < r < r3 ; 0 < z < z2 ; t ¼ 0;
(7)
where TN is the surrounding temperature and T0 is the initial temperature. r is the density, k(T) and c(T)are the temperaturedependent thermal conductivity and specific heat, respectively, while h is the convection heat transfer coefficient at the outer surfaces. q(r,t) in Eq. (6) is the frictional heat flux into the disc. If the contact pressure between the pad and disc is assumed constant, the frictional heat flux is a function of time and space variable r. This is because work done by friction force grows as r increases. However, if uniform wear contact is adopted, the frictional heat flux is just a function of time and is independent of the space variable r.
The sensitivity problem is obtained from the original direct problem defined by Eqs. (1)e(7) in the following manner: It is assumed that when q(r,t) undergoes a variation q(r,t), T(r,z,t), k(T), and c(T) are perturbed by DT, Dk, and Dc, respectively. Then replacing in the direct problem q(r,t) by qðr; tÞ þ Dqðr; tÞ, T by T þ DT, k(T) by k(T) þ Dk, and c(T) by c(T) þ Dc, subtracting from the resulting expressions the direct problem, using the fact that Dk ¼ ½dkðTÞ=dTDT and Dc ¼ ½dcðTÞ=dTDT, and neglecting the second-order terms, the following sensitivity problem for the sensitivity function DT can be obtained.
1 v v v2 r ½kðTÞDTðr; z; tÞ þ 2 ½kðTÞDTðr; z; tÞ r vr vr vz v ¼ ½rcðTÞDTðr; z; tÞ; in r1 < r < r3 ; 0 < z < z2 ; t > 0; vt
(9)
v ½kðTÞDTðr;z;tÞ ¼ hDTðr;z;tÞ; at r ¼ r1 ;0 < z < z2 ;t > 0; vr
(10)
v ½kðTÞDTðr;z;tÞ ¼ hDTðr;z;tÞ; at r ¼ r3 ;0 < z < z2 ;t > 0; vr
(11)
v ½kðTÞDTðr; z; tÞ ¼ 0; at z ¼ 0; r1 < r < r3 ; t > 0; vz
(12)
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
2443
v ½kðTÞDTðr;z;tÞ ¼ hDTðr;z;tÞ; at z ¼ z2 ;r1 < r < r2 ;t > 0; vz
(13)
kðTÞ
vlðr; z; tÞ ¼ hlðr; z; tÞ; at r ¼ r1 ; 0 < z < z2 ; t > 0; vr
(19)
v ½kðTÞDTðr;z;tÞ ¼ Dqðr;tÞ; at z ¼ z2 ;r2 < r < r3 ;t > 0; vz
(14)
kðTÞ
vlðr; z; tÞ ¼ hlðr; z; tÞ at r ¼ r3 ; 0 < z < z2 ; t > 0; vr
(20)
DTðr; z; 0Þ ¼ 0; for r1 < r < r3 ; 0 < z < z2 ; t ¼ 0:
(15)
kðTÞ
vlðr; z; tÞ ¼ 0 at z ¼ 0; r1 < r < r3 ; t > 0; vz
(21)
The sensitivity problem of Eqs. (9)e(15) can be solved by the same method as the direct problem of Eqs. (1)e(7).
kðTÞ
vlðr; z; tÞ ¼ hlðr; z; tÞ at z ¼ z2 ; r1 < r < r2 ; t > 0; vz
(22)
2.4. Adjoint problem and gradient equation
kðTÞ
vlðr; z; tÞ ¼ 0 at z ¼ z2 ; r2 < r < r3 ; t > 0; vz
(23)
To formulate the adjoint problem, Eq. (1) is multiplied by the Lagrange multiplier (or adjoint function) l, and the resulting expressions are integrated over the time and correspondent space domains. Then the results are added to the right hand side of Eq. (8) to yield the following expression for the functional J½qðr; tÞ:
Zz2
Ztf
Zr3 X M
J½qðr; tÞ ¼
½Ti ðr; z; tÞ
i¼1
t¼0 z¼0 r¼r1
lðr; z; tÞ ¼ 0 for r1 < r < r3 ; 0 < z < z2 ; t ¼ tf :
The adjoint problem is different from the standard initial value problem in that the final time condition at time t ¼ tf is specified instead of the customary initial condition at time t ¼ 0. However, this problem can be transformed to an initial value problem by the transformation of the time variable as s ¼ tf t. Then the adjoint problem can be solved by the same method as the direct problem. Finally the following integral term is left:
Yi ðr; z; tÞ2 dðr ri Þdðz zm Þdrdzdt Ztf
Zz2
DJ ¼
1 v vT r$lðr; z; tÞ$ kðTÞr r vr vr
Zr3
þ t¼0 z¼0 r¼r1
Zz2
Zr3 X M
t¼0 z¼0 r¼r1
(16)
DJ ¼
i¼1
Zr3
Ztf
rJ 0 ðr; tÞDqðr; tÞdrdt;
(26)
t¼0 r¼r2
where J0 (t) is the gradient of the functional J(q). A comparison of Eqs. (25) and (26) leads to the following form:
J 0 ðr; tÞ ¼ lðr; z2 ; tÞ:
(27)
The following iteration process based on the conjugate gradient method is now used for the estimation of q(r,t) by minimizing the above functional J[q(r,t)]:
Yi ðr; z; tÞDT dðr ri Þdðz zm Þdrdzdt ( Ztf Zz2 Zr3 v v l$ þ $r$ ½kðTÞDT vr vr
K qKþ1 ðr; tÞ ¼ qK ðr; tÞ b pK r; t ;
)
v2 v þ r$ 2 ½kðTÞDT r$ ½rcðTÞDT drdzdt vt vz
(17)
We can integrate the second triple-integral term in Eq. (17) by parts, utilizing the initial and boundary conditions of the sensitivity problem. Then DJ is allowed to go to zero. The vanishing of the integrands containing DT leads to the following adjoint problem for the determination of lðr; z; tÞ:
pK ðr; tÞ ¼ J 0 Kðr; tÞ þ gK pK1 r; t ;
gK ¼ (18)
(28)
(29)
which is conjugation of the gradient direction J 0 Kðr; tÞ at iteration K and the direction of descent pK1 ðr; tÞ at iteration K 1. The conjugate coefficient gK is determined from:
Ztf X M h
M 1 X þ $ 2½Ti ðr; z; tÞ Yi ðr; z; tÞdðr ri Þdðz zm Þ r i¼1
K ¼ 0; 1; 2; .;
where bK is the search step size in going from iteration K to iteration Kþ1, and pK(r,t) is the direction of descent (i.e., search direction) given by:
1 v vlðr; z; tÞ v2 lðr; z; tÞ vlðr; z; tÞ þ rcðTÞ $kðTÞ$ r þ kðTÞ$ r vr vr vt vz2
¼ 0; in r1 < r < r3 ; 0 < z < z2 ; t > 0;
(25)
2.5. Conjugate gradient method for minimization
2½Ti ðr; z; tÞ
t¼0 z¼0 r¼r1
r lðr; z2 ; tÞDqðr; tÞdrdt:
t¼0 r¼r2
where dð,Þ is the Dirac function. The variation DJ is derived by perturbing q(r,t) by qðr; tÞ þ Dqðr; tÞ, T by T þ DT, k(T) by k(T) þ Dk, and c(T) by c(T) þ Dc, respectively, in Eq. (16). Subtracting from the resulting expressions the original Eq. (16), using the fact that Dk ¼ ½dkðTÞ=dTDT and Dc ¼ ½dcðTÞ=dTDT, and neglecting the secondorder terms, we thus find:
Ztf
Zr3
Ztf
From the definition used in Ref. [7], we have:
v vT vT kðTÞ rcðTÞ drdzdt þ vz vz vt
DJ½qðr; tÞ ¼
(24)
t¼0
i2 J 0 Kðri ; tÞ dt
i¼1
Ztf X M h i2 J 0 K 1ðri ; tÞ dt t¼0
i¼1
with g0 ¼ 0:
(30)
2444
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
The convergence of the above iterative procedure in minimizing the functional J is proved in Ref. [7]. To perform the iterations according to Eq. (28), we need to compute the step size bK and the gradient of the functional J 0 Kðr; tÞ. The functional J½qKþ1 ðr; tÞ for iteration K þ 1 is obtained by rewriting Eq. (8) as:
Ztf X M h h i i2 K Ti qK b pK Yi ðri ;zm ;tÞ dt; J qKþ1 ðr;tÞ ¼
Inner disc diameter Outer disc diameter Mean sliding radius Sliding length Disc thickness Pad thickness Initial speed Braking time Deceleration Total energy Convection coefficient
(31)
i¼1
t¼0
where we replace qKþ1 by the expression given by Eq. (28). If K temperature Ti ðqK b pK Þ is linearized by a Taylor expansion, Eq. (31) takes the form:
h i J qKþ1 ðr; tÞ ¼
Table 1 Brake dimensions and operating conditions.
Ztf t¼0
M h X
K Ti ðqK Þ b DTi ðpK Þ
i¼1
i2 Yi ðri ; zm ; tÞ dtdt;
ð32Þ
properties are taken as constant. Furthermore, under the assumption of constant deceleration and uniform wear contact between the pad and the disc, the boundary heat flux in Eq. (6) is a function of time only and can be written as:
Ztf X M
bK ¼
t¼0
DTi pK
h i Ti qK Yi dt
i¼1
Z tf X M h
i2 DTi ðp Þ dt
t¼0
:
(33)
2.6. Stopping criterion If the problem contains no measurement errors, the traditional check condition specified as
(34)
where h is a small specified number, can be used as the stopping criterion. However, the observed temperature data contains measurement errors; as a result, the inverse solution will tend to approach the perturbed input data, and the solution will exhibit oscillatory behavior as the number of iteration is increased [26]. Computational experience has shown that it is advisable to use the discrepancy principle [23] for terminating the iteration process in the conjugate gradient method. Assuming Ti ðri ; zm ; tÞ Yi ðri ; zm ; tÞys, the stopping criteria h by the discrepancy principle can be obtained from Eq. (8) as
h ¼ Ms2 tf ;
(36)
mfpmax r2 u0 ;
(37)
with
q0 ¼
q0 2p
where pmax is the maximum pressure distributed in the pad and f is the heat partition coefficient between the disc and the pad defined as:
f¼
x d Sd ; xd Sd þ xp Sp
(38)
K
i¼1
J qKþ1 < h;
! t 1 ; tf
qðtÞ ¼ q0
ðqK Þ
where Ti is the solution of the direct problem at ðr; zÞ ¼ ðri ; zm Þ by using estimated qK(r,t) for exact q(r,t) at time t. The sensitivity function DTi ðpK Þ is taken as the solution of Eqs. (9)e(15) at the measured position ðr; zÞ ¼ ðri ; zm Þ by letting Dq ¼ pK [7]. The search step size bK is determined by minimizing the functional given by Eq. (32) with respect to bK. The following expression can be obtained:
0.132 (m) 0.227 (m) 0.094.5 (m) 0.037 (m) 0.011 (m) 0.010 (m) 100 (km h1) 3.96 (s) 7 (ms2) 165 (kJ) 60 (W m2 K1)
(35)
where s is the standard deviation of the measurement error. Then the stopping criterion is given by Eq. (34) with h determined from Eq. (35). The computational procedure is the same as in [20], hence there is no need to repeat here. 3. Results and discussion A bench-mark test case in [24] is first considered to validate the code used in this study. The details of the dimensions and operating conditions are given in Tables 1 and 2. In this test case, the thermal
where xd and xP are the thermal effusivities of the disc and the pad, and and Sd and Sp are frictional contact surfaces of the disc pffiffiffiffiffiffiffi ffi the pad, respectively. Thermal effusivity x is defined as x ¼ krc. The numerical procedure in this paper is based on the unstructured-mesh, fully collocated, finite-volume code, ‘USTREAM’ developed by W.-L. Chen. This is the descendent of the structuredmesh, multi-block code of ‘STREAM’ [27]. The computational mesh consists of 860 cells, where there are 20 and 43 cells, respectively, in the z- and r-directions. This mesh has been proved in a gridindependent test to achieve grid-independent solutions. The comparison of mean sliding radius surface temperature variations is given in Fig. 3. To find an appropriate time-step interval to achieve time-step-interval independent solutions, the results are computed by using three different numbers of time steps, namely 100, 200, and 300. It can be seen that the results by the three time-step numbers are almost identical and are in very good agreement with Newcomb’s analytical solution. This indicates that the time step number of 100, corresponding to a time-step interval of 0.0396 s, is good enough to ensure time-step-interval independent solutions in this problem. In addition, the test also versifies the correctness of the code used in this study. To test the effect of adopting temperature-dependent thermophysical properties, the disc is assumed made of AISI-1045 steel, which is commonly used in friction elements. Since the density of AISI-1045 steel has relatively weak temperature dependence Table 2 Thermo-physical properties of disc and pad.
Conductivity (W m1 K1) Density (kg m3) Specific heat (J kg1 K1)
Disc
Pad
43.5 7850 445
12 2500 900
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
2445
a
Fig. 3. Comparison disc mean sliding radius surface temperature variations.
over a wide temperature range [28], a constant average value of 7844 kg m3 is adopted. However, the specific heat and thermal diffusivity from the experimental measurements of Touloukian [29] are used here, and Fig. 4 shows the variations of specific heat and thermal conductivity with respect to temperature from 280 K to 1300 K. The mean sliding radius temperature variations computed with temperature-dependent thermo-physical properties are also shown in Fig. 3. Although the results are not significantly different from those adopting constant thermo-physical properties, the deviation, about 5 K at most places, is clearly present. To ensure high level of accuracy, temperature-dependent thermo-physical properties should be used, and in the following, all results were obtained using temperature-dependent thermo-physical properties. The objective of this article is to validate the current approach when used in estimating the unknown heat flux into the disc during a braking process accurately with no prior information on the functional form of the unknown quantities, a procedure called function estimation. In this study, we did not have an experiment setup to measure the temperatures on the disc. Instead, we substitute an assumed interface heat flux such as Eq. (36) into the direct problem of Eqs. (1)e(7) to calculate the temperatures at the locations where the temperature measurements are taken. The results are termed the computed temperature Yexact ðri ; zm ; tÞ. However, in reality, temperature measurements always contain some degree of error, whose magnitude depends upon the particular measuring method employed. To take measurement errors into account, a random error noise is added to the above computed temperature Yexact ðri ; zm ; tÞ to obtain the measured temperature Yi ðri ; zm ; tÞ. Hence, the measured temperature Yi ðri ; zm ; tÞ is expressed as:
Yi ðri ; zm ; tÞ ¼ Yexact ðri ; zm ; tÞ þ 6s;
(39)
where 6 is a random variable within 2.576 to 2.576 for a 99% confidence bounds, and s is the standard deviation of the measurement. The measured temperature Yi ðri ; zm ; tÞ generated in such way is called the simulated measurement temperature. At first, the temperature measurement is assumed located at zm ¼ 0.0054 m, equivalent to 0.1 mm inside the disc, and along the friction interface (see Fig. 2). The comparison between the inverse and exact heat flux, using Eq. (36), at three different time steps is given in Fig. 5. The inverse results were obtained with an initial guess of q0(r,t) ¼ 1.0 105 W m2 and measurement errors s ¼ 0.0.
b
Fig. 4. Measured AISI-1045 steel specific heat and thermal conductivity variations against temperature. (a) Specific heat variations, and (b) thermal conductivity variations.
Fig. 5 shows that for this simple heat flux function, the inverse estimation is almost identical to the exact value. Next, we adopt the uniform pressure assumption in Talati and Jalalifar [4], where the heat flux function is written as:
!
r t 1 ; qðr; tÞ ¼ q0 r2 tf
(40)
Eq. (40) indicates that the heat flux is now a function of radius and time. The results at three different time steps are shown in Fig. 6, and yet again very good agreement between the inverse and the exact values can be seen. Here the inverse results were obtained with an initial guess of q0(r,t) ¼ 1.0 105 W m2 and measurement errors s ¼ 0.0. Once the heat flux has been estimated correctly, the temperature evolution anywhere on the disc can be calculated accurately too. This is shown in Fig. 7, where the comparison of temperature evolution at four different locations on the disc surface is given. It is noticeable that the variations of calculated and exact temperature are almost identical.
2446
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
Fig. 5. Estimated q(t) (Eq. (37)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
Fig. 7. Estimated disc surface temperature variations at 4 different radial positions with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
In order to investigate the influence of measurement error upon the estimated results, Fig. 8 shows the estimated values of the unknown function q(r,t), obtained with the initial guesses q0(r,t) ¼ 1.0 105 W m2, and measurement error of deviation s ¼ 0.005. The deviation s ¼ 0.005 corresponds to 1.28% of measurement error, which is roughly in the range of 4e6 K over the entire disc surface temperature variation during the braking process. This magnitude of temperature measurement error is huge according to the accuracy standard of modern thermal couples, whose precision is normally within 1 K. Despite this unrealistically large measurement error, the results indicate that this temperature measurement error in general creates only about 2% of error on heat flux magnitude, demonstrating that the accuracy of the current method is not very sensitive to temperature measurement error. Since the direct problem is a transient heat transfer problem, there have been studies reporting the sensitivity of inverse accuracy on the relative distance between the estimated and measurement quantities in inverse transient heat conduction problems [30]. Under some
circumstances, the accuracy deteriorates rapidly as the measurement location moves away from the estimated quantity location. Therefore in the current inverse problem, this issue needs to be investigated. To examine the sensitivity of inverse method’s accuracy on measurement location, the thermal couple is assumed to be located respectively at three other positions, zm ¼ 0.005 m, 0.0035 m, and 0.0001 m. The estimated heat flux distributions, together with those of zm ¼ 0.0054 m, are plotted in Fig. 9. These results were all obtained with an initial guess of q0(r,t) ¼ 1.0 105 W m2 and measurement errors s ¼ 0.0, and they obviously show that the inverse accuracy deteriorates progressively as the measurement position moves away from the friction interface, especially at the vicinity of r ¼ r2. In the case of zm ¼ 0.0001 m, the inverse heat flux at t ¼ 3.66 s is in general half the magnitude of the exact value, resulting in an discrepancy even larger than that produced by taking measurement error into account.
Fig. 6. Estimated q(r,t) (Eq. (40)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
Fig. 8. Estimated q(r,t) (Eq. (40)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.005.
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
Fig. 9. Estimated q(r,t) (Eq. (40)) with four different measurement locations, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
2447
Fig. 11. Estimated q(r,t) (Eq. (41)) with zm ¼ 0.0054 m, the initial guess q0(r,t) ¼ 1.0 105 W m2, and measurement error s ¼ 0.0.
a In order to demonstrate the capability of the present method in obtaining an accurate estimation of arbitrary form of the unknown heat flux, we consider another case of q(r,t), based on the thermoelastic simulation in Choi and Lee [2], which is calculated by:
qðr; tÞ ¼ mr uðtÞpðr; tÞ;
b
(41)
where u(t) and p(r,t) are distributions taken from Figs. 4 and 6 in [2]. The heat generation during the braking stage is shown in Fig. 10(a). Please note that due to data digitizing and scaling factors, the absolute magnitude of heat flux shown here is not identical to that in [2] but the tendency of data variation is preserved. In their study, heat conduction and elastic equations are coupled through the heat generated at discepad interface calculated by Eq. (41). Therefore, the current inverse method can be applied to thermoelastic approach to estimate the heat flux into the disc and then evaluate the heat partition between the disc and pad at the friction interface. Fig. 10(b) shows the estimated results of q(r,t), obtained with zm ¼ 0.0054 m, the initial guesses q0(r,t) ¼ 1.0 105 W m2, and measurement error of deviation s ¼ 0.0. Excellent agreement between the exact and inverse heat flux can be seen by comparing Fig. 10(a) and (b). More detailed comparison on the heat flux distributions along the radial direction at three different time steps given in Fig. 11 further confirms the accuracy of the inverse estimation. The results indicate that the current inverse method performs very well even on this highly complicated heat flux function. 4. Conclusion
Fig. 10. The heat flux function taken from Choi and Lee [2]; (a) exact q(r,t) distributions, and (b) inverse q(r,t) distributions.
An inverse algorithm based on the conjugate gradient method and the discrepancy principle was successfully applied for the solution of the inverse problem to determine the unknown space- and time-dependent heat flux into the disc during a braking process, while knowing the temperature history at some measurement locations. Due to the wide range of temperature variation in a braking process, temperature-dependent specific heat and thermal conductivity are used here; hence the governing equation of the direct problem becomes nonlinear. Numerical results confirm that the method proposed herein can accurately estimate the space- and time-dependent heat flux and temperature distributions for the problem even involving the inevitable measurement errors. The
2448
Y.-C. Yang, W.-L. Chen / Applied Thermal Engineering 31 (2011) 2439e2448
proposed inverse algorithm does not require prior information for the functional form of the unknown quantities to perform the inverse calculation, and excellent estimated values can be obtained for the considered problem. However, the results also show that the inverse accuracy is very sensitive to the distance between the estimated and measurement quantities. The temperature measurement location should be put as close to the friction interface as possible because the accuracy of the inverse method deteriorates rapidly as the measurement location moves away from the friction interface. Acknowledgements This work was supported by the National Science Council, Taiwan, Republic of China, under the grant numbers NSC 98-2221E-168-035 and NSC 98-2221-E-168-034. References [1] P. Dufrenoy, Two-/three-dimensional hybrid model of the thermomechanical behavior of disc brakes, Proc. IMechE Part F: J. Rail Rapid Transit. 218 (2004) 17e30. [2] J.H. Choi, I. Lee, Finite element analysis of transient thermoelastic behaviors in disc brakes, Wear 257 (2004) 47e58. [3] C.H. Gao, X.Z. Lin, Transient temperature field analysis of a brake in a nonaxisymmetric three-dimensional model, J. Mater. Process. Technol. 129 (2002) 513e517. [4] F. Talati, S. Jalalifar, Analysis of heat conduction in a disk brake system, Heat Mass Transfer 45 (2009) 1047e1059. [5] N. Laraqi, Phenomene de constriction thermique dans les contacts glissants (thermal constriction phenomenon in sliding contacts), Int. J. Heat Mass Transfer 39 (1996) 3717e3724. [6] D. Majcherczak, P. Dufrenoy, M. Nait-Abdelaziz, Third body influence on thermal friction contact problems: application to braking, ASME J. Tribol. 127 (2005) 89e95. [7] O.M. Alifanov, Inverse Heat Transfer Problem. Springer-Verlag, New York, 1994. [8] C.H. Huang, C.C. Shih, A shape identification problem in estimating simultaneously two interfacial configurations in a multiple region domain, Appl. Therm. Eng. 26 (2006) 77e88. [9] Y.C. Yang, Simultaneously estimating the contact heat and mass transfer coefficients in a double-layer hollow cylinder with interface resistance, Appl. Therm. Eng. 27 (2007) 501e508. [10] H.T. Chen, H.C. Wang, Estimation of heat-transfer characteristics on a rectangular fin under wet conditions, Int. J. Heat Mass Transfer 51 (2008) 2123e2138. [11] W.L. Chen, Y.C. Yang, S.S. Chu, Estimation of heat generation at the interface of cylindrical bars during friction process, Appl. Therm. Eng. 29 (2009) 351e357.
[12] Y.C. Yang, W.L. Chen, An iterative regularization method in simultaneously estimating the inlet temperature and heat-transfer rate in a forced-convection pipe, Int. J. Heat Mass Transfer 52 (2009) 1928e1937. [13] W.L. Chen, Y.C. Yang, Simultaneous estimation of heat-transfer rate and coolant fluid velocity in a transpiration cooling process, Int. J. Therm. Sci. 49 (2010) 1407e1416. [14] C. Le Niliot, F. Lefevre, Multiple transient point heat sources identification in heat diffusion: application to numerical two- and three-dimensional problems, Numer. Heat Transfer Part B 39 (2001) 277e301. [15] M. Raudensky, K.A. Woodbury, J. Kral, Genetic algorithm in solution of inverse heat conduction problems, Numer. Heat Transfer Part B 28 (1995) 293e306. [16] C. Le Niliot, F. Lefevre, A parameter estimation approach to solve the inverse problem of point heat sources identification, Int. J. Heat Mass Transfer 47 (2004) 827e841. [17] C.H. Huang, H.C. Lo, A three-dimensional inverse problem in estimating the internal heat flux of housing for high speed motors, Appl. Therm. Eng. 26 (2006) 1515e1529. [18] H.L. Lee, W.J. Chang, W.L. Chen, Y.C. Yang, An inverse problem of estimating the heat source in tapered optical fibers for scanning near-field optical microscopy, Ultramicroscopy 107 (2007) 656e662. [19] B. Jin, L. Marin, The method of fundamental solutions for inverse source problems associated with the steady-state heat conduction, Int. J. Numer. Methods Eng. 69 (2007) 1570e1589. [20] W.L. Chen, Y.C. Yang, Inverse problem of estimating the heat flux at the roller/ workpiece interface during a rolling process, Appl. Therm. Eng. 30 (2010) 1247e1254. [21] J. Thevenet, M. Siroux, B. Desmet, Measurements of brake disc surface temperature and emissivity by two-color pyrometry, Appl. Therm. Eng. 30 (2010) 753e759. [22] P. Litos, M. Honner, V. Lang, J. Bartik, M. Hynek, A measuring system for experimental research on the thermomechanical coupling of disc brakes, Proc. IMechE, Part D: J. Automobile Eng. 222 (2008) 1247e1257. [23] O.M. Alifanov, E.A. Artyukhin, Regularized numerical solution of nonlinear inverse heat-conduction problem, J. Eng. Phys. 29 (1975) 934e938. [24] T.P. Newcomb, Transient temperatures attained in disk brakes, Br. J. Appl. Phys. 10 (1959) 339e340. [25] D. Majcherczak, P. Dufrenoy, Y. Berthier, Tribological, thermal and mechanical coupling aspects of the dry sliding contact, Tribol. Int. 40 (2007) 834e843. [26] O.M. Alifanov, Application of the regularization principle to the formulation of approximate solution of inverse heat conduction problem, J. Eng. Phys. 23 (1972) 1566e1571. [27] F.S. Lien, W.L. Chen, M.A. Leschziner, A multiblock implementation of a nonorthogonal, collocated finite volume algorithm for complex turbulent flows, Int. J. Numer. Methods Fluids 23 (1996) 567e588. [28] R.D. Pelhlke, A. Jeyarajan, H. Wada, Report No. NSF/MAE-82028, NSF Applied Research Division, University of Michigan (1982). [29] Y.S. Touloukian, Thermophysical Properties of High Temperature Solid Materials. Macmillan, New York, 1967. [30] W.L. Chen, Y.C. Yang, Inverse prediction of frictional heat flux and temperature in sliding contact with a protective strip by iterative regularization method, Appl. Math. Model. 35 (2011) 2874e2886.