Thermoelastic vibration analysis of disc brake based on multi-physical field coupling method

2.1. Theoretical analysis of coupling process

The thermodynamics characteristic of brake disc caused by friction is one kind of non-equilibrium problem [10]. The phenomenological equation can be expressed in space $V$ as follows:

where $\rho z$ is density of extensive quantity $Z$, $J_z$ is density of for extended flow, ${\sigma }_Z$ is generation rate of extensive quantity, $\tau$ is time. Assume that there are $i$ kinds of forces in the friction system. According to Newton’s second law and the symmetry of the stress tensor, the momentum equation of the brake disc can be deduced as follows:

where $U$ is the velocity vector of center of mass of brake disc, $\rho$ is density of the disc, $P$ is stress tensor, $F_i$ is the volume force of each function component.

Ignoring chemical reaction factor between the brake disc and brake pad, the entropy balance equation of brake disc can be deduced by the momentum equation, mass and energy conservation laws as followed:

where $s$ is the entropy density, $J_q$ is density of heat flow rate, $J_k$ and ${\mu }_k$ are respectively diffusion flow and chemical potential of the material, $T$ is temperature, $M$ is the molar mass, II is interfacial stress tensor. The entropy balance equation can fully express the fluid-thermal-structural coupling state of the brake disc.

According to the energy conservation equation, the three-dimensional transient heat transfer control equation of the brake disc can be obtained as follows:

where $a$ is the thermal diffusivity of brake disc, which is related to material quality, unit of m²/s. $T$ is temperature, unit of K. $\dot{\mathrm{\Phi }}$ is internal heat source, unit of J/m³.

Convective heat transfer is a ubiquitous physical phenomenon in nature, which is more complicated than pure heat conduction. Therefore, the static convective heat transfer coefficient is generally calculated in engineering, but the accuracy is limited. In order to obtained more accurate results, the dynamic change of convection heat transfer coefficient in space is considered in the paper.

The contact surface between the brake pad and brake disc is discrete. Therefore, the actual contact area consists of all the contact points and determines the size of the friction heat source. According to Uppal and Probert’s research, it can be known that the contact points will remain or even decrease along with the increase of the actual contact area if the brake pressure increases to the critical load. Thus, it can be seen that the braking torque will not always increase when the brake pressure is increasing.

In the braking process, the contact of friction pair has three statuses, including elastic, plastic and middle (elastic-plastic) forms, and the increase rate of the actual contact area in the elastic state is smaller than that of the plastic state [11]. According to the contact point of the spherical micro convex body theory, when $P_0=$2.76${\tau }_y$ (where $P_0$ is actual contact area center pressure, ${\tau }_y$ is ultimate tangential stress), plastic deformation will occur at the contact surface. The condition of the elastic-plastic deformation on the surface of the brake disc can be expressed by Tabor’s depth of press calculation method as follows:

where $a_c$ is the contact spot radius, $r$ is spherical micro convex body radius, ${\sigma }_T$ denotes the average shear stress, $k_n$ is a scaling factor (determined by ratio of the average shear stress and the effective yield stress).

The influence of surface profile on the elastic-plastic state can be expressed by Williamson and Greenwood index method as shown in Eq. (6):

where $\sigma$ is the standard deviation of the profile, $E\text{'}$ is the average elastic modulus of brake pad and brake disc, $H$ is hardness of the brake disc, ${\mu }_1$ and $E_1$, ${\mu }_2$ and $E_2$are respectively Poisson’s ratio and elastic modulus of brake disc and brake pad. In unit load, when $\varphi <\text{0.6}$, elastic state will remain, and when $\varphi >1$, the plastic state will occur.

In the frictional thermal analysis of disc brake, many scholars assume that the heat flux on the two contact surfaces is equally distributed, which is unreasonable in practice. The frictional heat flux on the interface caused by sliding is related to the thermal resistance of the two friction surfaces. At present, it is still quite difficult to quantitatively determine the thermal resistance between friction interfaces. According to the calculation principle of friction and wear, the reasonable thermal conduction boundary conditions between the friction pairs can be approximated as follows. In the real contact area, the temperature of each contact point pair composed of the contact point of the brake pad and the brake pad must be equal, but outside the real contact area, the surface temperature is not necessarily equal, that is, the temperature of the contact point pair is not necessarily equal:

where $T_1$ and $T_2$ is the characteristic temperature of contact point between the brake pad and brake disc, unit of K. $q_1$ and $q_2$ is heat flux density of contact surfaces between brake disc and brake disc , unit of J/(m²·s). $q$ is the total heat flux, unit of J/(m²·s).

The heat flux distribution coefficient of brake pad is:

where ${\lambda }_1$ and ${\lambda }_2$ are the thermal conductivity of the brake pad and the brake disc, unit of W/(m·K). $c_1$ and $c_2$ are the specific heat of the brake pad and the brake disc, unit of J/(kg·K). ${\rho }_1$ and ${\rho }_2$ are the density of the brake pad and the brake disc, unit of kg/m³.

According to Eq. (8), the brake pads of different materials have important influence on heat flow distribution. For example, the heat flow distribution coefficient of semi-metallic materials is 9.8, and the heat flow accounts for only 9.3 %.

According to Hooke’s law, the total strain component of the thermoelastic element in Cartesian coordinate system $Oxyz$ can be obtained as follows:

where $\alpha$ is the thermal expansion coefficient of brake disc, ${\epsilon }_x$, ${\epsilon }_y$, ${\epsilon }_z$ is positive strain in different direction, ${\gamma }_{yz}$, ${\gamma }_{zx}$, ${\gamma }_{xy}$ is shear strain, ${\sigma }_x$, ${\sigma }_y$, ${\sigma }_z$ is positive stress, ${\tau }_{yz}$, ${\tau }_{zx}$, ${\tau }_{xy}$ is shear stress, $\mu$ is Poisson ratio; $E$ is modulus of elasticity.

It can be obtained in two-dimensional $Oxy$ plane as followed:

The geometric equation of strain components can be expressed as:

The two-dimensional temperature-displacement coupling differential equation (expression of tight coupling) of brake disc can be derived as followed:

During the heat exchange calculation (loose coupling) between the brake disc and air, the convection heat transfer coefficient (also called surface heat transfer coefficient) is dynamic. Generally, the key physical quantities include air physical properties, flow rate, flow state, coupling surface shape, etc. In the paper, assume air velocity at the coupling point $\left(x,y\right)$ is $\left(u_a,v_a\right)$, temperature is $T_a$. Ignoring internal energy that produced by viscous dissipation, according to the conservation of energy law, the energy equation of convective heat transfer in two-dimensional space can be deduced as follows:

The mass conservation equation can be expressed as:

The momentum conservation equation can be expressed as:

where $T_a$, ${\rho }_a$, $c_p$, ${\lambda }_a$ are respectively temperature, density, specific heat, time, heat conductivity coefficient of air, $p_a$ stands for air pressure, $F_x$ and $F_y$ are respectively the volume force in $x$ direction and y direction.

It can be seen from Eqs. (14-17) that these equations are closed and there are only four unknown boundary parameters. But due to the complexity and nonlinear characteristics of the momentum equation, the mathematical computation cannot be completely obtained [12]. In order to solve the problem, boundary layer conditions are introduced, so that the loose coupling method of temperature field and air flow field calculation is proposed and achieved by node data sharing platform of MPCCI in the paper.

2.2. Coupling implementation scheme

In order to realize the fluid-thermal-structural coupling calculation of the brake disc, a numerical simulation method based on multi physics field coupling standard MPCCI, nonlinear structural analysis software ABAQUS and fluid analysis software FLUENT is presented as shown in Fig. 1. The principle of coupling calculation is that ABAQUS and FLUENT will calculate solid-thermal model and fluid-solid model [12, 13] respectively and simultaneously according to the FEA (Finite Element Analysis) method. The FEA models [14, 15] of solid and fluid (air flow direction is parallel to the disc transverse plane) should be established in the same coordinate system, and the coupling area should be defined in each software as shown in Fig. 2. In the coupling calculation process, ABAQUS performs the function of coupling data receiving and provides parameters such as node displacement and wall temperature. FLUENT performs the function of coupling data exchanging, and provides parameters such as the convection temperature and the convective heat transfer coefficient. Meanwhile, MPCCI will make sure the coupling data has synchronization, which means that the whole interpolation calculation step is the same, and the boundary conditions of the coupling region are continuously updated.

During the interpolation calculation, MPCCI will use the bucket pre-contact search algorithm to determine the grid matching attribute between the fluid model and the solid model. According to the physical quantity, the coupling data can be divided into non-conservation data (such as displacement, velocity, temperature, etc.) and conservation data (such as quality, flow, etc.). In the pre-processing of MPCCI, the coupling area and calculation time step will be defined. In order to ensure the data is shared properly, the fluid model is set as “data exchange”, and the solid model is set as “data receiving”.

Fig. 1Principle of coupling calculation

Fig. 2The FEA model of coupling calculation

2.3. Model pretreatment and analysis of ABAQUS

In order to ensure the computational efficiency and convergence, the element partition and local optimization method are used in hexahedral mesh division of both model as shown in Fig. 2. The mesh type of C3D8HT in ABAQUS model is chosen, of which cells number is 17370 and nodes number is 25821. Temperature effects are considered during the material properties setting, so that the results can be more accurate.

In the preprocessing of the model, the main parameters and boundary conditions are set as follows.

(1) The analysis step is established by means of temperature-displacement coupling. Considering the convergence and calculation efficiency, the time step is 0.001 s and the load step is 10000. The surface coupling reference point method is applied to the degree of freedom constraint of brake disc and brake pad. The radiant emissivity is defined as 0.28, the heat flux distribution coefficient is calculated by Eq. (8), and the interface thermal conductivity is defined as followed:

where $h_c$ is interface thermal conductivity between brake disc and brake pad, unit of K/W. $k$ is the minimum value of thermal conductivity of friction pairs, W/(m·K). $\mathrm{t}\mathrm{a}\mathrm{n}\theta$ is mean value of absolute slope of contour. $P$ is brake pressure, unit of Pa. $H$ is microhardness, kg/m². $\sigma$ is contact stress, unit of Pa.

(2) The brake pressure amplitude is set as 3 MPa, the brake force acting time is set as 1 s, and the brake disc speed amplitude is set as 46.4 rad/s. In order to match the test conditions, the initial temperature of the brake disc brake is defined as 100 °C, and set the surface of the brake disc as coupling surface.

2.4. Model pretreatment analysis of FLUENT

The fluid model includes the moving region and the static zone. The moving region is mainly used to simulate the rotation of the brake disk and calculate the flow field around the coupling surface. The static region is mainly used to provide the inlet and outlet of air flow. Considering the working load and calculation efficiency of the fluid model, the sliding grid method [16, 17] is used to realize the wall movement in the paper. In each time step, the flow variables on the interface nodes need to be transferred to realize the coupling solution of the two regions. This method can save a lot of calculation costs compared with the grid reconstruction method, and is conducive to improving the calculation efficiency of multi-physical field coupling.

As a typical near-wall mesh, boundary layer meshes have different grid requirements for different near-wall treatments [18]. In order to guarantee the quality of the boundary layer mesh, the dimensionless sublayer velocity $u^{+}$ and the distance to the wall $y^{+}$ should be defined:

where $u_a$ is air velocity, unit of m/s. ${\tau }_w$ is wall shear stress, Pa. $y_a$is normal distance along the wall, unit of m.

According to Plante’s theory, the relationship between dimensionless parameters $u^{+}$ and $y^{+}$ can be obtained by different test method as shown in Fig. 3. It can be seen that there are mainly three different boundary layer structures as follows.

Laminar sublayer: 0 $<y^{+}<$ 5, $u^{+}=y^{+}$.

Buffer layer: 5 $\le y^{+}\le$30, $u^{+}=5\mathrm{l}\mathrm{n}{y}^{+}-3.25$.

Turbulent core layer: $y^{+}>$30, $u^{+}=2.5\mathrm{l}\mathrm{n}{y}^{+}+5.5$.

According to the Reynolds number formula, it can be known that the braking process is a typical high turbulence model. Under the condition of multi-physical field coupling analysis, the wall function method is more reasonable. The initial conditions of the fluid model and the settings of the solver are shown in Table 2, and the high turbulence RNG $k$-$\epsilon$ model [19] is used.

The non-equilibrium wall function is used to calculate the turbulent flow energy of the near wall controlled body. Sublayer structures under the wall rule can be divided into laminar boundary layer (also known as linear or viscous sublayer) and turbulent boundary layer (turbulent core layer) as shown in Fig. 4. According to the average velocity wall method defined by FLUENT wall function, the expressions of dimensionless sublayer velocity $U^{*}$ and height $y^{*}$ in the first layer grid of node p are given as follows:

where $y_p$ is the distance from point $p$ to the coupling wall, $U_p$ is the average velocity of air at point, $\mu$ is dynamic viscosity coefficient of air, $k_p$ is turbulent kinetic energy at point $p$, $C_{\mu }$ is the coefficient of RNG turbulence model, value of 0.0845.

In the non-equilibrium wall function, the temperature wall rule is still consistent with the standard wall function, but the average velocity of turbulent core layer is more sensitive to the pressure gradient.

Fig. 3The distribution of the boundary layer

Fig. 4Sublayer structure under wall function

2.5. Convergence judgment and analysis of simulation

In this paper, the control of convergence in simulation mainly comes from grid optimization. Under the condition of near-wall function, the mesh generation in boundary layer is very important. The height of the first layer of mesh should be set in the logarithmic rate layer. This is because the laminar flow (linear) criterion is used when $y^{+}<$12.25 in CFD calculation, so the first layer mesh need not be too dense, because the wall function does not work in the viscous bottom layer. For standard wall function and non-equilibrium wall function, the y⁺ value of the first grid node is required to be 11.5-30 $<y^{+}<$ 200-400. According to the checking results, for the non-equilibrium wall function, when the $y^{+}$ value of the first grid node is between 30 and 60, good convergence can be obtained. However, the grid with good quality needs to be adjusted repeatedly. Whether $y^{+}$ meets the requirement in the result of trial calculation must be checked. If not, the height of the first layer of grid should be redefined. The value of $y^{+}$ in the radial direction of the coupled surface mesh is derived as shown in Fig. 5.

Fig. 5The y+ value of the coupling surface in radial direction

It can be seen that the value of $y^{+}$ in the first layer mesh node is 32.5, and the height and density of the mesh meet the requirements. Through the continuous adjustment and optimization, the number of elements in the fluid model is 618114, the number of nodes is 118445, and the maximum value of skewness is 0.85.

The residual curve of numerical simulation model is shown in Fig. 6. It can be seen that in the process of iteration, the residual value of energy is less than 10^-6, the residual value of other variables is less than 10^-3, and the mass, momentum and energy of the whole system are conserved. Therefore, it can be judged that the digital simulation is convergent well.

Fig. 6Residual curve of numerical simulation model

2.6. Experimental verification

The experimental verification for fluid-thermal-structural coupling calculation is realized by bench test method of Link3900 NVH platform. At present, the Link3900 NVH device is one of the most advanced tribology experimental machine for disc brake currently, of which detected parameters includes temperature, pressure, friction coefficient, rotate speed, humidity, torque, moment of inertia, noise, etc. The overall installation diagram of the disc brake and sensor is shown in Fig. 7, the high precision thermocouple is used in temperature detection. The drill hole is inserted into the brake disk and the conducting wire is directly connected to the test machine as shown in Fig. 8.

Fig. 7The overall installation diagram

Fig. 8The thermocouple installation method

The experimental parameters are set up in accordance with the simulation boundary conditions, in other words, the simulation models are built based on experimental conditions, especially friction coefficient [20, 21]. Then the temperature comparison between temperature sensor node and corresponding node of FEA model can be obtained as in Fig. 9. It can be seen that the simulation results of temperature match well with the experimental results, with an average deviation of 3.9 %. As the braking time continues, there is a small phase difference between the test results and the simulation results, which is due to the error accumulation of the constant speed control. All in all, the experiment proved that the coupling analysis scheme of fluid-thermal-structural multi-physical field of disc brake has high accuracy and reliability.

Fig. 9Comparison of simulation results and test results

a) Comparison of 0-1 s temperature changes

b) Comparison of 1-10 s temperature changes

3.1. Path definition

In order to study the variation of the convective heat transfer coefficient in different conditions and directions, two paths (lines) are designed for the structure of the coupling surface, including the path of PS in the radial direction of the end face of the brake disc and the path of PF in the radial direction of the bottom face of the ventilation channel. Path of PS and path of PF are defined in static coordinates, and their directions are perpendicular to the velocity direction of the air and parallel to the $x$-axis of the static coordinate system, as shown in Fig. 10(a). In order to study the characteristics and regularity of temperature and stress field imbalance, different paths (points) are defined for the end face of the brake disk (path of PA) and the side face of the fin (path of PB). The schematic diagram of the path definition is shown in Fig. 10(b).

Fig. 10Path definition in different face

a) Path definition of FLUENT model

b) Path definition of ABAQUS model

3.2. Fluctuation of flow field

At the time of 1.0 s, the velocity along the $Oxy$ plane and the convective heat transfer coefficient $h_w$ along the coupling surface are shown in Fig. 11 and Fig. 12, respectively.

As can be seen in Fig. 11, the air flow field is the result of the interaction of the inlet speed and the brake disk speed, and the streamline direction is obviously affected by the rotation of the brake disk. The wind speed and flow direction of air in different ventilation ducts are quite different, and the wind speed perpendicular to the windward surface is larger.

As can be seen from Fig. 12, the convective heat transfer coefficients at different positions on the surface of the disk vary greatly. The convective heat transfer coefficients near the windward surface are the largest and decrease gradually along the $x$-direction. When the air passes through the ventilation slot, the upwind surface appears again, and the convective heat transfer coefficient shows the same trend. It also can be seen that the greater air speed and pressure will come out greater heat transfer coefficient.

Fig. 11Velocity pathline of profile

Fig. 12Heat transfer coefficient of coupled faces

In order to study the influence of vehicle speed on the heat transfer coefficient of the coupling surface, the flow field characteristics were calculated under three wind speeds of 10 m/s (brake disc speed 23.2 rad/s), 20 m/s (brake disc speed 46.4 rad/s) and 30 m/s (brake disc speed 69.6 rad/s). Then the convective heat transfer coefficient at different speeds of path PS and path PF can be obtained as in Fig. 13. From Fig. 13, it can be seen that the heat transfer coefficients along both path show the same change characteristics under three different wind speeds. Along the $x$ direction, the values increases slowly and then increases rapidly.

Fig. 13Fluctuation of convective heat transfer coefficient in different path

a) Path of PS

b) Path of PF

3.3. Fluctuation of temperature field

The overall transient temperature field nephogram at the specific time of 1.0 s and 5.0 s is shown in Fig. 14. It can be seen that the “hot spot” phenomenon in the temperature field is quite obvious, but it would be gradually disappeared as heat transfer process went on. From the direction of temperature conduction, the temperature of the brake disc in the axial direction is greater than the radial conduction velocity.

The curves of temperature variation on different paths are shown in Fig. 15. From Fig. 15(a), it can be seen that during the braking process (0-1 s), the temperatures of the four nodes (nodes 1-4) in the friction plane all show a serrated upward trend. The temperature at the center of the brake pad (node 2) is the highest because the heat source area of the brake pad is large and the heat transfer rate at the center is relatively low compared with the heat generation rate. From Fig. 15(b), it can be seen that the radial temperature of the central node on the side of the fin always increases with time, but does not show zigzag characteristics. The farther away from the center of the friction surface, the lower the rate of temperature rise, and the more obvious the linear change.

Fig. 14Transient temperature field of disc brake

a) The time of 1.0 s

b) The time of 5.0 s

Fig. 15Fluctuation of transient temperature in different path

a) Path of PA

b) Path of PB

3.4. Fluctuation of stress field

The overall transient stress field nephogram at the specific time of 1.0 s and 5.0 s is shown in Fig. 16. It can be seen that the “hot spot” phenomenon is more obvious, especially in initial braking stage. It is mainly caused by the fin structure, which will result in unbalancing brake pressure on the contact surface. The brake pressure on the overlapping position of the brake pad and fin is larger than both sides, and the “hot spot” phenomenon in stress field does not last long because the elastic deformation is recovered very soon.

Fig. 16Transient stress field of disc brake

a) The time of 1.0 s

b) The time of 5.0 s

The curves of stress variation on different paths are shown in Fig. 17. From Fig. 17(a), it can be seen that the stress change in radial direction (nodes 1-4) of friction plane is consistent with the change of temperature in the braking process. The stress variation of node 5 shows more obvious fluctuation, but it remains relatively small because it is far from the friction surface. From Fig. 17(b), it can be seen that there is a significant difference between the law of stress change and the law of temperature change, especially in the braking stage. The phenomenon of stress mutation is very obvious because the brake pad has a direct compressive effect on the ribs. When the brake pad skips the joint, the stress decreases instantaneously and there is almost no time delay phenomenon.

Fig. 17Fluctuation of transient stress in different path

a) Path of PA

b) Path of PB

3.5. Fluctuation of thickness size and braking torque

Within the friction radius range, the fluctuations of the average brake disc thickness and the braking moment are extracted as shown in Fig. 18 and Fig. 19, respectively. It can be seen that with the increase of braking pressure, the fluctuation value of the brake disc thickness is negative, that is, the brake disc pressure thickness is thinner than the original brake disc, and the thickness will gradually recover when leaving the brake disc, so the reciprocating motion will fluctuate. At the beginning of braking, with the increase of static braking moment caused by pedal force, the braking moment also shows an upward trend. When the static braking moment reaches its maximum and remains constant, the change trend of braking moment and braking pressure tends to be gentle, but then the braking moment shows a slow downward trend. This is because the friction coefficient between the disc and the disc decreases with the increase of temperature.

Fig. 18Fluctuation of brake disc thickness size

Fig. 19Fluctuation of braking torque