1.一种ABAQUS中铁轨滚动接触疲劳裂纹的数值分析方法，其特征在于，按照以下步骤实施：步骤1：建立轨道有限元微观模型，

设车轮与轨道接触的裂纹分析区为6b*4b，b为轨道与车轮滚动接触的接触半宽，将其中的2b*b区定为中央区，使用Neper软件在轨道裂纹分析区建立Voronoi及内聚力单元镶嵌的轨道有限元微观模型；

步骤2：导入步骤1建立的轨道有限元微观模型，生成ABAQUS的轨道计算模型，利用ABAQUS的“Import”功能，将镶嵌Voronoi以及内聚力单元的轨道有限元微观模型文件mode.inp导入ABAQUS中；

步骤3：采集并初始化参数，设置轨道计算模型边界条件，根据轨道的材料初始化相关参数：晶粒的弹性模量E，泊松比γ，晶界内聚力单元初始所能承受的最大应力值T0，内聚力单元的初始刚度K0，断裂能GIC；并对疲劳退化损伤变量以及累积损伤变量初始赋值为0，即Ds＝0，Df＝0；

并且，根据轨道与车轮接触的工况，在ABAQUS中限制轨道底部的全部自由度；

步骤4：设置ABAQUS中的分析步和增量步；

步骤5：计算材料损伤性能，

利用疲劳退化损伤变量Df以及累积损伤变量Ds，修改材料内聚力单元刚度以及最大承载应力，计算得到内聚力单元退化后所能承受的应力Tf，计算得到损伤后的内聚力单元刚度K1；

步骤6：对轨道计算模型进行加载；

步骤7：分析内聚力单元的应力幅Δτ以及应变δ，

利用ABAQUS的静力学分析功能，得到轨道模型在滚动接触载荷作用下的单元最大应力幅Δτ和应变δ；

步骤8：分析内聚力单元的损伤，

在ABAQUS的job界面，添加P1.for子程序，将内聚力单元疲劳损伤失效的UMAT子程序P1.for关联到轨道分析模型，调用UMAT子程序P1.for计算轨道模型的内聚力单元损伤变量Df以及Ds；

步骤9：判断内聚力单元是否失效，

若Ds≠1，返回步骤5计算材料损伤性能，循环加载计算；

若Ds＝1，单元失效，进入步骤10；

步骤10：修订ABAQUS计算模型，

如果Ds＝1说明内聚力单元已失效，将该晶界的内聚力单元从模型中删除，并在模型中标记开始萌生裂纹；然后，返回步骤5计算材料损伤性能，循环加载计算；

在ABAQUS提交上述分析模型，ABAQUS将按设置的分析步以及增量步，循环轨道滚动接触疲劳裂纹萌生与扩展的数值分析算法中的步骤5到步骤10，从而实现对轨道裂纹的萌生与扩展分析，及时输出分析结果。

2.根据权利要求1所述的ABAQUS中铁轨滚动接触疲劳裂纹的数值分析方法，其特征在于：所述的步骤1中，使用基于Linux系统的开源软件Neper建模，具体过程是，

1.1)生成Voronoi晶粒位置文件，

利用Matlab的随机数产生函数rand()，在轨道的裂纹分析区6b*4b随机生成700个晶粒种子位置，将所有晶粒种子位置存贮在晶粒位置文件c1.txt中；

1.2)构建轨道Voronoi模型，

利用Neper的-T命令，由晶粒位置文件c1.txt，构建含有700个晶粒位置的轨道Voronoi模型文件mode_V.tess，其命令格式为：neper-T-n 700-dim 2-morphooptiini"coo:file(c1)"-domain"square(6b，4b)"-morpho gg-o mode_V；

1.3)划分镶嵌内聚力单元及有限元网格，

利用Neper的-M命令，对轨道Voronoi模型进行有限元网格划分，有限元网格单元类型为CPE3，并在晶粒间嵌入类型为COH2D4内聚力单元，生成镶嵌Voronoi以及内聚力单元的轨道有限元微观模型文件mode.inp，Neper的命令格式如下：neper-M mode_V.tess-dim all-nset 0，1，2，3-interface cohesive-format inp-rcl 0.5-o mode。

3.根据权利要求2所述的ABAQUS中铁轨滚动接触疲劳裂纹的数值分析方法，其特征在于：所述的步骤1.1)中，在裂纹分析区的中央区随机生成350个晶粒种子位置，在裂纹分析区的其它位置生成350个晶粒种子位置。

4.根据权利要求1所述的ABAQUS中铁轨滚动接触疲劳裂纹的数值分析方法，其特征在于：所述的步骤4中，对于轨道的滚动接触载荷，设置一个分析步，增量步的范围设为0.001～0.01。

5.根据权利要求1所述的ABAQUS中铁轨滚动接触疲劳裂纹的数值分析方法，其特征在于：所述的步骤6中，具体过程是，当车轮在轨道上滚动时，车轮与轨道间的作用力包括轨道表面所受的法向力p(x)与切向力t(x)，其中的法向力p(x)分布表达式如下式：式(6)中，xc为载荷中心位置横坐标；pmax为最大接触应力，单位为Mpa；b为轨道与车轮滚动接触的接触半宽，单位为mm；

切向力大小与牵引系数有关，其表达式为：

t(x)＝μp(x)(7)

式(7)中，μ为牵引系数。

6.根据权利要求1所述的ABAQUS中铁轨滚动接触疲劳裂纹的数值分析方法，其特征在于：所述的步骤8中，调用内聚力单元疲劳损伤失效的UMAT子程序P1.for，程序输入为：内聚力单元所受应力幅Δτ以及应变δ；

程序输出为：内聚力单元的损伤变量，包括疲劳退化损伤变量Df以及累积损伤变量Ds；

程序中材料的弹性模量E、泊松比γ、晶界内聚力单元初始所能承受的最大应力值T0、内聚力单元的初始刚度K0及断裂能GIC常量，该UMAT子程序的工作过程，按照以下步骤实施：

第一步：导入单元所受应力幅Δτ以及应变δ，

UMAT子程序通过数据交流接口从ABAQUS主程序中获取计算模型在外载荷作用下，内聚力单元所受的应力幅Δτ以及应变δ值；

第二步：判断内聚力单元是否开始损伤，

根据损伤起始条件，若Δτ≥Tf，判断内聚力单元开始损伤；若Δτ＜Tf，判断内聚力单元还未损伤，内聚力单元处于疲劳退化阶段；

第三步：计算内聚力单元的损伤变量Df以及DS，

若Δτ≥Tf，内聚力单元开始损伤，根据式(5)计算累积损伤变量DS；

若Δτ＜Tf，内聚力单元未损伤，计算疲劳退化损伤变量Df；

最后输出累积损伤变量Ds或疲劳退化损伤变量Df。