第18卷 第3期 重 庆 交 通 学 院 学 报1999年9月Vol.18 No.3 JOURNALOFCHONGQINGJIAOTONGINSTITUTESep.1999
文章编号:10012716(1999)0320014207
车桥耦合振动分析的数值方法Ξ单德山,李 乔(西南交通大学土木工程学院桥梁及结构工程系,四川成都610031)
摘要:车桥耦合振动问题是铁路和公路桥梁中十分重要的研究课题,而目前所采用的数值算法所需的时间比较长,为了减少计算机时,本文在对高速铁路曲线梁车2桥耦合振动研究中,建立了一种基于激励非线性振动的数值计算方法,并完成了计算程序BSNDS的编制,取得了较好的计算结果.并将其与其他模型进行比较,在保证精度的前提下,较大地节省了计算时间.
关 键 词:结构工程;耦合振动;数值方法中图分类号:U443234 文献标识码:A
对于车桥耦合振动分析这一类复杂问题,常用的算法有两种:时间序列的逐步积分法和频响函数法.时间序列的逐步积分法是将车辆和桥梁看作一个大的振动系统,建立该系统的运动微分方程并用直接积分法求解,得到各自由度上的位移、速度和加速度的时程[1];频响函数法是基于
随机振动的一种方法,该方法首先计算出车桥耦合系统的频响函数,用激励力的功率谱作为输入,求得系统在频域的响应[2].本文所介绍的方法是基于激励非线性振动的一种逐步积分法,在
计算中应用了求解非线性振动的Newmark预测2校正法[9],即在每一时段里预测桥梁的位移、速
度、加速度和车桥系统的耦合力,此时车桥系统的位移条件是协调的,以此作为迭代的开始进行计算,从而减少了迭代次数,进而减少了计算机时.
车桥耦合振动分析的困难在于寻找一种能处理车桥运动耦合的方法.在接触点处采用常规的运动方程的形式来描述车2桥系统的耦合振动M¨W+CW+KW=fcp(1)
式中,桥的特性由M(质量阵)、C(阻尼阵)、K(刚度阵)和W(位移)来描述.位移函数W是在t时刻接触点的位移;W、¨W分别表示其速度和加速度;点号(・)表示对时间求导;式(1)中fcp表示车桥间的耦合力,它可以看成是由桥上移动的车辆所施加的力.fcp是车辆运动的函数,它
还与桥梁的振动和路线的不平顺有关,这种相互关联的运动称为车2桥系统的运动耦合.当t时刻有两个或更多的车辆在桥上时,耦合力fcp还与桥上其它车辆有关.车与车之间的耦合通过桥
Ξ收稿日期:1998211220
基金项目:铁道部科技开发研究项目97G07
作者简介:单德山(1969-),男,四川大竹县人,西南交通大学讲师(博士),从事的研究是结构的空间行为.
© 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.梁的运动而包含在耦合力fcp中,这样就可考虑同时通过桥梁的所有车辆的影响.然而在每一时间步里,fcp事先并不知道,且是与未知量W及其导数、车的特性等有关的函数.显然,考虑运动耦合时,耦合力fcp是非线性的,式(1)为非线性运动方程,因而系统的响应也为非线性.如果忽略耦合运动,fcp仅是一与时间有关的荷载,此时该问题简化为移动荷载作用下的振动问题.
1 系统运动微分方程的离散桥梁结构承受动荷载的微分方程为:
M¨D+CD+KD=F(2)式中,M、C、K分别为桥梁的质量、阻尼和刚度矩阵;D=D(t)为离散系统的结点位移向量;
D和¨D分别为结点的速度和加速度向量;F为作用在结构上的力矢量,它是时间t的函数,F
主要由耦合力、车的重量和桥的重量等组成.
桥梁单元的质量和刚度矩阵是根据考虑剪切变形的空间直梁和考虑翘曲的曲线梁而推导出来的[3]、[6].空间直线梁单元由6个自由度组成,即由3个线位移和3个角位移组成;空间曲线梁单
元由7个自由度组成,即由3个线位移、3个角位移和截面翘曲位移组成.根据有限元法,将桥梁单元的有限单元矩阵和车桥耦合力的等效节点力组装起来就可得到公式(2)的离散矩阵方程.
如果单元的局部坐标与总体坐标方向不一致,应进行坐标变换.在车桥耦合振动中,还应考虑系统阻尼的影响;对于桥梁单元仍采用Raleigh阻尼,车辆的阻尼由其阻尼元件提供.
2 车桥耦合振动求解的数值算法为了能够求解方程(2),应在方程(2)中引入边界条件(包括位移边界条件和力的边界条件)
.
此时满足方程(2)的位移矢量D=D(t)即为方程的解.给定的初始条件为:D(0)=u
0;D
(0)=
u
0;u0和u0分别为在t=0时刻位移和速度的初始值
.
211 系统方程的时间离散根据Newmark法的要求,方程(2)尚需对时间进行离散.实际上时间离散仅是一个代数问题.
为了反映控制微分方程的这个特征,在以后的描述中桥梁的位移D(t)、速度D(t)、加速度¨D(t)
分别用Dn、Vn、An来表示;tn为第n时段(1≤n≤N);N为时间段离散的数目.同样dn、vn、an分别表示tn时段车辆的位移、速度、加速度,即dn=d(tn);vn=d(tn);an=¨d(tn).当求解下一时段时,即n+1时段车桥运动方程按如下表示:
MAn+1+CVn+1+KDn+1=Fn+1(3)
式中,Dn+1、Vn+1、An+1为t=tn+1时桥梁的位移、速度、加速度;F
n+1=Fcpn+1+Fspn+1,Fspn+1
包括了桥的恒载及规范所规定作用的荷载.Fcpn+1为车2桥间的耦合力,对于一维的运动质量(图
1),t
n+1
时刻的耦合力为:
Fcpn+1=k1(dcp-dr1)n+1+c1(vcp-vr1)n+1(4)
式中,dcp、vcp为接触点的位移和速度;vr1、dr1为车辆簧下质量ml的位移和速度;k1、c1为车轮的弹簧刚度和阻尼(图1).由于车在接触点的dcp、vcp、acp是与桥梁的运动有关的,那么矢量
51第3期 单德山等:车桥耦合振动分析的数值方法© 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.图1 一维运动质量可用以下的符号表示:
Fn+1=Fspn+1+Fcpn+1(Dcpn+1,Vcpn+1,Acpn+1,dn+1,vn+1,an+1)(5)
应用Newmark法,可将tn+1时刻的位移和速度表示如下:
Dn+1=Dn+ΔtVn+12Δt2[(1-2β)An+2βAn+1]Vn+1=Vn+Δt[(1-γ)An+γAn+1](6)
式中,Δt=t
n+1-tn;β,γ为保证计算精度和积分稳定性的控制参数.
将式(6)代入式(3),并用符号表示有:
MAn+1+CVn+1(An+1)+KDn+1(An+1)=Fn+1(An+1,dn+1,vn+1,an+1)
(7) 式(7)中并没有明显地含前一时段tn所确定的带有下标n的项.式(7)所表示的车桥耦合振动为非线性振动,因而最好采用迭代法求解.在这里采用一种与Newmark法相关的多次预测2校正过程进行求解.在求解过程中,当前迭代步的校正值作为下一迭代步的预测值,那么预测2校正的过程将不断进行直到达到所需的精度要求.
在tn+1时段,为将多次预测2校正过程初始化,必须假设式(7)中的初值A
n+1.
该值可设为
A0n+1且等于0,上标0表示预测2校正的迭代次数i=0.那么由式(6)可得
D0n+1=Dn+ΔtVn+12Δt2(1-2β)An V0n+1=Vn+Δt(1-γ)An(8)
由式(8)可知,A0n+1为0时,V0n+1和D0n+1由前一时段的值完全确定.已知A0n+1、V0n+1和D0n+1后,就可确定式(7)中的力矢量.至此可以开始预测2校正循环.而从计算实现的角度来看,
更新迭代次数i是非常方便的.再一次强调,车2桥系统的初始预测值不仅包括了A0n+1、V
0
n+1
和D0n+1,而且包括了F0n+1,对应于Ai-1n+1、Vi-1n+1和Di-1n+1的力矢量为F
i-1
n+1.如式(5)所示的那
样,Fi-1n+1是耦合力(Fcpn+1)i-1的函数,根据预测的桥梁位移和速度值,按如下方法确定F
i-1
n+1.
已知桥梁的Vi-1n+1和Di-1n+1,那么车在桥上的接触点在tn+1时段的dcp、vcp可通过包含接触点的桥梁单元的位移插值函数和截面的几何形状来得到.若要考虑到桥面的不平顺,那么不平顺的影响也应加到dcp、vcp中去.已知了dcp、vcp,那么车辆各自由度在tn+1时段的位移、速度和加速度可由下式得到:
mai-1n+1+cvi-1n+1+kdi-1n+1=fi-1n+1(9)
式(9)中,m、c和k表示车辆的质量矩阵、阻尼矩阵和刚度矩阵,它们的具体表达式祥见文[5]、[6].根据车辆动力学理论,求解方程(9)即可得到接触点处的耦合力(fcpn+1)i-1,随后根据有限元法得到耦合力的等效节点力,并由式(5)得到荷载力矢量,代入式(7)中即可进行求解.为了不中断当前的讨论,确定的dcp、vcp、dr、vr过程将在下一节中讨论.
完成了初始预测步骤后,就可进入下一相应的校正步骤.一般来说,预测的Ai-1n+1、Vi-1n+1、
Di-1n+1和力矢量Fi-1n+1并不能满足方程(7).对应于时段tn+1的Ai-1n+1、Vi-1n+1、Di-1n+1的力矢量
Fi-1n+1,系统残余力矢量的表达式为:ΔF=Fi-1n+1-(MAi-1n+1+CVi-1n+1+KDi-
1
n+1)(10)
61重庆交通学院学报 第18卷© 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.