Analytical-Numerical Model for Temperature Prediction of a Serpentine Belt Drive System

The serpentine belt drive system is used in the auto industry. To avoid thermal destruction inside the belt drive and improve the thermal fatigue life of pulley materials under a variety of operating conditions, the temperature information for each load case must be determined within only a few seconds. To this end, this paper proposes an advanced thermal model to calculate the temperature distribution of a serpentine belt drive at static state operating conditions in an efficient manner. In this model, using analytical and numerical methods, a set of equations is developed according to the thermal flows and heat exchanges occurring in the system. After calculating the thermal flows of each pulley and the belt temperature, the baseline numerical simulations are modified to output the temperature distribution for each pulley. In this manner, the time-consuming numerical calculations for each pulley are performed only once and then analytically modified to provide the temperature predictions for various designed load cases, which dramatically reduces the computational time while maintaining the accuracy. Furthermore, experiments were performed to obtain the temperature data, and the results exhibited a good agreement with the corresponding calculated results. The proposed model can thus be effectively utilized for several types of belt systems and the material development of pulleys.


Introduction
The serpentine belt drive system is used in a variety of applications, especially in the auto industry ( Figure 1). In this system, power is transmitted from a driver (DR) pulley to several driven (DN) pulleys through the belt. The advantages of this system are that is economical, efficient, and compact and has a good overload protection. However, an increasing number of pulleys used in high torque/speed power transmission systems are being fabricating using fiber-reinforced polymers (FRPs), resulting in an increased risk of thermal deterioration [1]. The heat from the torque/speed loss accumulates inside these components owing to the low thermal conductivity of the FRPs, which can damage the material and reduce its life expectancy [2]. To avoid thermal destruction inside the belt drive and improve the thermal fatigue life of pulley materials, a novel thermal model is required that can promptly predict the temperature distributions of the belt drives under different load cases. Moreover, this model must be independent of the component geometries because a belt drive consists of an arbitrary belt layout with many types of pulleys.
Considering this aspect, this paper proposes a model to determine the temperature distribution of the belt drives. The proposed approach, similar to any general thermal analysis model, consists of two major substeps: determining the heat generation within the system, and analyzing the heat transfer among the components. Considering this aspect, this paper proposes a model to determine the temperature distribution of the belt drives. The proposed approach, similar to any general thermal analysis model, consists of two major substeps: determining the heat generation within the system, and analyzing the heat transfer among the components.
The heat generation of the belt drive system is strongly related to the power loss during the operation because the power loss is converted to heat and dissipated into the environment through the surfaces of the components. Several well-developed theories and experimental research are available, focusing on the power loss of belt drive systems to improve the system efficiency. The most prevalent theory for flat and V-belt drive systems groups the losses into five main contributions: belt bending, belt radial compression, belt tension, belt shear, and belt-pulley slip [3]. The first four types of power losses occur because of the belt hysteresis, whereas the last type of power loss can be attributed to the belt frictional slippage against the pulley outer radius (OD) surfaces. Manin [4] expanded this theory to a poly-V belt transmission, which involves the combination of flat and Vbelts. Silva [5,6] focused on the power losses owing to the belt hysteresis caused by the dynamic bending, stretching, shear, flank and radial compression of the belt rubber, and they further considered other losses such as the bearing loss [7] in the belt transmission. Another theory classified the power losses into three types: frictional sliding, belt longitudinal and lateral material hysteresis, and frictional losses owing to the engagement/disengagement of the belt at the entrance and exit from the pulleys [8]. However, another existing study indicated that the slip loss between the belt and pulley is the main cause of power loss in belt drive transmissions [9]. Furthermore, the power loss is influenced by several factors, one of which is the belt pressure on the pulley. Some authors reported on an analytical model [10] for the friction damping of round clamp band joints, which can be used to predict the energy dissipation caused by the friction contact between the joint components. The belt pressure on the pulley also leads to a change in the axial force of the pulley, and some experiments indicated that this force influences the power transmission efficiency, although this influence is not significant [11]. In addition, Julio and Plante [12] developed a dynamic model to predict the transmission ratio time response under different conditions of the drive pulley axial force, load torque, and engine RPM.
In some studies, the heat transfer analysis of the belt drive system was considered after determining the power loss or heat generation within the belt drive system, because the belt is an important component of the system and belt failures may occur because of the high temperature. For example, Merghache and Ghernaout focused on the heat transfer of a synchronous belt through numerical and experimental methods [13]. A general theory related to the rubber frictional heat against solids with arbitrary thermal properties was developed [14]. Mackerle concluded finiteelement analyses, including sliding frictional simulations, applied to the rubber-like materials [15]. The heat generation of the belt drive system is strongly related to the power loss during the operation because the power loss is converted to heat and dissipated into the environment through the surfaces of the components. Several well-developed theories and experimental research are available, focusing on the power loss of belt drive systems to improve the system efficiency. The most prevalent theory for flat and V-belt drive systems groups the losses into five main contributions: belt bending, belt radial compression, belt tension, belt shear, and belt-pulley slip [3]. The first four types of power losses occur because of the belt hysteresis, whereas the last type of power loss can be attributed to the belt frictional slippage against the pulley outer radius (OD) surfaces. Manin [4] expanded this theory to a poly-V belt transmission, which involves the combination of flat and V-belts. Silva [5,6] focused on the power losses owing to the belt hysteresis caused by the dynamic bending, stretching, shear, flank and radial compression of the belt rubber, and they further considered other losses such as the bearing loss [7] in the belt transmission. Another theory classified the power losses into three types: frictional sliding, belt longitudinal and lateral material hysteresis, and frictional losses owing to the engagement/disengagement of the belt at the entrance and exit from the pulleys [8]. However, another existing study indicated that the slip loss between the belt and pulley is the main cause of power loss in belt drive transmissions [9]. Furthermore, the power loss is influenced by several factors, one of which is the belt pressure on the pulley. Some authors reported on an analytical model [10] for the friction damping of round clamp band joints, which can be used to predict the energy dissipation caused by the friction contact between the joint components. The belt pressure on the pulley also leads to a change in the axial force of the pulley, and some experiments indicated that this force influences the power transmission efficiency, although this influence is not significant [11]. In addition, Julio and Plante [12] developed a dynamic model to predict the transmission ratio time response under different conditions of the drive pulley axial force, load torque, and engine RPM.
In some studies, the heat transfer analysis of the belt drive system was considered after determining the power loss or heat generation within the belt drive system, because the belt is an important component of the system and belt failures may occur because of the high temperature. For example, Merghache and Ghernaout focused on the heat transfer of a synchronous belt through numerical and experimental methods [13]. A general theory related to the rubber frictional heat against solids with arbitrary thermal properties was developed [14]. Mackerle concluded finite-element analyses, including sliding frictional simulations, applied to the rubber-like materials [15]. Furthermore, some authors focused on the temperature profile of a cylindrical disk-like rotor or brake disks, providing guidelines regarding the thermal analysis of the pulley in this case. Phan and Kondyles performed a conjugate computational fluid dynamics (CFD) analysis to accurately predict the thermal coefficient and thermal fading cycle of a rotor [16]. A 3D numerical analysis was performed to simulate the carbon nanotubes flow pasting an inclined rotating disk [17]. In addition, Talati and Jalalifar focused on the heat conduction of a brake disk and provided useful suggestions to reduce the surface temperature of the disk [18]. The naphthalene sublimation technique was adopted for the experimental heat and mass transfer analysis for a rolling wheel [19]. Belhocine and Bouchetara [20] analyzed the transient thermal behavior of a full and ventilated brake disk using a numerical approach. And an analytical and experimental analysis was conducted for a brake rotor with fins [21]. Furthermore, Reddy, Mallikarjuna, and Ganesan applied the same approach to investigate the influence of a vane shape on the flow field and heat transfer characteristics of a brake disk under various configurations and speeds [22]. Yan investigated the effect of cross-drilled holes on the cooling efficiency of brake disks [23]. A mesh-free method is used for the thermoelastic analysis of finite-length functionally graded cylinders reinforced by carbon nanotubes [24,25]. The same method is also used for thermoelastic behaviors during the optimization of thermomechanical responses of nanocomposite sandwich plates reinforced by carbon nanotubes [26,27] and nanocomposite cylinders reinforced with graphene [28]. By comparing the temperature results between standard and cross-drilled back disks, the theoretical and experimental research suggested that the heat dissipations from the surfaces of the cross-drilled holes played a prominent role in the overall heat transfer enhancement. Some research also focused on the temperature prediction of a belt drive system. In the simplest analytical method, the pulley was simplified as a basic round disk, and the Bessel function was used to calculate the temperature plot. In addition, the belt was a considered as a constant medium, and the heat exchange between the pulleys and belt was determined. However, this method was only applicable for a two-pulley-one-belt drive system [3]. In the latest research, numerical simulations were conducted, and the updated multiple reference frame (MLF) approach was utilized to calculate the surface temperature distributions of pulleys. This method was then extended to determine the temperature distribution of a two-pulley-one-belt continuously variable transmission (CVT) system [29,30].
In particular, the existing research focused on the methods and techniques related to the power loss determination and thermal analysis of the belt drive and its components. However, only a few studies have attempted to accelerate the thermal analysis of the belt drive system. In the design and optimization process, specifically, in the design and research of the belt drive system, several pulleys are usually considered to create various belt drive layouts operating under different conditions. However, the belt drive layouts and operating conditions change considerably more frequently than the changes in the pulley and its structures. The belt drive analysis thus involves a large number of operating conditions to be simulated, and in the existing techniques, the thermal analysis on each pulley must be performed for each condition, which is extremely time-consuming.
To overcome this limitation, this paper proposes a thermal model to efficiently determine the temperature plots of a belt drive system under various conditions. This model combines the advantages of both analytical and numerical analyses and provides accurate and efficient temperature maps for any belt drive system. Furthermore, an experiment was performed to validate the output temperatures under various operating conditions. The remaining paper is organized as follows: Section 2 illustrates the development of the thermal model, Section 3 describes the validation experiment for the proposed model, and Section 4 presents a comparison between the calculated and experimental values of the temperature corresponding to a variety of conditions. Finally, Section 5 presents the conclusions of this research.

Thermal Model
A set of thermal balance equations are established according to the thermal behavior inside the belt drive system to formulate a thermal model, which, along with the heat generation information from the existing power loss theory, can be used to calculate the temperature distribution. Two approaches are considered to create two different groups of thermal balance equations. The first approach is to use the analytical method, which is suitable for establishing the equations for the belt Appl. Sci. 2020, 10, 2709 4 of 23 heat dissipation and belt-pulley heat exchanges, because the belt has a simple structure. The second approach corresponds to the pulley thermal analysis. In this analysis, the temperature distribution of the pulleys under various conditions is examined, and a mathematical algorithm is formulated to describe the temperature deviations at different operating parameters. Subsequently, the numerical method is used to perform a thermal analysis for the baseline temperature distributions. This step provides the coefficients and functions for the pulley thermal balance equations. The next step involves simultaneously solving all the thermal balance equations to determine the belt temperature and heat dissipation of each pulley. Finally, the proposed model outputs the temperature plot of each pulley under the designed operating conditions by calculating the temperature changes based on the baseline numerical simulations. At this point, the model can effectively provide the temperature distributions of a belt drive system under a large number of conditions in a reasonable processing time.

Establishment of Thermal Model
The first step in the thermal model is to determine the heat flux and locations within the pulley system during the operation. The power loss theory is used to calculate the heat generation within the system because all the power losses are converted into heat, according to the conservation of energy. A prevalent theory is considered to calculate the power losses of the belt drive systems [3], and according to the locations of these losses, two forms of energy losses exist: The first form is the belt internal power loss P h , which is caused by the belt bending, stretching, and shear radial compression. These phenomena generate heat at different sections of the belt. As the belt operates at high speeds, these losses are considered as a heat source uniformly distributed along the belt. The second form is the contact surface power loss P fn , which occurs at the nth belt-pulley contact surface. Figure 2 shows an example of these losses, in which P h is uniformly distributed inside the belt, as shown by the dashed area, and P f1 and f2 are located at the belt-pulley contact areas, highlighted by the orange surfaces.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 23 approaches are considered to create two different groups of thermal balance equations. The first approach is to use the analytical method, which is suitable for establishing the equations for the belt heat dissipation and belt-pulley heat exchanges, because the belt has a simple structure. The second approach corresponds to the pulley thermal analysis. In this analysis, the temperature distribution of the pulleys under various conditions is examined, and a mathematical algorithm is formulated to describe the temperature deviations at different operating parameters. Subsequently, the numerical method is used to perform a thermal analysis for the baseline temperature distributions. This step provides the coefficients and functions for the pulley thermal balance equations. The next step involves simultaneously solving all the thermal balance equations to determine the belt temperature and heat dissipation of each pulley. Finally, the proposed model outputs the temperature plot of each pulley under the designed operating conditions by calculating the temperature changes based on the baseline numerical simulations. At this point, the model can effectively provide the temperature distributions of a belt drive system under a large number of conditions in a reasonable processing time.

Establishment of Thermal Model
The first step in the thermal model is to determine the heat flux and locations within the pulley system during the operation. The power loss theory is used to calculate the heat generation within the system because all the power losses are converted into heat, according to the conservation of energy. A prevalent theory is considered to calculate the power losses of the belt drive systems [3], and according to the locations of these losses, two forms of energy losses exist: The first form is the belt internal power loss PR hR , which is caused by the belt bending, stretching, and shear radial compression. These phenomena generate heat at different sections of the belt. As the belt operates at high speeds, these losses are considered as a heat source uniformly distributed along the belt. The second form is the contact surface power loss PR fnR , which occurs at the nth belt-pulley contact surface. Figure 2 shows an example of these losses, in which PR hR is uniformly distributed inside the belt, as shown by the dashed area, and PR f1R and PR f2R are located at the belt-pulley contact areas, highlighted by the orange surfaces. The next step is to create a thermal model by analyzing the heat flux flow behaviors through the belt drive components and dissipations to the environment at the exposed surfaces of the pulleys and the belt. A general model is first introduced in this section, followed by a description of a detailed subanalysis on each component from Section 2.2 to Section 2.4.
The first set of thermal balance equations corresponds to pulleys. No heat source exists inside pulleys, and the source of the heat dissipated on the pulley surfaces is the frictional heat and pulleybelt exchanged heat, which is transferred through the belt-pulley engaged surfaces. In this study, a variable, ξR nR, is introduced to represent the ratio of the thermal flux flowing into a pulley to the frictional heat flux between the nth pulley and the belt, as shown in Equation (1)  The next step is to create a thermal model by analyzing the heat flux flow behaviors through the belt drive components and dissipations to the environment at the exposed surfaces of the pulleys and the belt. A general model is first introduced in this section, followed by a description of a detailed subanalysis on each component from Sections 2.2-2.4.
The first set of thermal balance equations corresponds to pulleys. No heat source exists inside pulleys, and the source of the heat dissipated on the pulley surfaces is the frictional heat and pulley-belt exchanged heat, which is transferred through the belt-pulley engaged surfaces. In this study, a variable, ξ n , is introduced to represent the ratio of the thermal flux flowing into a pulley to the frictional heat flux between the nth pulley and the belt, as shown in Equation (1) Appl. Sci. 2020, 10, 2709

of 23
In addition, the heat dissipation at the pulley surfaces is strongly related to the temperature profiles, ambient temperature T a , geometry dimensions, rotational speed, material thermal conductivity, and other parameters. Hence, Equation (1) is transformed to where the temperature plot of the nth pulley is represented by the matrix T pn , the pulley geometry dimensions are stored in the matrix D pn , the rotational speed of nth pulley can be represented as ω pn , and the material thermal conductivity of nth pulley is Λ pn . The second set of thermal balance equations is related to the heat exchange at the belt-pulley engaged surface. The heat flux flowing into the pulley is influenced by both the frictional heat and heat exchange at the belt-pulley engaged surface. The temperature difference between the belt and pulley causes a certain amount of heat exchange. This condition can be represented as in the following simple equation: ξ n P f n = g(T pn , T b , P f n ) where the contact temperature on the belt side is T b , and the contact temperature on the pulley side is T pn , which belongs to the matrix T pn . The third thermal balance equation corresponds to the heat analysis of the belt, which is in contact with all the pulleys. The heat dissipated from the belt surfaces is equal to the total heat generation minus the heat dissipated on the pulley surfaces. This heat is also strongly related to the ambient temperature, belt geometry, speed, and other parameters, which can be simply represented by a function, as follows: where the belt temperature plots belong to the matrix T b , belt geometry dimensions are stored in the matrix D b , belt velocity can be represented as V b , and the material thermal conductivity is Λ b . N is the total number of pulleys in the belt drive system. By combining Equations (2)-(4), a simple thermal model can be established, as follows: The subsequent subsections describe the analysis of the thermal flows inside the components and belt-pulley heat exchanges to provide clear information regarding the functions f (T pn , . . . ), g(T pn , . . . ) and k(T b , . . . ) and solve this set of equations.

Thermal Analysis of the Belt
The belt is a long and thin structure, and because of its approximately rectangular cross section, the analytical method is applicable to perform its heat transfer analysis. The temperature distribution of the belt can be considered uninformed because it has a thickness of as low as 0.05 m but a thermal conductivity of as high as 0.8 W/(m·K). The experiments described in Section 3 demonstrate that the two sides of the belt have similar temperatures, as shown in Figure 3a,b. The temperature distribution obtained using a simple numerical simulation also supports this assumption, as shown in Figure 3c.
Hence, the temperature distribution inside the belt can be replaced by a single value T b . All the values in matrix T b thus become T b . The heat dissipation on the belt surfaces can be easily calculated [31], thereby yielding function k(T b , . . . ). The thermal balance of the belt can be expressed as where A b is the total belt surface area exposed to the environment, h b is the heat transfer coefficient of the belt surfaces, and T a is the constant ambient temperature in this static analysis.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 24 conductivity of as high as 0.8 W/(m⋅K). The experiments described in Section 3 demonstrate that the two sides of the belt have similar temperatures, as shown in Figure 3a,b. The temperature distribution obtained using a simple numerical simulation also supports this assumption, as shown in Figure 3c. Hence, the temperature distribution inside the belt can be replaced by a single value TR bR. All the values in matrix TR bR thus become TR bR. The heat dissipation on the belt surfaces can be easily calculated [31], thereby yielding function k(TR bR ,…). The thermal balance of the belt can be expressed as where AR bR is the total belt surface area exposed to the environment, hR bR is the heat transfer coefficient of the belt surfaces, and TR aR is the constant ambient temperature in this static analysis. In Equation (6), hR bR varies with the local air speed. The relationship between hR bR and the local air velocity VR bR can be represented through a fourth-order polynomial created by the curve-fitting method: where kR 1R , kR 2R , kR 3R, kR 4R , and kR 5R denote five coefficients whose values related to the material properties of the belt. In this study, hR bR is measured from a small wind tunnel device, and the coefficients for the belt are as follows: kR 1R = −6.93 × 10P −5 P , kR 2R = 0.00714, kR 3R = 0.252, kR 4R = 4.3 and kR 5R = 1.7. This model focuses on the static condition, meaning that the belt speed is constant under each scenario. Equation (7) is used to calculate the corresponding hR bR under the designed belt speed, which is then substituted into Equation (6) for the model calculation. In Equation (6), h b varies with the local air speed. The relationship between h b and the local air velocity V b can be represented through a fourth-order polynomial created by the curve-fitting method: where k 1 , k 2 , k 3 , k 4 , and k 5 denote five coefficients whose values related to the material properties of the belt. In this study, h b is measured from a small wind tunnel device, and the coefficients for the belt are as follows: k 1 = −6.93 × 10 −5 , k 2 = 0.00714, k 3 = 0.252, k 4 = 4.3 and k 5 = 1.7. This model focuses on the static condition, meaning that the belt speed is constant under each scenario. Equation (7) is used to calculate the corresponding h b under the designed belt speed, which is then substituted into Equation (6) for the model calculation.

Thermal Analysis at Belt-Pulley Engaged Surfaces
The heat analysis at the belt-pulley engaged surfaces consists of two simultaneous phenomena, namely, the frictional heat and heat exchange. When the pulley temperature at the pulley-belt engaged surface T pcn is equal to T b , no heat exchange occurs between the belt and pulley; the pulley acquires half of the frictional heat, and thus ξ n = 0.5. When T b > T pcn and some heat flows from the belt to the pulley, the pulley acquires additional heat flux; thus, ξ n > 0.5. Thus, the following equation can be established: where Φ exn is the heat exchange flux per second at the nth belt-pulley engaged surface from the belt to the nth pulley. Moreover, ΦR exn R is also related to the pulley-belt temperature difference (T b − T pcn ), as follows: Combining Equations (8) and (9) yields Equation (10), which provides the definition of the function g(T b , T pcn , P fn ) expressed in Equation (5):

Pulley Internal Thermal Analysis
Formulating the relationship between Φ pn and the input parameters of function f (T pn , . . . ) for a general pulley is a dedicated task, as the pulley could have any shape, and its structure is likely complex. In such cases, the optimal analysis method involves numerical simulations because they can be applied to any pulley shape, enabling the proposed method to be free from the limitations of the complex pulley structure. Furthermore, the numerical method can be used to simulate the turbulence flows caused by the pulley rotation, providing the local air velocities for accurately calculating the heat dissipations, because the heat dissipation coefficient is mainly related to the local air velocity.
However, the numerical method is time-consuming, and it is nearly impossible to calculate the pulley temperature plot for each condition. Thus, in this work, the thermal flow within the pulley system is investigated to examine how the operating parameters influence the change in Φ pn and T pn . Through these observations, a mathematical algorithm can be formulated to describe the relationship between Φ pn and the parameters. Subsequently, thermal balance equations based on this relationship can be established, and the corresponding coefficient can be determined using the numerical simulations.

Pulley Mathematical Algorithm
To formulate a mathematical algorithm for a complex pulley, complex features such as grooved surfaces and some unimportant fillets are suppressed. Figure 4 shows the geometry transition from a general pulley to a simplified pulley. The mathematical model is established on the main round I-structure and the fins are later considered.  The main round I-structure can be divided into three layers ( Figure 5) to better analyze the thermal conduction inside the pulley. In each layer, the cross section is considered as a rectangular shape. The temperature distribution can be described using the Bessel function: The main round I-structure can be divided into three layers ( Figure 5) to better analyze the thermal conduction inside the pulley. In each layer, the cross section is considered as a rectangular shape. The temperature distribution can be described using the Bessel function: where T(r) is the temperature distribution along the radial direction of the pulley. The main round I-structure can be divided into three layers ( Figure 5) to better analyze the thermal conduction inside the pulley. In each layer, the cross section is considered as a rectangular shape. The temperature distribution can be described using the Bessel function: where T(r) is the temperature distribution along the radial direction of the pulley. However, Layer 2 includes fins ( Figure 6). In an infinitesimal element dr in Layer 2, the fins amplify the heat conduction and dissipation rate on the areas connected to the fins (AR fR in Figure 6) compared to those in the areas not connected to the fins (AR nR in Figure 6). A new coefficient ηR fR is used to represent the ratio of the heat dissipation rate on surfaces AR fR has been increased ηR fR times by fins comparing with the original rate. If no fins exist, the dissipation rates on AR fR and AR nR are equivalent, and the dissipation rate on one side of this infinitesimal element dr is 2πr⋅dr⋅hR pR. However, when the fins are considered in the analysis, the dissipation rate on one side changes to [NR fR ⋅WR fR⋅dr⋅ηRfR⋅hR pR +(2πr ⋅dr-NR fR ⋅WR fR⋅dr)⋅hR pR ]. Therefore, the heat dissipation rate is amplified by a factor of [1 + NR fR ⋅WR fR ⋅(ηR fRR R-1)/2πr]. It also can be considered that hR pR on the side surfaces of Layer 2 is increased by this ratio, and Equation (11) can be modified as However, Layer 2 includes fins ( Figure 6). In an infinitesimal element dr in Layer 2, the fins amplify the heat conduction and dissipation rate on the areas connected to the fins (A f Rin Figure 6) compared to those in the areas not connected to the fins (A n in Figure 6). A new coefficient η f is used to represent the ratio of the heat dissipation rate on surfaces A f has been increased η f times by fins comparing with the original rate. If no fins exist, the dissipation rates on A f and A n are equivalent, and the dissipation rate on one side of this infinitesimal element dr is 2πr·dr·h p . However, when the fins are considered in the analysis, the dissipation rate on one side changes to [N f ·W f ·dr·η f ·h p +(2πr·dr-N f ·W f ·dr)·h p ]. Therefore, the heat dissipation rate is amplified by a factor of [1 + N f ·W f ·(η f -1)/2πr]. It also can be considered that h p on the side surfaces of Layer 2 is increased by this ratio, and Equation (11) can be modified as The value of η f must be calculated based on the fin geometry during the solution of Equation (11). The fin section in a rounded infinitesimal element dr in Layer 2 is considered as a rectangle beam. The heat dissipation on the exposed surfaces of this beam is equal to the heat flux through the fin connected area A f . The general thermal balance equation for this long beam is [32] where G(x) is the temperature distribution of the beam along the pulley axial direction. G(0) is the temperature at A f , which is equal to T(r), the temperature at radius r on the main I-structure. Hence, the total heat dissipation for each fin can be calculated by integrating the surface temperature.
The first and second parts of this equation correspond to the heat dissipations on the two sides and the far-end side (A e in Figure 6), respectively. According to the definition, After the integration of G(x), which is a solution of an ODE, η f at the radius r can be calculated as η f can be considered as a constant, because Equation (16) shows that it is only a function of the fin geometry (L f , W f ), thermal conductivity (Λ p ) and heat transfer coefficient (h p ). The value of ηR fR must be calculated based on the fin geometry during the solution of Equation (11). The fin section in a rounded infinitesimal element dr in Layer 2 is considered as a rectangle beam. The heat dissipation on the exposed surfaces of this beam is equal to the heat flux through the fin connected area AR fR. The general thermal balance equation for this long beam is [32] where G(x) is the temperature distribution of the beam along the pulley axial direction. G(0) is the temperature at AR fR , which is equal to T(r), the temperature at radius r on the main I-structure. Hence, the total heat dissipation for each fin can be calculated by integrating the surface temperature.
The first and second parts of this equation correspond to the heat dissipations on the two sides and the far-end side (AR eR in Figure 6), respectively. According to the definition, After the integration of G(x), which is a solution of an ODE, ηR fR at the radius r can be calculated as Because η f is independent of r, Equation (17) can be regarded as Kummer's equation [33]. Its analytical solution T(r) provides the temperature distribution of one layer considered in this study, as follows: where the constant coefficients C 1 and C 2 can be derived from the thermal analysis boundary conditions considered in this study. When r = R y , value of function F 2 (R y ) in Equation (17) will become an infinity. However, the temperature at R y in real world cannot be an infinite high temperature. Therefore, C 2 has to be set to 0. Then, C 1 can be calculated as The temperature distribution for one layer can be obtained as Equation (19) shows that T(r) is related to the temperature at the outer radius of this layer, R x . The temperature distribution for a pulley can thus be calculated by using Equation (19) three times from Layers 1 to 3. The connection of any two layers exhibits the same temperature, and T(R y ) in one layer is thus equal to T(R x ) in the next layer. Thus, the three-layer temperature distribution of the pulley is related to the temperature at the outer radius of this pulley, T pcn .
Subsequently, the total heat dissipated per second at the side surfaces of all the three layers (surfaces B 1 , B 2 , and B 3 in Figure 7), Φ ps , is where R o and R i denote the outer and inner radii of the investigated pulley, and R 1 and R 2 denote the inner radii of Layers 1 and 2, as shown in Figure 5. The heat dissipation rate at the two side surfaces in one layer is Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 23 material and the surface at the inner diameter of the pulley (surface indicated as AR 3R in Figure 7). Therefore, the total heat dissipation rate of the investigated pulley is The above mathematic algorithm can then be used to formulate the mathematical relationship between ΦR pnR and the input parameters of function f(TR pnR,R…R ). The number of these input parameters is more than ten, and these parameters are thus classified into two groups to enable a clearer analysis. The first group includes the pulley related parameters, which are geometry related and stored in matrix DR pnR . In particular, this group includes NR fR , LR fR, WR fR, HR fR , AR 0R , AR 1R , AR 2R, AR 3R , AR fcR , WR 1R, WR 2R, WR 3R , RR oR , RR 1R, RR 2R , RR iR and the material related parameters ΛR pnR. The second group includes the operation related parameters hR pnR , ωR pnR and TR pcnR. hR pnR is a function of ωR pnR , which is determined experimentally.
In this case, ΦR pnR is a function of only two parameters ωR pnR and TR pcnR for a specific pulley, because the parameters in the first group are determined by the type of pulley and not the operating condition. In this work, a pulley thermal coefficient λR pnR is introduced to represent the ability of heat dissipation of the nth pulley in the belt drive system, which is the ratio of ΦR pnR to TR pcnR. And then each investigated pulley has its own thermal coefficient λR pR . The physical meaning of λR pnR is when TR pcnR remains same, the higher λR pnR is, the larger heat flux dissipates from pulley surfaces. In this thermal analysis, this coefficient can reduce the complexity of the thermal calculation.
( ) pn pcn a n fn The mathematical algorithm illustrates the changes in λR pnR for various parameters ωR pnR and TR pcnR. When a pulley operates at a constant rotating speed, the severity of the turbulent flow and Reynolds number remains the same, and hR pnR remains constant. In this case, ΦR pnR is directly proportion to TR pcnR for a pulley, and λR pnR does not change. However, if the pulley rotation becomes faster, the severity of the turbulent flow and Reynolds number increase dramatically [34], and hR pnR increases. Thus, λR pnR changes for different values of ωR pnR , and λR pnR is a function of ωR pnR for the nth pulley in the system. This function is defined as SR nR (ωR pnR ). Therefore, Equation (24) No analytical method can be used to analyze the turbulence flow. Thus, in this work, a numerical method is used to calculate the data points of λR pnR at different ωR pnR from 2000 RPM to 6000 RPM in intervals of 1000 RPM, and SR nR (ωR pnR ) is determined by the curve-fitting of these points. Moreover, the Furthermore, this algorithm includes the heat dissipated from the other surfaces, Φ pb : where A 0 , A 1 , A 2 , and A 3 denote the surface areas, as shown in Figure 7, α is the wrap angle, which is the degrees that the belt contacts the pulley, and h s is the heat transfer coefficient between the shaft material and the surface at the inner diameter of the pulley (surface indicated as A 3 in Figure 7). Therefore, the total heat dissipation rate of the investigated pulley is The above mathematic algorithm can then be used to formulate the mathematical relationship between Φ pn and the input parameters of function f (T pn , . . . ). The number of these input parameters is more than ten, and these parameters are thus classified into two groups to enable a clearer analysis. The first group includes the pulley related parameters, which are geometry related and stored in matrix D pn . In particular, this group includes N f , L f , W f , H f , A 0 , A 1 , A 2 , A 3 , A fc , W 1 , W 2 , W 3 , R o , R 1 , R 2 , R i and the material related parameters Λ pn . The second group includes the operation related parameters h pn , ω pn and T pcn . h pn is a function of ω pn , which is determined experimentally.
In this case, Φ pn is a function of only two parameters ω pn Rand T pcn for a specific pulley, because the parameters in the first group are determined by the type of pulley and not the operating condition. In this work, a pulley thermal coefficient λ pn is introduced to represent the ability of heat dissipation of the nth pulley in the belt drive system, which is the ratio of Φ pn to T pcn . And then each investigated pulley has its own thermal coefficient λ p . The physical meaning of λ pn is when T pcn remains same, the higher λ pn is, the larger heat flux dissipates from pulley surfaces. In this thermal analysis, this coefficient can reduce the complexity of the thermal calculation.
ξ n P f n = λ pn T pcn − T a (24) The mathematical algorithm illustrates the changes in λ pn for various parameters ω pn and T pcn . When a pulley operates at a constant rotating speed, the severity of the turbulent flow and Reynolds number remains the same, and h pn remains constant. In this case, Φ pn is directly proportion to T pcn for a pulley, and λ pn does not change. However, if the pulley rotation becomes faster, the severity of the turbulent flow and Reynolds number increase dramatically [34], and h pn increases. Thus, λ pn changes for different values of ω pn , and λ pn is a function of ω pn for the nth pulley in the system. This function is defined as S n (ωR pn ). Therefore, Equation (24) becomes No analytical method can be used to analyze the turbulence flow. Thus, in this work, a numerical method is used to calculate the data points of λ pn at different ω pn from 2000 RPM to 6000 RPM in intervals of 1000 RPM, and S n (ω pn ) is determined by the curve-fitting of these points. Moreover, the numerical simulations are also used to provide a baseline temperature database for the thermal model to output the thermal distributions of the pulleys considered. This aspect is discussed in more detail in Section 2.4.2.

Numerical Method
The numerical method uses a different approach from that of the mathematical algorithm to provide the λ pn for the nth pulley at different rotating conditions ω pn for each pulley. This section introduces the modeling, meshing, boundary conditions, and other general settings for the numerical calculation.

Modeling Structure
Both solid and fluid regions are created in this numerical study as shown in Figure 8. The solid region is created for the pulley and its connected shaft to simulate the thermal conduction inside the pulley and shaft. In addition, a rounded fluid region encloses the solid region. The diameter of the fluid region is three times larger than the pulley diameter, which provides a sufficiently large fluid domain to simulate the turbulence and small vortex flows around the pulley. The structure of the contacted belt is not considered. Furthermore, the frictional heat caused by the belt slipping against the pulley is applied through the boundary conditions, instead of modeling the slip motions to acquire the heat in the simulation.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 23 introduces the modeling, meshing, boundary conditions, and other general settings for the numerical calculation.

Modeling Structure
Both solid and fluid regions are created in this numerical study as shown in Figure 8. The solid region is created for the pulley and its connected shaft to simulate the thermal conduction inside the pulley and shaft. In addition, a rounded fluid region encloses the solid region. The diameter of the fluid region is three times larger than the pulley diameter, which provides a sufficiently large fluid domain to simulate the turbulence and small vortex flows around the pulley. The structure of the contacted belt is not considered. Furthermore, the frictional heat caused by the belt slipping against the pulley is applied through the boundary conditions, instead of modeling the slip motions to acquire the heat in the simulation.

Modeling Motions
A single reference frame (SRF) is adopted to model the relative rotation of the pulley and shaft against the surrounding fluid, which is considered as air in this work. The use of the SRF is suitable for such a static state application. The fluid domain is rotated against the stable solid domain to simulate the relative motion between the fluid and solid regions. In this manner, the fluid domain is used only to simulate the turbulence flow around pulley, and the solid domain is used only to simulate the heat conduction inside the pulley. The scale of the numerical simulation is decreased from a fluid-structure interaction to a static analysis, thereby reducing the required computational

Modeling Motions
A single reference frame (SRF) is adopted to model the relative rotation of the pulley and shaft against the surrounding fluid, which is considered as air in this work. The use of the SRF is suitable for such a static state application. The fluid domain is rotated against the stable solid domain to simulate the relative motion between the fluid and solid regions. In this manner, the fluid domain is used only to simulate the turbulence flow around pulley, and the solid domain is used only to simulate the heat conduction inside the pulley. The scale of the numerical simulation is decreased from a fluid-structure interaction to a static analysis, thereby reducing the required computational time and resources. Moreover, the SRF is applicable for the situations in which the fluid and solid domains do not change their shapes during the analysis.
Although several other options are available to supplement this analysis, they are not suitable to be applied for the analysis considered in this work. One option, known as the mixing plane approach, uses a circumferential averaging technique to compute a single radial profile, which in turn is used for all the cell faces at the interface. This technique is implemented to accommodate rotor/stator flows in centrifugal and axial flow processors. Another option is the so-called sliding mesh approach, which accounts for the unsteady interactions owing to the relative motions of the stationary and rotating components. The third option is the dynamic mesh approach, which can move boundaries and/or objects to adjust the mesh accordingly. These time-consuming options are suitable only for unsteady and/or transient analysis, which is in contrast to the considered analysis. Hence, the use of the SRF method was preferable for the considered analysis and used to conduct the numerical simulations.

Meshing Procedure
The solid and fluid regions are meshed simultaneously using the same technique. A triangular surface mesh is generated and converted to polyhedrons to reduce the number of cells. A minimum edge length of 0.2 mm with a growth rate of 1.1 for the mesh is considered as sufficient to resolve the given geometry. The fine mesh is focused on the surfaces at the pulley outer radius and gradually becomes coarser as it approaches the pulley inner radius (Figure 9). In addition, the interfaces for each surface are defined to ensure a conformal mesh at the solid/fluid connected boundaries. The conformal meshes share the same nodes and element edges at the boundaries, facilitating the numerical algorithms to calculate the heat exchange while enhancing the computational stability.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 23 becomes coarser as it approaches the pulley inner radius (Figure 9). In addition, the interfaces for each surface are defined to ensure a conformal mesh at the solid/fluid connected boundaries. The conformal meshes share the same nodes and element edges at the boundaries, facilitating the numerical algorithms to calculate the heat exchange while enhancing the computational stability.

Modeling the Heat Transfer
Establishing the settings for the heat transfer in this numerical simulation is challenging, as this setting requires an accurate turbulence model to not only simulate the airflow caused by the pulley rotation but also calculate the heat dissipation rate on the pulley surfaces. The model used in this numerical simulation is the Menter Shear Stress Transport (SST) k-ω turbulence model, which belongs to the Reynolds Averaged Navier-Stokes (RANS) regime. This two-equation eddy-viscosity model is ideal for this analysis because of its good near wall treatment, acceptable stability and reasonable computational time. The fluid medium is air, and the temperature must be below the melting temperature of FRPs which is 140 °C. The maximum rotational speed of the pulley is 7000 RPM, and the maximum diameter is 108 mm, resulting in subsonic analysis within the fluid domain. Hence, the density of air is computed by using the incompressible ideal gas equation and its change is negligible. The finite volume method used in this analysis solves the governing equations in integral form, as shown below [31]: j j ∂ Figure 9. Polyhedral mesh for the numerical analysis.

Modeling the Heat Transfer
Establishing the settings for the heat transfer in this numerical simulation is challenging, as this setting requires an accurate turbulence model to not only simulate the airflow caused by the pulley rotation but also calculate the heat dissipation rate on the pulley surfaces. The model used in this numerical simulation is the Menter Shear Stress Transport (SST) k-ω turbulence model, which belongs to the Reynolds Averaged Navier-Stokes (RANS) regime. This two-equation eddy-viscosity model is ideal for this analysis because of its good near wall treatment, acceptable stability and reasonable computational time. The fluid medium is air, and the temperature must be below the melting temperature of FRPs which is 140 • C. The maximum rotational speed of the pulley is 7000 RPM, and the maximum diameter is 108 mm, resulting in subsonic analysis within the fluid domain.
Hence, the density of air is computed by using the incompressible ideal gas equation and its change is negligible. The finite volume method used in this analysis solves the governing equations in integral form, as shown below [31]: where W is the vector of conserved quantities, which are mass, momentum and energy. J represents the inviscid terms. The vector of viscous terms is given by K. The vector H is the body forces. In these vectors, ρ denotes the density of air, υ is the vector of the air velocity, υR g R is the velocity vector of the grid, E represents the total energy per unit mass, p is the pressure, and . q is the occurring heat flux. The additional SST k-ω model transport equations, which are turbulence kinetic energy and specific dissipation rate equations, are also necessary [31]: where Γ k and Γ ω is the effective diffusivity of k and ω. G k and G ω represents the kinetic energy of k and ω. The dissipation of k and ω resulting from turbulence is given by Y k and Y ω . The extra term, D ω , is a blend function which combines the k-ω model for the robust near wall treatment and k-ε model for the free stream in the far field. The constants in these equations have been remained as the default. Moreover, a second-order upwind scheme is used to realize the spatial discretization of the governing equations. In addition, pseudo transient, a type of implicit under-relaxation for steady-state cases, is activated for better stabilizing the case and ensuring rapid convergence. In terms of the analysis parameters, the initial temperature is the ambient temperature. The thermal conductivity for the pulley material is measured according to ISO 22007-2 to accurately simulate the heat conduction.

Input/Output Setting
The input value of the thermal flux Φ pn flowing into the pulley model cannot be defined directly. In one full cycle of pulley rotation, the highlighted red surfaces at the crest of the grooved surface engage with the belt surfaces, frictional heat is generated, and a portion of this heat flows into the pulley. The flow direction is toward the pulley. Then, the surface disengages with the belt surfaces when the pulley rotates by a certain degree, and heat is dissipated through this surface. The flow direction at this instant is toward the outside of the pulley. It is impossible to change the flow direction based on the definition of the boundary conditions during the pulley high-speed revolutions. Hence, small solid sections (blue area in Figure 10) close to the surface are created to simulate the frictional heat generation. This section is close to the surface, and its influence on the temperature distribution at the pulley can be neglected. The highlighted red surfaces in Figure 10, together with the rest of the pulley surfaces, are set to dissipate the heat into the fluid domain.
together with the rest of the pulley surfaces, are set to dissipate the heat into the fluid domain.
The thermal flux of the pulley ΦR pnR is dissipated into the environment and the end of the shaft. The boundary of the fluid domain is an open space, and hence it is defined as a pressure outlet with zero-gauge pressure. The backflow turbulent intensity is 5%, and the backflow turbulent viscosity ratio is 10. The boundary for the end of the shaft is set as a constant temperature, which is the ambient temperature.  The thermal flux of the pulley Φ pn is dissipated into the environment and the end of the shaft. The boundary of the fluid domain is an open space, and hence it is defined as a pressure outlet with zero-gauge pressure. The backflow turbulent intensity is 5%, and the backflow turbulent viscosity ratio is 10. The boundary for the end of the shaft is set as a constant temperature, which is the ambient temperature.
The output results of the numerical method provide abundant information, for instance, the temperature distribution of the pulley T pn , area-weight average temperature on the pulley contact surfaces T pcn , and heat dissipation flux on the pulley exposed surfaces Φ pn . The primes on these parameters denote that this information is based on the baseline conditions and not the predicted temperatures and information under the designed operating conditions.

Parametric Setting
According to the definition of λ pn , the change of T pcn and ω pn influences the value of λ pn . Figure 11a presents the comparison of λ pn obtained from the numerical and mathematical methods for a pulley with a radius of 76.2 mm under various contacted temperatures of pulley T pcn at the speed of 6000 RPM. The diagram demonstrates that the numerical approach, which is the same as the mathematical approach, can prove that λ pn is irrelevant with T pcn . Hence, the following is to investigate the only impact of ω pn on λ pn . T , and heat dissipation flux on the pulley exposed surfaces ' pn Φ . The primes on these parameters denote that this information is based on the baseline conditions and not the predicted temperatures and information under the designed operating conditions.

Parametric Setting
According to the definition of λR pnR , the change of TR pcnR and ωR pnR influences the value of λR pnR. Figure 11a presents the comparison of λR pn Robtained from the numerical and mathematical methods for a pulley with a radius of 76.2 mm under various contacted temperatures of pulley TR pcnR at the speed of 6000 RPM. The diagram demonstrates that the numerical approach, which is the same as the mathematical approach, can prove that λR pn Ris irrelevant with TR pcnR . Hence, the following is to investigate the only impact of ωR pnR on λR pnR . This numerical simulation formulates the function SR nR (ωR pnR) or the relationship between λR pnR and ωR pnR through parametric studies in which ωR pnR is increased in intervals of 1000 RPM. In each condition of ωR pnR , the heat generation is input into the numerical model, and the output includes ' pn Φ and ' pcn T . Hence, a curve of SR nR (ωR pnR ) is plotted to describe this relationship. Figure 11b presents the comparison of λR pn Robtained from the numerical and mathematical methods for a pulley with a radius of 76.2 mm under various rotational speeds ranging from 2000-6000 RPM at TR pcnR is 50 °C. The λR pn Rvalues obtained using both the methods increase with the increase in ωR pnR, indicating that the heat dissipation on the pulley surfaces increases at higher speeds and then improving the heat dissipation ability of a pulley.
Meanwhile, the λR pn Rvalues obtained using the numerical algorithm are generally higher than those from the mathematical method because in the case of the mathematical algorithm, the grooved feature at the pulley outer surface was simplified to a flat surface (Figure 4), thus reducing the surface area AR 0R (Figure 7) and further decreasing the surface dissipation ΦR pnR through Equation (22) to cause λR pnR having a lower value than in the numerical approach. The heat dissipation from the surface area AR 0R (The first right term in Equation 22) is the primary contribution for ΦR pbR because T(RR 0R ) has the highest temperature within the pulley. Plus, the mathematical approach also does not consider the influence of the severity of turbulence with increased speeds, making the curve of SR nR (ωR pnR) from the mathematical approach less tilted than that from the numerical results. Therefore, the SR nR (ωR pnR ) from the numerical approach, instead of the mathematical approach, is adopted in this study. Therefore, function f(TR pnR ,R …R) can be solved by using this numerical simulation. SR nR(ωR pnR ) is determined by using the thermal balance equation, Equation (24), which is integrated into the thermal This numerical simulation formulates the function S n (ω pn ) or the relationship between λ pn and ω pn through parametric studies in which ω pn is increased in intervals of 1000 RPM. In each condition of ω pn , the heat generation is input into the numerical model, and the output includes Φ pn and T pcn . Hence, a curve of S n (ωR pn ) is plotted to describe this relationship. Figure 11b presents the comparison of λ pn obtained from the numerical and mathematical methods for a pulley with a radius of 76.2 mm under various rotational speeds ranging from 2000-6000 RPM at T pcn is 50 • C. The λ pn values obtained using both the methods increase with the increase in ω pn , indicating that the heat dissipation on the pulley surfaces increases at higher speeds and then improving the heat dissipation ability of a pulley.
Meanwhile, the λ pn values obtained using the numerical algorithm are generally higher than those from the mathematical method because in the case of the mathematical algorithm, the grooved feature at the pulley outer surface was simplified to a flat surface (Figure 4), thus reducing the surface area A 0 (Figure 7) and further decreasing the surface dissipation Φ pn through Equation (22) to cause λ pn having a lower value than in the numerical approach. The heat dissipation from the surface area A 0 (The first right term in Equation 22) is the primary contribution for Φ pb because T(R 0 ) has the highest temperature within the pulley. Plus, the mathematical approach also does not consider the influence of the severity of turbulence with increased speeds, making the curve of S n (ωR pn ) from the mathematical approach less tilted than that from the numerical results. Therefore, the S n (ωR pn ) from the numerical approach, instead of the mathematical approach, is adopted in this study.
Therefore, function f (T pn , . . . ) can be solved by using this numerical simulation. S n (ω pn ) is determined by using the thermal balance equation, Equation (24), which is integrated into the thermal model. A set of values {S n (ω pn ), T pn ,Φ pn } is created for each pulley rotating at speeds in increments of 1000 RPM, and a database is established for all the available pulleys. The model can then import S n (ωR pn ) from this database based on the pulley and designed ω pn . Moreover, T pn and Φ pn are stored for subsequent pulley thermal distribution calculations.

Overall Structure
Equations (6), (10), and (25) can be combined to yield the following thermal model for an arbitrary multi-pulley belt drive system: The belt drive system generally has a maximum of 10 pulleys (N ≤ 10) depending on the design and engine requirements. Equation (26) is a system of equations with a maximum of 21 unknown variables (T p1 , · · · , T pN , ξ 1 , · · · , ξ N and T b ). This study adopts the LU decomposition to solve Equation (29) to obtain the thermal flow information for each component and belt temperature. The calculation time is only a few seconds by using a computer with 12-core 3.4 GHz CPU.
In addition, the temperature distribution of the nth pulley in this system, T pn , can be determined by suitably modifying the referenced numerical solutions. First, the database stores the two sets of {T pn ,Φ pn } under the upper and lower speeds of ω pn . This model uses the linear interpolation to calculate Φ pn and T pn at the designed speed ω pn . Then, the temperature field of the nth pulley is

Experimental Setup
The reliability of the proposed thermal model in terms of the accuracy was validated experimentally. In particular, an experiment was performed to measure the temperature distribution of the entire belt drive and compare it with the corresponding results obtained from the thermal model under various operating conditions. Therefore, a rigorous experimental arrangement was established before validating the thermal algorithm with the multi-pulley configuration, as described in this section.

Load Cases
Considering a wide range of test operating conditions of the belt drive system is vital to validate the calculated results from the thermal model. The types of pulley for the serpentine belt drive system vary considerably in terms of the radius and materials. To cover the range of potential pulley variations, three radii and two materials were selected for the experiments. The engine operation varies significantly from idle to the maximum RPM conditions, and hence, three speeds for the DN pulley were selected. Table 1 summarizes the test configurations to fully consider the significant variation in the engine operating conditions. The transmitted speed and torque of the belt drive system were controlled by altering the revolution of the DR pulley and torque load of the DN pulley. The temperature distributions of the belt system were measured under the designed working conditions for validating the thermal model using the above devices and equipment as a test facility.

Belt Layouts
Two layouts were designed to consider the various changes in the parameters listed in Table 1. The DR pulley was larger than the DN pulley in a typical front-end accessory drive system; hence, the DR pulley had a radius of 148 mm, and the DN pulley had a radius of 76.2 mm and 108 mm in the two layouts. In addition, the DR pulley was made of steel alloy, and the DN pulleys were made of FRP. Moreover, the layout had to contain five pulleys instead of only the DN and DR pulleys. The extra three pulleys were used to measure and maintain the belt tension and speed, especially when the transmitted speed for a pulley was as high as 6000 RPM. Table 2 lists the pulley locations for the layouts of the two-belt drive systems for this experiment.

Equipment Setup
The designed five-pulley belt drive system was installed on an engine dynamometer ( Figure 12). The bottom left pulley was the DR pulley, marked as pulley 1, and the top middle pulley was the DN pulley, marked as pulley 4. The remaining pulleys were as follows: pulley 2 connected to a hub load sensor for belt tension measurement; pulley 3 connected to a speedometer for belt speed measurement; and pulley 5 connected to a tensioner that reduces the vibration of the belt and maintained a constant belt tension. Moreover, an insulated chamber equipped with a temperature control system is installed on this engine dynamometer to maintain the constant temperature. Two thermocouples were mounted at the two sides of the chamber to monitor the internal temperature. A case fan present inside the chamber ensured that the air temperature was uniformly distributed. In addition, thermal measurement devices were installed close to the measurement locations. Five infrared cameras were installed to measure the temperatures at the selected locations. A thermal imaging camera targeting the DN pulley was installed to capture the temperature distribution of this pulley. Furthermore, an ancillary data acquisition system was used to record the temperature values measured from the sensors. The thermal model used the abovementioned information related to the belt layout and operating conditions, along with the pulley geometries and thermal properties (hR pR , ΛR pR ,  ) measured using standard procedures, to calculate the temperature distributions of the belt drive system under the designed operating conditions. The comparison between the experimental and analytical results is presented in the following section.

Results and Discussion
To demonstrate the accuracy and reliability of this model, the experimentally obtained temperature and calculated temperature under various pulley radii, material types, and speeds, as described in Table 1, are compared. First, the reliability of λR pnR is validated by comparing the temperature distributions of the pulleys obtained using the numerical method with the experimental measurements. Second, the thermal model is validated by comparing the numerically and experimentally obtained temperatures of the two designed belt drives at different locations within the system. Finally, pulley temperature distribution results are validated by obtaining the experimental measurements on different locations of the pulleys.

Verification of Sn(ωpn)
Obtaining an accurate curve for SR nR (ωR pnR ) is critical during the calculation of the thermal model. According to the definition of λR pnR, the numerical and measured results of the heat dissipation flux of pulleys must be compared under the identical the temperatures at the contact area TR pcnR and rotating speed ωR pnR. However, no device can measure the heat dissipation of the pulley at its high rotational speeds. Hence, it is reasonable to compare the calculated and measured results of the temperature distribution of the pulley, because the temperature drop along the radial direction is strongly linked to the heat dissipation rate on pulley surfaces. Figure 13 presents that the temperature distributions of Layer 2 for the pulleys from the numerical model approach and compare this with the experimental temperature data at five different points along the radial direction under the condition of TR pcnR = 55 °C and rotating speed ωR pn R= 4000 RPM. In this diagram, the calculated results have a good agreement with the measured temperatures, indicating that the heat dissipation for a pulley is The thermal model used the abovementioned information related to the belt layout and operating conditions, along with the pulley geometries and thermal properties (h p , Λ p , · · · ) measured using standard procedures, to calculate the temperature distributions of the belt drive system under the designed operating conditions. The comparison between the experimental and analytical results is presented in the following section.

Results and Discussion
To demonstrate the accuracy and reliability of this model, the experimentally obtained temperature and calculated temperature under various pulley radii, material types, and speeds, as described in Table 1, are compared. First, the reliability of λ pn is validated by comparing the temperature distributions of the pulleys obtained using the numerical method with the experimental measurements. Second, the thermal model is validated by comparing the numerically and experimentally obtained temperatures of the two designed belt drives at different locations within the system. Finally, pulley temperature distribution results are validated by obtaining the experimental measurements on different locations of the pulleys.

Verification of S n (ω pn )
Obtaining an accurate curve for S n (ω pn ) is critical during the calculation of the thermal model. According to the definition of λ pn , the numerical and measured results of the heat dissipation flux of pulleys must be compared under the identical the temperatures at the contact area T pcn and rotating speed ω pn . However, no device can measure the heat dissipation of the pulley at its high rotational speeds. Hence, it is reasonable to compare the calculated and measured results of the temperature distribution of the pulley, because the temperature drop along the radial direction is strongly linked to the heat dissipation rate on pulley surfaces. Figure 13 presents that the temperature distributions of Layer 2 for the pulleys from the numerical model approach and compare this with the experimental temperature data at five different points along the radial direction under the condition of T pcn = 55 • C and rotating speed ω pn = 4000 RPM. In this diagram, the calculated results have a good agreement with the measured temperatures, indicating that the heat dissipation for a pulley is reliable. Moreover, to ensure that the heat dissipation of the pulley is reasonable, the simulation result of the airflow around the pulley is examined. Because the heat dissipation coefficient hR pR is strongly related to the local airspeed, realizing an accurate airflow simulation is critical to calculate the heat dissipation of the pulley. Figure 14 shows the numerical results of the air flows surrounding the pulley during a rotation speed of 4000 RPM. The adopted SST k-ω model uses a blending function to integrate the Wilcox k-ω and the k-ε models. The k-ω model is ideal for the flow simulation in the viscous sub-layer, and the graph exhibits the high velocity air flow near the pulley surfaces. The k-ε models is well suited for free flow in regions away from the wall, and the graph exhibits the vertex flow above the surfaces at the pulley outer radius. Moreover, this turbulent model used in this simulation is time-efficient. In particular, the CPU time for each parametric study is approximately 2-3 h on a desktop with a 12-core 3.4 GHz CPU and 64 GB memory. Therefore, the turbulent model selected in the numerical method can ensure sufficient accuracy of the turbulent flow. The curve for SR nR(ωR pnR ) for different rotation speeds is eligible for the subsequent thermal model calculation. Moreover, to ensure that the heat dissipation of the pulley is reasonable, the simulation result of the airflow around the pulley is examined. Because the heat dissipation coefficient h p is strongly related to the local airspeed, realizing an accurate airflow simulation is critical to calculate the heat dissipation of the pulley. Figure 14 shows the numerical results of the air flows surrounding the pulley during a rotation speed of 4000 RPM. The adopted SST k-ω model uses a blending function to integrate the Wilcox k-ω and the k-ε models. The k-ω model is ideal for the flow simulation in the viscous sub-layer, and the graph exhibits the high velocity air flow near the pulley surfaces. The k-ε models is well suited for free flow in regions away from the wall, and the graph exhibits the vertex flow above the surfaces at the pulley outer radius. Moreover, this turbulent model used in this simulation is time-efficient. In particular, the CPU time for each parametric study is approximately 2-3 h on a desktop with a 12-core 3.4 GHz CPU and 64 GB memory. Therefore, the turbulent model selected in the numerical method can ensure sufficient accuracy of the turbulent flow. The curve for S n (ω pn ) for different rotation speeds is eligible for the subsequent thermal model calculation. Figure 15 shows the curves of S n (ω pn ) for three different pulleys under different speeds. This figure also provides some guidelines regarding the temperature reduction and reliability improvement of the pulley. The pulley with a diameter of 108 mm has a higher λ p than pulleys with a diameter of 76.2 mm. The larger diameter of the pulley can lead to an increase in λ p because of the larger heat dissipation area. The steel pulley has a considerably higher λ p than that of the FRP pulleys. In other words, the high thermal conductivity of the pulley material can increase λ p efficiently because the pulley becomes conductive in terms of the thermal flow, thereby reducing the thermal concentration. models is well suited for free flow in regions away from the wall, and the graph exhibits the vertex flow above the surfaces at the pulley outer radius. Moreover, this turbulent model used in this simulation is time-efficient. In particular, the CPU time for each parametric study is approximately 2-3 h on a desktop with a 12-core 3.4 GHz CPU and 64 GB memory. Therefore, the turbulent model selected in the numerical method can ensure sufficient accuracy of the turbulent flow. The curve for SR nR(ωR pnR ) for different rotation speeds is eligible for the subsequent thermal model calculation.  Figure 15 shows the curves of SR nR(ωR pnR ) for three different pulleys under different speeds. This figure also provides some guidelines regarding the temperature reduction and reliability improvement of the pulley. The pulley with a diameter of 108 mm has a higher λR pR than pulleys with a diameter of 76.2 mm. The larger diameter of the pulley can lead to an increase in λR pR because of the larger heat dissipation area. The steel pulley has a considerably higher λR pR than that of the FRP pulleys. In other words, the high thermal conductivity of the pulley material can increase λR pR efficiently because the pulley becomes conductive in terms of the thermal flow, thereby reducing the thermal concentration.

Model Verification
The thermal model is validated by comparing the temperatures of the belt-side and pulley-side at the belt-pulley surface temperature with the experimental values. These values are critical to determine the thermal failures of materials, because belt-side contact surface temperature is the belt uniformed temperature and pulley-side temperature is the highest temperature in the pulley. The thermal model calculates the heat flux on each component and the temperatures on the contact surfaces, and it can provide the changes in the temperature distribution of the contact surfaces from 2000 RPM to 6000 RPM in intervals of 2000 RPM. Figure 16a,b shows the comparison of the measured and calculated temperatures for a two-belt drive system under various conditions. The model predicts the upward trends of the temperatures when the rotation speed of the belt drive system increases. The diagrams show that the calculated results are in good agreement with the experimental data. The maximum difference is less than 20% of total temperature rise, which is acceptable for industrial use or the subsequent thermal fatigue analysis.

Model Verification
The thermal model is validated by comparing the temperatures of the belt-side and pulley-side at the belt-pulley surface temperature with the experimental values. These values are critical to determine the thermal failures of materials, because belt-side contact surface temperature is the belt uniformed temperature and pulley-side temperature is the highest temperature in the pulley. The thermal model calculates the heat flux on each component and the temperatures on the contact surfaces, and it can provide the changes in the temperature distribution of the contact surfaces from 2000 RPM to 6000 RPM in intervals of 2000 RPM. Figure 16a,b shows the comparison of the measured and calculated temperatures for a two-belt drive system under various conditions. The model predicts the upward trends of the temperatures when the rotation speed of the belt drive system increases. The diagrams show that the calculated results are in good agreement with the experimental data. The maximum difference is less than 20% of total temperature rise, which is acceptable for industrial use or the subsequent thermal fatigue analysis.

Temperature Plot Verification
Another target output of this thermal model is the temperature distribution of a pulley under a certain load/speed condition. Layer2 has the greatest radial distance, as shown in Figure 4, this comparison focusses on the Layer2 along the outer surface of the pulley. Figure 17a,b shows the calculated temperature distributions on 76.2 mm and 108 mm pulleys, and five experimental temperatures at different points along the radial direction of each pulley are used to validate the calculated result. The calculated temperature plots are in good agreement with the measurement data, thereby demonstrating that the model can successfully predict the temperature decrease from the inner to the outer radial distances on Layer2. The maximum temperature difference is within 4 • C in these two plots. Therefore, the thermal model can successfully provide the reliable temperature distribution and heat dissipation for each pulley.
the belt-pulley surface temperature with the experimental values. These values are critical to determine the thermal failures of materials, because belt-side contact surface temperature is the belt uniformed temperature and pulley-side temperature is the highest temperature in the pulley. The thermal model calculates the heat flux on each component and the temperatures on the contact surfaces, and it can provide the changes in the temperature distribution of the contact surfaces from 2000 RPM to 6000 RPM in intervals of 2000 RPM. Figure 16a,b shows the comparison of the measured and calculated temperatures for a two-belt drive system under various conditions. The model predicts the upward trends of the temperatures when the rotation speed of the belt drive system increases. The diagrams show that the calculated results are in good agreement with the experimental data. The maximum difference is less than 20% of total temperature rise, which is acceptable for industrial use or the subsequent thermal fatigue analysis.  Appl. Sci. 2020, 10, x FOR PEER REVIEW 21 of 23

Temperature Plot Verification
Another target output of this thermal model is the temperature distribution of a pulley under a certain load/speed condition. Layer2 has the greatest radial distance, as shown in Figure 4, this comparison focusses on the Layer2 along the outer surface of the pulley. Figure 17a,b shows the calculated temperature distributions on 76.2 mm and 108 mm pulleys, and five experimental temperatures at different points along the radial direction of each pulley are used to validate the calculated result. The calculated temperature plots are in good agreement with the measurement data, thereby demonstrating that the model can successfully predict the temperature decrease from the inner to the outer radial distances on Layer2. The maximum temperature difference is within 4 °C in these two plots. Therefore, the thermal model can successfully provide the reliable temperature distribution and heat dissipation for each pulley. In the above comparisons, some differences can be noted between the calculated and experimental values. One reason for these differences is the measurement error from the devices, which cannot be eliminated completely even if the measurement data is considered as the average values of test results obtained three times for each condition. The second reason is that the turbulence modeling is slightly different compared to the real airflow, especially in terms of the airflow near the walls of the pulley. This aspect represents a balance between the accuracy and time consumption of the numerical method. Nevertheless, the above figures generally show good agreement between the predicted and experimental values.
Overall, all the predictions provided by the thermal model require considerably less computational time and resources to be obtained than those required by the existing approaches. The pulley temperature prediction on this system requires less than 5 s on a workstation with a 12-core 3.4 GHz CPU, except for the one-time computational cost of the baseline numerical simulations. In contrast, the current numerical method requires approximately 17 h [29] to realize the temperature prediction for a belt drive system under one working condition. The comparison thus demonstrates that the proposed thermal model can provide reliable and time-efficient temperature prediction. In the above comparisons, some differences can be noted between the calculated and experimental values. One reason for these differences is the measurement error from the devices, which cannot be eliminated completely even if the measurement data is considered as the average values of test results obtained three times for each condition. The second reason is that the turbulence modeling is slightly different compared to the real airflow, especially in terms of the airflow near the walls of the pulley. This aspect represents a balance between the accuracy and time consumption of the numerical method. Nevertheless, the above figures generally show good agreement between the predicted and experimental values.
Overall, all the predictions provided by the thermal model require considerably less computational time and resources to be obtained than those required by the existing approaches. The pulley temperature prediction on this system requires less than 5 s on a workstation with a 12-core 3.4 GHz CPU, except for the one-time computational cost of the baseline numerical simulations. In contrast, the current numerical method requires approximately 17 h [29] to realize the temperature prediction for a belt drive system under one working condition. The comparison thus demonstrates that the proposed thermal model can provide reliable and time-efficient temperature prediction.

Conclusions
This paper presents a novel thermal model to enable the temperature distribution prediction of a serpentine belt drive system in an accurate and efficient manner. The model uses an analytical-numerical method, which can promptly provide precise temperature results pertaining to pulley irregular structures under abundant belt drive load cases. In this study, a general set of equations was developed according to the flow behaviors within the belt drive system. In addition, the coefficient λ pn and function S n (ωR pn ) were investigated by performing numerical simulations. By solving this set of equations, the belt temperature and thermal flux flowing into each pulley was obtained. Then, utilizing this information, together with the baseline results obtained from the numerical simulations in the previous step, the model outputs the temperature distribution of each pulley on the designed operating conditions. To validate the proposed model, experiments were performed under various operating conditions, and the experimental and calculated results were compared. The good agreements noted in these comparisons prove the accuracy and reliability of the model. In summary, the proposed model can successfully predict the temperature distribution of a belt drive system under several operating conditions in only a few seconds, not accounting for the time taken to establish the database of the baseline numerical simulations in advance. Furthermore, the model provides considerable avenues for further development of the belt drive system. The time-efficient and accurate nature of this model can avoid the thermal failure of the belt drive system in the early design stages. In addition, the model can be applied to the thermal and fatigue analyses of the materials of pulleys and belts. Moreover, because this model involves a numerical method, it can be applied to any type of pulley geometry and belt transmission.