realyyy 发表于 2006-8-17 13:59:29

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 编辑 ]

tonnyw 发表于 2006-8-18 06:50:35

个人认为是没有根据的猜测。你那么多的系数,为什么不好好查一查是否出错?比如b4=gama*(delta*dt)就不对,应该是b4=gama/(delta*dt).

realyyy 发表于 2006-8-18 15:57:44

首先感谢tonnyw!
系数已核对并改正!结果没有太大变化啊。。。。。。郁闷中。。。。。。

tonnyw 发表于 2006-8-18 18:55:20

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);

速度公式不对,记住先用位移求加速度,用得到的加速度求速度。

realyyy 发表于 2006-8-21 12:36:26

把下面两句的顺序换一下对结果没有影响啊:
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 发表于 2006-8-21 22:13:56

请仔细检查一下,我说了你的速度公式不对,顺序不对。你现在调整了顺序,但公式还是不对。

realyyy 发表于 2006-8-23 14:03:32

感谢tonnyw的热心关注!
我检查了啊:用得到的加速度求速度的式子整理后和x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i)是一样的。

tonnyw 发表于 2006-8-23 15:15:51

x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i)你这里系数前面都是负号,很明显不对。仔细对照一下我贴的公式,就知道错在什么地方了。

realyyy 发表于 2006-8-24 14:36:03

再次感谢tonnyw!
我又认真、仔细的核对了一下所有的系数,发现了问题:我是按照《有限单元法》(薛守义编著)中的公式写的程序,系数与其书中所述完全一致,但与你所贴的就不一致。不知你用的是哪本书中的公式?而且更惊奇的是:我发现《机械振动基础》(胡海岩编著)一书中所述的公式与你我所用公式也不一样!!
到底大家对NEWMARK公式是怎么搞的吗?!!只能有一个是对的。。。。。。其他的全错!
更加郁闷中。。。。。。

realyyy 发表于 2006-8-24 15:25:07

最后经本人验证:tonnyw的贴图中的公式完全正确!程序运行后结果相当准确。
(还是外国的公式准!没办法。)
本讨论至此完美结束。再次感谢tonnyw的热心帮助!

kuki29 发表于 2006-8-24 19:21:00

求助随机有限元的程序

求助随机有限元的程序
本人一点头绪都没有
应该如何开始做啊

molen 发表于 2006-9-9 10:06:51

国内的大专家们写的书,
sign,,,

jacobi 发表于 2006-9-24 11:43:44

tonnyw 贴的是大牛K.J.Bathe的书,呵呵

ccclgz_cast 发表于 2006-9-26 14:54:35

不是吧!
国内的人写书来骗人啦
真是郁闷啦,把我们这些人给害掺了不是!

pp786195 发表于 2009-11-9 14:47:10

问一下,MATLAB计算的每一步结果,需要判断收敛吗,怎么判断?

孤叶 发表于 2011-1-16 18:49:08

静力分析不需要判断,如果存在材料几何的非线性的话要进行矩阵的更新和迭代分析 15# pp786195 jing

lucieze 发表于 2011-1-18 05:49:26

厉害啊 学习了

聆听随风 发表于 2011-5-1 12:16:23

学习了。。

xiaojia86230910 发表于 2011-5-14 10:13:07

都是牛人啊!学习!!!!

07mumu0923 发表于 2018-6-30 16:18:58

请问质量矩阵和刚度矩阵,阻尼矩阵是怎么得到的呢?采用模态分析还是数学建模计算?
页: [1]
查看完整版本: MATLAB之newmark法计算的结构响应与ANSYS的结果差异为何如此之大?