基于PIM法车桥耦合模型Duhamel项迭代格式选择

桂水荣， 万 水， 陈水生

(1.东南大学交通学院 南京, 210096) (2.华东交通大学土木建筑学院 南昌,330013)

10.16450/j.cnki.issn.1004-6801.2017.05.007

基于PIM法车桥耦合模型Duhamel项迭代格式选择

桂水荣1,2， 万 水1， 陈水生3

(1.东南大学交通学院 南京, 210096) (2.华东交通大学土木建筑学院 南昌,330013)

考虑荷载在积分步长内的协调分解和不变常量，结合模态综合叠加技术，建立移动弹簧质量车桥耦合振动模型，引入精细积分算法(precise intergration method,简称PIM)，并采用科茨及高斯积分格式展开精细积分中Duhamel非齐次项进行求解。以移动弹簧质量车模型作用于简支梁桥为例，分析积分步内荷载等效方法及Duhamel非齐次项展开方式对车桥耦合振动响应的影响。研究结果表明：荷载协调分解并结合科茨积分格式计算结果与解析解更接近，Newmark-β法积分步长内荷载分解形式对计算结果影响较小；荷载协调分解的高斯积分格式计算结果呈发散趋势，待定系数法计算结果与实际偏离较大；Newmark-β法能满足工程需要，但在保证相同计算精度时，需采用较小积分步长、同时耗费较多计算时间；与其他几种数值方法相比，积分步长内荷载协调分解并将科茨积分格式引入精细积分的算法(PIM-Cotes-Harmonize,简称PIM-C-H法)，具有快速收敛且计算快的优势。

车桥耦合振动模型； 精细积分法； Duhamel项； 迭代格式选择

引 言

公路桥梁与汽车荷载之间的相互作用，近些年一直是国内外学者研究桥梁问题的一个热点。数值模拟车桥耦合振动响应，运用模态综合叠加技术并结合隐式Newmark-β算法求解，一直是国内外学者广泛采用的方法[1]，分析中需通过保证积分步长来确保计算结果的收敛性[2]。钟万勰[3]提出的精细积分法，能给出齐次结构动力响应的高精度解答。对于非齐次结构动力方程引起的Duhamel非齐次项求解，因非齐次项的求解方式直接影响了非齐次动力方程计算精度，国内许多学者对此展开研究。谭述君等[4]利用Duhamel积分矩阵，导出非线性方程求解一般格式，结合Adams线性多步法，构造基于精细积分方法的算法。高强等[5]利用大规模动力系统矩阵的稀疏性和动力问题的物理特性，提出一种改进的快速精细积分方法。张文志等[6]将初值问题转化为边值问题，利用快速傅里叶变换(fast Fourier trans formation,简称FFT)进行求解。汪梦甫等[7]将高斯积分方法与精细积分法中的指数矩阵相结合，求解非齐次方程的特解，计算时不需对矩阵求逆。富明慧等[8]提出了一种广义精细积分方法，对非齐次项为多项式和指数函数，能实现特解的精细积分，但对任意形式荷载并不适用。同时，储德文等[9]对比几种数值算法指出，非齐次方程的特解求解需根据荷载的性质选择合适的积分方法。

根据非齐次方程的荷载形式，黄涛等[10]提出了一种求解非线性动力学方程的预测-校正数值算法。结合车桥耦合振动特点，杜宪亭等[11]构造了一种改进的高斯精细积分算法用于求解结构非线性问题。运用精细积分法求解非线性方程，引入指数矩阵后，加大积分步长并不影响齐次方程计算精度，但对于行驶在公路桥梁上的汽车荷载，积分步长内荷载作用位置不断发生变化，荷载项为时间、系统运动的非线性隐性函数。积分步长内荷载的等效方法及非齐次项的积分格式，对计算结果将有怎样的影响，需进一步研究。

根据车桥耦合振动特点，结合模态综合叠加技术，笔者建立移动弹簧质量车桥耦合振动模型，考虑节点荷载在积分步长内协调分解和不变常量两种格式，运用精细积分算法并结合科茨及高斯积分格式进行求解，研究移动弹簧质量车模型作用下非齐次项迭代格式对车桥耦合振动响应的影响。

1 车桥耦合精细积分模型

1.1 车桥耦合振动模型

将车辆模型简化为由两个质量体系组成的移动弹簧质量车模型。假定车辆与桥梁始终保持接触，车辆模型的质量由车轮质量m2及车体的簧载质量m1组成，车轮和悬架系统的弹簧刚度及阻尼分别等效为刚度k1、阻尼c1的弹簧-阻尼系统。假设简支梁桥静止时为平衡位置，车辆以速度v行驶，梁的动挠度为y(x,t)，簧载质量m1的动位移为z(t)，车轮m2始终与梁体保持接触不脱离，且位移与所在位置的桥梁挠度y(x,t)一致，移动弹簧质量车模型如图1所示。

图1 移动弹簧质量车桥耦合振动模型Fig.1 A beam subjected to a moving spring mass system

根据达朗贝尔原理，移动弹簧质量车辆模型振动方程为

(1)

车辆振动过程中，作用于桥梁的力为

(2)

运用有限元方法，将桥梁离散成欧拉梁单元，桥梁振动方程可以表示为

(3)

其中:Mb,Cb,Kb分别为桥梁系统的质量、阻尼和刚度矩阵；u为桥梁单元节点向量；Fb为作用于桥梁上的外荷载。

采用Rayleigh阻尼计算桥梁结构阻尼Cb=αMb+βKb，其中α和β为Rayleigh阻尼系数。桥梁节点模态空间取r阶，根据模态振型的正交性，不考虑各阶模态的相干性，经振型正交分解， 式(3)可以改写为

(4)

(5)

其中：Mbv，Cbv，Kbv为车桥耦合振动广义刚度、阻尼及质量矩阵；ybv为车桥耦合振动方程坐标向量；ybv=[q1,q2,…,qr,z]T；Fbv为车桥耦合振动广义荷载向量。

(6)

其中：N为单元插值函数；Q，S，A1，A2，B1，B2分别为车桥耦合振动引起的阻尼及刚度矩阵修正系数矩阵；MV，CV，KV分别为车辆模型质量、阻尼、刚度矩阵。

1.2 移动荷载节点等效分解

1.2.1 积分步长内荷载协调分解

如图2所示，假设简支梁第i单元长度为li，移动弹簧质量荷载在i单元内，自i端向j端以速度v匀速行驶。某一积分步长Δt内，在τ时刻弹簧质量系统由A点移动到B点，此时A，B两点均在i单元内。数值迭代时，任意tk时刻作用于A点，而tk+1时刻移动至B点，A点和B点分别距i单元起点的距离为x1和x2。假设t时刻，移动弹簧质量车模型距i单元起点距离为x，有

x=vtk-xk，ξ=x/l=ξ1t-ξ2

(7)

图2 移动弹簧质量车模型在i单元内移动Fig.2 A moving spring mass system moving in the i-th element

图3 i单元内，荷载节点等效分解Fig.3 A moving spring mass system moving in the i-th element

(8)

其中:Ni(ξ)(i=1～4)为梁单元的插值函数。

此时作用于桥梁上的外荷载等效到单元节点荷载可以表示成

(9)

将式(7),(8)代入(9)，并整理可得

Fv=(r0+r1t+r2t2+r3t3)[-ΦT(m1+m2)g]

(10)

其中:ri(i=0～3)为外荷载的常系数向量。

1.2.2 积分步长内荷载为常量

如图2所示，假设积分步长内荷载为常量，此时单元内荷载等效到节点荷载，积分步长内等效节点荷载为不变常量，任意tk时刻，节点i和节点j的等效节点荷载相等，式(8)改写为

(11)

此时，式(11)在积分步长内，荷载为定值，求解算法与Newmark-β法求解车桥耦合振动算法相同[12]。

1.3 精细积分算法分析

车桥耦合振动模型的二阶非齐次微分方程，运用精细积分法求解其通解，并对Duhamel非齐次项采用数值迭代积分格式分析。

将式(5)转化成状态空间方程

(12)

(13)

式(12)的解可以写成如下形式

(14)

将计算时间间隔Δt划分成若干步，则可以得到t+Δt时刻与t时刻系统响应之间关系

(15)

其中:exp为指数矩阵，第1项为通解，采用精细积分算法求解；第2项为非奇次方程的特解，Duhamel非齐次项运用科茨积分法及高斯积分法，分别进行数值迭代求解。

2 PIM法Duhamel项迭代格式分析

2.1 科茨积分格式

式(15)Duhamel非齐次项的解，可以采用科茨积分格式进行求解。对式(15)进行离散，令τ=tk+1-tk，τ为积分步长，则

(16)

2.2 高斯积分法

采用高斯积分公式展开式(15)中Duhamel项[11]，得

vt+Δt=exp(HΔt)+Δt/2×

(17)

其中：m为高斯积分点个数；θi,wi分别为高斯积分点位置、权重系数。

随着m增大积分精度提高，计算量随之增大。结合有限元等参元高斯积分的计算经验[11]，取m=3。式(17)中参数按如下参数取值

w1=5/9；w2=8/9；w3=5/9

(18)

式(17)离散后的非线性荷载项涉及t+θiΔt时刻的系统运动状态，可以通过以下方法确定。假定t到t+Δt之间任一时刻τ的位移，可由t,t+Δt时间步的位移、速度依据三次Hermite函数插值唯一确定

(19)

其中:N1(θ)=1-3θ2+2θ3；N2(θ)=3θ2-2θ3；N3(θ)=θ-2θ2+θ3；N4(θ)=θ3-θ2；θ=(τ-t)/Δt；0≤θ≤1。

(20)

3 算例分析

例1为研究PIM数值积分法计算车桥耦合振动响应的有效性及积分步长内荷载分解方法对车桥耦合振动响应的影响，如图4所示，以某一简支梁桥受移动常量力为例，对比分析几种常用数值算法、PIM数值算法与解析解的误差。图中：m为单位长度质量；EI为简支梁桥抗弯刚度。

图4 移动常量力作用于简支梁桥示意图Fig.4 Diagram of simple support beam subjected to moving constant force

定义ζ=T1/t，t=L/v，μ=(Rd-Rj)/Rj，其中：T1为结构的第1阶固有振动周期；t为移动荷载通过桥梁的时间；L为桥梁的计算跨径；Rd为跨中振动响应最大值；Rj为跨中最大静态位移；μ为位移冲击系数。

为对比积分步长内荷载等效分解方法对计算结果的影响，分析移动常量力以不同的ζ值通过简支梁时，跨中位移冲击系数，分别采用解析法、Runge-kutta法(解析方程数值解)、Newmark-β法(全自由度和模态综合叠加两种方法)、精细积分法(积分步长内荷载的协调分解及不变常量)求解。简支梁桥结构参数，L=30m，ρ=2 600kg/m3，A=1.062 2m2，I=0.509 2m4，E=3.5·1010N/m2，移动常量力P=3.278·1010N。

表1列出了PIM法与Newmark-β法计算跨中节点位移冲击系数与参考文献对比结果。其中：exact[1]结果为简支梁取前10阶模态的解析解；full-size[13]结果是采用全自由度直接耦合方法计算所得；而Newmark-β[12]结果是利用ANSYS将简支梁离散成有限元模型，提取前十阶模态向量计算所得；PIM-Cotes-Harmonize(PIM-C-H)法，荷载在积分步长内等效协调分解(Harmonize)，采用精细积分法(PIM)，结合科茨积分格式(Cotes)求解；PIM-Cotes-Same(PIM-C-S)法，积分步长内荷载为不变常量，采用精细积分法，结合科茨积分格式；PIM-Gauss-Same(PIM-G-S)法，积分步长内荷载为不变常量，采用精细积分法，结合高斯积分格式(Gauss)求解。

从表1可以看出，模态综合叠加法并结合精细积分算法，计算简支梁受移动常量力作用，跨中位移冲击系数与解析解最大相对误差为1.96%；当ζ=1.5时，各种数值算法计算的冲击系数均与文献[1]有较大偏差，其他ζ值下，各种算法的冲击系数最大绝对误差值为0.005 8；精细积分法的积分步长内荷载协调分解格式与积分步长内常量荷载迭代格式的误差最大值为0.05%。积分步长内荷载为不变常量，PIM-C-H计算结果比Newmark-β法更接近解析值，且不受积分步长限制，能更快收敛。

例2以图1所示移动弹簧质量车模型行驶在简支梁桥上为例，分析科茨积分、待定系数法及高斯积分三种数值迭代格式求解车桥耦合振动模型Duhamel非齐次项响应的差异。简支梁桥取例1中桥梁参数，车辆参数取m1=1 425 kg，m2=32 025 kg，k1=6.5×105N/m，c1=2.1×104(N·s)/m。

图5为采用PIM-C-H，PIM-C-S，Runge-Kutta法和Newmark-β法分别计算简支梁跨中位移响应。从图5中可以看出，Runge-Kutta法计算结果与PIM法及Newmark-β法偏差较大；积分步长较小时， PIM-C-H法与Newmark-β法曲线重合，PIM-C-H法和PIM-C-S存在一定误差，精细积分结合科茨积分格式，荷载分解方式对计算结果产生一定影响， Runge-Kutta法计算结果偏离其他三者较大。图6为分别采用PIM-C-H法、PIM-C-S法、PIM-G-S法求解车桥耦合振动响应曲线。图7为采用PIM-G-H法求解车桥耦合振动响应。从图6和图7中可以看出，积分步长内荷载为不变常量，精细积分法的特解采用Gauss或Cotes迭代格式，两者计算结果完全吻合；但考虑积分步长内荷载协调分解，Gauss算法发散，文献[11]针对车桥耦合的Gauss算法提出一种预测-校正的精细积分-高斯积分格式，三节点高斯积分格式求解任意荷载形式的车桥耦合振动特解需进行校正。采用待定系数法计算车桥耦合振动响应，因待定系数法需对矩阵求逆，计算结果虽不发散，但曲线明显与实际不符。

表1 移动集中力通过简支梁桥，跨中位移冲击系数对比表

Tab.1 Contrast table of Impact factors of the simple supported beam with mobile concentration force

T1/t①exact[1]②lin[14]③full-size[13]④Runge-Kuta[12]⑤Newmark-β[12]⑥PIM-C-S⑦PIM-C-H(⑦-①)/①(%)0.10.0500.0530.0650.0480.0510.0510.0511.600.50.2500.2520.2560.2550.2560.2550.2551.961.00.7070.7050.7030.7050.7020.7030.703-0.521.2340.7430.7300.7370.7320.7320.7330.733-1.401.50.7100.7040.6950.7020.6990.6990.699-1.512.00.5500.5500.5490.5500.5500.5500.550-0.04

(⑦-①)/①%为PIM-C-H法计算结果与精确值误差。

图5 不同算法分析跨中位移响应Fig.5 Displacement curve of dynamic response at mid span using different numerical algorithm

图6 精细积分法不同特解迭代格式跨中位移响应Fig.6 Dynamic displacement response of at mid span using different PIM integral formulation

根据上文分析，Newmark-β法、PIM-C-S法、PIM-G-S法、PIM-C-H法能有效计算车桥耦合振动响应，表2分析了这4种方法计算移动弹簧质量车模型行驶在简支梁桥上，积分步长对程序计算时间、收敛速度及跨中最大动挠度的影响。从表2可以看出，相同积分步长，Newmark-β法计算精度远小于PIM-C-H法，但相同积分次数，其运行时间较PIM-C-H法少；在保证相同计算精度条件下，PIM-C-H法能快速稳定。精细积分法，在单元节点荷载等效时，若不考虑荷载在积分步长内协调分解，采用高斯积分格式和科茨积分格式，计算误差均较Newmark-β法大，但考虑积分步长内荷载按协调分解，误差最小。同时，从解的稳定性来看，相同行车速度下，Newmark-β法需PIM-C-H法近10倍的积分步、3倍的计算时间方可达到相同的稳定解，PIM-C-H具有快速收敛，但相同积分步长，计算时间较长。

图7 PIM-G-H法分析跨中位移响应Fig.7 Dynamic displacement response of at mid span using PIM-G-H method

车速/(m·s-1)①Newmark-β法②PIM-Gause-Same③PIM-Cotes-Same④PIM-Cotes-Harmonize总积分步最大动位移/mm计算时间/s总积分步最大动位移/mm计算时间/s总积分步最大动位移/mm计算时间/s总积分步最大动位移/mm计算时间/s530060012006000120006000011.055611.255111.323311.343911.344211.34440.132760.263770.523492.533455.0332324.968523006001200600010.770610.778510.779210.78010.466340.903471.766798.729323006001200600010.764910.778410.779210.78010.447050.853221.715318.550273006001200600011.151711.337411.345011.34440.504381.007881.962179.7721015200400200040002000010.853910.863710.874410.874810.87500.097730.177870.848261.678378.203842004002000400010.690910.689310.691010.69100.300760.589762.966645.842892004002000400010.690910.689310.691010.69100.302890.585632.879125.66656200400200010.876210.874310.87500.346480.666613.3646530100200100020001000011.305511.250011.254511.253811.25360.055010.101210.434340.849154.234131002001000200011.121311.126111.128511.12850.155880.296851.479052.957111002001000200011.121311.126111.128511.12850.154360.294401.417612.931651002001000200011.254711.258511.253711.25560.179260.343891.685213.31891

4 结束语

结合模态综合叠加技术，建立了移动弹簧质量车桥耦合振动模型，考虑积分步长内荷载按协调分解和不变常量两种格式，等效引入精细积分法，研究Duhamel非齐次项积分格式对简支梁跨中位移响应的影响。研究结果表明：

1) 精细积分结合科茨积分格式计算移动常量力作用下车桥耦合振动响应，协调分解值与解析值更接近；积分步长较小时，Newmark-β法积分步内荷载可按常量处理。

2) 精细积分法结合高斯积分考虑荷载协调分解，求解车桥耦合振动响应呈发散趋势。

3) Newmark-β法能满足工程需要，但为保证其计算精度，需较小的积分步长、且耗费较多的计算时间。

4) 精细积分法求解车桥耦合振动响应，采用科茨积分格式，计算结果不受积分步长的影响，考虑积分步内荷载协调分解，计算精度高，PIM-C-H法具有快速收敛且计算较快的优势。

[1] Henchi K, Fafard M,Talbot M,et al. An efficient algorithm for dynamic analysis of bridges under moving vehicles using a coupled modal and physical components approach[J]. Journal of Sound and Vibration, 1998 ,212(4):663-683.

[2] 翟婉明. 车辆-轨道耦合动力学[M].北京：科学出版社，2007:1-8.

[3] 钟万勰. 结构动力方程的精细时程积分法[J].大连理工大学学报，1994, 32(2):131-136.

Zhong Wanxie. On precise time-integration method for structural dynamics [J]. Journal of Dalian University of Technology, 1994, 32(2):131-136.(in Chinese)

[4] 谭述君, 高强, 钟万勰. Duhamel 项的精细积分方法在非线性微分方程数值求解中的应用[J].计算力学学报,2010, 27(5):752-758.

Tan Shujun, Gao Qiang, Zhong Wanxie. Application of Duhamel term′s precise integration method in solving nonlinear differential equations[J]. Chinese Journal of Computational Mechanics, 2010, 27(5):752-758. (in Chinese)

[5] 高强, 吴锋, 张洪武,等. 大规模动力系统改进的快速精细积分方法[J].计算力学学报, 2011，28(4):493-498.

Gao Qiang, Wu Feng, Zhang Hongwu, et al. A fast precise intergration method for large-scale dynamic structures[J].Chinese Journal of Computational Mechanics, 2011，28(4):493-498. (in Chinese)

[6] 张文志, 富明慧, 刘祚秋. 结构动力方程的精细积分-FFT方法[J].中山大学学报：自然科学版,2008, 47(6):12-15.

Zhang Wenzhi, Fu Minghui, Liu Zuoqiu. Precise integration-FFT method for structural dynamics[J]. ACTA Scientiarum Naturalium Universitatis Sunyatseni, 2008,47(6):12-15. (in Chinese)

[7] 汪梦甫, 周锡元. 结构动力方程的高斯精细时程积分法[J]. 工程力学, 2004, 21(4):13-16.

Wang Mengfu, Zhou Xiyuan. Gauss precise time-integration of structural dynamic analysis[J]. Engineering Mechanics, 2004, 21(4):13-16. (in Chinese)

[8] 富明慧，刘祚秋, 林敬华. 一种广义精细积分法[J]. 力学学报, 2007, 39(5):72-77.

Fu Minghui, Liu Zuoqiu, Lin Jinhua. A generalized precise time step integration method[J]. Chinese Journal of Theoretical and Applied Mechanics, 2007, 39(5):72-77. (in Chinese)

[9] 储德文, 王元丰. 精细直接积分法的积分方法选择[J]. 工程力学, 2002，19(6):115-119.

Chu Dewen, Wang Yuanfeng. Integration formula selection for precise direct integration method[J]. Engineeying Mechanics, 2002,19(6): 115-119. (in Chinese)

[10] 黄涛, 樊建平, 何建平,等. 求解非线性动力学方程的预测-校正数值算法[J]. 固体力学学报, 2007,28(1):1-6.

Huang Tao, Fan Jianping, He Jianping, et al. A predict-correct numerical integration scheme for solving nonlinear dynamic equations[J]. Acta Mech Anica Solida Sinica, 2007,28(1):1-6. (in Chinese)

[11] 杜宪亭, 夏禾, 李慧乐,等. 基于改进高斯精细积分法的车桥耦合振动分析框架[J]. 工程力学, 2013，30(9):171-176.

Du Xianting, Xia He, Li Huile, et al. Dynamic analysis framework of train-bridge system based on improved gauss precise integration method[J]. Engineering Mechanics, 2013，30(9):171-176. (in Chinese)

[12] 王运金，桂水荣，陈水生.连续梁桥车桥耦合振动分析的数值解法[J]. 华东交通大学学报，2007,24(4):25-29.

Wang Yunjin，Gui Shuirong，Chen Shuisheng. An efficient algorithm for coupled vibration analysis of continuous bridges under moving vehicle[J]. Journal of East China Jiaotong University，2007,24(4):25-29. (in Chinese)

[13] 陈水生. 公路车桥耦合振动响应计算方法对比研究[J].华东交通大学学报，2011, 28(3):18-25.

Chen Shuisheng. Comparative research on the calculating methods for coupling highway vehicle-bridge system[J]. Journal of East China Jiaotong University, 2011, 28(3):18-25. (in Chinese)

[14] Lin Y H,Trethewey M W. Finite element analysis of elastic beams subjected to moving dynamic loads[J]. Journal of Sound and Vibration, 1990,136(2):323-342.

国家自然科学基金资助项目(51468018 ,51268013)；江西省教育厅科研资助项目(GJJ14384)；教育部工程研究中心资助项目(13TM02)

2015-07-17；

2016-04-18

U441.3; TH.3

桂水荣，女，1979年10月生，博士生、副教授。主要研究方向为公路桥梁车桥耦合振动。曾发表《系杆拱桥吊杆节点足尺模型承载性能试验研究》(《桥梁建设》2013年第43卷第2期)等论文。

E-mail:guishuirong@163.com