|
楼主 |
发表于 2011-5-26 17:50
|
显示全部楼层
clear all;
clc
global M rou A L E I k v c1 c2 delta1 g alpha1 miu1 alpha2 beta2 f delta2 alpha3 xi2 xi3 gama3 l k1 m
rou=7.2*10^3;
A=0.36;
E=1.995*10^8;
I=0.6;
M=2.5*10^3;
L=25;
v=20;
g=9.8;
c1=3.14*10^3;
c2=2*10^3;
k=3.4*10^4;
k1=3.14*10^2;
m=1;
omega=(pi*v)/L;
alpha2=2*k/(rou*A*L*omega^2);
beta2=2*c2/(rou*A*L*omega^2);
delta1=E*I/(rou*A*omega^2)*(pi/L)^4-alpha2/2;
alpha1=pi^4*E/(4*rou*A*omega^2*L^4);
miu1=c1/(rou*A*omega^2)-beta2/2;
f =2*M/(rou*A*L*omega^2)*g;
delta2= k/M/omega^2 ;
xi2=c2/(rou*A*omega^2);
alpha3=2*sqrt(2)*k1/(rou*A*L*omega^2);
xi3=0.1;
gama3=2*k1/m-sqrt(2)/2*alpha3;
l=8;
period=2*pi/omega;
tspan=0:period/100:50*period
u0=[0,0,0,0,0,0];
% tspan=0:0.01:1
figure(1)
[t,u]=ode45('CQ6',tspan,u0)
plot(t,u(:,1))
hold on
plot(t,u(:,3),'r-')
legend('vehicle','bridge'
子函数:function du=CQ6(t,U)
global delta1 delta2 alpha1 miu1 alpha2 xi2 beta2 omega f gama3 alpha xi3 alpha3 l g
omega =2.5133;
alpha2=0.1661;
beta2 =0.0098;
delta1 =1.7401;
alpha1 =0.7596;
miu1 =0.1869;
f =0.1197;
delta2 =2.1531;
xi2 = 0.1222;
alpha3 =0.0022;
xi3 =0.1000;
gama3 =627.9985;
l =8;
alpha=0.1;
g=9.8;
u=U(1);
x=U(2);
y=U(3);
z=U(4);
p=U(5);
q=U(6);
du=zeros(6,1);
du=[x;
-delta1*u-alpha1*u^3-miu1*x-alpha2*y*sin(omega*t)-alpha2/2*u*cos(2*omega*t)-beta2*z*sin(omega*t)-1/2*beta2*x*cos(2*omega*t)+f*sin(omega*t)+alpha3*p*(1-l/sqrt(alpha^2+p^2));
z;
-delta2*y-xi2*z+delta2*u*sin(omega*t)+xi2*x*sin(omega*t);
q;
gama3*p*(1-l/sqrt(alpha^2+p^2))-xi3*q-g+sqrt(2)/2*(delta1*u+alpha1*u^3+miu1*x+alpha2*y*sin(omega*t)+alpha2/2*u*cos(2*omega*t)+beta2*z*sin(omega*t)+1/2*beta2*x*cos(2*omega*t)-f*sin(omega*t))];
改成这个之后这个是可以运行的。但是就是想不通这和矩阵维数投什么关联??
|
|