一个关于平面有限元的问题（在线等）
平面有限元的，运行结果总是NAN，麻烦给看看！要交作业，所以着急
[ 本帖最后由 eight 于 2007-7-23 12:55 编辑 ] 1. 你的运行,及其出错信息没有讲清楚;
2. 不要灌水, 没人愿意回答的话,你一直重复也没有用. 如果没人理就考虑是不是你问问题的方式有问题，没有错误提示根本没提到你遇到的具体问题。整个传了一堆文件，摆明了让别人帮你做作业呀？这里是讨论版，要弄清楚哦
[ 本帖最后由 eight 于 2007-7-23 12:55 编辑 ] 不是，我自己编好了，能运行，不过结果里面有好多NAN。楼上的说的我都惭愧死了，虽说能力不怎么样，但也不是让你高手帮着做作业！
回复 #6 xyxyclq 的帖子
不好意思，言语不当之处我向你道歉，我完全没有你说的那个意思（无心之过）。不过问问题确实不该这样。既然是讨论版，问题直观就最好了，除非程序特别大才用附件，所以直接贴上来会更直观。别人帮你看程序，排查错误也更方便些 按花主任的说法我把程序贴出来了，帮看看问题出那边吧%mainprogram
global node element KE t mu
node_number=input('input node number:');
node=zeros(node_number,2);
node=input('input coordinates of the node');
element_number=input('input element number:');
element=zeros(element_number,4);
element=input('input element') ;
bc_number=input('input number of boundary condition');
bc=zeros(bc_number,3);
bc=input('input boundary condition:');
nf_number=input('input number of node force');
nf=zeros(nf_number,3);
nf=input('input node force');
E=input('input Yang’s modulus');
t=input('input thickness');
mu=input('input Poisson’s ratio');
K=sparse(node_number*2,node_number*2);
f=sparse(node_number*2,1);
for ie=1:element_number;
k=stiffnessmatrix(ie,E,t,mu);
assemblestiffnessmatrix(ie,k);
end;
for inf=1:nf_number;
n=nf(inf,1);
d=nf(inf,2);
f(((n-1)*2+d),1)=nf(inf,3);
end ;
for ibc=1:bc_number
n=bc(ibc,1)
d=bc(ibc,2)
m=(n-1)*2+d
f(m)=bc(ibc,3)*K(m,m)*1e15
K(m,m)=K(m,m)*1e15
end
delta=K\f
el=input('input element number you choose:');
d=zeros(3,6);
b=zeros(3,1);
c=zeros(3,1);
b(1)=node(element(el,2),2)-node(element(el,3),2);
c(1)=node(element(el,3),1)-node(element(el,2),1);
b(2)=node(element(el,3),2)-node(element(el,1),2);
c(2)=node(element(el,1),1)-node(element(el,3),1);
b(3)=node(element(el,1),2)-node(element(el,2),2);
c(3)=node(element(el,2),1)-node(element(el,1),1);
A=(b(1)*c(2)-b(2)*c(1))*0.5;
chang1=E*t/(2*(1-mu^2)*A);
dl=zeros(6,1);
dl(1)=delta(((element(el,1)-1)*2+1),1);
dl(2)=delta(((element(el,1)-1)*2+2),1);
dl(3)=delta(((element(el,2)-1)*2+1),1);
dl(4)=delta(((element(el,2)-1)*2+2),1);
dl(5)=delta(((element(el,3)-1)*2+1),1);
dl(6)=delta(((element(el,3)-1)*2+2),1);
for i=1:3
d(1,(i-1)*2+1)=b(i)*chang1;
d(1,(i-1)*2+2)=mu*c(i)*chang1;
d(2,(i-1)*2+1)=mu*b(i)*chang1;
d(2,(i-1)*2+2)=c(i)*chang1;
d(3,(i-1)*2+1)=0.5*(1-mu)*c(i)*chang1;
d(3,(i-1)*2+2)=0.5*(1-mu)*b(i)*chang1;
end
stress=d*dl
%刚度矩阵
function k=stiffnessmatrix(ie,E,t,mu)
global node element E t mu
k=zeros(6,6)
b=zeros(3,1)
c=zeros(3,1)
b(1)=node(element(ie,2),2)-node(element(ie,3),2)
c(1)=node(element(ie,3),1)-node(element(ie,2),1)
b(2)=node(element(ie,3),2)-node(element(ie,1),2)
c(2)=node(element(ie,1),1)-node(element(ie,3),1)
b(3)=node(element(ie,1),2)-node(element(ie,2),2)
c(3)=node(element(ie,2),1)-node(element(ie,1),1)
A=(b(1)*c(2)-b(2)*c(1))*0.5
chang=E*t/(4*(1-mu^2)*A)
for i=1:3
for j=1:3
m1=(i-1)*2+1
m2=(i-1)*2+2
n1=(j-1)*2+1
n2=(j-1)*2+2
k(m1,n1)=chang*(b(i)*b(j)+0.5*(1-mu)*c(i)*c(j))
k(m1,n2)=chang*(b(i)*c(j)*mu+0.5*(1-mu)*c(i)*b(j))
k(m2,n1)=chang*(b(j)*c(i)*mu+0.5*(1-mu)*c(j)*b(i))
k(m2,n2)=chang*(c(i)*c(j)+0.5*(1-mu)*b(i)*b(j))
end
end
return
%组合
function assemblestiffnessmatrix(ie,k)
global element K
for i=1:1:3
for j=1:1:3
for p=1:1:2
for q=1:1:2
m=(i-1)*2+p
n=(j-1)*2+q
M=(element(ie,i)-1)*2+p
N=(element(ie,j)-1)*2+q
K(M,N)=K(M,N)+k(m,n)
end
end
end
end
return
要输入的参数
input node number:6
input coordinates of the node
input element number:4
input element
input number of boundary condition3
input boundary condition:
input number of node force1
input node force
input Yang’s modulus1.0e+6
input thickness0.1
input Poisson’s ratio0.3
运行出现了如下的结果：
delta =
(1,1) NaN
(2,1) NaN
(3,1) NaN
(4,1) NaN
(5,1) NaN
(6,1) NaN
(7,1) NaN
(8,1) NaN
(9,1) NaN
(10,1) NaN
(11,1) NaN
(12,1) NaN
input element number you choose:1
Warning: Divide by zero.
(Type "warning off MATLAB:divideByZero" to suppress this warning.)
> In C:\Documents and Settings\small_bridge\桌面\有限元\mainpro3.m at line 48
stress =
NaN
NaN
NaN
[ 本帖最后由 xyxyclq 于 2007-7-22 13:30 编辑 ] 帮你看了下，没弄过有限元分析，也找不出什么问题:@L 提2点建议：
（1）刚开始的数据输入可以先改成直接初始化比较好，程序没问题后在用可以交互式的。这样读入数据容易出错。
（2）问题是否可以简化一下呢？先对一个简单些的问题按你的思路做一下，可以查看每一步的计算结果是否正常。这个程序参数太多，调试比较麻烦。
一家之言了，能力有限，能想到的也就是先简化问题。希望路过的高手可以帮你早日解决问题 谢谢花主任，不过有限元输入的是很多，也不知道改减那些，路过得高手帮忙看看吧 根据警告先检查一下A的值，再确认一下你的刚度矩阵和质量矩阵是否是正定的。 Ellips，里面的刚度矩阵是对的，看了f语言的和其他人都是这样编，郁闷了 原帖由 xyxyclq 于 2007-7-23 12:41 发表 http://www.chinavib.com/forum/images/common/back.gif
Ellips，里面的刚度矩阵是对的，看了f语言的和其他人都是这样编，郁闷了
难度别人这样做就不需要检查了吗？那别人做错了怎办（也许是笔误、或者有意隐藏真相）？建议先搞清楚，变成自己的东西，然后再作下一步的改进。就算你拿来就用，也不要紧，现在遇到的问题是涉及你的专业，所以，建议还是自己一步一步调试一下最好，否则究竟是错是对别人根本不知道
页:
[1]