Parameter Estimation of the Thermal Network Model of a Machine Tool Spindle by Self-made Bluetooth Temperature Sensor Module

Thermal characteristic analysis is essential for machine tool spindles because sudden failures may occur due to unexpected thermal issue. This article presents a lumped-parameter Thermal Network Model (TNM) and its parameter estimation scheme, including hardware and software, in order to characterize both the steady-state and transient thermal behavior of machine tool spindles. For the hardware, the authors develop a Bluetooth Temperature Sensor Module (BTSM) which accompanying with three types of temperature-sensing probes (magnetic, screw, and probe). Its specification, through experimental test, achieves to the precision ±(0.1 + 0.0029|t|) °C, resolution 0.00489 °C, power consumption 7 mW, and size Ø40 mm × 27 mm. For the software, the heat transfer characteristics of the machine tool spindle correlative to rotating speed are derived based on the theory of heat transfer and empirical formula. The predictive TNM of spindles was developed by grey-box estimation and experimental results. Even under such complicated operating conditions as various speeds and different initial conditions, the experiments validate that the present modeling methodology provides a robust and reliable tool for the temperature prediction with normalized mean square error of 99.5% agreement, and the present approach is transferable to the other spindles with a similar structure. For realizing the edge computing in smart manufacturing, a reduced-order TNM is constructed by Model Order Reduction (MOR) technique and implemented into the real-time embedded system.


Introduction
The thermal characteristics of machine tool spindles play an important role in the development of high-precision and high-speed machining. Many literatures about the machining error analysis had indicated that a large amount of the total machining errors attributed to the thermal issue [1,2]. In recent years, diagnostics and prognostics of sudden failure of machine tool spindle have drawn considerable attention. Hence, the model-based method is the key technology for predicting the thermal behavior of spindle. Bossmanns and Tu et al. developed a comprehensive thermo-mechanical model to characterize the heat transfer and quantitative power flow of a high-speed motorized spindle based on Finite Difference Method (FDM) [3][4][5]. On the other hand, Finite Element Method (FEM) has the advantage of fully simulating the thermo-dynamic-mechanical behavior of the spindle. These FEM models not only simulate the temperature distribution and the thermal deformation of machine tool, but also take the bearing stiffness and thermal contact resistance into consideration [6,7]. Experimental and simulation For angular contact ball bearing, the gyroscopic moment of the rolling elements necessarily leads to boring motion that would cause the friction torque on the contact surface, and thus also introduce heat generation. The friction torque related to spinning motion (M s ) is obtained by Equation (4) [25].
where F roll is the contact load, a represent the major axes of the elliptical contact area, e is the elliptical integral of the second kind, and µ s is the friction coefficient. The primary concern is to evaluate the heat source within the spindle thermal system associated with the spindle rotational speed (ω). Many researchers had indicated that the viscous friction has the influential contribution in total frictional torque within the bearing [12]. As for simplification, due to the power function trend of thermally induced preload influence, the polynomial-like function might have the ability to predict the heat generation of bearing. As mentioned above, substituting Equations (2)-(4) into Equation (1), the authors conclude that the heat source can be expressed as a function of the viscous-related term (ω 5/3 ) and the remaining term (ω), as described in Equation (5).
An analytical solution describes the forced convective heat transfer over rotating shaft [26]. The average Nusselt number (Nu), Reynolds number (Re), and Prandtl number (Pr) are determined by the following equations. Nu = 0.6366(RePr) 1/2 (6) Re = ωd 2 2ν air ,Pr = ν air α air Assume constant material properties, the rotational speed dominates the forced heat convection. This implies that the forced convection coefficient is a speed-related term being proportional to ω 0.5 . In this case, we take the rotational speed as the leading term and express the thermal resistance of forced convection around moving surface as a function of ω −0. 5 .
The forced heat convection over the surfaces of rotating annulus can be determined by the following Equations (10) [27]. Taylor number (Ta) and the geometry factor (F g , P) are obtained. Note that the geometry factors (Fg and P) depend on the inner and outer radius of annulus only; therefore, it can be viewed as a constant value.
, for 10 4 < Ta 2 F g < 10 7 (10) where and r i and r o are the inner and outer radius. Namely, the coefficient of forced heat convection can be derived in Equation (11), which is proportional to ω 0.482 . Thus, it is appropriate to formulate the thermal resistance of forced convection near annulus as a function of ω −0.482 in Equation (12).
h annulus ∝ Nu ∝ Ta 0.482 ∝ ω 0.482 (11) R annulus = 1 h annulus A annulus ∝ ω −0.482 (12) According to [28], free convection around a horizontal cylinder and stationary ambient can be estimated by the following Equation (13 Except for the material properties and geometry factor, the free convection coefficient is dependent on the spindle surface temperature (T s ) and ambient temperature (T a ) particularly. It would be difficult to simulate the free convective phenomenon according to the temperature difference. Alternatively, the rotational speed has strong correlation to the temperature difference. We assume that the free convective coefficient can be expressed as the function in Equation (14). Hence, the free convective thermal resistance can be derived as two terms, a constant term and a speed-related term. Especially, all of the assumptions must be verified by experimental results, and it will be further discussed in Section 5.
In general, the radiative thermal resistance can be expressed as Equation (15). It indicates that the resistance of thermal radiation depends on surface temperature (T s ) and ambient temperature (T a ). Due to the extremely tiny value of Stefan-Boltzmann constant (σ) and the emissivity (ε), the influence of surface temperature variation would be limited. Consequently, the radiative thermal resistance can be considered as a constant value, which would be further verified by experiment.
The thermal capacitance is directly determined by the element properties (volume, heat capacity, density) through Equation (16). It represents the store of energy that is required to increase the temperature. In other words, the thermal capacitance value is determined once the spindle structure is designed in practice. Thereby, it is assumed as a constant, when considering that the material properties remain relatively the same under different operating conditions. The estimated value of each thermal capacitance has the decisive influence on the thermal time constant.

System Identification Technique
The TNM of machine tool spindle is identified by system identification methodology based on physical concept, namely the so-called grey-box estimation. It is utilized to predict the temperature distribution of spindle. Figure 1 explains the system identification procedure for the TNM of machine tool spindle comprehensively. First, through the temperature data measured at steady state, solve for the thermal resistances of TNM at operating mode by means of the "lsqnonlin" function and the "trust-region-reflective" algorithm; the aforesaid function and algorithm are the nonlinear least square regression provided by MATLAB [29]. Second, solve for the thermal capacitance of TNM by means of the "greyest" function provided by MATLAB; that function is exactly the grey-box model identification [29]. Finally, the thermal resistances at natural cooling mode are obtained by the grey-box model identification as well. The least square system identification problem usually exists non-unique solutions satisfying the cost function (minimization). It is expected that some irrational solutions can be eliminated based on physical sense, for example, the convective thermal resistance is generally greater than the conductive one. Furthermore, it must be emphasized that while inversely estimate the heat generation, thermal resistance, and thermal capacitance at the same time, some analogous solutions might occur because different solutions with the same ratio among all of the parameters might introduce the same result. Fortunately, the total thermal capacitance (∑ 4 j=1 C j ) is a constant value; thus, the appropriate solution can be determined accordingly. Finally, it is necessary to verify the performance of the estimated TNM; hence, the self-validation and external validation are further investigated in Section 5.

Parameterization Strategy
The assumptions for modeling the TNM are listed as follows: 1) lumped element model assumption is valid, Bi << 0.1; 2) temperature distribution of spindle is axisymmetric; 3) heat generation and forced convective resistance are assumed according to the theoretical and empirical Equations (5), (9), and (12); 4) heat transfer through radiation and conduction are considered as a constant thermal resistance, including thermal contact resistance; and, 5) free convection coefficient is assumed as a function of ℎ , ℎ , According to previous research, the temperature distribution of spindle is almost axisymmetric. For this reason, the two-dimensional heat equation with internal heat generation can be implemented to describe the three-dimensional thermal behavior. Furthermore, based on the assumption of lumped element model, the TNM can characterize the heat transfer inside the spindle.

Estimated Thermal Network Model in State-Space
While the spindle operates at a constant rotational speed, the frictional heat generation are assumed to be a constant source, denoted as Qf1, Qf2. As a result, the heat balance equation with internal heat generation is expressed as

Parameterization Strategy
The assumptions for modeling the TNM are listed as follows: (1) lumped element model assumption is valid, Bi << 0.1; (2) temperature distribution of spindle is axisymmetric; (3) heat generation and forced convective resistance are assumed according to the theoretical and empirical Equations (5), (9), and (12); (4) heat transfer through radiation and conduction are considered as a constant thermal resistance, including thermal contact resistance; and, (5) free convection coefficient is assumed as a function of h f ree,c + h f ree,s ω −0.5 .
According to previous research, the temperature distribution of spindle is almost axisymmetric. For this reason, the two-dimensional heat equation with internal heat generation can be implemented to describe the three-dimensional thermal behavior. Furthermore, based on the assumption of lumped element model, the TNM can characterize the heat transfer inside the spindle.

Estimated Thermal Network Model in State-Space
While the spindle operates at a constant rotational speed, the frictional heat generation are assumed to be a constant source, denoted as Q f1 , Q f2 . As a result, the heat balance equation with internal heat generation is expressed as The authors define the temperature difference (θ) with respect to the ambient temperature.
Substituting it into Equation (17) gives the heat balance equations Consider the physical insight and geometry of the spindle; we construct the lumped-parameter TNM of the spindle with four nodes, illustrated in Figure 2, which are rear bearing A (T 1 ), midpoint of housing (T 2 ), front bearing D (T 3 ), midpoint of shaft (T 4 ), respectively. The specific locations will be illustrated latterly in Section 4.2. TNM contains four types of components, thermal resistance, speed-related thermal resistance, thermal capacitance, and heat source. The authors classify the thermal behavior of spindle into two modes, operating mode and natural cooling mode. At operating mode, the spindle is operated at a constant rotation speed and the initial temperature is equal to the ambient temperature, illustrated in Figure 2a. At natural cooling mode, the spindle is stopped and releases the heat that is stored during operating mode, as illustrated in Figure 2b. The state-space representation of the TNM at operating mode is expressed as .
Sensors 2018, 18, x FOR PEER REVIEW 7 of 18 The authors define the temperature difference (θ) with respect to the ambient temperature.
Substituting it into Equation (17) gives the heat balance equations Consider the physical insight and geometry of the spindle; we construct the lumped-parameter TNM of the spindle with four nodes, illustrated in Figure 2, which are rear bearing A(T1), midpoint of housing(T2), front bearing D(T3), midpoint of shaft(T4), respectively. The specific locations will be illustrated latterly in Section 4.2. TNM contains four types of components, thermal resistance, speed-related thermal resistance, thermal capacitance, and heat source. The authors classify the thermal behavior of spindle into two modes, operating mode and natural cooling mode. At operating mode, the spindle is operated at a constant rotation speed and the initial temperature is equal to the ambient temperature, illustrated in Figure 2a. At natural cooling mode, the spindle is stopped and releases the heat that is stored during operating mode, as illustrated in Figure 2b. The state-space representation of the TNM at operating mode is expressed as According to the assumptions, the convective resistance coefficient (ravj, r24) and the heat generation coefficient (qf1-qf4) are unknown parameters needed to be solved.
where ω0 is 5000 rpm, thus the unit of the coefficients can be meaningful. At steady state, the thermal According to the assumptions, the convective resistance coefficient (r avj , r 24 ) and the heat generation coefficient (q f1 -q f4 ) are unknown parameters needed to be solved.
where ω 0 is 5000 rpm, thus the unit of the coefficients can be meaningful. At steady state, the thermal capacitances can be ignored due to the temperature derivative terms decay at steady state. Therefore, rewriting the heat balance equation as: Equation (22) is a Single-Input-Multiple-Output (SIMO) system whose input is the spindle rotation speed (ω) while the spindle temperatures (θ) are the outputs. As the speed-related parameters varying with different spindle rotation speeds, the thermal response of TNM would change as well; therefore, these nonlinear behaviors make the TNM more complicated. Likewise, the state-space representation of the TNM at natural cooling mode is expressed as: The specific system parameter matrices are shown in Appendix A.

Self-Made Bluetooth Temperature Sensor Module
The performance (accuracy, resolution, stability . . . ) of the sensor is one of the key points directly affecting the quality of the system identification results. Due to the high stability and accuracy of RTD, the authors develop a BTSM that provides multi-channel temperature measurement with high precision. Figure 3a shows the circuit design of BTSM, it is implemented with a Bluno Beetle [30] as the microcontroller, a 24-bit AD7794 [31] as the analog-to-digital converter, a reference resistor with low temperature coefficient of resistance, 4-wire Kelvin measurement, and the Printed Circuit Board (PCB) Layout ( Figure 3b). After several versions of evolution, the BTSM (Figure 4a) has achieved with high accuracy of ±(0.1 + 0.0029|ϑ|) • C, the resolution of 0.00489 • C, and the module size is only Ø40 mm × 27 mm. The specification is listed in Table 1. As the front-end temperature acquisition element, three types of temperature probe (screw type, magnetic mount type, and probe type) are fabricated by platinum RTD (PT1000) and sealed with thermal conductive Gel (Figure 4a). Magnet type probe is made for surface temperature measurement, and the other types are made for the shaft, housing, bearing temperature measurement. The fabrication process of the temperature probe is shown in supplementary materials Figure S1. It should be mentioned here that, to improve the precision of the BTSM, the Kalman filter is used to filter the noise from the measured temperature signal.
types of temperature probe (screw type, magnetic mount type, and probe type) are fabricated by platinum RTD (PT1000) and sealed with thermal conductive Gel (Figure 4a). Magnet type probe is made for surface temperature measurement, and the other types are made for the shaft, housing, bearing temperature measurement. The fabrication process of the temperature probe is shown in supplementary materials Figure S1. It should be mentioned here that, to improve the precision of the BTSM, the Kalman filter is used to filter the noise from the measured temperature signal.

Experiment Setup
For validating the commercial application of the present TNM and BTSM, the authors purchased a customized, but without loss of generality, spindle warm-up system. Illustrated in Figure 4b, it contains a spindle driven by a 7.5 kW motor through belt, a controller accompanying with a PC being used to conduct the programmable control of motor, a commercial thermocouple device being used to verify the feasibility of the self-made BTSM, etc. The embedded locations of the temperature probes are designed based on the thermal characteristic of machine tool spindle. Specifically illustrated in Figure 5, the temperature sensors are strategically located at 12 nodes near the front and rear bearings, spindle housing, shaft, spindle surface, and ambient temperature, respectively.

Experiment Setup
For validating the commercial application of the present TNM and BTSM, the authors purchased a customized, but without loss of generality, spindle warm-up system. Illustrated in Figure 4b, it contains a spindle driven by a 7.5 kW motor through belt, a controller accompanying with a PC being used to conduct the programmable control of motor, a commercial thermocouple device being used to verify the feasibility of the self-made BTSM, etc. The embedded locations of the temperature probes are designed based on the thermal characteristic of machine tool spindle. Specifically illustrated in Figure 5, the temperature sensors are strategically located at 12 nodes near the front and rear bearings, spindle housing, shaft, spindle surface, and ambient temperature, respectively. contains a spindle driven by a 7.5 kW motor through belt, a controller accompanying with a PC being used to conduct the programmable control of motor, a commercial thermocouple device being used to verify the feasibility of the self-made BTSM, etc. The embedded locations of the temperature probes are designed based on the thermal characteristic of machine tool spindle. Specifically illustrated in Figure 5, the temperature sensors are strategically located at 12 nodes near the front and rear bearings, spindle housing, shaft, spindle surface, and ambient temperature, respectively.

Results
The transient temperature distribution of spindle is measured by BTSM at various rotation speeds. Figure 6a shows a typical temperature rise curves at the rotation speed of 6021 rpm. Note that, at operating mode is similar to the charging mode of R-C circuit, more specifically, the zero-state response. Likewise, the temperature fall behavior at the natural cooling mode is similar to the discharging mode of R-C circuit, more specifically, the zero-input response. Turn off the spindle right after reaching its steady state, meanwhile insert the temperature probe into the shaft.

Results
The transient temperature distribution of spindle is measured by BTSM at various rotation speeds. Figure 6a shows a typical temperature rise curves at the rotation speed of 6021 rpm. Note that, at operating mode is similar to the charging mode of R-C circuit, more specifically, the zero-state response. Likewise, the temperature fall behavior at the natural cooling mode is similar to the discharging mode of R-C circuit, more specifically, the zero-input response. Turn off the spindle right after reaching its steady state, meanwhile insert the temperature probe into the shaft. Consequently, the complete steady-state temperature distribution of spindle can be measured. Zoom in the region of the rectangular frame line of Figure 6a, namely the transition from operating mode to natural cooling mode, Figure 6b reveals that the uneven temperature distribution accumulated at operating mode will gradually reach the thermal equilibrium in nearly 1000 s. After that, the entire spindle behaves relatively as a homogeneous temperature element and cools through free heat convection and radiation heat transfer. It implies that the thermal convective resistances are greater than the thermal conductive resistances. As a result, it further proves that the lumped element assumption is valid. Following the system identification procedure, as mentioned in Section 2.2, the estimated parameters of the TNM are listed in Tables 2 and 3, and the value of the thermal parameters are calculated at the rotational speed of 5000 rpm. Consequently, the complete steady-state temperature distribution of spindle can be measured. Zoom in the region of the rectangular frame line of Figure 6a, namely the transition from operating mode to natural cooling mode, Figure 6b reveals that the uneven temperature distribution accumulated at operating mode will gradually reach the thermal equilibrium in nearly 1000 s. After that, the entire spindle behaves relatively as a homogeneous temperature element and cools through free heat convection and radiation heat transfer. It implies that the thermal convective resistances are greater than the thermal conductive resistances. As a result, it further proves that the lumped element assumption is valid. Following the system identification procedure, as mentioned in Section 2.2, the estimated parameters of the TNM are listed in Tables 2 and 3, and the value of the thermal parameters are calculated at the rotational speed of 5000 rpm.   Table 3. Estimated thermal parameters of TNM on natural cooling mode.  To verify the assumptions of the free convective and radiative heat transfer, the measured temperature of rear bearing (T 1 ) at steady state is used to simulate h free and h rad . Due to the temperature of rear bearing is the highest one among all of the other positions, the simulation results can indicate the significant influence. In Figure 7a, the approximate computation of R rad (the purple dashed line) from Equation (15) is almost a constant value. As for free convective heat transfer coefficient, the authors consider the spindle to be a homogeneous cylinder with the diameter (D) of 0.135 m and length (L) of 0.278 m, and its other material properties are based on 25 • C. Consequently, the approximately calculated R free (green dotted line) from Equation (13) can be curve-fitted with the formula 1/A h f ree,c + h f ree,s ω −0.5 . For this reason, it is appropriate to adopt these assumptions. Figure 7b shows the predicted heat generations of rear and front bearings. One can estimate the dissipated heat through radiation (Q rad.approx. ) and free convection (Q free.approx. ) by means of the preceding approximated h free and h rad . Furthermore, Figure 7a shows the speed-dependent parameters with respect to the rotational speed. It indicates that the nonlinear behavior is more significant under lower rotational speed. For this reason, it is necessary to take the speed-dependent thermal behavior into consideration. In other words, due to the varying of the thermal parameters, the poles of the estimated TNM would slightly change with different operating speeds, which will be discussed latterly in Section 5.4. . For this reason, it is appropriate to adopt these assumptions. Figure 7b shows the predicted heat generations of rear and front bearings. One can estimate the dissipated heat through radiation (Qrad.approx.) and free convection (Qfree.approx.) by means of the preceding approximated hfree and hrad. Furthermore, Figure 7a shows the speed-dependent parameters with respect to the rotational speed. It indicates that the nonlinear behavior is more significant under lower rotational speed. For this reason, it is necessary to take the speed-dependent thermal behavior into consideration. In other words, due to the varying of the thermal parameters, the poles of the estimated TNM would slightly change with different operating speeds, which will be discussed latterly in Section 5.4.

Steady State Self-Validation
After applying the nonlinear least square estimation to find the minimum 2-norm solution of the overdetermined system in state space representation, the global solutions at steady state can be found, which represent the sets of the optimized thermal resistance value. First, the model is self-validated at steady state through a mutual comparison of the predicted and measured

Steady State Self-Validation
After applying the nonlinear least square estimation to find the minimum 2-norm solution of the overdetermined system in state space representation, the global solutions at steady state can be found, which represent the sets of the optimized thermal resistance value. First, the model is self-validated at steady state through a mutual comparison of the predicted and measured temperature at five different rotation speeds, say 4006, 5013, 6021, 7023, and 8028 rpm, as shown in Figure 8a. The predicted and measured temperatures agree very well with each other, and the maximum mean average percentage error (MAPE) is only 5.54%. Moreover, it should be emphasized that the heat generation and steady state temperature are closely related, which is confirmed from the similarity of Figures 7b and 8b. (a) (b) Figure 7. (a) Simulated heat generation varying with rotational speeds; (b) Estimated convective resistance varying with rotational speeds.

Steady State Self-Validation
After applying the nonlinear least square estimation to find the minimum 2-norm solution of the overdetermined system in state space representation, the global solutions at steady state can be found, which represent the sets of the optimized thermal resistance value. First, the model is self-validated at steady state through a mutual comparison of the predicted and measured temperature at five different rotation speeds, say 4006, 5013, 6021, 7023, and 8028 rpm, as shown in Figure 8a. The predicted and measured temperatures agree very well with each other, and the maximum mean average percentage error (MAPE) is only 5.54%. Moreover, it should be emphasized that the heat generation and steady state temperature are closely related, which is confirmed from the similarity of Figures 7b and 8b. (a) (b) Figure 8. (a) Self-validation between predicted and measured steady state temperature, T1 to T4 are the positions of the rear bearing A, midpoint of inner housing, front bearing D, and midpoint of shaft, respectively; (b) Steady-state temperature validation.

Transient State Self-Validation
The grey-box model identification technique is implemented to estimate the remaining parameters of the TNM at both operating and natural cooling modes. The transient self-validations are obtained by comparing the transient temperatures rise estimated by TNM with the measured ones. Five transient states are validated; namely, the spindle runs from stationary to five different steady speeds, say 4006, 5013, 6021, 7023, and 8028 rpm. Figures 9 and 10 show promising validation with the minimum normalized mean square errors (NMSE) of 94.2% at operating mode and 95.8% at natural cooling mode. When the spindle operates at operating mode, the spindle shaft rotates at a constant speed; hence, the temperature probe is not able to contact the shaft. For this reason, the experimental temperature of spindle shaft (T 4 ) cannot be obtained, and the T 4 -est represents the estimated shaft temperature, as shown in Figure 9d.

External Validation
To verify the performance and robustness of the TNM, the external validation is necessary. The experiment under stepwise rotational speed (3001, 5018, 7028 rpm) is utilized especially. The results ( Figure 11) reveal that the predicted temperature is in satisfactory agreement with the experimental results. As a result, this TNM methodology is appropriate to implement to the machine tool spindle as a simplified thermal model. The grey-box model identification technique is implemented to estimate the remaining parameters of the TNM at both operating and natural cooling modes. The transient self-validations are obtained by comparing the transient temperatures rise estimated by TNM with the measured ones. Five transient states are validated; namely, the spindle runs from stationary to five different steady speeds, say 4006, 5013, 6021, 7023, and 8028 rpm. Figures 9 and 10 show promising validation with the minimum normalized mean square errors (NMSE) of 94.2% at operating mode and 95.8% at natural cooling mode. When the spindle operates at operating mode, the spindle shaft rotates at a constant speed; hence, the temperature probe is not able to contact the shaft. For this reason, the experimental temperature of spindle shaft (T4) cannot be obtained, and the T4-est represents the estimated shaft temperature, as shown in Figure 9d.

External Validation
To verify the performance and robustness of the TNM, the external validation is necessary. The experiment under stepwise rotational speed (3001, 5018, 7028 rpm) is utilized especially. The results ( Figure 11) reveal that the predicted temperature is in satisfactory agreement with the experimental results. As a result, this TNM methodology is appropriate to implement to the machine tool spindle as a simplified thermal model.

Model Order Reduction
To implement the system parameter estimation into a micro control unit for real-time onboard application on edge computing, the order of TNM model should be further simplified to reduce the load of computation and data storage. MOR technique is an excellent approach for reducing the complexity of the estimated TNM. Balanced realization method is proposed to reduce the 4th-order TNM into 1st-order truncated model, which is a simplified realization with equal controllability and observability. Figure 12a shows the Hankel singular values of the estimated TNM; it points out that there exists a leading state with the larger value than other states. For this reason, we remove the less three states (2, 3, 4) and apply the first dominant state to reconstruct the 1st-order truncated TNM. The pole location of the estimated TNM and truncated model are shown in Figure 12c. In addition, the Bode diagram illustrates the relationship between the input heat generation (Q f1 , Q f2 ) and the output temperature T 1 to T 4 at rotational speed of 6000 rpm. It indicates that the truncated model has equivalent performance to the original estimated TNM under low frequency (<10 −3 ) condition, as shown in Figure 13.

External Validation
To verify the performance and robustness of the TNM, the external validation is necessary. The experiment under stepwise rotational speed (3001, 5018, 7028 rpm) is utilized especially. The results ( Figure 10) reveal that the predicted temperature is in satisfactory agreement with the experimental results. As a result, this TNM methodology is appropriate to implement to the machine tool spindle as a simplified thermal model.

Model Order Reduction
To implement the system parameter estimation into a micro control unit for real-time onboard application on edge computing, the order of TNM model should be further simplified to reduce the load of computation and data storage. MOR technique is an excellent approach for reducing the complexity of the estimated TNM. Balanced realization method is proposed to reduce the 4th-order TNM into 1st-order truncated model, which is a simplified realization with equal controllability and observability. Figure 11a shows the Hankel singular values of the estimated TNM; it points out that there exists a leading state with the larger value than other states. For this reason, we remove the less three states (2, 3, 4) and apply the first dominant state to reconstruct the 1st-order truncated TNM. The pole location of the estimated TNM and truncated model are shown in Figure 11c. In addition, the Bode diagram illustrates the relationship between the input heat generation (Qf1, Qf2) and the output temperature T1 to T4 at rotational speed of 6000 rpm. It indicates that the truncated model has equivalent performance to the original estimated TNM under low frequency (<10 -3 ) condition, as shown in Figure 12.

Short Circuit Time Constant
The Short Circuit Time Constant (SCTC) method provides an approximation of the poles of the estimated TNM. Each pole can be calculated according to each thermal capacitance multiplied by the nearby equivalent parallel thermal resistance while the other capacitances are short circuited, refer to Figure 11b. As the aspect of system degradation, we can simply recognize the time-varying trend by examining the variation of the time constant instantly by applying SCTC method.

Conclusions
In this article, we investigate the possibility of simplifying the spindle thermal behaviors by

Short Circuit Time Constant
The Short Circuit Time Constant (SCTC) method provides an approximation of the poles of the estimated TNM. Each pole can be calculated according to each thermal capacitance multiplied by the nearby equivalent parallel thermal resistance while the other capacitances are short circuited, refer to Figure 12b. As the aspect of system degradation, we can simply recognize the time-varying trend by examining the variation of the time constant instantly by applying SCTC method.

Conclusions
In this article, we investigate the possibility of simplifying the spindle thermal behaviors by estimating the lumped-parameter TNM with system identification technique. We classify the spindle heat transfer condition into two modes, the operating mode and the natural cooling mode. The predicted thermal response of TNM agrees well with the experimental temperature results, verified by self-validation and external validation. In the aspect of the real-time onboard application, the reduced-order model by balanced realization has successfully lowered the computational complexity. In practice, the spindle warm-up procedure is essential to all of the machine tool spindles, especially for the new spindle. By operating at step-wise speeds, the lubricant within the spindle bearings can be evenly spread. The modeling methodology in this research can be appropriately integrated with the warm-up procedure and establish the estimated TNM for further operation. Figure 14 summarizes that the TNM methodology has the ability to predict the spindle thermal behavior; and the BTSM provides real-time monitoring the spindle temperature distribution. In that case, the vision of edge computing within a minimize temperature sensor module can be realized by integrated with the estimated TNM performed as the thermo-feature identification system (TID-sys), which can be customized for each machine tool spindle. Figure 15 demonstrates the real-time user interface of TID-sys on smartphone and PC. On the right-hand side, the TID-sys diagnosis the state (normal or alarm) of the machine tool spindle displayed through the terminal App [32] on smartphone. On the left-hand side, the measured temperature (black dot) of rear bearing (T 1 ), midpoint of housing (T 2 ), and front bearing (T 3 ) attempt to grow along with the curve predicted by the estimated TNM. The TID-sys ensure that the machine tool spindle runs under normal condition. Once the significant temperature variation or abnormal temperature rise is detected, the operator can slow down the spindle or even stop for early protection. The full demonstration video is shown in supplementary materials Video S1.

Supplementary Materials:
The following are available online at www.mdpi.com/link, Figure S1: Fabrication process of the temperature probe, Figure S2: Experimental setup; and Video S1: Demonstration of the TID-sys.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. System Parameter Matrices
The input, output and state vector of TNM and linear in state space representation: The heat balance equations of TNM: Figure 15. Demonstration of the thermo-feature identification system.

Supplementary Materials:
The following are available online at www.mdpi.com/1424-8220/18/2/656/s1, Figure S1: Fabrication process of the temperature probe, Figure S2: Experimental setup; and Video S1: Demonstration of the TID-sys. Author Contributions: This paper was performed in collaboration among the authors. Yuan-Chieh Lo was responsible for the parameter estimation of the thermal network model, experiment works, and writing the article. Yuh-Chung Hu was the research project instructor and responsible for advising the research works and article writing. Pei-Zen Chang is Lo's graduate advisor.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. System Parameter Matrices
The input, output and state vector of TNM and linear in state space representation: The heat balance equations of TNM: Sensors 2018, 18, 656 18 of 20 The TNM system matrices of state-space representation at operating mode: (A3) The TNM system matrices of state-space representation at steady state: q f 1 · · · q f 1 0 · · · 0 q f 2 · · · q f 2 0 · · · 0      ; E std =      q f 3 · · · q f 3 0 · · · 0 q f 4 · · · q f 4 0 · · · 0      ; (A4) The TNM system matrices of state-space representation at natural cooling mode: Normalized mean square error: x re f − x re f 2 (A7)