[求助]NewMark法计算结构响应的问题

realyyy

我按照书上的NEWMARK原理编了一段MATLAB程序来计算一个3自由度结构的响应，可是感觉有问题，下面是代码，请高手指点！

%用newmark法计算一个3自由度结构的响应：
%参数初始化：
clear
clc
N=3;
%N为自由度数。
M=[1 1 1];M=diag(M);
C=[0.3];
K=9*[1 1 1];K=diag(K);
d=5;
%间隔时间为1/d秒。
totaltime=6;tt=totaltime*d+2;
%totaltime为外力总持续时间，tt为划分的时间点数，划分精细程度由d决定。
u=zeros(N,tt);
du=zeros(N,tt);
ddu=zeros(N,tt);
%u、du、ddu分别表示位移、速度和加速度。
gama=0.5;
beta=0.25*(0.5+gama^2);
%gama和beta就是γ和β。
detat=1/d;
%按公式计算积分常数bi：
b1=1/(beta*detat^2);
b2=1/(beta*detat);
b3=-1/(2*beta)+1;
b4=gama*b1*detat;
b5=1+gama*b2*detat;
b6=(1+gama*b3-gama)*detat;
EK=K+b1*M+b4*C;
%EK表示有效刚度矩阵。
f=zeros(N,tt);
Ef=zeros(N,tt);
t1=1;
%用t1是为了保证矩阵索引号为整数。
for t=0:1/d:totaltime
    %因为detat=0.1，所以t=0:0.1:totalt
    f(:,t1)=[sin(t);0;0];
    t1=t1+1;
end
t1=1;
for t=0:1/d:totaltime
    %计算t+detat时刻的有效载荷：
    Ef(:,t1+1)=f(:,t1+1)+M*(b1*u(:,t1)-b2*du(:,t1)-b3*ddu(:,t1))+C*(b4*u(:,t1)-b5*du(:,t1)-b6*ddu(:,t1));
    %计算t+detat时刻的位移：
    u(:,t1+1)=inv(EK)*Ef(:,t1+1);
    %计算t+detat时刻的速度和加速度：
    du(:,t1+1)=b4*(u(:,t1+1)-u(:,t1))+b5*du(:,t1)+b6*ddu(:,t1);
    ddu(:,t1+1)=b1*(u(:,t1+1)-u(:,t1))+b2*du(:,t1)+b3*ddu(:,t1);
    t1=t1+1;
end
t=1:tt;
plot(t,Ef(1,1:tt))
%为什么Ef的曲线不是正弦形式？（位移u的也不是）
grid on
clear

[ 本帖最后由 realyyy 于 2006-8-13 19:42 编辑 ]

happy

这种问题最难解决了，只能自己慢慢查

guodan02312

我的也是，我想是不是激励力激起了一部分自由振动，从而两个波叠加在了一起。我用振型叠加法算出来的是正弦

sogooda

建议楼主找找动力学书上的经典例子验证一下程序。

yangfengdehaizi

有没有非线性动力方程的newmark法 计算程序啊

meiyongyuandeze

Newmark方法在数值计算的书中都有现成的代码，可以参考！用代码计算下duffing方法来验证算法！但该算法中的几个参数的确定是有讲究的，希望注意下！