高速铁路双块式无砟轨道计算模型约束方程的建立

卿启湘^{1, 2}，胡萍²，王永和²，尹汉锋¹， 张春顺²

 (1. 湖南大学 汽车车身先进制造技术国家重点实验室 航天技术研究所，湖南 长沙，410082；

2. 中南大学 土木建筑学院，湖南 长沙，410075)

摘  要：为合理处理高速铁路双块式无砟轨道结构动力学分析中梁-板与弹簧、梁-杆与板壳、板壳与实体等连接问题，基于Timoshenko板壳理论和直接法建立多点约束方程，构造双块式无砟轨道组合结构计算模型，将组合结构中相关节点的位移通过偏心关系与板壳理论的直线假设建立的多点约束方程引入Galerkin法的弱积分形式中，解决各种不同类型单元因自由度不同和偏移连接而对组合结构总体刚度矩阵的贡献不同等问题。从理论上分析各种不同类型单元相结合的多点约束方程构建问题，并在此基础上，建立高速铁路双块式无砟轨道结构动力学有限元模型。结果表明：基于Timoshenko板壳理论和直接法所建立的有限元模型是合理的，双块式无砟轨道不同结构层连接良好，能消除常规方法产生的不合理附加应力。

关键词：双块式无砟轨道；有限元；约束方程；直接引入法

中图分类号：TU416          文献标志码：A         文章编号：1672-7207(2010)04-1360-09

Constraints equations establishment for calculation model of double-block ballastless track in high speed railway

QING Qi-xiang^{1, 2}, HU Ping², WANG Yong-he², YIN Han-feng¹, ZHANG Chun-shun²

(1. State Key Laboratory of Advanced Design and Manufacture for Vehicle Body, Institute of Space Technology, Hunan University, Changsha 410082, China;

2. School of Civil Engineering and Architecture, Central South University, Changsha 410075, China)

Abstract: In order to appropriately deal with connection problems in kinetics analysis such as shell-beam-bar, superposition plates, beam axis-spring, shell-spring, shell-solid connections in high speed railway’s slab track, Timoshenko-type shell theory and direct introduction method were based for establishing multi-point constraint equations and calculation model of the composite slab track. In the composite structure, to counteract different type elements’ attribution to the composite structure’s total rigid matrix caused by the elements’ offset connections, the multi-point constraint equations on the basis of both the linear hypothesis in shell theory and eccentricity relations among relative nodes pairs in the composite structure were introduced to equilibrium equation. Multi-point constraint equations of different type elements were constructed theoretically, and followed by which the dynamic FEM model of the high speed railway’s slab track was thereby established. The reliability of the FEM model based on Timoshenko-type shell theory and direct introduction method is reliable by the comparison between numerical analysis and real measured results, and meanwhile the model produces good connections among different layers of the structure, eliminating irrational additional stress generated in conventional methods.

Key words: double-block ballastless track; FEM; constraint equation; direct introduction method

在工程实际问题中，经常遇到由连续体和薄壁板壳以及细长杆、梁等组成的结构。在数值计算中，利用实体单元分析杆件和板壳结构可以避免引入结构力学的简化。但是这样做在实际分析中将遇到困难，这是由于在用实体单元对结构进行离散时，若单元的2个方向或1个方向比其他方向小得多，则使得单元不同方向的刚度系数相差过大，从而导致求解方程出现病害或奇异，最后使解精度很低或失败。为避免上述问题，作者以双块式无砟轨道路基系统为例，对不同的结构引入不同的单元，如扣件及轨下胶垫系统采用弹簧-阻尼单元，而与它连接区域钢轨或双块式轨枕采用梁单元；整体混凝土道床采用板壳单元，而与它连接区域支承层采用实体单元。但是，这些单元节点的配制不可能一致，有的还产生位置偏移，为保证交界面上的位移协调，就要研究和解决不同类型单元的连接问题。为保证交界面上的位移协调性，基于弹性力学理论的实体单元对位移场函数未提出进一步假设，而基于板壳理论的板壳单元或梁单元则假设中面的法线变形后仍保持为直线，且忽略其长度的变化。对于板壳单元或梁单元与实体单元相连接的情况，恰当的提法是要求变形后实体单元与板壳单元或梁单元的交界面上各点位移沿垂直于该面方向(即平行于壳体或梁中面的方向)的分量相同，即两类单元在该面上仍保持贴合，但不能要求两类单元交界面上各点位移沿平行于该面方向(即垂直于壳体中面的方向)的分量应相同，因为板壳理论中已忽略该方向位移的变化。如果要求两类单元在该方向位移一致，那么就将板壳理论中的假设也引入了实体单元，也就是强迫实体单元在交界面上该方向的应变为0。这个不恰当的强制条件将会使实体单元在交界面附近产生相当大的附加应力，从而使仿真计算结果失真。Surana^[1-2]提出用于轴对称分析和三维应力分析的过渡单元以解决不同类型单元的连接问题，还有很多研究者对壳元与梁-杆单元的连接问题也建立了约束方程^[3-6]。本文作者首先采用直接引入法，建立不同类型单元之间的约束方程，然后，将约束方程引入Galerkin法的弱积分形式中，经转换后的实体部分的刚度矩阵和载荷向量可按通常的步骤进行，将板壳或梁部分的刚度矩阵或载荷向量集合成系统的刚度矩阵和载荷向量，从而得到系统的求解方程组。

1  多点约束方程的建立

双块式无砟轨道由下至上有3个界面：第1界面为垫层与下部支承结构；第2界面为混凝土灌筑层与垫层；第3界面为钢轨与轨枕，用扣件相联结。

1.1  板壳单元与梁-杆单元连接

在双块式无砟轨道中，双块式轨枕嵌入到道床板中，如图1所示。现以板壳的节点I和梁单元节点i为例，讨论建立两者交界面上的位移协调条件。在总体坐标系内节点I的位移参数是u_I, v_I, w_I, 和   ；梁单元节点i的位移参数是u_i, v_i, w_i, , 和。为了建立两类单元间的位移约束条件，首先将总体坐标系内的节点位移参数转换到局部坐标系。这里假设板壳单元在节点I与梁单元在节点i的局部坐标系方向是相同的，所以，可以统一采用梁单元节点i的局部坐标系x′, y′, z′。x′沿梁单元中心线，y′轴和z′轴分别沿截面的横向对称轴和纵向对称轴。经过坐标转换后，可以得到梁单元节点i在局部坐标系内的位移参数, , 和, , ；板壳单元节点I在局部坐标系内的位移参数, , 和, ，但不能得到。这是因为原来板壳单元节点I的位移参数和是绕与法线(即z′ 轴)相垂直的2个轴的转动，所以，经过坐标转换后只能得到和，而不能得到。为了进一步建立位移协调条件的需要，

可从中面内过节点I沿边界的切线的转动得到，即

                  (1)

实际上考虑到是可以用板壳单元的节点位移参数表示，即

               (2)

式中：分别为板壳单元节点I+1与I-1之间局部坐标系轴向距离。

图1  梁轴与板中面之间具有偏移的连接

Fig.1  Offset connection between beam axis and middle surface of shell element

关于板壳单元与梁单元在交界面上的位移和变形条件可以表述为：(1) 变形前板壳单元上通过节点I的中面法线与梁单元上通过节点i的纵向对称轴在同一条直线上，且板壳单元法线上的I_t点与梁单元纵向对称轴上的i_t点相连接；变形后，该直线仍保持为直线，即连线I-I_t-i_b-i仍是一条直线，且点I_t与点i_b仍保持连接。(2) 变形前板壳单元的包含该法线在内的边界横向截面的切平面与梁单元通过节点i的截平面是在同一个平面内；变形后它们仍保持在同一个平面内。但板壳单元的底面与梁单元的顶面在梁单元的宽度方向(即局部坐标系的方向)上，2个面应允许有一定的滑动。这是因为梁的理论中假设宽度b在变形过程中保持不变，而板壳单元底面在其自身平面内是可以变形的。若要求板壳单元底面与梁单元顶面之间在方向不能滑动，则给板壳单元施加了不合理的约束，将会产生很大的误差。因此，两类单元在交界面上的位移协调条件应是在局部坐标系, 和中，板壳单元的I_t点与梁单元的i_b点的3个位移分量相等，以及板壳节点I与梁单元节点i的3个转动分量相等，即

                 (3a)

                 (3b)

                (3c)

                (3d)

                (3e)

                (3f)

为了进一步利用板壳单元与梁单元的节点位移参数表示上述的位移约束条件，根据板壳单元与梁单元的几何方面的假设，可以表达为：

              (4a)

              (4b)

                 (4c)

              (4d)

              (4e)

                  (4f)

将式(4a)~(4f)代入式(3a)~(3f)，并考虑式(1)，就可以得到局部坐标系内板壳单元与梁单元节点位移参数间的约束方程，即：

       (5)

式中：H为梁厚度的1/2。

因为刚度集合以及求解方程是在整体坐标系中进行的，所以必须将节点位移参数转回到总体坐标系，这是通过坐标转换矩阵来实现的。设

            (6)

式(6)是两坐标系的关系矩阵，其中， 是与ox轴夹角的余弦，其余类推。由

              (7)

得到12×12阶的坐标转换矩阵，有

               (8)

因而最后可以建立整体坐标系内的板壳单元与梁单元节点之间的约束方程^[3]。

1.2  梁轴与弹簧之间具有偏移的连接

当车辆在轨道上行驶时，机车车辆的垂直轮载不是作用在钢轨中心线上，横向力也不是作用在钢轨截面的剪力中心上，因此，钢轨上作用有扭矩。在车辆正常运行情况下，钢轨的垂向、横向位移和转角相互之间无影响，可以线性叠加，因此，钢轨受力情况如图2所示。对于轨下结构的动力响应可以考虑上述轮载力的作用^[7]。这样，轨道与轨道板之间关系的集成模型如图3所示。

图2  钢轨受力的组成

Fig.2  Superposition of three forces on track

图3 双块式轨道-柔性基础的集成模型

Fig.3  Coupling model of double-block ballastless track-flex foundation

对于高速铁路双块式无砟轨道结构，钢轨的正下方为扣件和轨下胶垫，依实际受力情况来看，可假设钢轨与扣件和轨下胶垫的连接为铰接，钢轨可以转动，且其相接触处位移相等，如图4所示。

若扣件和轨下胶垫采用弹簧-阻尼单元模拟，钢轨采用梁单元模拟，并且梁单元可以自由转动，则根据上述分析，变形前梁单元上通过节点I的中面法线与弹簧单元上通过节点i的纵向对称轴z在同一条直线上，且梁单元法线上的I_b点与弹簧单元纵向对称轴上的i点相连接；变形后该直线仍保持为直线，即连线I-I_b-i仍是一条直线，且点I_b与点i仍保持连接。这样，当弹簧单元位于梁单元下方时，其多点约束方程为：

或       (9a)

或       (9b)

或           (9c)

而当弹簧单元位于梁单元上方时，其多点约束方程为：

或      (10a)

或      (10b)

或          (10c)

把式(9a)~(9c)和(10a)~(10c)写成矩阵形式，有

      (11)

图4  梁轴与弹簧之间具有偏移的连接

Fig.4  Offset connection between beam axis and spring element

1.3  板单元与三维实体单元的连接

为了增强双块式无砟轨道结构的弹性，在道床板下面还铺有一层水凝性混凝土稳定承载层。依实际工况，可假设道床板与混凝土稳定承载层的连接为刚接，且其相接触处位移相等。根据前面分析，基于Timoshenko梁理论假设，若要求板单元节点(I, J)位移很好地与实体单元节点(i, j)相衔接，则其对应节点位移之间的关系可由刚性杆相连接来模拟，如图5所示。其约束方程为：

          (12a)

           (12b)

               (12c)

因其连接为刚接，则和还应满足下列关系式

           (13a)

           (13b)

式中：和分别为实体单元节点(i+1, j)与(i-1, j)和(i, j+1)与(i, j-1)之间的整体坐标系x和y坐标之差。这样，由式(11)可得：

 (14a)

    (14b)

              (14c)

把上式写成矩阵形式，即

                 (15)

上面讨论了梁单元与实体单元、梁单元与弹簧单元以及板壳单元与梁-杆单元、板壳单元与板壳单元的连接问题。前者是板壳单元在其顶面(或底面)上与实体单元相连接，而后者则是板壳单元的顶面(或底面)与梁单元的底面(或顶面)相连接。研究不同类型单元连接的首要问题是正确地建立它们交界面上的约束方程。不能将仅适用于一种单元的力学简化强加于与其相连接的另一种单元，即不能将板壳单元中忽略法线长度变化的假设强加于实体单元，也不能将梁单元中截面形状不变的假设强加于板壳单元，否则将会由于过分的约束而导致相当大的误差。其次是引入约束方程的方法，罚函数近似是方法，而直接引入法是精确的方法。

至此，本文已经采用直接引入法，建立了各种不同类型单元之间的约束方程。最后就是将约束方程引入Galerkin法的弱积分形式中，形成刚度矩阵，再利用数值计算方法求解^[8-12]。

图5  板单元与实体单元的连接关系

Fig.5  Shell-solid elements contact relation

2  车辆载荷表述及材料参数

2.1  车辆载荷表述

在任一时刻t，设车辆匀速移动轮载为^[13-15]：

式中：p为轮载；x为轮载在t时刻的位置；δ为δ函数，即δ(x_i=vt_i)=1，δ(x_i≠vt)=0；x_i为双块式轨枕所在位置，x_i=；x₀为初始参考点，为双块式轨枕间距。选用我国高速铁路实际使用的动力分散式高速列车，动力车定距为11.4 m，轴距为3.0 m，车长为21.29 m，动力车的平均轴质量为19 435.22 kg，相应轮载力为125 N。拖车定距为18.0 m，平均轴质量为12 146.34 kg，相应轮载力为78 N。列车经过两相邻承轨槽所需的时间。其施加轮载算法见图6。

2.2  计算模型及材料参数

计算时，列车编组共取6辆车，车与车之间距离为0.95 m，计算时间等于第1个轮载进入模型至第24个轮载驶出模型所用时间加上列车驶出模型后所设定的延长时间，分析路基响应情况。

钢轨选用日本60轨，横截面为“工”字形梁；轨距为1 435 mm。扣件和轨下胶垫系统视为弹簧-阻尼结构，初始长度为38 mm。取其总等效刚度k=4.38×10⁷ N/m，阻尼系数c=4.5×10⁵ N·s/m。基床底层及路堤填高为3.0~5.1 m。其他计算参数见表1。

位移边界条件：模型左右两端面位移U_x=0；两侧面位移U_y=0；底面固定深度为40 m。

3  过渡段路基的动态响应特征

某客运专线有一双块式无砟轨道的隧-隧过渡段，涵洞为一框架涵，宽3.6 m，净高3.5 m，离路基面高5.1 m，涵洞中心线距离计算模型起始点(一端隧道口)19.8 m，距离计算模型终止点(另一端隧道口)26.4 m，计算长度共52.2 m，其中，过渡段长46.2 m，隧道内各取长3 m。基床表层0.4 m厚采用级配碎石掺加5%水泥填筑，隧道内挖除换填2.3 m厚，基坑回填C25混凝土(路堑地段基底挖除换填2.3 m厚的级配碎石掺5%水泥)；路堤地段基床表层以下采用级配碎石掺加5%水泥填筑，两侧为级配碎石；基底强风化层以上土层及全风化层全部挖除换填C15片石混凝土。考虑到对称性，本文只取路基一半作为研究对象，车速为350 km/h。有限元模型见图7。

图6  施加轮载算法

Fig.6  Wheel loads arithmetic

表1  计算所用参数

Table 1  Parameters used for calculation

图7  计算模型

Fig.7  Calculational model

3.1  动响应随时间的变化

图8所示为列车经过距离两隧道口各5.4 m时轨道面的动位移时程的变化，其中波峰为各轮载经过过渡段某横断面时轨面上的动位移随时间的变化，波谷为无轮载作用时轨面上的动位移的变化。这说明1个轮载对路基的作用分加载和卸载2个过程，对于路基结构，1个转向架就相当于2次加卸载。转向架的前轴通过后，动态响应降低至一定程度便开始增加，这是转向架前后轴轮载对路基作用的叠加结果。前轮经过时，动态响应小，后轮经过时，动态响应大，这与观察到的现象一致^{[7, 16-18]}。动力车产生的轮载力大，对路基影响也大；拖车产生的轮载力小，对路基影响也小。

1—离武汉向隧道口5.4 m; 2—离广州向隧道口5.4 m

图8  离两隧道口各5.4 m的轨道面动位移时程曲线

Fig.8  Track’s dynamic displacement changing with time for 5.4 m distant from two tunnels abutments

图9所示为列车经过离两隧道口各5.4 m时轨道正下方基床面的动应力时程的变化。由图9可见：基床面的竖向应力随时间的变化规律与动位移的变化规律非常类似。

1—离武汉向隧道口5.4 m的基床面上动应力;
2—离广州向隧道口5.4 m的基床面上动应力

图9  离两隧道口各5.4 m的基床面动应力时程曲线

Fig.9  Roadbed’s dynamic stress-time curve for 5.4 m distant from two tunnels abutments

图10所示为列车经过隧-隧过渡段最中央的轨道正下方离轨道面1.0和3.3 m深度处动应力时程的变化。由图10可见：动应力随时间的变化规律也与动位移的变化规律非常类似，离轨道面深度越浅，动响应越大。

1—离轨道面垂向距离1.0 m;
2—离轨道面垂向距离3.3 m

图10  隧-隧过渡段轨道正下方不同深度处动应力时程曲线

Fig.10  Roadbed’s dynamic stress-time curve at midpoint of two tunnels abutments at different depths below track

3.2  动响应随路程的变化

    图11所示为钢轨型面上的动位移仿真计算结果(有约束与无约束方程)，图12所示为轨道正下方路堤面上的动应力仿真计算结果(有约束与无约束方程)与实测值对比。从图11和图12可见：一方面，不考虑约束方程的仿真计算位移达85 mm(钢轨型面上)，而轨道以下结构的动态响应应力始终为0，这显然与实际情况不符，这证明无砟轨道路基结构的有限元模型若不考虑不同性质单元间的约束是完全错误的；另一方面，仿真计算值与实测值基本吻合，由此证明了有限元模型考虑不同性质单元间的连接约束是必需的。

1—考虑约束方程；2—不考虑约束方程

图11  钢轨型面上的动位移变化

Fig.11  Dynamic displacements-distance curves at surface of rail

1—考虑约束方程的计算值；2—实测值；
3—不考虑约束方程的计算值

图12  轨道正下方路堤面的动应力

Fig.12  Dynamic stress-distance curve at roadbed below track

4  结论

(1) 基于板壳理论所建立的板与梁接触面上点铰接的多点约束方程，保证了模型中板与梁的良好连接，使其具有较高的数值计算精度。

(2) 基于板壳理论所建立的梁与弹簧交界面上点铰接的多点约束方程，保持了弹簧单元分析的独立性，易于编程实现。

(3) 采用板壳理论所建立的壳元与体叠加面上点刚接的多点约束方程，不依赖于实体与壳连接的角度和位置，使其对实体与壳的叠加连接形式有较强的适应能力。

(4) 建立无砟轨道路基结构的有限元模型，必须考虑不同性质单元间的约束方程，否则，计算结果是错误的，这为类似的工程问题解决提供了一种较为精确的计算方案。

参考文献：

[1] Surana K S. Transition finite elements for axi-symmetric stress analysis[J]. Int J Num Meth Eng, 1980, 15: 809-832.

[2] Surana K S. Transition finite elements for three-dimensional stress analysis[J]. Int J Num Meth Eng, 1980, 15: 991-1020

[3] 王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003: 23-25.

WANG Xun-cheng. Finite element method[M]. Beijing: Tsinghua University Press, 2003: 23-25.

[4] Makoto S. Vibration behavior of roadbed on soft grounds under trainload[J]. Quartly Report of RTRI, 1990, 42(31): 29-35

[5] 卿启湘. 高速铁路无碴轨道-软岩路基系统动力特性研究[D]. 长沙: 中南大学土木学院, 2005: 125-150.

QING Qi-xiang. Study on dynamical characteristics in soft zone for ballastless tracks on high speed railway[D]. Changsha: Central South University. School of Civil Engineering and Architecture, 2005: 125-150.

[6] 文明才, 戴振东. 冲击荷载作用下Winkler地基上横观各向同性浅球壳的动力分析[J]. 中南大学学报: 自然科学版, 2008, 39(3): 629-634.

WEN Ming-cai, DAI Zhen-dong. Dynamic analysis of transversely isotropic shallow spherical shells of Winkler foundation subjected to impact force[J]. Journal of Central South University: Science and Technology, 2008, 39(3): 629-634.

[7] 翟婉明. 车辆-轨道棍合动力学[M]. 北京: 中国铁道出版社, 2002: 56-75.

ZHAI Wan-ming. Vehicle-track coupling dynamics[M]. Beijing: China Railway Press, 2002: 56-75.

[8] Tan L H, Bathe K J. Studies of finite element procedures—The conjugate gradient and GMRES methods in ADINA and ADINA-F[J]. Computers & Structures, 1991, 40: 440-449.

[9] LIU W K, Belytschko T, Zhang Y F. Implementation and accuracy of mixed-time implicit-explicit methods for structural dynamics[J]. Computers & Structures, 1984, 19: 521-530.

[10] HE Xu-hui, YU Zhi-wu, CHEN Zheng-qing. Finite element model updating of existing steel bridge based on structural health monitoring[J]. Journal of Central South University of Technology, 2008, 15(3): 399-403.

[11] Wills J R. Hertz contact of anisotropic bodies[J]. Journal of Mechanics and Physics of Solids, 1966, 14(2): 163-176.

[12] Chung L, Hulbert G. A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized α-method[J]. J Appl Mech, 1993, 60(4): 371-375.

[13] Wilson E L, Yuan M W, Dickens J M. Dynamic analysis by direct superposition of Ritz vectors[J]. Earth Eng & Struct Dynamic, 1982, 10: 813-821.

[14] Li D Q, Sellin E T. Wheel-track dynamic interaction: Track substructure perspective[J]. Vehicle System Dynamics Supplement, 1995, 24: 183-196.

[15] 梁波, 罗红, 孙常新. 高速铁路振动荷载的模拟研究[J]. 铁道学报, 2006, 28(4): 89-94.
LIANG Bo, LUO Hong, SUN Chang-xin. Simulated study on vibration load of high speed railway[J]. Journal of the China Railway Society, 2006, 28(4): 89-94.

[16] Monley G J, Wu T H. Tensile reinforcement effects on bridge-approach settlement[J]. Journal of Geotechnical Engineering, 1993, 119(4): 7749-7762.

[17] 练松良, 王继军, 杨文忠. 路桥过渡段轨道结构的动力性能试验研究[J]. 中国铁道科学, 2009, 30(4): 19-24.
LIAN Song-liang, WANG Ji-jun, YANG Wen-zhong. Test study on the dynamic performance of the track structure at the transition section between bridge and plain track[J]. China Railway Science, 2009, 30(4): 19-24.

[18] 雷晓燕, 张斌, 刘庆杰. 过渡段动力特性有限元分析[J]. 中国铁道科学, 2009, 30(5): 15-21.
LEI Xiao-yan, ZHANG Bin, LIU Qing-jie. Finite element analysis on the dynamic characteristics of the track transition[J]. China Railway Science, 2009, 30(5): 15-21.

收稿日期：2009-08-17；修回日期：2009-12-02

基金项目：国家自然科学基金资助项目(50678177)

通信作者：卿启湘(1964-)，男，湖南隆回人，博士，副教授，从事铁道工程与CAD设计研究；电话：0731-88829567；E-mail: qixiangcn@163.com

(编辑 陈爱华)