• 春华秋实*
    了解作者
  • matlab
    开发工具
  • 1KB
    文件大小
  • zip
    文件格式
  • 5
    收藏次数
  • 10 积分
    下载积分
  • 78
    下载次数
  • 2018-04-01 11:39
    上传日期
运用龙格库塔法-ode45求解一个多自由度的振动方程,得到振动的时域响应曲线
龙格库塔法解振动方程.zip
  • 龙哥库塔法解振动方程
  • xiao.m
    2.5KB
  • xiao2.m
    281B
内容介绍
%轨道垂向振动模型 %四轴车辆 %变步长龙格库塔法(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
评论
    相关推荐