Thermal Stress Coupling Analysis of Ventilated Disc Brake Based on Moving Heat Source

)ere is a thermomechanical coupling phenomenon in the braking process of vehicle disc brakes. Since the solution time of the direct coupling method is long and it is hard to converge, it is not suitable for analyzing heavy-load braking conditions with long braking time. Meanwhile, the sequential coupling method based on fixed heat source cannot achieve the rotation of friction heat source of the brake disc or accurately simulate the change of temperature.)erefore, Abaqus user subroutine is invoked to achieve the heat source rotation. )e results of the simulated temperature field are compared with the results of the complete thermodynamic coupling method, the sequential coupling method based on fixed heat source and dyno test. )e comparison shows that the new method combines the advantages of the first two methods and has a better engineering value.


Introduction
e uneven temperature field causes the thermoelastic deformation of the brake disc and the pad that changes the contact state and contact pressure between brake disc and pad, which greatly affects the important brake characteristics like vibration, noise, and friction [1][2][3].Friction-induced vibrations are known to affect many engineering applications, such as turbine blade joints, wheel/rail of mass transit systems, machine tool/work piece systems, and brake systems [4].Four mechanisms of friction-induced vibration in disc-brake systems have been proposed in the literature [5][6][7][8][9][10][11] as follows: (1) Modal coupling of the structure involving sliding parts (2) Stick-slip motion often seen as the main reason for friction squeal at low sound pressure level (3) Sprag-slip phenomenon proposed by Spurr (4) Negative friction-velocity slope Stribeck [5] found the dependence of the friction coefficient on the relative velocity.Within a certain range, velocity decreases and friction increases in turn, which is called the Stribeck model.Papangelo et al. [12,13] studied a chain of friction-excited oscillators with the nearestneighbour elastic coupling where friction was modelled with an exponentially decaying friction law.It is shown that, in a certain range of driving velocities, multiple stable spatially localized solutions exist, where the vibration localizes on few oscillators.e bifurcation diagram suggests an underlying snaking pattern, as has now been observed in many other fields of physics.
In general, the brake system belongs to weak coupling (thermomechanical coupling) and will experience frictioninduced vibrations.
erefore, the friction coefficient depends on velocity, which has the characteristics of a negative slope in a certain range of velocities.It is noteworthy that, in practical application, it may trigger unstable vibrations and lead to multiple solutions.As a result, thermomechanical coupling analysis becomes an indispensable part in brake design and development [14,15].
Numerous studies have been devoted to the thermal stress coupling analysis.Floquet and Dubourg [16] used the nonlinear finite element program to simulate the brake disc and the friction plate.ey analyzed the temperature and the equivalent stress distribution characteristics of the friction coupling.Lee et al. [17] introduced a finite element method to analyze the temperature field and stress field of a pressure plate.e results show that the stress concentration has great relationship with design parameters, especially the shapes of pressure plates and hub plates.Zagrodzki [18] thought that thermoelastic instability (TEI) can result in "hot spots" and has a great influence on the brake system.erefore, he proposed a novel method to simulate the thermoelastic instability process in the dry friction system.And the finite element spatial discretization method and modal aggregation principle were used.Yevtushenko and Kuciej [19] established a two-element model of tribosystem during braking.e dependences of temperature and stress on the boundary conditions of the ceramic-metal strip were studied.e results show that the obtained results can be applied to the model of a thermal cracking of the frictional elements in the braking stage.
Frictional heat generation is a common phenomenon during the process of friction brake, which leads to high temperature and variable material properties.Li and Barber [20] proposed solutions of the eigenvalue of thermoelastic instability.ey used a suitable interpolation procedure to achieve the simulation of transient thermoelastic contact problems.Burton et al. [21] also proposed a friction contact model, which is based on the nonuniform pressure distribution.e thermoelastic stability conditions between two sliding surfaces were discussed.At the same time, the viewpoint of critical velocity was presented; that is, the thermoelastic instability phenomenon only occurs when the relative velocity of two sliding surfaces exceed the critical velocity.After that, Lee from the University of Michigan and Barber [22] perfected Burton's model and established a sliding analysis model of a finite thickness layer in two half planes for automobile brake analysis, which is now widely used.
e research methods for the coupling problem are mainly theoretical and experimental studies.eoretical research is divided into analytical method and finite element method.With the development of finite element technology, the main way of studying the thermal structure of a disc brake is to simulate the temperature field and stress field of the friction pairs.Choi and Lee [23,24] combined the force equilibrium equation and heat conduction equation to solve the problem and analyzed the transient thermoelastic instability of the disc brake.However, the heat flow boundary condition was considered as axisymmetric and the surface heat source was not mobile, so this model is called the twodimensional axisymmetric model.It is difficult to show the transient situation of heat conduction and the deformation of the disc brake.
According to Avile et al. [25], the warpage of brake discs is the main reason for producing "hot spots" and causing thermoelastic instability, compared with experimental results.Altuzarra et al. [26] pointed out that the variation of the friction coefficient and thermoelastic instability caused the fluctuation of normal and tangential braking forces, which led to the low-frequency brake judder.
From previous studies, it can be concluded that temperature-stress full coupling and heat-stress sequential coupling are the main methods employed to conduct simulation.Actually, the contact area of the friction coupling will change continuously during braking and form a rotating heat source.e immobile heat source is applied to the disc surface in the shape of a complete circle with the traditional coupling method based on the fixed heat source, which disaccords with the actual braking conditions.erefore, in order to make up for the shortcomings of the above analysis methods, it is proposed that a rotating heat source is applied to the disc surface to simulate the friction heat source generated by the actual situation.e simulation is carried out in a typical working condition to verify its effectiveness.Furthermore, the comparison with the experimental results shows that the temperature error is within the allowable range, which further proves that the method is feasible.

Model and Theory
2.1.3D Model.Figure 1 shows the 3D model assembly drawing of the ventilated disc brake.e basic dimensions of the friction coupling are shown in Table 1.
e material of the brake disc is HT250, whose thermophysical properties [27] are shown in Table 2. e density is 7.28 × 10 −9 (t/mm 3 ) and Poisson's ratio is 0.3.
e material of the brake pad is resin matrix composite, whose thermophysical properties [28] are shown in Table 3.
During the braking stage, the heat generated by the friction pair increases the temperature of the contact area.e friction coefficient is not fixed, which is affected by the temperature.
erefore, one should take the change in friction coefficient into consideration for simulation (Table 5).

Mesh.
e hexahedral element has the advantages of fewer numbers, faster solution speed, and higher precision than the tetrahedral one.e hexahedral element is chosen for the brake components, while the tetrahedral element is used for the rigid surface.Given that thermostress coupling analysis is necessary, the tetrahedral unit type is set to C3D4, and the hexahedral unit type is set to C3D8RT.At the same time, the appropriate element size should be selected to reduce the calculation cost under the premise of ensuring the calculation precision.e mesh models of the brake disc and the pad are shown in Figure 2.

Heat Partition Coefficient.
e heat generated by the brake is allocated to the brake disc and the friction plate based on the heat partition coefficients.Figure 3 shows the one-dimensional heat conduction contact model of the friction pair.
According to references [30,31], heat partition coefficient is 2 Advances in Materials Science and Engineering where c denotes the heat partition coefficient, q denotes the heat flux density, k denotes the heat conduction coefficient, c denotes the specific heat, and ρ denotes the density.

Advances in Materials Science and Engineering
One can achieve the following with the above parameters: e ratio that the heat ux absorbed by the friction plate accounts for in the total heat ux is It can be seen that the proportion of heat absorbed by the friction plate is low, which is determined by the material properties of the friction pair, but that by the brake disc is up to 0.9.

Convective Heat Transfer Coe cient.
Figure 4 shows the diagram of convection heat transfer between brake disc and air.
e heat radiation of the brake disc surface is di erent from that of the heat rib.e convective heat transfer coe cient is di erent [32].e convective heat transfer coe cient of the disc surface is [33] where k a is the air thermal conductivity, D is the outer diameter of the disc, and Re is the Reynolds number denoting the dimensionless parameter of uid ow (Equation ( 5)).When Re ≤ 2.4 × 10 5 , the air ow is laminar ow; when Re > 2.4 × 10 5 , the air ow turns into turbulent ow: Re where ω is the brake disc angular velocity, R is the wheel rolling radius, ρ a is the air density, d 0 is the characteristic length, namely, the outer diameter, and μ a is the air viscosity.
From the related literature, k a 0.0276 W/(m• °C), ρ a 1.13 kg/m 3 , and μ a 1.91 × 10 5 kg/(m•s): Figure 5 shows the curve of convective heat transfer and angular velocity.

Disc Pad Pad
q p q d q d q p Substituting equation ( 7) into ( 6), the expression of the convective heat transfer coe cient and the braking time can be achieved: Only the tting convective heat transfer coe cient can be set in ABAQUS, so one must get a tting curve (Figure 6).

Thermal Stress Coupling Analysis Based on Moving Heat Source
In the analysis of heat-stress sequential coupled problems, it is considered that temperature eld is applied in the form of heat source in ABAQUS rather than directly by friction.erefore, only a single brake disc is selected as the research object while other components are without consideration.
e material property and mesh size of the brake disc are the same as those of temperature-stress full coupled methods.At the same time, the type of the brake disc mesh unit is changed to DCC3D8, and the analysis step type is changed to the heat transfer step.
During thermal analysis, it is necessary to impose the heat ow load on the disc surface by a certain shape.It is known from the contact area of the friction pair that the input shape of the friction heat ow is inside the diameter of 98 mm, outside the diameter of 146 mm, and wrapping angle of 66 °, as shown in Figure 7.
In order to realize the movement of heat source, the mass ow rate of the structure should be set in ABAQUS so that the heat source can be rotated by heat conduction analysis.However, the mass ow rate does not support CAE, which cannot be directly set up in ABAQUS interface.It is necessary to write user subroutine in Visual Studio.For this research object, all nodes of the brake disc are selected as the region of the mass ow rate and set to an aggregation.e inp le can be output after setting up the correct material, analysis step, boundary condition, load, and so on.After the analysis step, the key word for adding the mass ow rate is handwritten.e speci c calculation formula of the mass ow rate should be based on di erent braking conditions in Visual Studio to write di erent subroutines.Finally, these subroutines can be called when the calculation is submitted.
Figure 8 presents that thermal boundary conditions of the model are mainly convective heat transfer and heat conduction.e initial temperature of the brake disc model is de ned as 20 °C.
en thermal stress analysis is carried out.At this time, analysis step type should be changed to coupled temp-displacement type, and mesh unit type is changed to C3D8RT.
e temperature eld model obtained by the thermal analysis model is imported into the stress analysis model.For the boundary condition, the mechanical stress of the disc can be ignored when compared with thermal stress, so the mechanical stress is not considered.At the same time, six degrees of freedom of the disc must be restrained.
Generally, when adopting the temperature-stress full coupled method, heat is generated due to the frictional energy loss [34].According to references [34,35], the expression is as follows: q(t, r) ημpv(r, t), (9) where q(t, r) denotes the heat ux, η denotes the proportion of heat absorbed by the disc, μ denotes the friction   Advances in Materials Science and Engineering coefficient, p denotes the brake pressure, and v(r, t) denotes the velocity of a certain point on the disc.
As we have discussed in Introduction, the friction coefficient also depends on velocity.But equation ( 9) is more common about the friction energy generated.More importantly, once by applying the heat-stress sequential method, friction heat is directly exerted to the surface of the brake disc.

Model Accuracy Verification
e heat-stress sequential method based on moving heat source can simulate the actual heat source rotation and save calculation cost, which is based on the sequential method with immobile heat source.In order to verify its accuracy and feasibility, the heat-stress sequential coupling method is adopted in this section to simulate the temperature rise and stress changes of the brake disc under emergency braking conditions.e results will be compared with the dyno test results.

ermal Load in Braking Process.
e energy conversion method is used to calculate frictional heat generated during the braking process [36].e expression is as follows: where s is the slip rate, u a is the speed of the car under the condition of no wind, f is the rolling resistance coefficient, C D is the air resistance coefficient, A is the windward area, ρ a is the air density, and i is the slope.e distribution of frictional heat is calculated according to the braking force distribution coefficient β ′ .Because the research object is the front wheel brake and the heat generated by the brake is assigned to the brake disc and the friction plate according to the heat flow distribution coefficient, the heat flux on the single model working face is as follows: where A f represents the input area of friction heat flow, 0.00675 m 2 .

Emergency Braking Condition.
e initial speed of the vehicle under the emergency braking condition is 80 km/h.e braking deceleration is 5.8 m/s 2 .e slope is out of consideration because of the emergency braking condition.
e friction input heat flux density is calculated as follows: Figure 9 shows the linear fitting relationship curve between input heat flux and braking time, which is q in (t) � 5.923 × 10 6 − 1.546 × 10 6 t. ( Meanwhile, the subprogram needs to determine the change of the mass flow rate caused by the change of disc speed and along the radial brake disc.e expression of the mass flow rate at the initial braking time is as follows: where q v , equal to ω • h, is the volume flow rate through unit section area at the braking initial time.e mass flow rate of some point at any time is where t is the brake time.With the relevant data into the expression, the mass flow rate is en, the subroutine can be called in the job module.It is proved that during the whole braking process, the heat source is rotated on the brake disc surface.e rotating heat flux diagrams are depicted in Figure 10.
Figure 11 presents the distribution nephogram at the highest temperature of the brake disc.
e temperature rise curve of the maximum temperature node is shown as Figure 12, which presents the trend of rising and zigzag.
Compared with thermostress full coupling and thermostress sequential coupling based on the fixed heat source, thermostress sequential coupling based on rotated heat source shows the change in temperature more accurately (Figure 13).
Likewise, the stress distribution is basically the same as that of full coupling methods.e stress change curve of the maximum stress node is shown in Figure 14.rough comparison, it can be seen that the conclusion is the same as that of the temperature rise curve (Figure 15).
It can be seen from Table 6 that when adopting the same computer configuration for simulation calculation, the result errors between these methods are very slow but the calculation time is greatly reduced.

Brake Dynamometer Validation.
A brake dynamometer test is conducted to validate the model.e configuration of the brake dynamometer test is shown in Figure 16.e brake  Advances in Materials Science and Engineering disc is connected to the brake dynamometer rotating part.A thermocouple sensor is connected to the disc outboard surface and the pad friction surface.e disc is equipped with brake caliper and pads (Figure 17).e performance test is carried out using the QC/T 564-2008 Performance Requirements and Bench Test Methods for Passenger Car Brakes.e braking pressure is 3.56 MPa, and the temperature is 20 °C.
e analysis temperatures of the brake disc are compared with the experimental temperatures.e temperature results of the disc in the brake dynamometer test and the coupling analysis are shown in Figure 18. e maximum error is less than 10%, and the maximum temperature error is only 1.7%.It is concluded that the thermostress sequential coupling method based on rotating heat source is feasible.Advances in Materials Science and Engineering

Thermostress Analysis under Repeated Braking Condition
5.1.ermal Boundary Condition.It is not easy to realize the simulation by thermostress full coupling under repeated braking conditions that need long periods.However, the new method can not only signi cantly reduce the computation load but also satisfy the calculation precision.e test requirements are shown in Table 7 according to GB12676-1999.e whole braking process is divided into 15 cycles where the speed includes three stages: brake, acceleration, and steady speed.e input heat ux should be segmented.
ere is no heat generated during acceleration and steady speed stages, so the heat ux is zero.As a result, the heat ux of the brake stage is e corresponding heat ux density curve is shown in Figure 19.
Meanwhile, the mass ow rate should be rewritten when the subroutine is invoked.e expression of the mass ow rate for one cycle is as follows:

Temperature Field.
e temperature curve of the maximum temperature node under the cycle braking condition is shown in Figure 20.As shown in the picture, node temperature is uctuating upward until six or seven cycles where the maximum temperature can reach 172.9 °C.In the process of braking, the temperature rises rapidly.In the acceleration stage, the disc temperature will gradually decrease because of no input heat ux.Since the whole braking process is long and the disc has enough time to dissipate heat, it does not generate high temperature.

Stress Field.
With the same node as the temperature eld analysis, the stress curve under repeated braking conditions is shown in Figure 21.e variation curves of thermal stress are basically the same as those of temperature, which uctuate upward and tend to approach dynamic equilibrium gradually.

Conclusion
e solution time of the full coupling method is longer than that of the sequential coupling method, and it is hard to converge.Meanwhile, the traditional coupling method based 8 Advances in Materials Science and Engineering on the fixed heat source cannot accurately simulate the actual rotation of the heat source.erefore, the thermostress sequential coupling based on rotating heat source is proposed, which makes the heat source rotate by invoking the subroutine that includes the mass flow rate of the disc.e following conclusions can be drawn from this paper: (1) When in the emergency braking condition, not only it is easy for the new method to converge but also the    Advances in Materials Science and Engineering simulation results are consistent with the results from tests or full coupling method, thus verifying the reliability and practicability of the proposed method.(2) According to the simulation results of repeated braking conditions, it can be concluded that the disc surface temperature fluctuates periodically and tends to be stable, which has the same trend as thermal stress.

Figure 1 :
Figure 1: 3D model assembly drawing of the ventilated disc brake.

Figure 3 :Figure 4 :Figure 2 :
Figure 3: Heat conduction contact model of the friction pair.

Figure 6 :
Figure 6: e relation curve of the convective heat transfer coe cient and the braking time on the surface of the brake disc.

Figure 5 :
Figure 5: Change curve of the convective heat transfer coe cient with angular velocity.

Figure 10 :
Figure 10: Heat ux rotated on the disc brake.

Figure 11 :Figure 9 :Figure 12 :
Figure 11: Distribution nephogram at the highest temperature of the brake disc.

Figure 15 :Figure 14 :Figure 13 :
Figure 15: Comparison of three simulation methods for the stress change curve.

Figure 16 :
Figure 16: Configuration of the brake dynamometer test.

Figure 19 :Figure 20 :
Figure 19: e heat flux of a cycle braking period.

Figure 21 :
Figure 21: ermal stress change curve of disc surface in repeated braking.

Table 3 :
ermophysical properties of the brake pad.

Table 4 :
ermophysical properties of the brake backplate

Table 6 :
Comparison of simulation results and calculation time of three simulation methods.Full couplingSequential coupling based on rotated heat source Sequential coupling based on fixed heat source T max ( °C)

Table 7 :
Type M1 vehicle type I test requirements.