盘式制动器热-应力-磨损耦合行为的数值模拟*

张方宇，桂良进，范子杰

(清华大学，汽车安全与节能国家重点实验室，北京 100084)

前言

盘式制动器的制动性能主要取决于摩擦片和制动盘之间的摩擦状况。制动过程中，摩擦力做功生热，摩擦片被磨损，热、接触应力和磨损三者相互耦合，共同影响制动器的制动性能。

涉及磨损的理论有很多，如文献[1]中的线性磨损理论、文献[2]中提出的磨粒磨损微切削理论、文献[3]中提出的粘着理论、文献[4]中提出的疲劳磨损理论、文献[5]中提出的能量磨损理论、文献[6]中提出的剥层磨损理论、文献[7]中基于“摩擦疲劳”概念发展的相应的理论和文献[8]中提出的磨损率与磨损循环次数的幂指数关系法则。其中文献[1]中的线性磨损理论应用最为广泛。

目前对于制动热应力和磨损行为的研究主要采用的是有限单元法。文献[9]中研究了飞机正常着陆时，制动盘在制动过程中的热应力分布，以及热应力对C/C复合材料磨损表面形貌的影响，并利用有限元软件进行了热应力的仿真，但未对磨损行为进行有限元的模拟。文献[10]中利用有限元仿真软件Abaqus对织物衬垫自润滑关节轴承摩擦副间的应力分布进行了仿真，并根据Archard的磨损理论，计算了单步磨损量，以此循环，利用显式积分的方法计算总磨损量并描绘了应力分布情况的演变。文献[11]中使用了类似的方法对盘式制动器磨损行为进行了研究，其中还使用了模型结构矩阵(MSM)[12]的模块化建模方法，但它们在计算磨损时均未考虑摩擦热的影响。总之，当前对盘式制动器进行热-应力-磨损耦合行为的数值模拟，尚无文献报到。

本文中探索了盘式制动器热-应力-磨损耦合行为的数值模拟方法。首先建立了盘式制动器的有限元模型；然后利用Archard的磨损理论实现了摩擦片接触应力-磨损耦合分析的数值计算方法，并模拟了摩擦片和制动盘的磨合过程；最后探索了模型在热-应力耦合条件下摩擦片的磨损行为，并对汽车100km时速紧急制动工况下的摩擦片磨损过程进行了仿真。

1 有限元模型

盘式制动器组件的结构示意图如图1所示。本文中只关心摩擦片和制动盘接触界面的接触应力和磨损的分布以及摩擦热和热传导的过程，因此只对活塞、背板、摩擦片和制动盘进行有限元建模。制动盘两侧的制动衬片在几何形状、支撑形式和受力上虽有差异，但在与制动盘的作用形式方面还是比较接近，可近似地认为它们关于制动盘中面对称，因此可以只对一个制动衬片和半个制动盘建模，简化后的结构如图2所示。另外，由于活塞和制动盘材料较硬，制动过程中的变形对仿真结果影响很小，其磨损率相比摩擦片小很多，可以忽略，因此两者均可以当作刚体处理，只将制动衬片(包括背板和摩擦片)视为变形体来建模，且忽略制动盘的磨损，只考虑摩擦片的磨损。

整个有限元模型的爆炸图如图3所示。其中，活塞是解析刚体；背板是三维变形体，共包含1 467个节点，876个单元，单元类型为C3D8T；摩擦片是三维变形体，包含2 420个节点，1 662个单元，单元类型为C3D8T；制动盘是离散刚体，包含1 500个节点，960个单元，单元类型为C3D8T。背板的材料是普通钢铁，使用线弹性材料特性，弹性模量E=2.09×109Pa，泊松比μ=0.3，其热力学材料性能参数来自通用钢材的材料特性数据，其材料性能见表1。

热导率/(W/(m·K))密度/(kg/m3)弹性模量/GPa泊松比膨胀系数/K-1比热容/(J/(kg·K))4878002090311×10-5452

摩擦片的材料特性极其复杂，此处假定其力学性能为线弹性，弹性模量和泊松比参考文献[11]，热力学材料性能参数参考文献[13]。本文中边界条件、载荷和接触统称为接口(interface)。制动盘只有绕其轴线的转动自由度(i2)，其他各个方向的自由度全部为零(i1)。背板受到制动钳的约束，只能沿z轴方向移动，沿x和y轴被约束(i3)。摩擦片嵌合在背板上(i7)。活塞只能沿z轴移动(i4)，其参考点受到沿z轴负方向的集中载荷作用(i5)，提供制动力。摩擦片和制动盘之间为有摩擦接触约束(i6)并且包含摩擦生热特性，生热效率为100%，分配系数为0.5，接触时相互之间的热传导系数为1×109W/(m2·℃)。背板和活塞之间为无摩擦接触约束(i8)。

在Abaqus中，模型的散热边界条件通过背板、摩擦片和制动盘三者表面的surface film condition实现，其散热方程为

q=-h(θ-θ0)

(1)

式中：q为边界热流密度；h为膜层系数；θ为表面节点温度；θ0为冷源温度。

计算时，设定θ0=20℃，制动盘和背板的膜层系数设为100，摩擦片的膜层系数设为10。

2 应力-磨损耦合分析方法

在制动过程中，摩擦片和制动盘的接触伴随着摩擦和磨损，其程度受众多因素的影响，比如接触应力、滑动速率、材料特性等等。有很多磨损理论揭示了这种关系，其中Archard的摩擦定律相对简单，使用较为广泛，并且工程应用比较成熟，因此本文中也采用此理论模型。

Archard的磨损模型通常具有以下形式：

(2)

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

在无穷小的接触面积ΔA，无穷小时间增量dt下，将此形式的磨损模型推广到微分形式：

(3)

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

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

(4)

式中：FN/ΔA为接触微面积ΔA处的接触应力，假设其为pc；k/H为有量纲磨损系数，可以用kd表示，此时式(4)变为

dh=kdpcds

(5)

如果取有限小的时间增量为Δt，相应的磨损厚度增量为Δh，相对滑动位移增量为Δs，kd和pc假设在Δt内变化不大，仍取瞬时值，则式(5)改写为

Δh=kdpcΔs

(6)

在盘式制动器制动过程中，可以认为摩擦片是静止不动的，Δs为制动盘相对摩擦片滑动位移，Δh为Δt时间内的磨损深度。

不考虑热效应时，kd可认为是常数，Δh是pc的函数，但Δh在摩擦片的各个部位并不相同，从而对接触应力pc又产生影响，所以磨损过程是Δh和pc耦合的过程。为了求解这个耦合过程，本文中采用显式欧拉积分的方法。首先利用有限元软件计算得到接触应力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)

(7)

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

(8)

利用Abaqus Standard的ALE功能来实现磨损的仿真，其中磨损的数值算法通过用户子例程UMESHMOTION实现。ALE算法中设置磨损边界的约束类型为速度约束，将式(7)左右两侧同时除以第j个时间步长Δt(j)，得到如下关系：

(9)

(10)

3 摩擦片-制动盘界面磨合仿真

在全新的制动器进入稳定磨损之前，摩擦片-制动盘接触界面先要经历磨合过程。由于整个磨合期滑动位移和接触时间很长，现有的普通计算机的计算能力几乎无法模拟整个过程中温度场、接触应力场和磨损三者的耦合行为，因此在摩擦片-制动盘接触界面磨合过程的仿真中忽略了温度场的影响，只对接触应力-磨损两个场的耦合过程进行模拟。

磨合期摩擦片和制动盘之间的摩擦磨损可以看作一个准静态过程，其工况如图4(a)所示，对制动器的活塞施加沿制动盘轴线方向1.8kN载荷(相当于1MPa制动液压力)，制动盘做匀速运动。

此工况的仿真在Abaqus中采用Static General 算法来实现，计算用时约3h，计算得到的摩擦片接触应力分布和累积磨损量分布随制动盘转数的演变过程如图5和图6所示。

由图5和图6可见：

(1)制动盘刚开始转动(0转)时，摩擦片-制动盘接触界面在受活塞底部环面挤压的部位接触应力比较大，其他部位接触应力相对较小。但随制动盘转数的增加，应力大的区域磨损较大，使接触应力的分布又趋于均匀；

(2)制动盘刚开始转动(0转)时，摩擦片-制动盘接触界面接触应力的分布左右并不对称，这是由摩擦力的存在造成的。制动盘从摩擦片的左侧向右侧滑过，它们之间的摩擦力使摩擦片的左侧压紧制动盘，右侧脱离接触，因此右侧的接触应力为零；但是随着制动盘转数的增加，接触面中间应力较大，磨损较快，而右侧没有接触，也就没有磨损，这就导致摩擦片右侧逐渐进入接触，接触应力由小变大；

(3)当制动盘转数达到一定值后(10 000～20 000转)，摩擦片接触面上的接触应力分布达到稳定状态，几乎不再变化，其应力值沿制动盘半径方向逐渐减小，即从摩擦片内侧向外侧逐渐减小。这是由于在Archard磨损模型里磨损量与滑动距离和接触应力的乘积成正比，当磨损达到稳定时，摩擦片的整个接触面均匀地向背板方向磨损，也就是说接触面上各点的磨损速率相等，因此接触应力与滑动距离成反比。又因为滑动距离和制动盘上相应点的半径成正比，所以接触应力和半径成反比，这也就是接触面应力沿半径方向均匀减小的原因。

以上3点分析在一定程度上可以验证利用Abaqus中ALE算法模拟磨损过程的正确性。

4 摩擦片-制动盘界面紧急制动仿真

在制动盘转动20 000转之后，摩擦片和制动盘结束了磨合期，进入稳定磨损阶段。本节在上一节摩擦片和制动盘完成磨合期(20 000转)后的摩擦片轮廓基础上，研究一次紧急制动过程中热-应力-磨损三者的耦合行为。一次紧急制动工况如图4(b)所示，对制动盘的活塞施加沿制动盘轴线方向9kN载荷(相当于5MPa制动液压力)，制动盘做初速度为120rad/s(相当于汽车100km/h)的匀减速转动，3s后停止。

本文中利用Abaqus中Coupled temp-displacement (Transient) 算法来实现此工况，同时模拟了考虑和不考虑磨损两种情况，并对结果做出了比较。

考虑磨损时，计算用时约16h，计算得到的摩擦片接触应力分布、温度分布、累积磨损量分布和制动盘温度分布随时间的演变过程如图7～图10所示。

对图7～图10展示的结果做如下分析。

(1)一次紧急制动的制动力较大，所以开始制动 (t=0.001s) 时，摩擦片和制动盘接触界面在受活塞环挤压的部位接触应力较大，但随着摩擦力做功，接触界面温度升高，摩擦片向制动盘方向膨胀，这时摩擦片表面的中心部位与制动盘接触，接触应力的分布变为中心大，四周小。

(2)接触应力大的部位摩擦功率大，温度升高快，因此摩擦片温度分布和接触应力的分布大体一致。但是考虑到摩擦片与制动盘和空气的热交换，以及摩擦片内部的热传导都会使摩擦片的温度分布趋于均匀，因此摩擦片的温度分布又与应力分布不完全一致。

(3)刚开始制动时，制动盘转速快，摩擦功率大，摩擦片的温度迅速上升。但在制动过程中随着制动盘转速逐渐减小，摩擦功率也会减小，而随着温度的上升，散热功率增加，因此摩擦片的温升越来越慢，到1.5s左右摩擦片和制动盘的温度达到最高值，随后摩擦功率小于散热功率，摩擦片和制动盘的温度开始下降并趋于均匀。

以上3点分析在一定程度上说明本次模拟比较符合实际经验，具有准确性。

不考虑摩擦片磨损时，计算用时约12h，计算得到的摩擦片接触应力分布、温度分布和制动盘温度分布随时间的演变过程与考虑磨损的情况类似。

选取摩擦片接触应力最大值附近一点，节点编号为N2588，如图11(a)所示。在考虑磨损和不考虑磨损两种情况下将该点接触应力随时间变化的曲线进行比较，如图12(a)所示。由图可见，考虑磨损时接触应力值比不考虑磨损时的接触应力值小，其原因是接触应力大的地方磨损量大，使接触应力的分布相对均匀。

选取摩擦片温度最高值附近一点，节点编号为N2602，如图11(b)所示。在考虑磨损和不考虑磨损两种情况下将该点温度随时间变化的曲线进行比较，如图12(b)所示。从图中可以看到，考虑磨损时温度值比不考虑磨损时的温度值略小，其原因是考虑磨损时，接触应力相对均匀，摩擦生热的热流密度也相对均匀。但这种差别不大，原因是一次紧急制动磨损量很小，对温度分布的影响不显著。

5 结论

本文中探索了盘式制动器热-接触应力-磨损耦合行为的数值模拟方法。首先建立了盘式制动器的有限元模型；然后模拟了盘式制动器磨合期摩擦片的接触应力-磨损耦合过程；最后模拟了模型在一次紧急制动过程中热-接触应力-磨损三者的耦合情况，并将结果与不考虑磨损情况下的结果进行了比较。分析结果表明，所提出的模拟方法是有效的，可为盘式制动器的结构参数设计提供理论指导。

[1] Archard J F. Contact and Rubbing of Flat Surfaces [J]. Journal of Applied Physics，1953, 24: 981-988.

[2] 赫罗绍夫M M, 巴比契夫M A. 金属的磨损 [M]. 胡绍衣, 佘沪生，译.北京: 机械工业出版社, 1966.

[3] Bowden F P,Tabor D.The Friction and Lubrication of Solids [M].Oxford at the Clarenden Press, 1964.

[4] 克拉盖尔斯基иB，陀贝钦M H, 康巴洛夫B C. 摩擦磨损计算原理 [M]. 汪一麟,朱安仁，范明德，译.北京: 机械工业出版社, 1982.

[5] Fleisher G. Energishe Methode der Bestimung der Verschleiβes [J]. Schmierungstechnik, 1973(4): 9.

[6] Suh N P. An Overview of the Delamination Theory of Wear [J]. Wear, 1977, 44(1):1-16.

[7] СосновскийЛ А.Основы Трибофатики [M]. Гомель, 2003.

[8] Siniawski M T, Harris S J, Wang Q. A Universal Wear Law for Abrasion [J]. Wear, 2007, 262: 883-888.

[9] 徐惠娟, 易茂中, 熊翔，等. 制动过程中热应力对C/C复合材料磨损表面形貌的影响 [J]. 中国有色金属学报, 2011, 21(1): 131-137.

[10] Shen X,Cao L, Li R. Numerical Simulation of Sliding Wear Based on Archard Model[C]. Mechanic Automation and Control Engineering (MACE)，2010：325-329.

[11] Soderberg A, Andersson S. Simulation of Wear and Contact Pressure Distribution at the Pad-to-rotor Interface in a Disc Brake Using General Purpose Finite Element Analysis Software [J]. Wear, 2009, 267: 2243-2251.

[12] Sellgren U. Architecting Models of Technical Systems for Non-routine Simulations[C]. Proceedings of ICED 03-14th International Conference on Engineering Design, Stockholm, 19-21 August 2003.

[13] Day A J, Newcomb T P. Dissipation of Frictional Energy from the Interface of an Annular Disc Brake [C]// Proceedings of the Institution of Mechanical Engineers, Part D: Transport Engineering, 1984, 198(11): 201-209.