动车组闸片磨损分析及寿命预测

张辉,李斌,史立,尹剑,张生芳

(1.天津铁路信号有限责任公司,天津 300300;2.大连交通大学 机械工程学院,辽宁 大连 116028)

随着制动技术的发展,盘式制动器以其摩擦稳定性高、制动噪声小等众多优点,逐渐取代了摩擦制动方式中的闸瓦制动器,成为动车制动系统的重要形式[1]。盘式制动器在工作过程中需要在短时间内实现能量转换,产生的摩擦热会使制动盘/片处于高温高压环境中,在高温高压作用下,制动盘/片内会产生热应力、压应力以及微弱的变形。同时,由于摩擦作用,盘式制动摩擦副会产生不同程度的磨损。磨损的累积会改变摩擦副工作过程中的接触面积,进而影响制动稳定性,降低动车运行安全性。在动车组运维过程中,针对闸片的疲劳寿命监测多为预防性检修,需要经常对动车组闸片进行拆卸、检查、维护,因此需要对动车组闸片制动过程中的热、力、磨损及其寿命预测进行相关研究。

王磊等[2]利用Ansys建立了动车组用闸片循环对称模型,分析了闸片几何尺寸对热流密度分配、对流散热和辐射散热的影响。夏毅敏等[3]基于传热理论和有限元方法,考虑温度对材料特性的影响,利用Ansys分析了闸片温度场和应力场在制动过程的分布以及变化规律。张生芳等[4]用间接耦合法将热分析结果代入结构场,分析闸片的温度及应力变化。

摩擦必然伴随着磨损,磨损的计算方法种类繁多。Rajesh等[5]建立了圆盘与销轴接触分析的有限元模型,用3D表面到表面的体-地连接和旋转类型表示接触面的滑动行为。张方宇等[6]模拟了摩擦片和制动盘的磨合过程,证明了应力-磨损耦合分析方法的可行性。郝强[7]使用微元法改进了Archard磨粒磨损计算公式,利用有限元法计算探讨了宏观尺度下的闸片表面磨损分布规律。尹家宝等[8]提出一种求解摩擦块磨损的数值方法,利用Abaqus结合磨损子程序实现了磨损量及磨损量累计过程的模拟,最终通过Matlab二次开发直接提取了磨损数值。赵子衡等[9]对传统的Archard磨损模型进行改进,模拟了某鼓式制动器制动过程,结果表明完全耦合法能够更真实地模拟制动过程的热-力-磨损耦合。

基于统计分析的寿命预测主要通过大量采集初始数据建立统计模型,从而实现寿命的预测。曲线拟合、灰色预测和神经网络等方法能够很好地实现磨损预测。Zhang等[10]提出一种基于深度学习的寿命预测方法,组建了混合学习结构网络,通过多项式回归模型计算预测。Chang等[11]利用磨损模型和磨损进程的递归函数,通过原始算法推理了磨损随时间的变化。Fan等[12]利用遗传算法优化了长短时记忆网络,利用训练好的网络进行电机轴承寿命的预测。

本文基于热-力耦合有限元模型,通过Archard磨损理论进行二次开发,实现了闸片磨损的有限元计算,模拟分析了闸片的磨损演变分布规律。采用单一变量法研究了制动参数对闸片磨损的影响,并根据闸片磨损仿真,建立磨损与制动参数对应关系的数据库,通过BP算法调取数据建立磨损预测模型,建立动车组闸片磨损量与行驶里程的对应关系,最后设计了动车组闸片寿命预测系统。

1 制动盘热-力耦合磨损计算模型

1.1 热传导初始条件及边界条件

根据传热学理论,摩擦热会以热传导、热对流以及热辐射的形式传递到外部环境。

(1)热传导模型

将闸片看作无内热源的各向同性材料,根据传热理论,其热传导方程为:

(1)

式中:T为温度,K;t为时间,s;ρ为材料密度,kg/m3;C为材料比热容,J/(kg·K);λ为材料导热系数,W/(m·K)。

动车组盘式制动器主要有热传导、热对流及热辐射,模型初始边界条件为:在制动初始时刻t=0时,T(x,y,z)=T0,初始温度T0取20 ℃。

在闸片的所有换热界面上:

(2)

在摩擦副接触摩擦界面上:

(3)

式中:T0为环境温度;ni为各界面的法向单位向量;z为摩擦界面的法向单位向量;αi为各界面的对流换热系数;αz为摩擦界面的对流换热系数;εi、εz为辐射换热系数,W/(m2·K);σ为斯蒂芬-波尔兹曼常数;q为闸片吸收的热流密度,W/m2。

(2)对流换热模型

设强制对流换热系数为αc,对于普朗特数Pr>0.6的流体,当雷诺数Rel

(4)

当雷诺数Rel>Rexc时,对流换热系数为:

(5)

式中:Re为雷诺数;Pr为空气的普朗特数,取0.7;λ为空气的热导率,取λ=6.14×10-6W/(m·K);l为闸片的特征长度。

(3)辐射换热系数

制动过程中存在辐射换热现象。根据斯蒂芬-波尔兹曼定律,闸片与空气之间的辐射换热关系式为:

(6)

式中:ε为物体发射率;σ0为斯蒂芬-玻尔兹曼系数,σ0=5.67×10-8W/(m2·k4)。

1.2 动车组盘式制动器有限元模型的建立

动车组盘式制动系统在工作过程中发生了大量的能量转换,能量的转换又包含多种形式,根据研究条件,提出如下假设:

a.闸片及制动盘为各向同性材料,摩擦系数恒定为0.37。

b.模拟制动过程是匀减速。

c.忽略摩擦过程中磨屑带走的能量。

d.制动压力在闸片背面均布分布。

e.忽略制动过程中环境温度的变化,初始值设置为20 ℃。

f.忽略振动对结果的影响。

制动盘外半径为320 mm,内半径为180 mm,制动盘原型厚度为80 mm,制动盘厚度为40 mm;闸片外半径为310 mm,内半径为188 mm,厚度为20 mm,摩擦面积为28 500 mm2。

(1)摩擦副材料参数

摩擦副材料初始参数见表1,制动盘材料参数随温度变化见表2。

表1 摩擦副材料初始参数[4]

(2)有限元模型前处理

制动盘和闸片采用六面体网格,制动盘单元数量为2 400个,闸片单元数量为540个,采用减缩积分单元C3D8T,见图1。

图1 制动盘和闸片网格模型

表面热交换条件膜层散热系数为0.001 W/(m·K),热辐射发射率设置为0.2。图2为有限元模型边界条件示意图,闸片四周表面受x、y方向的转动和位移约束,RP点约束在x、y方向的位移,施加转速。制动盘非接触表面施加z方向的对称约束。

图2 有限元模型边界条件示意图

1.3 磨损计算方法

Archard的磨损模型[13]形式如下:

(7)

式中:V为磨损体积;s为磨损位移;k为无量纲的磨损系数;FN为法向载荷;H为磨损表面的硬度。

将闸片接触面划分为N个无穷小的接触面积,在单一无穷小的接触面积ΔA上,时间增量为dt,则磨损模型演变的微分形式为:

(8)

式中:dV为磨损体积增量;ds为滑动位移增量;FN为法向接触力。

假设ΔA在dt下的磨损厚度增量为dh,则dV=ΔAdh,代入式(7)得:

(9)

式中:FN/ΔA为无穷小接触面积处的接触应力,令FN/ΔA为pc;k/H为有量纲磨损系数,用kd表示,则式(9)表示为:

dh=kdpcΔs

(10)

取时间增量为Δt,磨损厚度增量为Δh,相对滑动位移增量为Δs,得:

Δh=kdpcΔs

(11)

首先计算接触应力pc。当Δs较小时,认为pc为不变量,求得对应的Δh,根据h变化之后的模型,获得接触应力pc。

设闸片和制动盘间接触点i=1,…,N,欧拉积分步数j=1,…,M,第i个接触点第j个积分步的接触应力为pc(i,j),磨损深度为Δh(i,j),j次磨损深度为h(i,j),制动盘转动位移为Δs(i,j)。则:

Δh(i,j)=kdpc(i,j)Δs(i,j)

(12)

h(i,j+1)=h(i,j)+Δh(i,j)

(13)

在ALE自适应网格选项中设置磨损边界为速度约束,将式(14)左右两侧同时除以第j个时间步长Δt(j),得:

(14)

(15)

2 闸片热-力耦合及磨损结果分析

2.1 闸片瞬态温度场应力场仿真结果分析

本节模拟计算制动初速度为120 km/h、制动压力为23 000 N时的制动工况。图3为制动过程中闸片温度场分布云图。制动开始后,摩擦产热远大于摩擦散热,闸片表面的温度不断上升,随着摩擦的持续,摩擦表面温度分布有明显的变化。在1 s时摩擦面靠近闸片外径处开始出现温度梯度;在13 s左右,由于转速不断降低,摩擦热流的输入与热量的耗散达到平衡,达到了制动过程中的最高温度;之后摩擦产热小于热量的耗散,温度开始下降。整个制动过程中,闸片表面最高温度达到了209.4 ℃,制动结束时最高温度为122.6 ℃,闸片的最高温度出现在摩擦半径273 mm处。

图3 制动过程中闸片温度分布云图

图4为制动过程中闸片应力分布云图。在制动过程中,由于摩擦作用,闸片摩擦表面累积了大量的热量,在摩擦表面存在相应的温度梯度,高温使得闸片的变形向径向扩展,闸片摩擦区域同时产生很大的热应力。制动初期,由于摩擦入口先发生滑动,闸片摩擦入口区域产生局部热应力;摩擦制动持续进行之后,摩擦出口区域应力逐渐大于其余位置,闸片径向的温度梯度同时影响闸片应力分布,在摩擦半径273 mm附近出现了最大应力。从制动开始13 s后闸片摩擦表面的径向应力分布呈明显的梯度分布,由于应力集中的存在,周向应力的分布存在较大的差异。

图4 制动过程中闸片应力分布云图

2.2 闸片磨损仿真结果分析

本节将基于上述有限元模型,设置制动初速度为160 km/h、制动压力为14 000 N、仿真时长为10 s,不同时刻闸片磨损分布云图见图5。闸片左侧为摩擦入口区域,右侧为摩擦出口区域。在制动初期,压力大、速度快,磨损值增加迅速,闸片中间部分磨损最小,随后在温度以及应力共同作用下,闸片的上下边缘部分磨损程度最小。制动初期摩擦入口区域磨损大,由于摩擦方向向右,随着相对运动的持续,逐渐演变成摩擦出口区域的磨损大。

图5 不同时刻闸片磨损分布云图

选择闸片不同周向以及径向节点,周向节点从左到右编号1～19,径向节点从下到上编号1～11,见图6。闸片表面磨损规律见图7, 闸片径向磨损呈

图6 闸片表面磨损分析节点图

(a) 周向磨损

(b) 径向磨损图7 闸片表面磨损规律

现先增大后减小的趋势,在制动过程中闸片内径的线速度相对较小,因此节点1磨损值最小,随着各节点摩擦半径的增大,磨损程度相比于内侧均有所增大。制动过程中温度及应力最大值区域位于摩擦半径273 mm(节点8)附近,温度与应力大小与磨损深度成正比,因此节点8处磨损最大,节点8向两侧呈现磨损梯度值。

2.3 制动初速度对闸片磨损的影响

将初始压力设定为23 kN,随着制动初速度的增大,闸片表面磨损深度不断增加,但闸片摩擦表面磨损分布趋势无明显变化。速度从80到180 km/h,最大磨损深度为1.40×10-2、1.69×10-2、2.01×10-2、2.27×10-2、2.53×10-2、2.80×10-2mm。制动过程中,初始速度直接决定摩擦副之间的相对滑动速度,在相同时间、相同制动压力下,初速度越大,摩擦副表面相对滑动距离越大,导致闸片表面磨损加剧,不同制动初速度下闸片磨损曲线见图8。

图8 不同制动初速度下闸片磨损曲线

从图8可以看出,在制动初期,不同转速下最大磨损量差值较小,随着制动时间的增加,其差值逐渐增大,且差值存在一定的差异。以80 km/h为起点,20 km/h为一个跨度,速度每增加20 km/h,磨损深度分别增加2.91×10-3、3.19×10-3、2.66×10-3、2.52×10-3、2.78×10-3mm。如图9所示,随着速度的增加,磨损深度差值的变化呈现先上升后下降再上升的趋势,但磨损深度的变化率逐渐减小,后期趋于平稳。速度增大,磨损深度变化较大,速度达到120 km/h以后,摩擦副温度较大。同时,由于闸片节点磨损,摩擦副接触面积发生变化,摩擦力大小发生一定的变化,导致磨损深度增加量减小,因此磨损变化率逐渐减小。由于速度不断上升,摩擦副相对滑动速度逐渐变大,滑动速度的增大将会加剧闸片接触表面的磨损,磨损深度进一步增大。

(a) 磨损深度差值

(b) 磨损变化率图9 磨损变化曲线

2.4 制动压力对闸片磨损的影响

将初始速度设定为120 km/h,闸片摩擦表面磨损分布趋势随着制动压力的增大无明显变化。图10为不同制动压力下闸片磨损曲线,不同压力下,磨损深度的变化趋势相同,摩擦面摩擦力随着制动压力的增大而增大,进一步增大摩擦表面接触应力,迫使闸片表面发生变形。随着制动压力的增大而增大,挤压摩擦力和剪切力也增大,进一步加剧了磨损程度。

图10 不同制动压力下闸片磨损曲线

在10 s时间内,制动压力为17 kN时最大磨损深度为1.52×10-2mm,随后压力每增加3 kN,最大磨损深度分别增加2.313 8×10-3、2.580 6×10-3、2.231 4×10-3、2.332 9×10-3、2.253 7×10-3mm,利用差值计算磨损变化率绘图,见图11。随着制动压力的增大,摩擦力线性增大,摩擦制动使得闸片温升变大,导致闸片发生热变形,在摩擦副相对滑动的同时,由于热变形的影响,摩擦接触面积一直变化,因此磨损差值会呈现出上下波动,但其磨损变化率呈现逐渐降低的趋势。

(a) 磨损深度差值

(b) 磨损变化率图11 磨损变化曲线

3 动车组闸片寿命预测系统设计

3.1 动车组闸片磨损预测模型的建立

数据预处理保证数据库完整性的同时还可以减小数据库的量,数据预处理包含数据去噪、特征选取以及特征数据处理3个过程。

(1)输入输出层结构

以闸片各节点的磨损值为输出参数,输入层为制动初速度、制动压力、制动时间。输出层神经元包括闸片209个节点的磨损仿真数值。

(2)隐藏层结构

本文采用单隐层、3层输入、多输出的BP神经网络建立预测模型。目前,隐层神经元数量有几种相关的经验公式可供参考,本文参考公式(16):

(16)

式中:m、n为输入层、输出层神经元数量;a∈[1,10]。

确定该网络结构的隐藏层神经元数量个数为15～24,该模型选择18个隐层神经元。

(3)激励函数

选择S型正切函数tansig作为隐藏层神经元的激励函数,选取S型对数函数logsig作为输出层神经元的激励函数。

(4)学习率

学习率大小直接影响学习精度、学习速度。通常情况下,学习率的范围大概取0.01～0.70即可,本模型学习率为0.01。

选择3种制动压力,分别为14、23、32 kN;初速度选择范围为60～200 km/h,以20 km/h为一个跨度;仿真时长20 s,以2 s为一个时间跨度选择数据进行训练,共提取出230组数据进行模型的训练、验证以及测试。

表3 为Progress区域训练结果统计表,Epoch表示实际训练次数(19次);Time表示训练总时间(31 min 9 s);Performance为性能指标,该模型训练均方误差为1.37×10-7;Gradient为梯度,进度条中显示当前的梯度值为8.37×10-8,右侧设定的梯度值为1×10-7,梯度值达到设定值则停止训练。

表3 网络训练结果数据统计表

3.2 磨损量-行驶里程对应关系的建立

目前复兴号动车组闸片厚度原型为20 mm,其余动车组闸片厚度原型为17 mm。其中毫米磨耗行驶里程=寿命万公里/(原型厚度-磨耗到限)。根据最大磨损深度计算失效磨耗量,以CRH380BL上的博科瑞BKR380FA为例,动车组闸片最大寿命为79.3 万km(取79 万km),拖车为44.2 万km(取44 万km)。闸片原型厚度为17 mm,允许磨耗到限厚度最薄为5 mm。对上述数据进行计算,闸片每磨耗1 mm,可行驶里程动车取n1=658 33 km,拖车取n2=36 666 km。根据动车与拖车行驶里程和单位磨损量的关系,当动车组闸片磨损x时,计算得到拖车闸片磨损x2=xn1/n2≈1.8x。综合以上选择数据,其寿命计算为:

y1=(12-x)×65 833

(17)

y2=(12-1.8x)×36 666

(18)

式中:y1、y2分别为动车组闸片的剩余行驶里程和拖车闸片的剩余行驶里程;x为动车组闸片磨损厚度。

3.3 系统软件界面设计结果及应用

寿命预测计算App基于Matlab中的图形用户界面设计,在用户与计算机之间建立交互操作。使用面板、图标按钮等功能框进行设计界面布局,通过编辑回调函数,逻辑选择的switch分支结构见图12。

图12 switch分支结构

图13为预测系统使用案例。可以看出,当输入制动初速度为160 km/h、制动压力为32 kN、制动时间为20 s时,磨损预测模型磨损深度预测结果最大值为0.030 76 mm,期望输出最大目标值为0.030 741 mm,误差值为1.9×10-5mm;磨损深度预测结果平均值为0.024 52 mm,期望输出目标的平均值为0.024 485 mm,误差值为3.5×10-5mm,可以看出磨损深度预测精度较高。随后针对寿命计算模块进行型号选择,闸片型号选择DC20,动车型号选择CRH380AL,选择型号之后动车、拖车磨损极限弹出初始值为12 mm,此时计算按钮的回调函数根据switch语句选择出CRH380AL车型上DC20闸片对应的寿命计算方式,随后根据输入磨损量计算剩余行驶里程以及剩余磨损厚度,并自动更新迭代至动车、拖车磨损极限。

图13 预测系统使用案例

4 结论

(1) 本文基于热-力耦合仿真模型,结合自适应网格技术进行二次开发,模拟计算了闸片磨损过程,分析了制动初始参数对于闸片磨损的影响。结果表明:该模型磨损仿真结果与试验闸片磨损分布趋势相同,闸片的摩擦出口区域出现了较严重的磨损,磨损在径向呈现梯度分布。制动初速度从80 km/h增大到180 km/h时,最大磨损深度增加了1.40×10-2mm。制动压力从17 kN增大到32 kN时,最大磨损深度增加了0.94×10-2mm。

(2)本文主要分析动车组闸片寿命预测的实现途径及实现过程,根据闸片磨损仿真结果,建立磨损与制动参数对应关系的数据库,基于BP神经网络模型,运用数据库数据对预测模型进行了训练、验证和测试,建立了闸片磨损预测模型。根据动车组运行数据建立了闸片磨损量-行驶里程的对应关系。通过编程语言调用训练好的模型,编写了各车型闸片的寿命计算逻辑,建立了集磨损预测、寿命计算为一体的动车组闸片寿命预测系统。