一个关于平面有限元的问题（在线等）

xyxyclq

平面有限元的，运行结果总是NAN，麻烦给看看！ assemblestiffnessmatrix.m (380 Bytes, 下载次数: 12)

2007-7-21 17:12 上传

点击文件名下载附件

mainpro3.m (1.9 KB, 下载次数: 11)

2007-7-21 17:12 上传

点击文件名下载附件

stiffnessmatrix.m (889 Bytes, 下载次数: 12)

2007-7-21 17:12 上传

点击文件名下载附件

输入的数据.txt (362 Bytes, 下载次数: 12)

2007-7-21 17:12 上传

点击文件名下载附件

要交作业，所以着急

[ 本帖最后由 eight 于 2007-7-23 12:55 编辑 ]

xjzuo

1. 你的运行,及其出错信息没有讲清楚;
2. 不要灌水, 没人愿意回答的话,你一直重复也没有用.

花如月

如果没人理就考虑是不是你问问题的方式有问题，没有错误提示根本没提到你遇到的具体问题。整个传了一堆文件，摆明了让别人帮你做作业呀？这里是讨论版，要弄清楚哦

[ 本帖最后由 eight 于 2007-7-23 12:55 编辑 ]

xyxyclq

不是，我自己编好了，能运行，不过结果里面有好多NAN。楼上的说的我都惭愧死了，虽说能力不怎么样，但也不是让你高手帮着做作业！

花如月

不好意思，言语不当之处我向你道歉，我完全没有你说的那个意思（无心之过）。不过问问题确实不该这样。既然是讨论版，问题直观就最好了，除非程序特别大才用附件，所以直接贴上来会更直观。别人帮你看程序，排查错误也更方便些

xyxyclq

按花主任的说法我把程序贴出来了，帮看看问题出那边吧
%mainprogram
global node element K  E 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[0 0;1 0;2 0;0 1;1 1;1 2]
input element number:4
input element[1 1 5 4;2 1 2 5;3 2 6 5;4 2 3 6]
input number of boundary condition3
input boundary condition:[1 1 0;1 2 0;4 1 0]
input number of node force1
input node force[6 1 1]
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）问题是否可以简化一下呢？先对一个简单些的问题按你的思路做一下，可以查看每一步的计算结果是否正常。这个程序参数太多，调试比较麻烦。
一家之言了，能力有限，能想到的也就是先简化问题。希望路过的高手可以帮你早日解决问题

xyxyclq

谢谢花主任，不过有限元输入的是很多，也不知道改减那些，路过得高手帮忙看看吧

w89986581

根据警告先检查一下A的值，再确认一下你的刚度矩阵和质量矩阵是否是正定的。

xyxyclq

Ellips，里面的刚度矩阵是对的，看了f语言的和其他人都是这样编，郁闷了

eight

原帖由 xyxyclq 于 2007-7-23 12:41 发表

                               
登录/注册后可看大图

Ellips，里面的刚度矩阵是对的，看了f语言的和其他人都是这样编，郁闷了

难度别人这样做就不需要检查了吗？那别人做错了怎办（也许是笔误、或者有意隐藏真相）？建议先搞清楚，变成自己的东西，然后再作下一步的改进。就算你拿来就用，也不要紧，现在遇到的问题是涉及你的专业，所以，建议还是自己一步一步调试一下最好，否则究竟是错是对别人根本不知道