请教newmark求解动力学方程遇到的问题

tianlei2005

下午看了看书，在网上找了找代码，试着用newmark法求解动力学方程，但是结果不对。好像以前在网上有人遇到过相似的问题，
请教达人指点！
%用newmark法计算一个3自由度结构的响应：
%参数初始化：
N=3;%N为自由度数。
M=[1 1 1];M=diag(M);
K=[3 -1 0;-1 2 -1;0 -1 3];
C=0.*(M+K);
d=10;%间隔时间为1/d秒。
totaltime=20;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就是γ和β。
dt=1/d;
%按公式计算积分常数bi：
alpha0=1/beta/dt^2;
alpha1=gama/beta/dt;
alpha2=1/beta/dt;
alpha3=1/2/beta-1;
alpha4=gama/beta-1;
alpha5=dt/2*(gama/beta-2);
alpha6=dt*(1-gama);
alpha7=gama*dt;
EK=K+alpha0*M+alpha1*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(3*t);0;0];
    t1=t1+1;
end
t1=1;
for t=0:1/d:totaltime
    %计算t+detat时刻的有效载荷：
    Ef(:,t1+1)=f(:,t1+1)+M*(alpha0*u(:,t1)+alpha2*du(:,t1)+alpha3*ddu(:,t1))+C*(alpha1*u(:,t1)+alpha4*du(:,t1)+alpha5*ddu(:,t1));
    %计算t+detat时刻的位移：
    u(:,t1+1)=inv(EK)*Ef(:,t1+1);
    %计算t+detat时刻的速度和加速度：
    ddu(:,t1+1)=alpha0*(u(:,t1+1)-u(:,t1))-alpha2*du(:,t1)-alpha3*ddu(:,t1);
    du(:,t1+1)=u(:,t1)+alpha6*du(:,t1)+alpha7*ddu(:,t1+1);
    t1=t1+1;
end
t=0:0.1:20+0.1;
% plot(t,Ef(1,1:tt))；
plot(t,u(1,:))；
grid on
clear

[ 本帖最后由 无水1324 于 2007-7-20 08:41 编辑 ]

无水1324

C=0.*(M+K);
这个阻尼有什么用？

tianlei2005

以前看到交大有个老师的讲义，例子是c＝0，给出了解析解和Rk方法的结果。所以我这里也让c＝0，验证我的代码是否正确

appleseed05

对于每一个时刻你这里面应该要做迭代阿

pp786195

弱弱的问一句，对于每步计算出来的结果需要进行收敛判断吗？

lin20080111

原帖由 pp786195 于 2009-11-9 14:49 发表
弱弱的问一句，对于每步计算出来的结果需要进行收敛判断吗？

只要满足newmark方法时间步长大小的要求，每步计算出来的结果不需要进行收敛判断

hew136

同意无水给一个合适的C值和合适的时间不长（要满足NEWMARK方法的时间步长要求）应该收敛的