基于显式有限元的无砟轨道系力学行为分析与程序开发

甘 愿，黄钰润

(西南交通大学土木工程学院，四川成都 610031)

无砟轨道在我国应用广泛，但随着我国高速铁路的快速发展，行车速度越来越快，各部件的变形失效明显加快，轮轨系统的动力作用也越来越明显，在列车长期循环荷载和混凝土桥梁收缩徐变等作用下产生的桥墩变形和梁体结构变形、砂浆损伤、板端上拱变形等典型病害引起无砟轨道系统出现局部脱空，严重影响了无砟轨道系统的动力特性。

由于脱空的存在，系统的力学行为因可能的碰撞而演变为高度非线性，处理此类问题时，显式有限元相对于隐式有限元具有先天的优势，其是否收敛只与时间步长有关，具有良好的稳定性，且可直接解耦，不形成大型矩阵，占用内存空间较小。目前常用的商业显式有限元软件有LS-DYNA、ABAQUS/Explicit等，但由于代码未开源，计算过程不透明，也不利于进行后续的参数化设计、结构优化、集成计算等工作，故基于显式有限元分析方法，自主开发了应用于有脱空无砟轨道应力分析的通用计算软件。

1 有限元模型建立

1.1 无砟轨道有限元离散

目前无砟轨道的理论计算方法中运用最多的是叠合梁模型、梁板模型和梁体模型，因为本文研究重点是轨道系统在有脱空状态下应力状态及程序应用可行性，故采用梁体模型。用减缩六面体实体单元离散轨道板，砂浆填充层，混凝土底座；钢轨采用弹性点支承梁模型，用大转动空间梁单元离散；考虑扣件的尺寸效应，用弹簧单元连接钢轨节点与轨道板扣件尺寸范围内节点，在不考虑地基不均匀沉降时，同样也用弹簧单元模拟地基的支撑作用。

1.2 接触问题显式有限元列式

显式有限元处理接触问题的有限元理论与基本求解方法可以参考文献[1-5]，在这里仅给出不考虑摩擦接触的基本公式。

假定存在2个物体A与B，或者一个物体的2个部分A和B，在t1时刻占据空间VA和VB，其边界分别为ΓA与ΓB，当至t2时刻时，系统的空间位置或者构型发生改变，2个物体发生了接触。其基本控制方程与不含接触系统的控制方程是一致的，但在接触界面上需增加额外的动力学与运动学条件，运动学约束为不可侵入条件，即对于ΓA中任意一点xA与ΓB中距离最近一点xB满足:

gn=(xA-xB)·n≥0

(1)

式中:n为xB处ΓB的外法线。同时，法向作用力只能为压力:

Pn≤0

(2)

根据虚功原理，系统的运动方程弱形式可以表述为:

(3)

pN=kNgN(gN<0)

(4)

对式(3)进行离散，可以得到系统的半离散:

MÜ=f

(5)

式中:f=fext-fint+fc。fint为系统内力，fext为系统外力，fc为接触力，M为质量矩阵，Ü为加速度向量。

上式以中心差分法直接积分，由于本程序中采用集中质量矩阵，因此可以将式(5)解耦。如已求得tn时刻的解，时刻tn+1的位移为:

(6)

2 有限元程序设计与开发

有限元程序的设计内容大致可以分为前处理，分析计算，后处理3个部分。其中分析计算部分由于有限元分析方法其需要描述真实世界，需要对实例进行多个层级的抽象与离散，天然的与面向对象程序语言的形式逻辑吻合；同时从节点到单元和材料、对于各类荷载与边界条件的处理，再到集总计算，这之间涉及到大量的数据交换与处理问题，使得程序不可避免的变得庞大易错。考虑到上述问题，本文选择面向对象的c++语言在Visual Studio 2019集成开发环境下进行设计开发。将各层级的问题抽象为类进行组织，从而通过对不同类的实例化来实现有限元分析，提高程序的可重用性、易维护性和可拓展性。程序框架如图 1所示。

图1 程序框架

2.1 前处理模块

用户在通过前处理模块完成与有限元分析系统的交互，控制一系列参数完成几何模型的建立，同时确定材料特性、单元划分、荷载和边界条件、计算终止时间，脱空位置等数据，再进行整理并传递给分析计算模块。

2.2 分析计算模块

此模块为有限元计算的核心，接收前处理模块整理后的数据，实例化节点类，再以节点类为基础，结合材料参数实例化单元类。其中节点类储存了节点的节点号、位置、惯性、速度、加速度、节点力等信息，单元类储存了此单元成员节点的索引、单元号、单元集中质量矩阵、单元体积(长度)、自由度等信息。再将各实例化后的节点与单元输入系统集总类，形成整体集中质量矩阵，遍历节点，将集总得到的惯性赋予节点相应的自由度，同时进行加载与边界条件的处理。最后再引入运算类进行迭代计算，其中包含节点内力计算、接触搜索、接触力计算与沙漏力计算等过程。分析计算模块程序执行流程如图2所示。

图2 分析计算模块程序执行流程

2.3 后处理模块

TECPLOT是一款功能强大的数据分析和可视化处理的软件，它提供了丰富的绘图格式，能后根据有限元分析输出的节点信息生成网格，支持云图绘制等有限元后处理过程中不可缺少的功能。

故本文使用TECPLOT软件进行可视化后处理，后处理模块整理计算结果，根据用户需求提取不同时间步结果，输出为结构化数据导入TECPLOT软件。

3 程序验证

3.1 算例建立

本文依据CRTS Ⅰ 型板式无砟轨道结构参数建立有限元模型作为算例，与通用显式有限元软件ANSYS/LS-DYNA进行对比验证。

在ANSYS/LS-DYNA软件中，用Beam161单元划分轨道，对应本文程序中的大转动空间梁单元；用Solid164单元划分轨道板、砂浆与混凝土底座板，对应于本文程序中的减缩六面体实体单元；用Spring-damping165单元模拟地基与扣件，对应于本文程序中的弹簧单元。

轨道采用CHN60钢轨；扣件间距为0.625 m，刚度取50 MN/m；轨道板长度为5 m、宽度为2.4 m、厚度0.19 m；轨道板选用C50混凝土，弹性模量取36 000 MPa；砂浆充填层的长宽与轨道板相同，厚度为0.03 m，弹性模量取300 MPa；混凝土底座板宽3 m，厚0.3 m，选用C40混凝土，弹性模量取32 500 MPa；路基基础面刚度取190 MPa/m。

为了避免边界条件对结果的影响，建立3块轨道板长度的模型，取中间块作为分析和研究的对象，模型如图3所示。同时为了更好的模拟出CA砂浆离缝脱空状态下轨道的力学行为，选择最不利位置进行加载。据参考文献[12]可知，当荷载加载于扣件附近的钢轨之上时，结构的受力最明显，故本文在模型纵向方向8 500 mm位置设置长1 500 mm，宽2 400 mm，高0.1 mm的离缝，再选择在离缝范围内中间位置的扣件上方作为最不利位置，按照单轴双轮荷载施加集中力P，示意见图4，其中P取150 kN，以随时间线性增大的方式加载于处于最不利位置的轨道上，以达到拟静态的效果。

图3 轨道模型

图4 荷载示意

3.2 结果分析

提取ANSYS/LS-DYNA与本文程序结果进行对比研究，发现：

(1) 考虑到有限元计算过程包含接触分析，故提取由于离缝脱空而与砂浆发生接触行为的轨道板部分计算结果进行研究分析。通过对比轨道板部分最大垂向扰度关于时间的发展曲线,如图5所示，ANSYS/LS-DYNA与本文程序结果符合良好，发展趋势一致，最大相对误差为2.6%，发生在0.19 s时，参考图6轨道板部分最大垂向加速度关于时间的发展曲线，可以看到本文程序与ANSYS/LS-DYNA计算结果的相对误差同样也是在0.19 s附近波动较大，而这正是发生接触的时刻，故推测误差可能是因为接触算法部分的区别造成的。

(2) 提取1.0 s时刻模型整体的垂向扰度计算结果绘制云图(图7为ANSYS/LS-DYNA结果,图8为本文程序结果)，可以看到垂向扰度分布规律一致，最大垂向扰度都出现在同一位置，位于加载位置扣件的下方。

图5 轨道板垂向位移

图6 轨道板垂向加速度

图7 ANSYS/LS-DYNA垂向位移云图

图8 程序垂向位移云图

4 结束语

本文基于显式有限元的基本思想，针对有脱空轨道板力学分析时难以处理的高度非线性接触问题，采用面向对象的c++语言设计开发了一套完整的有限元程序，并通过与通用商业有限元软件ANSYS/LS-DYNA的对比验证了程序的正确性与、实用性与高效性，为有脱空轨道板力学分析计算提供了参考，为未来更进一步的显式有限元软件开发工作打下基础。