随着国民经济的发展，高速铁路轻量化，大跨度铁路桥梁及桥上列车，易受湍流的脉动作用而发生影响行车安全的抖振响应[1-3]。复气动导纳作为阵风场和抖振力的幅值和相位角在频域内的传递函数，描述了抖振力的非定常特性，对于列车抖振力的精细化模拟至关重要[4-5]。因此，为准确评估列车行车安全和乘客舒适度，有必要研究阵风场中桥上移动列车的复气动导纳。从流体力学的基本控制方程推导列车复气动导纳的理论解非常困难，风洞模型试验仍是桥上移动列车复气动导纳研究的重要方法[6-9]。虽然桥上移动列车的气动参数与静止列车存在差异[4]，但由于试验设备、测试装置和研发成本的限制，常常使用静止列车模型来代替移动列车模型开展风洞试验[10-12]。马存明等[6]基于静止列车模型风洞试验，研究了横风作用下列车位置、风攻角对列车气动传递函数的影响；何旭辉等[7]基于节段模型风洞试验对双层桁架主梁断面上列车进行测力和测压，分析了静止列车在不同车道和桥层的气动特性。风洞试验的湍流场由格栅或尖劈被动生成的风场所模拟，并假定顺风向和竖向湍流分量对复气动导纳的贡献相同[13-14]，在该湍流场下测取的列车复气动导纳仅含幅值的传递特性，而湍流风场和抖振力之间相位角的频域特性则由于测量技术的限制无法获得。可见，现行的列车复气动导纳识别方法无法考虑列车抖振力和湍流场间的相位角和湍流分量的影响，基于复气动导纳的幅值计算的列车抖振阻力和升力有可能同时达到极大值，导致列车抖振响应的计算结果偏于保守[3-4]。复气动导纳能同时描述阵风和抖振力的幅值与相位角的传递特性，主动阵风生成装置的实现使得复气动导纳的试验识别成为现实。基于该装置，DIANA等[15-16]识别了流线型箱梁和分离式三箱梁的复气动导纳；韩艳[17]识别了薄平板断面的复气动导纳，识别结果与薄机翼理论解Sears函数趋近；LI等[18]识别了分离式双箱梁断面的复气动导纳并发现，其结果与薄机翼理论存在差异。可见，目前复气动导纳研究主要集中于桥梁断面。由于桥上移动列车试验的复杂性和昂贵成本，限制了主动阵风装置在列车复气动导纳的试验识别方面的应用。计算机技术的发展导致计算流体力学技术(CFD)在桥梁风工程领域焕发活力。利用CFD求解器Fluent中用户自定义函数和修改参数灵活的特点，可以模拟满足入流自保持能力的二维阵风场。基于CFD二维阵风模拟技术，研究人员识别了矩形断面[19]和典型桥梁断面[20]的复气动导纳，且识别结果与试验值吻合良好。需要指出，复气动导纳属于二维气动导纳，“片条假设”的成立是复气动导纳应用的前提条件，湍流场的三维特性和列车有限展长的特点使“片条假设”的适用性有待商榷。然而，KAVRAKOV等[20]指出，当钝体结构的展长与特征宽度的比值达到一定程度时，“片条假设”是成立的。上述钝体结构复气动导纳的数值识别和理论研究成果，为列车复气动导纳的数值识别建立了信心。本文针对阵风场中桥上移动列车，基于列车气动力理论建立列车复气动导纳的识别理论。利用CFD数值模拟技术，以列车通过箱桁组合梁断面为研究对象，识别桥上列车的复气动导纳，研究脉动风分量和列车位置对列车复气动导纳的影响。基于抖振力的三维谱张量理论，论证了“片条假设”的适用性，为复气动导纳用于桥上列车的抖振力模拟提供了理论支撑。研究结果可为桥上移动列车抖振力的精细化模拟和相关研究提供参考。
1 数值识别方法
1.1 流体控制方程
对于二维不可压缩流体，流场的控制方程为纳维-斯托克斯方程：
式中：U为速度矢量；v为空气分子黏度；vt为湍流黏度；
对于二维阵风场，其必须满足连续性方程和边界条件。根据LI等[19]的研究，可采用如下表达形式的二维阵风场：
式中：Aw和Au分别为竖向和顺风向阵风场的振幅，m/s；f为阵风的频率，Hz；H为流场的高度，m；n为涡旋数。
由方程(2)可知，当涡旋数n=1和n=2时，计算域的中部将分别产生竖向和顺风向谐波阵风场。根据RIBNER[22]的三维谱张量理论，多频谐波阵风场可视作多个单频谐波阵风场线性叠加的结果。因此，为减少计算耗时，将采用多频叠加的方法生成二维阵风场，该方法的有效性已在王鸿轩等[21]的研究中得到试验验证。此外，LI等[23]的风洞试验研究表明，湍流度对钝体结构的二维气动导纳几乎没有影响。因此，参考已有文献[21, 25]，数值模拟中，阵风场的幅值均设置为0.336 m/s，来流平均风速U=10 m/s，相应的湍流度为2.38%，这也符合SEARS[26]所要求的脉动风幅值要远小于来流风速的假设。此外，实际湍流风场的积分尺度通常等于甚至数倍于列车的宽度，对于本文数值模拟的阵风场而言，最大脉动频率为f=8 Hz，对应的二维阵风场的最小积分尺度为U/f=1.25 m，接近2.2节缩尺后列车宽度的30倍，这符合Sears函数推导的前提条件[26]。
1.2 计算域和边界条件
如图1所示，数值模拟的计算域为长18B，高6B的矩形，其中，B为主梁宽度，经试算，该计算域尺寸能够确保主梁的漩涡脱落不至于触及上下边界，避免了回流对气动力的影响。计算域左侧为速度入口，可由Fluent的自定义函数将速度函数加载至入口边界，右侧为压力出口。为了克服传统的周期边界造成的边界波动特征不协调，上下侧均为速度入口，Y方向速度随时间和空间位置的不同而呈简谐脉动变化，以确保波动特征与左侧边界相同[25]。数值模型所在位置为坐标系原点，模型表面采用无滑移壁面条件。
1.3 列车复气动导纳识别理论
假设桥上移动列车处于平均风速为U的湍流场中，此时，体轴坐标系下列车中车断面受到的气动力如图2所示。
根据钝体空气动力学理论，LI等[23]建立了考虑顺风向和竖向脉动风的移动列车抖振阻力和抖振升力的表达式，在此基础上，严乃杰[3]基于小风攻角的假设，推导得到体轴系下移动列车抖振阻力和抖振升力的表达式，分别为
式中：FS(t)和FL(t)分别为列车受到的抖振阻力和抖振升力时程；VR=(U2+v2)1/2为作用于移动列车的合成风速，U为来流平均风速，v为列车的移动速度；CS(α,φ)和CL(α,φ)分别为与风攻角α和风向角φ相关的阻力系数和升力系数；
假设主梁所在位置顺风向和竖向阵风场的表达形式分别为u(t)=Aueiωt和w(t)=Aweiωt，将其代入方程(3)并对方程的两边进行傅里叶变换，整理后得到列车复气动导纳的表达式为
其中，
现有的关于列车复气动导纳的研究，多是针对静止列车而言，且认为风向角φ=90°时为最不利情况，此时列车受到的合成风速VR=U[3]。实际列车在运行过程中，列车受到的合成风速和周围的绕流特性会发生改变，与静止列车存在差异[4]。然而，列车作为线状的结构，风沿结构的展向流动可以忽略[8, 11]，多数学者[7, 11-12]仍基于桥上静止的列车开展气动参数识别研究。因此，本文以桥上静止的列车中车断面为研究对象，开展桥上列车的复气动导纳数值识别。
2 数值模型及验证
2.1 工程背景
如图3所示，以某在建的主跨320 m的斜拉桥为研究对象，其主梁断面采用箱桁组合梁双层布置形式，上层为公路8车道布置、桁式结构，下层为双线高速铁路时速350 km+4车道公路、箱形断面。上层和下层桥面总宽度分别为33.5 m和39 m，桁高为13.2 m。
2.2 数值模型及网格划分
数值模拟中，列车模型采用CRH3，模型的缩尺比为1∶80，以2车交会的主梁断面为例给出了网格划分示意图，如图4所示。主梁和列车表面设置15层边界层，首层高度为0.000 01 m，网格增长率为1.2，以确保首层边界层网格节点位于黏性底层。参考唐煜等[25]的研究，为了使得阵风场的简谐速度脉动幅值的衰减率足够小，主梁周围采用非结构化三角形网格加密，网格增长率为1.15，最大网格尺寸为0.005 m，此时无量纲的时间步长等于1，速度脉动幅值的衰减率满足要求。计算域的其余部分均采用结构化的四边形网格。网格总数为583 126个。经试算，该网格划分策略下，主梁壁面的最大y+值小于1.5。
数值模拟时，均匀流场和二维阵风场下求解器的湍流模型均选择SST k-ω模型，这与已有研究[19, 21, 25]一致。时间推进格式为二阶隐式算法，压力和速度耦合方程的求解采用SIMPLEC算法，其余变量的空间离散均为二阶迎风格式，残差检测界限设为1×10-6，数值求解在求解器软件Fluent 17.0中完成。
2.3 工况设置
数值模拟的工况包括2个部分，一是均匀流场下列车的静力三分力系数，二是顺风向和竖向二维阵风场下列车的复气动导纳识别。工况设置如表1所示。需要指出，均匀流场下，图1所示计算域的上下边界为对称边界，左侧速度入口为仅含x方向恒定风速U的速度入口，其余设置保持不变。
| 工况 | 流场类型 | 列车位置 | 风攻角/(°) |
|---|---|---|---|
静力三分 力系数 | 均匀流 (U=10 m/s) | 迎风侧、背风侧和两车交会 | -5～+5 |
| 复气动导纳 | 阵风场(顺风向和竖向) | 迎风侧、背风侧和两车交会 | 0 |
2.4 网格和时间步长独立性及数值模型精度验证
静力三分力系数的定义为
(5(a),(b),(c))
式中：CD、CL和CM分别为阻力系数、升力系数和力矩系数；Ds、Ds和Ds分别为单位展长主梁上的阻力、升力和力矩的平均值。
为验证数值模型的网格独立性，如表2所示，以0°风攻角下主梁断面的静力三分力系数为例，给出了3种网格划分策略下数值模拟的结果，模拟的时间步长暂取为Δt=0.000 5 s。从表2可以看出，随着网格数量的增加，静力三分力系数逐渐收敛。综合考虑网格数量和计算耗时，选择58万网格作为计算方案。
| 网格数目/万 | 阻力系数 | 升力系数 | 力矩系数 |
|---|---|---|---|
| 22 | 0.784 | -0.186 | 0.054 |
| 58 | 0.829 | -0.219 | 0.059 |
| 86 | 0.831 | -0.220 | 0.061 |
以58万网格的数值模型为基准，基于升力系数进行时间步独立性验证，时间步长设置为0.000 1，0.000 2，0.000 5和0.001 s。分析发现，时间步长越小，精度越高，但计算耗时越长。相对于0.000 1 s的模拟值，当Δt=0.000 5 s时，相对误差仅-0.9%，综合考虑精度与耗时，取步长0.000 5 s。
为进一步验证数值模型的精度，如图5所示，开展节段模型风洞试验，测取了无车的主梁在-5°～+5°风攻角下的静力三分力系数。图6给出了风轴坐标系下三分力系数的模拟值与试验值的对比结果。从图中可以看出，在研究的风攻角范围内，静力三分力系数的模拟值与试验值吻合良好，这验证了数值模拟参数设置的准确性，数值模型可用于列车复气动导纳的识别。
基于试验验证的数值模型，根据表1的工况，图7给出了不同车道处列车静力三分力系数的数值模拟结果，以便于后续复气动导纳的计算。
3 列车复气动导纳识别结果与讨论
3.1 脉动风的影响
以单独迎风侧列车为例，图8和图9分别给出了顺风向和竖向阵风场下列车复气动导纳的识别结果。由于缺乏列车复气动导纳的理论解，为便于比较，将薄机翼的理论解(Sears函数[25])和钝体结构的经验函数(Davenport函数[26]，Vickery函数[27]和Homles函数)也绘于图中。
从图中可知，复气动导纳的幅值和相位角受脉动风分量的影响显著。折算频率较低时，顺风向阻力复气动导纳的幅值大于Davenport函数和Vickery函数，而竖向阻力复气动导纳的幅值则低于Davenport函数。低频处升力复气动导纳的幅值受脉动分量的影响不明显，且小于Sears函数和Holmes函数；当折算频率大于0.05时，竖向升力复气动导纳的幅值大于顺风向升力复气动导纳。可见，列车复气动导纳的幅值与理论解存在较大差异，这是由于来流流经列车表面时发生流动分离和再附着现象[6]，导致列车表面的气动特性不同于势流理论的结果。此外，桥梁的存在影响了列车周围的流场，列车与桥梁间的气动干扰会影响列车表面的压力分布[7]，这是导致列车复气动导纳区别于理论解的另一个原因。
此外，不同脉动风场下，阻力和升力复气动导纳的相位角与理论解存在差异。顺风向阻力复气动导纳的相位角在折算频率大于0.126时低于Sears函数，且随折算频率的增大呈先增大后减小的趋势；竖向阻力复气动导纳的相位角整体上大于Sears函数。随着折算频率的增加，顺风向升力复气动导纳的相位角逐渐减小且低于Sears函数，而竖向升力复气动导纳的相位角呈先增大后减小，在折算频率大于0.126后再增大，且大于Sears函数。上述分析表明，若用薄机翼的理论解模拟列车的抖振力，将无法准确地考虑列车抖振力和阵风场之间的气动传递关系，应将复气动导纳用于列车抖振力的模拟计算。
3.2 列车位置的影响
以顺风向阵风场为例，图10和图11分别给出了列车在不同车道位置时，其阻力和升力复气动导纳的识别结果。从图中可知，由于桥梁与列车之间的气动干扰效应，列车位置对阻力和升力复气动导纳的影响较大。2车交会时，背风侧列车的阻力复气动导纳幅值最大，这与马存明等[6]的研究结论相同，这可能是由于迎风侧列车上表面的流动分离和再附着改变了背风侧列车的周围风场所致，形成了“放大效应”，对行车安全不利；其余工况阻力复气动导纳的幅值和变化规律较一致。单独背风侧列车的升力复气动导纳幅值最大，2车交会时的背风侧列车次之，这与马存明等[6]的试验结果相反，可能是由于数值模拟的风场和列车展长与风洞试验不同导致的；其余工况的升力复气动导纳幅值和变化规律相似。
此外，图10和图11的结果显示，列车阻力和升力复气动导纳间的相位角随列车位置的不同而存在差异，其差异在数值上随折算频率的增大而呈现不同的变化规律。根据严乃杰[3]的研究，列车阻力和升力复气动导纳相位角间的差异，将导致列车抖振升力和抖振阻力互相关系数的降低，进而影响抖振响应的计算结果。以顺风向的计算结果为例，图12给出了顺风向升力和阻力复气动导纳相位差以及抖振阻力和抖振升力的互相关系数，不考虑相位差影响的“准定常理论”结果也绘于图中作比较。从图中可知，复气动导纳的相位差随折算频率的增大呈现先增大后减小的趋势，该趋势因列车位置的不同略有差异；考虑复气动导纳的相位角后，列车抖振阻力和升力的互相关系数整体上低于准定常值，随着折算频率的增大，互相关性的值先增大后减小，而后又增大，这表明列车抖振力模拟中，若忽略复气动导纳的相位角差异，将高估抖振力的互相关性，进而导致偏保守的列车耦合抖振响应结果。
4 适用性分析
复气动导纳属于二维气动导纳，仅与断面形状有关，将其应用于列车抖振力计算时，“片条假设”必须成立。然而，现实中湍流场的三维特性和列车有限展长的特点，使得桥上列车的气动导纳具有三维特性，“片条假设”的适用性有待商榷。
对于湍流场中有限展长的列车，根据抖振力的三维谱张量理论[22]，两波数的抖振阻力功率谱为
式中：
由于列车的
对方程(7)沿展向波数进行积分，得到抖振阻力的点谱表达式：
如果“片条假设”成立，方程(8)中气动导纳的修正项应为1，此时列车的抖振阻力点谱为
根据上述推导过程，“片条假设”的适用性越好，
为量化“片条假设”的适用性，根据LI等[29]的研究，引入无量纲的精度系数
将湍流积分尺度与列车宽度之比
根据文献[8]中列车抖振力的试验数据，基于非线性最小二乘法拟合气动导纳修正项中的系数，拟合时选取了多组初始值，以尽量排除初值对于拟合结果的影响。图13给出了不同折算频率下，精度系数
5 结论
1) 基于列车气动力理论所建立的复气动导纳数值识别方法，能够充分地考虑顺风向和竖向脉动风对列车升力和阻力复气动导纳的影响，识别的复气动导纳能够描述抖振力和阵风场的幅值和相位角在频域的传递关系，可以更准确地刻画湍流作用下桥上列车的抖振力。
2) 由于流动分离和车桥间的气动干扰，列车的复气动导纳不同于薄机翼理论结果，其幅值和相位角受脉动风分量和列车位置的影响较大，2车交会时，背风侧列车的顺风向阻力复气动导纳的幅值最大，单独背风侧列车的顺风向升力复气动导纳最大，在列车耦合抖振响应的计算中应加以重视。
3) 列车阻力和升力复气动导纳的相位角存在差异，这种差异随折算频率的增大而变化。考虑复气动导纳的相位角后，随着折算频率的增大，顺风向列车抖振阻力和抖振升力的互相关系数呈先增大后减小，而后又增大的趋势，且整体上低于不考虑相位差影响的准定常理论结果。
4) 列车复气动导纳的幅值受湍流积分尺度和特征宽度的比值
未来的研究将考虑把复气动导纳集成到列车-桥梁耦合振动模型，以实现列车-桥梁耦合抖振响应的精细化分析。此外，本文移动列车复气动导纳的数值识别方法理论上可以推广至任意主梁断面的桥上列车，未来可针对不同的桥梁类型(如悬索桥、斜拉桥)和主梁断面(如流线型箱梁、分体式双箱梁等)开展桥上列车复气动导纳的数值识别，以验证该方法的普适性。
分体式三箱梁的二维气动导纳研究
[J/OL]. 铁道标准设计, 2024: 1-11. (邰学洋,郑史雄,严正曦等.阵风场中桥上移动列车复气动导纳的数值识别[J].铁道科学与工程学报,2025,22(10):4726-4736.
TAI Xueyang,ZHENG Shixiong,YAN Zhengxi,et al.Numerical identification of complex aerodynamic admittance of vehicles on bridge under gusts[J].Journal of Railway Science and Engineering,2025,22(10):4726-4736.