基于宏微观摩擦接触模型的湿式摩擦元件温度场失效概率分析

Failure Probability Analysis of Temperature Field in Wet Friction Components Based on Macro-Micro Frictional Contact Model

摘要

为精确量化湿式摩擦元件的温度场失效概率问题，先建立热-机械耦合仿真与单对粗糙峰接触模型实时交互的宏微观摩擦模型获取温度，再构建温度场失效概率统计模型。利用核密度估计建立失效参数的概率密度函数，采用蒙特卡罗模拟进行了概率计算。试验结果表明，仿真与试验的数据高度吻合，所建立的统计模型能有效且精确计算湿式摩擦元件的失效概率。

Abstract

Aiming at accurately quantifying the failure probability of temperature fields for wet friction components， a macro-micro friction contact model， with the interaction of thermal-mechanical coupling simulation and single asperity contact model， was established to acquire temperature， and a statistical model for temperature field failure probability was constructed. The method utilized kernel density estimation to establish the probability density function of failure parameters and emploied Monte Carlo simulation for probability calculation. Experimental results show high agreement between simulation data and test data. The established statistical model may effectively and accurately calculate the failure probability of wet friction components.

Graphical abstract

关键词

湿式摩擦副 / 失效概率分析 / 宏微观交互 / 核密度函数 / 蒙特卡罗模拟

Key words

wet friction pair / failure probability analysis / macro-micro interaction / kernel density function / Monte Carlo simulation

0 引言

湿式摩擦元件是重型车辆、船舶等装备动力传输的关键部件，其失效会严重影响设备安全。工作时产生的高温易导致元件翘曲变形，引发传动故障^［1］，因此，研究其温度场失效对保障设备可靠具有重要意义。

早期的摩擦元件温度研究基于一维模型^［2］开展，即假设热量仅在沿接触体深度方向传播，通过推导一维瞬态温度方程计算接触体内部不同深度的温度变化。随着计算能力的增强和数值分析方法的发展，二维模型^［3］、三维模型^［4］因能更准确模拟和分析摩擦元件表面温度场而逐渐取代一维模型。李和言等^［5］建立摩擦热流密度与真实接触面积的动态滑摩过程分析模型，通过动态测温和静态比压试验验证了理论模型对接触比压分布研究的适用性。王立勇等^［6］采用有限元模拟和实验的方法对比分析了5种不同的径向非均布压力分布方式对湿式摩擦副工作过程中热机耦合的作用。郑良杰等^［7］建立湿式多片离合器分离过程动力学数值模型，研究不同花键摩擦因数下湿式多片离合器的分离过程。SHIN等^［8］提出一种考虑热动力学的湿式离合器摩擦建模方法。ABDULLAH等^［9］采用热分配法和总热产生法模拟汽车离合器的热模型以获得温度场。为更准确地研究温度场的变化，一些学者在此基础上建立宏观与微观交互的仿真模型进行温度场研究。FURMAŃSKI等^［10］提出一个宏微观导热模型来预测具有复杂微观结构的非均相介质的局部温度分布，获得微观温度场和宏观温度场之间的关系。MBEN等^［11］将微尺度流体动力学模型与热力学模型相结合，并使用所得算法推导淬火速度和温度之间的关系。研究温度场的变化规律有助于确定摩擦元件的潜在失效，因此需要建立宏微观摩擦接触模型，研究摩擦元件的温度场变化规律，以揭示温度对摩擦元件性能的影响，从而优化设计与材料选择，提高摩擦元件的可靠性和使用寿命。

尽管仿真模型是摩擦元件温度场研究中重要的定量分析工具，但其结果通常较为单一，存在一定的局限性，无法全面反映实际工况下的复杂性和不确定性。与之相比，温度场失效概率统计模型能结合温度梯度变化更清晰地表征特定失效概率下的摩擦元件是否接近失效，更符合实际失效规律，具备更高的客观性、鲁棒性。近年来，众多学者针对摩擦元件的失效展开了一系列的研究。吴健鹏等^［12］基于摩擦元件的实际工况设计加速寿命试验，根据摩擦因数变化规律和最大径向温差划分了滑摩特性的不同阶段，探究湿式摩擦副温度场变化规律和失效机理。李乐等^［13］构建了一种基于PSO-BP神经网络的湿式摩擦元件损伤预测模型。蔡运迪等^［14］建立盘式制动器三维热-机耦合模型来获取制动盘表面及内部的温度场与应力场的分布特征，并以此为基础分析了制动盘热疲劳失效的机理。GANESH等^［15］对非公路变速器中使用的两种离合器系统进行了现场故障原因分析。WANG等^［16］基于离合器故障的混合传播特性构建了系统故障诊断策略，在系统层面揭示离合器的失效机理。当前研究表明，摩擦元件温度场失效的区域性分布规律尚未得到系统揭示。局部高温会显著改变摩擦性能、加速磨损，并可能引发热翘曲变形，最终导致温度场失效。

为此，本文以湿式摩擦元件为研究对象，首先建立湿式摩擦元件宏微观摩擦接触模型，获取温度场数据，并通过全寿命周期试验对仿真模型进行验证。利用核密度函数建立温度场失效特征参数的概率密度函数，构建温度场失效极限状态函数，通过蒙特卡罗模拟得出摩擦元件温度场失效概率，分析摩擦元件的温度场失效规律，并通过试验数据对失效概率统计模型进行验证。最后分析了转速和压力对摩擦元件温度场失效概率的影响。

1 宏微观摩擦接触模型

1.1 摩擦副热固耦合仿真模型

湿式离合器中的摩擦副由钢片和摩擦片组成，其中，钢片的材料为65Mn，摩擦片的材料为铜基粉末冶金。摩擦片和钢片都有2个接触面并且两者都呈轴对称分布，鉴于摩擦元件的对称性，为简化计算，选取其轴向厚度的一半建立分析模型，并在两端面施加相应的边界条件。摩擦元件的几何尺寸和材料属性如表1所示。

钢片在不同温度有不同的热传导规律，其热导率随温度的升高呈一定下降趋势，见表2。

使用ABAQUS建立摩擦副三维有限元仿真模型（径向划分30个网格，周向划分90个网格）。为研究接合过程中接触表面的温度场分布特性，该模型考虑了摩擦副实际运转中的接触压力、相对旋转运动和摩擦接触等因素。使用该模型模拟不同转速和压力下的温升和温度场。压力1.5 MPa、转速1000 r/min的模拟钢片温度场见图1。

不同工况下，钢片的高温区集中在中径附近，低温区集中在内径和外径。压力1.5 MPa、转速1000 r/min工况下摩擦元件的最高温度为中外径区域的88.04 ℃，远远高于内径和外径的最高温度。这是因为内径和外径的散热条件优于中径处，且内径处的热流密度更低。

1.2 宏微观摩擦接触模型

塑性指数计算公式为

\Omega =\frac{E_{\mathrm{e}}}{H_{\mathrm{f}}}\sqrt[]{\frac{{\sigma }_{\mathrm{e}}}{\beta }}（1）

E_e=1/［（1-ν{}_{\mathrm{1}}^{\mathrm{2}}）/E₁+（1-ν{}_{\mathrm{2}}^{\mathrm{2}}）/E₂］

式中：E_e为等效弹性模量；ν₁、ν₂分别为钢片和摩擦片的泊松比；E₁、E₁分别为钢片和摩擦片的弹性模量；σ_e为粗糙度的均方根；H_f为摩擦片材料的硬度，H_f=2.8σ_y；σ_y为摩擦片的屈服强度；β为两粗糙表面接触模型的表面膜破坏分数，通常取0.5^［17］。

根据塑性指数判断接合过程中钢片上的接触状态，塑性指数Ω<0.6时，摩擦面处于弹性接触状态；0.6<Ω<1.4时，发生从纯弹性到弹塑性变形的转变；1.4<Ω<8时，弹性变形和塑性变形共同存在于粗糙表面的摩擦接触过程^［18］。表面材料为铜基粉末冶金时，摩擦盘在接合过程中通常会发生弹塑性变形。实际粗糙接触面积为

A_{\mathrm{c}}=k{\mathrm{\pi }}^{\mathrm{3}}N({\beta }_{\mathrm{1}}{\sigma }_{\mathrm{e}}{)}^{\mathrm{2}}(R_{\mathrm{o}}^{\mathrm{2}}-R_{\mathrm{i}}^{\mathrm{2}})A_{\mathrm{r}\mathrm{e}\mathrm{d}}{F}_{\mathrm{2}}(H)（2）

式中：k为塑性变形系数，k∈［1，2］；N为微凸峰的密度，N=3×10⁷ kg/m³；β₁为微凸峰的曲率半径；H为膜厚比，H=h_y；h_y为油膜厚度；\sigma为表面粗糙度的均方根；；A_red为非沟槽区域面积比；R_o为摩擦元件的外半径；R_i为摩擦元件的内半径。

粗糙峰最初与刚性平面接触时，粗糙峰顶点处的接触应力相对较低，应力分布区域较小。这意味着粗糙峰处于弹性变形阶段，塑性指数Ω<0.6，满足赫兹弹性接触理论。随后，确定粗糙峰顶点处的变形，变形与平均接触压力之间的关系表示为^［19］

\omega =\frac{\mathrm{9}{\mathrm{\pi }}^{\mathrm{2}}{p}_{\mathrm{m}\mathrm{e}}^{\mathrm{2}}{r}_{\mathrm{e}}}{\mathrm{16}{E}_{\mathrm{e}}^{\mathrm{2}}}（3）

式中：p_me为微凸体处于弹性变形阶段的平均接触压力；r_e为单微凸体顶点的曲率半径。

通过求解曲率半径的公式，可得粗糙峰顶点的曲率半径

r_{\mathrm{e}}={(\frac{{\mathrm{d}}^{\mathrm{2}}y(x)}{\mathrm{d}{x}^{\mathrm{2}}})}^{-\mathrm{1}}=\frac{b^{\mathrm{2}}}{\mathrm{2}{h}_{\mathrm{w}}}（4）

式中：b为微凸体椭圆截面的长半轴；h_{\mathrm{w}}为微凸体高度。

接触热导是通过界面的热通量密度与两个接触界面之间温差的比值。通过粗糙峰法向载荷可计算得到摩擦元件的接触热导。接触热导与粗糙峰法向载荷的关系如下：

K=\frac{\mathrm{2}{\lambda }_{\mathrm{s}}}{\Phi }{(\frac{\mathrm{3}{r}_{\mathrm{e}}p}{\mathrm{4}{E}_{\mathrm{e}}})}^{\mathrm{1}/\mathrm{3}}（5）

式中：K为单个微凸体接触热导；p为微凸体法向接触载荷；λ_s 为摩擦片热导率；Φ为接触热导因子。

摩擦副之间的对流传热系数的表达式为

h_{\mathrm{1}}=\mathrm{0.332}{\lambda }_o{(Re)}^{\mathrm{1}/\mathrm{2}}{(Pr)}^{\mathrm{1}/\mathrm{3}}/r（6）

Pr=z{c}_{\mathrm{1}}{\rho }_{\mathrm{1}}/{\lambda }_o（7）

式中，λ_o为油液热导率，λ_o=0.3 W/（m·K）；Re为雷诺数；Pr为普朗特系数；r为摩擦副径向半径；z为油液运动黏度；z=31×10{}^{-\mathrm{6}} m²/s；c₁为油液比热，c=2231 J/（kg·K）；ρ₁为油液密度，ρ₁=850 kg/m³。

为提高温度场仿真结果的准确性，本文建立了与宏观热-机械耦合仿真模型实时交互的单对粗糙峰接触仿真模型。粗糙峰的初始顶点曲率半径、热导率分别为8×10{}^{-\mathrm{4}} m和30 W/（m·K）。微观模型的曲率半径、高度与宏观模型的对流传热系数、摩擦因数等参数的设置充分考虑了摩擦元件真实表面的随机粗糙度。将微观模型与宏观模型进行实时交互以达到更准确的摩擦元件温度场仿真效果，宏微观摩擦接触模型交互示意如图2所示。

湿式摩擦元件的宏微观摩擦接触模型在ABAQUS平台上仿真。首先，根据式（3）计算粗糙峰弹性变形阶段的变形量ω。将此变形量分配给微观接触仿真模型，进行微观模型的仿真，单次仿真的时间步长为0.1 s。仿真完成后，提取钢片上粗糙峰顶点处的压力和温度。根据式（5）求解钢片粗糙峰的热导率，并将其分配给宏观接触仿真模型。更新宏观模型的参数后，开始下一轮仿真，确保每次仿真迭代选择的节点位置一致。压力0.3 MPa、转速500 r/min下，钢片的界面温度随滑摩时间的延长而增大，如图3所示。

分析摩擦界面的温度分布规律发现，钢片的高温区主要分布在中径和中外径处，低温区域主要分布在内径和外径。随着交互次数的增加，摩擦界面的整体温度呈现上升趋势。最高温度从第一次交互时的25.23 ℃上升至第十次交互时的104.4 ℃。

2 湿式离合器摩擦性能全寿命试验

为获取摩擦元件温度，研究表面温度对摩擦元件综合性能的影响，并验证仿真模型的可靠性和合理性，搭建了湿式离合器摩擦磨损试验台。在钢片埋入热电偶以测量温度，并制定试验方案。

2.1 湿式离合器试验台

如图4所示，试验台系统主要包括动力装置、惯量、离合器箱体、制动装置、液压系统、数据采集装置等。试验装置的润滑油压由泵站提供，控制油通过油路进入离合器包箱内油缸，推动活塞，使离合器摩擦片和钢片压紧接合；切断控制油路的压力油供给，回位弹簧将活塞推回，实现摩擦片和钢片的分离。

2.2 试验方案

在钢片埋入热电偶以测量18个点位的温度，如图5所示。摩擦元件的温升是随时间逐渐累积的过程，为保障试验的安全性，摩擦副不宜长时间处在高速重载环境下，因此交替进行多工况短时滑摩与单工况长时滑摩。多工况短时滑摩的目的是加快摩擦元件的损伤，使摩擦副在短时间的极限工况下近似受损。单工况长时滑摩试验旨在模拟摩擦副的实际失效过程，即通过在设定工况下使其发生损伤，以加速其寿命演化。基于此，本研究将从此项加速寿命试验中获取的单一工况长时滑摩损伤参数，作为该工况下的损伤特征参数进行分析。

多工况短时滑摩试验是不同压力与转速的正交试验，1.0 MPa、2.0 MPa与1000 r/min、2000 r/min的组合中，每次接合时间为2s，接合总时间为8 s。基于试验安全性，不宜设定较大的压力和转速。对此，单一工况长时滑摩阶段的压力为0.3、0.5、0.6、0.8、1.0 MPa，转速为200、400、500、600、800、1000、1500 r/min，每次接合时间为10 min。由此得出不同工况的摩擦元件温度，通过测量得到这些实验数据对宏微观摩擦接触模型进行验证。试验过程中出现异响和振动，且摩擦因数发生剧烈变化时停止实验。停止实验后发现摩擦元件出现翘曲变形，此时，可判定摩擦元件因温度场异常而发生失效。

2.3 宏微观摩擦接触模型验证

通过单一工况长时滑摩中各个工况的温度验证宏微观摩擦接触模型。1.0 MPa、1500 r/min工况的初始温度为20 ℃，宏观仿真模型仿真时间为1 s，宏微观摩擦接触模型的时间步长为0.1 s，循环10次后，提取钢片表面中径处的温度，对比试验相同位置的温度和宏观模型相应位置的仿真数据，仿真与试验结果对比如图6所示。为确保试验测量的温度准确可靠，对相同工况进行多次重复试验，剔除异常数据后取多次数据的均值。

由图6可看出，随着时间的增加，宏微观摩擦接触模型的仿真结果越来越逼近试验结果。10次循环交互中，宏微观摩擦接触模型的仿真结果始终优于单一宏观模型；宏微观摩擦接触模型在4次交互后的仿真结果大于真实值；随着交互次数的增加，宏微观摩擦接触模型的仿真结果不断逼近试验数据。10次交互完成后，宏微观摩擦接触模型与试验数据的温差仅为4.29 ℃，远远低于单一宏观模型的18.24 ℃。

3 温度场失效概率统计模型

通过宏微观摩擦接触模型构建湿式摩擦元件温度场失效概率统计模型（temperature-field failure probability model， TFPM）。以转速、压力、温度和滑摩时间为失效特征参数构建极限状态函数，通过蒙特卡罗模拟计算得到摩擦界面的失效概率，并通过试验数据验证TFPM的准确性。该模型的结构如图7所示。

3.1 核密度函数

研究失效概率时，首先需建立失效特征参数的概率密度函数。湿式摩擦元件的工作环境复杂多变，不同工况的失效特征参数变化规律并不遵循特定分布，因此，本研究采用核密度估计法构建失效特征参数的概率密度函数，以作为失效概率研究的可靠基础^［20］。

假设一维数据有n个样本，则该样本数据的累积分布函数为

（8）

f(x)=\underset{h\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{F(x+h_{\mathrm{d}})-F(x-h_{\mathrm{d}})}{\mathrm{2}{h}_{\mathrm{d}}}（9）

式中：f（x）为概率密度函数；x为失效特征参数；i为第i个样本；h{}_{\mathrm{d}}为带宽。

为简化表达，定义函数

{\mathrm{1}}_{y_i}=\left\{\begin{array}{cc}\mathrm{1}& y_i\mathrm{成立}\\ \mathrm{0}& y_i\mathrm{不成}\mathrm{立}\end{array}\right.

式中：y_i 为数学表达式。

引入累积分布函数的经验分布函数

F_n(T)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathrm{1}}_{x_i(T)\le T}（10）

式中：T为失效边界；n为样本总数。

使用n次观测中x_i ≤T出现的次数与n的比值近似描述P（x≤T）。将式（10）代入式（9）可得

f(x_i)=\underset{h_{\mathrm{d}}\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{1}}{\mathrm{2}n{h}_{\mathrm{d}}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathrm{1}}_{x_i-h_{\mathrm{d}}\le x_j\le x_i+h_{\mathrm{d}}}（11）

在f（x_i ）的计算过程中，必须给定h_d（核密度估计中的带宽）。确定带宽后，f（x）的表达式为

f(x)=\frac{\mathrm{1}}{\mathrm{2}n{h}_{\mathrm{d}}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathrm{1}}_{x-h_{\mathrm{d}}\le x_i\le x+h_{\mathrm{d}}}=\frac{\mathrm{1}}{\mathrm{2}n{h}_{\mathrm{d}}}{\displaystyle \sum _{i=\mathrm{1}}^n}K(\frac{|x-x_i|}{h_{\mathrm{d}}})（12）

对等式两边同时积分可以得到

\int f(x)\mathrm{d}x=\int \frac{\mathrm{1}}{\mathrm{2}}K(T)\mathrm{d}T（13）

令K₀（T）=0.5K（T），根据概率密度函数的定义，可以得到

\int K_{\mathrm{0}}(T)\mathrm{d}T=\mathrm{1}（14）

其中，K₀（T）称为核函数。0≤T≤1时，K₀（T）=0.5。

将K₀（T）代入式（12），此时f（x）的表达式为

f(x)=\frac{\mathrm{1}}{n{h}_{\mathrm{d}}}{\displaystyle \sum _{i=\mathrm{1}}^n}{K}_{\mathrm{0}}(\frac{|x-x_i|}{h_{\mathrm{d}}})（15）

对于二维数据，f（x， y）可以表示为

f(x,y)=\frac{\mathrm{1}}{n{h}_{\mathrm{d}}^{\mathrm{2}}}{\displaystyle \sum _{i=\mathrm{1}}^n}{K}_{\mathrm{0}}(\frac{D((x,y),(x_i,y_i))}{h_{\mathrm{d}}})（16）

式中：D（A，B）为A、B两点之间的距离。

以此可以得到任意数据的概率密度函数。

3.2 极限状态函数

为评估摩擦界面的温度场失效概率，需要确定摩擦界面的失效边界。失效判据的数学定义对材料或结构的寿命预测和可靠性评估至关重要^［21］。失效边界的判定公式^［22］为

Z=R(x)-S(x)（17）

R(x)=R=\mathrm{400}\mathrm{M}\mathrm{P}\mathrm{a}（18）

S(x)=S=\frac{E\alpha }{r^{\mathrm{2}}}(\frac{r^{\mathrm{2}}-r_{\mathrm{i}\mathrm{n}}^{\mathrm{2}}}{r_{\mathrm{o}\mathrm{u}\mathrm{t}}^{\mathrm{2}}-r_{\mathrm{i}\mathrm{n}}^{\mathrm{2}}}{\int }_{r_{\mathrm{i}\mathrm{n}}}^{r_{\mathrm{o}\mathrm{u}\mathrm{t}}}\mathrm{\Delta }Tr\mathrm{d}r-{\int }_{r_{\mathrm{i}\mathrm{n}}}^{r_{}}\mathrm{\Delta }Tr\mathrm{d}r)（19）

式中：R（x）、S（x）分别为当前区域的许用热应力和当前热应力；x为失效特征参数，x=x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_m；α为热膨胀系数；E为弹性模量；r_in、r_out分别为钢片的内径和外径；ΔT为温度梯度函数。

由式（17）可以判断出，Z>0表示摩擦界面处于安全状态，Z=0表示摩擦界面处于极限状态，Z0表示摩擦界面处于失效区域。

利用符号回归方法结合工况参数拟合ΔT的计算公式^［23］：

\begin{array}{c}\mathrm{\Delta }T=\mathrm{17.604}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{2}{X}_{\mathrm{3}}+X_{\mathrm{1}}+\mathrm{4.335}}{\mathrm{2.491}}+(X_{\mathrm{2}}+X_{\mathrm{4}})+\\ (\mathrm{l}\mathrm{n}(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{1.885})+\frac{\mathrm{6.920}+X_{\mathrm{1}}}{\sqrt[]{\mathrm{4.245}}})\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{l}\mathrm{n}\frac{\mathrm{s}\mathrm{i}\mathrm{n}(-\mathrm{6.920}-X_{\mathrm{1}})}{X_{\mathrm{3}}})\end{array}（20）

式中：X₁为转速；X₂为接合时间；X₃为半径；X₄为压力。

通过极限状态函数可确定每个微观区域是否处于失效状态。通过核密度函数可得Z的概率密度函数：

f_z(z)=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}{\sigma }_z}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{(z-{\mu }_z{)}^{\mathrm{2}}}{\mathrm{4}{\sigma }_z^{\mathrm{2}}})（21）

式中：σ{}_z为标准差；μ{}_z为均值。

摩擦界面的失效概率为

P_{\mathrm{f}}={\int }_{-\mathrm{\infty }}^{\mathrm{0}}{f}_Z(z)\mathrm{d}z=

\int \cdots {\int }_{\Omega }{f}_X(x_{\mathrm{1}},x_{\mathrm{2}}\cdots x_m)\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{2}}\cdots \mathrm{d}{x}_m（22）

式中：f_x （x₁，x₂，…，x_m ）为x的联合概率密度函数。

摩擦界面处于失效区域或临界区域时，Z≤0，因此积分区间Ω为（-∞，0）。就摩擦界面而言，可假设各失效特征参数是相互独立的。因此，联合概率密度函数可用单个密度函数代替，即

P_{\mathrm{f}}={\int }_{\Omega }{f}_X(x)\mathrm{d}x=

\int ...{\int }_{\Omega }{f}_{X_{\mathrm{1}}}(x_{\mathrm{1}})f_{X_{\mathrm{2}}}(x_{\mathrm{2}})...,f_{X_m}(x_m)\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{2}}\dots \mathrm{d}{x}_m（23）

式中：f_Xi （x_i ）为随机变量x_i 的概率密度函数。

一般来说，P_f的积分复杂且难以求解，因此本文引入蒙特卡洛模拟方法计算这个多重积分。

3.3 蒙特卡罗模拟

蒙特卡罗模拟基于各随机变量的概率分布生成样本，并通过计算系统的响应来评估其行为。其中，失效概率由失效样本的数量与总样本数之比进行估计，即通过统计失效频率来实现。具体而言，该方法首先依据概率密度函数f_x （x）抽取随机样本x。如果总采样数为N，则将采集到的样本值依次代入式（17）。如果Z<0，则模拟中发生了结构失效，失效次数N_f加1。根据伯努利大数定律，在N次独立重复试验中，随机事件Z<0的频率N_f /N收敛于发生概率P_f。因此，失效概率的预测值为^［24］

{\widehat{P}}_{\mathrm{f}}=N_{\mathrm{f}}/N（24）

根据式（22），界面的失效概率为

P_{\mathrm{f}}={\int }_{{\mathrm{\Omega }}_{\mathrm{f}}}{f}_X(x)\mathrm{d}x=

{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}I[g_X(x)]f_X(x)\mathrm{d}x=E\{I[g_X(x)]\}（25）

（26）

利用指标函数I（α）保证Z>0时的积分为0，从而将积分域扩展到整个定义域。根据式（23），失效概率预测值为

{\widehat{P}}_{\mathrm{f}}=\frac{\mathrm{1}}{N}\{I[g_X(x_i)]\}（27）

式中：I[g_X(x_i)]为指标函数生成的样本。

根据上述蒙特卡罗模拟原理，建立温度场失效概率统计模型。假设有m个随机变量 x=（x₁，x₂，…，x_m ）^T，且所有随机变量互相独立，则通过随机变量的统计分布产生N个不同向量：

{\mathit{X}}_i=(X_{i\mathrm{1}},X_{i\mathrm{2}},\cdots ,X_{im})i=1，2，…，N（28）

向量 X_i 代表一个样本，假设失效区域F=｛x：Z≤0｝，落在失效区域的样本数为N_f，则失效概率P_f=N_f/N。

蒙特卡罗模拟基于各参数的概率密度分布，通过大量随机抽样，统计系统失效的样本比例，从而近似估计微观失效概率。该方法依赖大量局部失效数据的统计分析，以得出系统整体的失效可能性。

3.4 模型验证

为验证温度场失效概率统计模型的准确性，引入标准差σ、微观失效概率的均值μ、失效概率大于i的区域的面积S_i 与整个摩擦界面的面积S之比β_i。

σ可体现摩擦界面径向失效概率的差异。设每条半径上所有P_f的平均值为x_i，计算所有半径对应的x_i 之间的σ。σ的计算公式为

\sigma =\sqrt[]{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}(x_i{-a)}^{\mathrm{2}}}（29）

式中：x_i 为第i条半径上所有P_f 的平均值；a为x_i 的均值。

摩擦元件整个界面所有微观失效概率的均值为

\mu =\frac{\mathrm{1}}{M_{\mathrm{1}}{N}_{\mathrm{1}}}{\displaystyle \sum _{i=\mathrm{1}}^{M_{\mathrm{1}}}}{\displaystyle \sum _{j=\mathrm{1}}^{N_{\mathrm{1}}}}{a_i}_j（30）

式中：a_ij 为第i行、第j列的失效概率；M₁为总行数；N₁为总列数。

通过多个工况的温度来验证TFPM。通过宏微观摩擦接触模型对摩擦元件在转速500 r/min、压力0.5 MPa工况下进行全寿命周期仿真，获取湿式摩擦元件在2、4、6、8 h的温度场，并计算得到不同时间节点下摩擦元件的微观失效概率分布图。图8所示为四个时间节点下μ和真实退化曲线的对比。图9所示为4个时间节点的σ和β_0.5。

由图8可看出，2、4、6、8 h的μ分别为0.254、0.505、0.759和0.951。摩擦元件的失效概率大体呈线性增长。试验过程中没有较大的冲击载荷，且转速和压力也为恒定值，因此可认为试验过程中的摩擦元件的退化趋势为线性退化^［25］，与TFPM的计算结果一致。图9中，β_0.5上升而σ下降。结合图8易发现，随着滑摩时间的延长，摩擦元件表面温度升高，摩擦元件的径向温差逐渐减小。中径偏外区域（中外径）的温度显然最高，而随着热量的积累，摩擦元件内径区域和外径区域的温度逐渐升高，导致摩擦元件径向上的温度差异降低，与试验结果一致，验证了TFPM具有较高的准确性和合理性。根据温度场数据分析得出，β_0.5=0.9时，摩擦元件表面产生较大的温差，表面热应力较大，摩擦因数的稳定性剧烈下降，摩擦元件处于失效边界。不同工况下的摩擦元件温度分布不同，因此不同工况的失效边界对应的β_i 及i的大小也应根据具体工况数据分析。

3.5 不同模型对比

为评估不同模型的失效概率计算结果，选取均方根误差RMSE和相关指数R²作为评价指标，相关数学公式可以表示为^［26］

I_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt[]{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}(X_{\mathrm{p}\mathrm{r}\mathrm{e},i}-X_{\mathrm{o}\mathrm{b}\mathrm{s},i}{)}^{\mathrm{2}}}（31）

I_{R^{\mathrm{2}}}=\mathrm{1}-{\displaystyle \sum _{i=\mathrm{1}}^n}(X_{\mathrm{p}\mathrm{r}\mathrm{e},i}-X_{\mathrm{o}\mathrm{b}\mathrm{s},i}{)}^{\mathrm{2}}/{\displaystyle \sum _{i=\mathrm{1}}^n}(\overline{X}-X_{\mathrm{o}\mathrm{b}\mathrm{s},i}{)}^{\mathrm{2}}（32）

式中：X_pre，i 为预测值；X_obs，i 为真实值；\overline{X}为真实值的均值。

图10a所示为基于宏微观摩擦接触模型的TFPM（I-TFPM）、基于单一宏观模型的TFPM（S-TFPM）与另外9种神经网络模型μ的计算结果。由图10a可看出，I-TFPM的预测精度最接近线性退化曲线，单一宏观模型得出的失效特征参数用于失效概率的计算误差较大。这是由于宏微观摩擦接触模型仿真得到的温度更接近试验数据。这进一步验证了宏微观摩擦接触模型对温度场失效研究的必要性。相较于另外9种深度学习模型，I-TFPM显然具有更好的拟合效果，对温度场失效概率的计算能力更加优越，与试验结果更接近。

4 影响因素分析

考虑转速和压力对摩擦元件温度场失效的影响，利用μ、σ、β_0.5分析不同工况下的摩擦元件温度场分布规律，以验证I-TFPM的分析结果和泛化能力。

4.1 转速

研究转速对摩擦元件温度场失效的影响，取压力为0.65 MPa，转速分别为300、800、1300、1800 r/min，对比四组工况下摩擦元件的温度场和失效概率，结果如图11、图12所示。

由图11可看出，转速增大时，σ、μ和β_0.5都有着不同程度的增大；随着转速的增大，微凸体温度升高。这是由于随着转速的增大，钢片和摩擦片之间的相对滑动速度增大，摩擦元件表面微凸体之间的摩擦接触更加紧密，摩擦界面温升加快，使摩擦界面高温区分布更广，中径区域无法有效散热，导致径向温差和周向温差增大，可能加剧摩擦元件产生热翘曲变形，使摩擦元件产生失效现象。

较高的转速影响润滑油的流动和表面压力的分布，导致表面温度急剧升高，而中径的散热条件要远远差于内径和外径，因此中径区域的热量难以散发，热应力也随之增大，进而导致中径的失效概率增大。由图12可以看出，摩擦界面温度场的失效概率沿径向呈不对称分布，其中中径区域最高，而内外径区域相对较低。这主要是因为，在摩擦过程中，热流密度随半径增大而升高，导致内径处的热流密度最低，因而失效概率也最低。对于外径区域，尽管其热流密度最高，但由于其散热条件优于中径处，更多的热量被及时耗散，故其失效概率也低于中径处。摩擦元件不同径向位置的材料热导率不同。中径区域由于材料结构和冷却条件不同，无法有效散热，导致热量在这一位置积累。内径和外径的散热条件相对中径处更优，导致这两个位置发生温度场失效的概率相对较低。

4.2 压力

研究压力对摩擦元件温度场失效的影响，转速850 r/min，压力分别为0.2、0.4、0.7、1.0 MPa的摩擦元件温度场和失效概率如图13、图14所示。

由图13可看出，压力增大时，σ、μ和β_0.5都有着不同程度的增大。随着转速的增大，微凸体的温度升高。这是由于随着压力的增大，摩擦片和钢片表面的微凸体产生弹塑性变形，摩擦副之间的空隙减小，表面微凸体之间的接触更紧密，增大了实际接触面积加剧了摩擦元件之间的摩擦和磨损，微凸体之间的摩擦加剧，导致界面温度上升。中径区域的散热条件远差于内径和外径，因此热量在中径位置大量积累，导致中径处的径向温差和周向温差增大，加剧摩擦元件的热翘曲变形，从而改变其表面的热流分布与接触状态，最终导致温度场失效。

较大的压力会破坏润滑油膜，降低润滑冷却效果。高压力下，摩擦元件的材料发生热膨胀，导致形变或接触不均匀，加剧局部热量的累积。这种不均匀的温度分布可能会加速材料的损伤和失效。由图14可看出，摩擦界面内径和外径的温度场失效概率较低，中外径的温度场失效概率较高。湿式摩擦元件中，摩擦力通常在接触区域的分布是不均匀的。较大的摩擦力会产生更多的摩擦热，使中径的温度较高。

5 结论

1）宏微观摩擦接触模型仿真结果与试验结果吻合度高，模型有效且具有很高的准确性。10次循环交互中，宏微观摩擦接触模型仿真效果始终优于单一宏观模型，随着交互次数的增加，宏微观摩擦接触模型的仿真结果不断逼近试验数据。10次交互完成后，宏微观摩擦接触模型与试验数据的温差仅为4.29 ℃，远远低于单一宏观模型的18.24 ℃。

2）基于宏微观摩擦接触模型的TFPM（I-TFPM）能精确评估湿式摩擦元件的温度场失效概率。I-TFPM的失效概率计算结果与试验数据基本一致，大体呈线性增长趋势。取摩擦元件全寿命周期中的2 h、4 h、6 h和8 h，摩擦元件的微观失效概率的均值分别为0.254、0.505、0.759和0.951。I-TFPM的相关系数和均方根误差分别为0.991和0.025，比其他10种模型的预测精度高。

3）摩擦元件的温度场失效概率随转速或压力的增大而增大，300、800、1300、1800 r/min下微观失效概率均值分别为0.279、0.307、0.326、0.359，0.2、0.4、0.7、1 MPa时微观失效概率均值分别为0.274、0.306、0.338、0.363。径向失效概率标准差和失效概率大于0.5区域的面积占总面积的比值随转速和压力的增大而增大。

参考文献

[1] 彭程， 陈俐， 傅圣来. 基于离合器滑摩扭矩估计的并联式气电混合动力船舶模式切换控制［J］. 中国机械工程， 2024， 35（8）： 1489-1497.

[2] PENG Cheng， CHEN Li， FU Shenglai. Mode Transition Control of Parallel Gas-electric Hybrid Ships Based on Clutch Slipping Torque Estimation［J］. China Mechanical Engineering， 2024， 35（8）： 1489-1497.

[3] WANG Baolin， MAI Y W. Transient One-dimensional Heat Conduction Problems Solved by Finite Element［J］. International Journal of Mechanical Sciences， 2005， 47（2）： 303-317.

[4] LIU Yuwei， SUN Yuanzhi， GAO Ziyin， et al. Investigation of Thermal Behavior of a Novel Multi-conical Friction Pair during Wet Clutch's Engagement Based on Finite Difference Method［J］. International Communications in Heat and Mass Transfer， 2024， 155： 107555.

[5] ZHAO Qin， MA Biao， YU Liang， et al. Numerical and Experimental Studies on the Thermodynamic Characteristics of Post-buckling Separate Plate in the Clutch［J］. Tribology International， 2024， 195： 109607.

[6] 李和言， 王宇森， 熊涔博， 等. 离合器配对摩擦副径向温度梯度对接触比压的影响［J］. 机械工程学报， 2018， 54（1）： 136-143.

[7] LI Heyan， WANG Yusen， XIONG Cenbo， et al. Effect of Radial Temperature Gradient on Interface Pressure of Clutch Friction Pair［J］. Journal of Mechanical Engineering， 2018， 54（1）： 136-143.

[8] 王立勇， 吴瑾， 李乐， 等. 径向非均布压力分布对湿式摩擦副热-机耦合影响［J］. 北京理工大学学报， 2021， 41（6）： 588-596.

[9] WANG Liyong， WU Jin， LI Le， et al. Influence of Radial Non-uniform Pressure Distribution on Thermo Mechanical Coupling of Wet Friction Pairs［J］. Transactions of Beijing Institute of Technology， 2021， 41（6）： 588-596.

[10] 郑良杰， 马彪， 陈漫， 等. 花键摩擦对湿式多片离合器分离过程影响［J］. 哈尔滨工业大学学报， 2022， 54（7）： 89-95.

[11] ZHENG Liangjie， MA Biao， CHEN Man， et al. Influence of Spline Friction Coefficient on a Wet Multi-disc Clutch in the Disengaging Process［J］. Journal of Harbin Institute of Technology， 2022， 54（7）： 89-95.

[12] SHIN S， CHOI S B. Control-oriented Modeling of Wet Clutch Friction Considering Thermal Dynamics［J］. Mechatronics， 2024， 99： 103146.

[13] ABDULLAH O I， SCHLATTMANN J， MAJEED M H， et al. The Temperatures Distributions of a Single-disc Clutches Using Heat Partitioning and Total Heat Generated Approaches［J］. Case Studies in Thermal Engineering， 2018， 11： 43-54.

[14] FURMAŃSKI P， ŁAPKA P. Micro-macro Heat Conduction Model for the Prediction of Local， Transient Temperature in Composite Media［J］. International Journal of Thermal Sciences， 2020， 154： 106401.

[15] MBEN DAVID， ZIMMELS Y， ZVIRIN Y. Determination of the Quench Velocity and Rewetting Temperature of Hot Surfaces. Part II： an Integrated Approach of Hydrodynamic-Thermodynamic Micro- and Macro-scale Models［J］. International Journal of Heat and Mass Transfer， 2001， 44（7）： 1343-1371.

[16] 吴健鹏， 马彪， 李和言， 等. 加速寿命过程的湿式摩擦副滑摩温升特性［J］. 哈尔滨工业大学学报， 2019， 51（7）： 103-111.

[17] WU Jianpeng， MA Biao， LI Heyan， et al. Temperature Rise of Wet Friction Pair during Sliding in the Accelerated Life Test［J］. Journal of Harbin Institute of Technology， 2019， 51（7）： 103-111.

[18] 李乐， 舒越超， 吴健鹏， 等. 基于PSO-BP神经网络湿式摩擦元件损伤预测模型［J］. 北京理工大学学报， 2022， 42（12）： 1246-1255.

[19] LI Le， SHU Yuechao， WU Jianpeng， et al. A Damage Prediction Model of Wet Friction Elements Based on PSO-BP Neural Network［J］. Transactions of Beijing Institute of Technology， 2022， 42（12）： 1246-1255.

[20] 蔡运迪， 唐文献， 黄秋芸， 等. 水冷盘式制动器热疲劳失效有限元分析［J］. 中国机械工程， 2012， 23（22）： 2726-2731.

[21] CAI Yundi， TANG Wenxian， HUANG Qiuyun， et al. Finite Element Analysis（FEM） on Thermal Fatigue Failure of Water-cooling Disc Brake［J］. China Mechanical Engineering， 2012， 23（22）： 2726-2731.

[22] GANESH K C， DHANASEKARAN K， SARAVANAN P， et al. Failure Analysis of Multi Plate Wet-clutch Used in Off-highway Transmission［J］. Engineering Failure Analysis， 2021， 127： 105534.

[23] WANG Shaohua， HUANG Huanming， SHI Dehua， et al. Hybrid Propagation Modeling Based Clutch Fault Diagnosis of Multi-mode Electromechanical Transmission System Using Particle Filter［J］. Measurement， 2025， 243： 116385.

[24] STOLARSKI T A. A System for Wear Prediction in Lubricated Sliding Contacts［J］. Lubrication Science， 1996， 8（4）： 315-351.

[25] KOGUT L， ETSION I. A Finite Element Based Elastic-plastic Model for the Contact of Rough Surfaces［J］. Tribology Transactions， 2003， 46（3）： 383-390.

[26] ZHANG Xin， ZHANG Lei， YANG Zhi， et al. Tracing the Evolution of Contact Area during the Running-in Process of AgCu-MoS₂ Composite［J］. Tribology International， 2022， 174： 107739.

[27] CHEN Jiaqi， HE Yulin， CHENG Yingchao， et al. A Multiple Kernel-based Kernel Density Estimator for Multimodal Probability Density Functions［J］. Engineering Applications of Artificial Intelligence， 2024， 132： 107979.

[28] GUO Jinwei， YANG Li， ZHOU Yichun， et al. Reliability Assessment on Interfacial Failure of Thermal Barrier Coatings［J］. Acta Mechanica Sinica， 2016， 32（5）： 915-924.

[29] XIONG Cenbo， MA Biao， LI Heyan， et al. Experimental Study and Thermal Analysis on the Buckling of Friction Components in Multi-disc Clutch［J］. Journal of Thermal Stresses， 2015， 38（11）： 1323-1343.

[30] YU Liang， MA Biao， LI Heyan， et al. Numerical and Experimental Studies of a Wet Multidisc Clutch on Temperature and Stress Fields Excited by the Concentrated Load［J］. Tribology Transactions， 2019， 62（1）： 8-21.

[31] XIAO Y Q， LIU Z Y， ZHU W， et al. Reliability Assessment and Lifetime Prediction of TBCS on Gas Turbine Blades Considering Thermal Mismatch and Interfacial Oxidation［J］. Surface and Coatings Technology， 2021， 423： 127572.

[32] WU Jianpeng， LI Pengpeng， GUO Liuyang， et al. Fuzzy Comprehensive Evaluation of Hydroviscous Drive Performance Degradation Based on Improved Game Theory Combination Weighting［J］. IEEE Sensors Journal， 2024， 24（15）： 23575-23584.

[33] 刘会永， 张松， 李剑峰， 等. 采用改进CNN-BiLSTM模型的刀具磨损状态监测［J］. 中国机械工程， 2022， 33（16）： 1940-1947.

[34] LIU Huiyong， ZHANG Song， LI Jianfeng， et al. Tool Wear Detection Based on Improved CNN-BiLSTM Model［J］. China Mechanical Engineering， 2022， 33（16）： 1940-1947.

基金资助

国家自然科学基金(52105084)

深圳市基础研究专项（自然科学基金）(JCYJ20240813113329038)