舒 瑶，蒋忠城，张 俊，张 波，刘国云，杨新文

(1. 西南交通大学牵引动力国家重点实验室，四川，成都 610031；2. 大功率交流传动电力机车系统集成国家重点实验室，湖南，株洲 412001；3. 中车株洲电力机车有限公司，湖南，株洲 412001；4. 同济大学道路与交通工程教育部重点实验室，上海 201804)

无砟轨道结构损伤将会降低结构的承载能力，恶化轮轨关系，威胁运营安全，且轨道结构一旦出现伤损病害，不仅维修困难，而且维修费用昂贵。无砟轨道基础结构主要为多孔性、多物相、极限应变小的水泥基准脆性材料，在车辆动荷载重复作用下，结构逐渐损伤劣化甚至断裂失效。轨道结构损伤缺陷的动力演化及长时劣化，必然引起结构服役状态与动态性能的持续改变[1 − 3]，导致轮轨关系恶化，引起列车振动加剧，反过来又使得车辆系统对轨道的动力破坏作用增大，轮轨关系恶化与轨道结构损伤是一个长期的恶性循环过程。

国内外许多学者都曾围绕混凝土结构损伤及其相关问题进行了研究[4 − 12]，无砟轨道结构为典型的混凝土结构，为了研究无砟轨道结构的动力损伤问题，许多学者将混凝土结构分析理论引入无砟轨道伤损病害问题的研究，使无砟轨道结构的耐久性及伤损病害分析也大量采用数值计算同现代破坏力学理论相结合的分析方法。断裂力学[3, 13 − 14]、损伤力学[15 − 16]、疲劳分析[16 − 17]等现代破坏力学理论，以及界面力学[18 − 19]理论和方法都被用于含损伤缺陷轨道的精细化分析。该类模型不再局限于线性的材料本构关系，而是引入断裂力学、损伤力学等混凝土非线性力学行为的分析方法，建立精细化的数值模型，分析无砟轨道结构或其局部力学性能的变化。利用精细化数值模型进行动力仿真分析时，计算十分费时，且该类模型为了消除数值模型人工截断边界的影响，而将实际模型范围延拓，往往导致计算量巨大。为深入、细致地研究轨道结构伤损的产生与动态演化，而不得不追本溯源地考虑轮轨之间的动态相互作用。现有可以考虑轨道结构损伤缺陷发展演化的精细化数值模型因其计算量过大，难以实现车-轨耦合动力学分析，这也是目前很少见到同时考虑材料非线性及接触非线性的车辆-轨道耦合动力分析的原因之一。

笔者曾在文献[20 − 21]中提出了一种能考虑轨道结构损伤发展的非线性轨道动力学模型，并利用拟谱法(Pseudo-spectral method，PSM)实现了模型的高效快速求解。为了探究无砟轨道结构损伤与车轨动态相互作用的相互影响机制，本文进一步基于该考虑轨道结构损伤的动力学模型，将车辆简化为多刚体系统，且假定随机性的轨面不平顺以单节车长为周期重现，利用Hertz 非线性接触实现车辆系统与含损伤无砟轨道系统的垂向传力耦合，从而建立含损伤无砟轨道-车辆垂向耦合周期性动力模型；为了加速该同时包含材料非线性和接触非线性的动力模型求解的收敛速度，采用隐式动力预测-校正算法和轨道-车辆系统交叉迭代的求解策略，实现了含损伤无砟轨道-车辆垂向耦合动力模型的隐式快速求解。

1 模型及运动方程

图1 所示为含损伤无砟轨道-车辆垂向耦合系统，包含上部车辆子系统和下部含损伤无砟轨道子系统。将车辆简化为包括车体、转向架、轮对在内的多刚体系统，具有车体沉浮和点头运动、前后转向架沉浮和点头运动、4 个轮对垂向运动等10 个自由度[22 − 23]，单节高速车辆运动方程为：

图1 含损伤无砟轨道-车辆垂向耦合系统Fig.1 Vertical coupling system of ballastless track and vehicle with damage

层状的无砟轨道体系在竖向上具有极大的材料不均匀性和差异性，自上而下的钢轨、轨道板、垫层和支承层的变形性能和损伤破坏模式差异巨大，在线弹性范围内评估轨道结构的动力响应，忽视了不同结构材料的损伤特性和及其承载能力的差异，因而难以精确表征轨道系统的动力行为特性。无砟轨道的轨道板、支承层均为受拉易开裂的水泥基准脆性材料，研究表明，该类材料基本上都是带缺陷或损伤服役，而且在列车荷载、温度力等荷载或环境因素长期作用下，无砟轨道结构的损伤会不断发展而影响其服役性能。真实材料的变形往往既有刚度退化，又有塑性变形存在，无砟轨道混凝土基础结构在列车动荷载循环加卸载作用下动应变不大，一般忽略塑性变形，仅考虑弹性应变，轨道结构在服役过程中由于损伤累积而导致刚度降低的现象，则利用混凝土弹性损伤本构来描述。故从宏观唯象的角度考虑混凝土受拉损伤的影响，采用《混凝土结构设计规范》(GB 50010−2010，2015 年版)[24]推荐的混凝土单轴受拉的弹性损伤本构模型，该弹性损伤本构的数学关系式可以表达为：

该弹性损伤本构应力-损伤-应变曲线如图2 所示，无砟轨道轨道板和支承层的初始损伤场往往难以给定，本文采用的是方法是在荷载作用下，利用式(2)计算生成，轨道板C55 损伤因子d0为0.175 891，支承层C15 损伤因子d0为0.062 208。

图2 弹性损伤本构应力-损伤-应变曲线Fig.2 Stress-damage-strain curves for elastic damage constitutive

如图1 所示，截取周期长度为车长LC的一段无砟轨道结构，钢轨、轨道板及支承层均采用欧拉梁模拟。钢轨、轨道板和支承层的振动方程[20 − 21]分别见式(3)～式(5)：

车辆运动方程式(1)与轨道系统振动方程通过轮轨接触关系互相耦合，轮轨间Hertz 非线性接触轮轨力pj(t)为[22 − 23, 25]：

式中：相对位移δZ=Zwj(x,t)−[zr(xpj,t)+zr0(xpj,t)]，轮轨接触常数G=3.86R−0.115×10−8m·N−2/3(磨耗型踏面车轮)；Zwj(x,t)为t 时刻第j 位轮对的垂向位移；zr(xpj,t)为t 时刻第j 位轮对下钢轨的垂向位移；zr0(xpj,t)为t 时刻第j 位轮对下轨道的不平顺。

如图3 所示，假定随机性的轨面不平顺以单节车长为周期重现，进行车辆运动仿真计算时，先假设Car A 完全占据单个周期长度的轨道节段，如图3(a)所示，随着列车向前行进，邻车Car B也逐渐进入该轨道节段，如图3(b)～图3(c)所示，直至邻车Car B 完全占据该轨道节段，图3(a)～图3(d)恰好为单节车辆通过周期轨道节段的全过程。仿真过程中涉及到2 节车辆系统的平衡，利用模型的周期性进行简化，当Car A 超出周期性轨道模型的边界时，将邻车Car B 前转向架下的车轮所受到的轮轨力作为Car A 前转向架下的车轮所受到的轮轨力，轨道系统的车辆荷载也需进行对应处理。

图3 含损伤无砟轨道-车辆系统耦合仿真计算示意图Fig.3 Schematic diagram of coupling simulation of ballastless track vehicle system with damage

2 模型求解

2.1 含损伤无砟轨道振动方程空间域离散

文献[20 − 21]应用周期性谱方法缩减模型求解规模，利用基函数的周期性处理模型周期性的边界条件，并给出了含损伤无砟轨道振动方程的离散过程。谱方法(Spectral method)起源于Ritz-Galerkin 方法，以逼近性质良好、且便于利用快速变换来计算的正交多项式(Fourier 多项式、Chebyshev多项式、Legendre 多项式等)作为基函数，又分为Galerkin 谱方法、Tau 方法或拟谱方法(Pseudospectral method，又称谱配点法)。谱方法利用正交基函数系将待求函数展开作为其某种意义上的逼近函数，然后通过人为截断，把无限维问题简化为有限维问题，将截断的逼近函数展式代入原偏微分方程，利用Galerkin 方法将偏微分方程转化为常微分方程，而后求解常微分方程的得到逼近函数的系数，从而求得原方程解。

钢轨受到移动轮载作用，应用周期性Fourier-Galerkin 谱方法(简记为FGSM)比较方便，将钢轨位移写成Fourier 级数，利用Galerkin 法将钢轨振动偏微分方程离散成常微分方程组。

进而应用拟谱方法对轨道板和支承层的振动方程进行空间离散，令轨道板和支承层的振动方程在每一个计算配点满足加权残值为0，从而将轨道板和支承层的振动方程离散为常微分方程组。最后将离散后的常微分方程组整合为：

图4 计算配点和应力点示意图Fig.4 Schematic diagram of calculation collocation points and stress points

2.2 含损伤无砟轨道-车辆耦合动力模型的时间域求解

车辆-轨道耦合系统中，车辆系统被简化为由质量-弹簧-阻尼集成的线性刚柔体系，其运动方程的求解采用常规多刚体运动学求解方法即可；而轨道系统被简化为考虑材料损伤非线性的弹性地基上的叠合梁，其振动方程的求解往往需要采用增量迭代的方法，结构应力-应变的更新采用映射回退算法。而轮轨间Hertz 非线性接触进一步引入了车辆-轨道耦合系统运动方程边界条件的非线性，在每一时间步又必须采用迭代的方法求解。为了加速该同时包含材料非线性和接触非线性的动力模型求解的收敛速度，采用隐式动力预测-校正算法和轨道-车辆系统交叉迭代的求解策略，以期实现含损伤无砟轨道-车辆垂向耦合动力模型的隐式快速求解。

2.2.1 隐式动力预测-校正算法

为了加速考虑损伤效应的无砟轨道振动方程的求解，在隐式求解方法[21]的基础上，进一步引入预测-校正求解策略。动力方程的增量形式[21]可写为：

由预估值式(10)计算得到不平衡的残余力后，继续利用残余力计算得到广义位移增量，并利用其对预估值进行校正，如此重复操作，直至达到前后两个迭代子步广义位移的计算结果之差满足预设容差，此时认为迭代收敛而转入下一时间步继续计算。t+∆t时步内迭代计算过程概要如下：

① 令迭代计数变量i=1，开始预估阶段，令：

③ 用下式计算等效广义刚度矩阵 K∗：

⑤ 进入修正阶段：

⑥ 如果第i 迭代步不是初始迭代步，判断第i 步与第i−1 步两个迭代子步轨道系统动位移的计算结果之差是否满足预设容差，不满足则令i=i+1，并转到第②步，满足则转向⑦；

⑦ 利用式(16)得到的修正值作为t+∆t时刻轨道的最终计算结果，并进一步更新轨道的应变及损伤等其他物理量，而后开始下一时间步的计算。

2.2.2 轨道-车辆系统隐式动力交叉迭代求解

为便于根据方程特性针对性地选用求解方法，分别对车辆系统及轨道系统交叉求解，如图4所示，即在每一时间步，首先通过假设车辆系统、轨道系统的动位移，分别提取车轮位移、轮轨接触点的轨面高低不平顺及钢轨位移，试算Hertz 非线性接触轮轨力，进而将其输入轨道系统求解轨道结构的动力响应而得到新的钢轨位移，更新Hertz非线性接触轮轨力后将其代入车辆系统求解车辆的动力响应而得到新的车轮位移，如此便完成了一个非线性接触的循环迭代过程，判断是否已经满足迭代收敛条件，不满足则再次更新轮轨力继续迭代计算，反之则进入下一时间步的计算。

严格来说，考虑材料损伤非线性的轨道系统运动方程的求解，包含两重意义的非线性问题的迭代，如图5 所示，分别是轨道系统的平衡迭代和本构关系的迭代，其动力响应的计算需采用增量迭代的方法，另外车辆-轨道耦合的Hertz 非线性接触也需要迭代，再加上时域求解固有的最外层时间步循环，故整个车辆-轨道耦合系统运动方程的时域数值求解的完整过程包含4 层循环，且其中3 层是非线性迭代。为简化求解过程，尝试将这三个层面的非线性迭代计算合并，不再需要在每一次试算轮轨力作用下的轨道系统都达到完全的迭代平衡，而是在Hertz 非线性接触的迭代计算中不断根据残余力修正轨道系统的位移，以期达到材料本构、轨道系统及车辆-轨道耦合系统三个层面非线性迭代的平衡。如此处理虽然简化了求解流程，但要达到迭代收敛往往会增加迭代步数。

图5 每一时步迭代计算示意图Fig.5 Schematic diagram of iterative calculation for each time step

为了加速迭代收敛，借鉴Newton-Raphson 法的思想，将隐式动力预测-校正算法嵌入轮轨非线性接触的循环迭代过程中，提出一种同时求解材料非线性和接触非线性的预测-校正迭代算法。在时刻t+∆t，车辆系统和轨道系统动力方程的典型形式分别为：

同理，t+∆t 时刻车辆系统的位移预估值为：

进而利用轨道、车辆系统的位移预估值，试算Hertz 非线性接触轮轨力，然后展开t+∆t时步各迭代子步的计算，直到达到前后两个迭代子步轨道系统动位移的计算结果之差满足预设容差，此时认为迭代平衡而转入下一时间步继续计算，t+∆t时步内迭代计算过程概要如下：

1) 令迭代计数变量i=1，开始预估阶段，令轨道系统的位移、速度和加速度为式(11)，令车辆系统的位移为式(19)；

3) 按照第2.2.1 节所述预测-校正算法中t+∆t 时步迭代计算过程②～⑤，计算轨道系统的位移，更新过程变量及Hertz 非线性接触轮轨力，进而将其输入车辆系统求解车辆的动力响应而得到新的车轮位移；

4) 如果第i 步不是初始迭代步，判断第i 步与第i−1 步两个迭代子步轨道系统动位移的计算结果之差是否满足预设容差，不满足则令i=i+1并转到第2)步，满足则转向第5)步；

5) 利用步骤4)得到的轮轨力生成轨道系统广义荷载向量，再按②～⑤修正轨道系统的广义位移、速度和加速度，以便得到与车辆系统位移同步更新的轨道系统位移，将其作为t+∆t时刻轨道系统的最终计算结果，并进一步更新轨道系统的应变及损伤等其他物理量，而后开始下一时间步的计算。

轨道-车辆系统隐式动力交叉迭代求解流程见图6。

3 模型验证与讨论

3.1 预测-校正算法的验证与讨论

为了验证预测-校正算法的正确性，将其计算结果与隐式动力N-R 法[21]的计算结果进行比较，时间步长Δt =1×10−4s。图7 给出了预校正算法与隐式动力N-R 法轨道板和支承层位移计算结果的对比，图7 中x 表示纵向坐标，t 表示时间，ΔZs表示轨道板垂向位移绝对差，ΔZb表示支承层垂向位移绝对差。由图可以看出，两种方法的轨道板和支承层位移计算结果的绝对差均非常小，可以认为两种方法的计算结果一致，预校正算法的位移计算结果准确可靠。

图6 交叉迭代求解示意图Fig.6 Schematic diagram of cross iteration solution

进一步提取预测-校正算法与隐式动力N-R 法轨道板和支承层最终损伤计算结果进行对比，如图8 所示，图中x 表示纵向坐标，t 表示时间，Δds表示计算的轨道板损伤因子绝对差，Δdb表示计算的支承层损伤因子绝对差。两种方法的轨道板和支承层最终损伤计算结果的绝对差也都非常小，可以认为两种方法最终损伤的计算结果一致，进一步验证了预校正算法计算结果的正确性，也相互应证了两种计算方法的可靠性。

3.2 含损伤无砟轨道-车辆垂向耦合动力模型的验证与讨论

以和谐号高速动车组CRH3 为例[25]，进行含损伤无砟轨道-车辆耦合动力分析，CRH3 高速动车组一般为8 节编组，最高运营速度达到350 km/h，仿真所用单节车辆系统参数见表1，轨道系统参数见表2。

图7 轨道板和支承层位移计算结果对比Fig.7 Comparison of displacement calculation results of slab and bed

轨面高低不平顺采用中国高速铁路无砟轨道轨面高低不平顺频谱[27]，计算所采用的高低不平顺样本见图9，假定随机性的轨面不平顺以单节车长为周期重现时，周期性模型偏向于考虑波长小于单节车长的不平顺成分，如需考虑波长成分更为丰富的轨面不平顺，可以将周期性模型的周期长度拓展为整数倍单节车辆长度。

为验证含损伤无砟轨道-车辆垂向耦合动力模型计算的正确性，特将无砟轨道系统退化到弹性无损阶段，采用FGSM 法[20 − 21]计算了钢轨为离散点支承的无砟轨道-车辆系统的动力响应，简记为FGSM-e，也采用PSM 法[20 − 21]计算了钢轨为连续均匀支承的弹性无损阶段无砟轨道-车辆系统的动力响应，简记为PSM-e；而将含损伤无砟轨道-车辆垂向耦合动力模型简记为PSM-d，三者均采用Hertz 非线性接触模拟车辆-轨道耦合系统的动力相互作用。进而分别提取轮轨力F、车体加速度及轨道结构位移等物理量进行对比，从而验证含损伤无砟轨道-车辆垂向耦合动力模型计算的正确性。图10 分别给出了FGSM-e、PSM-e 和PSM-d 三种模型计算结果的对比，图10 中x 表示纵向坐标，t 表示时间，F 表示轮轨力，由图可以看出，FGSM-e、PSM-e 和PSM-d 三种模型轮轨力F、车体加速度及轨道结构位移等物理量的计算结果基本一致，相互应证了各自模型计算结果的正确性。

图8 轨道板和支承层损伤计算结果对比Fig.8 Comparison of damage factors calculation results of slab and bed

表1 车辆系统参数Table 1 Parameters of vehicle system

表2 轨道系统参数Table 2 Parameters of track system

表3 给出了FGSM-e、PSM-e 和PSM-d 三种模型计算耗时对比，由表可知，对于弹性阶段的无砟轨道-车辆垂向耦合动力模型的求解，PSM 法求解速度最快，仅耗时67.1 s，而FGSM 法耗时100.8 s；考虑损伤效应后，利用PSM 法求解含损伤无砟轨道-车辆垂向耦合动力模型需耗时125.8 s，每一时间步长内最大迭代11 次收敛，整个仿真时间内每一时间步长平均迭代6 次，计算仍然非常高效，方便应用于多参数、多工况的仿真计算。

图10 计算结果对比Fig.10 Comparison of calculation results

表3 计算耗时对比Table 3 Calculation time comparison

图11 轨道板纵剖面损伤分布云图Fig.11 Damage distribution of longitudinal section of track slab

图12 轨道板上、下表面应变与损伤因子变化曲线Fig.12 Time history curves of strain and damage factor on slab surface

图11 给出了单节车辆通过后轨道板纵剖面损伤分布云图，图中x 表示纵向坐标，z 表示垂向坐标，以轨道板中心线为坐标0 位置，由该图可以看出损伤更多主要集中在轨道板的上、下表面。进一步提取轨道板跨中上、下表面应变与损伤因子，图12 给出了轨道板跨中上、下表面应变与损伤因子时程曲线，t 表示时间，由图可以看出，在车辆动荷载作用下，相对于初始的损伤因子0.175 891，含初始损伤的轨道板的损伤因子会有所增加，在当前的车辆动荷载水平下，轨道板在列车动载作用下弯曲时，如图12(b)所示，轨道板底部最大拉应变出现在车辆轮轴作用时刻，最大压应变出现在相邻两节车辆转向架作用的中间时刻，由最大负弯矩引起；如图12(a)所示，轨道板顶部最大压应变出现在车辆轮轴作用时刻，最大拉应变出现在相邻两节车辆转向架作用的中间时刻，也由最大负弯矩引起，并导致了轨道板顶部受拉损伤的进一步增加，故在进行轨道结构设计时，负弯矩对轨道板顶部受拉损伤的影响不容忽视。

4 结论

本文运用车辆-轨道耦合动力学理论，利用混凝土弹性损伤本构考虑无砟轨道结构的损伤效应，将车辆简化为多刚体系统，建立考虑轨道结构损伤效应的无砟轨道-车辆垂向耦合周期性动力模型，并采用隐式动力预测-校正算法和轨道-车辆系统交叉迭代的求解策略，实现了同时包含材料非线性及接触非线性的含损伤无砟轨道-车辆垂向耦合动力模型的隐式快速求解。所建模型不仅计算结果正确，而且具有非常高的求解效率，方便应用于多参数、多工况的仿真计算，可用于分析和研究无砟轨道结构损伤与车轨动态相互作用的相互影响规律。