车辆系统动力学计算机仿真 - 2014
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
END
一、无外加激扰（自由衰减振动）
M1加速度
10 5 0 -5 -10
k m3
2l m1 Iy c z3 m2 z2 z1
Acceleration m/s
2
0
1
2
Time/s
3
4
M1速度
m1
2l Iy c m3 z3 m2 z2 z1
8 4
k
Velocity m/s
0 -4 -8 0.0
DIMENSION DY(ND),DV(ND),DA(ND),D1A(ND),CY(ND),CV(ND)
VST2=VST*VST DO 700 I=1,ND,1
CY(I)=DY(I)+DV(I)*VST+(F1KC*DA(I)-FKC*D1A(I))*VST2
CV(I)=DV(I)+(FI1*DA(I)-FI*D1A(I))*VST 700 CONTINUE RETURN END
END
SUBROUTINE EXCH2(ZZC,VZC,AZC,DZZC,DVZC,DAZC,D1AZC) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DZZC=ZZC DVZC=VZC D1AZC=DAZC DAZC=AZC
ZZC=0.
VZC=0. AZC=0.
RETURN
车辆系统动力学计算机仿真
假设Mb为车体质量，Mc 为附加质量块，列车运行 速度为V，如何设置相关 参数使得车体振动完全静 止？
本题实际上是一个反共振问题，动力吸振器包括质量
m2 、弹簧 k 2 ， m1 与 k1 是主系统的质量及弹簧刚度，偏
心质量引起的机械垂向振动力为 m 2 e sin(t ) 。 由此 不难得到系统的强迫振动方程为：
zt2=0.05*sin(4.*pi*tt)
zt3=0.05*sin(4.*pi*tt) CALL EXP1(dzb1,dvb1,dab1,d1ab1,zb1,vb1,nd) ! 积分
CALL EXP2(dzb2,dvb2,dab2,d1ab2,zb2,vb2)
CALL EXP2(dzb3,dvb3,dab3,d1ab3,zb3,vb3)
COMMON /EXPL/ step,F1KC,FKC,FI1,FI dimension zb1(nd),vb1(nd),ab1(nd)
dimension dzb1(nd)源自文库dvb1(nd),dab1(nd),d1ab1(nd)
open(11,file='rusta.dat',status='unknown') open(12,file='rustd.dat',status='unknown')
cmk=1.0e+6;dmp=2.0e+4
cbk=1.0e+7;dbp=2.0e+4
bm1=3250.0 bm2=2100.0 bim=3020.0 bl=2.5
F1KC=1.0 FKC=0.5 fi1=1.5 fi=0.5 step=1.0e-4
do i=1,nn,1 tt=dble(i)*step
k
m1 m3
Iy c z3 m2 z2
z1
20 10
Y Axis Title
0 -10 -20
0
1
2
3
4
5
6
加速度（M2）
k
2l m1 m3 Iy c z3 m2 z2 z1
40 20
Y Axis Title
0 -20 -40
0
2
4
6
M1位移
k
2l m1 m3 Iy c z3 m2 z2 z1
! 积分
! 积分
! 求弹簧悬挂点质量 1的位移与速度
dm2=zb1(1)+0.5*bl*zb1(2) vm2=vb1(1)+0.5*bl*vb1(2) dm3=zb1(1)-0.5*bl*zb1(2) vm3=vb1(1)-0.5*bl*vb1(2)
fm2=cmk*(dm2-zb2)+dmp*(vm2-vb2) fm3=cmk*(dm3-zb3)+dmp*(vm3-vb3)
1
2
3
4
5
6
加速度（质量1）
40 20
Y Axis Title
0 -20 -40
0
1
2
3
4
5
6
加速度（质量2）
80 40
Y Axis Title
0 -40 -80
0
1
2
3
4
5
6
加速度（质量3）
100 50
Y Axis Title
0 -50 -100
0
1
2
3
4
5
6
位移（质量1）
60 30
Y Axis Title
fb2=cbk*(zb2-zt2)+dbp*vb2
fb3=cbk*(zb3-zt3)+dbp*vb3
ab1(1)=(-fm2-fm3)/bm1-gg
ab1(2)=0.5*bl*(fm3-fm2)/bim
ab2=(-fb2+fm2)/bm2-gg
ab3=(-fb3+fm3)/bm2-gg
if(mod(nn,800).eq.0)then write(11,100) tt,ab1(1),ab2,ab3,zt2 write(12,100) tt,zb1(1)*th,zb2*th,zb3*th,zt3 end if CALL EXCH1(zb1,vb1,ab1,dzb1,dvb1,dab1,d1ab1,ND) CALL EXCH2(zb2,vb2,ab2,dzb2,dvb2,dab2,d1ab2)
三、激扰
x2 (t ) 0.05 sin(4t ) x3 (t ) 0.05 sin(6t )
implicit double precision (a-h,o-z) parameter(nn=40000,nd=2,gg=9.81,pi=3.14159)
parameter(th=1000.)
60
displacement /mm
0
-60
-120
0
2 Time/s
4
6
第四章 车辆系统动力学数值分析方法
位移（质量2）
k
2l m1 m3 Iy c z3 m2 z2 z1
60 30
Y Axis Title
0 -30 -60
0
2
4
6
三、外加激扰（不同频率）
x2 (t ) 0.05 sin(4t ) x3 (t ) 0.05 sin(8t )
由此，当 k2 / m2 时，主系统的振幅为 B1=0，即主系统不再振动， 因此使建筑物 稳态振动振幅最小吸振系统应满足的设计 条件为： 2 k2 m2 ，即使得吸振器的固 有频率等于主系统的固有频率。
车辆系统的振动微分方程可以化为如下形式： } [C ]{X } [ K ]{X } {P} [M ]{X (1)
（轮轨间具有Hertz非线性接触、非线性悬挂系 统、轮轨接触几何关系非线性）
• 目前，求解这类问题，只能采用直接数值积 分法
s s0 v0 t 1 / 2a t 2 v v0 a t
sn1 sn vn t 1 / 2an t 2 vn 1 vn an t sn 1 sn vn t (1 / 2 )an t 2 an1 t 2 vn1 vn (1 )an t - an-1 t
时间积分步长不便于选取、稳定性和精度一般较差；
计算过程简单、速度效率高。
一、模型
2l m1 k m3 z3 m2 Iy c z2 z1
二、参数
m1 3250kg m 2100kg 2 2 I y 3020kg m k 0.3M N/m c 20kN s/m l 1.25m
0.5
1.0 Time/s
1.5
2.0
2l
0
M1位移
k m3
m1
Iy c z3 m2 z2
z1
displacement /mm
-10
-20
-30 0.0
0.5
1.0 Time/s
1.5
2.0
2l
M2位移
0
k m3
m1
Iy c z3 m2 z2
z1
displacement /mm
-2
-4
-6 0.0
0.5
1.0 Time/s
1.5
2.0
二、外加激扰（相同频率）
x2 (t ) 0.05 sin(4t ) x3 (t ) 0.05 sin(4t )
0.06 0.03
Y Axis Title
0.00 -0.03 -0.06
0
1
2
3
X Axis Title
4
5
6
激扰（1）
2l
加速度（M1）
由此可以得到稳态响应的振幅为
k1 k 2 m 1 2 B1 k 2 k1a m 2e k 2 m 2 2 ( ) 2 B k k m k2 0 2 2 2 2
1
2 () (k1 k2 m1 2 )(k2 m2 2 ) k2
{X } { X }(或{ V })
[M ] [C ]
[K ]
系统的广义位移矢量 系统的广义速度矢量 系统的广义加速度矢量 质量矩阵 阻尼矩阵 刚度矩阵
X { }(或{ A })
} [C ]{X } [ K ]{X } {P} [M ]{X
(1)
{P}—系统的广义载荷矢量，一般难于显式表示， 往往是与广义位移和广义速度有非线性过程量。 • 这是一个非线性动力学微分方程。
大型非线性动力学数值求解现状
Newmark β Houbolt 隐式法 Wilson θ Hiber Hghes 的α法和 法 积分法 Park 显式法四阶Runge Kutta法 中心差分法
隐式法特点： 无论对线性还是非线性问题，其数值稳定性一般较好； 时间积分步长便于选取（积分步长可以取得较大）； 计算速度较慢（尤其是对于大型非线性问题）。 显式法特点： 经典的线性多步法中，不存在无条件稳定的显式法；
0.06 0.04 0.02 0.00 -0.02 -0.04 -0.06 0 1 2 3 4 5 6
激扰
0.06 0.04
Y Axis Title
0.02 0.00 -0.02 -0.04 -0.06 0 1 2 3 4 5 6
激扰（8）
60 30
Y Axis Title
0 -30 -60
0
m1 0 1 c c x 1 k1 k 2 0 x c c x k m2 x 2 2 2 k 2 x1 m 2 e sin(t ) x k2 0 2
CV=DV+(FI1*DA-FI*D1A)*VST
RETURN END
SUBROUTINE EXCH1(ZZC,VZC,AZC,DZZC,DVZC,DAZC,D1AZC,ND)
IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION ZZC(ND),VZC(ND),AZC(ND) DIMENSION DZZC(ND),DVZC(ND),DAZC(ND),D1AZC(ND) DO 5200 I=1,ND,1 DZZC(I)=ZZC(I) DVZC(I)=VZC(I) D1AZC(I)=DAZC(I) DAZC(I)=AZC(I) ZZC(I)=0. VZC(I)=0. AZC(I)=0. 5200 CONTINUE RETURN
CALL EXCH2(zb3,vb3,ab3,dzb3,dvb3,dab3,d1ab3)
end do 100 format(1x,f12.5,2x,f12.5,2x,f12.5,2x,f12.5,2x,f12.5) stop end
SUBROUTINE EXP1(DY,DV,DA,D1A,CY,CV,ND) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /EXPL/ VST,F1KC,FKC,FI1,FI
SUBROUTINE EXP2(DY,DV,DA,D1A,CY,CV) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /EXPL/ VST,F1KC,FKC,FI1,FI VST2=VST*VST CY=DY+DV*VST+(F1KC*DA-FKC*D1A)*VST2
一、无外加激扰（自由衰减振动）
M1加速度
10 5 0 -5 -10
k m3
2l m1 Iy c z3 m2 z2 z1
Acceleration m/s
2
0
1
2
Time/s
3
4
M1速度
m1
2l Iy c m3 z3 m2 z2 z1
8 4
k
Velocity m/s
0 -4 -8 0.0
DIMENSION DY(ND),DV(ND),DA(ND),D1A(ND),CY(ND),CV(ND)
VST2=VST*VST DO 700 I=1,ND,1
CY(I)=DY(I)+DV(I)*VST+(F1KC*DA(I)-FKC*D1A(I))*VST2
CV(I)=DV(I)+(FI1*DA(I)-FI*D1A(I))*VST 700 CONTINUE RETURN END
END
SUBROUTINE EXCH2(ZZC,VZC,AZC,DZZC,DVZC,DAZC,D1AZC) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DZZC=ZZC DVZC=VZC D1AZC=DAZC DAZC=AZC
ZZC=0.
VZC=0. AZC=0.
RETURN
车辆系统动力学计算机仿真
假设Mb为车体质量，Mc 为附加质量块，列车运行 速度为V，如何设置相关 参数使得车体振动完全静 止？
本题实际上是一个反共振问题，动力吸振器包括质量
m2 、弹簧 k 2 ， m1 与 k1 是主系统的质量及弹簧刚度，偏
心质量引起的机械垂向振动力为 m 2 e sin(t ) 。 由此 不难得到系统的强迫振动方程为：
zt2=0.05*sin(4.*pi*tt)
zt3=0.05*sin(4.*pi*tt) CALL EXP1(dzb1,dvb1,dab1,d1ab1,zb1,vb1,nd) ! 积分
CALL EXP2(dzb2,dvb2,dab2,d1ab2,zb2,vb2)
CALL EXP2(dzb3,dvb3,dab3,d1ab3,zb3,vb3)
COMMON /EXPL/ step,F1KC,FKC,FI1,FI dimension zb1(nd),vb1(nd),ab1(nd)
dimension dzb1(nd)源自文库dvb1(nd),dab1(nd),d1ab1(nd)
open(11,file='rusta.dat',status='unknown') open(12,file='rustd.dat',status='unknown')
cmk=1.0e+6;dmp=2.0e+4
cbk=1.0e+7;dbp=2.0e+4
bm1=3250.0 bm2=2100.0 bim=3020.0 bl=2.5
F1KC=1.0 FKC=0.5 fi1=1.5 fi=0.5 step=1.0e-4
do i=1,nn,1 tt=dble(i)*step
k
m1 m3
Iy c z3 m2 z2
z1
20 10
Y Axis Title
0 -10 -20
0
1
2
3
4
5
6
加速度（M2）
k
2l m1 m3 Iy c z3 m2 z2 z1
40 20
Y Axis Title
0 -20 -40
0
2
4
6
M1位移
k
2l m1 m3 Iy c z3 m2 z2 z1
! 积分
! 积分
! 求弹簧悬挂点质量 1的位移与速度
dm2=zb1(1)+0.5*bl*zb1(2) vm2=vb1(1)+0.5*bl*vb1(2) dm3=zb1(1)-0.5*bl*zb1(2) vm3=vb1(1)-0.5*bl*vb1(2)
fm2=cmk*(dm2-zb2)+dmp*(vm2-vb2) fm3=cmk*(dm3-zb3)+dmp*(vm3-vb3)
1
2
3
4
5
6
加速度（质量1）
40 20
Y Axis Title
0 -20 -40
0
1
2
3
4
5
6
加速度（质量2）
80 40
Y Axis Title
0 -40 -80
0
1
2
3
4
5
6
加速度（质量3）
100 50
Y Axis Title
0 -50 -100
0
1
2
3
4
5
6
位移（质量1）
60 30
Y Axis Title
fb2=cbk*(zb2-zt2)+dbp*vb2
fb3=cbk*(zb3-zt3)+dbp*vb3
ab1(1)=(-fm2-fm3)/bm1-gg
ab1(2)=0.5*bl*(fm3-fm2)/bim
ab2=(-fb2+fm2)/bm2-gg
ab3=(-fb3+fm3)/bm2-gg
if(mod(nn,800).eq.0)then write(11,100) tt,ab1(1),ab2,ab3,zt2 write(12,100) tt,zb1(1)*th,zb2*th,zb3*th,zt3 end if CALL EXCH1(zb1,vb1,ab1,dzb1,dvb1,dab1,d1ab1,ND) CALL EXCH2(zb2,vb2,ab2,dzb2,dvb2,dab2,d1ab2)
三、激扰
x2 (t ) 0.05 sin(4t ) x3 (t ) 0.05 sin(6t )
implicit double precision (a-h,o-z) parameter(nn=40000,nd=2,gg=9.81,pi=3.14159)
parameter(th=1000.)
60
displacement /mm
0
-60
-120
0
2 Time/s
4
6
第四章 车辆系统动力学数值分析方法
位移（质量2）
k
2l m1 m3 Iy c z3 m2 z2 z1
60 30
Y Axis Title
0 -30 -60
0
2
4
6
三、外加激扰（不同频率）
x2 (t ) 0.05 sin(4t ) x3 (t ) 0.05 sin(8t )
由此，当 k2 / m2 时，主系统的振幅为 B1=0，即主系统不再振动， 因此使建筑物 稳态振动振幅最小吸振系统应满足的设计 条件为： 2 k2 m2 ，即使得吸振器的固 有频率等于主系统的固有频率。
车辆系统的振动微分方程可以化为如下形式： } [C ]{X } [ K ]{X } {P} [M ]{X (1)
（轮轨间具有Hertz非线性接触、非线性悬挂系 统、轮轨接触几何关系非线性）
• 目前，求解这类问题，只能采用直接数值积 分法
s s0 v0 t 1 / 2a t 2 v v0 a t
sn1 sn vn t 1 / 2an t 2 vn 1 vn an t sn 1 sn vn t (1 / 2 )an t 2 an1 t 2 vn1 vn (1 )an t - an-1 t
时间积分步长不便于选取、稳定性和精度一般较差；
计算过程简单、速度效率高。
一、模型
2l m1 k m3 z3 m2 Iy c z2 z1
二、参数
m1 3250kg m 2100kg 2 2 I y 3020kg m k 0.3M N/m c 20kN s/m l 1.25m
0.5
1.0 Time/s
1.5
2.0
2l
0
M1位移
k m3
m1
Iy c z3 m2 z2
z1
displacement /mm
-10
-20
-30 0.0
0.5
1.0 Time/s
1.5
2.0
2l
M2位移
0
k m3
m1
Iy c z3 m2 z2
z1
displacement /mm
-2
-4
-6 0.0
0.5
1.0 Time/s
1.5
2.0
二、外加激扰（相同频率）
x2 (t ) 0.05 sin(4t ) x3 (t ) 0.05 sin(4t )
0.06 0.03
Y Axis Title
0.00 -0.03 -0.06
0
1
2
3
X Axis Title
4
5
6
激扰（1）
2l
加速度（M1）
由此可以得到稳态响应的振幅为
k1 k 2 m 1 2 B1 k 2 k1a m 2e k 2 m 2 2 ( ) 2 B k k m k2 0 2 2 2 2
1
2 () (k1 k2 m1 2 )(k2 m2 2 ) k2
{X } { X }(或{ V })
[M ] [C ]
[K ]
系统的广义位移矢量 系统的广义速度矢量 系统的广义加速度矢量 质量矩阵 阻尼矩阵 刚度矩阵
X { }(或{ A })
} [C ]{X } [ K ]{X } {P} [M ]{X
(1)
{P}—系统的广义载荷矢量，一般难于显式表示， 往往是与广义位移和广义速度有非线性过程量。 • 这是一个非线性动力学微分方程。
大型非线性动力学数值求解现状
Newmark β Houbolt 隐式法 Wilson θ Hiber Hghes 的α法和 法 积分法 Park 显式法四阶Runge Kutta法 中心差分法
隐式法特点： 无论对线性还是非线性问题，其数值稳定性一般较好； 时间积分步长便于选取（积分步长可以取得较大）； 计算速度较慢（尤其是对于大型非线性问题）。 显式法特点： 经典的线性多步法中，不存在无条件稳定的显式法；
0.06 0.04 0.02 0.00 -0.02 -0.04 -0.06 0 1 2 3 4 5 6
激扰
0.06 0.04
Y Axis Title
0.02 0.00 -0.02 -0.04 -0.06 0 1 2 3 4 5 6
激扰（8）
60 30
Y Axis Title
0 -30 -60
0
m1 0 1 c c x 1 k1 k 2 0 x c c x k m2 x 2 2 2 k 2 x1 m 2 e sin(t ) x k2 0 2
CV=DV+(FI1*DA-FI*D1A)*VST
RETURN END
SUBROUTINE EXCH1(ZZC,VZC,AZC,DZZC,DVZC,DAZC,D1AZC,ND)
IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION ZZC(ND),VZC(ND),AZC(ND) DIMENSION DZZC(ND),DVZC(ND),DAZC(ND),D1AZC(ND) DO 5200 I=1,ND,1 DZZC(I)=ZZC(I) DVZC(I)=VZC(I) D1AZC(I)=DAZC(I) DAZC(I)=AZC(I) ZZC(I)=0. VZC(I)=0. AZC(I)=0. 5200 CONTINUE RETURN
CALL EXCH2(zb3,vb3,ab3,dzb3,dvb3,dab3,d1ab3)
end do 100 format(1x,f12.5,2x,f12.5,2x,f12.5,2x,f12.5,2x,f12.5) stop end
SUBROUTINE EXP1(DY,DV,DA,D1A,CY,CV,ND) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /EXPL/ VST,F1KC,FKC,FI1,FI
SUBROUTINE EXP2(DY,DV,DA,D1A,CY,CV) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /EXPL/ VST,F1KC,FKC,FI1,FI VST2=VST*VST CY=DY+DV*VST+(F1KC*DA-FKC*D1A)*VST2