ABAQUS---轮轨瞬态滚动接触有限元模型(直线半轮对)

      目前，轮轨瞬态滚动接触有限元模型日渐成熟，尤其针对直线半轮对情况。利用该模型已经详细开展了大量的轮轨滚动仿真，比如：1）轮轨不平顺（钢轨波磨、焊接接头、硌伤、隐伤；车轮多边形、擦伤、凹磨）；2）道岔瞬态冲击振动；3）单点-两点接触；4）轮轨低黏着；5）热机耦合，并分析了各种情形下的轮轨滚动接触力学行为、磨耗和疲劳损伤问题。然而，该成熟的模型大多都是基于ANSYS软件建立，而ABAQUS软件本身在模拟强非线性接触、材料塑性本构、CAE界面操作等方面具有显著的优势，但是当下基于ABAQUS软件建立的轮轨瞬态滚动接触模型仍存在很多问题，比如：轮轨力不稳定、车轮网格沙漏引起畸变、牵引/制动模拟困难、一系耦合约束和扣件模拟不当等，使得该模型推广受阻。本文旨在从作者经验角度，分享轮轨滚动接触有限元建模时可能面临的问题，如有不当，还欢迎批评指正。

       轮轨瞬态滚动接触有限元模型中，由于车轮具有较高的滚动速度，使得车轮瞬态滚动时对系统激扰较大，轮轨接触力稳定困难。因此，采用隐式-显式方法模拟瞬态滚动接触行为，其中隐式模型可得到车轮在重力场下的轮轨静态位移和应力场结果，然后将其导入至显式模型中，再在显式模型中模拟车轮滚动。以下分别介绍这两个模型及其之间的关联。

        注意：模型中单位制为mm、t、N、MPa。1为横向x，2为纵向y，3为垂向z

1. 隐式-静态（模型名称：standard.inp）

（1）Part

    模型包含三维的车轮、钢轨和道床板，等效模拟的一系悬挂和扣件。

（2）Property

    分别对三维实体结构赋予材料属性，包括：密度、弹性参数（弹性模量、泊松比）、塑性参数（屈服强度、硬化参数）

（3）Assembly

    将车轮、钢轨和道床板进行装配。注意：由于是半轮对模型，轮轨一般是对中接触，所以装配时保证轮对内侧距和轨距正确，保证轮对中心线和轨道中心线对齐即可。如果需要施加车轮横移、侧滚和冲角，可由动力学模型计算得到。

（4）Step

    1）建立分析步：Step-1 1s Static, General Nlgeom: On

*Step, name=Step-1, nlgeom=YES, inc=100000
*Static
0.05, 1., 1e-08, 1.

    2）在standard模型属性的edit keywords中的“*End step”前加入语句*contact controls, stabilize，收敛速度可大大提高。

    3）重启动输出 Restart requests：选择Frequency=1，每个增量步都输出一个重启动文件；选择Intervals=1，只输出最后一个增量步结果。这里输出后，以供显式模型调用。

    4）场输出 Field output：E,S/U,V/CDISP,CFORCE,STRESS（every n increment, n=1）

    5）历史输出 History output：ALLAE,ALLSE,ETOTAL（every n increment, n=1）

（5）Interaction

    1）主从面选择：在initial分析步中，选择车轮踏面为主面，钢轨踏面为从面，滑移计算公式为有限滑移，面-面接触，不调整接触距离，无接触控制。

    注意：最好将轮轨接触间隙调整为10^-4-10^-5mm左右，如果网格大和间隙比较大，可适当设置接触间隙。

** Interaction: Int-1
*Contact Pair, interaction=IntProp-1, type=SURFACE TO SURFACE
Surf-rail, Surf-wheel

    2）法切向接触属性：法向接触选择硬接触，切向接触选择Penalty，摩擦系数为0.1-0.6（湿润-干燥），out-of-plane surface thickness or cross-sectional area(standard):1

*Surface Interaction, name=IntProp-1
1.,
*Friction, slip tolerance=0.005
 0.5,
*Surface Behavior, pressure-overclosure=HARD

    3）建立一系质量：由于构架、车轮等部件影响轮轨系统低频特性，而瞬态滚动接触模型主要关注高频特性，因此这里忽略其实体结构建模，而是进行等效建模。建立四个参考点：RP1、RP2、RP3、RP4，分别赋予其质量Inertia=1.80 t，相当于轮重为7.2t，轴重为14.4t。也可以施加转动惯量，这里忽略。

*Element, type=MASS, elset=Set-RP1_Inertia-1_
1, 1
*Mass, elset=Set-RP1_Inertia-1_
1.80,

    4）选择车轴中心网格节点，创建四个参考点：RP5、RP6、RP7、RP8。

    5）建立一系耦合关系：将参考点RP5、RP6、RP7、RP8分别与车轴表面进行运动耦合Kinematic，所有自由度U1,U2,U3,UR1,UR2,UR3均进行约束。

*Coupling, constraint name=Constraint-1, ref node=Set-RP5, surface=Mass1
*Kinematic

    6）建立一系悬挂：由RP1和RP5、RP2和RP6、RP3和RP7、RP4和RP8两参考点之间分别建立4条Spring A弹簧，以模拟一系悬挂，仅考虑其垂向特性时，设置单个垂向刚度为200 N/mm，单个垂向阻尼为2.00N s/mm。

*Spring, elset=Springs/Dashpots-yixi-z-spring
3, 3
200.0
*Dashpot, elset=Springs/Dashpots-yixi-z-dashpot
3, 3
2.00
*Element, type=SpringA, elset=Springs/Dashpots-yixi-z-spring
7, 1, 5
9, 2, 6
11, 3, 7
13, 4, 8
*Element, type=DashpotA, elset=Springs/Dashpots-yixi-z-dashpot
8, 1, 5
10, 2, 6
12, 3, 7
14, 4, 8

另外，将RP1和RP5、RP2和RP6、RP3和RP7、RP4和RP8设置为Kinematic运动耦合关系，自由度约束U1,U2,UR1,UR2,UR3。

*Coupling, constraint name=Constraint-5, ref node=Set-RP1, surface=Set-RP5_CNS_
*Kinematic
1, 1
2, 2
4, 4
5, 5
6, 6

    注意：这里也可以选用axial或者connector连接器模拟，后续在显式模型中会进行介绍。如果考虑三向特性，可以选择1,2,3三个方向，分别设置三个刚度和阻尼值，一般直线半轮对模型仅考虑垂向即可。另外，上下两两参考点之间的距离为320 mm。

    7）建立扣件：选择钢轨和轨道板对应的两个点建立Spring A弹簧，以模拟扣件，仅考虑其垂向特性时，设置单个垂向刚度为1050 N/mm, 单个垂向阻尼为2.33N s/mm。这里同样也可以选用axial或者connector连接器进行模拟。

*Spring, elset=Springs/Dashpots-koujian-z-spring
3, 3
1050.0
*Dashpot, elset=Springs/Dashpots-koujian-z-dashpot
3, 3
2.33

（6）Load

    1）载荷：Step-1施加系统重力，垂向，9810，无幅值曲线

    2）边界条件：在Initial 分析步

        ★ RP1、RP2、RP3、RP4：U1，U2约束，一系有压缩，放开U3。（U1=U2=UR1=UR2=UR3=0）

        ★ RP5、RP6、RP7、RP8：U1，U2约束，一系有压缩，放开U3。（U1=U2=UR1=UR2=UR3=0）

        ★ 由于扣件仅考虑了垂向，因此需要将钢轨底面或者仅扣件位置的横向U1约束（U1=0），以及钢轨两端面纵向对称约束（U2=UR1=UR3=0）。

        ★ 车轮中心面横向对称约束（U1=UR2=UR3=0），车轮另一端面U1、U2约束（U1=U2=0）

        ★ 轨道板端面纵向对称约束（U2=UR1=UR3=0），中心面横向对称约束（U1=UR2=UR3=0），底面固定约束（U1=U2=U3=0）。这里若不用固定约束，也可以设置接地弹簧。

（7）Mesh

    1）网格划分：车轮、钢轨、道床板三维实体结构，采用sweep扫掠网格进行划分。轮轨接触区域网格需细化为1 mm，其余区域可进行过渡划分，最大尺寸8 mm。

    注意：目前比较习惯用hypermesh划分网格，优点是网格数量少，质量高，计算速度比ABAQUS快很多，缺点是重新划分局部网格后，还需要导入ABAQUS中重新设置所有，比较麻烦。

    2）网格类型：standard，3Dstress，Linear，Hex，Reduced integration，C3D8R，八节点六面体缩减积分单元。

    3）网格畸变控制：Distortion control: Length ratio=0.5，Hourglass control: Enhanced

    Scaling factors: Displacement hourglass=1, Linear bulk viscosity=1, Quadratic bulk viscosity=1

    *Solid Section, elset=Set-rail, controls=EC-1, material=Material-rail

,

*Hourglass Stiffness
0.6, , 0., 0.
*End Part

---------

** ELEMENT CONTROLS
** 
*Section Controls, name=EC-1, DISTORTION CONTROL=YES, length ratio=0.5, hourglass=ENHANCED
1., 1., 1.

    注意：为了保证计算沙漏力一致，只需要在隐式模型中设置enhanced增量应变方法控制沙漏问题，在显式里CAE中以及inp文件里设置其参数都是没有用的。另外，enhanced方法比stiffness刚度系数方法来控制沙漏问题效果更明显。

    4）检查网格质量：包括分析检查，稳定时间增量估计，长宽比（小于10比较困难，50即可），角度（小于10，大于160）等。

2. 显式-动态（模型名称：explicit.inp）

    目的：将静态位移场结果以数据传递的方式导入显式动态分析模型中；方法：复制隐式模型，重新定义。

    数据传递时注意事项：

    ①对于需要传递数据的部件，其名称需要在原始分析和后续分析中保持一致。

    ②该方法在CAE中仅能导入原始分析中的网格和材料状态（应力和等效应变PEEQ等）。

    ③流体单元、无限单元、弹簧单元、质量单元、阻尼单元以及旋转惯性单元不能进行数据传递。

    ④已经定义了初始状态场的部件实体、网格及其相应的部件，都不能再在CAE中进行修改，不允许为其赋予新的截面属性、材料方向、定义质量和惯性矩，在后续分析中会自动替换为原始分析中的材料。

    ⑤在standard和explicit之间进行数据传递时，要保证沙漏的计算方法一致。

（1）Assembly

    复制隐式CAE模型，模型中部件名字千万不能更改替换。

    注意：取掉隐式模型中在edit keywords中加入的*contact controls, stabilize，该关键词不适用于显式，否则报错。

（2）Step

    1）删掉之前的step，重新创建动态，显式分析步

*Step, name=Step-1, nlgeom=YES
*Dynamic, Explicit
, 0.12
*Bulk Viscosity
0.06, 1.2

    2）建立分析步：Step-1 Dynamic, explicit 0.15s Nlgeom: On Automatic Time scaling factor=1

    3）重启动输出 Restart requests：选择选择Intervals=1，只输出最后一个增量步结果。不选择Overlay和Time marks

4）场输出 Field output：E,S/U,V,A/STRESS,CFORCE,FSLIPR,FSLIP,PPRESS （Evenly spaced time intervals, interval=1500）

*Output, field, number interval=1500
*Element Output, directions=YES
E, S

5）历史输出 History output

    ALLAE,ALLSE,ETOTAL 所有模型

    CFNM,CFN1,CFN2,CFN3,CFSM,CFS1,CFS2,CFS3,CFTM,CFT1,CFT2,CFT3,CAREA 接触对：int-rail-wheel

    U1,U3,A1,A3 集合：rail-1,2,3

    U1,U2,U3,V2,A1,A2,A3 集合：wheel-zhouxiang

    输出数据间隔，0.15s内总共输出1500个数据，采样间隔为1×10-4s。Evenly spaced time intervals, interval=1500

*Output, history, time interval=0.0001
*Contact Output, cpset=Int-rail-wheel
CAREA, CFN3, CFS1, CFS2