龙格库塔法解振动方程.zip

%轨道垂向振动模型 %四轴车辆 %变步长龙格库塔法（ode45） clear clear global clc close all %% 车辆参数 global A B F vc lc lt Ir ktz mc=41800;%车体质量 mt=7480;%构架质量 mw=1800;%轮对质量(拖车) Icy=2.446e6;%车体点头惯量（绕y轴） Ity=3605;%构架点头惯量（绕y轴） kpz=2.14e6;%一系悬挂刚度（垂向） ksz=2.535e6;%二系悬挂刚度（垂向） ktz=1e9;% 轮轨接触刚度 cpz=4.9e4;%一系悬挂阻尼（垂向） csz=1.95e5;%二系悬挂阻尼（垂向） lc=8.4;%车辆定距之半 lt=1.25;%车辆轴距之半 vc=60/3.6;%车辆运行速度 g=9.8;%重力加速度 Vf=6;%车辆自由度 %% 钢轨不平顺样本 load Irr1 Ir=Irr(1:2,:)*1.5; %% 车辆参数矩阵的生成 %惯量矩阵 M=[ mc , 0 , 0 , 0 , 0 ,0; 0 , Icy , 0 , 0 , 0, 0; 0 , 0 , mt, 0 , 0 ,0; 0 , 0 , 0 , Ity , 0 ,0; 0 , 0 , 0 , 0 , mw ,0; 0 , 0 , 0 , 0 , 0, mw]; %刚度矩阵 K=[ksz , 0 , -ksz , 0 , 0 , 0 ; 0 , ksz*lc^2 , ksz*lc , 0 , 0 , 0 ; -ksz , ksz*lc , 2*kpz+ksz , 0 , -kpz , -kpz ; 0 , 0 , 0 , 2*kpz*lt^2 , kpz*lt , -kpz*lt ; 0 , 0 , -kpz , kpz*lt , kpz+2*ktz, 0 ; 0 , 0 , -kpz , -kpz*lt , 0 ,kpz+2*ktz ]; %阻尼矩阵 C=[csz , 0 , -csz , 0 , 0 , 0 ; 0 , csz*lc^2 , csz*lc , 0 , 0 , 0 ; -csz , csz*lc , 2*cpz+csz , 0 , -cpz , -cpz ; 0 , 0 , 0 , 2*cpz*lt^2 , cpz*lt , -cpz*lt ; 0 , 0 , -cpz , cpz*lt , cpz , 0 ; 0 , 0 , -cpz , -cpz*lt , 0 , cpz ]; % 力向量 Q=[mc*g;0;mt*g;0;mw*g;mw*g]; %Q=[0;0;0;0;0;0;0;0;0;0]; %计算矩阵 O1=zeros(Vf,Vf);O2=zeros(Vf,1); A=[C,M;M,O1]; B=[K,O1;O1,-M]; F=[Q;O2]; %% 数值求解 z0=zeros(2*Vf,1);% 系统初值 fs=1e4;% 计算步长 tx=8;% 计算时间 tspan=0:1/fs:tx;% 时间序列 tic [t,z]=ode45(@xiao2,tspan,z0,[]); toc %% 数据分类存储 Zc=z(:,1:Vf); Vc=z(:,Vf+1:2*Vf); Ac=diff(z(:,Vf+1:2*Vf))*fs; %% 结果绘图 figure(1) for i=1:Vf subplot(1,6,i) plot(tspan,Zc(:,i),'b','linewidth',1.5);xlabel('时间[s]') legend('未加非线性刚度','加入非线性刚度'); end figure(2) for i=1:Vf subplot(1,6,i) plot(tspan,Vc(:,i),'b','linewidth',1.5);xlabel('时间[s]') end figure(3) for i=1:Vf subplot(1,6,i) plot(tspan(1,2:end),Ac(:,i),'b','linewidth',1.5);xlabel('时间[s]');ylabel('加速度[m/s^2]') end