tianlei2005 发表于 2007-7-19 21:49

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

下午看了看书,在网上找了找代码,试着用newmark法求解动力学方程,但是结果不对。好像以前在网上有人遇到过相似的问题,
请教达人指点!
%用newmark法计算一个3自由度结构的响应:
%参数初始化:
N=3;%N为自由度数。
M=;M=diag(M);
K=;
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)=;
    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 发表于 2007-7-20 08:41

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

tianlei2005 发表于 2007-7-20 11:01

回复 #2 无水1324 的帖子

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

appleseed05 发表于 2007-7-20 12:48

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

pp786195 发表于 2009-11-9 14:49

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

lin20080111 发表于 2009-11-9 16:36

本帖最后由 VibInfo 于 2016-5-17 13:22 编辑

原帖由 pp786195 于 2009-11-9 14:49 发表
弱弱的问一句,对于每步计算出来的结果需要进行收敛判断吗?
只要满足newmark方法时间步长大小的要求,每步计算出来的结果不需要进行收敛判断

hew136 发表于 2009-11-9 17:23

回复 楼主 tianlei2005 的帖子

同意无水给一个合适的C值和合适的时间不长(要满足NEWMARK方法的时间步长要求)应该收敛的
页: [1]
查看完整版本: 请教newmark求解动力学方程遇到的问题