Co-simulation-based digital twin for thermal characteristics of motorized spindle

To improve the accuracy of thermal characteristics analysis of motorized spindle, an online correction model of thermal boundary conditions is proposed based on BP neural network (BPNN), the experimental data and simulation results are used to construct the BPNN model to correct the thermal boundary conditions of motorized spindle. Based on the co-simulation of Ansys, Matlab, and LabVIEW, a digital twin system for thermal characteristics is built to precisely predict the temperature field and thermal deformation of a motorized spindle under varied operating conditions. The experimental results show that the prediction accuracy of temperature field is greater than 98%, and the prediction accuracy of thermal deformation is greater than 96%, which effectively improves the simulation accuracy and robustness of thermal characteristics, and provides the foundation for the error compensation and thermal optimization design.


Introduction
Intelligent manufacturing closely combines industrialization and informatization, in a manner that significantly enhances production efficiency and quality [1,2]. As one of the primary technologies of intelligent manufacturing, digital twin simulates the behavior of real-world physical things using actual production data. As shown in Fig. 1, its interaction [3] with the Internet of Things (IoT) and cyber-physical system (CPS) is as follows: the physical model of the real world collects data in real time through IoT acquisition system, which is then mapped to the virtual world via CPS to reflect the whole life cycle process of the corresponding physical entity. Liu et al. [4] combined digital twin with finite element method (FEM) and proposed a digital twin-driven approach for traceability and dynamic control of processing quality, which enabled interactive integration of physical and virtual workshops and was applied to the processing of diesel engine connecting rod. Combining computer simulation, algorithm analysis, and experience assessment, Lu et al. [5] constructed a multi-dimensional digital twin model devoted to the product life cycle of the constant velocity joint of a car. As one of the key technologies of intelligent manufacturing, digital twin has been applied in many fields such as product design and manufacturing. However, the application of digital twin in thermal characteristics of motorized spindle is relatively scarce.
As the core component of high-speed computer numerical control machine tools, the thermal deformation of the motorized spindle generated by high-speed operation has a significant impact on the precision of machine tools. Consequently, it is crucial to precisely analyze the thermal characteristics of motorized spindle. In recent years, a significant deal of study has been conducted on the thermal characteristics and thermal design of motorized spindles [6], which has substantially improved their performance. Fan et al. [7] used partial least squares (PLS) to predict the thermal deformation of motorized spindle under varied operating conditions, and analyzed the correlation between multiple temperature variables and three-dimensional thermal deformation, which increased the accuracy of the model. Uhlmanna et al. [8] proposed a 3D finite element (FE) prediction model for high precision thermal characteristics of motorized spindle under complicated boundary conditions such as heat source, convection and contact. Zhang et al. [9] used the least square to optimize the heat transfer coefficient of the motorized spindle and took the optimized heat transfer coefficient as the boundary of the FE model to increase the accuracy of the temperature rise prediction model. Zhang et al. [10] also proposed a method for optimizing the heat transfer coefficient based on the biogeography optimization algorithm and proposed a thermal deformation prediction model of an intelligent motorized spindle. Compared to the traditional thermal deformation prediction model, the proposed model is more precise.
With the increasing speed of motorized spindle, the analysis and improvement of its thermal characteristics have become a hot spot. At present, the thermal analysis method of motorized spindle is mainly FE thermal simulation. By establishing FE model and applying thermal boundary conditions, the temperature field and thermal deformation of motorized spindle are obtained by FE analysis software. The key is identifying thermal boundary conditions precisely. Currently, the most prevalent identification approaches consist of theoretical calculation method [11], experimental measurement method [12], intelligent optimization method, and inverse method [13], etc. The existing identification methods can obtain the thermal boundary conditions under specific conditions, providing a foundation for thermal characteristics analysis and thermal optimization design. However, due to the complexity of thermal boundary conditions and experimental limitations, there is a certain deviation between the identification accuracy and the actual value [6]; Furthermore, due to the shape, location and material of the thermal boundary, as well as other factors, it is difficult to determine the exact value [14]. Therefore, it is necessary to carry out in-depth research on accurate identification methods and correction techniques of thermal boundary conditions.
In order to increase the accuracy of thermal analysis of motorized spindle, a modified thermal boundary condition model based on BP neural network (BPNN) is proposed to identify more accurate thermal boundary conditions under varied operating conditions. The digital twin system of motorized spindle thermal characteristics based on co-simulation is intended to increase the accuracy of temperature field and thermal deformation predictions for motorized spindles.

BPNN-based correction mechanism for heat generation of bearing
During high-speed rotation, the friction between the bearing rollers and the inner and outer rings is the primary cause of bearing heat generation [15,16]. The higher the bearing speed is, the more intense the friction is and the greater the heat generation is. In finite element simulation, the bearing roller is usually regarded as the heat source of the bearing, and its heat generation can be calculated according to Palmgren formula and empirical correction formula [17] as, where, H f is the initial heat generation of the bearing in W∕m 3 ; n is the spindle speed in rpm; M is the bearing friction torque in N ⋅ mm , V m is the bearing volume in m 3 ; M 1 is the viscous friction torque of lubricating oil; M 2 is the load torque. f 0 is a constant related to bearing type and lubrication mode; v 0 (Δt) is the kinematic viscosity of lubricating oil changing with temperature in mm 2 ∕s ; D m is the average diameter of the bearing in mm; f 1 is a coefficient related to Fig. 1 Relationship between digital twin, CPS, and IoT bearing type and load; P 1 (Δt) is the bearing load changing with temperature in Newton. It can be seen from Eq. (1) that for an assembled spindle, the bearing heat generation is proportional to the speed n of the motorized spindle, the bearing load P 1 (Δt) and the kinematic viscosity of the lubricating oil v 0 (Δt) . Since it is difficult to quantify the bearing load precisely and its magnitude fluctuates with the bearing temperature, it is difficult to accurately calculate the heat generation of the bearing. In order to increase the identification accuracy of bearing heat generation, Reference [15] proposed the correction model of bearing heat generation as, where, H ′ f is the corrected bearing heat in W∕m 3 ; t m is the measured temperature of the thermal key points in °C; t f is the finite element analysis temperature of the thermal key points in °C; t ∞ is ambient temperature in °C.
Using Eq. (2), the identification accuracy of heat generation can be increased greatly. To further increase the identification accuracy of heat generation, a BPNN-based correction method using the instantaneous temperature rise ratio is proposed in order to further correct the heat generation calculated by Eq. (2). If the measured temperature is greater than the simulated temperature, heat transfer theory dictates that the identified heat generation of the bearing is smaller than the actual heat generation of the bearing. According to the relationship between the heat and temperature rise Q = c·m·Δt, the heat flux is proportional to the temperature rise Δt. Therefore, the instantaneous temperature rise ratio is introduced as a correction factor in this study. Using the ratio of measured and simulated instantaneous temperature increase as a correction factor for bearing heat generation, the correction function can be established as, where, C H is the correction factor; Δt m andΔt f are the measured and simulated instantaneous temperature rise, respectively; H ′′ f is the bearing heat generation corrected by correction factor in W∕m 3 .
Due to the nonlinearity of temperature rise, BP neural network is used to predict the heat generation in this study. According to Eqs. (1,2,3), the spindle speed, ambient temperature, measured, and simulated temperatures at thermal key points of bearing, operating duration, and initial heat generation can be regarded as the input of neural network; The corrected bearing heat generation can be regarded as the output of the neural network. According to the analysis abovementioned, the neural network prediction model of front and rear bearing heat generation is composed of a 6-neuron input layer, a 10-neuron hidden layer using the purelin transfer function, and a single output neuron with a logsig transfer function respectively.
In order to avoid too small adjustment range of neural network weight and threshold, the elastic gradient descent algorithm (trainrp) is used to train BP network. Table 1 shows some samples of the heat generation correction of the front bearing of the motorized spindle.
The front bearing's heat generation varies nonlinearly, as shown in Table 1 This is mostly due to the thermal expansion of the ball during the bearing's initial operation, which increases the bearing load. As the inner and outer rings of the bearing expand, the bearing load progressively reduces, resulting in the decrease of the heat generation. With the expansion of the main shaft, bearing sleeve and other components, the overall expansion of the bearing also changes continuously. Finally, when the heat balance is reached, the heat generation of the bearing is constant.

Correction mechanism of thermal contact resistance using neural network
Because the bonding surface is microscopically uneven, i.e., the two surfaces are not totally in touch, and the gap between the non-contact interface is filled with medium to form thermal contact resistance [17][18][19]. In ANSYS thermal simulation, it is necessary to accurately define the thermal contact resistance to ensure the accuracy of simulation. Therefore, it is significant to identify the thermal contact resistance through calculation and combined with the correction formula. The initial thermal contact resistance of the bonding surface can be calculated by the following semiempirical formula as, where, R is the initial thermal contact resistance in m 2 /K·W; r is the root mean square roughness in Ra; k is the tuned average thermal conductivity in W/m·K; P is the contact pressure in Mpa; E is the elastic modulus in Gpa. It can be seen from Eq. (4) that the thermal contact resistance of the bonding surface is proportional to the contact pressure. For the assembled motorized spindle, it is difficult to quantify the bonding surface's contact pressure. Therefore, the correction formula of thermal contact resistance is proposed in Reference [17] as, where, R' is the corrected thermal contact resistance in m 2 /K·W; Δt is the temperature change of thermal key points in °C; α is the coefficient of thermal expansion; A is the area of contact surface in m 2 .
If the measured temperature is greater than the simulated temperature, heat transfer theory dictates that the identified thermal contact resistance is greater than the actual value of the bonding surface, and the thermal contact resistance is directly proportional to the interface temperature difference. Therefore, this study introduces a correction factor, i.e., the ratio of the measured temperature rise to the simulated temperature rise as a correction factor of the thermal contact resistance, and establishes a correction function as, where, C C is the correction factor; R" is the thermal contact resistance corrected by correction factor in m 2 /K·W.
In Eq. (6), the instantaneous temperature rise is related to the heat generation which is related to the spindle speed. Therefore, the spindle speed, the measured and simulated temperatures of thermal key points, the ambient temperature, operating duration, and the initial thermal contact resistance can be regarded as the input of the neural network; The corrected thermal contact resistance can be regarded as the output of the neural network. The thermal contact resistance neural network model is composed of a 6-neuron input layer, a 12-neuron hidden layer, and 1-neuron output layer. The transfer function and network training algorithm are consistent with the neural network for predicting bearing heat generation.

Sample collection and training
The accuracy of a neural network's predictions relies on the precision of its samples. In Sections 2.1 and 2.2, input and output layer samples are collected by a combination of experiment, simulation, and calculation. The particular steps are as follows: The first step is to conduct 20 thermal characteristic tests of motorized spindle rotating at various speeds for 30 min, and store the spindle speed, ambient temperature, and the measured temperature of thermal key points at different running time in the database as input layer samples.
In the second step, Eqs. (1) and (4) are applied to calculate the initial heat generation and thermal contact resistance of the motorized spindle as the input layer samples of the neural network.
The third step is to calculate other boundary conditions of thermal characteristic analysis. The motor's heat generation is calculated using the electromagnetic loss formula, while the convective heat transfer coefficient is calculated using the Baz formula.
In the fourth step, the thermal characteristics of the motorized spindle are simulated under the same working conditions as the tests, and the simulated temperatures of thermal key points at the same running time were stored in the database as the input layer samples.
In the fifth step, the measured and simulated temperature of thermal key points are fitted a function of time to calculate the correction parameters of C H and C C . Then, the bearing heat generation and thermal contact resistance with correction factors can be calculated by Eqs. (3) and (6) as the training samples of output layers of the neural network. The first corrected heat generation H' and thermal contact resistance R' are obtained by online simulation.
To ensure sufficient convergence, the training parameters are set as 5000 epochs, the learning speed is set as 0.05, and the error is set as 10 −6 , which is calculated by the mean square error function as, where, k is the number of neurons in the output layer; p is the number of training samples; ŷ kj is the expected output of the neural network; y kj is the output of neural network.
To improve the training accuracy, the training samples are standardized to the range of (0, 1) through Eq. (8), where, S is the standardized training sample; X is the initial training samples; n is the maximum order of magnitude of training samples, n = 6. Using Eq. (8) to normalize the training samples of neural network can avoid the issue of computer rounding error and obtain more accurate training results. For the analysis effect of neural network algorithm regression, the R squared index as shown in Eq. (9) is utilized; the closer it is to 1, the better.
where Y is the experimental value; O is the prediction value. The samples are divided into training set and test set according to the ratio of 8:2, and the thermal boundary condition correction model is obtained by training. Figure 2 is the prediction results of partial thermal boundary conditions. It can be seen from Fig. 2 that all the model evaluation indexes are greater than 0.99. It demonstrates that the prediction effect of motorized spindle's thermal boundary conditions based on BP neural network is good, which can effectively increase the accuracy of finite element analysis. Fig. 2 a Prediction of thermal contact resistance between front bearing and bearing sleeve. b Prediction of thermal contact resistance between rear bearing and bearing sleeve. c Prediction of heat generation of front bearing. d Prediction of heat generation of rear bearing

Digital twin for thermal characteristics of motorized spindle
Digital twin is to establish the digital image of the actual physical system in the virtual system, expand new capabilities for the planning and design, state monitoring and motion control of the physical entity, and map the corresponding entity reality through virtual space. As shown in Fig. 3, the digital twin principle of thermal characteristics of motorized spindle is to use the IoT data acquisition system to collect the temperature of thermal key points inside the motorized spindle. After data preprocessing, the thermal boundary conditions are predicted by BP neural network and applied to the virtual model of the motorized spindle. ANSYS parametric design language (APDL) is used to analyze the thermal characteristics of the motorized spindle. The analysis results are mapped to the virtual space to realize the digital twin of the thermal characteristics of the motorized spindle. Figure 4 shows the architecture of thermal characteristic digital twin system [20]. The digital twin system includes three parts: physical space, twin space, and virtual space. The twin space is connected with physical space and virtual space through IoT data acquisition system. Combined with the parameters such as ambient temperature and motorized spindle speed input from the user system interface, the real-time thermal boundary conditions of motorized spindle are predicted through BP neural network and output to APDL command stream file for thermal characteristic simulation. The simulation results are then mapped to virtual space to realize digital twin of thermal characteristics of motorized spindle.

Physical space development
The physical space is composed of motorized spindle, sensors, IoT data acquisition system and computer, as shown in Fig. 5. The IoT data acquisition system collects the temperature and displacement data of thermal key points of the motorized spindle measured by the sensors. IoT acquisition system realizes data interaction with host computer database through IP address and port number. Figure 6 shows the embedded position of the temperature sensors in the motorized spindle. Five temperature sensors are used to measure the temperature of thermal key points during motorized spindle running [21], in which T1 and T2 are installed at the front bearing outer ring and front bearing sleeve to measure the temperature of the front bearing outer ring and bearing sleeve; T3 and T4 are installed at the rear bearing outer ring and rear bearing sleeve to measure the temperature of the rear bearing outer ring and bearing sleeve; T5 is pasted on the motorized spindle housing to measure the temperature of the motorized spindle shell. The information acquired by each temperature sensor is utilized to correct the thermal boundary conditions of the motorized spindle in real time.

Twin space development
Twin space mainly includes data access module and graphic modeling module [22]. The data access module is mainly the data interaction between MATLAB, LabVIEW, and IoT data acquisition system. The communication software based on TCP / IP protocol is developed by using LabVIEW to import the experimental data into the database in real time. Then, combined with the ambient temperature and motorized spindle speed input by the user in the graphical interface as the input layer of the neural network model, the  Fig. 7.
The client and the server realize information transmission through sockets. After establishing the connection, the server analyzes the request message according to the protocol specification after receiving the request message, and then returns the response message to the client in order to acquire real-time data.
The graphical modeling module is the data interaction between Matlab GUI and APDL command stream file. Due to the fact that using ANSYS for thermal characteristics simulation analysis requires a certain amount of calculation time, in order to avoid mutual error, the predicted thermal boundary conditions are transmitted to the APDL command stream file every 60 s for Ansys batch startup for thermal-structural coupling simulation. The APDL post-processing content includes the export of motorized spindle temperature field, thermal deformation field,   simulation temperature of thermal key points, and correction data. The post-processing results are saved to the database, where the simulation temperature of thermal key points and the operating duration of motorized spindle are called to the neural network model as the input layer. The specific data interaction process of twin space is shown in Fig. 8.

Virtual space development
As shown in Fig. 9, a digital twin system based on Matlab GUI, APDL, and LabVIEW is developed to study the thermal characteristics of motorized spindle. The application can realize the starting, stopping, and speed change of motorized spindle, real-time data interaction, result visualization, monitoring correction accuracy, etc. In the system interface, the user only needs to input the ambient temperature and motorized spindle speed. In the background, the input parameters are transmitted to the twin space every 60 s through the I/O flow to complete the transient thermal-structure coupling simulation with a step length of 4 s and a step end time of 60 s, and then the simulated temperature field and thermal deformation of the motorized spindle are displayed in real time via UIAxes component in Matlab GUI. The real-time temperature of the five key points, the thermal deformation of the core spindle along X-, Y-, and Z-directions, and the monitoring accuracy are displayed in real time through the UITable component. Finally, the historical data is stored in the database to provide a basis for error compensation and thermal optimization design.

Experimental validation
For a type of high-speed motorized spindle, the physical space is arranged as shown in Figs. 5 and 6. To verify the accuracy and the robustness of the model built in this paper, two verification experiments named experiment I and experiment II have been carried out on the motorized spindle.
In experiment I, the thermal characteristics of the motorized spindle operating at the speed given in Table 2 Fig. 9 User interface of digital twin system for 30 min are tested by using the thermal characteristics digital twin system. Three eddy current displacement sensors are used to measure the thermal deformation of core spindle along the X-, Y-, and Z-directions, and the laboratory temperature was 25.34 °C. Figure 10 shows the digital twin diagram of thermal characteristics of motorized spindle when it rotates at 9000r/min for 30 min. The temperature field of the motorized spindle is shown in Fig. 10a. The areas with high temperature of the motorized spindle are concentrated at the front and rear bearings, stator, and rotor, of which the rotor temperature is the highest which is 57.017 °C, followed by the front and rear bearings and stator. The main reason is that the bearings and stator are assembled inside the spindle and the heat dissipation conditions are poor. However, a portion of the heat from the front and rear bearings and stator is transferred to the cooling system and taken away by the coolant, whereas the rotor can only dissipate heat through internal air convection and radiation. Figure 10a also shows that the temperature field of the motorized spindle is ununiform, which will lead to inconsistent thermal deformation of the components and thus affect the thermal boundary conditions. Figure 10b shows the thermal deformation of the motorized spindle, and the maximum thermal deformation is 8.47 μm which occurs at the front end of the core spindle.
The comparison of thermal properties between the digital twin and real measurement data at 7800 rpm and 11,400 rpm is shown in Fig. 11. It demonstrates that the digital twin of thermal characteristic based on neural network correction of thermal boundary can better map the thermal characteristic of the motorized spindle, can realize the whole process mapping of thermal characteristic of the motorized spindle, and has good following performance. Table 2 shows the digital twin accuracy of thermal characteristics of motorized spindle under various test speeds. The minimum prediction accuracy of five key temperature measurement points is 98.5%, and the average accuracy is 98.62%. The minimum prediction accuracy of thermal deformation is 95.4%, and the average accuracy is 96.06%. It proves that the digital twin for thermal characteristics of motorized spindle based on BP neural network can effectively increase the accuracy of thermal characteristic analysis of the prediction model.
Based on the comprehensive analysis of the digital twin results of the thermal characteristics of the motorized spindle, when the temperature rise rate of the motorized spindle changes rapidly, the following performance needs to be improved: As shown in Fig. 11, the temperature accuracy of the key temperature measurement point T2 near 200 s is lower than at other times due to the fact that the neural network training samples are collected at evenly distributed time intervals, while the temperature rise rate changes rapidly when it tends to thermal equilibrium. The prediction accuracy can be increased by collecting more samples throughout this period.
Experiment II is to verify whether the motorized spindle can still have good monitoring accuracy and robustness when considering temperature rise and fall. During operation, the motorized spindle may be working in different status, including the operating status, the halt status, and the power-off status. The existence of the heat generated by the internal heat source and the status of the cooling system under each spindle status is shown in Table 3. In experiment II, the spindle rotates at 9600 rpm for 20 min at first and subsequently remains in the halt status for 20 min; then, the spindle rotates at 9600 rpm for 20 min, and finally remains in the power-off status for 20 min. Figure 12 shows the comparison of thermal properties between the digital twin and real measurement data. Experimental results show that the prediction accuracy of five key temperature measurement points is 98.47%, and the prediction accuracy of thermal deformation is 96.13%, which proved that the method has good accuracy, robustness and stability under different spindle states.

Conclusion
Accurate thermal characteristic analysis of motorized spindle is the premise of thermal optimization design and error compensation. The precision of the thermal characteristic analysis of a motorized spindle is hampered due to the intricacy of the boundary conditions. To improve the accuracy  of thermal characteristic analysis, the correction model of thermal contact resistance of joint surface and bearing heat generation is constructed based on BP neural network, and the digital twin system for thermal characteristics of motorized spindle is developed based on ANSYS, MATLAB, and LabVIEW to realize the digital twin for thermal characteristics of motorized spindle. According to the test results, the following conclusions are drawn. The training accuracy of BP neural network is highly related to the samples. The samples of the input layer can be accurately obtained by sensors, while the samples of the output layer need to be obtained by combining test, simulation, and correction function. Constructing an accurate correction function can effectively improve the prediction accuracy of neural network. The collection of training samples also has a great impact on the training accuracy of the network. The sample collection volume can be increased in the area with large changes in thermal characteristics to increase the prediction accuracy. The test results show that the prediction accuracy of digital twin temperature and thermal deformation of motorized spindle based on BP neural network correcting thermal boundary conditions can reach more than 98% and 96% respectively under varied operating conditions, which can effectively increase the accuracy of finite element analysis of thermal characteristics and the robustness of the prediction model, and lay a foundation for error compensation and performance optimization of motorized spindle.
Author contribution Fan pointed out the research direction and the feasibility of the experiment. Yi operated all the experiments, programmed, and wrote the manuscript. Fan finally polished the manuscript.
Data availability All data generated or material during this study are included in this published article.
Code availability All the Matlab, APDL, and LabView code are available.

Declarations
Ethics approval The authors confirm that this work is an original work and has not been published in other places, and all the data is true.

Consent to participate
The authors consent to participate.

Consent for publication
The authors consent to publish the paper.

Competing interests
The authors declare no competing interests.