高铁车辆和轨道耦合模型

车辆模型参数

g=9.8;                     %重力加速度(m/s^2)
Lc=9.0;                    %车辆定距之半(m)
Lt=1.25;                   %转向架固定轴距之半(m)
Kpz=1.87E6;                %一系悬挂刚度(N/m)
Cpz=5E5;                   %一系悬挂阻尼(N.s/m)
Ksz=1.72E6;                %二系悬挂刚度(N/m)
Csz=1.96E5;                %二系悬挂阻尼(N.s/m)
Mc=30000;                  %车体质量(Kg)
Mt=3000;                   %前后转向架构架质量(Kg)
Mw=1000;                   %各轮对质量(Kg)
Jc=2.0E6;                 %车体点头惯量(Kg.m^2)
Jt=2000;                   %前后转向架点头惯量(Kg.m^2)
R=0.4575;                  %车轮半径(m)

高速线路轨道模型参数

E=2.059E11;                %钢轨弹性模量(Pa/m)
I=3.217E-5;                %钢轨截面惯性矩(m^4)
Mr=60.64;                  %钢轨每延米质量(Kg)
Ms=170;                    %每根轨枕质量(Kg)
Kp=6.0E7;                  %轨下垫层刚度(N/m)
Cp=7.5E4;                  %轨下垫层阻尼(N.s/m)
Ls=0.6;                    %轨枕间距(m)
Le=1.175;                  %半轨枕有效支承长度
Lb=0.277;                  %轨枕底面宽度(m)
Rho=1.9E3;                 %道床密度(Kg/m^3)
Eb=1.2E8;                  %道床弹性模量(Pa/m)
Cb=6.0E4;                  %道床阻尼(N.s/m)
Kw=7.84E7;                 %道床剪切刚度(N/m)       
Cw=8.0E4;                  %道床剪切阻尼(N.s/m)
a=35;                      %轨下压力扩散角(°)
hb=0.35;                   %道床厚度(m)
Ef=1.9E8;                  %路基弹性模量(Pa/m)
Cf=1.0E5;                  %路基阻尼(N.s/m)

计算耦合系统中的各矩阵

M=blkdiag(Mv,MT);
C=blkdiag(Cv,CT);
K=blkdiag(Kv,KT);

计算各车轮的运动坐标

xw(:,t+1)=[2*(Lc+Lt)+v*t.*dt;2*Lc+v*t.*dt;2*Lt+v*t.*dt;v*t.*dt];

for j=1:4
    Zrw(j,t+1)=0;
    Vrw(j,t+1)=0;
    Arw(j,t+1)=0;
    for k=1:NM
  Zrw(j,t+1)=Zrw(j,t+1)+(2/(Mr*L))^0.5*sin(k*pi*xw(j,t+1)).*X(10+k,t+1);
  Vrw(j,t+1)=Vrw(j,t+1)+(2/(Mr*L))^0.5*sin(k*pi*xw(j,t+1)).*V(10+k,t+1);
  Arw(j,t+1)=Arw(j,t+1)+(2/(Mr*L))^0.5*sin(k*pi*xw(j,t+1)).*A(10+k,t);
    end
    Z0(:,t+1)=0.005*sind(20.*t/1000)+0.003*cosd(10.*t/1000+30);
    Deta_Z(j,t+1)=X(j+6,t+1)-Zrw(j,t+1)-Z0(:,t+1);
    pw(j,t+1)=(1/G.*Deta_Z(j,t+1))^1.5;
    if  pw(j,t+1)<0
        pw(j,t+1)=0;
    else pw(j,t+1)=pw(j,t+1);
    end
    sum_pw=sum_pw+pw(j,2);
    sum_deta_pw=sum_deta_pw+abs(pw(j,2)-pw(j,1));
end
if sum_deta_pw./sum_pw>1e-2
    pw(1:4,2)=1e-2.*sum_pw;
end