*Article* **Investigation into Multiaxial Character of Thermomechanical Fatigue Damage on High-Speed Railway Brake Disc**

**Chun Lu 1,2, Jiliang Mo 1,2,\*, Ruixue Sun 1,2, Yuanke Wu 1,2 and Zhiyong Fan 1,2**


**Abstract:** The multiaxial character of high-speed railway brake disc thermomechanical fatigue damage is studied in this work. Although the amplitudes and distributions of temperature, strain and stress are similar with uniform and rotating loading methods, the multiaxial behavior and out-of-phase failure status can only be revealed by the latter one. With the help of a multiaxial fatigue model, fatigue damage evaluation and fatigue life prediction are implemented, the contribution of a uniaxial fatigue parameter, multiaxial fatigue parameter and out-of-phase failure parameter to the total damage is discussed, and it is found that using the amplitude and distribution of temperature, stress and strain for fatigue evaluation will lead to an underestimation of brake disc thermomechanical fatigue damage. The results indicate that the brake disc thermomechanical fatigue damage belongs to a type of multiaxial fatigue. Using a uniaxial fatigue parameter causes around 14% underestimation of fatigue damage, while employing a multiaxial fatigue parameter without the consideration of out-of-phase failure will lead to an underestimation of about 5%. This work explains the importance of studying the thermomechanical fatigue damage of the brake disc from the perspective of multiaxial fatigue.

**Keywords:** thermomechanical fatigue; multiaxial fatigue; out-of-phase failure; brake disc; life prediction; high-speed railway; damage mechanism

### **1. Introduction**

Fatigue is a common damage type in the engineering field, which is resulted from cyclic load [1]. The characteristics of fatigue damage are different from static damage, i.e., fatigue failure may occur when the cyclic stress is much smaller than the static strength limit. In railway engineering, the cycle of the braking process can cause serious fatigue damage of the brake disc.

Due to the large tonnage and high speed of a vehicle, a tremendous amount of kinetic energy is transformed into thermal energy during the actions of implementing braking, precision stopping and stabling speeding [2]. This transformation arouses significant temperature increasement and exposes the brake disc to severe thermomechanical loading, resulted from the frictional contact between the brake disc and pad [3]. Because of this, the thermomechanical fatigue damage is considerably severe [4,5] (see Figure 1), and it could be more serious in the long slope line. This problem sets a restriction against the further speed increasement, economic saving and security improvement of railway transportation, and even could cause catastrophic accidents. Therefore, the study of thermomechanical fatigue damage on the brake disc is of great importance from both academic and engineering perspectives.

**Citation:** Lu, C.; Mo, J.; Sun, R.; Wu, Y.; Fan, Z. Investigation into Multiaxial Character of Thermomechanical Fatigue Damage on High-Speed Railway Brake Disc. *Vehicles* **2021**, *3*, 287–299. https:// doi.org/10.3390/vehicles3020018

Academic Editor: Ralf Stetter

Received: 27 March 2021 Accepted: 28 May 2021 Published: 1 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Thermomechanical fatigue of brake disc [6],[7]. **Figure 1.** Thermomechanical fatigue of brake disc [6,7].

There are different factors affecting fatigue failure [8],[9] and several fatigue damage types can be detected on the brake disc (Figure 1), i.e., heat checking (crackle), radial crack, circumferential crack and hot spot. Kasem et al. [10] analyze the subsurface damage beneath the hot spots, plastic deformation and microcrack propagation are detected under the rubbing surface. Collignon et al. [11] study the damage of brake discs and it is found that residual hoop tensile stress is relevant to radial microcracks. Honner et al. [12] study the thermoelastic instability of the brake system, which is verified to be relevant for nonuniform temperature distribution. Yang et al. [13] analyze the fatigue mechanism on the brake disc, fatigue cracking is observed in the interior and border of hot spots. Panier et al. [14] propose a classification for hot spots and suggest a scenario of their occurrence by experimental methods. Kim et al. [6] find that the maximum thermal stress under a variable pressure distribution is quite relevant with fatigue cracking. Wu et al. [15] study the peak temperature and crack geometry evolution with the help of an inserted semi elliptical surface crack. Adamowica and Grzes [16] study the brake disc temperature distribution, and it is found that the influence of convection cooling is significant in the brake release period. A review of the thermomechanical performance of brake discs with numerical and experimental approaches can be found in the work of Afzal and Mujeebu [17]. According to previous research, the commonly used method for brake disc thermomechanical fatigue investigation is adopting the distribution and amplitude of tempera-There are different factors affecting fatigue failure [8,9] and several fatigue damage types can be detected on the brake disc (Figure 1), i.e., heat checking (crackle), radial crack, circumferential crack and hot spot. Kasem et al. [10] analyze the subsurface damage beneath the hot spots, plastic deformation and microcrack propagation are detected under the rubbing surface. Collignon et al. [11] study the damage of brake discs and it is found that residual hoop tensile stress is relevant to radial microcracks. Honner et al. [12] study the thermoelastic instability of the brake system, which is verified to be relevant for non-uniform temperature distribution. Yang et al. [13] analyze the fatigue mechanism on the brake disc, fatigue cracking is observed in the interior and border of hot spots. Panier et al. [14] propose a classification for hot spots and suggest a scenario of their occurrence by experimental methods. Kim et al. [6] find that the maximum thermal stress under a variable pressure distribution is quite relevant with fatigue cracking. Wu et al. [15] study the peak temperature and crack geometry evolution with the help of an inserted semi elliptical surface crack. Adamowica and Grzes [16] study the brake disc temperature distribution, and it is found that the influence of convection cooling is significant in the brake release period. A review of the thermomechanical performance of brake discs with numerical and experimental approaches can be found in the work of Afzal and Mujeebu [17].

ture, stress and strain items. However, due to the rotating movement, the thermomechanical fatigue should belong to multiaxial fatigue damage [18], which is complex with numerous influencing factors [19]. Regrettably, the multiaxial feature of brake disc thermomechanical fatigue is ignored at this moment; this simplification might lead to an underestimation in the thermomechanical fatigue damage and a blur in the thermomechanical fatigue mechanism reveal. Based on the analysis above, the study of multiaxial thermomechanical fatigue damage on the high-speed railway brake disc is insufficient. In this study, the thermomechan-According to previous research, the commonly used method for brake disc thermomechanical fatigue investigation is adopting the distribution and amplitude of temperature, stress and strain items. However, due to the rotating movement, the thermomechanical fatigue should belong to multiaxial fatigue damage [18], which is complex with numerous influencing factors [19]. Regrettably, the multiaxial feature of brake disc thermomechanical fatigue is ignored at this moment; this simplification might lead to an underestimation in the thermomechanical fatigue damage and a blur in the thermomechanical fatigue mechanism reveal.

ical fatigue damage of the brake disc will be studied with the help of a multiaxial fatigue model instead of using the distribution and amplitude of temperature, stress and strain items as fatigue indicators. A finite element simulation will be used to reproduce the multiaxial material response state under thermomechanical load. Afterward, a multiaxial fatigue model, recently proposed by the first author, will be adopted to analyze the multiaxial thermomechanical fatigue damage and out-of-phase failure status. Finally, the methodology for multiaxial thermomechanical fatigue evaluation will be proposed and the fatigue life prediction will be performed. **2. Finite Element Model** Based on the analysis above, the study of multiaxial thermomechanical fatigue damage on the high-speed railway brake disc is insufficient. In this study, the thermomechanical fatigue damage of the brake disc will be studied with the help of a multiaxial fatigue model instead of using the distribution and amplitude of temperature, stress and strain items as fatigue indicators. A finite element simulation will be used to reproduce the multiaxial material response state under thermomechanical load. Afterward, a multiaxial fatigue model, recently proposed by the first author, will be adopted to analyze the multiaxial thermomechanical fatigue damage and out-of-phase failure status. Finally, the methodology for multiaxial thermomechanical fatigue evaluation will be proposed and the fatigue life prediction will be performed.

#### 2.1.1. Thermal formulation **2. Finite Element Model**

#### Ignoring the influence of wheel–rail friction and aerodynamic resistance, most vehi-*2.1. Theoretical Background*

#### cle kinetic energy is transferred into heat by friction contact between the brake disc and 2.1.1. Thermal Formulation

pad when mechanical braking is implemented. Therefore, the total thermal energy flowed Ignoring the influence of wheel–rail friction and aerodynamic resistance, most vehicle kinetic energy is transferred into heat by friction contact between the brake disc and pad when mechanical braking is implemented. Therefore, the total thermal energy flowed into

*2.1. Theoretical Background* 

by Equation (3):

(see Equation (4)):

*2.2. FE Model Description* 

loading method.

(5):

() is described in Equation (1):

transfer () is demonstrated in Equation (2):

the braking system can be calculated by friction energy. The corresponding heat flux (*q <sup>f</sup>* ) is described in Equation (1): () is described in Equation (1): ∙ ∙ into the braking system can be calculated by friction energy. The corresponding heat flux

into the braking system can be calculated by friction energy. The corresponding heat flux

*Vehicles* **2021**, *3*, FOR PEER REVIEW 3

$$q\_f = \frac{\eta \cdot \mathbf{P} \cdot \mu}{A} \cdot (v\_0 + a \cdot t) \tag{1}$$
 
$$\dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots \qquad \dots$$

where *η* is the conversion rate between friction energy and heat, *P* denotes the contact force and *µ* stands for the friction coefficient, *A* is the braking contact area, *v*<sup>0</sup> is the initial relative frictional translation speed between the brake disc and pad, and *a* and *t* represent the acceleration and time, respectively. where is the conversion rate between friction energy and heat, denotes the contact force and stands for the friction coefficient, *A* is the braking contact area, <sup>0</sup> is the initial relative frictional translation speed between the brake disc and pad, and and represent the acceleration and time, respectively. = where is the conversion rate between friction energy and heat, denotes the contact force and stands for the friction coefficient, *A* is the braking contact area, <sup>0</sup> is the initial relative frictional translation speed between the brake disc and pad, and and rep-

Apart from the heat input, there are heat exchanges between the braking system and environment through convection heat transfer and heat radiation. The convection heat transfer (*q <sup>f</sup>* ) is demonstrated in Equation (2): Apart from the heat input, there are heat exchanges between the braking system and environment through convection heat transfer and heat radiation. The convection heat transfer () is demonstrated in Equation (2): resent the acceleration and time, respectively. Apart from the heat input, there are heat exchanges between the braking system and environment through convection heat transfer and heat radiation. The convection heat

$$q\_f = h \cdot (T - T\_0) \tag{2}$$

where *h* is the film coefficient, and *T* and *T*<sup>0</sup> represent the temperature of braking component and surrounding environment, respectively. The heat radiation (*q<sup>h</sup>* ) is calculated by Equation (3): where ℎ is the film coefficient, and and <sup>0</sup> represent the temperature of braking component and surrounding environment, respectively. The heat radiation (ℎ) is calculated by Equation (3): = ℎ ∙ ( − 0) (2) where ℎ is the film coefficient, and and <sup>0</sup> represent the temperature of braking component and surrounding environment, respectively. The heat radiation (ℎ) is calculated

$$q\_h = \mathbb{X} \cdot (T - T\_\*)^4 \tag{3}$$

where ж is the Stefan-Boltzmann constant, <sup>∗</sup> stands for the sink temperature. where <sup>ℎ</sup> = ж ∙ ( − <sup>∗</sup> ) 4 (3) is the Stefan-Boltzmann constant, *T*<sup>∗</sup> stands for the sink temperature.

#### 2.1.2. Thermal–mechanical solution where ж is the Stefan-Boltzmann constant, <sup>∗</sup> stands for the sink temperature. 2.1.2. Thermal–Mechanical Solution

∙ ∙

With regard to the solution technology for the coupled thermal–mechanical problem, the full Newton method is usually implemented with the enonsymmetric Jacobian matrix (see Equation (4)): 2.1.2. Thermal–mechanical solution With regard to the solution technology for the coupled thermal–mechanical problem, With regard to the solution technology for the coupled thermal–mechanical problem, the full Newton method is usually implemented with the enonsymmetric Jacobian matrix (see Equation (4)):

$$\begin{aligned} \left[ \begin{array}{c} \text{K}\_{\text{df}} \\ \text{K}\_{\text{fd}} \end{array} \right] \left\{ \begin{array}{c} \text{K}\_{\text{fd}} & \text{K}\_{\text{fd}} \\ \text{K}\_{\text{fd}} & \text{K}\_{\text{fd}} \end{array} \right] \left\{ \begin{array}{c} \Delta d \\ \Delta t \end{array} \right\} = \left\{ \begin{array}{c} \text{R}\_{\text{d}} \\ \text{R}\_{t} \end{array} \right\} \end{aligned} \tag{4}$$

where ∆ and ∆ are increments of displacement and temperature, respectively. is the submatrices of fully coupled Jacobian matrix, and stand for the mechanical and thermal residual vectors. Particularly, with respect to the issues of the coupled thermal–mechanical problem in the brake disc, the coupling effect between the mechanical and thermal solutions is weak. Therefore, the off-diagonal submatrices are assumed to be ] { ∆ ∆ } = { } (4) where ∆ and ∆ are increments of displacement and temperature, respectively. is the submatrices of fully coupled Jacobian matrix, and stand for the mechanical and thermal residual vectors. Particularly, with respect to the issues of the coupled thermal–mechanical problem in the brake disc, the coupling effect between the mechanical where ∆*d* and ∆*t* are increments of displacement and temperature, respectively. *Kij* is the submatrices of fully coupled Jacobian matrix, *R<sup>d</sup>* and *R<sup>t</sup>* stand for the mechanical and thermal residual vectors. Particularly, with respect to the issues of the coupled thermal– mechanical problem in the brake disc, the coupling effect between the mechanical and thermal solutions is weak. Therefore, the off-diagonal submatrices are assumed to be zero in order to have a less costly solution. The solution equation is described in Equation (5):

zero in order to have a less costly solution. The solution equation is described in Equation

] {

∆

} = {

[

0

$$
\begin{bmatrix}
\begin{array}{cc}
\mathbf{K}\_{dd} & \mathbf{0} \\
\mathbf{0} & \mathbf{K}\_{lt}
\end{array}
\end{bmatrix}
\begin{Bmatrix}
\Delta d \\
\Delta t
\end{Bmatrix} = \begin{Bmatrix}
\begin{array}{cc}
\mathbf{R}\_d \\
\mathbf{R}\_t
\end{array}
\end{Bmatrix}
\tag{5}
$$

} (5)

#### 0 ∆ *2.2. FE Model Description*

(5):

[

[

*2.2. FE Model Description*  There are two common methods for the loading implementation on the brake disc; one is the heat flux functions method and the other one is the mechanical friction method. Particularly, the heat flux functions method can be subdivided into the rotating loading method and the uniform loading method. The advantage of the uniform loading method is the calculation efficiency, but the multiaxial material response state and relevant out-0 ] { ∆ } = { } (5) There are two common methods for the loading implementation on the brake disc; one is the heat flux functions method and the other one is the mechanical friction method. Particularly, the heat flux functions method can be subdivided into the rotating loading There are two common methods for the loading implementation on the brake disc; one is the heat flux functions method and the other one is the mechanical friction method. Particularly, the heat flux functions method can be subdivided into the rotating loading method and the uniform loading method. The advantage of the uniform loading method is the calculation efficiency, but the multiaxial material response state and relevant outof-phase failure feature during the braking process have to be revealed by the rotating loading method.

of-phase failure feature during the braking process have to be revealed by the rotating loading method. The three-dimensional geometry model of the investigated disc brake system is shown in Figure 2a. It can be seen that the disc brake system is complex, comprising the brake disc, brake pad, backplate, pad bracket and brake clamp, etc. In order to improve the calculation efficiency, the brake disc is built with some simplifications and only half of it is established. The finite element model of the brake disc can be seen Figure 2b. The description of the brake disc dimension can be seen in Figure 2c. The red solid area and method and the uniform loading method. The advantage of the uniform loading method is the calculation efficiency, but the multiaxial material response state and relevant outof-phase failure feature during the braking process have to be revealed by the rotating The three-dimensional geometry model of the investigated disc brake system is shown in Figure 2a. It can be seen that the disc brake system is complex, comprising the brake disc, brake pad, backplate, pad bracket and brake clamp, etc. In order to improve the calculation efficiency, the brake disc is built with some simplifications and only half of it is established. The finite element model of the brake disc can be seen Figure 2b. The description of the brake disc dimension can be seen in Figure 2c. The red solid area and The three-dimensional geometry model of the investigated disc brake system is shown in Figure 2a. It can be seen that the disc brake system is complex, comprising the brake disc, brake pad, backplate, pad bracket and brake clamp, etc. In order to improve the calculation efficiency, the brake disc is built with some simplifications and only half of it is established. The finite element model of the brake disc can be seen Figure 2b. The description of the brake disc dimension can be seen in Figure 2c. The red solid area and red dotted area in Figure 2c are the regions for rotating loading and uniform loading, respectively. rL1, rL2 and r<sup>M</sup> indicate the inner, outer and average friction radii. The geometry dimension can be found in Table 1. A total of 29 elements (denoted as ER1-ER29) are evenly distributed along radial direction (between r<sup>2</sup> and r3) and 3 elements (represented by EA1-EA3) are uniformly

arranged from the contact surface into the brake disc with a distance of 2 mm between each other. by EA1-EA3) are uniformly arranged from the contact surface into the brake disc with a distance of 2 mm between each other.

red dotted area in Figure 2c are the regions for rotating loading and uniform loading, respectively. rL1, rL2 and r<sup>M</sup> indicate the inner, outer and average friction radii. The geometry dimension can be found in Table 1. A total of 29 elements (denoted as ER1-ER29) are evenly distributed along radial direction (between r<sup>2</sup> and r3) and 3 elements (represented

**Figure 2.** (**a**) Three-dimensional geometry model of the investigated disc brake system, (**b**) simplified finite element model of the brake disc, (**c**) description of brake disc dimension. **Figure 2.** (**a**) Three-dimensional geometry model of the investigated disc brake system, (**b**) simplified finite element model of the brake disc, (**c**) description of brake disc dimension.

**Table 1.** Geometry dimension of brake disc. **Table 1.** Geometry dimension of brake disc.


The finite element model is established in ABAQUS [20] in order to ensure the calculation accuracy and improve the calculation efficiency. An 8-node three-dimensional linear coupled temperature-displacement reduced-integration brick element is used to discretize the model. The material properties are listed in Table 2. The loading (heat flux, normal load and tangential load) is realized by subroutines in ABAQUS [20]. The heat flux can be calculated by Equation (1), which changes with time. The normal load is resulted from the brake force, and the tangential load is resulted from the friction force. The inner ring where the brake disc is in contact with the axle is fully constrained. A global Cartesian coordinate system is employed, and the original point is the center of the brake disc. The z direction is the rotation axis and the rotation speed can be calculated according to the adopted braking condition. The transient coupled temperature-displacement calculation is performed. The finite element model is established in ABAQUS [20] in order to ensure the calculation accuracy and improve the calculation efficiency. An 8-node three-dimensional linear coupled temperature-displacement reduced-integration brick element is used to discretize the model. The material properties are listed in Table 2. The loading (heat flux, normal load and tangential load) is realized by subroutines in ABAQUS [20]. The heat flux can be calculated by Equation (1), which changes with time. The normal load is resulted from the brake force, and the tangential load is resulted from the friction force. The inner ring where the brake disc is in contact with the axle is fully constrained. A global Cartesian coordinate system is employed, and the original point is the center of the brake disc. The z direction is the rotation axis and the rotation speed can be calculated according to the adopted braking condition. The transient coupled temperature-displacement calculation is performed.

**Table 2.** Material properties of brake disc.


The braking operation parameters are listed in Table 3.: the initial braking speed is 250 km/h, the wheel radius is 430 mm, the transfer efficiency is 0.9, the frictional coefficient between brake disc and pad is assumed as 0.4, a uniform acceleration −1.389 m/s<sup>2</sup> is employed under a brake force of 17,000 N, the corresponding brake time is 50 s and the total brake distance is 1736.1 m. The environment temperature is 20 °C. The braking operation parameters are listed in Table 3: the initial braking speed is 250 km/h, the wheel radius is 430 mm, the transfer efficiency *η* is 0.9, the frictional coefficient between brake disc and pad is assumed as 0.4, a uniform acceleration <sup>−</sup>1.389 m/s<sup>2</sup> isemployed under a brake force of 17,000 N, the corresponding brake time is 50 s and the total brake distance is 1736.1 m. The environment temperature is 20 ◦C.

**Table 3.** The braking operation parameters.


### **3. Multiaxial Fatigue Model**

A fatigue model called the multiaxial fatigue space (referred as MFS criterion for convenience), recently proposed by the first author, is adopted to assess the multiaxial thermomechanical fatigue damage of the brake disc. The MFS criterion is validated by the fatigue test data of different materials under various loading conditions and displays a pleasant prediction capability, compared with other commonly used multiaxial fatigue models [21]. A brief introduction is presented, and detailed information can be found in the corresponding reference [21].

The fundamental assumption in the MFS criterion is that the changing of the material response state is the source of the fatigue damage. As shown in Figure 3, subscripts T1 and T2 indicate different time points. σ and τ are normal and shear stresses, subscripts x and y represent the components following the global coordinates, and subscripts 1 and 2 mean the components following the principal directions. σ<sup>1</sup> T10 and σ<sup>2</sup> T10 indicate the response status at moment T2 following the principal directions of T1. C and R are the center and radius of Mohr circle, and A is twice the angle between the loading plane and the principal plane. It can be seen that the changing of the material state can be described by the changing of C, R and A. Therefore, the material response status is described by three parameters (*P1*, *P2* and *P3*), which are the equivalent shear parameter, equivalent normal parameter and out-of-phase parameter, respectively [21]. The changing of them (three fatigue basic units, i.e., *f* 1(*P*1), *f* 2(*P*2) and *f* 3(*P*3)) can be used for fatigue damage evaluation, and the multiaxial fatigue space can be established by fatigue basic units (see Figure 4) [21]. The multiaxial fatigue parameter is described in Equation (6):

$$FP = f0(f1(P1), f2(P2), f3(P3))\tag{6}$$

with

$$P1 = \frac{\sqrt{\left(N\_{\rm C1} - N\_{\rm C2}\right)^2 + \left(S\_{\rm C1} - S\_{\rm C2}\right)^2}}{2} \tag{7}$$

$$P2 = \frac{N\_{\rm C1} + N\_{\rm C2}}{2} \tag{8}$$

$$P\mathfrak{3} = P1 \frac{\theta^2}{\theta\_{\text{max}}} \tag{9}$$

where *fn* indicates a certain function, *P1* and *P2* are in-phase parameters and *P3* is out-ofphase parameter [21], *N* and *S* indicate the normal and shear components on the critical plane (a couple of two orthogonal planes, noted by subscripts *C*1 and *C*2), stress, strain and energy components can be used as the adopted items, *θ*/2 is the out-of-phase angle, which represents the relative distortion deformation, and *θmax*/2 is the limit value of the out-of-phase angle [21]. It can be noticed that the MFS criterion belongs to the type of critical plane criterion. *Vehicles* **2021**, *3*, FOR PEER REVIEW 6

According to the MFS criterion, the amplitude of parameter variation (changing of *P1*, *P2* and *P3*) is the foundational component in multiaxial fatigue space, which is consistent with the definition of fatigue damage. The influence of the mean scales the vector *OP* in multiaxial fatigue space [21]. Therefore, the multiaxial fatigue parameter in Equation (6) can be expressed by amplitude parameter () and mean parameter (). See

= (1 + 11) ∙ (1 + 22) ∙ (1 + 33) (11)

)

) <sup>ℎ</sup>]

<sup>ℎ</sup> <sup>+</sup> ( <sup>∙</sup>

<sup>ℎ</sup> + ( ∙ 2

= 3

+ ()

where subscripts *a* and *m* indicate the amplitude and mean of corresponding parameter, subscripts *in* and *out* denote the in-phase and out-of-phase failure conditions, *k* and *h* are material-dependent coefficients, and is the mean influence coefficient. Afterward, the predicted fatigue life () can be obtained by Equations (15) or (16), depending on the

= ∙ (10)

<sup>ℎ</sup> ]

1⁄ℎ (12)

1⁄ℎ (13)

(140

(15)

**Figure 3.** Material response states at different time instants. **Figure 3.** Material response states at different time instants.

**Figure 4.** Multiaxial fatigue space.

= [(

)

)

= A()

= [(1

Equation (10):

fatigue curve pattern:

with

**Figure 3.** Material response states at different time instants.

*Vehicles* **2021**, *3*, FOR PEER REVIEW 6

**Figure 4.** Multiaxial fatigue space. **Figure 4.** Multiaxial fatigue space.

According to the MFS criterion, the amplitude of parameter variation (changing of *P1*, *P2* and *P3*) is the foundational component in multiaxial fatigue space, which is consistent with the definition of fatigue damage. The influence of the mean scales the vector *OP* in multiaxial fatigue space [21]. Therefore, the multiaxial fatigue parameter in Equation (6) can be expressed by amplitude parameter () and mean parameter (). See According to the MFS criterion, the amplitude of parameter variation (changing of *P1*, *P2* and *P3*) is the foundational component in multiaxial fatigue space, which is consistent with the definition of fatigue damage. The influence of the mean scales the vector *OP* in multiaxial fatigue space [21]. Therefore, the multiaxial fatigue parameter in Equation (6) can be expressed by amplitude parameter (*FPa*) and mean parameter (*FPm*). See Equation (10):

$$FP = FP\_a \cdot FP\_m \tag{10}$$

with with

Equation (10):

$$FP\_m = (1 + \chi\_{P1} P1\_m) \cdot (1 + \chi\_{P2} P2\_m) \cdot (1 + \chi\_{P3} P3\_m) \tag{11}$$

$$FP\_a = \left[ \left( FP\_a^{in} \right)^{h\_{out}} + \left( k\_{out} \cdot FP\_a^{out} \right)^{h\_{out}} \right]^{1/h\_{out}} \tag{12}$$

$$\mathbf{F}P\_a^{in} = \left[ (P\mathbf{1}\_a)^{h\_{in}} + (k\_{in} \cdot \mathbf{P2}\_a)^{h\_{in}} \right]^{1/h\_{in}} \tag{13}$$

$$FP\_a^{out} = P\mathfrak{B}\_a \tag{14}$$

 = [(1 ) <sup>ℎ</sup> + ( ∙ 2 ) <sup>ℎ</sup>] 1⁄ℎ (13) = 3 (140 where subscripts *a* and *m* indicate the amplitude and mean of corresponding parameter, where subscripts *a* and *m* indicate the amplitude and mean of corresponding parameter, subscripts *in* and *out* denote the in-phase and out-of-phase failure conditions, *k* and *h* are material-dependent coefficients, and *χPx* is the mean influence coefficient. Afterward, the predicted fatigue life (*N<sup>f</sup>* ) can be obtained by Equations (15) or (16), depending on the fatigue curve pattern:

$$FP = \mathbf{A} \left(\mathbf{N}\_f\right)^B + \mathbf{C} \left(\mathbf{N}\_f\right)^D \tag{15}$$

$$FP = \mathbf{A} \left(\mathbf{N}\_f\right)^B \tag{16}$$

fatigue curve pattern: = A() + () (15) where *A*, *B*, *C* and *D* are material-dependent coefficients with solid physical meaning in multiaxial fatigue space [21].

It should be noticed that the original MFS criterion is established and validated in a 2-dimensional condition. With regard to 3-dimensional cases, the general formula is the same; the difference is in the determination of the critical plane. In this work, the directions of the maximum and the second maximum principal strain range are employed for critical plane determination. That is to say, the direction of the maximum principal strain range will be determined firstly. Afterward, the directions perpendicular to the maximum principal strain range will be checked to obtain the direction of the second maximum principal strain range. In this way, the issues of the 3-dimensional condition can be transformed into 2-dimensional cases and the MFS criterion can be used smoothly.

About the material dependent coefficients, normally the fatigue test is performed under mechanical fatigue loading instead of thermomechanical fatigue loading. The thermomechanical fatigue-dependent coefficients of the brake disc are not available at

this moment. Thus, it is assumed that they are identical with the fatigue dependent coefficients of hot-rolled 45 steel [21]. In order to present the methodology, this is a common strategy when the fatigue-related information is unavailable [22,23]. The adopted fatiguedependent coefficients are listed in Table 4. An accurate multiaxial thermomechanical fatigue evaluation of the brake disc can be performed following the proposed approach when the thermomechanical fatigue-dependent coefficients are collected in the future.

**Table 4.** The adopted material-dependent coefficients.


### **4. Results and Discussion**

Considering the calculation efficiency, uniform loading is performed first, the calculation time is set as 500 s to have an overall response status, the highest temperature 454 ◦C occurs at 35 s, and the temperature distribution on the brake disc at four moments before the highest temperature time is demonstrated in Figure 5. As discussed above, the multiaxial characters cannot be revealed with uniform loading. Rotating loading is needed for the multiaxial thermomechanical fatigue evaluation. According to the results with uniform loading, it can be seen that the temperature increases significantly and reaches above 400 ◦C in the first 20 s. Thus, the calculation time is set as 20 s with rotating loading so as to reduce the simulation time and result file size. The comparison between the results with rotating loading and uniform loading can be seen in Figure 6. *Vehicles* **2021**, *3*, FOR PEER REVIEW 8

**Figure 5.** Temperature distribution with uniform loading at different moments (highest temperature 454 °C occurs at t4 = 35 s). **Figure 5.** Temperature distribution with uniform loading at different moments (highest temperature 454 ◦C occurs at t4 = 35 s).

(**a**) Temperature in radial direction (**b**) Huber–Mises in radial direction It can be seen in Figure 6 that the results manifestation is reasonable [7,17], which verifies the correctness of the calculation to a certain extent. The difference between the results with uniform and rotating loading is negligible. It can also be noticed that the heating time (50 s) is much shorter than the cooling time (much greater than 450 s) and the results variation during the heating process (see the changing of Huber–Mises results in the first 50 s in Figure 6) is much more violent than that during the cooling process (see the Huber–Mises results from 50 s to 500 s in Figure 6). Thus, as the changing of the material state during the heating process is used for fatigue evaluation, the calculation time is greatly reduced without losing much accuracy. To have a clear view about the multiaxial character of thermomechanical fatigue damage of the brake disc, the evolutions of the strain items with rotating loading can be seen in Figure 7. It can be noticed clearly that the changing of shear strain is non-proportional, which indicates that the principal strain direction is changed during the braking process. It means that the thermomechanical fatigue damage of the brake disc not only belongs to the type of multiaxial fatigue damage, but also the out-of-phase failure status. The out-of-phase failure of the brake disc could be resulted from the mechanical interaction between the brake disc and the brake pad,

or could be resulted from the influence of the brake disc structure. The reasons will be discussed in the future.

**Figure 6.** Comparison between the results with uniform loading and rotating loading. **Figure 6.** Comparison between the results with uniform loading and rotating loading.

thermomechanical fatigue damage of the brake disc not only belongs to the type of multiaxial fatigue damage, but also the out-of-phase failure status. The out-of-phase failure of the brake disc could be resulted from the mechanical interaction between the brake disc and the brake pad, or could be resulted from the influence of the brake disc **Figure 7.** (**a**) The changing of strain items following xx, yy and zz directions with rotating loading, (**b**) the changing of strain items following xy, xz and yz directions with rotating loading. **Figure 7.** (**a**) The changing of strain items following xx, yy and zz directions with rotating loading, (**b**) the changing of strain items following xy, xz and yz directions with rotating loading.

structure. The reasons will be discussed in the future. It is known that the changing of the material response state is critical for multiaxial fatigue analysis; the time increment needs to be small enough to capture the details of this changing. This makes the realization with rotating loading unfavorable because of the lengthy calculation time and huge resulting file size. Thus, an equivalent accelerated method (referred to as the equivalent method) is proposed. An amplification factor is used to scale the heat flux in order to reach the target temperature with only one rotation. In It is known that the changing of the material response state is critical for multiaxial fatigue analysis; the time increment needs to be small enough to capture the details of this changing. This makes the realization with rotating loading unfavorable because of the lengthy calculation time and huge resulting file size. Thus, an equivalent accelerated method (referred to as the equivalent method) is proposed. An amplification factor is used to scale the heat flux in order to reach the target temperature with only one rotation. In

**Figure 8.** Temperature distribution with equivalent method at different times in one rotation.

are used for the fatigue evaluation with the MFS criterion.

order to validate equivalent method, the temperature with rotating loading at 20 s is used as the target temperature (406.73 °C). The accepted tolerance for the obtained temperature

It can be seen that the highest temperature by the equivalent method (406.38 °C) is close to that with rotating loading (406.73 °C). The position of the maximum temperature by the equivalent method is located at the average friction radius position (ER15), which is the same as the rotating loading calculation. The comparison of the results among uniform loading, rotating loading and equivalent methods at the maximum temperature position is shown in Table 5. It can be seen that the temperature, thermal strain and Huber–Mises stress are approximately equal with each other. The results with rotating loading is used as reference. The temperature and strain with the equivalent method are quite satisfactory and the relative error is less than 0.1%. The difference of Huber–Mises stress is around 5% because of the history effect and constraint from the extrusion by the surrounding elements. In general, the results with the equivalent method are acceptable and strain items order to validate equivalent method, the temperature with rotating loading at 20 s is used as the target temperature (406.73 ◦C). The accepted tolerance for the obtained temperature by equivalent method is 1%. The temperature distribution with equivalent method at different times in one rotation is shown in Figure 8; the highest temperature is 406.38 ◦C. to scale the heat flux in order to reach the target temperature with only one rotation. In order to validate equivalent method, the temperature with rotating loading at 20 s is used as the target temperature (406.73 °C). The accepted tolerance for the obtained temperature by equivalent method is 1%. The temperature distribution with equivalent method at different times in one rotation is shown in Figure 8; the highest temperature is 406.38 °C.

It is known that the changing of the material response state is critical for multiaxial fatigue analysis; the time increment needs to be small enough to capture the details of this changing. This makes the realization with rotating loading unfavorable because of the lengthy calculation time and huge resulting file size. Thus, an equivalent accelerated method (referred to as the equivalent method) is proposed. An amplification factor is used

(**a**) (**b**) **Figure 7.** (**a**) The changing of strain items following xx, yy and zz directions with rotating loading, (**b**) the changing of

strain items following xy, xz and yz directions with rotating loading.

*Vehicles* **2021**, *3*, FOR PEER REVIEW 9

**Figure 8.** Temperature distribution with equivalent method at different times in one rotation. **Figure 8.** Temperature distribution with equivalent method at different times in one rotation.

It can be seen that the highest temperature by the equivalent method (406.38 °C) is close to that with rotating loading (406.73 °C). The position of the maximum temperature by the equivalent method is located at the average friction radius position (ER15), which is the same as the rotating loading calculation. The comparison of the results among uniform loading, rotating loading and equivalent methods at the maximum temperature position is shown in Table 5. It can be seen that the temperature, thermal strain and Huber–Mises stress are approximately equal with each other. The results with rotating loading is used as reference. The temperature and strain with the equivalent method are quite satisfactory and the relative error is less than 0.1%. The difference of Huber–Mises stress is around 5% because of the history effect and constraint from the extrusion by the surrounding elements. In general, the results with the equivalent method are acceptable and strain items are used for the fatigue evaluation with the MFS criterion. It can be seen that the highest temperature by the equivalent method (406.38 ◦C) is close to that with rotating loading (406.73 ◦C). The position of the maximum temperature by the equivalent method is located at the average friction radius position (ER15), which is the same as the rotating loading calculation. The comparison of the results among uniform loading, rotating loading and equivalent methods at the maximum temperature position is shown in Table 5. It can be seen that the temperature, thermal strain and Huber–Mises stress are approximately equal with each other. The results with rotating loading is used as reference. The temperature and strain with the equivalent method are quite satisfactory and the relative error is less than 0.1%. The difference of Huber–Mises stress is around 5% because of the history effect and constraint from the extrusion by the surrounding elements. In general, the results with the equivalent method are acceptable and strain items are used for the fatigue evaluation with the MFS criterion.


**Table 5.** Comparison of the results among uniform, rotating and equivalent methods.

Regarding the multiaxial fatigue prediction capability with the equivalent method, the comparison with rotating loading is performed at different temperatures. Because the results of the finite element calculation are described in the global coordinate (coordinate xyz), coordinate transformation is executed with the help of Euler angles to get the material response status in coordinate XYZ (see Figure 9). Here, zXZ rotation order is adopted for the coordinate transformation. Afterward, the MFS criterion can be adopted. The calculated fatigue parameter and relative error (defined in Equation (17)) are demonstrated in Figure 10.

$$E\_r = \left(FP\_{rotation} - FP\_{equivalent}\right) / FP\_{rotation} \tag{17}$$

It can be seen that with the increasement of temperature, the relative error between the equivalent method and the rotating loading method is increased slowly, following a linear pattern, approximately. Based on this, the relationship function between the relative error (*Er*) and temperature (*T*) can be written as Equation (18):

$$E\_r = 0.00032637 \times T - 0.017 \tag{18}$$

with error correction (see Figure 10, red point).

in Figure 10.

error (

*Vehicles* **2021**, *3*, FOR PEER REVIEW 10

**Table 5.** Comparison of the results among uniform, rotating and equivalent methods.

**Item Rotating Uniform Equivalent** 

Temperature (°C) 406.73 417.82 −2.73 406.38 0.086 Thermal strain (%) 0.40220 0.41373 −2.87 0.40184 0.090 Huber–Mises (MPa) 1063.5 1082.8 −1.81 1117.8 −5.1

Regarding the multiaxial fatigue prediction capability with the equivalent method,

It can be seen that with the increasement of temperature, the relative error between

Now the multiaxial thermomechanical fatigue of the brake disc under the investi-

gated braking operation condition can be assessed. Due to the fact that the calculation time with rotating loading is quite long, the maximum temperature with uniform loading is used as the target temperature (454 °C, see Figure 6); the value of the fatigue parameter with the equivalent method can be obtained with this target temperature (see Figure 10, blue point). Hereafter, the fatigue parameter value with rotating loading can be obtained

the equivalent method and the rotating loading method is increased slowly, following a linear pattern, approximately. Based on this, the relationship function between the relative

) and temperature () can be written as Equation (18):

= ( − )⁄ (17)

= 0.00032637 × − 0.017 (18)

the comparison with rotating loading is performed at different temperatures. Because the results of the finite element calculation are described in the global coordinate (coordinate xyz), coordinate transformation is executed with the help of Euler angles to get the material response status in coordinate XYZ (see Figure 9). Here, zXZ rotation order is adopted for the coordinate transformation. Afterward, the MFS criterion can be adopted. The calculated fatigue parameter and relative error (defined in Equation (17)) are demonstrated

**Value Error % Value Error %**

**Figure 9. Figure 9.**  Coordinate transformation with Euler angles. Coordinate transformation with Euler angles.

**Figure 10.** Fatigue parameter comparison between rotating loading and equivalent methods. **Figure 10.** Fatigue parameter comparison between rotating loading and equivalent methods.

It is discussed above that the amplitudes and distributions of temperature, stress and strain items are not enough for accurate thermomechanical fatigue evaluation of the brake disc. Ignoring the influence of the multiaxial character and out-of-phase feature leads to an underestimation of the thermomechanical fatigue assessment for the brake disc. Based Now the multiaxial thermomechanical fatigue of the brake disc under the investigated braking operation condition can be assessed. Due to the fact that the calculation time with rotating loading is quite long, the maximum temperature with uniform loading is used as the target temperature (454 ◦C, see Figure 6); the value of the fatigue parameter with the equivalent method can be obtained with this target temperature (see Figure 10, blue point). Hereafter, the fatigue parameter value with rotating loading can be obtained with error correction (see Figure 10, red point).

on the result with the rotating loading calculation, the contribution of the equivalent shear parameter, equivalent normal parameter and out-of-phase parameter (represented by *P1*, *P2* and *P3* in the MFS criterion) to the total multiaxial thermomechanical fatigue parameter is demonstrated in Figure 11. It can be seen that just using a uniaxial fatigue parameter causes at least 14% underestimation of the fatigue damage, while employing a multiaxial fatigue parameter without the consideration of the out-of-phase failure leads to an underestimation of 5%. It should be noticed that the relationship between the fatigue parameter and fatigue life is exponential, which means that an insignificant deviation in the fatigue It is discussed above that the amplitudes and distributions of temperature, stress and strain items are not enough for accurate thermomechanical fatigue evaluation of the brake disc. Ignoring the influence of the multiaxial character and out-of-phase feature leads to an underestimation of the thermomechanical fatigue assessment for the brake disc. Based on the result with the rotating loading calculation, the contribution of the equivalent shear parameter, equivalent normal parameter and out-of-phase parameter (represented by *P1*, *P2* and *P3* in the MFS criterion) to the total multiaxial thermomechanical fatigue parameter is demonstrated in Figure 11. It can be seen that just using a uniaxial fatigue parameter causes at least 14% underestimation of the fatigue damage, while employing a multiaxial fatigue parameter without the consideration of the out-of-phase failure leads to an underestimation of 5%. It should be noticed that the relationship between the fatigue

parameter assessment leads to a significant predicted error of fatigue life.

**Figure 11.** The contribution of different parameters on multiaxial thermomechanical fatigue of the

with temperature [24]–[26]. The material character under high temperature is not available at this moment. As there is an approximate linear relationship between the yield stress and fatigue limit [27], the yield stress decreases exponentially with the increasing of the temperature, approximately (see Figure 12a) [15]. Thus, an exponential relationship between the fatigue limit and temperature is assumed in order to present the methodology in this work (see Figure 12b). The corresponding fatigue curves at different temperatures can be obtained using extrapolation method with the help of the yield stress at different temperatures. The multiaxial thermomechanical fatigue life of the investigated brake disc under the employed braking operation condition can be obtained, which is about 2500 braking cycles. The accurate multiaxial thermomechanical fatigue life prediction can be obtained with the proposed methodology when the temperature related material fatigue

It is known that the material mechanical and fatigue characters change significantly

characters become available in the future.

brake disc.

parameter and fatigue life is exponential, which means that an insignificant deviation in the fatigue parameter assessment leads to a significant predicted error of fatigue life. and fatigue life is exponential, which means that an insignificant deviation in the fatigue parameter assessment leads to a significant predicted error of fatigue life.

**Figure 10.** Fatigue parameter comparison between rotating loading and equivalent methods.

It is discussed above that the amplitudes and distributions of temperature, stress and

strain items are not enough for accurate thermomechanical fatigue evaluation of the brake disc. Ignoring the influence of the multiaxial character and out-of-phase feature leads to an underestimation of the thermomechanical fatigue assessment for the brake disc. Based on the result with the rotating loading calculation, the contribution of the equivalent shear parameter, equivalent normal parameter and out-of-phase parameter (represented by *P1*, *P2* and *P3* in the MFS criterion) to the total multiaxial thermomechanical fatigue parameter is demonstrated in Figure 11. It can be seen that just using a uniaxial fatigue parameter causes at least 14% underestimation of the fatigue damage, while employing a multiaxial fatigue parameter without the consideration of the out-of-phase failure leads to an underestimation of 5%. It should be noticed that the relationship between the fatigue parameter

*Vehicles* **2021**, *3*, FOR PEER REVIEW 11

**Figure 11.** The contribution of different parameters on multiaxial thermomechanical fatigue of the brake disc. **Figure 11.** The contribution of different parameters on multiaxial thermomechanical fatigue of the brake disc.

It is known that the material mechanical and fatigue characters change significantly with temperature [24]–[26]. The material character under high temperature is not available at this moment. As there is an approximate linear relationship between the yield stress and fatigue limit [27], the yield stress decreases exponentially with the increasing of the temperature, approximately (see Figure 12a) [15]. Thus, an exponential relationship between the fatigue limit and temperature is assumed in order to present the methodology in this work (see Figure 12b). The corresponding fatigue curves at different temperatures can be obtained using extrapolation method with the help of the yield stress at different temperatures. The multiaxial thermomechanical fatigue life of the investigated brake disc It is known that the material mechanical and fatigue characters change significantly with temperature [24–26]. The material character under high temperature is not available at this moment. As there is an approximate linear relationship between the yield stress and fatigue limit [27], the yield stress decreases exponentially with the increasing of the temperature, approximately (see Figure 12a) [15]. Thus, an exponential relationship between the fatigue limit and temperature is assumed in order to present the methodology in this work (see Figure 12b). The corresponding fatigue curves at different temperatures can be obtained using extrapolation method with the help of the yield stress at different temperatures. The multiaxial thermomechanical fatigue life of the investigated brake disc under the employed braking operation condition can be obtained, which is about 2500 braking cycles. The accurate multiaxial thermomechanical fatigue life prediction can be obtained with the proposed methodology when the temperature related material fatigue characters become available in the future. *Vehicles* **2021**, *3*, FOR PEER REVIEW 12

> In this study, the multiaxial characters of thermomechanical fatigue damage on the high-speed railway brake disc is investigated. Through the revealing of the multiaxial re-

> type of multiaxial thermomechanical fatigue. With the help of a multiaxial fatigue model, the methodology for multiaxial thermomechanical fatigue damage evaluation is proposed and multiaxial thermomechanical fatigue life prediction is implemented. In addition, the contribution of the uniaxial fatigue parameter, multiaxial fatigue parameter and out-ofphase failure parameter to the total multiaxial thermomechanical fatigue damage of the brake disc is discussed. The results indicate that, using the amplitude and distribution of temperature, stress and strain items for fatigue evaluation lead to an underestimation of the brake disc thermomechanical fatigue damage. This work explains the importance of studying the thermomechanical fatigue damage of brake disc from the perspective of mul-

> **Author Contributions:** Conceptualization, C.L. and J.M.; methodology, C.L. and Z.F.; software, C.L. and Y.W.; validation, C.L. and J.M.; formal analysis, C.L.; investigation, C.L. and R.S.; resources, C.L. and J.M.; data curation, C.L. and Y.W.; writing—original draft preparation, C.L., R.S., Y.W. and Z.F.; writing—review and editing, J.M.; visualization, C.L. and R.S.; supervision, J.M.; project administration, C.L. and J.M.; funding acquisition, C.L. and J.M. All authors have read and agreed to

> **Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51822508; Sichuan Science and Technology Program, grant number 2020JDTD0012; the Fun-

damental Research Funds for the Central Universities, grant number 2682021CX028.

**Figure 12.** (**a**) The relationship between the yield stress and temperature [15], (**b**) extrapolation of fatigue curve. **Figure 12.** (**a**) The relationship between the yield stress and temperature [15], (**b**) extrapolation of fatigue curve.

**5. Conclusions**

tiaxial fatigue.

*Int. J. Fatigue* **2017**, *104*, 99–111, doi:10.1016/j.ijfatigue.2017.07.018.

*Appl. Mech. Lett.* **2019**, *9*, 308–311, doi:10.1016/j.taml.2019.06.008.

by simulatin. *J. Mech. Eng.* **2011**, *47*, 126–131.

**References**

the published version of the manuscript.

**Institutional Review Board Statement:** Not applicable.

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

1. Lu, C.; Melendez, J.; Martínez-Esnaola, J. Fatigue damage prediction in multiaxial loading using a new energy-based parameter.

2. Yevtushenko, A.; Kuciej, M.; Wasilewski, P. Experimental study on the temperature evolution in the railway brake disc. *Theor.* 

3. Zhou, S.; Yang, Y.; Xie, J. Analysis of transient temperature and thermal stress distribution on the high-speed strain brake disk

**Informed Consent Statement:** Not applicable.

### **5. Conclusions**

In this study, the multiaxial characters of thermomechanical fatigue damage on the high-speed railway brake disc is investigated. Through the revealing of the multiaxial response state and the out-of-phase failure feature, it is proved that under mechanical friction braking conditions, the damage of the high-speed railway brake disc belongs to the type of multiaxial thermomechanical fatigue. With the help of a multiaxial fatigue model, the methodology for multiaxial thermomechanical fatigue damage evaluation is proposed and multiaxial thermomechanical fatigue life prediction is implemented. In addition, the contribution of the uniaxial fatigue parameter, multiaxial fatigue parameter and out-ofphase failure parameter to the total multiaxial thermomechanical fatigue damage of the brake disc is discussed. The results indicate that, using the amplitude and distribution of temperature, stress and strain items for fatigue evaluation lead to an underestimation of the brake disc thermomechanical fatigue damage. This work explains the importance of studying the thermomechanical fatigue damage of brake disc from the perspective of multiaxial fatigue.

**Author Contributions:** Conceptualization, C.L. and J.M.; methodology, C.L. and Z.F.; software, C.L. and Y.W.; validation, C.L. and J.M.; formal analysis, C.L.; investigation, C.L. and R.S.; resources, C.L. and J.M.; data curation, C.L. and Y.W.; writing—original draft preparation, C.L., R.S., Y.W. and Z.F.; writing—review and editing, J.M.; visualization, C.L. and R.S.; supervision, J.M.; project administration, C.L. and J.M.; funding acquisition, C.L. and J.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51822508; Sichuan Science and Technology Program, grant number 2020JDTD0012; the Fundamental Research Funds for the Central Universities, grant number 2682021CX028.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

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

### **References**

