前几周帮某个小课题搭程序要分析一座简支梁桥在列车通过时的竖向动力响应。第一版程序图省事用的是移动荷载模型把列车简化成一组按速度移动的集中力。结果一出来就傻眼跨中位移和实测差了差不多两成墩顶位置的时程曲线对不上频谱里还多出一段根本不该出现的振动分量。后来才反应过来问题出在把列车当成“扰动源”而不是“结构”。车辆本身是有质量、有悬挂系统的多自由度振动体系轮轨接触面之间还有轨道不平顺在不断地产生强迫激励。这些效应只要速度一上来、桥跨一敏感就会被移动荷载模型完完全全漏掉。于是换成了完整的三维车桥耦合分析车辆子系统、无砟轨道与桥梁子系统共同组成一个大耦合系统用Newmark法做时域数值积分轨道不平顺通过轮轨接触点的几何位移导入。这条路线很成熟但真要把MATLAB程序跑通中间有不少细节不是书上写的那样明明白白。这篇文章就把完整思路、方程、代码结构和调参过程梳理一遍给正在做同类程序的兄弟们一个能直接参考的底稿。1. 为什么非要用车桥耦合模型移动荷载方案在什么条件下根本不可用1.1 桥上跑车时真正发生的是什么先别急着写矩阵。得把物理过程搞清楚。列车过桥时桥上承受的并不只是轮子的重力。车辆通过弹簧悬挂压在轮对上轮对又通过轮轨接触压在轨道上轨道通过扣件、道床传给桥梁结构。这个传力路径里最要命的是车辆本身在振动。车体上下颠、轮对上下跳簧上质量和簧下质量之间是弹簧和阻尼器连接。车辆振动频率一旦和桥跨固有频率接近就会发生强烈的动力放大。移动荷载模型把车辆质量全忽略了也就是说不存在“车桥质量耦合”。对于跨度小、质量轻的桥车辆质量可能占到桥上总质量的相当比例忽略质量就等于把一个关键的动力参数砍掉了。再考虑一个因素车轮驶过不平顺轨道时轮轨之间会产生高频冲击力这个力通过轨道传给桥梁。这属于强迫振动的输入移动荷载模型压根没法统一表述。1.2 无砟轨道把问题变成了“三个子系统”题目里写的是“无砟轨道-桥梁”耦合。和有砟轨道不同无砟轨道的钢轨是直接固定在轨道板上的轨道板下面用砂浆层和底座板联结底座板再和桥面连接。它的整体刚度比有砟轨道高一个量级几何位置稳定所以在动力计算里有一个好处轨道板可以当作连续梁元或板元参与建模而不必像有砟道砟那样处理成大量离散颗粒。建模时可以有两种处理方式。第一种是精细的三层建模钢轨、轨道板、底座板各用一层梁单元扣件、砂浆层用弹簧阻尼单元。这种模型自由度多计算量大但能细致反映轮轨力向桥梁传递的过程。第二种是简化处理把无砟轨道看成桥面的一部分直接对桥梁截面做截面刚度叠加轨道板和梁体共同工作。这种方式对于大跨度桥的总体响应计算足够用但对局部轮轨力分布的分辨率不足。我这次程序采用的是折中方案钢轨单独作为一层梁它的节点与轨道板节点之间用扣件弹簧连接轨道板与桥梁主梁节点之间用等效刚度层连接。这样既保留了钢轨的局部振动特征又不至于让自由度爆炸。三层无砟轨道-桥梁子系统的自由度组成大致是钢轨梁节点竖向位移、转角轨道板节点竖向位移、转角桥梁主梁节点竖向位移、转角轮轨接触力作用在钢轨节点上通过弹簧层向轨道板、桥梁传递。这个弹簧层参数的取值很关键后面专门讲。2. 三个子系统的方程怎么写进同一个程序2.1 车辆子系统四分之一车模型的自由度与参数车辆模型最常用的入门级选择是四分之一车模型也就是把单个轮对位置上的车辆动力学简化成两个自由度簧上质量车体那部分质量和簧下质量轮对质量中间用悬挂弹簧和阻尼连接。方程是标准的二自由度振动方程[ m_s \ddot{y}_s c_s (\dot{y}_s - \dot{y}_w) k_s (y_s - y_w) 0 ][ m_w \ddot{y}_w - c_s (\dot{y}_s - \dot{y}_w) - k_s (y_s - y_w) k_h (y_w - y_r - r) 0 ]其中 ( y_s ) 是簧上质量位移( y_w ) 是轮对位移( y_r ) 是钢轨位移( r ) 是轨道不平顺幅值( k_h ) 是轮轨接触刚度( c_s ) 和 ( k_s ) 是车辆悬挂参数。[ M_v \begin{bmatrix} m_s 0 \ 0 m_w \end{bmatrix},\qquad C_v \begin{bmatrix} c_s -c_s \ -c_s c_s \end{bmatrix},\qquad K_v \begin{bmatrix} k_s -k_s \ -k_s k_s k_h \end{bmatrix} ]车辆对桥梁产生的力本质上是轮轨接触点上的 ( k_h (y_w - y_r - r) )这个力同时反向作用于桥梁子系统。这个项就是耦合的核心。如果有整车模型的需求可以在这个基础上扩展到前后两个轮对或者左右两侧轮对也就是半车模型或整车模型自由度会增加到 4个、6个甚至更多。四分之一车模型对单轮激励下的竖向响应分析足够我在程序里先用它验证整体逻辑。2.2 桥梁子系统有限元离散的思路与矩阵形式桥梁采用传统有限元梁离散。每个节点两个自由度竖向位移和转角位移。单元用欧拉-伯努利梁单元质量矩阵用一致质量矩阵刚度矩阵是标准单元刚阵。系统中还包括钢轨和轨道板它们同样用梁单元描述。整个桥梁子系统的矩阵组装过程就是把钢轨单元、轨道板单元、桥梁单元分别形成单元刚阵按节点编号次序叠加到总体坐标里。扣件刚度和砂浆刚度也加进去。特别注意钢轨和轨道板之间的扣件刚度如果直接以集中弹簧单元形式加进去应该在其两个端点节点上分别添加对应刚度项。例如在钢轨节点 ( i ) 和轨道板节点 ( j ) 之间刚度值 ( k_f ) 加入总体刚度矩阵的位置是[ K(i,i) k_f,\quad K(j,j) k_f,\quad K(i,j) - k_f,\quad K(j,i) - k_f ]阻尼矩阵在这个问题里比较麻烦。桥梁结构的阻尼通常采用瑞利阻尼假设( C_b \alpha M_b \beta K_b )系数由结构前两阶自振频率与阻尼比反推。钢轨和轨道板的阻尼可以考虑得简单一些给一个较低的模态阻尼比。2.3 轮轨接触与不平顺位移怎么作为边界条件进方程轮轨接触处理是车桥耦合里最关键的一步。在考虑不平顺的情况下轮轨间的弹性压缩量不只是轮对位移和钢轨位移的差还要减去不平顺幅值 ( r )。当车轮滚过钢轨上某一个点如果钢轨该位置本身因为不平顺高了 0.5mm轮对即使不动接触力也会增加。接触力的表达式[ F_c k_h (y_w - y_r - r) ]其中 ( r ) 是随时间和轮对位置变化的不平顺值。这个力同时出现在车辆方程和桥梁方程中方向相反。于是可以把它直接整合进系统矩阵[ M \ddot{X} C \dot{X} K X F(t) ]右端项 ( F(t) ) 中包含了由于不平顺造成的附加激励 ( k_h r(t) )而系统矩阵中已经包含了接触刚度 ( k_h ) 的贡献。这种做法的好处是无需在时间循环里反复做接触力迭代。传统“分离迭代法”在每个时间步需要在车辆和桥梁之间来回传力直到收敛。直接耦合法一次组装、一次求解收敛性天然有保障程序结构也清晰得多。3. 不平顺时域样本的生成功率谱密度转位移的完整做法3.1 不平顺功率谱密度函数怎么取轨道不平顺属于随机激励工程上通常用功率谱密度来描述它的频域特征。国内外都发展过对应的谱函数形式上多为分段幂函数或多项式的组合。这里用一个常见形式的竖向轨面不平顺谱[ S(\Omega) \frac{A}{\Omega^k B} ]其中 ( \Omega ) 为空间角频率单位是 rad/m( A )、( B )、( k ) 由轨道等级决定。等级越高不平顺越小。MATLAB里生成不平顺样本时需要先把谱函数从频域转换为空间域波形。最常用的方法是三角级数叠加法[ r(x) \sum_{n1}^{N} \sqrt{2 S(\Omega_n) \Delta \Omega} \cos(\Omega_n x \varphi_n) ]其中 ( \varphi_n ) 是 ([0, 2\pi)) 范围内的随机相位( \Delta \Omega ) 是空间频率离散间隔( \Omega_n ) 是离散频率点。3.2 MATLAB里用三角级数叠加法生成空间不平顺直接给一个能跑的MATLAB函数function [r] genTrackIrr(L, dx, A, k, B, seed) % 生成轨道竖向不平顺序列 % L: 轨道总长度, m % dx: 空间采样间隔, m % A, B, k: 不平顺谱参数 % seed: 随机种子 x (0:dx:L); nPoints length(x); Omega_min 2*pi/L; % 最低空间频率 Omega_max pi/dx; % 最高空间频率采样定理 nFreq 2000; Omega linspace(Omega_min, Omega_max, nFreq); dOmega Omega(2) - Omega(1); S A ./ (Omega.^k B); % 功率谱密度 rng(seed); phase 2*pi*rand(nFreq, 1); Amp sqrt(2 * S * dOmega); r zeros(nPoints, 1); for i 1:nPoints r(i) sum( Amp .* cos(Omega * x(i) phase) ); end % 去均值避免初始位移偏移 r r - mean(r); end注意几个关键点最低空间频率由样本总长决定。如果轨道长度只有几十米长波成分根本表现不出来样本的整体形状会偏“干瘪”。最高频率由空间采样间隔决定也就是 ( \Omega_{max} \pi / dx )。如果程序里忽略了这一点直接拿一个随意设置的高频上限生成出来的样本会在数值上失真。相位必须随机。同一组谱参数下换一组随机相位就等于换一条轨道。要做参数对比时应该固定同一个种子否则响应差异会把轨道随机性的影响和目标参数的影响混在一起。随机样本的统计特性需要多次试验平均才能和理论谱吻合。单次样本的功率谱会有较大波动这属于正常现象。3.3 你会踩的坑频率上限、两端处理和速度换算踩坑经验说几个第一车轮移动过程中轮对所在位置的 ( x ) 坐标随时间变化( x v \cdot t )不平顺值也要按这个位置实时插值。程序中不是把整个不平顺数组一次性赋给系统而是在每个时间步根据轮对坐标从生成的样本里插值取值。因此生成样本时轨道长度必须足够覆盖列车过桥的整个路程通常取桥梁长度加上一个过渡段。第二桥梁两端的不平顺样本最好做衰减处理。实际轮对是从引桥上驶入主桥的如果样本在桥端突然截断或者存在过大的初相位程序启动时会给系统一个虚假的初始冲击。更稳妥的做法是在计算区域内正常施加不平顺在桥端范围各留出一个过渡区用窗函数对不平顺做幅值渐变。第三时间采样频率与空间采样频率要匹配。车速为 ( v )轮对经过空间采样间隔 ( dx ) 所需时间为 ( dx/v )Newmark积分的时间步长应该取这个间隔的若干分之一。比如 ( dx 0.05m )车速 ( v 80 m/s )那么轮对每 0.000625s 经过一个采样点时间步长至少小于这个值的 1/5程序才有足够分辨率。4. Newmark法不是“黑盒”参数、迭代和收敛逐个说透4.1 Newmark公式与无条件稳定参数的选择逻辑Newmark法是一种隐式逐步积分法基本公式是[ u_{n1} u_n \Delta t \dot{u}_n \frac{\Delta t^2}{2} \left[ (1-2\beta)\ddot{u}n 2\beta \ddot{u}{n1} \right] ][ \dot{u}_{n1} \dot{u}_n \Delta t \left[ (1-\gamma)\ddot{u}n \gamma \ddot{u}{n1} \right] ]参数 ( \beta ) 和 ( \gamma ) 决定了算法的稳定性和精度。最常用的取值是 ( \beta 0.25 )、( \gamma 0.5 )对应常数平均加速度法其特点是无条件稳定且能量不会因时间积分放大。无条件稳定意味着时间步长可以取到比结构最高频率对应周期大一些的数值程序也能算下去不会发散。但注意数值不发散不等于结果准确。对于高频成分过大的时间步长相当于给系统加了一个低通滤波真实的高频轮轨力会被严重抹平。如果只关心桥梁低频响应( \Delta t 1/(100 f_1) ) 可能就够用其中 ( f_1 ) 是桥梁基频。但如果要分析轮轨力的高频冲击最高关心频率假设是 200Hz时间步长至少要小于 0.0005s。4.2 在MATLAB里实现多自由度Newmark的循环架构在多自由度系统中Newmark法通常需要把动力学方程改写成等效静力形式业内一般称为“等效刚度法”。先定义常数a0 1/(beta*dT^2); a1 gamma/(beta*dT); a2 1/(beta*dT); a3 1/(2*beta) - 1; a4 gamma/beta - 1; a5 dT/2*(gamma/beta - 2); a6 dT*(1 - gamma); a7 gamma*dT;在每个时间步需要形成等效刚度矩阵并解方程Keff K a0*M a1*C; for n 1:nSteps F_eff F(n,:) M*(a0*u(:,n) a2*v(:,n) a3*a(:,n)) ... C*(a1*u(:,n) a4*v(:,n) a5*a(:,n)); u_next Keff \ F_eff; a_next a0*(u_next - u(:,n)) - a2*v(:,n) - a3*a(:,n); v_next v(:,n) a6*a(:,n) a7*a_next; u(:,n1) u_next; v(:,n1) v_next; a(:,n1) a_next; end这里有几个效率优化点( K_{eff} ) 组装一次就保持不变可以在循环外预先分解Keff_factored decomposition(Keff)或[L,U,p] lu(Keff)循环里直接用换元法回代避免每步重复做一次矩阵分解。对于上千自由度的系统这一项优化能把总耗时减少一半以上。右端力向量 ( F(n,:) ) 每步需要更新一次因为轮对位置在变不平顺激励也在变。矩阵规模如果超过两三千自由度且时间步数很多MATLAB 的循环效率可能撑不住。建议把每若干个时间步的结果临时存下来而不是每一帧都写入磁盘。4.3 耦合迭代的收敛判据与接触力的平滑处理如果采用分离迭代的方式来处理车辆和桥梁的耦合那么每一次时间步内需要进行迭代假设当前轮轨力求解车辆子系统得到轮对位移。把轮对位移代入桥梁子系统求解钢轨接触点位移。用新的位移差值计算轮轨力和上一轮迭代的结果比较。收敛条件设为接触力相对变化小于某一阈值例如 0.5%。直接耦合法省去了这个迭代过程但也有一个代价系统的总体矩阵维度扩大了车辆自由度与桥梁自由度合在一起。如果车辆模型复杂整车多轮对维度扩大的比例会相当可观。权衡之下我给这个程序选择了直接耦合法换来的是稳定性和程序结构清晰。接触力的平滑处理是另一个容易被忽略的点。轮轨力本身随不平顺剧烈波动如果提取出来直接画图波形上会有大量锯齿。这不一定是数值错误可能是真实的随机激励响应。但如果锯齿非常密集呈明显的“麻点状”就要检查时间步长是否过大、轮轨赫兹刚度取值是否过于刚硬以及不平顺高频分量是否需要做低通滤波。5. 程序跑起来的全流程、调试技巧和后期处理5.1 从零搭一个可运行脚本的步骤顺序整个程序建议按六个模块组织每个模块单独写成函数或脚本段落参数定义区桥梁参数、轨道参数、车辆参数、车速、时间步长、总时长。子系统矩阵生成区车辆M/C/K、桥梁M/C/K、耦合矩阵组装。不平顺生成区按照第3节的函数生成轨道样本。初始条件设置区位移、速度、加速度初值。时间循环求解区Newmark主循环。结果后处理区提取跨中位移、轮对位移、轮轨力、加速度绘图分析。模块化最大的好处是便于定位错误。程序一旦算错如果你是几万行挤在一个脚本里光找错就能耗掉半天分模块的话可以单独验证桥梁自振频率单独检查车辆自由衰减确保子系统没毛病再联调。5.2 慢速起步几个特别容易出错的地方根据实测调试经验最常见的错误集中在下面这几类问题现象可能原因处理办法系统发散位移随时间指数增大时间步长过大矩阵组装错误阻尼为负缩小步长打印矩阵对角线检查正定性响应波形出现规则尖峰不平顺插值采样导致激励突变检查轮对位置与时间步的匹配关系桥梁自振频率和理论解对不上边界条件加载错误刚阵/质量阵组装遗漏先跑模态分析验证再跑动力响应车辆弹簧刚度参与耦合后矩阵奇异接触刚度与车辆刚度叠加位置冲突逐项检查刚度矩阵对角元素其中矩阵组装错误是最隐蔽的。在MATLAB里拼装大矩阵时节点编号排序错了会导致矩阵本身不错但物理意义错乱。建议组装完成后做一个自检运行一次模态分析比对桥梁前三阶频率是否与理论值一致。这个自检必须在施加耦合外力之前做不然分不清误差来源。5.3 计算结果怎么取位移、加速度与轮轨力的输出和复核算完之后先别急着画图。有几个步骤值得做。输出跨中竖向位移时程用快速傅里叶变换看频谱桥梁基频上应该有一条明显的峰。如果谱线最高峰出现在零频附近说明结果被某种趋势项污染通常是因为初始条件没处理好或者某阶不参与振动的高频成分没有滤掉。轮轨力是判断耦合模型是否合理的重要依据。静态情况下轮轨力应当等于车辆重力分配到该轮对上的荷载。加上动力效应后轮轨力围绕静轮重上下波动。如果程序跑出来的轮轨力长期偏离静轮重且不回落多半是接触刚度或不平顺均值处理出了问题。加速度信号的输出需要特别谨慎。Newmark法对加速度的求解是后验的即先求位移再通过位移差除以 ( \Delta t^2 ) 推回来因此对时间步长更敏感。如果时间步长不够小加速度时程上会出现人工振荡而位移时程看起来还正常。经验法则是加速度结果的可信程度至少要求时间步长小于关心频率周期的 1/30。6. 个人经验车速、不平顺和时间步长这三者怎么联动6.1 三者之间的物理关系车桥耦合程序跑出来的结果好不好很大程度上在调参阶段就决定了。车速影响的是激励频率。不平顺在空间上是固定的波形车辆以速度 ( v ) 通过时就把它转化成了时间激励激励频率 ( f v \Omega / 2\pi )。同样的不平顺波长车速越高激励频率越高。如果某一段不平顺的波长是 2m车速是 60m/s对应的激励频率就是 30Hz。这一频率如果和桥梁某一阶模态频率接近就会产生共振响应会异常放大。我在程序调试中见过这种情形把车速从 50m/s 调到 55m/s跨中位移峰值突然上涨了 30%。一开始以为是程序出错了后来确认就是车速接近了某阶模态的共振区。时间步长的选择必须同时满足驱动频率分辨率和数值稳定性的要求。我的经验是分两步先按车辆最高关心频率设定步长上限再按不平顺空间采样点的时间间隔进行校验取两者中的较小值作为最终步长。6.2 给新手的建议先跑通小模型再上大模型如果你刚接触这类程序我强烈建议先做一个 10 个桥梁单元加四分之一车模型的小算例。节点少矩阵维度小单步循环时间短你可以清楚地观察到每一步更新过程。等这个小模型的结果跟文献算例对上了再扩大自由度、增加轮对、改成整车模型。我用小算例做验证时发现过一个有意思的事车辆悬挂刚度对桥梁动力响应的影响并没有想象中那么大反而是阻尼比的影响更明显。这是因为桥梁的响应主要由低频成分主导车辆悬挂的频率相对桥梁更高两者的耦合主要存在于中高频段。如果你做的项目只关心跨中位移那四分之一车模型完全够用但如果要关心轮轨力和钢轨高频振动就必须考虑更精细的车辆模型和轨道模型。最后的实际体会是车桥耦合程序调试过程中最难的不是算法或公式而是数据的一致性和边界条件的合理性。一个符号错了可能跑半天才暴露出来一个单位错了可能整批结果报废。建议在程序里统一采用国际单位制并设一个参数自检段把所有关键参数打印到命令行再进入主循环。程序跑通之后再把轨道不平顺的种子换几组做统计平均你会发现单次响应和多次平均响应之间差距非常大——不定量的随机分析永远只能看到现象的一面。
