找回密码
 注册
Simdroid-非首页
查看: 940|回复: 19

[H. 有限元编程] MATLAB之newmark法计算的结构响应与ANSYS的结果差异为何如此之大?

[复制链接]
发表于 2006-8-17 13:59:29 | 显示全部楼层 |阅读模式 来自 陕西西安
请问高手:我用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=[0 0 0 0 0 0]';
v=[0 0 0 0 0 0]';
a=[10 0 0 0 0 0]';
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)=[f,0,0,0,0,0]';
      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=[x(1:1,t)';x(2:2,t)';x(3:3,t)';x(4:4,t)';x(5:5,t)';x(6:6,t)']
%x1respons=[x1(1:1,t)';x1(2:2,t)';x1(3:3,t)';x1(4:4,t)';x1(5:5,t)';x1(6:6,t)']
%x2respons=[x2(1:1,t)';x2(2:2,t)';x2(3:3,t)';x2(4:4,t)';x2(5:5,t)';x2(6:6,t)']
%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 编辑 ]

评分

1

查看全部评分

发表于 2006-8-18 06:50:35 | 显示全部楼层 来自 美国
Simdroid开发平台
个人认为是没有根据的猜测。你那么多的系数,为什么不好好查一查是否出错?比如b4=gama*(delta*dt)就不对,应该是b4=gama/(delta*dt).

评分

1

查看全部评分

 楼主| 发表于 2006-8-18 15:57:44 | 显示全部楼层 来自 陕西西安
首先感谢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);

速度公式不对,记住先用位移求加速度,用得到的加速度求速度。
 楼主| 发表于 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);
........................
期盼高人啊!
发表于 2006-8-21 22:13:56 | 显示全部楼层 来自 美国
请仔细检查一下,我说了你的速度公式不对,顺序不对。你现在调整了顺序,但公式还是不对。
 楼主| 发表于 2006-8-23 14:03:32 | 显示全部楼层 来自 陕西西安
感谢tonnyw的热心关注!
我检查了啊:用得到的加速度求速度的式子整理后和x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i)是一样的。
发表于 2006-8-23 15:15:51 | 显示全部楼层 来自 美国
x1(:,i+1)=b4*(x(:,i+1)-x(:,i))-b5*x1(:,i)-b6*x2(:,i)你这里系数前面都是负号,很明显不对。仔细对照一下我贴的公式,就知道错在什么地方了。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
 楼主| 发表于 2006-8-24 14:36:03 | 显示全部楼层 来自 陕西西安
再次感谢tonnyw!
我又认真、仔细的核对了一下所有的系数,发现了问题:我是按照《有限单元法》(薛守义编著)中的公式写的程序,系数与其书中所述完全一致,但与你所贴的就不一致。不知你用的是哪本书中的公式?而且更惊奇的是:我发现《机械振动基础》(胡海岩编著)一书中所述的公式与你我所用公式也不一样!!
到底大家对NEWMARK公式是怎么搞的吗?!!只能有一个是对的。。。。。。其他的全错!
更加郁闷中。。。。。。
 楼主| 发表于 2006-8-24 15:25:07 | 显示全部楼层 来自 陕西西安
最后经本人验证:tonnyw的贴图中的公式完全正确!程序运行后结果相当准确。
(还是外国的公式准!没办法。)
本讨论至此完美结束。再次感谢tonnyw的热心帮助!
发表于 2006-8-24 19:21:00 | 显示全部楼层 来自 辽宁锦州

求助随机有限元的程序

求助随机有限元的程序
本人一点头绪都没有
应该如何开始做啊
发表于 2006-9-9 10:06:51 | 显示全部楼层 来自 加拿大
国内的大专家们写的书,
sign,,,
发表于 2006-9-24 11:43:44 | 显示全部楼层 来自 上海
tonnyw 贴的是大牛K.J.Bathe的书,呵呵
发表于 2006-9-26 14:54:35 | 显示全部楼层 来自 北京
不是吧!
国内的人写书来骗人啦
真是郁闷啦,把我们这些人给害掺了不是!
发表于 2009-11-9 14:47:10 | 显示全部楼层 来自 江苏苏州
问一下,MATLAB计算的每一步结果,需要判断收敛吗,怎么判断?
回复 不支持

使用道具 举报

发表于 2011-1-16 18:49:08 | 显示全部楼层 来自 上海
静力分析不需要判断,如果存在材料几何的非线性的话要进行矩阵的更新和迭代分析 15# pp786195 jing
回复 不支持

使用道具 举报

发表于 2011-1-18 05:49:26 | 显示全部楼层 来自 法国
厉害啊 学习了
回复 不支持

使用道具 举报

发表于 2011-5-1 12:16:23 | 显示全部楼层 来自 福建福州
学习了。。
回复 不支持

使用道具 举报

发表于 2011-5-14 10:13:07 | 显示全部楼层 来自 湖北武汉
都是牛人啊!学习!!!!
回复 不支持

使用道具 举报

发表于 2018-6-30 16:18:58 | 显示全部楼层 来自 湖南长沙
请问质量矩阵和刚度矩阵,阻尼矩阵是怎么得到的呢?采用模态分析还是数学建模计算?
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-4-26 00:39 , Processed in 0.053747 second(s), 17 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表