solar617 发表于 2009-3-10 20:05

求matlab三维有限元patch test程序(六面体单元)

急~
求有限元单元法求构件中个点位移、应变、应力。
采用六面体等参单元。:@)

ZXinNCHU 发表于 2014-4-17 09:52

我最关心怎么求应力,是不是也要通过高斯积分

牛小贱 发表于 2014-4-17 21:38

给楼主一个关于“三节点平面应力单元有限元分析,计算出了其节点力、位移和应力”的程序。程序代码:%清空数据
clear;
%取杨氏弹性模量为210000,泊松比为0.3
E=210000;
mu=0.3;
%板厚:
h=0.003;
%此为平面应力问题,弹性矩阵为:
D=E/(1-mu^2)*;
%所有节点的坐标
XY=zeros(121,2);
for i=1:11
    for j=1:11
      XY(11*(i-1)+j,:)=[-0.75+0.15*(j-1),-1+0.2*(i-1)];
    end
end
%所有单元内部的节点排序
EPI=zeros(200,3);
for i=1:10
    for j=1:10
      EPI(20*(i-1)+2*j-1,:)=;
      EPI(20*(i-1)+2*j,:)=;
    end
end
%受力节点的信息:
FPI=zeros(11,3);
for i=1:11
    FPI(i,:)=;
end
%受约束节点的信息:
BPI=zeros(11,3);
for i=1:11
    BPI(i,:)=;
end
%画出网格图形,用星号标记受力点,用叉号标记受约束点
hold on
for i=1:size(EPI,1)
    plot(, ...
         );
end
for i=1:size(FPI,1);
    plot(XY(FPI(i,1),1),XY(FPI(i,1),2),'r*','MarkerSize',15,'LineWidth',3);
end
for i=1:size(BPI,1);
    plot(XY(BPI(i,1),1),XY(BPI(i,1),2),'kx','MarkerSize',15,'LineWidth',3);
end
axis equal
hold off
%下面整个for循环计算其总刚度矩阵
K=zeros(2*size(XY,1));
for i=1:size(EPI,1)
    x1=XY(EPI(i,1),1);
    y1=XY(EPI(i,1),2);
    x2=XY(EPI(i,2),1);
    y2=XY(EPI(i,2),2);
    x3=XY(EPI(i,3),1);
    y3=XY(EPI(i,3),2);
    A=det()/2;
    syms x y
    N1=((x2*y3-x3*y2)+(y2-y3)*x+(x3-x2)*y)/2/A;
    N2=((x3*y1-x1*y3)+(y3-y1)*x+(x1-x3)*y)/2/A;
    N3=((x1*y2-x2*y1)+(y1-y2)*x+(x2-x1)*y)/2/A;
    %三角形单元的几何矩阵
    B=[];
    B=[diff(N1,'x') 0 diff(N2,'x') 0 diff(N3,'x') 0
       0 diff(N1,'y') 0 diff(N2,'y') 0 diff(N3,'y')
       diff(N1,'y') diff(N1,'x') diff(N2,'y') ...
       diff(N2,'x') diff(N3,'y') diff(N3,'x')];
   B=double(B);
   BB(3*i-2:3*i,1:6)=B;
   ke=zeros(6);
   ke=h*transpose(B)*D*B*A;
   P=;
   for j=1:3
       for m=1:3
         for o=1:2
               for p=1:2
                   K(2*(P(j)-1)+o,2*(P(m)-1)+p)=K(2*(P(j)-1)+o, ...
                   2*(P(m)-1)+p)+ke(2*(j-1)+o,2*(m-1)+p);
               end
         end
       end
   end
end

t=1;
j=1;
l=1;
a=zeros(1,2*size(XY,1)-2*size(BPI,1));
for i=1:size(XY,1)
    if i==BPI(j,1)
      j=j+1;
      if j>size(BPI,1)
            j=j-1;
      end
    elseif i==FPI(l,1)
            fp(l)=t;
            for k=1:2
                a(t)=2*(i-1)+k;
                t=t+1;
            end
            l=l+1;
            if l>size(FPI,1)
                l=l-1;
            end
    else
      for k=1:2
            a(t)=2*(i-1)+k;
            t=t+1;
      end
    end
end
k=zeros(length(a));
for i=1:length(a)
    for j=1:length(a)
      k(i,j)=K(a(i),a(j));
    end
end
f=zeros(length(a),1);
for i=1:l
    f(fp(i))=FPI(i,2);
    f(fp(i)+1)=FPI(i,3);
end
u=k\f;
U=zeros(2*size(XY,1),1);
for i=1:length(a)
    U(a(i),1)=u(i);
end
F=K*U;
%计算xy方向的应力和主应力
for i=1:size(EPI,1)
    P=;
    ue=[U(2*P(1)-1);U(2*P(1))
      U(2*P(2)-1);U(2*P(2))
      U(2*P(3)-1);U(2*P(3))];
    sigma(i,:)=D*BB((3*i-2):3*i,1:6)*ue;
    R=(sigma(i,1)+sigma(i,2))/2;
    Q=((sigma(i,1)-sigma(i,2))/2)^2+sigma(i,3)*sigma(i,3);
    M=2*sigma(i,3)/(sigma(i,1)-sigma(i,2));
    s1=R+sqrt(Q);
    s2=R-sqrt(Q);
    theta=(atan(M)/2)*180/pi;
    yl(i,:)=;
end
%计算最右上角节点的位移:
ys=
%计算第60个节点(-0.15,0)上的位移和应力:
wy=
yl=(sigma(87,:)+sigma(88,:)+sigma(89,:)+sigma(108,:)+sigma(109,:)+sigma(110,:))/6希望对LZ有所帮助!!



页: [1]
查看完整版本: 求matlab三维有限元patch test程序(六面体单元)