U75V钢轨微观晶粒结构对滚动接触疲劳损伤行为影响

Effect of the Microscopic Grain Structure of U75V Rail on Rolling Contact Fatigue Damage Behavior

摘要

钢轨材料晶粒的拓扑结构、尺寸和取向等微观结构特性，对钢轨早期滚动接触疲劳（Rolling Contact Fatigue，RCF）裂纹演化影响显著。为分析晶粒结构对RCF裂纹的影响，基于晶体塑性理论和维诺图（Voronoi）原理建立反映U75V钢轨材料晶粒几何结构、取向等特征的微观模型；通过修正的内聚力单元构建材料微观结构与疲劳损伤之间的联系，实现在微观尺度模拟RCF裂纹演化，建立钢轨微观RCF裂纹演化模型；基于该模型，研究循环荷载作用下钢轨微观组织结构对表面RCF裂纹萌生位置和扩展路径的影响。结果表明：循环荷载作用下，RCF裂纹大多萌生于晶界处，微裂纹萌生的位置和寿命具有随机性，但萌生位置在特定区域内呈现集中趋势；RCF裂纹在不同阶段的扩展速率不同，随扩展长度的增加整体呈早期慢后期逐渐变快的趋势，可分为2个阶段；Ⅰ阶段为起步阶段，微裂纹扩展长度在该阶段为几个晶粒尺寸，且该阶段裂纹的萌生扩展行为主要受晶粒取向和晶界间相互作用支配，裂纹扩展速率在0.002 5 μm · r^-1左右，远低于裂纹平均扩展速率0.010 3 μm · r^-1；当萌生位置间隔较近的微裂纹在扩展过程中相互吸引聚合，形成1条新裂纹时，进入Ⅱ阶段；新裂纹的扩展行为主要受外部荷载的支配，裂纹扩展速率明显增加，最快可达0.018 7 μm · r^-1；模拟结果与双轮试验观测的裂纹形态吻合良好。

Abstract

Microstructural properties such as grain topology, size, and orientation of rail materials have a significant effect on early rolling contact fatigue (RCF) crack evolution. To analyze the influence of grain structure on RCF cracks, a microscopic model reflecting the geometric structure and orientation of grains in U75V rail materials was established based on crystal plasticity theory and the Voronoi principle. The model was further enhanced by modified cohesive elements to forge a link between the microstructure of the material and fatigue damage, facilitating the simulation of RCF crack evolution at the microscopic scale and establishing a microscopic RCF crack evolution model for rails. Based on this model, the effects of rail microstructure on the initiation locations and propagation paths of surface RCF cracks under cyclic loading were investigated. The results show that under cyclic loading, RCF cracks mostly initiate at grain boundaries. While the microcrack initiation locations and lifetimes exhibit randomness, the initiation locations show a tendency to be concentrated in specific regions. The growth rate of RCF cracks varies across different stages, with the increase of growth length, the overall trend of early slow late gradually become faster, which can be divided into stages I and II. Stage I is the initiation stage, where microcracks propagate over a length equivalent to several grain sizes. The crack growth behavior during this stage is mainly governed by the grain orientation and the interaction between grain boundaries, and the crack growth rate is around 0.002 5 μm · r^-1, which is lower than the crack growth expansion rate of 0.010 3 μm · r^-1. When the microcracks that are closely spaced at the initiation position are attracted to each other and polymerize during the growth process to form a new crack, they enter into stage II. The growth behavior of new cracks is mainly governed by external loading, and the crack growth rate increases significantly to 0.018 7 μm · r^-1. The simulated crack morphology matches the structure observed in the twin-disc test.

Graphical abstract

关键词

U75V钢轨 / 滚动接触疲劳 / 微观结构 / 晶体塑性 / 裂纹扩展 / 微观模型

Key words

U75V Rail / Rolling contact fatigue (RCF) / Microstructure / Crystal plasticity / Crack propagation / Microscopic model

钢轨表面RCF和磨耗一直是影响铁路钢轨使用寿命的主要问题^［1-3］。随着U75V和U78CrV等硬质钢轨在我国重载铁路的普遍使用，钢轨RCF裂纹呈现萌生早、扩展快等特征^［4-5］。在轮轨接触斑重复应力作用下，RCF裂纹可能在接触表面或亚表面萌生，并快速扩展，最终导致钢轨表面剥离掉块，影响行车安全。既有的试验研究证实了微观晶粒结构对金属材料疲劳行为的影响，RCF裂纹高度依赖微观结构^［6-9］。

晶体塑性理论是在微观结构尺度上研究材料变形的常用方法^［10］，该方法以位错滑移为主要变形机制，可以在微观尺度上准确描述材料的弹塑性变形行为。然而，当前学者们主要关注对材料本身变形行为的描述，提出了一系列针对研究材料的本构模型，对材料在不同荷载下的变形行为进行了大量研究^［11-13］，但对晶体塑性理论与RCF损伤相结合的研究较少。Pandkar等^［14］研究了金属材料微观结构对RCF损伤积累在空间和时间速率方面的影响；Vijay等^［15-16］提出了一种采用损伤力学-内聚单元耦合的方法（Damage Mechanics-Cohesive Element Method，DM-CEM）模拟轴承载荷作用下多晶材料的RCF裂纹萌生和扩展，进而研究轴承的失效特征；Ghodrati等^［17-18］建立了一种基于晶体塑性有限元研究RCF裂纹的仿真框架，分析了牵引系数、温度变化、最大接触压力和晶粒尺寸等不同参数对RCF裂纹的影响；Wang等^［19］在微观层面上评估了齿轮RCF性能，并研究了夹杂物对损伤积累的影响。这些研究主要集中于滚动轴承和齿轮领域，通过考虑材料的晶粒尺寸和取向的随机分布，建立材料晶粒模型，然后使用不同的疲劳模型分析晶粒模型中的疲劳损伤演化。然而金相观测表明^［9］，钢轨材料在经过热处理后，材料的晶粒取向和尺寸等这些参数不再满足随机分布，而这些参数又会影响材料的局部应力场，对RCF裂纹的早期演化行为起着关键作用，仍需一种方法能综合反映真实材料的晶粒尺寸、取向以及晶界等因素对RCF裂纹演化的影响。

本文首先基于晶体塑性理论，模拟钢轨材料微观晶粒的变形行为；然后对不同循环荷载次数下的钢轨材料晶粒演变特征进行统计分析，应用Voronoi原理提出了一种基于电子背散衍射（Electron Back Scatter Diffraction，EBSD）观测结果的材料真实微观结构建模方法；最后通过在晶粒内部和晶界处引入零厚度内聚力单元模拟材料微观结构对裂纹的影响，并结合多轴疲劳准则（Jiang-Sehitoglu，J-S）模型模拟内聚力单元的疲劳损伤累计，实现RCF裂纹萌生扩展与材料微观结构的相互作用。研究结果为理解钢轨表面RCF的机制，预测RCF裂纹演化提供参考。

1 微观晶粒描述理论及几何模型构建

1.1　BCC金属晶体塑性理论

采用Hill等^［21-22］、Peirce等^［23］、Asaro等^［24］及Havner等^［25］建立的率相关塑性本构模型，研究滚动接触下钢轨BCC金属材料的变形响应。总变形梯度\mathit{F}可以描述为

\mathit{F}={\mathit{F}}_{\mathrm{e}}{\mathit{F}}_{\mathrm{p}}（1）

式中：{\mathit{F}}_{\mathrm{e}}为弹性变形梯度；{\mathit{F}}_{\mathrm{p}}为位错滑移引起的塑性变形梯度。

晶体塑性理论认为塑性变形是由滑移面上的位错运动产生的。U75V钢轨材料为BCC金属，共有12个主滑移系，对12个滑移系做分解，定义{\mathit{m}}_{{\alpha }_{\eta }}和{\mathit{n}}_{{\alpha }_{\eta }}分别为第\eta个滑移系{\alpha }_{\eta }的滑移方向向量和法向向量，它们满足正交性。则塑性变形梯度张量{\mathit{L}}_{\mathrm{p}}为所有滑移系上滑移剪切率的叠加，即

{\mathit{L}}_{\mathrm{p}}={\displaystyle \sum _{\eta =\mathrm{1}}^{\mathrm{12}}}{\dot{\gamma }}_{{\alpha }_{\eta }}({\mathit{m}}_{{\alpha }_{\eta }}{\mathit{n}}_{{\alpha }_{\eta }})（2）

式中：{\dot{\gamma }}_{{\alpha }_{\eta }}为滑移系{\alpha }_{\eta }引起的滑移剪切率。

根据Schmid定理，若晶粒结构不改变，任意滑移系{\alpha }_{\eta }的剪切滑移率{\dot{\gamma }}_{{\alpha }_{\eta }}取决于该滑移系上的Schmid应力{\mathit{\tau }}_{{\alpha }_{\eta }}，则{\dot{\gamma }}_{{\alpha }_{\eta }}可描述为

{\dot{\gamma }}_{{\alpha }_{\eta }}={〈\frac{{\mathit{\tau }}_{{\alpha }_{\eta }}-{\mathit{\lambda }}_{{\alpha }_{\eta }}}{g_{{\alpha }_{\eta }}}〉}^n\mathrm{s}\mathrm{g}\mathrm{n}({\mathit{\tau }}_{{\alpha }_{\eta }}-{\mathit{\lambda }}_{{\alpha }_{\eta }})（3）

式中：g_{{\alpha }_{\eta }}为滑移系{\alpha }_{\eta }的变形阻力，表示晶体硬化状态；{\mathit{\lambda }}_{{\alpha }_{\eta }}为滑移系{\alpha }_{\eta }的背应力，表示材料的包辛格效应；< >为Macauley运算符；n为滑移系{\alpha }_{\eta }的应变速率敏感指数，主要控制材料的黏性，当n趋于无穷大时，即为率无关材料。

Hill^［26］的研究指出某一滑移系的硬化会受到其余滑移系的影响，则某一滑移系{\alpha }_{\eta }的变形阻力率{\dot{g}}_{{\alpha }_{\eta }}表示为

{\dot{g}}_{{\alpha }_{\eta }}={\displaystyle \sum _{\zeta }^{\mathrm{12}}}{h}_{{\alpha }_{\eta }{\beta }_{\zeta }}{\dot{\gamma }}_{{\beta }_{\zeta }}（4）

其中，

h_{{\alpha }_{\eta }{\beta }_{\zeta }}=q{h}_{{\alpha }_{\eta }{\alpha }_{\eta }}

式中：{\beta }_{\zeta }为与{\alpha }_{\eta }不同的滑移系，当\eta =\zeta时，即表示相同滑移系；h_{{\alpha }_{\eta }{\beta }_{\zeta }}为混合硬化参数，当\eta =\zeta时称为自硬化参数（也可写为h_{{\alpha }_{\eta }{\alpha }_{\eta }}），否则为潜硬化参数；q为潜硬化参数与自硬化参数的比值。

目前常用幂函数^［23］描述材料的自硬化参数，即

h_{{\alpha }_{\eta }{\alpha }_{\eta }}=h_{\mathrm{0}}\mathrm{s}\mathrm{e}\mathrm{c}{\mathrm{h}}^{\mathrm{2}}\left|\frac{h_{\mathrm{0}}}{{\tau }_{\mathrm{s}}-{\tau }_{\mathrm{0}}}\right|（5）

式中：h_{\mathrm{0}}为初始硬化模量；{\tau }_{\mathrm{s}}为滑移系初始应力；{\tau }_{\mathrm{0}}为出现大变形时的饱和应力。

基于上述晶体塑性本构模型，借助有限元软件Abaqus的材料子程序（User Material Subroutine，UMAT），并基于Fortran语言将晶体塑性本构方程离散化编写为用户子程序以实现模拟^［27-29］。

1.2　晶粒微观损伤演化准则

通过在晶粒内部和晶界嵌入零厚度内聚力单元模拟材料损伤演化，内聚力单元的静态损伤采用双线性内聚力模型，并基于Jiang和Sehitoglu提出的多轴疲劳准则（J-S模型^［30］）修正内聚力单元的疲劳损伤累计。

1.2.1　混合加载模式下牵引-分离准则

采用基于等效位移的双线性内聚力模型进行静态损伤仿真，内聚力单元牵引-分离定律的本构方程^［31］这里不再赘述。

在断裂力学中将裂纹扩展分为Ⅰ型、Ⅱ型和Ⅲ型3种模式，Ⅰ型主要承受拉应力，Ⅱ型和Ⅲ型主要承受剪切应力。Ⅰ型、Ⅱ型和Ⅲ型这3种单一加载模式下内聚力单元应力与位移之间的关系可表示为

（6）

其中，

D_i=\frac{{\delta }_{\mathrm{f}i}({\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}i}-{\delta }_{\mathrm{0}i})}{{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}i}({\delta }_{\mathrm{f}i}-{\delta }_{\mathrm{0}i})}

式中：i表示方向，其中n表示法向，s和t表示2个切向；{\mathrm{\sigma }}_i为对应位移方向的应力；k_i为对应位移方向的单元刚度；{\delta }_i为加载位移；D_i为不同加载模式对应的静态损伤参数；{\delta }_{\mathrm{f}i}为断裂时对应方向的位移；{\delta }_{\mathrm{0}i}为初始损伤时对应方向的位移；{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}i}为对应方向的最大加载位移。

至此可以得到3种单一加载模式下内聚力单元的牵引-失效准则。但真实轮轨接触下，钢轨RCF裂纹通常是多种加载模式共同作用的结果，内聚力单元的损伤应该由3个方向的应力共同作用^［32］，因此，采用二次名义应力准则作为初始损伤判定。

一般认为2个切向方向的内聚力刚度和断裂能相等，因此，混合加载模式下的混合位移比\Gamma可以表示为

\Gamma =\frac{{\delta }_{\mathrm{L}}}{{\delta }_{\mathrm{n}}}{\delta }_{\mathrm{n}}>\mathrm{0}（7）

式中：{\delta }_{\mathrm{L}}为2个切向位移的模。

则混合加载模式下损伤初始对应的位移{\delta }_{\mathrm{0}\mathrm{m}}可以通过式（8）得到。

{\delta }_{\mathrm{0}\mathrm{m}}=\left\{\begin{array}{l}{\delta }_{\mathrm{0}\mathrm{n}}{\delta }_{\mathrm{0}\mathrm{s}}\sqrt[]{\frac{\mathrm{1}+{\Gamma }^{\mathrm{2}}}{{\delta }_{\mathrm{0}\mathrm{s}}^{\mathrm{2}}+(\Gamma {\delta }_{\mathrm{0}\mathrm{n}}{)}^{\mathrm{2}}}}{\delta }_{\mathrm{n}}>\mathrm{0}\\ {\delta }_{\mathrm{0}\mathrm{s}}{\delta }_{\mathrm{n}}\le \mathrm{0}\end{array}\right.（8）

最后采用BK准则^［33］，就可以得到混合加载模式下临界断裂能G_{\mathrm{c}\mathrm{m}}为

G_{\mathrm{c}\mathrm{m}}=G_{\mathrm{c}\mathrm{n}}+(G_{\mathrm{c}\mathrm{s}}-G_{\mathrm{c}\mathrm{n}}){\left(\frac{G_{\mathrm{s}}+G_{\mathrm{t}}}{G_{\mathrm{n}}+G_{\mathrm{s}}+G_{\mathrm{t}}}\right)}^{\eta }（9）

式中：G_{\mathrm{c}\mathrm{n}}和G_{\mathrm{c}\mathrm{s}}分别为法向和切向加载模式下临界断裂能；G_{\mathrm{n}}，G_{\mathrm{s}}和G_{\mathrm{t}}分别为内聚力单元的法向、2个剪切方向上所做的功。

1.2.2　疲劳损伤累积模型

前节考虑了混合加载模式下内聚力单元的静态损伤，但钢轨RCF裂纹多为在滚动循环荷载下的疲劳损伤累积造成的^［34］。因此，将J-S模型^［30］结合到内聚力单元中，并假设疲劳损伤累积仅发生在内聚力单元中。J-S模型中疲劳参数F_{\mathrm{P}}可表示为

F_{\mathrm{P}}=〈{\mathrm{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{x}}〉\frac{∆\epsilon }{\mathrm{2}}+J∆\tau ∆r（10）

其中，

〈{\mathrm{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{x}}〉=\mathrm{0.5}\left(\left|{\mathrm{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right|+{\mathrm{\sigma }}_{max}\right)

式中：{\mathrm{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{x}}为裂纹所在平面内的最大拉伸应力；∆\epsilon为垂直于裂纹表面的拉伸应变；∆\tau和∆r均为裂纹所在平面内的切应变；J为断裂韧性常数，取0.2。

根据式（10）可得疲劳参数存在最大值F_{{\mathrm{P}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}，当时，疲劳不会产生，则结合Manson-Coffin公式^［35-36］，可得钢轨疲劳寿命N_{\mathrm{f}}与F_{{\mathrm{P}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}的关系即拉伸型和剪切型疲劳裂纹萌生寿命预测经验式，钢轨疲劳裂纹损伤累积的具体推导过程，可参考文献［37］。

则根据前节静态损伤概念，提出疲劳损伤参数D_{\mathrm{f}}，定义值在0~1之间变化，按照Miner累积法则有

\frac{\mathrm{d}{D}_{\mathrm{f}}}{\mathrm{d}{N}_{\mathrm{f}}}=\frac{\mathrm{1}}{N_{\mathrm{f}}}=\frac{(F_{\mathrm{P}}-F_{{\mathrm{P}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{)}^H}{C}（11）

式中：N为循环加载次数；C和H为材料控制常数。

在每次循环加载结束前计算1次疲劳参数F_{\mathrm{P}}。

根据式（11）计算并更新疲劳损伤参数D_{\mathrm{f}}。在计算当前循环总损伤D_{\mathrm{a}}时，将提出的疲劳损伤参数D_{\mathrm{f}}与前节的静态损伤参数D_{\mathrm{m}}相结合，描述内聚力单元材料的弱化过程^［38-39］，定义D_{\mathrm{a}}为

D_{\mathrm{a}}=\int \mathrm{m}\mathrm{a}\mathrm{x}\left({\dot{D}}_{\mathrm{m}},{\dot{D}}_{\mathrm{f}}\right)\mathrm{d}t（12）

内聚力单元受牵引-分离准则和疲劳损伤累积准则的共同控制，计算时每个单元积分点有2个损伤参数，对应于疲劳损伤参数D_{\mathrm{f}}和静态损伤参数D_{\mathrm{m}}。在每个循环中，选择损伤变化率更大的值存储为当前循环总损伤D_{\mathrm{a}}，随着循环次数的累计，直到D_{\mathrm{a}}达到1，则表明单元完全失效，定义该单元为裂纹，并将其从模拟中删除。在晶粒和晶界中嵌入的内聚力单元，作为裂纹的潜在萌生扩展路径，即可模拟钢轨材料晶界或者晶粒内部的RCF裂纹演化行为。

1.3　微观晶粒几何模型

1.3.1　几何结构

应用Voronoi原理，提出基于EBSD观测结果的材料微观晶粒几何建模方法。Voronoi原理为在特定区域内利用1组种子点构建Voronoi多边形，使每个多边形内任意点到与之对应的种子点的距离比到其他种子点都近^［40］，这与晶粒的结晶生长过程极为类似，每个多边形即是晶粒。

二维微观晶粒几何模型建立的主要步骤为：①利用Python软件中Scipy库函数提供的Voronoi类函数生成Voronoi图的种子点并迭代，使最终生成的多边形图尺寸分布符合EBSD观测结果（以生成的晶粒尺寸统计结果的期望\mu和方差\mathrm{\sigma }是否与EBSD观测结果一致为判据）；②提取步骤①中完整的Voronoi图中多边形各顶点信息，构造每个晶粒的边和面；③将步骤②中构造好的数据按照点、线、面的顺序输入到txt文件中，方便后续Abaqus软件调用；④在Abaqus/CAE软件中通过Python脚本接口读取步骤③中txt数据，按照点、线、面的规则存放于多级列表中；⑤在Part模块调用步骤④中数据，利用WirePloyLine命令生成晶粒的边，然后用AddFaces命令将生成的边组成面；⑥在Assembly模块中把步骤⑤中生成的所有面通过布尔运算相加，组装成实例，即可在Abaqus软件中得到二维微观晶粒几何模型如图1所示。

1.3.2　晶粒取向信息

除了晶粒尺寸存在差异之外，每个晶粒的晶体学取向也不同，因此需要对每个晶粒取向进行赋予。由于子程序计算时采用的是Miller指数方式进行输入，但实际EBSD观测中得到的晶粒取向数据为Euler表述，因此，首先需要从EBSD观测的ctf文件中提取3组Euler角参数，然后通过取向矩阵\mathit{g}将其转化为Miller指数，Euler角取向与取向矩阵\mathit{g}满足关系

\mathit{g}={\mathit{g}}_{{\phi }_{\mathrm{2}}}{\mathit{g}}_{\varphi }{\mathit{g}}_{{\phi }_{\mathrm{1}}}（13）

其中，

\begin{array}{l}{\mathit{g}}_{{\phi }_{\mathrm{1}}}=\left(\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{1}}& \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{1}}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right)\\ {\mathit{g}}_{\varphi }=\left(\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\varphi & \mathrm{s}\mathrm{i}\mathrm{n}\varphi \\ \mathrm{0}& -\mathrm{s}\mathrm{i}\mathrm{n}\varphi & \mathrm{c}\mathrm{o}\mathrm{s}\varphi \end{array}\right)\\ {\mathit{g}}_{{\phi }_{\mathrm{2}}}=\left(\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}& \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right)\end{array}

式中：{\phi }_{\mathrm{1}}，{\phi }_{\mathrm{1}}和{\phi }_{\mathrm{2}}分别为晶粒绕局部坐标系3个方向轴的转动欧拉（Euler）角。

将转化的Miller指数保存为txt文本文件，采用Python脚本实现对每个晶粒取向的赋予：首先遍历区域（Part）中的所有晶粒，将每个晶粒建立为单独的部分（Set）；然后通过f.read函数命令读取txt文件中每个晶粒的取向数据；最后由For循环对Set进行取向赋予，结果如图2所示。图中：不同颜色代表不同取向的晶粒。

综上所述，首先基于晶体塑性理论，模拟钢轨材料微观结构的变形行为；然后基于Voronoi原理，建立了基于EBSD观测结果的材料真实微观晶粒模型，最后将损伤累积法则引入内聚力模型，并通过在晶粒和晶界处引入内聚力单元模拟材料损伤演化，最终实现在微观尺度模拟RCF裂纹演化。

2 RCF裂纹演化模型建立

以U75V钢轨为研究对象，这种珠光体钢具有体心立方（Body-Centered Cubic，BCC）晶体结构，主要用于国内重载铁路。U75V钢轨的化学成分、单调拉伸的应力-应变曲线来自文献［20］、在滚动接触荷载下的磨耗、表面塑性特征及疲劳损伤来自文献［9］。钢轨RCF裂纹演化模型示意图如图3所示。

为了得到前节建立的微观晶粒模型的建模参数，并定量分析U75V钢轨材料在滚动接触条件下的演化行为，基于建立的RCF裂纹演化模型，开展小比例双轮滚动疲劳试验。首先基于EBSD对钢轨滚轮的微观结构特征进行表征，为建立微观晶粒模型提供计算依据；然后建立双轮有限元宏观模型，为微观模型计算提供相应的边界条件与激励；最后将建立的微观模型与EBSD观测结果进行对比验证。

2.1　U75V钢轨材料微观晶粒特征

对小比例双轮疲劳试验研究中^［9］的钢轨滚轮进行EBSD观测，对晶粒尺寸分布、晶体取向等信息进行统计分析。用EBSD表征的U75V钢轨材料结果如图4所示。

对U75V钢轨滚轮试样的晶粒取向和尺寸进行统计（取向差角>15°，为大角度晶界；反之为小角度晶界），结果如图5所示。图中：d₉₀=35.6 μm，表示90%的晶粒尺寸都集中在35.6 μm以下；d₅₀和d₁₀同理。

从图5可以看出：U75V钢轨材料中主要为小角度晶界，占60%左右，同时伴有少量的大角度晶界，呈对数正态分布；U75V钢轨材料中晶粒尺寸主要集中在10~15 μm左右，占比高达55.2%，只有10%的晶粒尺寸小于10 μm；晶粒尺寸分布规律和取向差角类似，符合对数正态分布，期望μ=11.46，方差σ=1.07。

2.2　RCF微观损伤模型

根据文献［9］中小比例双轮试验工况，建立双轮接触有限元宏观模型，并对双轮接触区处的局部网格进行细化，如图6所示。

建立的模型共包括不参与损伤计算的双轮试验有限元局部模型和二维微观损伤子模型2部分。

局部模型主要为损伤子模型提供合理边界，切割自双轮试验宏观模型，边界条件继承自宏观模型。根据Bomidi等^［41］的相似理论建立的二维局部模型整体尺寸应大于10a×4a（a为接触半宽），该尺寸可以在保证计算精度的基础上有效减少计算量，同时确保模型边界不会影响荷载分布，设置为5 000 μm×2 000 μm。

二维微观损伤子模型采用上文方法建立，损伤子模型尺寸为1 000 μm×500 μm（根据文献［9］中的金相观测，该尺寸足以包括接触表面所有损伤发生的范围，且远离局部模型边界），边界条件继承自局部模型，损伤仅允许在此区域发生；接触荷载选择提取双轮试验接触斑纵向中心位置的垂向和纵向应力，并假设接触过程中荷载大小不发生变化，接触荷载从损伤子模型一端移动到另一端为1次循环。二维RCF裂纹演化模型建立示意图如图7所示。

对建立的U75V钢轨二维RCF裂纹演化模型进行极图与反极图绘制^［42］，并与真实材料的EBSD极图结果进行对比，如图8所示。钢轨材料晶粒尺寸分布对比如图9所示。

从图8和图9可以看出：建立的微观晶粒模型与真实钢轨材料具有相同的晶粒取向分布；将建立的微观仿真模型晶粒尺寸分布与钢轨材料晶粒EBSD实际观测数据进行对比，2者具有相同的正态分布特征，可以认为建立的微观仿真模型晶粒特征在统计学意义上与真实钢轨材料是一致的。

3 钢轨微观晶粒结构中裂纹扩展特性

钢轨在累积20 000次循环后RCF裂纹扩展情况仿真结果和实测结果如图10所示。图中：1#—9#为9条裂纹编号；a，b，c，d为4条裂纹编号；另外，图中隐去了晶粒形状和网格单元，方便观测裂纹扩展路径。

从图10可以看出：钢轨材料晶粒结构相对均匀，裂纹均为与接触表面有一定角度的斜裂纹，与双轮试验实测结果一致；裂纹的出现会改变材料内部的应力状态，在其扩展路径周边有应力集中现象，裂纹尖端会出现较大的应力集中；仿真和实测的裂纹扩展角度均在10°~20°左右，但裂纹扩展程度并不一致；部分裂纹扩展长度已超出100 μm（1#，2#，4#和6#），另一部分裂纹处于刚萌生，或者仅扩展十几微米左右（3#，5#，7#，8#和9#），这是因为仿真考虑了微观晶粒结构的影响，不同晶粒之间的取向差较大，晶界间相互作用也不一致，这导致在相同外部荷载条件下，裂纹的具体扩展路径和扩展程度不尽相同，这与实测结果相符；对仿真结果中已经成形的主要裂纹扩展长度进行统计（命名为主裂纹，长度>100 μm），平均扩展长度为126.7 μm，而试验观测的裂纹扩展长度均值为102.3 μm，仿真裂纹扩展长度大于试验观测结果，误差为19.3%；考虑双轮试验在20 000次循环时的垂直磨耗深度为32.04 μm^［9］，在仿真结果中减去垂直磨耗，3#和7#裂纹会直接被磨损，其余主裂纹扩展长度也会大幅度缩减，此时仿真裂纹平均扩展长度为94.5 μm，略小于实际观测长度，误差为7.6%。

仿真的裂纹扩展路径并不平滑，为明确循环荷载过程中裂纹扩展行为，裂纹在不同循环次数下（每4 000次）的扩展情况如图11所示。

从图11可以看出以下结论。

（1）在循环早期，接触表面只有零星的应力集中现象，未出现裂纹萌生，当循环4 000次时，仅有4#裂纹在晶粒内部萌生；随着循环次数的累积，4#裂纹沿滚动方向平行扩展，其他裂纹也陆续萌生和扩展，在相同外部荷载下，裂纹萌生寿命和扩展速率并不一致。

（2）图10中，将循环结束时的主裂纹编为9条，但在循环过程中，实际裂纹萌生数量远多于9条。萌生位置间隔较近的微裂纹互相吸引，并在循环荷载作用下互相连接，最终形成循环结束时的9条主裂纹。

（3）裂纹大多萌生于晶界处，萌生后大部分沿晶界继续扩展，少部分直接向晶粒内部穿晶扩展，裂纹扩展过程中，沿晶和穿晶扩展一直伴随其中。

4#裂纹扩展路径非标准斜裂纹，路径呈半“Z”字型，在实际观测中也发现了类似裂纹走向（图10（b）中的c裂纹和d裂纹）。为探究这一路径的形成过程，对4#裂纹区域进行放大跟踪观测，其在不同循环次数下的扩展情况如图12所示。图中隐藏网格、取向配色和应力结果。

从图12可以看出：在4 000次循环时，微裂纹在晶粒内部首次出现，萌生并扩展约1个晶粒尺寸，在6 000次循环时，在附近又萌生了新的微裂纹，微裂纹继续扩展，在8 000次循环时相互连接成为1个整体，暂且命名为4^*裂纹；在16 000次循环时，4^*裂纹的右下角萌生了多个与滚动方向形成一定角度的微裂纹；在18 000次循环时，4^*裂纹右下角萌生的微裂纹已经基本连成一片；当循环达到20 000次时，4^*裂纹与右下角萌生的裂纹已经互相连接成为1个整体，形成循环结束时的4#裂纹。

将循环荷载作用下裂纹演化归为Ⅰ阶段和Ⅱ阶段，Ⅰ阶段包括微裂纹萌生和扩展的过程，Ⅱ阶段为微裂纹连接成片的过程，示意图如图13所示。图中：黑线表示晶界；黄线表示微裂纹；红色虚线表示均质材料时，在该外部荷载工况下裂纹扩展方向理论计算值。

从图13可以看出：微裂纹在荷载作用下大多萌生于晶界处，单个微裂纹的萌生位置和寿命具有随机性，但萌生位置在特定区域内呈现集中趋势，区域内萌生位置间隔较近的微裂纹在扩展过程中会相互吸引聚合，在较高循环次数后形成主裂纹；Ⅰ阶段为起步阶段，微裂纹扩展长度在该阶段通常为几个晶粒尺寸，该阶段裂纹的萌生扩展行为主要受晶粒取向和晶界间相互作用支配；主裂纹成形后进入Ⅱ阶段，此时多发的微裂纹连接形成1条新裂纹，新裂纹的扩展行为主要受外部荷载的支配，材料微观晶粒结构的影响大幅减弱，从而形成了在循环荷载结束时各条主裂纹扩展路径及角度高度相似的现象。

根据最不利原则，选取模型中最早萌生的微裂纹代表该模型的裂纹萌生寿命，选取模型中扩展最长的裂纹代表该模型中裂纹扩展速率，则选择4#裂纹代表（裂纹长度统计包括4#裂纹形成过程中的各微裂纹长度之和），裂纹萌生寿命为3 472次以单元疲劳累积到1，单元失效并删除时对应的荷载循环次数定义为裂纹萌生），裂纹在循环荷载期间的扩展速率和不同循环次数的相对裂纹扩展速率如图14所示。

从图14可以看出：平均裂纹扩展速率为0.010 3 μm · r^-1；在循环初期，裂纹此时属于Ⅰ阶段，在相同循环次数下的裂纹扩展长度较短，该阶段裂纹扩展速率为0.002 5 μm · r^-1左右，远低于裂纹平均扩展速率；随着裂纹继续扩展，相同循环次数下裂纹的扩展长度逐渐加快，在进入Ⅱ阶段后裂纹扩展速率明显增加，最高可达0.018 7 μm · r^-1。裂纹在不同阶段的扩展速率不同，随扩展长度的增加，整体呈早期慢后期逐渐变快的趋势。

4 结 论

（1）基于晶体塑性理论和Voronoi原理建立的反映U75V钢轨材料晶粒几何结构、取向等特征的微观模型，并通过修正的内聚力单元构建的材料微观结构与疲劳损伤之间的联系，实现了在微观尺度模拟滚动接触疲劳裂纹的萌生与扩展。

（2）微裂纹在荷载作用下大多萌生于晶界处，单个微裂纹的萌生位置和寿命具有随机性；但萌生位置在特定区域内呈现集中趋势，且该区域分布与宏观累积塑性应变的分布高度相关，区域内萌生位置间隔较近的微裂纹在扩展过程中会相互吸引聚合，在较高循环次数下形成新的裂纹。

（3）循环荷载作用下，接触表面裂纹的演化分为2个阶段：Ⅰ阶段为起步阶段，微裂纹扩展长度在该阶段为几个晶粒尺寸，该阶段裂纹的萌生扩展行为主要受晶粒取向和晶界间相互作用支配；当萌生位置间隔较近的微裂纹在扩展过程中相互吸引聚合，形成1条新裂纹时，进入Ⅱ阶段，新裂纹的扩展行为主要受外部荷载的支配，此时材料晶粒结构的影响大幅减弱，从而形成在循环结束时各条裂纹扩展路径及角度高度相似的现象。

（4）滚动接触疲劳裂纹在不同阶段的扩展速率不同，随扩展长度的增加呈早期慢后期逐渐变快的趋势。在循环初期，相同循环次数下的裂纹扩展长度较短，裂纹扩展速率在0.002 5 μm · r^-1左右，远低于裂纹平均扩展速率0.010 3 μm · r^-1；随着裂纹继续扩展，相同循环次数下裂纹的扩展长度逐渐加快，裂纹扩展速率明显增加，最高可达0.018 7 μm · r^-1。

参考文献

[1] WEN Z F， TAO G Q， ZHAO X， et al. Wear and RCF Problems of Metro Wheel/Rail Systems： Phenomena， Causes and Countermeasures in China ［J］. Wear， 2023， 534： 205118.

[2] WANG W J， LI S J， DING H H， et al. Wheel/Rail Adhesion and Damage under Different Contact Conditions and Application Parameters of Friction Modifier ［J］. Wear， 2023， 523： 204870.

[3] 赵鑫，温泽峰，王衡禹，等.中国轨道交通轮轨滚动接触疲劳研究进展［J］.交通运输工程学报，2021，21（1）：1-35.

[4] ZHAO Xin， WEN Zefeng， WANG Hengyu， et al. Research Progress on Wheel/Rail Rolling Contact Fatigue of Rail Transit in China ［J］. Journal of Traffic and Transportation Engineering， 2021， 21 （1）： 1-35. in Chinese

[5] 张银花，周清跃，陈朝阳，等.重载铁路钢轨现状及发展方向探讨［C］//2014年铁路重载运输技术交流会.北京：中国铁道学会， 2014： 381-386.

[6] ZHANG Yinhua， ZHOU Qinyue， CHEN Chaoyang， et al. Discussion on Current Situation and Development Direction of Heavy Haul Railway Rail ［C］// 2014 Heavy Haul Railway Technology Conference. Beijing： China Railway Society， 2014： 381-386. in Chinese）

[7] 许鑫，徐超，张倩，等.60N U71MnG钢轨轨头裂纹成因分析［J］.中国铁道科学，2023，44（5）：47-57.

[8] XU Xin， XU Chao， ZHANG Qian， et al. Cause Analysis of Cracks in the 60N U71MnG Rail Head ［J］. China Railway Science， 2023， 44 （5）： 47-57. in Chinese

[9] GUILHEM Y， BASSEVILLE S， CURTIT F， et al. Investigation of the Effect of Grain Clusters on Fatigue Crack Initiation in Polycrystals ［J］. International Journal of Fatigue， 2010， 32 （11）： 1748-1763.

[10] WALVEKAR A A， SADEGHI F. Rolling Contact Fatigue of Case Carburized Steels ［J］. International Journal of Fatigue， 2017， 95： 264-281.

[11] PEREIRA H B， ECHEVERRI E A A， ALVES L H D， et al. Influence of HAZ Microstructure on RCF under Twin-Disc Test of a Flash-Butt Welded Rail ［J］. Wear， 2024， 546： 205324.

[12] LI J P， ZHOU Y， CHENG Z N， et al. Evolution of Microscopic Characteristics of Rolling Contact Fatigue Coexisting with Wear of U75V Rail ［J］. Tribology International， 2023， 190： 109041.

[13] KUMAR A， DAWSON P R. Modeling Crystallographic Texture Evolution with Finite Elements over Neo-Eulerian Orientation Spaces ［J］. Computer Methods in Applied Mechanics and Engineering， 1998， 153 （3/4）： 259-302.

[14] RAABE D， ROTERS F. Using Texture Components in Crystal Plasticity Finite Element Simulations ［J］. International Journal of Plasticity， 2004， 20 （3）： 339-361.

[15] ZHU J C， BETTAIEB M B， ZHOU S， et al. Ductility Limit Prediction for Polycrystalline Aggregates Using a CPFEM-Based Multiscale Framework ［J］. International Journal of Plasticity， 2023， 167： 103671.

[16] YAN B Y， LIN P， YANG L， et al. Revealing Size Effect for Plastic Deformation of Nanocrystalline NiTi SMA at Moderate Temperature by Combining CPFEM with Experiment ［J］. Journal of Materials Research and Technology， 2024， 30： 7869-7884.

[17] PANDKAR A S， ARAKERE N， SUBHASH G. Microstructure-Sensitive Accumulation of Plastic Strain Due to Ratcheting in Bearing Steels Subject to Rolling Contact Fatigue ［J］. International Journal of Fatigue， 2014， 63： 191-202.

[18] VIJAY A， SADEGHI F. A Continuum Damage Mechanics Framework for Modeling the Effect of Crystalline Anisotropy on Rolling Contact Fatigue ［J］. Tribology International， 2019， 140： 105845.

[19] VIJAY A， SADEGHI F. Rolling Contact Fatigue of Coupled EHL and Anisotropic Polycrystalline Materials ［J］. Tribology International， 2022， 169： 107479.

[20] GHODRATI M， AHMADIAN M， MIRZAEIFAR R. Modeling of Rolling Contact Fatigue in Rails at the Microstructural Level ［J］. Wear， 2018， 406： 205-217.

[21] GHODRATI M， AHMADIAN M， MIRZAEIFAR R. Three-Dimensional Study of Rolling Contact Fatigue Using Crystal Plasticity and Cohesive Zone Method ［J］. International Journal of Fatigue， 2019， 128： 105208.

[22] WANG W， LIU H J， ZHU C C， et al. Micromechanical Analysis of Gear Fatigue-Ratcheting Damage Considering the Phase State and Inclusion ［J］. Tribology International， 2019， 136： 182-195.

[23] FAN Y L， KAN Q H， ZHAO J Z， et al. Experimental Study on Ratcheting Behavior of Different Rail Steels ［J］. Journal of Mechanical Engineering， 2020， 56： 35.

[24] HILL R. Generalized Constitutive Relations for Incremental Deformation of Metal Crystals by Multislip ［J］. Journal of the Mechanics and Physics of Solids， 1966， 14 （2）： 95-102.

[25] HILL R， RICE J R. Constitutive Analysis of Elastic-Plastic Crystals at Arbitrary Strain ［J］. Journal of the Mechanics and Physics of Solids， 1972， 20 （6）： 401-413.

[26] PEIRCE D， ASARO R J， NEEDLEMAN A. An Analysis of Nonuniform and Localized Deformation in Ductile Single Crystals ［J］. Acta Metallurgica， 1982， 30 （6）： 1087-1119.

[27] ASARO R J， RICE J R. Strain Localization in Ductile Single Crystals ［J］. Journal of the Mechanics and Physics of Solids， 1977， 25 （5）： 309-338.

[28] HAVNER K S， BAKER G S， VAUSE R F. Theoretical Latent Hardening in Crystals — I. General Equations for Tension and Compression with Application to F.C.C. Crystals in Tension ［J］. Journal of the Mechanics and Physics of Solids， 1979， 27 （1）： 33-50.

[29] HILL R. The Elastic Behaviour of a Crystalline Aggregate ［J］. Proceedings of the Physical Society， Section A， 1952， 65 （5）： 349-354.

[30] HUANG Y G. A User-Material Subroutine Incorporating Single Crystal Plasticity in the ABAQUS Finite Element Program ［D］. Cambridge： Harvard University， 1991.

[31] 罗娟，康国政，董亚伟，等.多晶铜和铝单轴棘轮行为的循环晶体塑性本构模拟［J］.四川大学学报（工程科学版），2012（增2）：80-83.

[32] LUO Juan， KANG Guozheng， DONG Yawei， et al. Modeling for Ratcheting of Polycrystalline Copper and Aluminum by Cyclic Crystal Plasticity ［J］. Journal of Sichuan University （Engineering Science Edition）， 2012 （）： 80-83. in Chinese

[33] 阚前华，康国政，徐祥.非线性本构关系在ABAQUS中的实现［M］.北京：科学出版社， 2019.

[34] KAN Qianhua， KANG Guozheng， XU Xiang， et al. Implementations of Nonlinear Constitutive Relations in ABAQUS ［M］. Beijing： Science Press， 2019. in Chinese

[35] JIANG Y Y， SEHITOGLU H. A Model for Rolling Contact Failure ［J］. Wear， 1999， 224 （1）： 38-49.

[36] SIMULIA. Abaqus 2016 Documentation ［EB/OL］. 2015.

[37] DONZELLA G， FACCOLI M， GHIDINI A， et al. The Competitive Role of Wear and RCF in a Rail Steel ［J］. Engineering Fracture Mechanics， 2005， 72 （2）： 287-308.

[38] BENZEGGAGH M L， KENANE M. Measurement of Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus ［J］. Composites Science and Technology， 1996， 56 （4）： 439-449.

[39] JIANG Y Y， SEHITOGLU H. Rolling Contact Stress Analysis with the Application of a New Plasticity Model ［J］. Wear， 1996， 191 （1/2）： 35-44.

[40] MANSON S S. Behavior of Materials under Conditions of Thermal Stress ［M］. Washington， D.C： National Advisory Committee for Aeronautics， 1953.

[41] COFFIN L F. A Study of the Effects of Cyclic Thermal Stresses on a Ductile Metal ［J］. Journal of Fluids Engineering， 1954， 76 （6）： 931-949.

[42] LI J P， WENG Z Y， ZHOU Y， et al. Prediction Model of Rail Head Check Initiation Considering Wear and Nonlinear Fatigue Damage Accumulation ［J］. Fatigue & Fracture of Engineering Materials & Structures， 2023， 46 （12）： 4591-4605.

[43] ROE K L， SIEGMUND T. An Irreversible Cohesive Zone Model for Interface Fatigue Crack Growth Simulation ［J］. Engineering Fracture Mechanics， 2003， 70 （2）： 209-232.

[44] ZHANG J W， ZHU S Y， CAI C B， et al. Cohesive Zone Modeling of Fatigue Crack Propagation in Slab Track Interface under Cyclic Temperature Load ［J］. Engineering Failure Analysis， 2022， 134： 106028.

[45] SCHNEIDER Y， WEBER U， WASSERBÄCH W， et al. A Numerical Method for the Generation of Hierarchical Poisson Voronoi Microstructures Applied in Micromechanical Finite Element Simulations， Part I： Method ［J］. Computational Mechanics， 2020， 66 （3）： 651-667.

[46] BOMIDI J A R， WEINZAPFEL N， SADEGHI F， et al. An Improved Approach for 3D Rolling Contact Fatigue Simulations with Microstructure Topology ［J］. Tribology Transactions， 2013， 56 （3）： 385-399.

[47] HIELSCHER R， SCHAEBEN H. A Novel Pole Figure Inversion Method： Specification of the MTEX Algorithm ［J］. Journal of Applied Crystallography， 2008， 41 （6）： 1024-1037.

基金资助

中央高校基本科研业务费项目(2022-5-ZD-04)

上海市科委基金资助项目(20dz1203100)

浙江省交通运输厅科技计划项目(2023024)

中国国家铁路集团有限公司科技研究开发项目(N2022G011)

中国国家铁路集团有限公司科技研究开发项目(K2023G013)

中国铁道科学研究院集团有限公司院基金课题(2023YJ052)

安徽省城市轨道交通安全与应急管理重点实验室开放课题(2024GD0008)

2024年度合肥大学人才科研基金项目(24RC24)