MATLAB之newmark法计算的结构响应与ANSYS的结果差异为何如此之大?
请问高手:我用MATLAB编了一个newmark法计算结构响应的小程序,又用ANSYS将相同的模型计算了一下,结果两者的结果差异很大(数量级差10000倍),不知道为什么?附MARLAB源程序及ANSYS命令流:
==========================
MATLAB源程序:
clc;clear;close;
M =1.0e+003 * [5.2000 0 0 1.3000 0 0
0 5.7943 0 0 1.0029 -0.2414
0 0 0.1486 0 0.2414 -0.0557
1.3000 0 0 5.2000 0 0
0 1.0029 0.2414 0 5.7943 0
0 -0.2414 -0.0557 0 0 0.1486];
C =1.0e+006 * [6.8267 0 0 -3.3955 0 0
0 6.8267 0 0 -3.3955 1.6989
0 0 2.2674 0 -1.6989 0.5664
-3.3955 0 0 6.8267 0 0
0 -3.3955 -1.6989 0 6.8267 0
0 1.6989 0.5664 0 0 2.2674];
K =1.0e+010 * [4.0017 0 0 -2.0008 0 0
0 4.0001 0 0 -2.0000 1.0000
0 0 1.3334 0 -1.0000 0.3333
-2.0008 0 0 4.0017 0 0
0 -2.0000 -1.0000 0 4.0001 0
0 1.0000 0.3333 0 0 1.3334];
u=';
v=';
a=';
t(1)=0; %时间
x(:,1)=u; %位移
x1(:,1)=v; %速度
x2(:,1)=a; %加速度
%newmark参数
gama=0.5;
dt=2e-1;
delta=0.25;
b1=1/(delta*dt^2);
b2=1/(delta*dt);
b3=1/(2*delta-1);
b4=gama/(delta*dt);
b5=gama/delta-1;
b6=dt*gama/(2*delta)-dt;
b7=1/(2*delta)-1;
%等效刚度矩阵
K1=K+b1*M+b4*C;
t_max=12; %计算时间总长
i=1;
t(1)=0.1;
q=zeros(6,1);
while t(i)<t_max
f=1*sin(15*t(i)*0.2);
%f=sin(t(i));
Q(:,i)=';
q(:,i+1)=Q(:,i)+M*(b1*x(:,i)+b2*x1(:,i)+b3*x2(:,i))+C*(b4*x(:,i)+b5*x1(:,i)+b6*x2(:,i));
x(:,i+1)=inv(K1)*q(:,i+1);
x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i);
x2(:,i+1)=b1*(x(:,i+1)-x(:,i))-b2*x1(:,i)-b7*x2(:,i);
i=i+1;
t(i)=t(i-1)+dt;
end
format long e
t=10:50;
%xrespons=
%x1respons=
%x2respons=
%plot(t,x(:,t),t,x1(:,t),t,x2(:,t)/10000000)
plot(t,x(:,t))
=================================
ANSYS命令流:
!施加正弦荷载进行瞬态分析
!用beam3计算
finish
/clear
/title,the transient analysis of a beam use of beam3
/FILNAME,frame3,1
/CWD,'C:\ansyswork\frame2'
save
/prep7
et,1,beam3
R,1,1,1/12,1, , , ,
mp,ex,1,2e11
mp,dens,1,7800
mp,prxy,1,0.2
k,1
k,2,0,1
k,3,1,1
k,4,1
lstr,1,2
lstr,2,3
lstr,3,4
lesize,all,,,1
lmesh,all
/SOLU
DK,1, , , ,0,ALL, , , , , ,
DK,4, , , ,0,ALL, , , , , ,
/solu
antype,4
trnopt,full
!dk,1,all
!dk,4,all
time,0
autots,0
deltim,0.2, , ,1
kbc,0
alphad,1.0926529, !由计算得知
betad,0.0001737,
outpr,basic,all,
*do,i,1,60,1 !计算到6秒
time,i*0.2
autots,0
deltim,0.2, , ,1
kbc,0
fn=sin(15*i*0.2)
fk,2,fx,fn
solve
*enddo
===========================
因为这个模型极其简单,所以我实在想不出理由。个人认为可能是MATLAB中采用的M、C、K与ANSYS模型不符,那应该是怎样的呢?肯请高人指点!!先行谢过!!
[ 本帖最后由 realyyy 于 2006-8-18 15:51 编辑 ] 个人认为是没有根据的猜测。你那么多的系数,为什么不好好查一查是否出错?比如b4=gama*(delta*dt)就不对,应该是b4=gama/(delta*dt). 首先感谢tonnyw!
系数已核对并改正!结果没有太大变化啊。。。。。。郁闷中。。。。。。 x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i);
x2(:,i+1)=b1*(x(:,i+1)-x(:,i))-b2*x1(:,i)-b7*x2(:,i);
速度公式不对,记住先用位移求加速度,用得到的加速度求速度。 把下面两句的顺序换一下对结果没有影响啊:
x2(:,i+1)=b1*(x(:,i+1)-x(:,i))-b2*x1(:,i)-b7*x2(:,i);
x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i);
........................
期盼高人啊! 请仔细检查一下,我说了你的速度公式不对,顺序不对。你现在调整了顺序,但公式还是不对。 感谢tonnyw的热心关注!
我检查了啊:用得到的加速度求速度的式子整理后和x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i)是一样的。 x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i)你这里系数前面都是负号,很明显不对。仔细对照一下我贴的公式,就知道错在什么地方了。 再次感谢tonnyw!
我又认真、仔细的核对了一下所有的系数,发现了问题:我是按照《有限单元法》(薛守义编著)中的公式写的程序,系数与其书中所述完全一致,但与你所贴的就不一致。不知你用的是哪本书中的公式?而且更惊奇的是:我发现《机械振动基础》(胡海岩编著)一书中所述的公式与你我所用公式也不一样!!
到底大家对NEWMARK公式是怎么搞的吗?!!只能有一个是对的。。。。。。其他的全错!
更加郁闷中。。。。。。 最后经本人验证:tonnyw的贴图中的公式完全正确!程序运行后结果相当准确。
(还是外国的公式准!没办法。)
本讨论至此完美结束。再次感谢tonnyw的热心帮助!
求助随机有限元的程序
求助随机有限元的程序本人一点头绪都没有
应该如何开始做啊 国内的大专家们写的书,
sign,,, tonnyw 贴的是大牛K.J.Bathe的书,呵呵 不是吧!
国内的人写书来骗人啦
真是郁闷啦,把我们这些人给害掺了不是! 问一下,MATLAB计算的每一步结果,需要判断收敛吗,怎么判断? 静力分析不需要判断,如果存在材料几何的非线性的话要进行矩阵的更新和迭代分析 15# pp786195 jing 厉害啊 学习了 学习了。。 都是牛人啊!学习!!!! 请问质量矩阵和刚度矩阵,阻尼矩阵是怎么得到的呢?采用模态分析还是数学建模计算?
页:
[1]