|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
%机车质量矩阵
Mv1=[mv 0;0 Iv];
Mv2=[m1 0;0 m2];
%机车阻尼、刚度矩阵元素
Cv_11=[Cs1+Cs2 (-Cs1*a1+Cs2*a2)*S;(-Cs1*a1+Cs2*a2)*S (Cs1*a1^2+Cs2*a2^2)*S^2];
Cv_12=[-Cs1 -Cs2;Cs1*a1*S -Cs2*a2*S];
Cv_21=[-Cs1 Cs1*a1*S;-Cs2 -Cs2*a2*S];
Cv_22=[Cs1 0;0 Cs2];
Kv_11=[Ks1+Ks2 (-Ks1*a1+Ks2*a2)*S;(-Ks1*a1+Ks2*a2)*S (Ks1*a1^2+Ks2*a2^2)*S^2];
Kv_12=[-Ks1 -Ks2;Ks1*a1*S -Ks2*a2*S];
Kv_21=[-Ks1 Ks1*a1*S;-Ks2 -Ks2*a2*S];
Kv_22=[Ks1 0;0 Ks2];
…
M=[Mb zeros(2*(j+1),2) Hc*Mv2;zeros(2,2*(j+1)) Mv1 zeros(2)];
C=[Cb Hc*Cv_21 Hc*Cv_22;zeros(2,2*(j+1)) Cv_11 Cv_12];
K=[Kb Hc*Kv_21 Hc*Kv_22;zeros(2,2*(j+1)) Kv_11 Kv_12];
g=9.8;
Ms=[(m1+a2*mv)*g;(m2+a1*mv)*g];%机车静载
F=[Hc*Ms;zeros(2,1)];
%简支梁边界条件
%theta1=0,Omega1=0;Omega(j+1)=0,theta'(j+1)=0
function xdot=B_v_system(time,x,M,C,K,F,j)
NN=2*(j+1);%自由度数
xdot(1:(NN+4),1)=x((NN+5):2*(NN+4),1);
xdot((NN+5):2*(NN+4),1)=M\(F-K*x(1:(NN+4),1)-C*x((NN+5):2*(NN+4),1));
%车-桥系统
%X=[Omega1;theta1;Omega2;theta2;……Omega(j+1);theta(j+1);yv;thetav;y1;y2];系统坐标
NN=2*j;
%仿真程序
t0=0;
h=0.01;
tfinal=5;
tspan=[t0 tfinal];
x0=zeros(2*NN+8,1);
[time,x]=ode45('B_v_system',tspan,x0,options);
H为三次样条插值形函数,Omega1;theta1…为梁单元节点转角与挠度,因为M,C,K并非方阵,因此方程组中未知数多与方程数(少两个方程),是否应该加入简支梁边界条件,怎样加??谢谢
|
|