基于残余应力模拟的高速铁路路基动力安定分析算法

毕宗琦 叶阳升 蔡德钩 安再展

1.北京铁科特种工程技术有限公司，北京 100081；2.中国铁道科学研究院集团有限公司铁道建筑研究所，北京 100081；3.中国铁道科学研究院集团有限公司高速铁路轨道技术国家重点实验室，北京 100081

建立科学有效的路基结构耐久性评估方法，是实现高速铁路路基高精度变形控制和长久服役须攻克的基础性理论难题。考虑到试验数据量巨大、计算的累积误差以及参数获取困难，传统方法很难直接评估预测路基长时间尺度下的耐久稳定状态。

作为一种有别于常规弹塑性增量分析的临界状态方法，安定理论方法为预测结构的安全性、评估其长期变形趋势和稳定性提供了有效的分析途径。安定理论相关定理和算法最早是由Melan等［1］提出，并由Maier、Ceradini、Polizzotto等［2-4］进一步发展。与传统方法相比，安定性分析不需要精确的荷载历史和时间历程，只需要对外部荷载的边界包络进行确定。因此，可以有效地确定承受大量荷载循环的给定路基结构是否会产生塑性应变的逐渐累积和最终失稳破坏，从而给出高速铁路路基长期服役耐久性的合理判断。近年来，安定理论在岩土工程中的应用越来越广泛，如路面车辙变形［5］、海洋平台地基承载能力［6］、边坡稳定性［7］以及无砟轨道结构［8］的相关评估和分析。目前安定荷载已逐渐成为近代工程强度设计规范中塑性失效准则的重要理论判断依据，为临界动应力的选取、结构安全性评估等提供了有力的分析工具。

在安定分析算法方面，国内外学者分别针对其研究的应用场景和假定条件，给出了相应的分析格式与求解流程。Yu等［9］提出了三维移动交通荷载下黏性摩擦材料安定性的严格下限解，并通过优化程序计算下限安定极限。Qian、Dai等［10-11］通过计算车轮和路面之间滚动和滑动接触引起的弹性应力响应，扩展了柔性路面的安定性分析。Alves等［12］针对板式无砟轨道相关的土体安定分析提供了一种基于2.5D数值模型的方法，分析了列车速度与安定极限载荷的关系。Connolly等［13］通过使用混合模型研究了高速铁路的动力学和安定问题，给出了弹性安定极限用于确定导致轨道破坏的最小荷载。Connolly等［14］将下限安定理论与基于强度的模拟方法相结合，实现了弹性应力生成和安定计算之间的解耦。

在列车反复荷载作用下的路基安定分析中，动力效应的影响对应力应变响应和累积滞回耗能起着重要作用，因此分析中应适当考虑振动频率和阻尼特性的影响。然而，目前大多数安定分析算法都使用了假设准静态载荷条件，动力安定分析算法的相关研究与应用还有待进一步完善。本文基于动力安定分析理论框架，提出一种基于残余应力模拟的动力安定分析算法，并通过典型多层铁路路基数值算例验证该方法的有效性，为路基结构的耐久性评估以及维修养护提供理论支撑。

1 基本格式

根据动力安定定理［15］，结构安定的充要条件是，对于所有初值条件∈ΩI泛函空间和动力作用∈ΩL泛函空间作用下的虚拟弹性动力响应，可以找到一个与时间无关的残余应力分布和时间t*，使结构内部各处均不违反屈服条件，则真实结构将安定于该残余应力状态。具体而言，在获取结构在特定初始条件和外载荷下产生的弹性动力响应后，需在时域内搜索获取任何有可能使结构达到临界安定的应力状态进行安定性检验。据此给出以下动力安定分析格式

式中：λsd为安定荷载乘子；λ为荷载乘子为结构内的应力关于位置坐标求偏导；nj为方向余弦。

该问题最终归结为一个最优化求解的数学规划格式，式中s.t.表示约束条件。其目标函数为荷载乘子λ，对应于结构所受外荷载的大小。求解时从数学规划问题的可行域内部开始搜索最优值。通过构建自平衡残余应力场，对应于载荷域ΩL的相关弹性域ΩE随着载荷的增加而扩展。最终，当λ达到最大值、即安定荷载乘子λsd时，可以确定安定荷载域的极限，并完成动力安定分析。

2 算法构建

2.1 残余应力场的模拟

动力安定分析求解的核心是寻找一个与时间无关的自平衡残余应力场采用应力模拟法，通过构造虚设的自平衡温度参数应力场对高速铁路路基结构的残余应力场进行模拟［16］。以此为基础，构造以节点温度T为变量的残余应力场

设物体的热膨胀系数为α，由温度载荷作用产生的热应变为ε0=αT。对于物体中存在初应变的情况下，引入弹性矩阵D，并利用几何矩阵B及结构位移矩阵U表示应变ε=BU，应力应变关系表示为

结构应变能表示为

式中：K为刚度矩阵，K=BDB；R为常数项，R=为温度应变引起的载荷项，Q=

根据最小势能原理，对式（3）进行变分，并引入应力矩阵S、关系矩阵G和H，可得以节点温度T为自变量表示的温度参数残余应力场为

将式（4）代入式（1）可得

2.2 屈服条件的线性化处理

在动力安定分析格式中，将非线性屈服函数f进行线性化处理，从而将大规模非线性规划问题转化为线性规划问题，降低问题求解难度。采用包含土体摩擦强度参数的Mohr-Coulomb屈服准则进行模型的安定性判别。在2τxy-（σx-σy）应力空间坐标系中将一点的屈服条件转化为l个线性方程，屈服条件不等式表示为

式中：σx、σy、τxy分别为x、y方向的正应力以及剪应力；Ak、Bk、Ck、Dk为系数，Ak=cos（2πk/l）+sinφcos（π/l），Bk=sinφcos（π/l）-cos（2πk/l），Ck=2sin（2πk/l），Dk=2ccosφcos（π/l）。

令

式中：N为线性化矩阵，由各线性屈服面的外法线向量组成；r为线性化后各屈服面至坐标原点的距离。

屈服条件转化为

代入（5），并考虑与变载对应的弹性应力σde（t）、与恒载对应的弹性应力σse，得

2.3 变量数目的控制

变量数目对求解的效率有关键性影响。对于利用大量节点参数描述残余应力的场变量，适当减少优化变量的个数、缩减规划问题的规模将有效提升计算速度。本算法考虑在模型中利用非均匀有理B样条插值手段，合理选取若干控制点赋予其新的变量，表征整个模型的残余应力场，从而减少变量数目，提高算法的运行效率。具体实施步骤如下：

1）给定n×r维控制点Pi，j，Pi，j定义在笛卡尔坐标系中，并给出对应于Pi，j的二维权系数wi，j。

2）给定在两个维度各自的插值次数p和q，依据定义，控制点Pi，j在两个维度的向量长度分别为n+1和r+1，按归一化B样条插值基函数的性质可知两个维度上节点向量最大编号分别为m1=n+p+1和m2=r+q+1，节点向量长度分别为m1+1和m2+1。

3）根据两个维度给出节点向量U=｛α，…，α，up+1，…，um1-p-1，β，…，β｝、V=｛α，…，α，vq+1，…，vm2-q-1，β，…，β｝，其中α和β的重复度分别为p+1和q+1。

4）给出NURBS有理多分式矢函数S(u，v)=其中u、v为坐标系中两个方向的参变量，Ni，p（u）、Nj，q（v）分别为u、v方向上的归一化B样条插值基函数。插值基函数通常可利用Cox-de Boor递归公式获取。

5）在矢函数S（u，v）的基础上，构建坐标参数与拟合目标参数的关系。每个控制点Pi，j对应一个变量Ti，j，分量形式表示为P（x，y，T）=（Px，Py，PT），相应的矢函数S（u，v）=［Sx（u，v），Sy（u，v），ST（u，v）］。生成拟合场后，利用反距离加权IDW法获取模型网格节点的具体变量值。

3 求解流程

求解时，首先利用动力响应分析模块，分别求出外荷载各基准荷载分量P（x，t）作用下的弹性动应力场σ*（x，t）；随后在模型空间内部设置控制点，其节点温度参数变化Tc作为自变量，生成整个模型区域的自平衡残余应力场σt（x）；将路基结构各点位动力、残余应力场叠加后，根据屈服函数设置约束条件，构建最优化问题分析格式并迭代求解，最终获取安定荷载乘子的最大值λsd，完成高速铁路路基结构安定分析。

图1 动力安定分析算法求解流程

4 算例应用

铁路路基动力安定分析的数值模型是根据文献［17］的试验断面建立，如图2所示。路基总高度为6.86 m。模型中的构件主要采用4节点平面应变单元离散化，如图3所示，同时给出了基于NURBS的残余应力场的控制点的分布。图4给出了在迭代中不断更新对的残余应力场插值曲面，其形态由各节点的参数变量和材料性质决定，反映了自重和外部列车荷载加卸载后土体中残余应力的可能分布形式。

图2 一高速铁路路基结构断面（单位：m）

图3 高铁路基结构模型与控制点设置

图4 残余应力场模拟插值曲面

安定分析中模型材料参数见表1，其中Ed为动弹性模量；ρ为密度；μ为泊松比；ξ为阻尼比；c为材料的黏聚力；φ为材料内摩擦角。假设路基本体为各向同性、均质材料。混凝土板视为线弹性材料。

表1 高速铁路路基安定分析模型材料参数

动力分析得到的路基动应力衰减变化曲线以及现场试验得到的实测数据见图5。可知，模型中的动应力衰减与试验结果吻合较好，验证了模型得到的动态响应基本规律的可靠性。

图5 动力响应分析结果

通过本文给出的算法流程，路基结构的安定极限可通过计算无量纲参数λsd获得。动力安定分析迭代收敛曲线见图6。可知，经过迭代分析，获取了最优值λsd=1.47。根据现场试验的测量值，结构表面所受列车荷载实测值P0约为92 kPa，将其视为基准荷载参考值，其与λsd的乘积表示外部荷载在参考值的不同倍数下的振幅。当λ=λsd时，外载荷幅值Psd=λsdP0=135.24 kPa，即为算法获取的路基的安定阈值。

图6 动力安定分析迭代收敛曲线

安定阈值与路基结构的变形累积发展趋势有密切的关联。根据λsd的预测值选择一系列荷载乘数λ1，施加振幅P1=λP0执行荷载循环，提取等效塑性应变的发展，对不同外荷载循环下的高速铁路路基塑性应变发展情况进行研究，见图7。可知：根据等效塑性应变随时间的不同发展趋势，主要分安定和连续累积两组。载荷振幅低于安定阈值时（P1=92、120、130 kPa），等效塑性应变的值在初始阶段后不再增加，材料的应变-应力行为最终趋于线性响应，这意味着模型达到安定状态。相反，当荷载超过安定阈值时（P1=133、135、140、150 kPa），等效塑性应变在每个循环中继续累积。因此，将外荷载作用下的路基动力响应控制在安定阈值内，是控制路基结构不发生无限制变形累积，维持耐久稳定状态的重要前提。

图7 路基结构塑性区等效塑性应变发展趋势及安定状态判别

5 结语

本文提出了一种基于残余应力模拟的高速铁路路基动力安定分析算法和数值求解程序，用于研究铁路路基的安定极限。通过在结构中使用控制点并结合NURBS插值，构建了与时间无关的自平衡残余应力场和相应的算法格式。通过应力模拟、屈服面线性化处理以及变量数目优化，实现了算法求解效率的提升。

高速铁路路基的动力安定极限无量纲参数λsd=1.47，与现场试验荷载基准值92 kPa对应的安定阈值为135.24 kPa。等效塑性应变的发展曲线给出了路基是否能够达到稳定状态或变形持续累积状态的分类标准，给出了安定和连续累积两大类累积趋势。

本文方法可以进一步扩展到考虑材料硬化的影响、使用不同的应力-应变屈服函数，或更复杂精细的路基结构形式。