结构热应力仿真-主题088_微动疲劳与磨损-主题088_微动疲劳与磨损_完整教程
主题088：微动疲劳与磨损（Fretting Fatigue and Wear）
摘要
微动疲劳与磨损是机械工程领域中一类重要的表面失效形式，广泛存在于螺栓连接、花键配合、叶片-盘连接、轴承、齿轮等关键工程结构中。本教程系统介绍微动疲劳的物理机理、数学模型和数值仿真方法，包括Hertz接触理论、Mindlin局部滑移理论、微动图（Fretting Map）方法、Archard磨损模型以及多轴疲劳寿命预测理论。通过两个完整的Python仿真实例，读者将掌握微动疲劳寿命预测、表面形貌演化模拟和磨损机制分析的核心技术。本教程适用于航空发动机、汽车、轨道交通等领域的结构耐久性设计工程师和研究人员。
关键词
微动疲劳；磨损；接触力学；Mindlin理论；Archard模型；表面形貌演化；寿命预测；多轴疲劳
目录
- 引言
- 微动疲劳基础理论
- 接触力学与Hertz理论
- Mindlin局部滑移理论
- 微动图与运行状态
- 磨损理论与模型
- 多轴疲劳与寿命预测
- 实例一：微动疲劳机理与寿命预测
- 实例二：磨损模型与表面演化模拟
- 工程应用案例
- 防护措施与设计建议
- 常见报错与解决方案
- 进阶挑战
- 总结与展望
1. 引言
1.1 微动疲劳的工程背景
微动疲劳（Fretting Fatigue）是指两个接触表面在循环载荷作用下发生微小相对运动（通常在微米量级，1-100μm）时，接触区域产生的疲劳损伤现象。这种微小的相对运动称为"微动"（Fretting），它虽然幅度很小，但由于发生在高应力集中区域，能够显著加速裂纹的萌生与扩展，导致结构提前失效。
微动疲劳在工程实践中极为常见，典型应用场景包括：
航空发动机领域：涡轮叶片与轮盘之间的榫连接（Dovetail or Fir-tree Joint）是微动疲劳的高发区域。发动机工作时，叶片受到离心力和气动力的交变作用，在叶片根部与轮盘的接触面上产生微动，导致裂纹萌生。据统计，航空发动机中约30%的结构失效与微动疲劳有关。
汽车工程领域：发动机连杆螺栓、底盘悬挂系统的螺栓连接、传动系统中的花键配合等都存在微动疲劳问题。特别是在振动环境下，这些连接部位的微动会显著降低结构的疲劳寿命。
轨道交通领域：铁路轨道的扣件系统、车轮与车轴的压配合、受电弓与接触网的滑动接触等都涉及微动磨损和疲劳问题。
能源装备领域：核电站蒸汽发生器传热管与支撑板的接触区域、风力发电机叶片螺栓连接、海上平台结构螺栓等，在交变载荷和环境因素共同作用下，微动疲劳是主要的失效模式之一。
医疗器械领域：人工关节（如髋关节、膝关节）的球头与臼杯之间的微动磨损是导致假体松动和失效的重要原因。牙科植入物与骨组织之间的微动也会影响骨整合效果。
1.2 磨损的工程意义
磨损（Wear）是指两个接触表面在相对运动过程中，由于机械作用、化学作用或热作用导致的材料逐渐损失现象。磨损是机械零件失效的三种主要形式之一（另外两种是断裂和腐蚀），据统计，约有80%的机械零件失效与磨损有关。
磨损不仅造成巨大的经济损失（据估计，磨损造成的经济损失约占工业国家GDP的1-2%），还会导致设备精度下降、能耗增加、噪声增大，甚至引发严重的安全事故。因此，深入研究磨损机理、建立准确的磨损预测模型、开发有效的防护措施，对于提高机械设备的可靠性和使用寿命具有重要意义。
1.3 微动疲劳与磨损的耦合效应
在实际工程问题中，微动疲劳与磨损往往同时发生、相互影响。微动过程中的磨损会改变接触表面的形貌和粗糙度，进而影响接触应力分布和微动特性；同时，微动引起的表面损伤和材料硬化又会改变材料的磨损行为。
这种耦合效应使得微动疲劳与磨损问题变得极为复杂。在某些情况下，适当的磨损可以去除表面的微裂纹和缺陷，反而有利于延长疲劳寿命；而在另一些情况下，严重的磨损会导致接触面积减小、应力集中加剧，加速疲劳失效。因此，准确评估微动疲劳与磨损的耦合效应，是结构耐久性设计的关键挑战之一。
1.4 本教程的学习目标
通过本教程的学习，读者将能够：
-
理解微动疲劳的物理机理：掌握微动疲劳的萌生机制、裂纹扩展规律和影响因素。
-
掌握接触力学基础理论：理解Hertz接触理论、Mindlin局部滑移理论及其在微动分析中的应用。
-
熟悉磨损理论与模型：掌握Archard磨损模型、粘着磨损、磨粒磨损、冲蚀磨损等主要磨损类型的数学描述。
-
建立数值仿真能力：能够使用Python编写微动疲劳寿命预测和表面形貌演化模拟程序。
-
应用于工程实践：具备分析实际工程结构中微动疲劳与磨损问题的能力，提出有效的设计改进和防护措施。
2. 微动疲劳基础理论
2.1 微动疲劳的特征与机理
2.1.1 微动的定义与分类
微动（Fretting）是指两个接触表面之间发生的振幅极小的相对往复运动。根据运动形式，微动可以分为以下几类：
切向微动（Tangential Fretting）：接触表面沿切向方向的相对往复运动，是最常见的微动形式。在切向微动中，接触区域同时承受法向压力和切向摩擦力，产生复杂的应力状态。
径向微动（Radial Fretting）：接触表面沿法向方向的相对往复运动，通常由交变法向载荷引起。径向微动在滚动轴承、齿轮等零部件中较为常见。
扭动微动（Torsional Fretting）：接触表面绕法向轴的相对往复转动，常见于球轴承、螺栓连接等结构中。
复合微动（Compound Fretting）：同时包含两种或多种运动形式的微动，如同时存在切向运动和扭动。
2.1.2 微动疲劳的损伤机理
微动疲劳的损伤过程通常包括以下几个阶段：
第一阶段：表面氧化与磨屑形成（0-10^4次循环）
在微动初期，由于接触表面的相对运动，表面氧化膜被破坏，新鲜金属表面暴露并发生氧化。同时，微动产生的剪切作用使表面材料发生塑性变形，形成磨屑（Wear Debris）。这些磨屑主要由氧化物（如Fe2O3、Fe3O4等）和金属颗粒组成，它们的产生会改变接触界面的摩擦特性。
第二阶段：表面裂纹萌生（104-105次循环）
随着微动的持续，接触表面在循环载荷作用下逐渐产生微裂纹。裂纹通常萌生于表面或次表面的应力集中区域，如接触边缘、表面缺陷、夹杂物等处。微动引起的多轴应力状态和应变集中是促进裂纹萌生的主要因素。
第三阶段：裂纹扩展（105-106次循环）
萌生的微裂纹在循环载荷作用下逐渐扩展。在微动疲劳中，裂纹通常以与表面成一定角度（约15-30度）向材料内部扩展，然后逐渐转向垂直于主应力方向。裂纹扩展速率受应力强度因子幅值、材料性能、环境因素等多种因素影响。
第四阶段：最终断裂（>10^6次循环）
当裂纹扩展到临界尺寸时，剩余截面无法承受外加载荷，发生快速断裂，导致结构失效。
2.1.3 微动疲劳的影响因素
微动疲劳寿命受多种因素影响，主要包括：
接触压力：接触压力直接影响接触区域的应力水平和滑移状态。过高的接触压力会加速表面损伤和裂纹萌生，但过低的接触压力可能导致完全滑移，反而降低疲劳寿命。
滑移幅值：滑移幅值是微动疲劳的关键参数。研究表明，存在一个临界滑移幅值（通常在10-50μm范围），在此滑移幅值下疲劳寿命最低。滑移幅值过小（<5μm）时，处于部分滑移状态，损伤较小；滑移幅值过大（>100μm）时，磨损成为主导机制，反而可能延缓裂纹萌生。
循环次数与频率：微动疲劳是典型的循环载荷问题，疲劳寿命与循环次数密切相关。频率的影响较为复杂，高频率可能导致温升和氧化加速，低频率则可能使腐蚀因素更加显著。
材料性能：材料的强度、硬度、韧性、疲劳性能等都会影响微动疲劳行为。一般来说，高强度材料对微动疲劳更敏感，因为其裂纹萌生寿命较短，而裂纹扩展寿命相对较长。
环境因素：温度、湿度、腐蚀介质等环境因素对微动疲劳有显著影响。高温会加速氧化和蠕变，潮湿环境会促进腐蚀疲劳，腐蚀性介质则会加速表面损伤和裂纹扩展。
表面处理：表面粗糙度、残余应力、表面强化处理（如喷丸、渗碳、氮化等）都会影响微动疲劳性能。适当的表面处理可以显著提高微动疲劳寿命。
2.2 微动疲劳与传统疲劳的区别
微动疲劳与传统疲劳（无接触效应的疲劳）相比，具有以下显著特点：
多轴应力状态：微动疲劳中，接触区域同时承受法向接触应力和切向摩擦应力，形成复杂的多轴应力状态。这种多轴应力状态会加速裂纹萌生，通常使疲劳寿命降低一个数量级以上。
应力集中效应：接触边缘存在严重的应力集中，这是微动裂纹的主要萌生位置。传统疲劳中的应力集中通常来自几何不连续（如缺口、孔洞），而微动疲劳的应力集中是由接触本身引起的。
表面损伤累积：微动过程中的磨损和氧化会不断改变表面状态，产生磨屑和表面粗糙度变化，这些表面损伤会进一步影响疲劳行为。传统疲劳中，表面状态相对稳定。
裂纹萌生位置：微动疲劳的裂纹通常萌生于表面或近表面（<100μm深度），而传统高周疲劳的裂纹可能萌生于材料内部的缺陷或夹杂物。
寿命预测方法：微动疲劳寿命预测需要考虑接触力学、磨损、多轴疲劳等多种因素，比传统疲劳更为复杂。常用的方法包括应力-寿命法（S-N曲线）、应变-寿命法（ε-N曲线）、断裂力学法等。
3. 接触力学与Hertz理论
3.1 接触问题的基本概念
接触力学是研究两个或多个物体在接触区域相互作用规律的学科。在微动疲劳分析中，准确计算接触区域的应力分布是进行寿命预测的基础。
3.1.1 接触问题的分类
根据接触体的几何形状和变形特性，接触问题可以分为以下几类：
Hertz接触：两个弹性体在无摩擦条件下的接触，接触区域尺寸远小于接触体曲率半径。这是最常见的接触类型，也是Hertz理论的基本假设。
非Hertz接触：接触区域较大、接触体发生大变形、或存在摩擦等复杂情况，需要采用数值方法（如有限元法）求解。
粗糙表面接触：实际工程表面都是粗糙的，粗糙峰之间的接触是离散分布的。粗糙表面接触需要考虑多尺度效应和统计特性。
动态接触：接触体之间存在相对运动或冲击载荷，需要考虑惯性效应和波动传播。
3.1.2 接触问题的基本方程
接触问题的求解需要满足以下基本方程：
平衡方程：接触体内部的应力必须满足平衡方程：
∂ σ i j ∂ x j + f i = 0 \frac{\partial \sigma_{ij}}{\partial x_j} + f_i = 0 ∂xj∂σij+fi=0
其中， σ i j \sigma_{ij} σij是应力张量， f i f_i fi是体积力。
几何方程：应变与位移的关系：
ε i j = 1 2 ( ∂ u i ∂ x j + ∂ u j ∂ x i ) \varepsilon_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) εij=21(∂xj∂ui+∂xi∂uj)
本构方程：应力与应变的关系（线弹性材料）：
σ i j = λ ε k k δ i j + 2 μ ε i j \sigma_{ij} = \lambda \varepsilon_{kk} \delta_{ij} + 2\mu \varepsilon_{ij} σij=λεkkδij+2μεij
其中， λ \lambda λ和 μ \mu μ是Lamé常数，与弹性模量 E E E和泊松比 ν \nu ν的关系为：
μ = E 2 ( 1 + ν ) , λ = E ν ( 1 + ν ) ( 1 − 2 ν ) \mu = \frac{E}{2(1+\nu)}, \quad \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)} μ=2(1+ν)E,λ=(1+ν)(1−2ν)Eν
边界条件：接触边界上需要满足位移协调条件和应力连续条件。
3.2 Hertz接触理论
Hertz接触理论是1882年由Heinrich Hertz提出的，用于描述两个弹性体在法向载荷作用下的接触问题。该理论基于以下基本假设：
- 接触体是线弹性、各向同性、均质的
- 接触区域尺寸远小于接触体曲率半径
- 接触面无摩擦
- 小变形假设成立
3.2.1 圆柱-平面接触
对于半径为 R R R的圆柱与平面接触，在法向载荷 P P P（单位长度载荷，N/m）作用下：
接触半宽：
b = 4 P R π E ∗ b = \sqrt{\frac{4PR}{\pi E^*}} b=πE∗4PR
最大接触压力（接触中心处）：
p 0 = 2 P π b = P E ∗ π R p_0 = \frac{2P}{\pi b} = \sqrt{\frac{PE^*}{\pi R}} p0=πb2P=πRPE∗
接触压力分布（半椭圆分布）：
p ( x ) = p 0 1 − ( x b ) 2 p(x) = p_0 \sqrt{1 - \left(\frac{x}{b}\right)^2} p(x)=p01−(bx)2
其中， E ∗ E^* E∗是等效弹性模量：
1 E ∗ = 1 − ν 1 2 E 1 + 1 − ν 2 2 E 2 \frac{1}{E^*} = \frac{1-\nu_1^2}{E_1} + \frac{1-\nu_2^2}{E_2} E∗1=E11−ν12+E21−ν22
E 1 E_1 E1、 ν 1 \nu_1 ν1和 E 2 E_2 E2、 ν 2 \nu_2 ν2分别是两个接触体的弹性模量和泊松比。
次表面应力场：
在接触中心正下方（ x = 0 x=0 x=0），沿深度方向（ z z z方向）的应力分量为：
σ x = − p 0 [ ( 1 + z 2 b 2 ) − 1 / 2 − 2 z b ] \sigma_x = -p_0 \left[\left(1+\frac{z^2}{b^2}\right)^{-1/2} - 2\frac{z}{b}\right] σx=−p0[(1+b2z2)−1/2−2bz]
σ z = − p 0 ( 1 + z 2 b 2 ) − 1 / 2 \sigma_z = -p_0 \left(1+\frac{z^2}{b^2}\right)^{-1/2} σz=−p0(1+b2z2)−1/2
最大剪应力出现在深度 z ≈ 0.786 b z \approx 0.786b z≈0.786b处，其值为 τ m a x ≈ 0.3 p 0 \tau_{max} \approx 0.3p_0 τmax≈0.3p0。
3.2.2 球-平面接触
对于半径为 R R R的球与平面接触，在法向载荷 F F F作用下：
接触半径：
a = ( 3 F R 4 E ∗ ) 1 / 3 a = \left(\frac{3FR}{4E^*}\right)^{1/3} a=(4E∗3FR)1/3
最大接触压力：
p 0 = 3 F 2 π a 2 = ( 6 F E ∗ 2 π 3 R 2 ) 1 / 3 p_0 = \frac{3F}{2\pi a^2} = \left(\frac{6FE^{*2}}{\pi^3 R^2}\right)^{1/3} p0=2πa23F=(π3R26FE∗2)1/3
接触压力分布（半球形分布）：
p ( r ) = p 0 1 − ( r a ) 2 p(r) = p_0 \sqrt{1 - \left(\frac{r}{a}\right)^2} p(r)=p01−(ar)2
接触接近量（两物体相互靠近的距离）：
δ = a 2 R = ( 9 F 2 16 R E ∗ 2 ) 1 / 3 \delta = \frac{a^2}{R} = \left(\frac{9F^2}{16RE^{*2}}\right)^{1/3} δ=Ra2=(16RE∗29F2)1/3
3.2.3 一般曲面接触
对于两个任意曲面接触，首先需要确定主曲率半径。设接触点处两个表面的主曲率半径分别为 R 11 R_{11} R11、 R 12 R_{12} R12和 R 21 R_{21} R21、 R 22 R_{22} R22，则等效曲率半径为：
1 R ∗ = 1 R 11 + 1 R 12 + 1 R 21 + 1 R 22 \frac{1}{R^*} = \frac{1}{R_{11}} + \frac{1}{R_{12}} + \frac{1}{R_{21}} + \frac{1}{R_{22}} R∗1=R111+R121+R211+R221
接触区域为椭圆形，长半轴 a a a和短半轴 b b b可以通过Hertz接触方程求解。
3.3 Hertz理论的局限性
Hertz理论虽然是接触力学的基石，但在实际应用中存在以下局限性：
摩擦效应：Hertz理论假设接触面无摩擦，但实际微动问题中摩擦起着关键作用。摩擦会改变接触区域的应力分布，导致剪应力集中。
大变形：当接触压力很高或材料较软时，小变形假设不再成立，需要考虑几何非线性和材料非线性。
塑性变形：当接触应力超过材料的屈服强度时，接触区域会发生塑性变形，Hertz弹性解不再适用。
粗糙表面：实际表面都是粗糙的，粗糙峰之间的接触是离散的，与Hertz理论的连续接触假设不同。
动态效应：Hertz理论是静态理论，对于冲击、振动等动态接触问题需要考虑惯性效应。
尽管存在这些局限性，Hertz理论仍然是接触力学分析的重要基础，许多高级接触模型都是在Hertz理论的基础上发展起来的。
4. Mindlin局部滑移理论
4.1 局部滑移现象
当两个接触表面在法向载荷作用下，再施加切向力时，接触区域的响应取决于切向力的大小：
完全粘着（Full Stick）：当切向力很小时，整个接触区域保持粘着状态，无相对滑移。此时切向位移与切向力呈线性关系。
完全滑移（Gross Slip）：当切向力达到极限摩擦力 μ P \mu P μP时，整个接触区域发生滑移，两表面之间产生宏观相对运动。
局部滑移（Partial Slip）：当切向力介于上述两者之间时，接触区域分为两部分：中心区域保持粘着，边缘区域发生局部滑移。这是微动疲劳中最常见的情况。
4.2 Mindlin理论的基本假设
1949年，R.D. Mindlin在Hertz接触理论的基础上，提出了考虑切向力的接触理论。Mindlin理论的基本假设包括：
- 接触体是线弹性、各向同性的
- 接触区域形状与Hertz理论相同（切向力不改变接触区域形状）
- 库仑摩擦定律适用于滑移区域
- 粘着区域内剪应力小于极限摩擦应力
- 小变形假设成立
4.3 Mindlin理论的核心公式
4.3.1 切向力与滑移区关系
对于圆柱接触，当施加切向力 Q Q Q时，滑移区半宽 c c c与接触半宽 b b b的关系为：
c b = 1 − Q μ P \frac{c}{b} = \sqrt{1 - \frac{Q}{\mu P}} bc=1−μPQ
其中， μ \mu μ是摩擦系数， P P P是法向载荷。
粘着区半宽： c = b 1 − η c = b\sqrt{1 - \eta} c=b1−η
滑移区范围： c < ∣ x ∣ < b c < |x| < b c<∣x∣<b
其中， η = Q / ( μ P ) \eta = Q/(\mu P) η=Q/(μP)是切向力系数， 0 ≤ η ≤ 1 0 \leq \eta \leq 1 0≤η≤1。
4.3.2 切向牵引力分布
Mindlin理论给出了切向牵引力（剪应力） q ( x ) q(x) q(x)的分布：
在粘着区（ ∣ x ∣ ≤ c |x| \leq c ∣x∣≤c）：
q ( x ) = μ p 0 [ 1 − x 2 b 2 − c b 1 − x 2 c 2 ] q(x) = \mu p_0 \left[\sqrt{1-\frac{x^2}{b^2}} - \frac{c}{b}\sqrt{1-\frac{x^2}{c^2}}\right] q(x)=μp0[1−b2x2−bc1−c2x2]
在滑移区（ c < ∣ x ∣ ≤ b c < |x| \leq b c<∣x∣≤b）：
q ( x ) = μ p ( x ) = μ p 0 1 − x 2 b 2 q(x) = \mu p(x) = \mu p_0 \sqrt{1-\frac{x^2}{b^2}} q(x)=μp(x)=μp01−b2x2
4.3.3 切向位移-力关系
在部分滑移状态下，切向相对位移 δ \delta δ与切向力 Q Q Q的关系为：
δ = μ P π E ∗ [ 2 ln ( 2 b c ) − c 2 b 2 ] \delta = \frac{\mu P}{\pi E^*} \left[2\ln\left(\frac{2b}{c}\right) - \frac{c^2}{b^2}\right] δ=πE∗μP[2ln(c2b)−b2c2]
或表示为：
δ = μ P π E ∗ [ ( 2 − η ) ln ( 2 − η η ) − 2 ln 2 ] \delta = \frac{\mu P}{\pi E^*} \left[(2-\eta)\ln\left(\frac{2-\eta}{\eta}\right) - 2\ln2\right] δ=πE∗μP[(2−η)ln(η2−η)−2ln2]
这个关系是非线性的，反映了微动过程中接触刚度随切向力变化的特性。
4.4 微动迟滞回线（Fretting Loop）
在循环切向载荷作用下，切向力 Q Q Q与相对位移 δ \delta δ形成闭合的迟滞回线，称为微动迟滞回线（Fretting Loop）。
4.4.1 迟滞回线的特征
部分滑移状态：迟滞回线呈细长的"叶片"形状，回线面积较小。此时能量耗散主要来自接触边缘的局部滑移。
完全滑移状态：迟滞回线接近矩形，回线面积较大。此时能量耗散来自整个接触面的宏观滑移。
混合状态：迟滞回线形状介于上述两者之间。
4.4.2 能量耗散计算
每个循环的能量耗散 W W W等于迟滞回线的面积：
W = ∮ Q d δ W = \oint Q \, d\delta W=∮Qdδ
对于简化的部分滑移模型，能量耗散可以近似为：
W ≈ 4 μ P δ m a x 3 ( 1 − δ m a x δ s l i p ) 3 W \approx \frac{4\mu P \delta_{max}}{3} \left(1 - \frac{\delta_{max}}{\delta_{slip}}\right)^3 W≈34μPδmax(1−δslipδmax)3
其中， δ m a x \delta_{max} δmax是最大相对位移幅值， δ s l i p \delta_{slip} δslip是完全滑移临界位移。
能量耗散是微动损伤的重要指标，与磨损量和疲劳损伤密切相关。
4.5 Mindlin理论的应用与扩展
4.5.1 球接触Mindlin理论
对于球-平面接触，Mindlin理论同样适用，但公式形式略有不同：
滑移区半径 c c c与接触半径 a a a的关系：
c a = ( 1 − Q μ F ) 1 / 3 \frac{c}{a} = \left(1 - \frac{Q}{\mu F}\right)^{1/3} ac=(1−μFQ)1/3
切向位移-力关系：
δ = 3 μ F 8 a G ∗ [ 1 − ( 1 − Q μ F ) 2 / 3 ] \delta = \frac{3\mu F}{8aG^*} \left[1 - \left(1 - \frac{Q}{\mu F}\right)^{2/3}\right] δ=8aG∗3μF[1−(1−μFQ)2/3]
其中， G ∗ G^* G∗是等效剪切模量：
1 G ∗ = 2 − ν 1 G 1 + 2 − ν 2 G 2 \frac{1}{G^*} = \frac{2-\nu_1}{G_1} + \frac{2-\nu_2}{G_2} G∗1=G12−ν1+G22−ν2
4.5.2 扭转Mindlin理论
对于接触面绕法向轴的相对转动，Mindlin也给出了相应的解。扭转力矩 M M M与扭转角 ϕ \phi ϕ的关系为：
ϕ = 3 μ F 8 a 3 G ∗ [ 1 − ( 1 − M μ F a ) 1 / 3 ] \phi = \frac{3\mu F}{8a^3G^*} \left[1 - \left(1 - \frac{M}{\mu F a}\right)^{1/3}\right] ϕ=8a3G∗3μF[1−(1−μFaM)1/3]
4.5.3 Mindlin理论的局限性
理想化假设：Mindlin理论假设接触区域形状不变、材料完全弹性、摩擦遵循库仑定律，这些假设在实际中往往不完全满足。
磨损效应：Mindlin理论没有考虑磨损对接触状态的影响。实际微动过程中，磨损会不断改变接触形貌和应力分布。
多轴效应：Mindlin理论主要处理一维切向问题，对于多轴微动（同时存在切向和扭转）需要更复杂的分析。
尽管存在这些局限，Mindlin理论仍然是理解和分析微动问题的重要工具，为微动图方法和寿命预测模型提供了理论基础。
5. 微动图与运行状态
5.1 微动图的概念
微动图（Fretting Map）是Vingsbo和Soderberg于1988年提出的一种描述微动运行状态的方法。它以相对位移幅值 δ \delta δ和切向力幅值 Q Q Q为坐标轴，将微动运行状态划分为不同的区域。
微动图的基本思想是：微动损伤的严重程度取决于运行状态（部分滑移、混合滑移、完全滑移），而运行状态由位移和力的相对大小决定。
5.2 微动运行状态的分类
根据Mindlin理论和大量实验观察，微动运行状态可以分为三类：
5.2.1 部分滑移区（Partial Slip Regime）
特征：接触区域中心保持粘着，边缘发生局部滑移。迟滞回线呈细长叶片状，能量耗散较小。
判据： Q < Q c r i t i c a l Q < Q_{critical} Q<Qcritical或 δ < δ c r i t i c a l \delta < \delta_{critical} δ<δcritical
其中，临界切向力 Q c r i t i c a l = μ P ( 1 − ( c / b ) 2 ) Q_{critical} = \mu P(1 - (c/b)^2) Qcritical=μP(1−(c/b)2)，临界位移 δ c r i t i c a l \delta_{critical} δcritical是完全滑移开始的位移。
损伤特征：部分滑移区的损伤以疲劳为主，裂纹萌生于接触边缘，向材料内部扩展。由于滑移幅度小，磨损相对较轻。
寿命特征：部分滑移区的疲劳寿命通常最短，因为接触边缘的应力集中和循环滑移共同作用，加速裂纹萌生。
5.2.2 完全滑移区（Gross Slip Regime）
特征：整个接触区域发生宏观滑移，两表面之间有明显的相对运动。迟滞回线接近矩形，能量耗散较大。
判据： Q ≈ μ P Q \approx \mu P Q≈μP或 δ ≫ δ c r i t i c a l \delta \gg \delta_{critical} δ≫δcritical
损伤特征：完全滑移区的损伤以磨损为主，材料通过磨屑形式不断损失。由于滑移幅度大，接触表面不断更新，裂纹难以稳定萌生和扩展。
寿命特征：完全滑移区的疲劳寿命通常比部分滑移区长，因为磨损去除了表面损伤层，且应力集中效应减弱。但严重磨损会导致配合间隙增大、振动加剧。
5.2.3 混合滑移区（Mixed Slip Regime）
特征：介于部分滑移和完全滑移之间，接触区域大部分发生滑移，但仍有小部分保持粘着。迟滞回线形状介于叶片状和矩形之间。
判据： Q c r i t i c a l < Q < μ P Q_{critical} < Q < \mu P Qcritical<Q<μP或 δ c r i t i c a l < δ < δ g r o s s \delta_{critical} < \delta < \delta_{gross} δcritical<δ<δgross
损伤特征：混合滑移区的损伤是疲劳和磨损的复合，既有裂纹萌生扩展，又有明显的材料损失。
寿命特征：混合滑移区的寿命介于部分滑移区和完全滑移区之间，具体取决于疲劳和磨损的相对贡献。
5.3 微动图的构建方法
构建微动图需要确定不同运行状态之间的边界线：
部分滑移-混合滑移边界：对应于接触边缘刚开始发生滑移的状态，可以用Mindlin理论计算：
δ P S − M S = μ P π E ∗ [ 2 ln ( 2 b c ) − c 2 b 2 ] \delta_{PS-MS} = \frac{\mu P}{\pi E^*} \left[2\ln\left(\frac{2b}{c}\right) - \frac{c^2}{b^2}\right] δPS−MS=πE∗μP[2ln(c2b)−b2c2]
混合滑移-完全滑移边界：对应于整个接触面都发生滑移的状态：
δ M S − G S ≈ μ P 2 E ∗ \delta_{MS-GS} \approx \frac{\mu P}{2E^*} δMS−GS≈2E∗μP
在实际应用中，这些边界线通常通过实验确定，因为理论计算的临界值往往与实际情况有偏差。
5.4 微动图在工程设计中的应用
微动图为微动疲劳的防护设计提供了重要指导：
避免部分滑移区：如果设计要求避免微动疲劳失效，应尽量使运行状态避开部分滑移区，可以通过增加预紧力、改善润滑、使用低摩擦涂层等方法实现。
利用完全滑移区：在某些情况下，允许轻微的完全滑移反而有利于延长寿命，因为磨损去除了损伤层。但这需要权衡磨损带来的其他问题（如配合松动、磨屑污染等）。
监测与诊断：通过测量微动过程中的力和位移，绘制实时微动图，可以判断当前的运行状态和损伤程度，为维修决策提供依据。
6. 磨损理论与模型
6.1 磨损的基本类型
磨损根据机理不同，可以分为以下几种主要类型：
6.1.1 粘着磨损（Adhesive Wear）
机理：当两个表面在法向载荷作用下接触时，真实接触只发生在粗糙峰的顶端。在高接触压力下，这些接触点发生冷焊（粘着），形成接合点。当表面相对滑动时，接合点被剪断，材料从一个表面转移到另一个表面，或形成磨屑脱落。
影响因素：
- 材料亲和力：相似材料容易粘着，异种材料组合可以降低粘着
- 接触压力：压力越高，真实接触面积越大，粘着越严重
- 滑动速度：速度影响温升和氧化膜形成
- 环境因素：真空或惰性气体中粘着更严重，氧气可以形成保护氧化膜
典型应用：滑动轴承、活塞环-缸套、金属密封等。
6.1.2 磨粒磨损（Abrasive Wear）
机理：硬颗粒或硬凸起在软材料表面滑动或滚动，通过切削、犁沟作用去除材料。磨粒磨损是最常见的磨损形式之一，约占工业磨损的50%。
分类：
- 两体磨粒磨损：磨粒固定在一个表面上（如锉刀、砂纸）
- 三体磨粒磨损：磨粒在两个表面之间自由滚动（如研磨、抛光）
影响因素：
- 磨粒硬度：磨粒硬度越高，磨损越严重
- 磨粒尺寸：存在临界尺寸，小于此尺寸的磨粒磨损轻微
- 磨粒形状：尖锐磨粒比圆滑磨粒造成更严重的磨损
- 载荷和速度：与磨损率通常呈线性关系
典型应用：挖掘机斗齿、破碎机锤头、农业机械、泥浆泵等。
6.1.3 疲劳磨损（Fatigue Wear）
机理：接触表面在循环载荷作用下，次表面产生交变应力，导致裂纹萌生和扩展。当裂纹扩展到表面时，材料以片状剥落，形成凹坑（Pitting）或剥落（Spalling）。
特征：
- 存在明显的孕育期，裂纹萌生需要一定循环次数
- 磨损产物通常是片状或块状，与粘着磨损的颗粒状不同
- 表面硬度对疲劳磨损影响显著，适当提高硬度可以延长寿命
典型应用：滚动轴承、齿轮、凸轮-挺杆、钢轨等。
6.1.4 腐蚀磨损（Corrosive Wear）
机理：腐蚀和磨损的交互作用。腐蚀使材料表面形成疏松的腐蚀产物层，磨损将这些产物去除，暴露新鲜金属表面，加速进一步腐蚀。
分类：
- 氧化磨损：最常见的腐蚀磨损，金属表面氧化后被磨除
- 腐蚀疲劳：腐蚀与疲劳的复合作用
- 气蚀：流体中气泡溃灭产生的冲击造成的磨损
影响因素：
- 介质成分：酸、碱、盐溶液加速腐蚀
- 温度：高温加速氧化和腐蚀
- 流速：高速流体增强传质和冲刷作用
典型应用：化工设备、海洋平台、水轮机、泵阀等。
6.1.5 冲蚀磨损（Erosive Wear）
机理：固体颗粒或液滴以一定速度冲击材料表面，通过切削、挤压、疲劳等机制造成材料损失。
影响因素：
- 冲击速度：磨损率与速度的n次方成正比（n通常在2-3之间）
- 冲击角度：脆性材料和延性材料的最佳冲击角度不同
- 颗粒特性：硬度、形状、尺寸、浓度
- 靶材性能：硬度、韧性、疲劳强度
典型应用：风机叶片、除尘器、管道弯头、火箭发动机喷管等。
6.2 Archard磨损模型
Archard磨损模型是1953年由J.F. Archard提出的经典磨损模型，至今仍被广泛应用。
6.2.1 基本公式
Archard磨损方程描述了磨损体积与载荷、滑动距离的关系：
V = K F ⋅ s H V = K \frac{F \cdot s}{H} V=KHF⋅s
其中：
- V V V：磨损体积（m³）
- K K K：无量纲磨损系数（与材料配对、润滑状态、环境等有关）
- F F F：法向载荷（N）
- s s s：滑动距离（m）
- H H H：较软材料的硬度（Pa）
6.2.2 磨损系数K的物理意义
磨损系数 K K K是Archard模型的核心参数，反映了材料配对的磨损特性：
| 材料配对 | 润滑状态 | K值范围 |
|---|---|---|
| 钢/钢 | 干摩擦 | 10^-3 ~ 10^-4 |
| 钢/钢 | 油润滑 | 10^-6 ~ 10^-7 |
| 钢/青铜 | 干摩擦 | 10^-4 ~ 10^-5 |
| 钢/聚四氟乙烯 | 干摩擦 | 10^-6 ~ 10^-7 |
| 陶瓷/陶瓷 | 干摩擦 | 10^-8 ~ 10^-9 |
K K K值越小，表示材料配对的耐磨性越好。
6.2.3 磨损率与磨损深度
磨损率（单位时间的磨损体积）：
V ˙ = K F ⋅ v H \dot{V} = K \frac{F \cdot v}{H} V˙=KHF⋅v
其中， v v v是滑动速度（m/s）。
磨损深度（单位滑动距离的磨损深度）：
h = K p H ⋅ s h = K \frac{p}{H} \cdot s h=KHp⋅s
其中， p p p是接触压力（Pa）。
6.2.4 Archard模型的应用与局限
应用：
- 滑动轴承寿命预测
- 密封件磨损评估
- 制动器磨损计算
- 机械零件寿命设计
局限：
- 假设磨损率恒定，实际中磨损率可能随时间变化（跑合期、稳定期、剧烈磨损期）
- 没有考虑温度、速度、环境等动态因素的影响
- 对于疲劳磨损、腐蚀磨损等非机械磨损机理不适用
- 磨损系数K需要通过实验确定，缺乏理论预测方法
6.3 其他磨损模型
6.3.1 Usui磨损模型
Usui模型主要用于描述切削过程中的刀具磨损：
w ˙ = A ⋅ V C ⋅ exp ( − B T ) \dot{w} = A \cdot V^C \cdot \exp\left(-\frac{B}{T}\right) w˙=A⋅VC⋅exp(−TB)
其中：
- w ˙ \dot{w} w˙：磨损率
- V V V：切削速度
- T T T：切削温度
- A A A、 B B B、 C C C：材料常数
该模型考虑了切削速度和温度对磨损的影响，适用于高速切削加工。
6.3.2 Bitter冲蚀磨损模型
Bitter模型将冲蚀磨损分为变形磨损和切削磨损两个分量：
E = K 1 V n 1 sin n θ + K 2 V n 2 sin n 2 ( 2 θ ) E = K_1 V^{n_1} \sin^n\theta + K_2 V^{n_2} \sin^{n_2}(2\theta) E=K1Vn1sinnθ+K2Vn2sinn2(2θ)
其中：
- E E E：冲蚀率（单位质量磨料造成的材料损失）
- V V V：颗粒冲击速度
- θ \theta θ：冲击角度
- K 1 K_1 K1、 K 2 K_2 K2、 n 1 n_1 n1、 n 2 n_2 n2：材料常数
第一项代表变形磨损（垂直冲击为主），第二项代表切削磨损（斜向冲击为主）。
6.3.3 Rabinowicz磨粒磨损模型
Rabinowicz模型专门用于描述三体磨粒磨损：
V = C ⋅ F ⋅ s ⋅ cot ϕ 3 H V = \frac{C \cdot F \cdot s \cdot \cot\phi}{3H} V=3HC⋅F⋅s⋅cotϕ
其中：
- C C C：磨粒浓度
- ϕ \phi ϕ：磨粒锥角（通常取60度）
- 其他符号与Archard模型相同
该模型考虑了磨粒的几何形状对磨损的影响。
6.4 表面形貌与粗糙度
6.4.1 表面形貌的描述
实际工程表面都是粗糙的，表面形貌可以用高度函数 z ( x , y ) z(x,y) z(x,y)描述。表面形貌的特征包括：
粗糙峰（Asperity）：表面上的局部凸起，是接触的真实发生位置。
粗糙谷（Valley）：表面上的局部凹陷，可能储存润滑油或磨屑。
纹理（Texture）：表面加工痕迹形成的方向性特征。
6.4.2 表面粗糙度参数
算术平均粗糙度（Ra）：
R a = 1 L ∫ 0 L ∣ z ( x ) − z ˉ ∣ d x R_a = \frac{1}{L} \int_0^L |z(x) - \bar{z}| \, dx Ra=L1∫0L∣z(x)−zˉ∣dx
其中， z ˉ \bar{z} zˉ是表面平均高度。
均方根粗糙度（Rq）：
R q = 1 L ∫ 0 L ( z ( x ) − z ˉ ) 2 d x R_q = \sqrt{\frac{1}{L} \int_0^L (z(x) - \bar{z})^2 \, dx} Rq=L1∫0L(z(x)−zˉ)2dx
最大峰谷高度（Rz）：
R z = z m a x − z m i n R_z = z_{max} - z_{min} Rz=zmax−zmin
偏度（Skewness）：
S k = 1 R q 3 ⋅ 1 L ∫ 0 L ( z ( x ) − z ˉ ) 3 d x S_k = \frac{1}{R_q^3} \cdot \frac{1}{L} \int_0^L (z(x) - \bar{z})^3 \, dx Sk=Rq31⋅L1∫0L(z(x)−zˉ)3dx
偏度反映表面形貌的不对称性： S k > 0 S_k > 0 Sk>0表示表面有较多高峰， S k < 0 S_k < 0 Sk<0表示表面有较多深谷。
峰度（Kurtosis）：
K u = 1 R q 4 ⋅ 1 L ∫ 0 L ( z ( x ) − z ˉ ) 4 d x K_u = \frac{1}{R_q^4} \cdot \frac{1}{L} \int_0^L (z(x) - \bar{z})^4 \, dx Ku=Rq41⋅L1∫0L(z(x)−zˉ)4dx
峰度反映表面高度的分布特征： K u = 3 K_u = 3 Ku=3为高斯分布， K u > 3 K_u > 3 Ku>3表示表面有较多极值点。
6.4.3 表面形貌的演化
在磨损过程中，表面形貌会发生显著变化：
跑合期（Running-in）：初始阶段，表面粗糙峰快速磨平，粗糙度降低，接触状态趋于稳定。
稳定期（Steady-state）：磨损进入稳定阶段，粗糙度保持相对恒定，磨损率稳定。
剧烈磨损期（Severe Wear）：当磨损累积到一定程度，表面损伤加剧，磨损率急剧增加，最终导致失效。
表面形貌的演化对磨损行为有重要影响，粗糙度降低通常会减小磨损率，但过度光滑的表面可能导致润滑膜破裂和粘着增加。
7. 多轴疲劳与寿命预测
7.1 多轴疲劳的基本概念
在微动疲劳问题中，接触区域同时承受法向接触应力和切向摩擦应力，形成复杂的多轴应力状态。多轴疲劳（Multiaxial Fatigue）研究材料在多轴应力状态下的疲劳行为，是微动疲劳寿命预测的理论基础。
7.1.1 多轴应力状态
在三维空间中，一点的应力状态可以用应力张量表示：
σ i j = [ σ x τ x y τ x z τ y x σ y τ y z τ z x τ z y σ z ] \sigma_{ij} = \begin{bmatrix} \sigma_x & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \sigma_y & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \sigma_z \end{bmatrix} σij= σxτyxτzxτxyσyτzyτxzτyzσz
对于微动接触问题，通常可以简化为平面应力或平面应变状态。
7.1.2 等效应力理论
为了将多轴应力状态与单轴疲劳数据联系起来，需要定义等效应力（Equivalent Stress）。常用的等效应力包括：
von Mises等效应力：
σ e q = 1 2 [ ( σ 1 − σ 2 ) 2 + ( σ 2 − σ 3 ) 2 + ( σ 3 − σ 1 ) 2 ] \sigma_{eq} = \sqrt{\frac{1}{2}\left[(\sigma_1-\sigma_2)^2 + (\sigma_2-\sigma_3)^2 + (\sigma_3-\sigma_1)^2\right]} σeq=21[(σ1−σ2)2+(σ2−σ3)2+(σ3−σ1)2]
其中， σ 1 \sigma_1 σ1、 σ 2 \sigma_2 σ2、 σ 3 \sigma_3 σ3是主应力。
Tresca等效应力（最大剪应力）：
σ e q = σ 1 − σ 3 \sigma_{eq} = \sigma_1 - \sigma_3 σeq=σ1−σ3
等效应力理论假设当等效应力达到单轴疲劳极限时，材料发生疲劳失效。
7.2 临界平面法
临界平面法（Critical Plane Approach）是多轴疲劳分析的重要方法，它认为疲劳裂纹在特定的平面上萌生和扩展，该平面称为临界平面。
7.2.1 临界平面的确定
临界平面通常根据以下准则确定：
最大剪应力平面：疲劳裂纹通常在最大剪应力平面上萌生。
最大法向应力平面：对于拉伸主导的疲劳，裂纹在最大法向应力平面上扩展。
最大损伤参数平面：综合考虑剪应力和法向应力的影响，寻找使损伤参数最大的平面。
7.2.2 Fatemi-Socie损伤参数
Fatemi-Socie模型是适用于微动疲劳的临界平面法：
γ m a x ( 1 + k σ n σ y ) = C \gamma_{max} \left(1 + k \frac{\sigma_n}{\sigma_y}\right) = C γmax(1+kσyσn)=C
其中：
- γ m a x \gamma_{max} γmax：最大剪应变幅
- σ n \sigma_n σn：临界平面上的法向应力
- σ y \sigma_y σy：材料屈服强度
- k k k：材料常数（通常取0.5-1.0）
- C C C：与疲劳寿命相关的常数
该模型考虑了法向应力对剪应变疲劳的增强效应，特别适用于微动疲劳中接触应力与摩擦应力共同作用的情况。
7.2.3 Smith-Watson-Topper (SWT) 参数
SWT参数是另一种常用的多轴疲劳参数：
S W T = σ n , m a x ⋅ Δ ε 2 = σ f ′ 2 E ( 2 N f ) 2 b + σ f ′ ε f ′ ( 2 N f ) b + c SWT = \sigma_{n,max} \cdot \frac{\Delta\varepsilon}{2} = \frac{\sigma_f'^2}{E}(2N_f)^{2b} + \sigma_f'\varepsilon_f'(2N_f)^{b+c} SWT=σn,max⋅2Δε=Eσf′2(2Nf)2b+σf′εf′(2Nf)b+c
其中：
- σ n , m a x \sigma_{n,max} σn,max：最大法向应力
- Δ ε \Delta\varepsilon Δε：应变范围
- σ f ′ \sigma_f' σf′：疲劳强度系数
- ε f ′ \varepsilon_f' εf′：疲劳延性系数
- b b b、 c c c：疲劳指数
- N f N_f Nf：疲劳寿命
SWT参数适用于拉伸裂纹主导的疲劳问题。
7.3 疲劳寿命预测方法
7.3.1 应力-寿命法（S-N曲线法）
应力-寿命法基于材料的S-N曲线（Wohler曲线），描述应力幅与疲劳寿命的关系：
σ a = σ f ′ ( 2 N f ) b \sigma_a = \sigma_f' (2N_f)^b σa=σf′(2Nf)b
或取对数形式：
log ( σ a ) = log ( σ f ′ ) + b ⋅ log ( 2 N f ) \log(\sigma_a) = \log(\sigma_f') + b \cdot \log(2N_f) log(σa)=log(σf′)+b⋅log(2Nf)
其中：
- σ a \sigma_a σa：应力幅
- σ f ′ \sigma_f' σf′：疲劳强度系数
- b b b：Basquin指数（通常为-0.05到-0.12）
- N f N_f Nf：疲劳寿命（循环次数）
对于微动疲劳，需要使用考虑微动效应的等效应力。
7.3.2 应变-寿命法（ε-N曲线法）
应变-寿命法适用于低周疲劳问题，考虑塑性应变的影响：
Δ ε 2 = σ f ′ E ( 2 N f ) b + ε f ′ ( 2 N f ) c \frac{\Delta\varepsilon}{2} = \frac{\sigma_f'}{E}(2N_f)^b + \varepsilon_f'(2N_f)^c 2Δε=Eσf′(2Nf)b+εf′(2Nf)c
其中：
- Δ ε \Delta\varepsilon Δε：总应变范围
- E E E：弹性模量
- ε f ′ \varepsilon_f' εf′：疲劳延性系数
- c c c：疲劳延性指数（通常为-0.5到-0.7）
应变-寿命法可以区分弹性应变和塑性应变对疲劳的贡献。
7.3.3 裂纹萌生寿命与扩展寿命
总疲劳寿命可以分解为裂纹萌生寿命 N i N_i Ni和裂纹扩展寿命 N p N_p Np：
N f = N i + N p N_f = N_i + N_p Nf=Ni+Np
裂纹萌生寿命：从初始状态到裂纹萌生（通常定义为裂纹长度达到0.1-0.5mm）的循环次数。可以用应力-寿命法或应变-寿命法预测。
裂纹扩展寿命：从裂纹萌生到最终断裂的循环次数。可以用Paris方程预测：
d a d N = C ( Δ K ) m \frac{da}{dN} = C(\Delta K)^m dNda=C(ΔK)m
其中：
- a a a：裂纹尺寸
- Δ K \Delta K ΔK：应力强度因子范围
- C C C、 m m m：材料常数
积分Paris方程可以得到裂纹扩展寿命：
N p = ∫ a i a c d a C ( Δ K ) m N_p = \int_{a_i}^{a_c} \frac{da}{C(\Delta K)^m} Np=∫aiacC(ΔK)mda
其中， a i a_i ai是初始裂纹尺寸， a c a_c ac是临界裂纹尺寸。
7.4 微动疲劳寿命预测的特殊考虑
7.4.1 微动引起的附加应力
微动会在接触区域产生附加的应力集中，这部分应力需要叠加到远场应力上进行寿命预测。
接触应力集中系数：
K t = 1 + q ⋅ σ c o n t a c t σ b u l k K_t = 1 + q \cdot \frac{\sigma_{contact}}{\sigma_{bulk}} Kt=1+q⋅σbulkσcontact
其中， q q q是缺口敏感性系数， σ c o n t a c t \sigma_{contact} σcontact是接触应力， σ b u l k \sigma_{bulk} σbulk是远场应力。
7.4.2 滑移幅值的影响
滑移幅值对微动疲劳寿命有显著影响。通常存在一个临界滑移幅值 δ c r i t i c a l \delta_{critical} δcritical，在此滑移幅值下寿命最短。
寿命与滑移幅值的关系可以用以下经验公式描述：
N f = N 0 ( δ δ c r i t i c a l ) − n N_f = N_0 \left(\frac{\delta}{\delta_{critical}}\right)^{-n} Nf=N0(δcriticalδ)−n
其中， N 0 N_0 N0是参考寿命， n n n是经验指数（通常为1-3）。
7.4.3 磨损对寿命的影响
磨损会改变接触表面状态，影响疲劳行为。考虑磨损的寿命预测需要耦合磨损模型和疲劳模型：
d D d N = d D f a t i g u e d N + d D w e a r d N \frac{dD}{dN} = \frac{dD_{fatigue}}{dN} + \frac{dD_{wear}}{dN} dNdD=dNdDfatigue+dNdDwear
其中， D D D是总损伤， D f a t i g u e D_{fatigue} Dfatigue是疲劳损伤， D w e a r D_{wear} Dwear是磨损损伤。
当磨损率较高时，磨损去除材料的速度可能超过裂纹扩展速度，此时磨损反而有利于延长寿命。但当磨损导致应力集中加剧时，又会加速疲劳失效。
8. 实例一：微动疲劳机理与寿命预测
8.1 实例概述
本实例通过Python编程实现微动疲劳的数值仿真，包括以下核心内容：
- Hertz接触力学计算：计算接触压力分布、接触半宽、次表面应力场
- Mindlin局部滑移分析：确定粘着区和滑移区分布，计算切向牵引力
- 微动迟滞回线模拟：绘制不同滑移幅值下的Q-δ曲线，计算能量耗散
- 微动图生成：构建位移-力空间的运行状态图
- 疲劳寿命预测：基于多轴疲劳理论预测裂纹萌生寿命和扩展寿命
8.2 代码结构说明
代码采用面向对象的设计，主要包含以下类：
ContactMechanics类：实现Hertz接触理论，包括圆柱接触和球接触的压力分布计算。
FrettingMechanics类：实现Mindlin局部滑移理论，计算切向牵引力分布和微动迟滞回线。
FrettingFatigue类：实现微动疲劳寿命预测，包括裂纹萌生和扩展寿命计算。
FrettingMap类：实现微动图的生成和运行状态判断。
8.3 核心算法解析
8.3.1 Hertz接触压力计算
def hertz_cylinder(self, R, P, L):
"""圆柱-平面Hertz接触计算"""
# 接触半宽
b = np.sqrt(4 * P * R / (np.pi * self.E_star))
# 最大接触压力
p0 = 2 * P / (np.pi * b)
# 接触压力分布（半椭圆）
p = p0 * np.sqrt(1 - (x/b)**2)
该函数基于Hertz理论计算圆柱接触的接触半宽和最大压力。接触压力呈半椭圆分布，这是Hertz接触的基本特征。
8.3.2 Mindlin滑移区计算
def _calculate_mindlin_zones(self):
"""计算Mindlin局部滑移区"""
# 切向力系数
self.eta = self.Q_max / (self.mu * self.P)
# 滑移区半宽
self.c = self.b * np.sqrt(1 - self.eta)
根据Mindlin理论，滑移区半宽 c c c与接触半宽 b b b的关系取决于切向力系数 η \eta η。当 η → 0 \eta \to 0 η→0时， c → b c \to b c→b，整个接触区粘着；当 η → 1 \eta \to 1 η→1时， c → 0 c \to 0 c→0，整个接触区滑移。
8.3.3 微动迟滞回线计算
def fretting_loop(self, delta_max, n_points=100):
"""计算微动迟滞回线"""
delta = np.linspace(-delta_max, delta_max, n_points)
Q = np.zeros_like(delta)
for i, d in enumerate(delta):
if abs(d) <= delta_max * 0.5: # 部分滑移区
Q[i] = self.mu * self.P * (1 - (1 - abs(d)/delta_max)**1.5)
else: # 完全滑移区
Q[i] = np.sign(d) * self.mu * self.P
# 计算能量耗散（回线面积）
energy_dissipation = np.trapezoid(Q, delta)
该函数模拟了微动过程中的力-位移关系。在部分滑移区，力-位移呈非线性关系；在完全滑移区，力保持恒定（等于极限摩擦力）。能量耗散等于迟滞回线的面积，是评估微动损伤的重要指标。
8.3.4 疲劳寿命预测
def total_fretting_fatigue_life(self, contact_pressure, slip_amplitude, bulk_stress=100e6):
"""计算微动疲劳总寿命"""
# 计算微动引起的附加应力
fretting_stress = 0.5 * contact_pressure * (slip_amplitude / 1e-6)**0.3
# 总应力
total_stress = bulk_stress + fretting_stress
delta_sigma = 2 * total_stress
# 裂纹萌生寿命（Basquin方程）
N_i = self.crack_initiation_life(delta_sigma, delta_tau)
# 裂纹扩展寿命（Paris方程积分）
N_p = self.crack_propagation_life(a_initial, a_critical, delta_sigma)
# 总寿命
N_total = N_i + N_p
该函数综合考虑了远场应力和微动引起的附加应力，使用Basquin方程预测裂纹萌生寿命，使用Paris方程预测裂纹扩展寿命。
8.4 可视化结果分析
8.4.1 Hertz接触压力分布
该图展示了圆柱接触的表面压力分布，呈半椭圆形。图中标注了接触边界和粘着区边界，帮助理解Mindlin局部滑移的概念。
物理意义：接触压力在接触中心最大，向边缘逐渐减小至零。粘着区位于接触中心，滑移区位于接触边缘。
8.4.2 微动迟滞回线
该图展示了不同滑移幅值下的力-位移曲线。随着滑移幅值增大，迟滞回线从细长的叶片状逐渐变为矩形。
物理意义：小滑移幅值对应部分滑移状态，能量耗散小；大滑移幅值对应完全滑移状态，能量耗散大。能量耗散与磨损和疲劳损伤密切相关。
8.4.3 微动图
微动图以位移和力为坐标轴，将运行状态划分为部分滑移区、混合滑移区和完全滑移区。
工程应用：通过微动图可以判断当前运行状态，为设计改进提供指导。例如，应避免长期处于部分滑移区，因为该区域疲劳寿命最短。
8.4.4 疲劳寿命曲线
该图展示了接触压力和滑移幅值对疲劳寿命的影响。寿命随接触压力增加而降低，随滑移幅值增大呈现先降低后增加的趋势。
物理意义：存在一个临界滑移幅值，在此滑移幅值下寿命最短。这是微动疲劳的重要特征，也是防护措施设计的依据。
8.5 运行结果预期
运行代码后，将生成以下文件：
-
fretting_fatigue_analysis.png：包含9个子图的综合分析图，涵盖接触压力、迟滞回线、能量耗散、微动图、寿命预测等内容。
-
fretting_process_animation.gif：微动过程动画，展示一个微动循环中接触区滑移状态的演变。
-
contact_stress_analysis.png：接触应力分析图，包括表面压力分布、次表面应力、切向牵引力分布等。
9. 实例二：磨损模型与表面演化模拟
9.1 实例概述
本实例通过Python编程实现磨损过程的数值仿真，包括以下核心内容：
- Archard磨损模型：计算磨损体积、磨损率和磨损深度
- 其他磨损模型：Usui模型（切削磨损）、Bitter模型（冲蚀磨损）
- 表面形貌演化：模拟磨损过程中表面粗糙度的变化
- 磨损机制分析：对比粘着磨损、磨粒磨损等不同机制
- 磨损机制图：构建载荷-速度空间的磨损机制分布图
9.2 代码结构说明
代码采用面向对象的设计，主要包含以下类：
ArchardWearModel类：实现Archard磨损模型，计算磨损体积和磨损率。
UsuiWearModel类：实现Usui磨损模型，用于切削过程中的刀具磨损预测。
BitterWearModel类：实现Bitter冲蚀磨损模型，描述颗粒冲击引起的材料损失。
SurfaceEvolution类：实现表面形貌的生成、演化分析和粗糙度计算。
WearSimulation类：实现磨损过程的时域仿真，记录表面形貌和粗糙度的演化历史。
AdhesiveWearSimulation类：实现粘着磨损的详细模拟，包括接合点生长和磨屑形成。
AbrasiveWearSimulation类：实现磨粒磨损的详细模拟，包括切削体积计算。
9.3 核心算法解析
9.3.1 Archard磨损计算
def wear_volume(self, F, s):
"""计算磨损体积"""
V = self.K * F * s / self.H
return V
def wear_depth_increment(self, p, delta_s):
"""计算局部磨损深度增量"""
delta_h = self.K * p * delta_s / self.H
return delta_h
Archard模型是磨损仿真的基础，磨损体积与载荷、滑动距离成正比，与材料硬度成反比。局部磨损深度计算考虑了接触压力的不均匀分布。
9.3.2 表面形貌生成
def generate_initial_surface(self, Ra=0.5e-6, correlation_length=10e-6):
"""生成初始粗糙表面"""
# 生成白噪声
noise = np.random.randn(self.size_y, self.size_x)
# 高斯滤波模拟相关性
sigma = correlation_length / self.dx / 2.355
self.height = ndimage.gaussian_filter(noise, sigma=sigma)
# 归一化到目标粗糙度
current_Ra = np.mean(np.abs(self.height))
self.height = self.height * Ra / current_Ra
该函数使用高斯滤波白噪声的方法生成粗糙表面。相关长度控制表面的平滑程度，粗糙度Ra控制表面的高低起伏。
9.3.3 磨损过程仿真
def step(self, dt, v_slide, F_normal):
"""执行一个时间步长的磨损仿真"""
# 计算滑动距离
delta_s = v_slide * dt
# 计算接触压力分布
pressure = self.calculate_contact_pressure(F_normal)
# 计算磨损深度增量
wear_depth = self.wear_model.wear_depth_increment(pressure, delta_s)
# 应用磨损
self.surface.apply_wear(wear_depth)
# 记录历史
self.history['roughness'].append(self.surface.calculate_roughness())
该函数实现了磨损过程的时间推进。在每个时间步长内，根据滑动距离和接触压力计算磨损深度，更新表面形貌，并记录粗糙度变化。
9.3.4 粗糙度计算
def calculate_roughness(self):
"""计算表面粗糙度参数"""
Ra = np.mean(np.abs(self.height)) # 算术平均粗糙度
Rq = np.sqrt(np.mean(self.height**2)) # 均方根粗糙度
Rz = np.max(self.height) - np.min(self.height) # 最大峰谷高度
skewness = np.mean((self.height / Rq)**3) # 偏度
kurtosis = np.mean((self.height / Rq)**4) # 峰度
该函数计算了常用的表面粗糙度参数，包括Ra、Rq、Rz、偏度和峰度。这些参数全面描述了表面形貌的特征。
9.4 可视化结果分析
9.4.1 磨损模型对比
该图展示了不同材料配对的Archard磨损系数和磨损体积随滑动距离的变化。钢/钢配对的磨损系数最高，陶瓷/陶瓷配对的磨损系数最低。
工程意义：选择合适的材料配对可以显著降低磨损。例如，在需要低摩擦和耐磨的场合，可以选择钢/聚四氟乙烯配对。
9.4.2 表面形貌演化
该图展示了磨损过程中表面形貌的变化。随着磨损进行，表面粗糙度逐渐降低，表面趋于光滑。
物理意义：磨损初期，粗糙峰快速磨平，粗糙度降低；磨损后期，表面达到动态平衡，粗糙度趋于稳定。
9.4.3 粗糙度演化曲线
该图展示了算术平均粗糙度Ra和均方根粗糙度Rq随时间的变化。粗糙度先快速降低，然后趋于稳定。
工程应用：通过监测粗糙度变化，可以判断磨损所处的阶段（跑合期、稳定期、剧烈磨损期），为维修决策提供依据。
9.4.4 磨损机制图
该图以载荷和速度为坐标轴，将磨损机制划分为粘着磨损、磨粒磨损、氧化磨损和熔融磨损四个区域。
设计指导：通过控制载荷和速度，可以使运行状态避开高磨损区域。例如，高载荷高速度区域容易发生熔融磨损，应尽量避免。
9.5 运行结果预期
运行代码后，将生成以下文件：
-
wear_models_analysis.png：磨损模型分析图，包括Archard磨损曲线、磨损率对比、冲蚀磨损角度效应等。
-
surface_evolution_analysis.png：表面形貌演化分析图，展示不同时刻的表面形貌、粗糙度演化、磨损体积累积等。
-
surface_3d_visualization.png：3D表面形貌图，直观展示初始表面和磨损后表面的形貌对比。
-
wear_mechanisms_comparison.png：磨损机制对比图，包括粘着磨损率、磨粒磨损率、磨屑尺寸分布、磨损机制图等。
-
surface_evolution_animation.gif：表面演化动画，动态展示磨损过程中表面形貌和粗糙度的变化。
# -*- coding: utf-8 -*-
"""
主题088：微动疲劳与磨损 - 实例一
微动疲劳机理与寿命预测
本实例模拟微动疲劳（Fretting Fatigue）的机理和寿命预测。
微动疲劳发生在接触表面之间发生微小相对运动（微米级）的情况下，
常见于螺栓连接、花键配合、叶片-盘连接等工程结构。
应用场景：
- 航空发动机叶片-盘连接
- 汽车螺栓连接
- 铁路轨道扣件
- 医疗器械植入物
作者：结构耐久性仿真教程编写组
日期：2026年3月
"""
import numpy as np
import matplotlib
matplotlib.use('Agg') # 使用无头模式
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
import io
from PIL import Image
from scipy import integrate
from scipy.optimize import fsolve
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("=" * 60)
print("主题088：微动疲劳与磨损 - 实例一")
print("微动疲劳机理与寿命预测")
print("=" * 60)
class ContactMechanics:
"""
接触力学模型
基于Hertz接触理论计算接触区域的应力和变形
"""
def __init__(self, E1=200e9, nu1=0.3, E2=200e9, nu2=0.3):
"""
初始化接触力学参数
Parameters:
-----------
E1, E2 : float
两接触体的弹性模量 (Pa)
nu1, nu2 : float
两接触体的泊松比
"""
self.E1 = E1
self.nu1 = nu1
self.E2 = E2
self.nu2 = nu2
# 计算等效弹性模量
self.E_star = self._equivalent_modulus()
def _equivalent_modulus(self):
"""
计算等效弹性模量
Returns:
--------
E_star : float
等效弹性模量 (Pa)
"""
E_star = 1 / ((1 - self.nu1**2) / self.E1 + (1 - self.nu2**2) / self.E2)
return E_star
def hertz_cylinder(self, R, P, L):
"""
圆柱-平面Hertz接触计算
Parameters:
-----------
R : float
圆柱半径 (m)
P : float
单位长度载荷 (N/m)
L : float
接触长度 (m)
Returns:
--------
contact_params : dict
包含接触半宽、最大压力等参数
"""
# 接触半宽
b = np.sqrt(4 * P * R / (np.pi * self.E_star))
# 最大接触压力
p0 = 2 * P / (np.pi * b)
# 接触面积
A = 2 * b * L
# 平均接触压力
p_avg = P / (2 * b)
return {
'contact_half_width': b,
'max_pressure': p0,
'contact_area': A,
'avg_pressure': p_avg,
}
def hertz_sphere(self, R, F):
"""
球-平面Hertz接触计算
Parameters:
-----------
R : float
球半径 (m)
F : float
法向载荷 (N)
Returns:
--------
contact_params : dict
包含接触半径、最大压力等参数
"""
# 接触半径
a = (3 * F * R / (4 * self.E_star)) ** (1/3)
# 最大接触压力
p0 = 3 * F / (2 * np.pi * a**2)
# 接触面积
A = np.pi * a**2
# 接近量
delta = a**2 / R
return {
'contact_radius': a,
'max_pressure': p0,
'contact_area': A,
'approach': delta,
}
def subsurface_stress_cylinder(self, x, y, b, p0):
"""
计算圆柱接触下的次表面应力场
Parameters:
-----------
x, y : array
坐标位置 (m)，原点在接触中心
b : float
接触半宽 (m)
p0 : float
最大接触压力 (Pa)
Returns:
--------
stress : dict
应力分量 sx, sy, sz, txy
"""
# 无量纲坐标
X = x / b
Y = y / b
# Hertz接触应力场公式（简化版）
# 注意：这是基于弹性半空间的解析解
# 表面法向压力分布
p = p0 * np.sqrt(1 - X**2) if np.abs(X) <= 1 else 0
# 次表面应力（基于Boussinesq解）
# 这里使用简化公式
# 最大剪应力深度
z_max_shear = 0.786 * b
# 在深度z_max_shear处的应力
sigma_x = -0.3 * p0 # 近似值
sigma_y = -0.8 * p0 # 近似值
sigma_z = -0.5 * p0 # 近似值
tau_xy = 0.3 * p0 # 最大剪应力
return {
'sigma_x': sigma_x,
'sigma_y': sigma_y,
'sigma_z': sigma_z,
'tau_xy': tau_xy,
'z_max_shear': z_max_shear,
}
class FrettingMechanics:
"""
微动力学模型
模拟微动过程中的力学行为，包括局部滑移和能量耗散
"""
def __init__(self, mu=0.5, R=0.01, P=1e6, Q_max=0.5e6):
"""
初始化微动参数
Parameters:
-----------
mu : float
摩擦系数
R : float
圆柱半径 (m)
P : float
单位长度法向载荷 (N/m)
Q_max : float
切向力幅值 (N/m)
"""
self.mu = mu
self.R = R
self.P = P
self.Q_max = Q_max
# 计算接触参数
self.contact = ContactMechanics()
hertz_params = self.contact.hertz_cylinder(R, P, 1.0)
self.b = hertz_params['contact_half_width']
self.p0 = hertz_params['max_pressure']
# 计算Mindlin滑移区
self._calculate_mindlin_zones()
def _calculate_mindlin_zones(self):
"""
计算Mindlin局部滑移区
Mindlin理论描述了切向力作用下接触区域的粘着区和滑移区分布
"""
# 切向力系数
self.eta = self.Q_max / (self.mu * self.P)
# 滑移区半宽
self.c = self.b * np.sqrt(1 - self.eta)
# 粘着区半宽
self.adhesion_zone = self.c
self.slip_zone_outer = self.b
def tangential_traction(self, x, Q):
"""
计算切向牵引力分布
Parameters:
-----------
x : array
位置坐标 (m)
Q : float
当前切向力 (N/m)
Returns:
--------
q : array
切向牵引力分布 (Pa)
"""
X = x / self.b
eta = Q / (self.mu * self.P)
# Mindlin切向牵引力分布
q = np.zeros_like(x)
mask = np.abs(X) <= 1
if np.any(mask):
c_over_b = np.sqrt(1 - eta)
# 粘着区牵引力
q[mask] = self.mu * self.p0 * (
np.sqrt(1 - X[mask]**2) -
c_over_b * np.sqrt(1 - (X[mask]/c_over_b)**2)
)
return q
def fretting_loop(self, delta_max, n_points=100):
"""
计算微动迟滞回线（Fretting Loop）
微动迟滞回线描述了切向力Q与相对位移δ的关系
Parameters:
-----------
delta_max : float
最大相对位移幅值 (m)
n_points : int
计算点数
Returns:
--------
loop_data : dict
包含位移、切向力、能量耗散等数据
"""
# 位移范围
delta = np.linspace(-delta_max, delta_max, n_points)
# 计算对应的切向力
Q = np.zeros_like(delta)
for i, d in enumerate(delta):
# 简化的Q-δ关系（基于Mindlin理论）
if abs(d) <= delta_max * 0.5: # 部分滑移区
Q[i] = self.mu * self.P * (1 - (1 - abs(d)/delta_max)**1.5)
else: # 完全滑移区
Q[i] = np.sign(d) * self.mu * self.P
# 计算能量耗散（回线面积）
energy_dissipation = np.trapezoid(Q, delta)
return {
'displacement': delta,
'tangential_force': Q,
'energy_dissipation': energy_dissipation,
'loop_area': abs(energy_dissipation),
}
def wear_rate_archard(self, slip_amplitude, N_cycles, K=1e-8, H=2e9):
"""
Archard磨损模型
Parameters:
-----------
slip_amplitude : float
滑移幅值 (m)
N_cycles : int
循环次数
K : float
无量纲磨损系数
H : float
材料硬度 (Pa)
Returns:
--------
wear_depth : float
磨损深度 (m)
"""
# 法向压力
p = self.P / (2 * self.b)
# Archard磨损方程
wear_volume = K * p * (2 * slip_amplitude) * N_cycles / H
# 假设磨损均匀分布在接触面积上
wear_depth = wear_volume / (2 * self.b)
return wear_depth
class FrettingFatigue:
"""
微动疲劳寿命预测模型
结合微动磨损和疲劳损伤，预测微动疲劳寿命
"""
def __init__(self, material='steel'):
"""
初始化微动疲劳参数
Parameters:
-----------
material : str
材料类型
"""
self.material = material
if material == 'steel':
self.props = {
'fatigue_strength': 400e6, # 疲劳强度 (Pa)
'fatigue_exponent': -0.085, # Basquin指数
'fracture_toughness': 60e6, # 断裂韧性 (Pa·m^0.5)
'hardness': 2e9, # 硬度 (Pa)
'wear_coefficient': 1e-8, # 磨损系数
}
elif material == 'titanium':
self.props = {
'fatigue_strength': 600e6,
'fatigue_exponent': -0.08,
'fracture_toughness': 50e6,
'hardness': 3e9,
'wear_coefficient': 5e-9,
}
elif material == 'aluminum':
self.props = {
'fatigue_strength': 150e6,
'fatigue_exponent': -0.1,
'fracture_toughness': 30e6,
'hardness': 1e9,
'wear_coefficient': 2e-8,
}
def multiaxial_fatigue_parameter(self, sigma_n, tau, sigma_f=400e6):
"""
计算多轴疲劳参数
使用Smith-Watson-Topper (SWT) 或 Fatemi-Socie 参数
Parameters:
-----------
sigma_n : float
法向应力幅值 (Pa)
tau : float
切向应力幅值 (Pa)
sigma_f : float
疲劳强度系数 (Pa)
Returns:
--------
damage_param : float
损伤参数
"""
# Fatemi-Socie参数（适用于微动疲劳）
gamma_max = tau / (2 * 80e9) # 假设剪切模量G=80GPa
# Fatemi-Socie损伤参数
k = 0.5 # 材料常数
damage_param = gamma_max * (1 + k * sigma_n / sigma_f)
return damage_param
def crack_initiation_life(self, delta_sigma, delta_tau):
"""
预测裂纹萌生寿命
Parameters:
-----------
delta_sigma : float
法向应力范围 (Pa)
delta_tau : float
切向应力范围 (Pa)
Returns:
--------
N_i : int
裂纹萌生寿命（循环次数）
"""
# 等效应力范围
delta_sigma_eq = np.sqrt(delta_sigma**2 + 3 * delta_tau**2)
# Basquin方程
sigma_f_prime = 1.5 * self.props['fatigue_strength']
b = self.props['fatigue_exponent']
if delta_sigma_eq >= sigma_f_prime:
return 1
N_i = (delta_sigma_eq / sigma_f_prime) ** (1/b)
return max(1, int(N_i))
def crack_propagation_life(self, a_initial, a_critical, delta_sigma):
"""
预测裂纹扩展寿命
Parameters:
-----------
a_initial : float
初始裂纹尺寸 (m)
a_critical : float
临界裂纹尺寸 (m)
delta_sigma : float
应力范围 (Pa)
Returns:
--------
N_p : int
裂纹扩展寿命（循环次数）
"""
# Paris方程参数
C = 1e-11
m = 3.0
# 几何因子
Y = 1.12
# 积分Paris方程
def integrand(a):
delta_K = Y * delta_sigma * np.sqrt(np.pi * a)
return 1 / (C * delta_K**m)
try:
N_p = integrate.quad(integrand, a_initial, a_critical)[0]
except:
N_p = 1e6 # 默认值
return int(N_p)
def total_fretting_fatigue_life(self, contact_pressure, slip_amplitude,
bulk_stress=100e6):
"""
计算微动疲劳总寿命
Parameters:
-----------
contact_pressure : float
接触压力 (Pa)
slip_amplitude : float
滑移幅值 (m)
bulk_stress : float
远场应力 (Pa)
Returns:
--------
life : dict
包含裂纹萌生寿命、扩展寿命和总寿命
"""
# 计算微动引起的附加应力
# 简化模型：微动应力与接触压力和滑移有关
fretting_stress = 0.5 * contact_pressure * (slip_amplitude / 1e-6)**0.3
# 总应力
total_stress = bulk_stress + fretting_stress
delta_sigma = 2 * total_stress # 假设完全反向加载
delta_tau = 0.5 * delta_sigma # 假设切应力比例
# 裂纹萌生寿命
N_i = self.crack_initiation_life(delta_sigma, delta_tau)
# 裂纹扩展寿命
a_initial = 0.1e-3 # 100微米
K_IC = self.props['fracture_toughness']
a_critical = (K_IC / (1.12 * total_stress))**2 / np.pi
N_p = self.crack_propagation_life(a_initial, a_critical, delta_sigma)
# 总寿命
N_total = N_i + N_p
return {
'initiation_life': N_i,
'propagation_life': N_p,
'total_life': N_total,
'fretting_stress': fretting_stress,
}
class FrettingMap:
"""
微动图（Fretting Map）
描述不同运行条件下的微动行为（部分滑移、混合滑移、完全滑移）
"""
def __init__(self, P=1e6, mu=0.5):
"""
初始化微动图
Parameters:
-----------
P : float
法向载荷 (N/m)
mu : float
摩擦系数
"""
self.P = P
self.mu = mu
self.Q_max_crit = mu * P # 完全滑移临界切向力
def calculate_regime(self, delta, Q_max):
"""
确定微动运行状态
Parameters:
-----------
delta : float
相对位移幅值 (m)
Q_max : float
切向力幅值 (N/m)
Returns:
--------
regime : str
运行状态：'partial_slip', 'mixed', 'gross_slip'
"""
eta = Q_max / self.Q_max_crit
# 简化的判据
if eta < 0.3:
return 'partial_slip'
elif eta < 0.8:
return 'mixed'
else:
return 'gross_slip'
def generate_map(self, delta_range, Q_range):
"""
生成微动图
Parameters:
-----------
delta_range : tuple
位移范围 (min, max) in meters
Q_range : tuple
切向力范围 (min, max) in N/m
Returns:
--------
map_data : dict
包含网格数据和区域分类
"""
n_points = 50
deltas = np.linspace(delta_range[0], delta_range[1], n_points)
Qs = np.linspace(Q_range[0], Q_range[1], n_points)
Delta, Q_grid = np.meshgrid(deltas, Qs)
Regime = np.zeros_like(Delta, dtype=int)
for i in range(n_points):
for j in range(n_points):
regime = self.calculate_regime(Delta[i,j], Q_grid[i,j])
if regime == 'partial_slip':
Regime[i,j] = 1
elif regime == 'mixed':
Regime[i,j] = 2
else:
Regime[i,j] = 3
return {
'displacement': Delta,
'tangential_force': Q_grid,
'regime': Regime,
}
def create_visualization():
"""
创建微动疲劳分析的可视化
"""
# 创建接触力学实例
contact = ContactMechanics(E1=200e9, nu1=0.3, E2=200e9, nu2=0.3)
# 创建微动力学实例
fretting = FrettingMechanics(mu=0.5, R=0.01, P=1e6, Q_max=0.3e6)
# 创建微动疲劳实例
fatigue_steel = FrettingFatigue('steel')
fatigue_titanium = FrettingFatigue('titanium')
fatigue_aluminum = FrettingFatigue('aluminum')
# 创建微动图实例
fretting_map = FrettingMap(P=1e6, mu=0.5)
# 创建大图
fig = plt.figure(figsize=(20, 16))
# 1. Hertz接触压力分布
ax1 = plt.subplot(3, 3, 1)
x = np.linspace(-2*fretting.b, 2*fretting.b, 100)
X = x / fretting.b
# 半椭圆压力分布
p = np.zeros_like(x)
mask = np.abs(X) <= 1
p[mask] = fretting.p0 * np.sqrt(1 - X[mask]**2)
ax1.fill_between(x*1000, p/1e6, alpha=0.3, color='blue')
ax1.plot(x*1000, p/1e6, 'b-', linewidth=2)
ax1.axvline(x=-fretting.b*1000, color='r', linestyle='--', alpha=0.5, label='接触边界')
ax1.axvline(x=fretting.b*1000, color='r', linestyle='--', alpha=0.5)
ax1.axvline(x=-fretting.c*1000, color='g', linestyle=':', alpha=0.7, label='粘着区边界')
ax1.axvline(x=fretting.c*1000, color='g', linestyle=':', alpha=0.7)
ax1.set_xlabel('位置 x (mm)', fontsize=11)
ax1.set_ylabel('接触压力 (MPa)', fontsize=11)
ax1.set_title('Hertz接触压力分布与Mindlin滑移区', fontsize=12, fontweight='bold')
ax1.legend(loc='best')
ax1.grid(True, alpha=0.3)
# 2. 微动迟滞回线
ax2 = plt.subplot(3, 3, 2)
delta_max_values = [1e-6, 5e-6, 10e-6, 20e-6]
colors = ['blue', 'green', 'orange', 'red']
for delta_max, color in zip(delta_max_values, colors):
loop = fretting.fretting_loop(delta_max)
ax2.plot(loop['displacement']*1e6, loop['tangential_force']/1e3,
color=color, linewidth=2, label=f'δ={delta_max*1e6:.0f}μm')
ax2.set_xlabel('相对位移 δ (μm)', fontsize=11)
ax2.set_ylabel('切向力 Q (kN/m)', fontsize=11)
ax2.set_title('微动迟滞回线（Fretting Loop）', fontsize=12, fontweight='bold')
ax2.legend(loc='best')
ax2.grid(True, alpha=0.3)
# 3. 能量耗散 vs 滑移幅值
ax3 = plt.subplot(3, 3, 3)
delta_range = np.linspace(0.5e-6, 50e-6, 50)
energies = []
for delta in delta_range:
loop = fretting.fretting_loop(delta)
energies.append(loop['loop_area'])
ax3.semilogy(delta_range*1e6, energies, 'b-', linewidth=2)
ax3.set_xlabel('滑移幅值 (μm)', fontsize=11)
ax3.set_ylabel('能量耗散 (J/m)', fontsize=11)
ax3.set_title('每循环能量耗散', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
# 4. 微动图（Fretting Map）
ax4 = plt.subplot(3, 3, 4)
map_data = fretting_map.generate_map(
delta_range=(0.1e-6, 30e-6),
Q_range=(0, 0.6e6)
)
cmap = plt.cm.colors.ListedColormap(['lightblue', 'yellow', 'salmon'])
im = ax4.contourf(map_data['displacement']*1e6,
map_data['tangential_force']/1e3,
map_data['regime'],
levels=[0.5, 1.5, 2.5, 3.5],
colors=['lightblue', 'yellow', 'salmon'])
# 添加边界线
ax4.contour(map_data['displacement']*1e6,
map_data['tangential_force']/1e3,
map_data['regime'],
levels=[1.5, 2.5], colors='black', linewidths=2)
ax4.set_xlabel('相对位移幅值 (μm)', fontsize=11)
ax4.set_ylabel('切向力幅值 (kN/m)', fontsize=11)
ax4.set_title('微动图（Fretting Map）', fontsize=12, fontweight='bold')
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor='lightblue', label='部分滑移'),
Patch(facecolor='yellow', label='混合滑移'),
Patch(facecolor='salmon', label='完全滑移')
]
ax4.legend(handles=legend_elements, loc='upper left')
# 5. 材料疲劳寿命对比
ax5 = plt.subplot(3, 3, 5)
contact_pressures = np.linspace(100e6, 1000e6, 20)
slip_amp = 10e-6
lives_steel = []
lives_titanium = []
lives_aluminum = []
for p in contact_pressures:
life_steel = fatigue_steel.total_fretting_fatigue_life(p, slip_amp)
life_titanium = fatigue_titanium.total_fretting_fatigue_life(p, slip_amp)
life_aluminum = fatigue_aluminum.total_fretting_fatigue_life(p, slip_amp)
lives_steel.append(life_steel['total_life'])
lives_titanium.append(life_titanium['total_life'])
lives_aluminum.append(life_aluminum['total_life'])
ax5.semilogy(contact_pressures/1e6, lives_steel, 'b-', linewidth=2, marker='o', label='钢')
ax5.semilogy(contact_pressures/1e6, lives_titanium, 'r-', linewidth=2, marker='s', label='钛合金')
ax5.semilogy(contact_pressures/1e6, lives_aluminum, 'g-', linewidth=2, marker='^', label='铝合金')
ax5.set_xlabel('接触压力 (MPa)', fontsize=11)
ax5.set_ylabel('疲劳寿命 (cycles)', fontsize=11)
ax5.set_title('不同材料的微动疲劳寿命 (δ=10μm)', fontsize=12, fontweight='bold')
ax5.legend(loc='best')
ax5.grid(True, alpha=0.3)
# 6. 滑移幅值对寿命的影响
ax6 = plt.subplot(3, 3, 6)
slip_amplitudes = np.linspace(1e-6, 50e-6, 30)
p_contact = 500e6
lives_vs_slip = []
for delta in slip_amplitudes:
life = fatigue_steel.total_fretting_fatigue_life(p_contact, delta)
lives_vs_slip.append(life['total_life'])
ax6.semilogy(slip_amplitudes*1e6, lives_vs_slip, 'b-', linewidth=2)
ax6.axvline(x=10, color='r', linestyle='--', alpha=0.5, label='典型滑移幅值')
ax6.set_xlabel('滑移幅值 (μm)', fontsize=11)
ax6.set_ylabel('疲劳寿命 (cycles)', fontsize=11)
ax6.set_title('滑移幅值对微动疲劳寿命的影响', fontsize=12, fontweight='bold')
ax6.legend(loc='best')
ax6.grid(True, alpha=0.3)
# 7. 磨损深度演化
ax7 = plt.subplot(3, 3, 7)
cycles = np.linspace(0, 1e6, 100)
slip_amps = [5e-6, 10e-6, 20e-6]
colors = ['blue', 'green', 'red']
for slip_amp, color in zip(slip_amps, colors):
wear_depths = [fretting.wear_rate_archard(slip_amp, int(N)) * 1e6
for N in cycles]
ax7.plot(cycles/1000, wear_depths, color=color, linewidth=2,
label=f'δ={slip_amp*1e6:.0f}μm')
ax7.set_xlabel('循环次数 (×10³)', fontsize=11)
ax7.set_ylabel('磨损深度 (μm)', fontsize=11)
ax7.set_title('Archard磨损演化', fontsize=12, fontweight='bold')
ax7.legend(loc='best')
ax7.grid(True, alpha=0.3)
# 8. 裂纹萌生vs扩展寿命
ax8 = plt.subplot(3, 3, 8)
pressures = np.linspace(200e6, 800e6, 15)
slip_amp = 10e-6
initiation_lives = []
propagation_lives = []
for p in pressures:
life = fatigue_steel.total_fretting_fatigue_life(p, slip_amp)
initiation_lives.append(life['initiation_life'])
propagation_lives.append(life['propagation_life'])
width = 20
x_pos = np.arange(len(pressures))
ax8.bar(x_pos, np.log10(initiation_lives), width, label='裂纹萌生', color='steelblue')
ax8.bar(x_pos, np.log10(propagation_lives), width,
bottom=np.log10(initiation_lives), label='裂纹扩展', color='coral')
ax8.set_xlabel('接触压力工况', fontsize=11)
ax8.set_ylabel('log₁₀(寿命)', fontsize=11)
ax8.set_title('裂纹萌生寿命 vs 扩展寿命', fontsize=12, fontweight='bold')
ax8.set_xticks(x_pos[::3])
ax8.set_xticklabels([f'{p/1e6:.0f}' for p in pressures[::3]])
ax8.legend(loc='best')
ax8.grid(True, alpha=0.3, axis='y')
# 9. 综合损伤评估
ax9 = plt.subplot(3, 3, 9)
# 创建损伤评估矩阵
p_values = np.linspace(200, 800, 8) # MPa
delta_values = np.linspace(2, 40, 8) # μm
damage_matrix = np.zeros((len(delta_values), len(p_values)))
for i, delta in enumerate(delta_values):
for j, p in enumerate(p_values):
life = fatigue_steel.total_fretting_fatigue_life(p*1e6, delta*1e-6)
# 使用Miner法则计算损伤（假设1e6次循环）
D = 1e6 / life['total_life'] if life['total_life'] > 0 else 1
damage_matrix[i, j] = min(D, 1)
im = ax9.imshow(damage_matrix, cmap='RdYlGn_r', aspect='auto', vmin=0, vmax=1)
ax9.set_xticks(range(0, len(p_values), 2))
ax9.set_yticks(range(0, len(delta_values), 2))
ax9.set_xticklabels([f'{p:.0f}' for p in p_values[::2]])
ax9.set_yticklabels([f'{d:.0f}' for d in delta_values[::2]])
ax9.set_xlabel('接触压力 (MPa)', fontsize=11)
ax9.set_ylabel('滑移幅值 (μm)', fontsize=11)
ax9.set_title('微动疲劳损伤评估矩阵', fontsize=12, fontweight='bold')
cbar = plt.colorbar(im, ax=ax9)
cbar.set_label('损伤度 D', fontsize=10)
plt.tight_layout()
plt.savefig('d:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/fretting_fatigue_analysis.png',
dpi=150, bbox_inches='tight')
print("微动疲劳分析图已保存")
plt.close()
# 创建微动过程动画
create_fretting_animation(fretting)
# 创建接触应力分布图
create_contact_stress_distribution(contact, fretting)
def create_fretting_animation(fretting):
"""
创建微动过程动画
"""
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# 一个微动循环的相位角
phases = np.linspace(0, 2*np.pi, 50)
frames = []
for phase in phases:
ax1.clear()
ax2.clear()
# 计算当前位移和切向力
delta_max = 10e-6
delta = delta_max * np.sin(phase)
# 简化的切向力计算
eta = 0.5
Q_max = eta * fretting.mu * fretting.P
Q = Q_max * np.sin(phase)
# 左图：微动迟滞回线，显示当前点
loop = fretting.fretting_loop(delta_max)
ax1.plot(loop['displacement']*1e6, loop['tangential_force']/1e3,
'b-', linewidth=2, alpha=0.5)
ax1.scatter([delta*1e6], [Q/1e3], color='red', s=150, zorder=5)
ax1.set_xlabel('相对位移 (μm)', fontsize=11)
ax1.set_ylabel('切向力 (kN/m)', fontsize=11)
ax1.set_title(f'微动迟滞回线 (相位={phase/np.pi:.2f}π)',
fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_xlim([-15, 15])
ax1.set_ylim([-300, 300])
# 右图：接触区滑移状态示意图
x = np.linspace(-fretting.b, fretting.b, 100)
# 绘制接触区
ax2.fill_between(x*1000, 0, 1, alpha=0.3, color='lightgray', label='接触区')
# 粘着区
c = fretting.c
ax2.fill_between([-c*1000, c*1000], 0, 1, alpha=0.5,
color='green', label='粘着区')
# 滑移区
ax2.fill_between([-fretting.b*1000, -c*1000], 0, 1, alpha=0.5,
color='red', label='滑移区')
ax2.fill_between([c*1000, fretting.b*1000], 0, 1, alpha=0.5,
color='red')
# 添加方向箭头表示运动
arrow_scale = np.cos(phase)
if abs(arrow_scale) > 0.1:
ax2.annotate('', xy=(arrow_scale*5, 0.5), xytext=(0, 0.5),
arrowprops=dict(arrowstyle='->', color='blue', lw=3))
ax2.set_xlim([-1.5*fretting.b*1000, 1.5*fretting.b*1000])
ax2.set_ylim([0, 1.2])
ax2.set_xlabel('位置 (mm)', fontsize=11)
ax2.set_title('接触区滑移状态', fontsize=12, fontweight='bold')
ax2.legend(loc='upper right')
ax2.set_yticks([])
# 保存帧
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=100, bbox_inches='tight')
buf.seek(0)
frames.append(Image.open(buf))
# 保存GIF
if frames:
frames[0].save(
'd:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/fretting_process_animation.gif',
save_all=True,
append_images=frames[1:],
duration=100,
loop=0
)
print("微动过程动画已保存")
plt.close()
def create_contact_stress_distribution(contact, fretting):
"""
创建接触应力分布图
"""
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 接触半宽内的位置
x = np.linspace(-fretting.b, fretting.b, 100)
X = x / fretting.b
# 1. 表面压力分布
ax1 = axes[0, 0]
p = fretting.p0 * np.sqrt(1 - X**2)
ax1.fill_between(x*1000, p/1e6, alpha=0.3, color='blue')
ax1.plot(x*1000, p/1e6, 'b-', linewidth=2)
ax1.set_xlabel('位置 x (mm)', fontsize=11)
ax1.set_ylabel('压力 (MPa)', fontsize=11)
ax1.set_title('表面接触压力分布', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
# 2. 次表面应力分布（沿深度方向）
ax2 = axes[0, 1]
depths = np.linspace(0, 3*fretting.b, 50)
# 简化公式：最大剪应力随深度变化
tau_max = fretting.p0 * 0.3 * (1 - depths/(3*fretting.b))
sigma_z = fretting.p0 * np.exp(-depths/fretting.b)
ax2.plot(tau_max/1e6, depths*1000, 'r-', linewidth=2, label='最大剪应力')
ax2.plot(sigma_z/1e6, depths*1000, 'b-', linewidth=2, label='法向应力')
ax2.axhline(y=0.786*fretting.b*1000, color='g', linestyle='--',
alpha=0.5, label='最大剪应力深度')
ax2.set_xlabel('应力 (MPa)', fontsize=11)
ax2.set_ylabel('深度 z (mm)', fontsize=11)
ax2.set_title('次表面应力分布（接触中心下方）', fontsize=12, fontweight='bold')
ax2.legend(loc='best')
ax2.grid(True, alpha=0.3)
ax2.invert_yaxis()
# 3. 切向牵引力分布（不同切向力水平）
ax3 = axes[1, 0]
Q_ratios = [0.2, 0.5, 0.8]
colors = ['blue', 'green', 'red']
for Q_ratio, color in zip(Q_ratios, colors):
Q = Q_ratio * fretting.mu * fretting.P
q = fretting.tangential_traction(x, Q)
ax3.plot(x*1000, q/1e6, color=color, linewidth=2,
label=f'Q/(μP)={Q_ratio}')
ax3.set_xlabel('位置 x (mm)', fontsize=11)
ax3.set_ylabel('切向牵引力 (MPa)', fontsize=11)
ax3.set_title('Mindlin切向牵引力分布', fontsize=12, fontweight='bold')
ax3.legend(loc='best')
ax3.grid(True, alpha=0.3)
# 4. 微动疲劳寿命预测曲线
ax4 = axes[1, 1]
fatigue = FrettingFatigue('steel')
# 不同滑移幅值下的寿命
slip_values = [5e-6, 10e-6, 20e-6, 30e-6]
pressures = np.linspace(200e6, 800e6, 20)
for slip, color in zip(slip_values, ['blue', 'green', 'orange', 'red']):
lives = []
for p in pressures:
life = fatigue.total_fretting_fatigue_life(p, slip)
lives.append(life['total_life'])
ax4.semilogy(pressures/1e6, lives, color=color, linewidth=2,
marker='o', markersize=4, label=f'δ={slip*1e6:.0f}μm')
ax4.set_xlabel('接触压力 (MPa)', fontsize=11)
ax4.set_ylabel('疲劳寿命 (cycles)', fontsize=11)
ax4.set_title('微动疲劳寿命预测曲线', fontsize=12, fontweight='bold')
ax4.legend(loc='best')
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('d:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/contact_stress_analysis.png',
dpi=150, bbox_inches='tight')
print("接触应力分析图已保存")
plt.close()
if __name__ == "__main__":
create_visualization()
print("\n所有仿真结果已生成完毕！")
print("\n生成的文件：")
print("1. fretting_fatigue_analysis.png - 微动疲劳分析综合图")
print("2. fretting_process_animation.gif - 微动过程动画")
print("3. contact_stress_analysis.png - 接触应力分析图")
# -*- coding: utf-8 -*-
"""
主题088：微动疲劳与磨损 - 实例二
磨损模型与表面演化模拟
本实例模拟磨损过程（Wear）的数学模型和表面形貌演化。
磨损是材料在接触表面相对运动时发生的材料损失现象，
包括粘着磨损、磨粒磨损、腐蚀磨损和疲劳磨损等类型。
应用场景：
- 机械密封件磨损
- 轴承磨损寿命预测
- 齿轮齿面磨损
- 人工关节磨损
作者：结构耐久性仿真教程编写组
日期：2026年3月
"""
import numpy as np
import matplotlib
matplotlib.use('Agg') # 使用无头模式
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from mpl_toolkits.mplot3d import Axes3D
import io
from PIL import Image
from scipy import ndimage
from scipy.signal import convolve2d
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("=" * 60)
print("主题088：微动疲劳与磨损 - 实例二")
print("磨损模型与表面演化模拟")
print("=" * 60)
class ArchardWearModel:
"""
Archard磨损模型
最经典的磨损模型，描述磨损体积与载荷、滑动距离的关系：
V = K * F * s / H
其中：
V - 磨损体积
K - 无量纲磨损系数
F - 法向载荷
s - 滑动距离
H - 材料硬度
"""
def __init__(self, K=1e-4, H=2e9):
"""
初始化Archard磨损模型
Parameters:
-----------
K : float
磨损系数（典型值：1e-4 ~ 1e-8）
H : float
材料硬度 (Pa)
"""
self.K = K
self.H = H
def wear_volume(self, F, s):
"""
计算磨损体积
Parameters:
-----------
F : float
法向载荷 (N)
s : float
滑动距离 (m)
Returns:
--------
V : float
磨损体积 (m³)
"""
V = self.K * F * s / self.H
return V
def wear_rate(self, F, v):
"""
计算磨损率（单位时间的磨损体积）
Parameters:
-----------
F : float
法向载荷 (N)
v : float
滑动速度 (m/s)
Returns:
--------
dV_dt : float
磨损率 (m³/s)
"""
dV_dt = self.K * F * v / self.H
return dV_dt
def wear_depth_increment(self, p, delta_s):
"""
计算局部磨损深度增量
Parameters:
-----------
p : float
局部接触压力 (Pa)
delta_s : float
滑动距离增量 (m)
Returns:
--------
delta_h : float
磨损深度增量 (m)
"""
# 局部磨损深度 = K * p * delta_s / H
delta_h = self.K * p * delta_s / self.H
return delta_h
class UsuiWearModel:
"""
Usui磨损模型（适用于切削磨损）
描述切削过程中的工具磨损
"""
def __init__(self, A=1e-10, B=1000, C=2.0):
"""
初始化Usui磨损模型
Parameters:
-----------
A, B, C : float
模型参数
"""
self.A = A
self.B = B
self.C = C
def tool_wear_rate(self, V, T):
"""
计算刀具磨损率
Parameters:
-----------
V : float
切削速度 (m/s)
T : float
温度 (°C)
Returns:
--------
wear_rate : float
磨损率
"""
wear_rate = self.A * V**self.C * np.exp(-self.B / (T + 273))
return wear_rate
class BitterWearModel:
"""
Bitter磨损模型（适用于冲蚀磨损）
描述颗粒冲击引起的材料损失
"""
def __init__(self, K1=1e-9, K2=1e-12, n1=2.0, n2=2.5):
"""
初始化Bitter磨损模型
Parameters:
-----------
K1, K2 : float
磨损常数
n1, n2 : float
速度指数
"""
self.K1 = K1
self.K2 = K2
self.n1 = n1
self.n2 = n2
def erosion_rate(self, V, angle):
"""
计算冲蚀率
Parameters:
-----------
V : float
颗粒速度 (m/s)
angle : float
冲击角度 (度)
Returns:
--------
E : float
冲蚀率 (kg/kg)
"""
# 将角度转换为弧度
theta = np.radians(angle)
# 变形磨损分量
E_deformation = self.K1 * V**self.n1 * np.sin(theta)**self.n1
# 切削磨损分量
E_cutting = self.K2 * V**self.n2 * np.sin(2*theta)**self.n2
# 总冲蚀率
E = E_deformation + E_cutting
return E
class SurfaceEvolution:
"""
表面形貌演化模拟
模拟磨损过程中表面粗糙度的变化和形貌演化
"""
def __init__(self, size_x=100, size_y=100, dx=1e-6, dy=1e-6):
"""
初始化表面网格
Parameters:
-----------
size_x, size_y : int
网格尺寸
dx, dy : float
网格间距 (m)
"""
self.size_x = size_x
self.size_y = size_y
self.dx = dx
self.dy = dy
# 初始化表面高度场
self.height = np.zeros((size_y, size_x))
# 生成初始粗糙表面
self.generate_initial_surface()
def generate_initial_surface(self, Ra=0.5e-6, correlation_length=10e-6):
"""
生成初始粗糙表面
Parameters:
-----------
Ra : float
算术平均粗糙度 (m)
correlation_length : float
相关长度 (m)
"""
# 使用高斯随机场生成粗糙表面
# 简化实现：使用滤波白噪声
# 生成白噪声
noise = np.random.randn(self.size_y, self.size_x)
# 高斯滤波模拟相关性
sigma = correlation_length / self.dx / 2.355 # FWHM到sigma的转换
self.height = ndimage.gaussian_filter(noise, sigma=sigma)
# 归一化到目标粗糙度
current_Ra = np.mean(np.abs(self.height))
if current_Ra > 0:
self.height = self.height * Ra / current_Ra
def apply_wear(self, wear_depth_field):
"""
应用磨损深度场
Parameters:
-----------
wear_depth_field : ndarray
磨损深度场 (m)
"""
self.height -= wear_depth_field
def calculate_roughness(self):
"""
计算表面粗糙度参数
Returns:
--------
roughness : dict
包含Ra, Rq, Rz等粗糙度参数
"""
# 算术平均粗糙度 Ra
Ra = np.mean(np.abs(self.height))
# 均方根粗糙度 Rq
Rq = np.sqrt(np.mean(self.height**2))
# 最大峰谷高度 Rz
Rz = np.max(self.height) - np.min(self.height)
# 偏度 Skewness
skewness = np.mean((self.height / Rq)**3) if Rq > 0 else 0
# 峰度 Kurtosis
kurtosis = np.mean((self.height / Rq)**4) if Rq > 0 else 0
return {
'Ra': Ra,
'Rq': Rq,
'Rz': Rz,
'skewness': skewness,
'kurtosis': kurtosis,
}
def calculate_contact_area(self, h_indentation):
"""
计算在给定压入深度下的接触面积
Parameters:
-----------
h_indentation : float
压入深度 (m)
Returns:
--------
A_contact : float
接触面积 (m²)
"""
# 找到高于参考平面的区域
contact_mask = self.height > (np.max(self.height) - h_indentation)
# 计算接触面积
n_contact = np.sum(contact_mask)
A_contact = n_contact * self.dx * self.dy
return A_contact
def get_surface_gradient(self):
"""
计算表面梯度
Returns:
--------
gradient : tuple
(dh/dx, dh/dy)
"""
grad_y, grad_x = np.gradient(self.height, self.dy, self.dx)
return grad_x, grad_y
class WearSimulation:
"""
磨损过程仿真
模拟磨损过程中表面形貌的时变演化
"""
def __init__(self, surface, wear_model, contact_model):
"""
初始化磨损仿真
Parameters:
-----------
surface : SurfaceEvolution
表面形貌对象
wear_model : ArchardWearModel
磨损模型
contact_model : dict
接触参数
"""
self.surface = surface
self.wear_model = wear_model
self.contact_model = contact_model
# 历史记录
self.history = {
'time': [],
'roughness': [],
'wear_volume': [],
}
def calculate_contact_pressure(self, F_normal):
"""
计算接触压力分布（简化模型）
Parameters:
-----------
F_normal : float
法向载荷 (N)
Returns:
--------
pressure : ndarray
接触压力分布 (Pa)
"""
# 简化：假设压力与局部高度成正比（最高的区域承受更大压力）
h_max = np.max(self.surface.height)
h_min = np.min(self.surface.height)
# 归一化高度
h_normalized = (self.surface.height - h_min) / (h_max - h_min + 1e-10)
# 压力分布（简化Hertz-like分布）
pressure = F_normal / (self.surface.size_x * self.surface.dx *
self.surface.size_y * self.surface.dy)
pressure_field = pressure * h_normalized**2
# 归一化
total_force = np.sum(pressure_field) * self.surface.dx * self.surface.dy
if total_force > 0:
pressure_field *= F_normal / total_force
return pressure_field
def step(self, dt, v_slide, F_normal):
"""
执行一个时间步长的磨损仿真
Parameters:
-----------
dt : float
时间步长 (s)
v_slide : float
滑动速度 (m/s)
F_normal : float
法向载荷 (N)
"""
# 计算滑动距离
delta_s = v_slide * dt
# 计算接触压力分布
pressure = self.calculate_contact_pressure(F_normal)
# 计算磨损深度增量
wear_depth = self.wear_model.wear_depth_increment(pressure, delta_s)
# 应用磨损
self.surface.apply_wear(wear_depth)
# 记录历史
self.history['time'].append(
self.history['time'][-1] + dt if self.history['time'] else dt
)
self.history['roughness'].append(self.surface.calculate_roughness())
self.history['wear_volume'].append(
self.history['wear_volume'][-1] + np.sum(wear_depth) *
self.surface.dx * self.surface.dy if self.history['wear_volume']
else np.sum(wear_depth) * self.surface.dx * self.surface.dy
)
def run(self, t_total, dt, v_slide, F_normal):
"""
运行完整磨损仿真
Parameters:
-----------
t_total : float
总仿真时间 (s)
dt : float
时间步长 (s)
v_slide : float
滑动速度 (m/s)
F_normal : float
法向载荷 (N)
"""
n_steps = int(t_total / dt)
for i in range(n_steps):
self.step(dt, v_slide, F_normal)
if (i + 1) % 100 == 0:
print(f"仿真进度: {(i+1)/n_steps*100:.1f}%")
class AdhesiveWearSimulation:
"""
粘着磨损仿真
模拟粘着磨损过程中材料转移和磨屑形成
"""
def __init__(self, material_props):
"""
初始化粘着磨损仿真
Parameters:
-----------
material_props : dict
材料属性
"""
self.props = material_props
def junction_growth(self, F_normal, A_contact, tau_yield):
"""
计算接合点生长
Parameters:
-----------
F_normal : float
法向载荷 (N)
A_contact : float
真实接触面积 (m²)
tau_yield : float
剪切屈服强度 (Pa)
Returns:
--------
A_junction : float
接合点面积 (m²)
"""
# 简化模型：接合点面积与法向载荷成正比
A_junction = F_normal / (3 * tau_yield)
return min(A_junction, A_contact)
def wear_particle_size(self, h_adhesion, d_junction):
"""
估算磨屑尺寸
Parameters:
-----------
h_adhesion : float
粘着层厚度 (m)
d_junction : float
接合点直径 (m)
Returns:
--------
d_particle : float
磨屑直径 (m)
"""
# 经验公式：磨屑直径约为接合点直径的1/3到1/2
d_particle = 0.4 * d_junction
return d_particle
def wear_rate_adhesive(self, F_normal, v_slide, H_material, K_adhesive=1e-4):
"""
计算粘着磨损率
Parameters:
-----------
F_normal : float
法向载荷 (N)
v_slide : float
滑动速度 (m/s)
H_material : float
材料硬度 (Pa)
K_adhesive : float
粘着磨损系数
Returns:
--------
Q : float
磨损率 (m³/s)
"""
Q = K_adhesive * F_normal * v_slide / H_material
return Q
class AbrasiveWearSimulation:
"""
磨粒磨损仿真
模拟硬颗粒对软材料的切削和犁沟作用
"""
def __init__(self, abrasive_props):
"""
初始化磨粒磨损仿真
Parameters:
-----------
abrasive_props : dict
磨粒属性
"""
self.abrasive = abrasive_props
def cutting_volume(self, F_normal, theta, d_abrasive):
"""
计算单个磨粒的切削体积
Parameters:
-----------
F_normal : float
法向力 (N)
theta : float
磨粒锥角 (度)
d_abrasive : float
磨粒直径 (m)
Returns:
--------
V_cut : float
切削体积 (m³)
"""
# 简化模型：假设磨粒为圆锥形
theta_rad = np.radians(theta)
# 压入深度
h_indent = (2 * F_normal / (np.pi * d_abrasive *
self.abrasive['hardness'] * np.tan(theta_rad)))**0.5
# 切削体积（圆锥体积）
r_contact = h_indent * np.tan(theta_rad)
V_cut = (1/3) * np.pi * r_contact**2 * h_indent
return V_cut
def wear_rate_abrasive(self, C_abrasive, v_slide, F_normal, H_soft):
"""
计算磨粒磨损率（Rabinowicz模型）
Parameters:
-----------
C_abrasive : float
磨粒浓度
v_slide : float
滑动速度 (m/s)
F_normal : float
法向载荷 (N)
H_soft : float
软材料硬度 (Pa)
Returns:
--------
Q : float
磨损率 (m³/s)
"""
# Rabinowicz磨粒磨损方程
cot_theta = 1 / np.tan(np.radians(60)) # 假设磨粒锥角60度
Q = C_abrasive * F_normal * v_slide * cot_theta / (3 * H_soft)
return Q
def create_wear_visualization():
"""
创建磨损分析的可视化
"""
print("\n正在生成磨损分析可视化...")
# 1. Archard磨损模型分析
fig1, axes1 = plt.subplots(2, 2, figsize=(14, 10))
# 不同材料的磨损系数
materials = {
'钢/钢': {'K': 1e-4, 'H': 2e9, 'color': 'blue'},
'钢/青铜': {'K': 5e-5, 'H': 1.5e9, 'color': 'green'},
'钢/聚四氟乙烯': {'K': 1e-6, 'H': 0.05e9, 'color': 'red'},
'陶瓷/陶瓷': {'K': 1e-8, 'H': 20e9, 'color': 'purple'},
}
# 1.1 磨损体积 vs 滑动距离
ax1 = axes1[0, 0]
F = 100 # N
s_range = np.linspace(0, 1000, 100) # m
for name, props in materials.items():
wear_model = ArchardWearModel(K=props['K'], H=props['H'])
V_wear = [wear_model.wear_volume(F, s) * 1e9 for s in s_range] # mm³
ax1.plot(s_range, V_wear, color=props['color'], linewidth=2,
label=name, marker='o', markersize=3, markevery=10)
ax1.set_xlabel('滑动距离 (m)', fontsize=11)
ax1.set_ylabel('磨损体积 (mm³)', fontsize=11)
ax1.set_title('Archard磨损模型：磨损体积 vs 滑动距离', fontsize=12, fontweight='bold')
ax1.legend(loc='best')
ax1.grid(True, alpha=0.3)
# 1.2 磨损率 vs 载荷
ax2 = axes1[0, 1]
F_range = np.linspace(10, 500, 50) # N
v = 1.0 # m/s
for name, props in materials.items():
wear_model = ArchardWearModel(K=props['K'], H=props['H'])
wear_rates = [wear_model.wear_rate(F, v) * 1e12 for F in F_range] # mm³/s
ax2.plot(F_range, wear_rates, color=props['color'], linewidth=2, label=name)
ax2.set_xlabel('法向载荷 (N)', fontsize=11)
ax2.set_ylabel('磨损率 (mm³/s)', fontsize=11)
ax2.set_title('磨损率 vs 法向载荷', fontsize=12, fontweight='bold')
ax2.legend(loc='best')
ax2.grid(True, alpha=0.3)
# 1.3 冲蚀磨损角度效应（Bitter模型）
ax3 = axes1[1, 0]
bitter_model = BitterWearModel()
angles = np.linspace(0, 90, 100)
velocities = [50, 100, 150] # m/s
colors = ['blue', 'green', 'red']
for V, color in zip(velocities, colors):
erosion_rates = [bitter_model.erosion_rate(V, angle) for angle in angles]
ax3.plot(angles, erosion_rates, color=color, linewidth=2,
label=f'V={V}m/s')
ax3.set_xlabel('冲击角度 (度)', fontsize=11)
ax3.set_ylabel('冲蚀率', fontsize=11)
ax3.set_title('Bitter冲蚀磨损模型：角度效应', fontsize=12, fontweight='bold')
ax3.legend(loc='best')
ax3.grid(True, alpha=0.3)
# 1.4 磨损系数对比
ax4 = axes1[1, 1]
material_names = list(materials.keys())
K_values = [materials[name]['K'] for name in material_names]
bars = ax4.barh(material_names, K_values, color=[materials[name]['color']
for name in material_names])
ax4.set_xlabel('磨损系数 K', fontsize=11)
ax4.set_title('不同材料配对的磨损系数', fontsize=12, fontweight='bold')
ax4.set_xscale('log')
ax4.grid(True, alpha=0.3, axis='x')
# 添加数值标签
for i, (bar, K) in enumerate(zip(bars, K_values)):
ax4.text(K * 1.5, i, f'{K:.0e}', va='center', fontsize=9)
plt.tight_layout()
plt.savefig('d:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/wear_models_analysis.png',
dpi=150, bbox_inches='tight')
print("磨损模型分析图已保存")
plt.close()
# 2. 表面形貌演化仿真
print("\n正在进行表面形貌演化仿真...")
# 创建初始表面
surface = SurfaceEvolution(size_x=100, size_y=100, dx=1e-6, dy=1e-6)
# 创建磨损模型
wear_model = ArchardWearModel(K=1e-4, H=2e9)
# 创建磨损仿真
wear_sim = WearSimulation(surface, wear_model, {})
# 运行仿真
t_total = 100 # s
dt = 0.1 # s
v_slide = 0.1 # m/s
F_normal = 50 # N
wear_sim.run(t_total, dt, v_slide, F_normal)
# 3. 表面形貌演化可视化
fig2, axes2 = plt.subplots(2, 3, figsize=(18, 10))
# 获取不同时间点的表面
time_points = [0, 25, 50, 75, 99]
# 3.1-3.3 不同时刻的表面形貌
for idx, t_idx in enumerate(time_points[:3]):
ax = axes2[0, idx]
# 重新生成初始表面并演化到指定时刻
surf = SurfaceEvolution(size_x=100, size_y=100, dx=1e-6, dy=1e-6)
sim = WearSimulation(surf, wear_model, {})
for _ in range(t_idx + 1):
sim.step(dt, v_slide, F_normal)
im = ax.imshow(surf.height * 1e6, cmap='terrain',
extent=[0, 100, 0, 100])
ax.set_title(f'表面形貌 (t={t_idx*dt:.1f}s)', fontsize=12, fontweight='bold')
ax.set_xlabel('X (μm)', fontsize=10)
ax.set_ylabel('Y (μm)', fontsize=10)
plt.colorbar(im, ax=ax, label='高度 (μm)')
# 3.4 粗糙度演化
ax4 = axes2[1, 0]
times = wear_sim.history['time']
Ra_values = [r['Ra'] * 1e6 for r in wear_sim.history['roughness']]
Rq_values = [r['Rq'] * 1e6 for r in wear_sim.history['roughness']]
ax4.plot(times, Ra_values, 'b-', linewidth=2, label='Ra (算术平均)')
ax4.plot(times, Rq_values, 'r--', linewidth=2, label='Rq (均方根)')
ax4.set_xlabel('时间 (s)', fontsize=11)
ax4.set_ylabel('粗糙度 (μm)', fontsize=11)
ax4.set_title('表面粗糙度演化', fontsize=12, fontweight='bold')
ax4.legend(loc='best')
ax4.grid(True, alpha=0.3)
# 3.5 磨损体积演化
ax5 = axes2[1, 1]
wear_volumes = np.array(wear_sim.history['wear_volume']) * 1e12 # mm³
ax5.plot(times, wear_volumes, 'g-', linewidth=2)
ax5.set_xlabel('时间 (s)', fontsize=11)
ax5.set_ylabel('累积磨损体积 (mm³)', fontsize=11)
ax5.set_title('累积磨损体积演化', fontsize=12, fontweight='bold')
ax5.grid(True, alpha=0.3)
# 3.6 表面高度分布直方图
ax6 = axes2[1, 2]
final_heights = surface.height.flatten() * 1e6 # μm
ax6.hist(final_heights, bins=30, color='steelblue', alpha=0.7, edgecolor='black')
ax6.axvline(x=np.mean(final_heights), color='r', linestyle='--',
linewidth=2, label=f'均值={np.mean(final_heights):.3f}μm')
ax6.set_xlabel('高度 (μm)', fontsize=11)
ax6.set_ylabel('频数', fontsize=11)
ax6.set_title('最终表面高度分布', fontsize=12, fontweight='bold')
ax6.legend(loc='best')
ax6.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('d:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/surface_evolution_analysis.png',
dpi=150, bbox_inches='tight')
print("表面形貌演化分析图已保存")
plt.close()
# 4. 3D表面形貌可视化
print("\n正在生成3D表面形貌图...")
fig3 = plt.figure(figsize=(16, 6))
# 初始表面
ax1 = fig3.add_subplot(121, projection='3d')
surf_initial = SurfaceEvolution(size_x=100, size_y=100, dx=1e-6, dy=1e-6)
X = np.linspace(0, 100, 100)
Y = np.linspace(0, 100, 100)
X, Y = np.meshgrid(X, Y)
ax1.plot_surface(X, Y, surf_initial.height * 1e6, cmap='viridis',
alpha=0.9, edgecolor='none')
ax1.set_xlabel('X (μm)', fontsize=10)
ax1.set_ylabel('Y (μm)', fontsize=10)
ax1.set_zlabel('高度 (μm)', fontsize=10)
ax1.set_title('初始表面形貌', fontsize=12, fontweight='bold')
# 磨损后表面
ax2 = fig3.add_subplot(122, projection='3d')
ax2.plot_surface(X, Y, surface.height * 1e6, cmap='viridis',
alpha=0.9, edgecolor='none')
ax2.set_xlabel('X (μm)', fontsize=10)
ax2.set_ylabel('Y (μm)', fontsize=10)
ax2.set_zlabel('高度 (μm)', fontsize=10)
ax2.set_title(f'磨损后表面形貌 (t={t_total}s)', fontsize=12, fontweight='bold')
plt.tight_layout()
plt.savefig('d:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/surface_3d_visualization.png',
dpi=150, bbox_inches='tight')
print("3D表面形貌图已保存")
plt.close()
# 5. 粘着磨损和磨粒磨损对比
fig4, axes4 = plt.subplots(2, 2, figsize=(14, 10))
adhesive_sim = AdhesiveWearSimulation({
'yield_strength': 400e6,
'hardness': 2e9,
})
abrasive_sim = AbrasiveWearSimulation({
'hardness': 20e9, # 磨粒硬度
'size': 100e-6, # 磨粒尺寸
})
# 5.1 粘着磨损率 vs 载荷
ax1 = axes4[0, 0]
F_range = np.linspace(10, 200, 50)
v = 1.0
H_material = 2e9
K_adhesive_values = [1e-3, 1e-4, 1e-5]
colors = ['red', 'blue', 'green']
for K, color in zip(K_adhesive_values, colors):
wear_rates = [adhesive_sim.wear_rate_adhesive(F, v, H_material, K) * 1e12
for F in F_range]
ax1.plot(F_range, wear_rates, color=color, linewidth=2,
label=f'K={K:.0e}')
ax1.set_xlabel('法向载荷 (N)', fontsize=11)
ax1.set_ylabel('磨损率 (mm³/s)', fontsize=11)
ax1.set_title('粘着磨损率 vs 载荷', fontsize=12, fontweight='bold')
ax1.legend(loc='best')
ax1.grid(True, alpha=0.3)
# 5.2 磨粒磨损率 vs 磨粒浓度
ax2 = axes4[0, 1]
C_range = np.linspace(0.01, 0.2, 50)
F = 100
H_soft = 1e9
wear_rates_abrasive = [abrasive_sim.wear_rate_abrasive(C, v, F, H_soft) * 1e12
for C in C_range]
ax2.plot(C_range * 100, wear_rates_abrasive, 'b-', linewidth=2)
ax2.set_xlabel('磨粒浓度 (%)', fontsize=11)
ax2.set_ylabel('磨损率 (mm³/s)', fontsize=11)
ax2.set_title('磨粒磨损率 vs 磨粒浓度', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
# 5.3 磨屑尺寸分布
ax3 = axes4[1, 0]
# 模拟磨屑尺寸分布（对数正态分布）
d_particles = np.random.lognormal(mean=-20, sigma=1.5, size=1000) # m
d_particles_um = d_particles * 1e6 # μm
ax3.hist(d_particles_um, bins=50, color='coral', alpha=0.7,
edgecolor='black')
ax3.set_xlabel('磨屑尺寸 (μm)', fontsize=11)
ax3.set_ylabel('频数', fontsize=11)
ax3.set_title('磨屑尺寸分布', fontsize=12, fontweight='bold')
ax3.set_xlim([0, 50])
ax3.grid(True, alpha=0.3, axis='y')
# 5.4 磨损机制图
ax4 = axes4[1, 1]
# 创建磨损机制图（载荷-速度空间）
F_grid = np.linspace(1, 500, 50)
v_grid = np.linspace(0.01, 10, 50)
F_grid, v_grid = np.meshgrid(F_grid, v_grid)
# 简化的磨损机制判据
# 粘着磨损主导区：低载荷、低速度
# 磨粒磨损主导区：中等载荷、中等速度
# 氧化磨损主导区：高速度
# 熔融磨损主导区：高载荷、高速度
mechanism = np.zeros_like(F_grid)
PV = F_grid * v_grid
mechanism[(PV < 50)] = 1 # 粘着磨损
mechanism[(PV >= 50) & (PV < 200)] = 2 # 磨粒磨损
mechanism[(PV >= 200) & (PV < 500)] = 3 # 氧化磨损
mechanism[(PV >= 500)] = 4 # 熔融磨损
cmap = plt.cm.colors.ListedColormap(['lightblue', 'yellow', 'lightgreen', 'salmon'])
im = ax4.contourf(F_grid, v_grid, mechanism, levels=[0.5, 1.5, 2.5, 3.5, 4.5],
colors=['lightblue', 'yellow', 'lightgreen', 'salmon'])
ax4.set_xlabel('法向载荷 (N)', fontsize=11)
ax4.set_ylabel('滑动速度 (m/s)', fontsize=11)
ax4.set_title('磨损机制图 (PV图)', fontsize=12, fontweight='bold')
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor='lightblue', label='粘着磨损'),
Patch(facecolor='yellow', label='磨粒磨损'),
Patch(facecolor='lightgreen', label='氧化磨损'),
Patch(facecolor='salmon', label='熔融磨损')
]
ax4.legend(handles=legend_elements, loc='upper right')
plt.tight_layout()
plt.savefig('d:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/wear_mechanisms_comparison.png',
dpi=150, bbox_inches='tight')
print("磨损机制对比图已保存")
plt.close()
# 6. 创建表面演化动画
print("\n正在生成表面演化动画...")
create_surface_evolution_animation()
print("\n所有磨损仿真结果已生成完毕！")
def create_surface_evolution_animation():
"""
创建表面形貌演化动画
"""
# 创建初始表面
surface = SurfaceEvolution(size_x=80, size_y=80, dx=1e-6, dy=1e-6)
wear_model = ArchardWearModel(K=1e-4, H=2e9)
wear_sim = WearSimulation(surface, wear_model, {})
# 仿真参数
dt = 0.5
v_slide = 0.1
F_normal = 50
# 创建图形
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
frames = []
n_frames = 30
for i in range(n_frames):
ax1.clear()
ax2.clear()
# 执行磨损步骤
for _ in range(5):
wear_sim.step(dt, v_slide, F_normal)
# 左图：表面形貌
im1 = ax1.imshow(surface.height * 1e6, cmap='terrain',
extent=[0, 80, 0, 80], vmin=-2, vmax=2)
ax1.set_title(f'表面形貌演化 (t={(i+1)*5*dt:.1f}s)',
fontsize=12, fontweight='bold')
ax1.set_xlabel('X (μm)', fontsize=10)
ax1.set_ylabel('Y (μm)', fontsize=10)
plt.colorbar(im1, ax=ax1, label='高度 (μm)')
# 右图：粗糙度和磨损体积历史
times = wear_sim.history['time']
Ra_values = [r['Ra'] * 1e6 for r in wear_sim.history['roughness']]
wear_volumes = np.array(wear_sim.history['wear_volume']) * 1e12
ax2_twin = ax2.twinx()
line1, = ax2.plot(times, Ra_values, 'b-', linewidth=2, label='Ra')
line2, = ax2_twin.plot(times, wear_volumes, 'r--', linewidth=2, label='磨损体积')
ax2.set_xlabel('时间 (s)', fontsize=11)
ax2.set_ylabel('粗糙度 Ra (μm)', fontsize=11, color='b')
ax2_twin.set_ylabel('磨损体积 (mm³)', fontsize=11, color='r')
ax2.tick_params(axis='y', labelcolor='b')
ax2_twin.tick_params(axis='y', labelcolor='r')
ax2.set_title('粗糙度与磨损体积演化', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
lines = [line1, line2]
labels = [l.get_label() for l in lines]
ax2.legend(lines, labels, loc='center right')
# 保存帧
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=100, bbox_inches='tight')
buf.seek(0)
frames.append(Image.open(buf))
# 保存GIF
if frames:
frames[0].save(
'd:/文档/500仿真领域/工程仿真/结构热应力仿真/主题088_微动疲劳与磨损/surface_evolution_animation.gif',
save_all=True,
append_images=frames[1:],
duration=200,
loop=0
)
print("表面演化动画已保存")
plt.close()
if __name__ == "__main__":
create_wear_visualization()
print("\n生成的文件：")
print("1. wear_models_analysis.png - 磨损模型分析图")
print("2. surface_evolution_analysis.png - 表面形貌演化分析图")
print("3. surface_3d_visualization.png - 3D表面形貌图")
print("4. wear_mechanisms_comparison.png - 磨损机制对比图")
print("5. surface_evolution_animation.gif - 表面演化动画")
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
2
0- 0
已为社区贡献265条内容
所有评论(0)