汽车制动器系统稳定性响应的不确定性 和相关性分析

吕 辉，杨 超，上官文斌，于德介，赵克刚; Lü Hui; YANG Chao; SHANGGUAN Wen-bin; YU De-jie; ZHAO Ke-gang

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

汽车制动器系统稳定性响应的不确定性和相关性分析 PDF

- ORCID：
吕辉 ^1,2
✉
- ORCID：
杨超 ¹
- ORCID：
上官文斌 ¹
- ORCID：
于德介 ²
- ORCID：
赵克刚 ¹
✉

1. 华南理工大学机械与汽车工程学院,广东广州 510641； 2. 湖南大学整车先进设计制造技术全国重点实验室,湖南长沙 410082

中图分类号： U463.51

最近更新：2024-10-11

DOI：10.16385/j.cnki.issn.1004-4523.2024.09.011

摘要

在制动噪声现象中，汽车制动器系统参数不可避免地存在着不确定性和相关性，使得系统响应亦可能同时存在一定的不确定性和相关性，针对该问题开展了制动器系统稳定性响应的不确定性和相关性分析研究。采用多椭球凸模型描述系统参数的不确定性和相关性，以不稳定模态阻尼比表征系统稳定性响应。将多椭球凸模型分别与蒙特卡罗仿真、一阶摄动法和二阶摄动法相结合，提出了三种系统稳定性响应的不确定性分析方法；结合蒙特卡罗仿真和一阶摄动法，分别提出了两种系统不确定响应的相关性分析方法；基于不确定性和相关性分析方法，提出了建立系统响应椭球域的组合方法。通过算例分析验证了方法的有效性。分析结果表明，所提出的方法可有效地求得系统稳定性响应的边界区间、相关系数和椭球域，并且该方法具有较高的计算精度和效率。

关键词

汽车制动器系统; 多椭球凸模型; 不确定性分析; 相关性分析; 制动噪声

引言

汽车制动噪声与制动器的振动稳定性密切相关，众多学者主要通过研究系统稳定性来改善制动噪声问题。受工作条件变化、制造误差、材料老化等因素的影响，制动器系统的诸多参数往往存在不确定性^［

1］。因此，研究考虑不确定性的制动器系统稳定性具有重要的工程意义。

近年来，基于不确定性模型的制动器系统稳定性研究受到了越来越多的关注。张立军等^［

2］通过考虑制动盘参数的随机不确定性，采用稳健性设计方法研究了制动尖叫问题。SARROUY等^{［参考文献 3‑4}3‑4］基于混沌多项式展开提出了一种不确定性分析方法，用于研究线性制动器系统稳定性问题。NOBARI等^{［参考文献 5

百度学术}5］引入代理模型研究制动噪声问题，大大减少了计算时间和成本。黄晓婷等^{［参考文献 6

百度学术}6］基于模糊模型对制动器的参数不确定性进行建模，提出了一种面向制动噪声控制的模糊不确定性分析方法。吕辉等^{［参考文献 7‑8}7‑8］基于证据理论提出了非精确的概率不确定性分析方法，用于分析制动尖叫问题。

可以看出，基于不确定性模型的制动器系统稳定性研究已经取得较多研究成果。但现有研究还存在如下两点不足：第一，上述对于制动器稳定性的研究均将系统参数视为相互独立的不确定变量，没有考虑不确定参数之间的相关性。然而，在工程实际中，系统不确定参数之间往往存在一定的相关性^［

9‑10］。例如系统制动压力和摩擦系数之间往往存在相关性。第二，对于实际的工程结构，不仅输入参数间存在相关性，输出响应之间也可能存在相关性。当一个响应发生变化，可能会影响其他响应的变化，通过响应相关性分析可以获得此类影响规律，进而更好地指导后续的优化设计。然而，现有研究很少涉及系统响应的相关性分析。

针对上述可能存在的问题，本文基于多椭球凸模型开展了制动器稳定性响应的不确定性和相关性分析，有针对性地研究了系统不确定参数存在相关性的情形。首先，采用多椭球凸模型描述系统参数的不确定性和相关性；然后，结合蒙特卡罗法、一阶摄动法和二阶摄动法开展了系统稳定性响应的不确定性分析；接着，结合一阶摄动法求解了系统响应的相关性；最后，通过算例验证了方法的有效性。

1 分析模型

1.1 盘式制动器模型

图1为一汽车浮动钳盘式制动器模型^［

8］，其主要由制动卡钳、定位销、支架、制动块、法兰盘、制动盘等部件组成。其中，制动盘刚性地连接在轮毂上随着车轮转动。制动块则由衬片底板和摩擦材料组成。系统主要通过制动盘和制动块间产生的摩擦阻力矩使得汽车减速或停止运动。制动过程中，如果系统处于不稳定状态，则有可能引发制动噪声^{［参考文献 8

百度学术}8］。

图1 盘式制动器模型

Fig.1 The model of a disc brake

制动器系统的振动方程可表示为^［

1］：

M \ddot{v} + C \dot{v} + (K - J) v = 0

（1）

式中 $M$ ， $C$ 和 $K$ 分别为无摩擦制动器系统的质量、阻尼和刚度矩阵； $J$ 为摩擦接触矩阵； $v$ 为系统振动的广义位移量，由欧拉公式可得：

v = φ e^{λ t}

（2）

式中 $λ$ 为系统特征值； $φ$ 为振型矩阵。

结合式（1）和（2）可得：

(λ^{2} M + λ C + K - J) φ = 0

（3）

系统第 $i$ 阶的特征值可表示为 $λ_{i} = σ_{i} + j ω_{i}$ ，其中， $σ_{i}$ 和 $ω_{i}$ 分别表示特征值的实部和虚部。

系统第i阶特征值对应的模态阻尼比定义为：

ζ_{i} = - \frac{σ_{i}}{π ω_{i}}

（4）

当 $ζ_{i}$ 为负时，对应的复特征值实部为正，此时系统不稳定，有可能引发制动噪声。因此，模态阻尼比可以作为评判制动器系统稳定性的指标。

1.2 多椭球凸模型

制动器的诸多不确定参数之间可能存在着一定的相关性。针对该复杂情形，本文引入多椭球凸模型描述系统参数的不确定性和相关性，即将系统参数分为若干组，每组参数用一个单椭球模型描述，组与组之间相互独立。

假定系统存在n个不确定参数x=（x₁ … x_i … x_n）^T，对于任一 $x_{i}$ ，记其上界为 $x_{i}^{U}$ ，下界为 $x_{i}^{L}$ ，则其中心 $x_{i}^{M}$ 和半径 $x_{i}^{R}$ 定义为：

x_{i}^{M} = \frac{x_{i}^{U} + x_{i}^{L}}{2}, x_{i}^{R} = \frac{x_{i}^{U} - x_{i}^{L}}{2}

（5）

$x_{i}$ 的不确定度表示为 $u_{i} = \frac{x_{i}^{R}}{x_{i}^{M}} \times 100 %$ ，方差表示为 $D (x_{i}) = (x_{i}^{R})^{2}$ 。

将n个不确定参数分成N组，即 $x = {x_{1}^{T} \dots x_{i}^{T}$ $\dots x_{N}^{T}}^{T}$ ， $i = 1,2, \dots, N$ 。

则不确定参数的中心值也会被分成N组： $x^{M} = {\{(x_{1}^{M})^{T} \dots (x_{i}^{M})^{T} \dots (x_{N}^{M})^{T}\}}^{T}, i = 1,2, \dots, N$ 。

对于 $x_{i}^{T}$ ，其单椭球凸模型描述为：

{(x_{i} - x_{i}^{M})}^{T} C_{i}^{- 1} (x_{i} - x_{i}^{M}) \leq 1, i = 1,2, \dots, N

（6）

式中 $C_{i}$ 为椭球凸模型的协方差矩阵。

任意两个不确定变量 $x_{α}$ 和 $x_{β}$ 的协方差定义为： $C o v (x_{α}, x_{β}) = \frac{t a n θ}{1 - t a n^{2} θ} [D (x_{α}) - D (x_{β})]$ ，其中， $θ$ 为椭圆的旋转角度。

为剔除参数级大小不同的影响，引入相关系数描述 $x_{α}$ 和 $x_{β}$ 的相关性：

ρ (x_{α}, x_{β}) = \frac{C o v (x_{α}, x_{β})}{\sqrt[]{D (x_{α})} \cdot \sqrt[]{D (x_{β})}}

（7）

不确定参数 $x$ 的协方差矩阵可表示为：

C_{i} = [\begin{matrix} C o v (x_{1}, x_{1}) & C o v (x_{1}, x_{2}) & \dots & C o v (x_{1}, x_{n}) \\ C o v (x_{2}, x_{1}) & C o v (x_{2}, x_{2}) & \dots & C o v (x_{2}, x_{n}) \\ ⋮ & ⋮ & ⋮ \\ C o v (x_{n}, x_{1}) & C o v (x_{n}, x_{2}) & \dots & C o v (x_{n}, x_{n}) \end{matrix}]

（8）

不确定参数 $x$ 的相关系数矩阵可表示为：

ρ = [\begin{matrix} ρ (x_{1}, x_{1}) & ρ (x_{1}, x_{2}) & \dots & ρ (x_{1}, x_{n}) \\ ρ (x_{2}, x_{1}) & ρ (x_{2}, x_{2}) & \dots & ρ (x_{2}, x_{n}) \\ ⋮ & ⋮ & ⋮ \\ ρ (x_{n}, x_{1}) & ρ (x_{n}, x_{2}) & \dots & ρ (x_{n}, x_{n}) \end{matrix}]

（9）

令 $δ_{i} = x_{i} - x_{i}^{M}$ ，单椭球凸模型也可描述为：

δ_{i}^{T} C_{i}^{- 1} δ_{i} \leq 1, i = 1,2, \dots, N

（10）

此时，N组单椭球模型构成的多椭球凸模型为：

\begin{array}{l} Γ_{x} = \\ \{δ_{1}^{T} C_{1}^{- 1} δ_{1} \leq 1 \dots δ_{i}^{T} C_{i}^{- 1} δ_{i} \leq 1 \dots δ_{N}^{T} C_{N}^{- 1} δ_{Ν} \leq 1\} \end{array}

（11）

式中 $Γ_{x}$ 即为描述制动器系统参数 $x$ 的N组单椭球凸模型构成的不确定模型。

2 响应不确定性分析

开展制动器系统模态阻尼比响应的不确定性分析，能获得响应的上、下边界。将响应边界控制在给定的设计范围内，能有效降低制动噪声的产生倾向^［

7］。以

ζ_{j} (x)

表示系统第j个待研究的模态阻尼比，

x

表示系统参数并且采用多椭球凸模型描述。下面给出三种方法用于求解

ζ_{j} (x)

。

2.1 MCUA方法

首先，基于蒙特卡罗仿真^［

11］提出一种求解

ζ_{j} (x)

的蒙特卡罗不确定性分析（Monte Carlo Uncertainty Analysis，MCUA）方法。MCUA法的主要步骤如下：

（1）根据参数x的相关性，将不确定参数分为N组， $x = {\{x_{1}^{T} \dots x_{i}^{T} \dots x_{N}^{T}\}}^{T} (i = 1,2, \dots, N)$ ，然后构建多椭球凸模型 $Γ_{x}$ 。

（2）假设不确定参数x相互独立，在 $[x^{L}, x^{U}]$ 范围内随机抽取样本。

（3）将上述样本代入多椭球凸模型，如满足模型方程，则将样本保存；若不满足，则舍去。

（4）重复步骤（2）和（3），直至获得 $P$ 组满足要求的样本。

（5）计算每一组满足要求的样本对应的 $ζ_{j} (x)$ 值，得到 $P$ 个 $ζ_{j} (x)$ 响应值。

（6）选取 $P$ 个响应值中的最大和最小值作为 $ζ_{j} (x)$ 的上、下界，分别记为 $ζ_{j}^{U} (x)$ 和 $ζ_{_{j}}^{L} (x)$ ，并得到 $ζ_{j} (x)$ 的响应区间 $[ζ_{_{j}}^{L} (x), ζ_{_{j}}^{U} (x)]$ 。

MCUA的计算精度随着样本数量的增大而提高，当抽取样本数量足够大时，能获得相当精确的结果。因此，MCUA可作为参考方法验证其他分析方法的有效性。

2.2 FPUA方法

MCUA的计算精度严重依赖于样本数量，计算效率往往较低。基于一阶摄动法、拉格朗日乘子法和中心差分法提出一种求解 $ζ_{j} (x)$ 的方法。

$ζ_{j} (x)$ 在 $x^{M}$ 处的一阶泰勒展开为：

\begin{array}{l} ζ_{j} (x) = ζ_{j} (x^{M}) + \sum_{p = 1}^{n} \frac{\partial ζ_{j} (x)}{\partial x_{p}} |_{x = x^{M}} (x_{p} - x^{M}) = \\ ζ_{j} (x^{M}) + {(g_{j}^{M})}^{T} δ \end{array}

（12）

式中 $g_{j}^{M}$ 为 $ζ_{i} (x)$ 的一阶偏导数，可表示为： $g_{j}^{M} = {[\frac{\partial ζ_{j} (x)}{\partial x_{1}} |x = x^{M} \dots \frac{\partial ζ_{j} (x)}{\partial x_{p}} |x = x^{M} \dots \frac{\partial ζ_{j} (x)}{\partial x_{n}} |x = x^{M}]}^{T}$ 。

采用中心差分法求解偏导数，有：

\begin{matrix} \frac{\partial ζ_{j} (x)}{\partial x_{p}} |x = x^{M} = \frac{ζ_{j} (x^{M} + h_{p}) - ζ_{j} (x^{M} - h_{p})}{2 h_{p}}, \\ p = 1,2, \dots, n \end{matrix}

（13）

式中 $h_{p}$ 为差分步长，且 $h_{p} = [0 \dots h_{p} {\dots 0]}^{T}$ 。

$δ$ ， $g_{j}^{M}$ 也可被分成N组， $δ = {\{δ_{1}^{T} \dots δ_{i}^{T} \dots δ_{N}^{T}\}}^{T}$ ， $g_{j}^{M} = {\{{(g_{j, 1}^{M})}^{T} \dots {(g_{j, i}^{M})}^{T} \dots {(g_{j, N}^{M})}^{T}\}}^{T}$ ，i=1，2，…，N。

定义 $ζ_{j} (x)$ 的拉格朗日方程为：

\begin{matrix} \begin{array}{l} L_{j} (δ) = ζ_{j} (x^{M}) + {(g_{j}^{M})}^{T} δ + \sum_{i = 1}^{N} l_{i} (δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1) = \\ ζ_{j} (x^{M}) + \sum_{i = 1}^{N} [{(g_{j, i}^{M})}^{T} δ_{i} + l_{i} (δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1)] \end{array} \end{matrix}

（14）

式中 $l_{1}, \dots, l_{i}, \dots, l_{N}$ 为拉格朗日乘子；结合拉格朗日函数取得极值的必要条件 $\frac{\partial L_{j} (δ)}{\partial δ_{i}}$ =0以及约束条件 $δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1 = 0$ ，可求得 $ζ_{j} (x)$ 的上、下界分别为：

ζ_{_{j}}^{U} (x) = ζ_{j} (x^{M}) + \sum_{i = 1}^{N} \sqrt[]{{(g_{j, i}^{M})}^{T} C_{i} g_{j, i}^{M}}

（15）

ζ_{_{j}}^{L} (x) = ζ_{j} (x^{M}) - \sum_{i = 1}^{N} \sqrt[]{{(g_{j, i}^{M})}^{T} C_{i} g_{j, i}^{M}}

（16）

为表述方便，上述方法称为一阶摄动不确定性分析（First‑order Perturbation Uncertainty Analysis，FPUA）方法。

2.3 SPUA方法

FPUA适用于不确定度较小或者非线性较弱的情况。为进一步提高计算精度，基于二阶摄动法提出一种求解 $ζ_{j} (x)$ 的方法。

基于二阶泰勒展开， $ζ_{j} (x)$ 可表示为：

ζ_{j} (x) = ζ_{j} (x^{M}) + {(g_{j}^{M})}^{T} δ + \frac{1}{2} δ^{T} H_{j}^{M} δ

（17）

式中 $H_{j}^{M}$ 为海森矩阵，可表示为：

H_{j}^{M} = [\begin{matrix} \frac{\partial^{2} ζ_{j} (x)}{\partial x_{1} \partial x_{1}} |_{x = x^{M}} & \frac{\partial^{2} ζ_{j} (x)}{\partial x_{1} \partial x_{2}} |_{x = x^{M}} & \dots & \frac{\partial^{2} ζ_{j} (x)}{\partial x_{1} \partial x_{n}} |_{x = x^{M}} \\ \frac{\partial^{2} ζ_{j} (x)}{\partial x_{2} \partial x_{1}} |_{x = x^{M}} & \frac{\partial^{2} ζ_{j} (x)}{\partial x_{2} \partial x_{2}} |_{x = x^{M}} & \dots & \frac{\partial^{2} ζ_{j} (x)}{\partial x_{2} \partial x_{n}} |_{x = x^{M}} \\ ⋮ & ⋮ & ⋮ \\ \frac{\partial^{2} ζ_{j} (x)}{\partial x_{n} \partial x_{1}} |_{x = x^{M}} & \frac{\partial^{2} ζ_{j} (x)}{\partial x_{n} \partial x_{2}} |_{x = x^{M}} & \dots & \frac{\partial^{2} ζ_{j} (x)}{\partial x_{n} \partial x_{n}} |_{x = x^{M}} \end{matrix}]

（18）

忽略 $H_{j}^{M}$ 的非对角元素，可简化为：

{\bar{H}}_{j}^{M} = d i a g [\frac{\partial^{2} ζ_{j} (x)}{\partial x_{1} \partial x_{1}} |x = x^{M} \dots \frac{\partial^{2} ζ_{j} (x)}{\partial x_{n} \partial x_{n}} |x = x^{M}]

。

采用中心差分法求二阶偏导有：

\begin{array}{l} \frac{\partial^{2} ζ_{j} (x)}{\partial x_{p} \partial x_{p}} |x = x^{M} = \\ \frac{1}{2 h_{p}} [\frac{\partial ζ_{j} (x + h_{p})}{\partial x_{p}} |x = x^{M} \frac{\partial ζ_{j} (x - h_{p})}{\partial x_{p}} |x = x^{M}] = \\ \frac{1}{4 h_{p}^{2}} [ζ_{j} (x^{M} + 2 h_{p}) - 2 ζ_{j} (x^{M}) + ζ_{j} (x^{M} - 2 h_{p})], \\ p = 1,2, \dots, n \end{array}

（19）

根据参数的相关性， $H_{j}^{M}$ 也可以分为N组， ${\bar{H}}_{j}^{M} = d i a g {{\bar{H}}_{j, 1}^{M} {\bar{H}}_{j, 2}^{M} \dots {\bar{H}}_{j, N}^{M}}$ ，其中， ${\bar{H}}_{j, i}^{M} = d i a g [\frac{\partial^{2} ζ_{j} (x_{i})}{\partial x_{i}^{(1)} \partial x_{i}^{(1)}} |x_{i} = x_{i}^{M} \dots \frac{\partial^{2} ζ_{j} (x_{i})}{\partial x_{i}^{(r)} \partial x_{i}^{(r)}} |x_{i} = x_{i}^{M}]$ ，r表示第i组的第r个参数。

$ζ_{j} (x)$ 的表达式可改写为：

\begin{array}{l} ζ_{j} (x) = ζ_{j} (x^{M}) + {(g_{j}^{M})}^{T} δ + \frac{1}{2} δ^{T} {\bar{H}}_{j}^{M} δ = \\ ζ_{j} (x^{M}) + \sum_{i = 1}^{N} [{(g_{j, i}^{M})}^{T} δ_{i} + \frac{1}{2} δ_{i}^{T} {\bar{H}}_{j, i}^{M} δ_{i}] \end{array}

（20）

以多椭球凸模型为约束，定义第i个拉格朗日函数为：

\begin{array}{l} L_{j} (δ) = ζ_{j} (x^{M}) + {(g_{j}^{M})}^{T} δ + \frac{1}{2} δ^{T} {\bar{H}}_{j, i}^{M} δ + \\ \sum_{i = 1}^{N} l_{i} (δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1) = ζ_{j} (x^{M}) + \\ \sum_{i = 1}^{N} [{(g_{j, i}^{M})}^{T} δ_{i} + \frac{1}{2} δ_{i}^{T} {\bar{H}}_{j, i}^{M} δ_{i} + l_{i} (δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1)] \end{array}

（21）

$L_{j} (δ)$ 存在多个方程累加，考虑第i个椭球凸模型 $L_{j, i} (δ) = {(g_{j}^{M})}^{T} δ + \frac{1}{2} δ^{T} {\bar{H}}_{j, 1}^{M} δ + \sum_{i = 1}^{N} l_{i} (δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1)$ ， $(i = 1,2, \dots, N)$ 为非线性方程，此时根据Kuhn‑Tucker条件，需要讨论以下两种情况。

（1）极值点出现在椭球内部时有：

\{\begin{matrix} \frac{\partial L_{j, i} (δ)}{\partial δ_{i}} = g_{j, i}^{M} + {\bar{H}}_{j, i}^{M} δ_{i} = 0 \\ s . t . δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1 < 0 \end{matrix}

（22）

求解式（22）可得 ${\tilde{δ}}_{i}^{(1)} = - {({\bar{H}}_{j, i}^{M})}^{- 1} g_{j, i}^{M}$ 。

（2）极值点出现在椭球边界时有：

\{\begin{array}{l} \frac{\partial L_{j, i} (δ)}{\partial δ_{i}} = g_{j, i}^{M} + {\bar{H}}_{j, i}^{M} δ_{i} + 2 l_{i} C_{i}^{- 1} δ_{i} = 0 \\ s . t . δ_{i}^{T} C_{i}^{- 1} δ_{i} - 1 = 0 \end{array}

（23）

求解式（23）可得第二种情形的解 ${\tilde{δ}}_{i}^{(2)}$ 。综合以上两种情形的解 ${\tilde{δ}}_{i}^{(1)}$ 和 ${\tilde{δ}}_{i}^{(2)}$ ，记为解集 ${\tilde{δ}}_{i}$ ，则存在N组解集 ${\tilde{δ}}_{1}, {\tilde{δ}}_{2}, \dots, {\tilde{δ}}_{N}$ 。对它们进行排列组合，得最终的解集为 $\tilde{δ} = \{{\tilde{δ}}^{(1)}, {\tilde{δ}}^{(2)}, \dots, {\tilde{δ}}^{(k)}, \dots\}$ 。

将这些解代入式（17）得到一系列的 $ζ_{j} (x)$ 值，则 $ζ_{j} (x)$ 的上、下界分别为：

ζ_{j}^{U} (x) = m a x \{ζ_{j} ({\tilde{δ}}^{(1)}), \dots, ζ_{j} ({\tilde{δ}}^{(k)}), \dots\},

ζ_{j}^{L} (x) = m i n \{ζ_{j} ({\tilde{δ}}^{(1)}), \dots, ζ_{j} ({\tilde{δ}}^{(k)}), \dots\}

（24）

为表述方便，上述方法称为二阶摄动不确定性分析（Second‑order Perturbation Uncertainty Analysis，SPUA）方法。

3 响应相关性及椭球域分析

引入相关系数 $ρ (ζ_{α}, ζ_{β})$ 描述系统任意两个待研究的模态阻尼比响应 $ζ_{α} (x)$ 和 $ζ_{β} (x)$ 间的相关性。下面将给出两种方法用于开展系统响应相关性分析，即求解 $ρ (ζ_{α}, ζ_{β})$ 。

3.1 MCCA方法

首先提出一种求解 $ρ (ζ_{α}, ζ_{β})$ 的蒙特卡罗相关性分析（Monte Carlo Correlation Analysis，MCCA）方法。MCCA法的主要步骤如下：

（1）执行MCUA方法的前4个步骤。

（2）计算每一组满足要求的样本对应的 $ζ_{α} (x)$ 和 $ζ_{β} (x)$ 值，得到各响应的P个值。

（3）基于各响应的 $P$ 个值，根据定义可以求得各响应的方差 $D (ζ_{α})$ 和 $D (ζ_{β})$ ，以及两个响应间的协方差 $C o v (ζ_{α}^{}, ζ_{β}^{})$ 。进而，可求得两个响应间的相关系数 $ρ (ζ_{α}, ζ_{β})$ 。

MCCA法计算过程主要基于蒙特卡罗抽样，计算效率较低。

3.2 FPCA方法

为进一步提高相关性分析效率，基于一阶摄动法提出一种求解 $ρ (ζ_{α}, ζ_{β})$ 的方法。

由式（6）可知第i个椭球方程为：

\begin{matrix} Ω_{x_{i}} = \{x_{i} |{(x_{i} - x_{i}^{M})}^{T} {C_{i}}^{- 1} (x_{i} - x_{i}^{M}) \leq 1\}, \\ i = 1,2, \dots, N \end{matrix}

（25）

对椭球方程进行处理，得到：

Ω_{x_{i}} = \{x_{i} |(x_{i} - x_{i}^{M}) {(x_{i} - x_{i}^{M})}^{T} C_{i}^{- 1} (x_{i} - x_{i}^{M}) \cdot {(x_{i} - x_{i}^{M})}^{T} \leq (x_{i} - x_{i}^{M}) {(x_{i} - x_{i}^{M})}^{T}\}

（26）

仅考虑椭球表面时：

\begin{array}{l} {\hat{Ω}}_{x_{i}} = \{x_{i} |(x_{i} - x_{i}^{M}) {(x_{i} - x_{i}^{M})}^{T} \cdot \\ ({C_{i}}^{- 1} (x_{i} - x_{i}^{M}) {(x_{i} - x_{i}^{M})}^{T} - E) = 0\} \end{array}

（27）

式中 $E$ 为单位矩阵。因此，有 ${C_{i}}^{- 1} (x_{i} - x_{i}^{M})$ ${(x_{i} - x_{i}^{M})}^{T} - E = 0$ 。即协方差矩阵可表示为：

C_{i} = (x_{i} - x_{i}^{M}) {(x_{i} - x_{i}^{M})}^{T}

（28）

得到任意两个响应 $ζ_{α}, ζ_{β}$ 的相关系数表达式为：

\begin{array}{l} C o v (ζ_{α}, ζ_{β}) = (g_{α}^{M})^{T} (x - x^{M}) {((g_{β}^{M})^{T} (x - x^{M}))}^{T} = \\ \sum_{i = 1}^{N} (g_{α, i}^{M})^{T} C_{i}^{T} (g_{β, i}^{M})^{T} + \\ \sum_{i = 1}^{N} \sum_{j = 1, j \neq i}^{N} (g_{α, i}^{M})^{T} (x_{i} - x_{i}^{M}) {(x_{j} - x_{j}^{M})}^{T} g_{β, i}^{M} \end{array}

（29）

等式两边同时在不确定域内积分，且N个椭球域都是关于中心对称的，可得：

\begin{array}{l} \int \overset{n}{\underset{Ω x_{1} + \dots + Ω x_{n}}{. . .}} \int C o v (ζ_{α}, ζ_{β}) d x_{1} d x_{2} \dots d x_{n} = \\ \int \overset{n}{\underset{Ω x_{1} + \dots + Ω x_{n}}{. . .}} (g_{α, i}^{M})^{T} (x_{i} - x_{i}^{M}) {(x_{j} - x_{j}^{M})}^{T} g_{β, i}^{M} \end{array}

因此，两个响应间的协方差为：

C o v (ζ_{α}, ζ_{β}) = \sum_{i = 1}^{N} (g_{α, i}^{M})^{T} C_{i}^{T} g_{β, i}^{M}

（30）

最后可求得相关系数为：

ρ (ζ_{α}, ζ_{β}) = \frac{C o v (ζ_{α}, ζ_{β})}{\sqrt[]{C o v (ζ_{α}, ζ_{α})} \sqrt[]{C o v (ζ_{β}, ζ_{β})}} \approx \frac{\sum_{i = 1}^{N} (g_{α, i}^{M})^{T} C_{i}^{T} g_{β, i}^{M}}{\sqrt[]{\sum_{i = 1}^{N} (g_{α, i}^{M})^{T} C_{i}^{T} g_{α, i}^{M}} \sqrt[]{\sum_{i = 1}^{N} (g_{β, i}^{M})^{T} C_{i}^{T} g_{β, i}^{M}}}

（31）

为表述方便，上述方法称为一阶摄动相关性分析（First‑order Perturbation Correlation Analysis，FPCA）方法。

3.3 椭球域分析

得到任意两个响应的上、下边界以及相关系数后，可建立一个椭球域用于量化 $ζ_{α} (x)$ 和 $ζ_{β} (x)$ 的不确定域，其数学表示为：

{[\begin{array}{l} ζ_{α} (x) - ζ_{α} (x^{M}) \\ ζ_{β} (x) - ζ_{β} (x^{M}) \end{array}]}^{T} C_{}^{- 1} (ζ_{α}, ζ_{β}) [\begin{array}{l} ζ_{α} (x) - ζ_{α} (x^{M}) \\ ζ_{β} (x) - ζ_{β} (x^{M}) \end{array}] \leq 1

（32）

式中 $C (ζ_{α}, ζ_{β})$ 为协方差矩阵：

C (ζ_{α}, ζ_{β}) = [\begin{matrix} D (ζ_{α} (x)) & C o v (ζ_{α} (x), ζ_{β} (x)) \\ C o v (ζ_{α} (x), ζ_{β} (x)) & D (ζ_{β} (x)) \end{matrix}]

，其中，

D (ζ_{α} (x)) = {(\frac{ζ_{α}^{U} {(x)}^{} - ζ_{α}^{L} (x)}{2})}^{2}

，

D (ζ_{β} (x)) = {(\frac{ζ_{β}^{U} {(x)}^{} - ζ_{β}^{L} (x)}{2})}^{2}

，

C o v (ζ_{α} (x), ζ_{β} (x)) = ρ (ζ_{α}, ζ_{β}) \sqrt[]{D (ζ_{α} (x))} \sqrt[]{D (ζ_{β} (x))}

。

根据MCUA和MCCA，可求得 $ζ_{α} (x)$ 和 $ζ_{β} (x)$ 边界和相关系数 $ρ (ζ_{α}, ζ_{β})$ ，进而建立响应的椭球域。将结合MCUA和MCCA求响应椭球域的方法简称为MCUA‑MCCA法。类似地，可得到FPUA‑FPCA法和SPUA‑FPCA法。

4 算例分析

4.1 研究模型

以图1所示的制动器系统为研究对象，由文献［

8］可知该系统在制动过程中存在的制动噪声问题。结合文献［8，12］，选取系统研究参数为：制动盘与制动块之间的摩擦系数

f

、制动压力

p

、摩擦材料密度

ρ_{l}

和摩擦材料的杨氏模量

e_{l}

。这些参数与摩擦接触特性紧密相关，在工程实际中表现出较强的不确定性。故将这些参数均视为不确定参数，其取值如表1所示。同时，根据文献［13‑14］可知摩擦系数

f

和制动压力

p

存在相关性，摩擦材料密度

ρ_{l}

和杨氏模量

e_{l}

亦存在一定相关性。

表1 参数不确定取值

Tab.1 The uncertain values of parameters

取值	$f$	$p / M P a$	$ρ_{l} / (k g \cdot m^{- 3})$	$e_{l} / G P a$
标称值	0.5	1.2	30	6.4
最小值	0.44	0.9	27	5.8
最大值	0.56	1.5	33	7.6

结合文献［

1］可知，在1~16 kHz的范围内，该模型在7.28和12.9 kHz左右的不稳定模态的阻尼比绝对值较大。因此，选取这两个模态阻尼比作为主要响应来反映系统的不稳定性。基于有限元分析和响应面法，可建立两个不稳定模态的阻尼比响应函数，分别记为

ζ_{1} (x)

和

ζ_{2} (x)

：

\begin{array}{l} ζ_{1} (x) = ζ_{1} (f, p, ρ_{l}, e_{l}) = \\ 10.09273 - 24.84553 f + 0.71694 p - \\ 0.48950 ρ_{l} - 1.12423 e_{l} - 0.21561 f p + \\ 1.29826 f ρ_{l} + 0.36183 f e_{l} - 0.23732 p ρ_{l} + \\ 0.01364 p e_{l} - 0.05786 ρ_{l} e_{l} + 16.98764 f^{2} - \\ 0.00558 p^{2} + 0.04551 ρ_{l}^{2} + 0.09857 e_{l}^{2} \end{array}

（33）

\begin{array}{l} ζ_{2} (x) = ζ_{2} (f, p, ρ_{l}, e_{l}) = \\ 1.8835 + 7.1091 f - 0.4801 p - \\ 0.02224 ρ_{l} - 1.1163 e_{l} + 1.0553 f p - \\ 0.6360 f ρ_{l} + 1.0196 f e_{l} + 1.0196 p ρ_{l} - \\ 0.0816 p e_{l} + 0.0619 ρ_{l} e_{l} - 12.86144 f^{2} + \\ 0.05915 p^{2} - 0.0304 ρ_{l}^{2} + 0.0441 e_{l}^{2} \end{array}

（34）

式中 $x = {[f p ρ_{l} e_{l}]}^{T}$ ，并采用多椭球凸模型描述其不确定性和相关性。

两个响应面模型的决定系数 $R^{2}$ 分别为0.9837和0.8932，模型 $F$ 值分别为64.83和17.73，两模型与有限元模型逼近程度较高，能较好地满足预测精度要求。

4.2 不确定性分析

考虑系统不确定参数间的相关性，设定相关系数 $ρ (f, p) = ρ (ρ_{l}, e_{l}) = 0.5$ 。为分析系统参数不确定性对响应不确定性的影响，以及分析FPUA和SPUA适用的不确定度范围，给定一系列参数的不确定度，使用不同方法计算响应的边界，表2和3给出了不同情形下MCUA，FPUA和SPUA的计算结果。其中，MCUA的样本数为10⁵。

表2 三种方法求得的

ζ_{1} (x)

范围

Tab.2 The ranges of

ζ_{1} (x)

obtained by the three methods

不确定度u	MCUA		FPUA		SPUA
不确定度u	下界	上界	下界	上界	下界	上界
±1%	-0.35008	-0.31279	-0.35079	-0.31325	-0.35013	-0.31257
±3%	-0.38237	-0.27119	-0.38832	-0.27571	-0.38256	-0.26952
±5%	-0.41006	-0.22418	-0.42586	-0.23817	-0.41025	-0.22062
±7%	-0.43351	-0.17289	-0.46340	-0.20063	-0.43352	-0.16563
±9%	-0.45201	-0.11715	-0.50094	-0.16309	-0.45264	-0.10437

以表2和3为基础，绘制出MCUA，FPUA和SPUA的计算结果如图2所示。计算FPUA和SPUA相较于MCUA的相对误差，如图3所示。

图2 三种方法求解 $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 的结果

Fig.2 The results of $ζ_{1} (x)$ and $ζ_{2} (x)$ solved by the three methods

图3 两种方法计算的相对误差

Fig.3 The relative errors calculated by two methods

由图2可知， $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 也存在不确定性，且随参数不确定度增大， $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 的不确定区间逐渐增大。此外，当 $u \leq \pm 3 %$ 时，三种方法的曲线几乎是重合的，说明此时FPUA和SPUA方法的计算精度均很高；而随着参数不确定度增大，FPUA和SPUA方法求得的曲线逐渐偏离MCUA求得的参考值曲线，其中FPUA的曲线差异更明显。在求解 $ζ_{1} (x)$ 的响应时，SPUA方法求得的曲线与参考曲线基本重合。这说明SPUA的计算精度整体上优于FPUA。

表3 三种方法求得的

ζ_{2} (x)

范围

Tab.3 The ranges of

ζ_{2} (x)

obtained by the three methods

不确定度u	MCUA		FPUA		SPUA
不确定度u	下界	上界	下界	上界	下界	上界
±1%	-0.05665	-0.04956	-0.05688	-0.04974	-0.05702	-0.04988
±3%	-0.06305	-0.04071	-0.06401	-0.04261	-0.06534	-0.04249
±5%	-0.07301	-0.02943	-0.07114	-0.03547	-0.07488	-0.03367
±7%	-0.08648	-0.01583	-0.07828	-0.02834	-0.08588	-0.02341
±9%	-0.10297	-0.00040	-0.08541	-0.02121	-0.09898	-0.01170

由图3可知，对于 $ζ_{1} (x)$ 下界，当 $u \geq \pm 5 %$ 时，FPUA方法计算的相对误差大于5%，但SPUA方法计算的相对误差均在5%以内；对于 $ζ_{1} (x)$ 上界，当 $u \geq \pm 5 %$ 时，FPUA方法计算的相对误差超过5%，而SPUA方法计算的相对误差均在5%以内；对于 $ζ_{2} (x)$ 下界，当 $u \geq \pm 7 %$ 时，FPUA的相对误差超过5%，但SPUA的相对误差均在5%以内；对于 $ζ_{2} (x)$ 上界，由于其值接近0，两种方法计算的相对误差均较大，FPUA方法和SPUA方法均在 $u \geq \pm 3 %$ 后相对误差超过5%。总的来说，SPUA方法的计算精度高于FPUA方法，且适用的不确定度范围更大。

在计算效率方面，在求解上述响应结果时，MCUA用时48 s，FPUA用时0.25 s，SPUA用时0.4 s。可见FPUA和SPUA在求解系统响应时均具有较高的计算效率。与FPUA相比，SPUA方法的计算效率略低，但大大提高了计算精度。

4.3 相关性分析

令参数不确定度 $u = \pm 3 %$ ，分5种情况考虑两组参数的相关性，即 $ρ_{0} = ρ (f, p) = ρ (ρ_{l}, e_{l}) =$ 0，0.3，0.5，0.7和0.9。然后分别采用MCCA方法和FPCA方法求解 $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 之间的相关系数，其中，MCCA的样本数为10⁵，结果如表4所示。以MCCA方法为参考，给出FPCA方法计算的相对误差，结果如图4所示。

表4 相关系数求解结果

Tab.4 The results of correlation coefficients

相关系数 $ρ_{0}$	$ρ (ζ_{1}, ζ_{2})$
相关系数 $ρ_{0}$	MCCA	FPCA
0	0.4752	0.4958
0.3	0.3973	0.4102
0.5	0.3372	0.3470
0.7	0.2708	0.2769
0.9	0.1941	0.1967

图4 FPCA方法计算的相对误差

Fig.4 The relative errors calculated by FPCA

由表4和图4可知，以MCCA为参考，FPCA求解响应相关系数的相对误差均在4.5%以下，两种方法的计算结果较为接近。此外，还可看出，参数相关性越大时，响应相关性的计算误差越小。

绘制响应的相关系数与参数的相关系数曲线图，如图5所示，由图可见，两种方法曲线的重合度较高，也体现了FPCA具有较高的计算精度。此外，随着参数的相关系数逐渐增大， $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 的相关系数逐渐减小，即两响应的关联程度逐渐减弱。

图5 响应的相关系数与参数的相关系数曲线图

Fig.5 Correlation coefficient of response and correlation coefficient curve of parameter

特别地，当系统参数的相关系数为0，即不考虑参数的相关性时，系统响应之间的相关系数并不为0。这表明响应间的相关性不是完全来自系统参数，当参数独立时，响应间也可能存在相关性。

在计算效率方面，MCCA用时5 s，FPCA用时0.3 s。可见FPCA在求解系统响应相关性时具有很高的计算效率。

综上，FPCA方法在求解系统响应相关性时具有较高的计算精度和计算效率。

4.4 椭球域分析

获得响应的不确定边界和相关性后，可建立响应的椭球域来直观反映响应的不确定性和相关性。分别采用MCUA‑MCCA，FPUA‑FPCA和SPUA‑ FPCA三种组合方法研究以下两种情形的响应椭球域：（1）保持参数不确定度 $u = \pm 5 %$ 不变，改变参数的相关系数取值为 $ρ_{0} =$ 0.3，0.5，0.7和0.9，分析结果如图6所示；（2）保持参数的相关系数 $ρ_{0} = 0.5$ 不变，改变参数的不确定度为 $u = \pm 3 %$ ， $\pm 5 %$ ， $\pm 7 %$ 和 $\pm 9 %$ ，分析结果如图7所示。

图6 不同参数相关系数下的椭圆域

Fig.6 The elliptical domains under different parameter correlation coefficients

图7 不同参数不确定度下的椭圆域

Fig.7 The elliptical domains under different parameter uncertainties

由图6可知，以MCUA‑MCCA求得的结果为参考，SPUA‑FPCA求得的椭圆域比FPUA‑FPCA求得的椭圆域更接近参考值，表明SPUA‑FPCA的计算结果更为精确。同时，随着参数相关系数增大，椭圆旋转角度变小且保持正值，椭圆域逐渐变圆。即响应 $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 的相关系数在减小，正相关程度在降低，这与前文分析一致。

由图7可知，随着参数不确定度增大，三种方法求得的椭圆域均增大，椭圆旋转角度基本不变。说明响应 $ζ_{1} (x)$ 和 $ζ_{2} (x)$ 的不确定范围在扩大，但正相关性变化不大。此外，FPUA‑FPCA和SPUA‑FPCA求得的椭圆域随着参数不确定度增大逐渐偏离参考值，但SPUA‑FPCA求得的结果始终接近参考值。这也说明SPUA‑FPCA的求解精度更高。

5 结论

随着参数不确定度增大，制动器系统稳定性响应的不确定域变大；以MCUA方法作为参考，SPUA方法计算的响应不确定边界具有更高的精度，且适用于更大的参数不确定度范围，FPUA方法仅适用于参数不确定度小的情形；FPUA和SPUA方法均具有较高的计算效率。

在求解制动器系统稳定性响应的相关性方面，以MCCA方法为参考，所提出的FPCA方法求得的响应相关系数是有效的，且计算效率更高。

随着不确定参数相关性增大，所研究的制动器系统两个稳定性指标响应之间的相关系数减小，正相关性减弱。本文方法亦适用于盘式制动器外的其他类型制动器系统。

参考文献

LÜ H， YU D J. Brake squeal reduction of vehicle disc brake system with interval parameters by uncertain optimization［J］. Journal of Sound and Vibration， 2014， 333（26）： 7313-7325. [百度学术]

张立军，庞明，孟德建，等.面向制动尖叫抑制的制动块稳健性设计［J］.汽车工程，2016，38（1）：65-71. [百度学术]

ZHANG Lijun， PANG Ming， MENG Dejian， et al. Robust design of brake pad for brake squeal suppression［J］. Automobile Engineering， 2016， 38（1）： 65-71. [百度学术]

SARROUY E， DESSOMBZ O， SINOU J J. Piecewise polynomial chaos expansion with an application to brake squeal of a linear brake system［J］. Journal of Sound and Vibration， 2013，332（3）： 577-594. [百度学术]

SARROUY E， DESSOMBZ O， SINOU J J. Stochastic study of a non-linear self-excited system with friction［J］. European Journal of Mechanics-A/Solids， 2013， 40： 1-10. [百度学术]

NOBARI A， OUYANG H， BANNISTER P. Uncertainty quantification of squeal instability via surrogate modelling［J］. Mechanical Systems and Signal Processing， 2015， 60： 887-908. [百度学术]

黄晓婷，李沛航，吕辉. 含模糊不确定性的汽车盘式制动器稳定性研究［J］. 重庆理工大学学报（自然科学）， 2021， 35（8）： 48-55. [百度学术]

HUANG Xiaoting， LI Peihang， LÜ Hui. Research on the stability of automotive disc brakes with fuzzy uncertainty［J］. Journal of Chongqing University of Technology（Natural Science）， 2021， 35（8）： 48-55. [百度学术]

吕辉，上官文斌，于德介. 基于证据理论的汽车制动器系统稳定性分析［J］. 华南理工大学学报（自然科学版）， 2019， 47（3）： 53-60. [百度学术]

LÜ Hui， SHANGGUAN Wenbin， YU Dejie. Stability analysis of automotive brake systems based on evidence theory［J］. Journal of South China University of Technology（Natural Science Edition）， 2019， 47（3）： 53-60. [百度学术]

LÜ H， SHANGGUAN W B， YU D J. An imprecise probability approach for squeal instability analysis based on evidence theory［J］. Journal of Sound and Vibration， 2017， 387（20）： 96-113. [百度学术]

贾爱芹，陈建军，徐亚兰.基于摄动法的不确定性汽车悬架振动控制特征值的凸模型分析［J］.中南大学学报（自然科学版），2012，43（4）： 121-125. [百度学术]

JIA Aiqin， CHEN Jianjun， XU Yalan. Convex model analysis of vibration control eigenvalues of vehicle suspension system based on perturbation method［J］. Journal of Central South University（Science and Technology）， 2012，43（4）： 121-125. [百度学术]

王攀，臧朝平. 改进的平行六面体凸模型识别动力学不确定参数区间的方法［J］. 振动工程学报，2019，32（1）： 97-106. [百度学术]

WANG Pan， ZANG Chaoping. Method of identifying dynamic uncertain parameter intervals with improved parallelepiped convex model［J］. Journal of Vibration Engineering， 2019， 32（1）： 97-106. [百度学术]

CAI B H， SHANGGUAN W B， LÜ H， et al. Hybrid uncertainties-based analysis and optimization design of powertrain mounting systems［J］. Science China Technological Sciences， 2020， 63（5）： 838-850. [百度学术]

RENAULT A， MASSA F， LALLEMAND B， et al. Experimental investigations for uncertainty quantification in brake squeal analysis［J］. Journal of Sound and Vibration， 2016， 367： 37-55. [百度学术]

BIJWE J， KUMAR M. Optimization of steel wool contents in non-asbestos organic （NAO） friction composites for best combination of thermal conductivity and tribo-performance［J］. Wear， 2007， 263（7-12）： 1243-1248. [百度学术]

ASHBY M F， CEBON D. Materials selection in mechanical design［J］. Le Journal de Physique Ⅳ， 1993， 3（C7）： C7-1-C7-9. [百度学术]