Dynamic response and stability analysis of brake chatter with stochastic parameter perturbation

Abstract

The parameter perturbation makes it difficult to analyze the dynamic mechanism of the brake chatter with the help of an ideal deterministic model. Hence, a stochastic dynamic model of brake chatter with the uncertainty of the torsional stiffness of the brake disc is established, the Stribeck model is used to describe the friction characteristic between the brake pad and brake disc. Then, the Itô stochastic differential equation is solved by means of the stochastic averaging method, the boundary type of the one-dimensional energy process is identified to discuss the stochastic stability of the brake system. Furthermore, the Fokker Planck Kolmogorov equation is derived, the dynamic bifurcation and phenomenological bifurcation are analytically proven, the influences of parameters, such as noise intensity, brake pressure, and friction coefficient, on bifurcation characteristics are revealed. The results show that the noise intensity, brake pressure, and friction coefficient difference need to be reduced to improve the system stability. Finally, the probability density function of the brake system is used to validate the analytical analysis. The relevant results provide a theoretical basis for better suppressing the brake chatter.

Keywords

disc brake brake chatter stability analysis stochastic bifurcation probability density

Introduction

Disc brake, with excellent heat dissipation, stability, and response speed, has become the mainstream choice for the brake system, its dynamic performance directly affects the driving safety of vehicles (Chen et al., 2025; Li et al., 2025; İbrahim Cana et al., 2025). However, the brake chatter of the disc brake usually occurs during vehicle braking, typically manifested as the shaking of the brake pedal and vehicle body. This phenomenon not only affects the driving comfort, braking performance, and components lifespan of the vehicle, but may also pose serious risks to the driving safety, many scholars have focused on the nonlinear dynamics of the brake chatter (Cantoni et al., 2009; Chatelet et al., 2008).

Since the disc brake depends on the friction to achieve vehicle braking, the brake chatter is manifested as the friction-induced vibration (Gu et al., 2023, 2024; Lu and Zhao, 2022; Meehan and Leslie, 2021), the influence of the friction characteristics between the brake disc and brake pad on the brake chatter has been studied (Marques et al., 2021; Zhang et al., 2024). Rhee et al. (1991) studied the friction film formation of the brake rotor and pad, the results indicated that the brake vibration and noise are caused by the transfer film destruction. Thomsen and Fidlin (2003) derived an approximate analytical expression of the brake system, the influence of the difference between static friction and dynamic friction on the system response was analyzed by means of time series diagrams, phase diagrams, and amplitude response diagrams. With the help of the analytical analysis and experimental data, Crowther and Singh (2008) analyzed the influence of the friction nonlinear coupling brake and rotor on the brake groan phenomenon. Manish et al. (2005) pointed out that the growth of the friction layer has a significant impact on the contact stiffness, thereby affecting the system stability. For the dynamic coupling between the coefficient of friction, contact temperature and sliding velocity of the disc brake, Grzes and Kuciej (2025) proposed a conversion method of the measured time profiles of them into analytical formulation of two variables.

In order to reveal the dynamic mechanism of the brake chatter of the disc brake, some scholars simplify the disc brake as a single degree-of-freedom (DOF) mass-on-moving-belt system, which is used to investigate the dynamic characteristics of the braking system. For example, Andreaus and Casini (2001) focused on a single DOF model with Coulomb friction, solved the closed-form solution of the system. Then, the influence of the initial speed on the system response was discussed. Hetzler et al. (2007) analytically proved that the subcritical Hopf bifurcation leads to an interesting phenomenon, that is, the system response changes from an unstable fixed-point to an unstable limit cycle. Chen and Xi (2014) proposed a single DOF torsional model of driveline with wedge brake, the influences of initial actuation force, excitation frequency, and wedge angle on the system stability were investigated.

Due to the dynamic coupling between the brake pad and brake disc through frictional characteristics, dynamic models with different degrees-of-freedom (DOFs) of the brake system have also established (Akif and Osman, 2023; Shin et al., 2002; Zeng and Luo, 2011). Considering both the vertical motion of the brake unit and the pitch motion of the bogie frame of the disc brake, Zeng and Luo (2011) established a two DOFs model of disc brake, the results indicated that the critical speed depends on the mass of the brake unit, the vertical damping of the suspension, the friction coefficient, and the brake normal force between the lining and disc. Wei et al. (2019) developed the three DOFs of a new brake system with double-layer pad, the impact of the pad mass discussion and the connection stiffness on the dynamic responses of the system was investigated by the numerical analysis. Du et al. (2024) derived the three DOFs and four DOFs models of the three-body brake system for different brake conditions, the effects of the road conditions and particles on the dynamic characteristics of the brake system pairs were examined. Wang et al. (2024) considered the nonlinear characteristics of the brake system, established the vehicle-track coupled model with the disc brake, which can be used to analyze the brake squeal and friction-induced vibration.

In addition, the brake system is a complex dynamic system, and its dynamic performances can be influenced by multiple factors. Hence, the influence of various factors on the brake chatter has also been discussed. Shin et al. (2002) found that the damping of the brake disc is as important as that of the brake pad, the nonlinear analysis was performed to obtain the unstable boundary and demonstrated limit cycles in the phase space. Kang et al. (2009) pointed out that the frictional mode-coupling is a major mechanism of the system instability. However, the rotational motion of the disc has little influence on flutter mode. With the help of the finite element method and eigenvalue analysis, Fischer et al. (2021) found that the flexibilities in the suspension bushings and the elastic deformations of the suspension part and tire exacerbated the self-excited vibration of the brake system. Song et al. (2025) explained the effects of the lateral external forces, the primary longitudinal stiffness, and the conicity on the optimization effect. Wang and Lei (2024) pointed out that the presence of short wavelength irregularity aggravates the stick-slip vibration.

The achievements mentioned above are based on the deterministic assumption of system parameters to study the dynamic mechanism of the brake chatter, which provides a theoretical basis for further researches on the nonlinear dynamics of the brake chatter. In particular, the vehicle is a typical nonlinear system, the model simplification, measurement error, material aging, and other factors can lead to the stochastic parameter perturbation, this can significantly affect the dynamic responses of the vehicle (Hu et al., 2025; Lü et al., 2017; Sarrouy et al., 2013; Zhao et al., 2024). Therefore, this paper deals with the dynamic response and stability analysis of the brake chatter with the stochastic parameter perturbation.

The rest of this paper is arranged as follows. Firstly, a two DOFs mathematical model considering stochastic parameter perturbation of torsional stiffness is derived. Then, the Itô stochastic differential equation is solved, the boundary type of the one-dimensional energy process is identified to discuss the stochastic stability of the brake system. Furthermore, the bifurcation characteristics of the brake system with the stochastic parameter perturbation are explored, numerical verifications are carried out. Finally, conclusions are summarized.

Mathematical equations

Mechanical model

Disc brake is mainly composed of the brake caliper and brake disc. To derive mathematical equations of the disc brake, as shown in Figure 1, the coordinate system XYZ is established with the center of mass (COM) of the brake disc as the origin and $ω$ is the initial angular velocity of the brake disc. The brake pad and brake disc are simplified into the spring-damping system, the mechanical model of the disc brake is shown in Figure 2.

Obviously, there are two DOFs in the mechanical system, namely, the displacement $y_{b}$ of the brake pad along the Y-axis and the rotation angle $θ_{r}$ of the brake disc around the Z-axis. In Figure 2, $m_{b}$ is the mass of the brake pad, $c_{b}$ is the equivalent damping of the brake pad, $k_{b}$ is the equivalent stiffness of the brake pad, $I_{r}$ is the moment of inertia of the brake disc with respect to the Z-axis, $c_{r}$ is the equivalent angular damping of the brake disc, $k_{r}$ is the torsional stiffness of the brake disc, and r_b is the distance from the COM of the brake disc to the COM of the brake pad.

Hence, the energy of the mechanical system can be calculated by equations (1) to (3).

T_{s} = \frac{1}{2} m_{b} {\dot{y}}_{b}^{2} + \frac{1}{2} I_{r} {\dot{θ}}_{r}^{2}

(1)where

T_{s}

is the kinetic energy.

V_{s} = \frac{1}{2} k_{b} y_{b}^{2} + \frac{1}{2} k_{r} θ_{r}^{2}

(2)where

V_{s}

is the potential energy.

D_{s} = \frac{1}{2} c_{b} {\dot{y}}_{b}^{2} + \frac{1}{2} c_{r} {\dot{θ}}_{r}^{2}

(3)where

D_{s}

is the dissipative function.

The generalized forces with respect to the displacement $y_{b}$ and rotation angle $θ_{r}$ are

Q_{b} = F_{b}

(4)

Q_{r} = T_{0} - T_{r}

(5)where

F_{b}

is the frictional force acting on the brake pad;

T_{0}

is the driving torque acting on the brake disc; and

T_{r}

is the frictional torque acting on the brake disc,

T_{r} = F_{b} r_{b}

For equations (1) ∼ (5), with the help of the second Lagrange equation, the mathematical equations of the brake system are derived as

m_{b} {\ddot{y}}_{b} + c_{b} {\dot{y}}_{b} + k_{b} y_{b} = F_{b}

(6)

I_{r} {\ddot{θ}}_{r} + c_{r} {\dot{θ}}_{r} + k_{r} θ_{r} = T_{0} - T_{r}

(7)

Based on the sensitivity analysis (Meehan and Leslie, 2021), the stiffness of the brake disc is a sensitive parameter for brake stability. Considering the model simplification and measurement error in the brake system, there is uncertainty in the torsional stiffness of the brake disc. Therefore, a Gaussian random process $ξ (t)$ with noise intensity 2D is used to simulate the uncertainty of the torsional stiffness of the brake disc, that is, $ρ ξ (t)$ . Then, equation (7) can be written as

I_{r} {\ddot{θ}}_{r} + c_{r} {\dot{θ}}_{r} + [k_{r} + ρ ξ (t)] θ_{r} = T_{0} - T_{r}

(8)where

ρ

is the uncertainty coefficient of the torsional stiffness of the brake disc.

Friction model

The friction characteristic between the brake pad and brake disc is an important nonlinear factor in the brake system, which is closely related to the brake chatter. To describe the negative friction-velocity gradient, the mathematical expression of the Stribeck model is used.

F_{b} = μ (v_{r}) F_{N}

(9)where

F_{N}

is the braking pressure acting on the brake pad and

μ (v_{r})

is a function of friction coefficient and relative velocity, which can be calculated by

μ (v_{r}) = [μ_{k} + (μ_{s} - μ_{k}) e^{- a | v_{r} |}] sgn (v_{r})

(10)where

μ_{k}

is the dynamic friction coefficient,

μ_{s}

is the static friction coefficient, a is the attenuation coefficient, which determines the degree of negative gradient, their values can be found in Table 1,

sgn (\cdot)

is a sign function, and

v_{r}

is the relative velocity between the brake pad and brake disc, which can be calculated by

v_{r} = (ω + {\dot{θ}}_{r}) r_{b} - {\dot{y}}_{b}

(11)

Table 1.

Coefficients of the Stribeck model.

$μ_{k}$	$μ_{s}$	a	$d_{1}$ (s/m)	$d_{2}$ (s³/m³)
0.50	0.30	4.6	−0.8598	0.2638

Since the relative velocity between the brake pad and brake disc is always positive, the Stribeck model can be derived as

F_{b} = μ (v_{r}) F_{N} = [μ_{k} + (μ_{s} - μ_{k}) e^{- a v_{r}}] F_{N}

(12)

The Stribeck model is difficult to analytically solve due to its strong nonlinearity. Since the relative velocity between the brake pad and brake disc is usually small (Wei et al., 2019), it is necessary to simplify the Stribeck model by means of polynomials to facilitate the stochastic dynamics analysis.

F_{b} = \{μ_{k} + Δ μ e^{- a ω r_{b}} [1 + d_{1} a ({\dot{θ}}_{r} r_{b} - {\dot{y}}_{b}) + d_{2} a^{3} {({\dot{θ}}_{r} r_{b} - {\dot{y}}_{b})}^{3}]\} F_{N} = β_{0} + β_{1} ({\dot{θ}}_{r} r_{b} - {\dot{y}}_{b}) + β_{3} {({\dot{θ}}_{r} r_{b} - {\dot{y}}_{b})}^{3}

(13)where

Δ μ = μ_{s} - μ_{k}

is the friction coefficient difference,

β_{0} = μ_{k} F_{N} + Δ μ F_{N} e^{- a ω r_{b}}

β_{1} = d_{1} a Δ μ F_{N} e^{- a ω r_{b}}

β_{2} = \frac{a^{2}}{2} Δ μ F_{N} - \frac{a^{3}}{2} Δ μ F_{N} ω r_{b}

β_{3} = d_{2} a^{3} Δ μ F_{N} e^{- a ω r_{b}}

, and d₁ and d₂ are fitting coefficients, their specific values are listed in Table 1.

Stability analysis

The singular boundary theory of the stochastic diffusion process is an effective method to explore the stochastic stability of the nonlinear system. Hence, the stochastic averaging method is used to solve the one-dimensional stochastic diffusion process, the stochastic stability of the brake system is analyzed for different interesting parameters (Huan and Zhu, 2009; Ling et al., 2015; Zhu and Cai, 2017; Zhu et al., 2002).

Itô stochastic differential equation

Due to the presence of the constant input in the brake system, the static solution of non-homogeneous differential equations (6) and (8) can be solved by setting $\ddot{y_{b}}$ = 0, $\ddot{θ_{r}}$ = 0, ${\dot{y}}_{b}$ = 0, ${\dot{θ}}_{r}$ = 0.

{\begin{cases} y_{0} = \frac{β_{0}}{k_{b}} \\ θ_{0} = \frac{- β_{0} r_{b}}{k_{r}} \end{cases}

(14)

Let the state variables ${\begin{cases} q_{1} = y_{b} - x_{0} \\ q_{2} = (θ - θ_{0}) r_{b} - (y_{b} - x_{0}) \end{cases}$ and ${\begin{cases} p_{1} = \dot{q_{1}} \\ p_{2} = \dot{q_{2}} \end{cases}$ , equations (6) and (8) can be transformed into

{\begin{cases} \dot{q_{1}} = p_{1} \\ \dot{p_{1}} = - k_{11} q_{1} - k_{12} q_{2} - k_{13} p_{1} - k_{14} p_{2} - δ_{1} p_{1}^{3} + \frac{ρ ξ (t)}{I_{r}} q_{1} \\ \dot{q_{2}} = p_{2} \\ \dot{p_{2}} = - k_{21} q_{2} - k_{22} p_{1} - k_{23} p_{2} - δ_{2} p_{1}^{3} \end{cases}

(15)where

k_{11} = \frac{k_{r}}{I_{r}}

k_{12} = \frac{k_{r}}{I_{r}} - \frac{k_{b}}{m_{b}}

k_{13} = \frac{β_{1}}{m_{b}} + \frac{β_{1} r_{b}^{2} + c_{r}}{I_{r}}

k_{14} = \frac{c_{r}}{I_{r}} - \frac{c_{b}}{m_{b}}

k_{21} = \frac{k_{b}}{m_{b}}

k_{22} = - \frac{β_{1}}{m_{b}}

k_{23} = \frac{c_{b}}{m_{b}}

δ_{1} = \frac{β_{3}}{m_{b}} + \frac{β_{3} r_{b}^{2}}{I_{r}}

δ_{2} = - \frac{β_{3}}{m_{b}}

Based on the kinetic and potential energy of the brake system, the Hamiltonian function of the brake system is set to

H = \frac{1}{2} p_{1}^{2} + \frac{1}{2} p_{2}^{2} + U

(16)where

U = \frac{1}{2} k_{11} q_{1}^{2} + k_{12} q_{1} q_{2} + \frac{1}{2} k_{21} q_{2}^{2}

With the help of the stochastic averaging method (Gu et al., 2021; Hu et al., 2023), the equivalent It $\hat{o}$ stochastic differential equation of equation (15) can be transformed as

{\begin{cases} d q_{1} = \frac{\partial H}{\partial p_{1}} d t \\ d p_{1} = - [\frac{\partial H}{\partial q_{1}} + m_{11} \frac{\partial H}{\partial p_{1}} + m_{12} \frac{\partial H}{\partial p_{2}}] d t + σ_{11} d B (t) \\ d q_{2} = \frac{\partial H}{\partial p_{2}} d t \\ d p_{2} = - [\frac{\partial H}{\partial q_{2}} + m_{21} \frac{\partial H}{\partial p_{1}} + m_{22} \frac{\partial H}{\partial p_{2}}] d t \end{cases}

(17)where

m_{11} = k_{13} + δ_{1} p_{1}^{2}

m_{12} = k_{14}

σ_{11} = \frac{- ρ q_{1}}{I_{r}} \sqrt{2 D}

m_{21} = k_{22} + δ_{2} p_{1}^{2}

m_{22} = k_{23} - \frac{k_{12} q_{1}}{p_{2}}

B (t)

is the Weiner process.

The stochastic brake system belongs to a quasi non integrable Hamiltonian system. Based on the properties of the quasi non integrable Hamiltonian system, equation (17) is transformed into the Itô stochastic differential equation of the one-dimensional energy process

d H = m (H) d t + σ (H) d B (t)

(18)where

m (H)

is the averaged drift coefficient and

σ (H)

is the averaged diffusion coefficient, which can be solved by

{\begin{cases} m (H) = \frac{1}{T (H)} \int_{Ω_{1}} (- m_{11} p_{1} - m_{12} p_{2} - m_{21} p_{2} - \frac{m_{22} p_{2}^{2}}{p_{1}} + \frac{σ_{11}^{2}}{2 p_{1}} + \frac{σ_{12}^{2}}{2 p_{1}}) d q_{1} d q_{2} d p_{2} \\ σ^{2} (H) = \frac{1}{T (H)} \int_{Ω_{1}} (σ_{11}^{2} p_{1} + σ_{12}^{2} p_{1}) d q_{1} d q_{2} d p_{2} \\ T (H) = \int_{Ω_{1}} \frac{1}{p_{1}} d q_{1} d q_{2} d p_{2} \end{cases}

(19)where

Ω_{1} = {(q_{1}, q_{2}, p_{2}) | \frac{1}{2} p_{2}^{2} + \frac{1}{2} k_{11} q_{1}^{2} + k_{12} q_{1} q_{2} + \frac{1}{2} k_{21} q_{2}^{2} \leq H}

According to equation (16), we can obtain

p_{1} = \sqrt{2 (H - U) - p_{2}^{2}}

(20)

Then, equation (21) can be obtained by performing the coordinate transformation on equation (20).

{\begin{cases} p_{1} = \sqrt{2 (H - U)} \cos (φ) \\ p_{2} = \sqrt{2 (H - U)} \sin (φ) \\ d p_{2} = \sqrt{2 (H - U)} \cos (φ) d φ \end{cases}

(21)where

- \frac{π}{2} \leq φ \leq \frac{π}{2}

By substituting equation (21) into equation (19) and integrating the state variable $p_{2}$ , we can obtain

{\begin{cases} m (H) = \frac{1}{T (H)} \int_{Ω_{2}} [(- (k_{13} π + k_{23} π) (H - \frac{1}{2} k_{11} q_{1}^{2} - k_{12} q_{1} q_{2} - \frac{1}{2} k_{21} q_{2}^{2}) \\ - \frac{3 π δ_{1}}{2} {(H - \frac{1}{2} k_{11} q_{1}^{2} - k_{12} q_{1} q_{2} - \frac{1}{2} k_{21} q_{2}^{2})}^{2} + \frac{π}{I_{r}^{2}} ρ^{2} D q_{1}^{2}] d q_{1} d q_{2} \\ σ^{2} (H) = \frac{1}{T (H)} \int_{Ω_{2}} [π (H - \frac{1}{2} k_{11} q_{1}^{2} - k_{12} q_{1} q_{2} - \frac{1}{2} k_{21} q_{2}^{2}) σ_{11}^{2}] d q_{1} d q_{2} \\ T (H) = \int_{Ω_{2}} π d q_{1} d q_{2} \end{cases}

(22)where

Ω_{2} = {(q_{1}, q_{2}) | \frac{1}{2} k_{11} q_{1}^{2} + k_{12} q_{1} q_{2} + \frac{1}{2} k_{21} q_{2}^{2} \leq H}

To solve the double integral mentioned above, the variables $q_{1}$ and $q_{2}$ are transformed into polar coordinates.

{\begin{cases} q_{1} = \frac{r \cos (α)}{\sqrt{k_{11}}} \\ q_{2} = \frac{r \sin (α)}{\sqrt{k_{21}}} \\ d q_{1} d q_{2} = \frac{r}{\sqrt{k_{11} k_{21}}} d r d α \end{cases}

(23)

Equation (24) is derived by substituting equation (23) into equation (22).

{\begin{cases} m (H) = \frac{1}{T (H) \sqrt{k_{11} k_{21}}} \int_{0}^{2 π} d α \int_{0}^{R} [- (k_{13} π + k_{23} π) (H - \frac{1}{2} r^{2} \\ - \frac{\cos (α) \sin (α) k_{12}}{\sqrt{k_{11}} \sqrt{k_{21}}} r^{2}) - \frac{3 π δ_{1}}{2} {(H - \frac{1}{2} r^{2} - \frac{\cos (α) \sin (α) k_{12}}{\sqrt{k_{11}} \sqrt{k_{21}}} r^{2})}^{2} \\ + \frac{π ρ^{2} D r^{2} \cos^{2} (α)}{I_{r}^{2} k_{11}}] r d r \\ σ^{2} (H) = \frac{1}{T (H) \sqrt{k_{11} k_{21}}} \int_{0}^{2 π} d α \int_{0}^{R} [\frac{2 D π ρ^{2} r^{2} \cos^{2} (α)}{I_{r}^{2} k_{11}} (H - \frac{1}{2} r^{2} \\ - \frac{\cos (α) \sin (α) k_{12}}{\sqrt{k_{11}} \sqrt{k_{21}}} r^{2})] r d r \\ T (H) = \int_{0}^{2 π} d α \int_{0}^{R} \frac{π r}{2 \sqrt{k_{11} k_{21}}} d r \end{cases}

(24)

Then, equation (25) can be calculated by integrating equation (24).

{\begin{cases} m (H) = \frac{1}{T (H) \sqrt{k_{11} k_{21}}} [- 2 π (k_{13} π + k_{23} π) (\frac{H}{2} R^{2} - \frac{1}{8} R^{4}) - \\ \frac{3 π δ_{1}}{2} (π H^{2} R^{2} + \frac{π}{12} R^{6} - \frac{π}{2} H R^{4} + \frac{π k_{12}^{2}}{24 k_{11} k_{21}} R^{6}) + \frac{π^{2} ρ^{2} D R^{4}}{I_{r}^{2} 4 k_{11}}] \\ σ^{2} (H) = \frac{1}{T (H) \sqrt{k_{11} k_{21}}} [\frac{D ρ^{2} π}{I_{r}^{2}} (\frac{H R^{4} π}{2 k_{11}} - \frac{R^{6} π}{6 k_{11}})] \\ T (H) = \frac{π^{2} R^{2}}{\sqrt{k_{11} k_{21}}} \end{cases}

(25)where R is the solution of

H = \frac{1}{2} R^{2} + \frac{k_{12} R^{2} \cos (α) \sin (α)}{\sqrt{k_{11} k_{21}}}

The expected value of R with respect to $α$ is solved as

H = E [\frac{1}{2} R^{2} + \frac{k_{12} R^{2} \cos (α) \sin (α)}{\sqrt{k_{11} k_{21}}}] = \frac{1}{2} R^{2}

(26)

Since $T (H) = \frac{π^{2} R^{2}}{\sqrt{k_{11} k_{21}}} = \frac{2 π^{2} H}{\sqrt{k_{11} k_{21}}}$ , the averaged drift coefficient and averaged diffusion coefficient are obtained.

{\begin{cases} m (H) = [\frac{ρ^{2} D}{2 k_{11} I_{r}^{2}} - \frac{(k_{13} + k_{23})}{2}] H - (\frac{δ_{1}}{2} + \frac{k_{12}^{2} δ_{1}}{4 k_{11} k_{21}}) H^{2} \\ σ^{2} (H) = \frac{D ρ^{2}}{3 I_{r}^{2} k_{11}} H^{2} \end{cases}

(27)

Boundary analysis

The stochastic stability of the brake system depends on the boundary type of the one-dimensional energy process. If $x_{s}$ satisfies the diffusion coefficient $σ (x_{s}) = 0$ , $x_{s}$ is called the first type singular boundary. If $x_{s}$ satisfies the diffusion coefficient $m (x_{s}) = \infty$ , $x_{s}$ is called a second type singular boundary. The boundary type of the one-dimensional energy process can be distinguished by the drift index, diffusion index, and character value (Zhu and Cai, 2017).

For the one-dimensional diffusion process H(t), it can be divided into two boundaries, namely, the left boundary (i.e., $H \to 0$ ) and the right boundary (i.e., $H \to \infty$ ). When the left boundary (i.e., $H \to 0$ ), the averaged drift coefficient and averaged diffusion coefficient are related to low-order terms, equation (27) can be rewritten as

{\begin{cases} m (H) = [\frac{ρ^{2} D}{2 k_{11} I_{r}^{2}} - \frac{(k_{13} + k_{23})}{2}] H + O (H) \\ σ^{2} (H) = \frac{D ρ^{2}}{3 I_{r}^{2} k_{11}} H^{2} + O (H^{2}) \end{cases}

(28)

Then, the drift index $β_{L}$ , diffusion index $α_{L}$ , and character value $c_{L}$ of the left boundary is solved.

{\begin{cases} σ^{2} (H) = O (H^{α_{L}}) = O (H^{2}) \\ m (H) = O (H^{β_{L}}) = O (H^{1}) \\ c_{L} = \lim_{H \to 0^{+}} \frac{2 m (H) H^{α_{L} - β_{L}}}{σ^{2} (H)} = \frac{3 ρ^{2} D - 3 I_{r}^{2} k_{11} (k_{13} + k_{23})}{ρ^{2} D} \end{cases}

(29)

According to equation (29), we can obtain $α_{L} = β_{L} + 1$ , the boundary type depends on the value of the character value $c_{L}$ . When the character value $c_{L} < 1$ , the trajectory of the system solution converges to the left boundary, causing the system energy to gradually decrease. When the character value $c_{L} > 1$ , the left boundary is a repulsive natural boundary, the trajectory of the system solution enters the system from the left boundary, causing the system energy to continuously increase.

In addition, when the right boundary (i.e., $H \to \infty$ ), the averaged drift coefficient and diffusion coefficient are related to higher-order terms, equation (27) is rewritten as

{\begin{cases} m (H) = - (\frac{δ_{1}}{2} + \frac{k_{12}^{2} δ_{1}}{4 k_{11} k_{21}}) H^{2} + O (H^{2}) \\ σ^{2} (H) = \frac{D ρ^{2}}{3 I_{r}^{2} k_{11}} H^{2} + O (H^{2}) \end{cases}

(30)

Then, the drift index $β_{R}$ , diffusion index $α_{R}$ , and character value $c_{R}$ of the right boundary can be solved.

{\begin{cases} σ^{2} (H) = O (H^{α_{R}}) = O (H^{2}) \\ m (H) = O (H^{β_{R}}) = O (H^{2}) \\ c_{R} = - \lim_{H \to \infty} \frac{2 m (H) H^{α_{R} - β_{R}}}{σ^{2} (H)} \end{cases}

(31)

The parameters of the brake system are listed in Table 2; according to equation (31), we can obtain

β_{R} > α_{R} - 1

m (+ \infty) < 0

, and

β_{R} \geq 1

, the right boundary is the entry boundary (Zhu and Cai, 2017), the trajectory of the system solution enters the interior from the right boundary, and the system energy will not diffuse out from the right boundary.

Table 2.

Parameters of the brake system.

Parameters	Symbols	Values (unit)
Moment of inertia of the brake disc with respect to the Z-axis $I_{r}$	$I_{r}$	0.65 (kg·m²)
Mass of the brake pad $m_{b}$	$m_{b}$	1.50 (kg)
Torsional stiffness of the brake disc $k_{r}$	$k_{r}$	300.00 (kN·m/rad)
Equivalent stiffness of the brake pad $k_{b}$	$k_{b}$	300.00 (kN/m)
Equivalent angular damping of the brake disc $c_{r}$	$c_{r}$	1.00 (N·m·s/rad)
Equivalent damping of brake pad $c_{b}$	$c_{b}$	1.00 (N·s/m)
Distance from the COM of the brake disc to the COM of the brake pad $r_{b}$	$r_{b}$	0.12 (m)
Uncertainty coefficient of the torsional stiffness of the brake disc $ρ$	$ρ$	1.72 (kN·m/rad)
Noise intensity D	D	0.15

Obviously, the property of the left boundary of the one-dimensional energy process depends on the character value $α_{L}$ , the right boundary of the one-dimensional energy process is always the entry boundary. Hence, when the character value $c_{L} < 1$ , the system energy converges to zero, the stationary solution of the one-dimensional diffusion process is stable. When the character value $c_{L} > 1$ , the left boundary is the repulsive natural boundary, the system energy will converge to a certain steady-state energy, at this point, the stationary solution of the one-dimensional diffusion process is unstable.

Discussions

First of all, the initial angular velocity of the brake disc is selected as an interesting parameter, the character values of the brake system for different initial angular velocities can be solved and plotted in Figure 3. Obviously, as the initial angular velocity increases, the character value decreases. When the initial angular velocity $ω = 11.97 rad / s$ , the character value $c_{L} = 1$ . Taking the rolling radius of the wheel R = 0.35 m as an example, the critical angular velocity of the brake disc can be converted into the critical vehicle speed v = ωR ≈ 15.08 km/h, which is consistent with the occurrence speed of the brake chatter in engineering (Wei et al., 2019). When the initial angular velocity further increases, the stationary solution of the one-dimensional energy process will change from unstable to stable, the brake chatter can be suppressed.

To analyze the influence of parameters, such as the noise intensity D, brake pressure F_N, and friction coefficient difference $Δ μ$ , on the stochastic stability of the brake system, the character values of the brake system for different interesting parameters are solved, as shown in Figure 4. Clearly, the character value increases with the increase of the noise intensity, brake pressure, and friction coefficient difference, the stochastic stability of the brake system is weakened. Hence, in order to improve the stochastic stability of the brake system for different initial angular velocities, the braking pressure should be reduced while ensuring the brake safety. Meanwhile, the difference between the static friction coefficient and the dynamic friction coefficient of the brake system should also be reduced to enhance the stochastic stability. Therefore, effective measures such as using non-asbestos organic friction material, optimizing the brake disc structure to ensure heat dissipation, and implementing standardized running-in of the brake disc and brake pad, can be used to reduce the friction coefficient difference.

Stochastic bifurcation

As depicted in Figure 4, when the system parameters cross the character value $c_{L} = 1$ , the stochastic stability of the brake system will change, which is closely related to the stochastic bifurcation. Therefore, the bifurcation characteristics of the brake system with the stochastic parameter perturbation are further explored to reveal the dynamic mechanism of the brake chatter.

Bifurcation characteristic

To analyze the stochastic bifurcation characteristics, it is necessary to solve the stationary probability density of the one-dimensional energy process of the brake system, which can be obtained from the Fokker Planck Kolmogorov (FPK) equation (32).

\frac{\partial P_{s t}}{\partial t} = - \frac{\partial [m (H) P_{s t}]}{\partial H} + \frac{1}{2} \frac{\partial^{2} [σ^{2} (H) P_{s t}]}{\partial H^{2}}

(32)where

P_{s t}

is the stationary probability density.

By integrating the FPK equation (32), the stationary probability density of the brake system can be obtained as

P_{s t} (H) = C \exp [- \int_{0}^{H} \frac{\frac{d σ^{2} (u)}{d u} - 2 m (u)}{σ^{2} (u)} d u] = C H^{η} \exp (κ H)

(33)where C is the normalization coefficient, which can be calculated from

C^{- 1} = \int_{0}^{\infty} H^{η} \exp (κ H) d H

η = 1 - \frac{3 I_{r}^{2} k_{11} (k_{13} + k_{23})}{ρ^{2} D}

κ = - (\frac{3 I_{r}^{2} k_{11} δ_{1}}{ρ^{2} D} + \frac{3 I_{r}^{2} k_{12}^{2} δ_{1}}{2 ρ^{2} D k_{21}})

Since the stochastic stability of the brake system depends on the property of the left boundary of the one-dimensional energy process, when the system energy $H \to 0$ , substituting equation (28) into equation (33), equation (33) can be rewritten as

\begin{array}{l} P_{s t} (H) = O {\exp [- \int_{0}^{H} (\frac{α_{L} u^{α_{L} - 1}}{u^{α_{L}}} - c_{L} u^{β_{L} - α_{L}}) d u]} \\ = O [H^{- α_{L}} \exp (c_{L} \ln H)] = O (H^{ζ}) \end{array}

(34)where

ζ = β_{L} - α_{L}

can be defined as the stochastic bifurcation coefficient.

When the stochastic bifurcation coefficient $ζ < - 1$ , the stationary probability density $P_{s t} (H)$ is a δ-function, the stationary solution of the brake system is absolutely stable. When the stochastic bifurcation coefficient $- 1 < ζ < 0$ , the stationary probability density $P_{s t} (H)$ is a decreasing function, its maximum value is at the system energy H = 0, and the stationary solution of the system is probabilistically stable. When the stochastic bifurcation coefficient $ζ > 0$ , the maximum value of the stationary probability density $P_{s t} (H)$ is at the system energy H ≠ 0, the stationary solution of the system is probabilistically unstable. Hence, the dynamic bifurcation (D bifurcation) occurs at the stochastic bifurcation coefficient $ζ = - 1$ , the phenomenological bifurcation (P bifurcation) occurs at the stochastic bifurcation coefficient $ζ = 0$ (Zhu and Cai, 2017).

Figure 5 shows the stochastic bifurcation coefficient of the brake system for different initial angular velocities. Obviously, there are the stochastic D bifurcation and stochastic P bifurcation in the brake system, with critical initial angular velocities of $ω = 11.97 rad / s$ and $ω = 10.46 rad / s$ , respectively. When the initial angular velocities of $ω > 11.97 rad / s$ , the stationary solution of the system is absolutely stable, which is consistent with the results shown in Figure 3, verifying the analysis results of the stochastic stability. When the initial angular velocities $10.46 ra d / s < ω < 11.97 rad / s$ , the stationary solution of the system is probabilistically stable. However, when the initial angular velocities $ω < 10.46 rad / s$ , the stationary solution is probabilistically unstable.

In addition, the stochastic bifurcation diagrams of the brake system for different initial angular velocities and interesting parameters are shown in Figures 6 –8, the blue dashed line is the D bifurcation boundary, the red dashed line is the P bifurcation boundary. According to the stochastic bifurcation coefficient, there are three regions in the bifurcation diagram: the absolutely stable region (i.e., $ζ < - 1$ ), the probability stable region (i.e., $- 1 < ζ < 0$ ), and the probability unstable region (i.e., $0 < ζ$ ).

As shown in Figure 6, when the initial angular velocity is fixed and the noise intensity increases, the brake system firstly undergoes the D bifurcation, the stationary solution changes from absolutely stable to probabilistically stable. As the noise intensity further increases, the brake system will undergo the P bifurcation, the stationary solution becomes probabilistically unstable. Similarly, it can be observed from Figures 7 and 8 that as the brake pressure and friction coefficient difference increases, the brake system firstly undergoes the D bifurcation, then the P bifurcation, causing the stationary solution to change from absolutely stable to probabilistically stable, and finally to probabilistically unstable. Hence, the noise intensity, brake pressure, and friction coefficient difference all need to be reduced to improve the stochastic stability of the brake system.

Numerical verifications

To verify the analysis results of the stochastic stability and bifurcation characteristics, substituting equation (16) into equation (33), the joint probability density function of the brake system can be obtained.

P_{s t} (p, q) = C {(\frac{1}{2} p_{1}^{2} + \frac{1}{2} p_{2}^{2} + \frac{1}{2} k_{11} q_{1}^{2} + k_{12} q_{1} q_{2} + \frac{1}{2} k_{21} q_{2}^{2})}^{η} \exp [κ (\frac{1}{2} p_{1}^{2} + \frac{1}{2} p_{2}^{2} + \frac{1}{2} k_{11} q_{1}^{2} + k_{12} q_{1} q_{2} + \frac{1}{2} k_{21} q_{2}^{2})]

(35)

The stationary probability density of the brake pad is selected as an example for the numerical analysis, the stationary probability density with the respect to state variables q₂ and p₂ is obtained by integrating equation (35).

P_{s t} (p_{2}, q_{2}) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} P_{s t} (p, q) d p_{1} d q_{1}

(36)

Then, the cases of different bifurcation regions in Figure 6(b) are solved by using the integral2 algorithm in Matlab software, as shown in Figure 9 the parameter settings of integral2 algorithm are as follows: the absolute error is 10^–6, the integral boundary of p₁ is [–0.6 m/s, 0.6 m/s], the integral boundary of q₁ is [–0.002 m, 0.002 m]. For case A₁ in the absolutely stable region, the stationary probability density is a δ-function, as shown in Figure 9(a). The value of the stationary probability density tends to infinity at the origin, while the value of the stationary probability density at other coordinates is zero. Hence, the system response always converges to a stable equilibrium point. For case A₂ in the probability stable region, as shown in Figure 9(b), the value of the stationary probability density at the origin is no longer infinite, the value of the stationary probability density values at positions away from the origin decreases monotonically. This indicates that the probability of the system response converging to the equilibrium point is higher, while the probability of converging to the non-equilibrium point is lower. Note that for case A₃ in the probability unstable region, as shown in Figure 9(c), the central attraction domain of the brake system becomes the attraction domain of the probability limit cycle, the system response no longer converges to the equilibrium point but converges to the probability stable limit cycle.

Moreover, Figures 10 and 11 show different cases from Figures 7 and 8. Similarly, when the brake system undergoes the D bifurcation, the system response converges to the equilibrium point with a high probability. However, when the brake system undergoes the P bifurcation, the system response converges to a limit cycle with a certain probability. It is evident that the numerical results of the stationary probability density agree well with bifurcation analysis, which verifies the accuracy of theoretical analysis.

Conclusions

This paper focuses on the dynamic response and stability analysis of the vehicle disc brake with the stochastic parameter perturbation. A Gaussian random process is used to simulate the uncertainty of the torsional stiffness of the brake disc; the stochastic dynamic model of the brake chatter is established. With the help of the stochastic averaging method, the boundary type of the one-dimensional energy process is discussed. The system stability depends on the character value of the left boundary of the one-dimensional energy process. When the character value is less than one, the system energy converges to zero, the trajectory of the system solution converges to the left boundary. When the character value is greater than one, the left boundary is a repulsive natural boundary, the system energy will converge to a certain steady-state energy, at this point, the trajectory of the system solution enters the system from the left boundary. When the initial angular velocity is less than the critical value, the stationary solution of the one-dimensional energy process will change from stable to unstable, that is, the brake chatter occurs.

Furthermore, the bifurcation characteristics of the brake system with the stochastic parameter perturbation are analyzed to reveal the chatter mechanism. The D bifurcation and P bifurcation are analytically proven by means of the FPK equation. As the noise intensity, brake pressure and friction coefficient difference increases, the brake system will firstly undergo the D bifurcation, then the P bifurcation, causing the stationary solution to change from absolutely stable to probabilistically stable, and finally to probabilistically unstable. Hence, the noise intensity, brake pressure, and friction coefficient difference all need to be reduced to improve the stochastic stability of the brake system. Numerical results of the stationary probability density of the brake system for different initial angular velocities and interesting parameters are in good agreement with the analytical results, which helps to improve the system stability and suppress the brake chatter of the vehicle.

Footnotes

ORCID iDs

Heng Wei

Jialin Li

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China [Grant number 52402477, 52275100] and the Science Fund for Distinguished Young Scholars of Anhui Province [Grant number 2408085J034].

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Akif

Osman

(2023) Investigation of the pad stiffness effects on the initiation of brake squeal phenomenon. Journal of Theoretical and Applied Mechanics 61(2): 219–232.

Andreaus

Casini

(2001) Dynamics of friction oscillators excited by a moving base and/or driving force. Journal of Sound and Vibration 245(4): 685–699.

Cantoni

Cesarini

Mastinu

, et al. (2009) Brake comfort-A review. Vehicle System Dynamics 47(8): 901–947.

Chatelet

Michon

Manin

, et al. (2008) Stick/slip phenomena in dynamics: choice of contact model. Numerical predictions & experiments. Mechanism and Machine Theory 43: 1211–1224.

Chen

(2014) Stability and response of a self-amplified braking system under velocity-dependent actuation force. Nonlinear Dynamics 78(4): 2459–2477.

Chen

Zhang

Neighborgall

, et al. (2025) Effect of drum and disc brakes on roll stability of double-trailer trucks during braking-in-turn manoeuvers. Vehicle System Dynamics 1–16. Available at: https://doi.org/10.1080/00423114.2025.2460493

Crowther

Singh

(2008) Identification and quantification of stick-slip induced brake groan events using experimental and analytical investigations. Noise Control Engineering Journal 56(4): 235–255.

Zhang

Wei

, et al. (2024) Dynamic characteristics of three-body braking system considering tire-road friction. Journal of Computational and Nonlinear Dynamics 19(1): 011006.

Fischer

Pürscher

Huemer-Kals

, et al. (2021) Characterization of brake creep groan vibrations. SAE International Journal of Advances and Current Practices in Mobility 3(2): 1049–1058.

10.

Grzes

Kuciej

(2025) Conversion of images of 3D friction maps to study the coupling between coefficient of friction, velocity and contact temperature of the disc brake. Mechanical Systems and Signal Processing 223: 111853.

11.

Jia

Deng

, et al. (2021) Stochastic response of nonlinear viscoelastic systems with time-delayed feedback control force and bounded noise excitation. International Journal of Structural Stability and Dynamics 21(13): 2150181.

12.

Zhao

Deng

, et al. (2023) Approximate analytical response of nonlinear functionally graded beams subjected to harmonic and random excitations. International Journal of Non-Linear Mechanics 148(Suppl C): 104269.

13.

Zhang

Mughal

, et al. (2024) Stochastic responses of nonlinear inclined cables with an attached damper and random excitations. Nonlinear Dynamics 112: 15969–15986.

14.

Hetzler

Schwarzer

Seemann

(2007) Analytical investigation of steady-state stability and Hopf-bifurcations occurring in sliding friction oscillators with application to low-frequency disc brake noise. Communications in Nonlinear Science and Numerical Simulation 12(1): 83–99.

15.

Zhang

, et al. (2023) Stochastic stabilization of multi-degree-of-freedom nonlinear random-time-delay controlled systems. International Journal of Robust and Nonlinear Control 33(3): 2288–2303.

16.

Zeng

Zhou

, et al. (2025) Stochastic dynamics in power systems excited by discrete-continuous random disturbances. Probabilistic Engineering Mechanics 82: 103858.

17.

Huan

Zhu

(2009) Stochastic optimal control of quasi integrable Hamiltonian systems subject to actuator saturation. Journal of Vibration and Control 15(1): 85–99.

18.

İbrahim Cana

Özgünb

Özgün Cema

(2025) Modelling and control of electromechanical disc brake for heavy-duty vehicles. Proceedings of the Institution of Mechanical Engineers - Part D: Journal of Automobile Engineering 239(5): 1717–1732.

19.

Kang

Krousgrill

Sadeghi

(2009) Comprehensive stability analysis of disc brake vibrations including gyroscopic, negative friction slope and mode-coupling mechanisms. Journal of Sound and Vibration 324(1-2): 387–407.

20.

Chen

, et al. (2025) A review of braking performance and dynamic performance of brake discs. Journal of Tribology 147(12): 120801.

21.

Ling

Jin

Wang

, et al. (2015) Stochastic stabilization of quasi-generalized Hamiltonian systems. Journal of Vibration and Control 21(1): 47–59.

22.

Zhao

(2022) Observer-based friction force estimation for brake system under low-frequency vibration. Journal of Vibration and Control 28(13-14): 1553–1564.

23.

Lü

Shangguan

(2017) A unified approach for squeal instability analysis of disc brakes with two types of random-fuzzy uncertainties. Mechanical Systems and Signal Processing 93(1): 281–298.

24.

Manish

Don

Chu

, et al. (2005) Noise and vibration analysis of a disc-brake system using a stick-slip friction model involving coupling stiffness. Journal of Sound and Vibration 282(3-5): 1273–1284.

25.

Marques

Woliński

Wojtyra

, et al. (2021) An investigation of a novel lugre-based friction force model. Mechanism and Machine Theory 166: 104493.

26.

Meehan

Leslie

(2021) On the mechanisms, growth, amplitude and mitigation of brake squeal noise. Mechanical Systems and Signal Processing 152: 107469.

27.

Rhee

Jacko

Tsang

PHS

(1991) The role of friction film in friction, wear, and noise of automotive brakes. Wear 146(1): 89–97.

28.

Sarrouy

Dessombz

Sinou

(2013) Stochastic study of a non-linear self-excited system with friction. European Journal of Mechanics - A: Solids 40(0): 1–10.

29.

Shin

Brennan

, et al. (2002) Analysis of disc brake noise using a two-degree-of freedom model. Journal of Sound and Vibration 254(5): 837–848.

30.

Song

Zhang

Zhou

, et al. (2025) Analytical study and optimisation of wheel-rail wear on curved track under traction/braking conditions. Vehicle System Dynamics 1–19. Available at: https://doi.org/10.1080/00423114.2025.2484459

31.

Thomsen

Fidlin

(2003) Analytical approximations for stick-slip vibration amplitudes. International Journal of Non-Linear Mechanics 38(3): 389–403.

32.

Wang

Lei

(2024) Trend analysis of rail corrugation on the metro based on wheel-rail stick-slip characteristics. Vehicle System Dynamics 62(1): 147–176.

33.

Wang

, et al. (2024) A novel vehicle-track coupled dynamics model with disc brake systems: modelling and validation. Vehicle System Dynamics 62(4): 789–812.

34.

Wei

Song

Nan

, et al. (2019) Analysis of the stick-slip vibration of a new brake pad with double-layer structure in automobile brake system. Mechanical Systems and Signal Processing 118(0): 305–316.

35.

Zeng

Luo

(2011) Non-linear analysis of disc brake-induced vibrations for railway vehicles. Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail and Rapid Transit 225(1): 48–56.

36.

Zhang

Zhou

, et al. (2024) Frictional vibration analysis of train braking system considering wheel-rail attachment and multi-body friction. International Journal of Non-Linear Mechanics 162: 104715.

37.

Zhao

Zhu

, et al. (2024) Dynamic characteristics and sensitivity analysis of a nonlinear vehicle suspension system with stochastic uncertainties. Nonlinear Dynamics 112(24): 21605–21626.

38.

Zhu

Cai

(2017) Introduction to Stochastic Dynamics. Beijing: Science Press.

39.

Zhu

Huang

Suzuki

(2002) Stochastic averaging and Lyapunov exponent of quasi partially integrable Hamiltonian systems. International Journal of Non-Linear Mechanics 37(3): 419–437.