用matlab计算运动方程的问题

shu_yl

我用matlab编程求解移动荷载作用下简支梁桥的振动分析，用newmark数值分析方法对运动微分方程求解，得出来的梁中点位移的数量级是十的－6次方，而手算大概这个数量级应该在－3。我用wilson发重新算了下，还是－6次方，我就郁闷了，哪位大哥帮分析下，谢啦！

fanghuikeer

程序传上来,这样怎么分析呢??

shu_yl

参数belt,gama,time,delt,n,v，我分别取0.25，0.5，6，0.01，5，5
function [s,w,w1,Fm,Mm,Km,Xm,Vm,Am,xp]=newmark2(belt,gama,time,delt,n,v)
N=time/delt;
i=1; L=30;%梁的长度
Xm=zeros(n,N);%设定存储位移矩阵
Vm=zeros(n,N);%设定存储速度矩阵
Am=zeros(n,N);%设定存储加速度矩阵
for t=delt:delt:time
    xp=v*t;%移动荷载的速度
    [Fm,Mm,Km]=sovle_matrix2(xp,n);
    KKm=Km+Mm/(belt*delt^2);
    if i==1
        Xm(:,i)=inv(KKm)*Fm;
        Vm(:,i)=Xm(:,i)*gama/(belt*delt);
        Am(:,i)=Xm(:,i)/(belt*delt^2);
    else
        FFm=Fm+(Mm/(belt*delt^2))*Xm(:,i-1)+(Mm/(belt*delt))*Vm(:,i-1)+(Mm*(1/(2*belt)-1))*Am(:,i-1);
        Xm(:,i)=inv(KKm)*FFm;
        Vm(:,i)=(Xm(:,i)-Xm(:,i-1))*gama/(belt*delt)-Vm(:,i-1)*(gama/belt-1)-Am(:,i-1)*delt*(gama/belt-2)/2;
        Am(:,i)=(Xm(:,i)-Xm(:,i-1))/(belt*delt^2)-Vm(:,i-1)/(belt*delt)-Am(:,i-1)*(1/(2*belt)-1);
    end
    i=i+1;
end
t=1:N;t=t*delt;w=zeros(1,N);w1=zeros(n,N);
for i=1:n
    w1(i,:)=Xm(i,:)*sin(i*pi*15/L);
    w=w+w1(i,:);
end
s=v*t;
%画出梁中点位移的曲线，在荷载移动过程中
plot(s,w)


*********************************************************************************

function [Fm,Mm,Km]=sovle_matrix2(xp,n)
Fm=zeros(n,1);%振型力向量
Mm=zeros(n);%质量矩阵
Km=zeros(n);%刚度矩阵
L=30;mb=3.6e+4;EI=2.943*8.72e+10;M=14000;g=9.8;
%形成质量和刚度矩阵
for i=1:n
    a(i)=L*mb/2;
    b(i)=i^2*pi^2*EI/(2*L);
end
m=diag(a);
k=diag(b);
for i=1:n
    for j=1:n
        Km(i,j)=k(i,j);
        Mm(i,j)=m(i,j);
    end
    Fm(i)=M*g*sin(i*pi*xp/L);
end

[ 本帖最后由 shu_yl 于 2006-10-25 16:27 编辑 ]

shu_yl

顶一哈！

happy

可能是单位之类的有问题，先把方程做一下无量纲化看看

shalingfeng

顶一哈！