声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2792|回复: 2

[编程技巧] 求matlab三维有限元patch test程序(六面体单元)

[复制链接]
发表于 2009-3-10 20:05 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

使用道具 举报

发表于 2014-4-17 09:52 | 显示全部楼层
我最关心怎么求应力,是不是也要通过高斯积分
发表于 2014-4-17 21:38 | 显示全部楼层
给楼主一个关于“三节点平面应力单元有限元分析,计算出了其节点力、位移和应力”的程序。程序代码:
  1. %清空数据
  2. clear;
  3. %取杨氏弹性模量为210000,泊松比为0.3
  4. E=210000;
  5. mu=0.3;
  6. %板厚:
  7. h=0.003;
  8. %此为平面应力问题,弹性矩阵为:
  9. D=E/(1-mu^2)*[1 mu 0;mu 1 0;0 0 (1-mu)/2];
  10. %所有节点的坐标
  11. XY=zeros(121,2);
  12. for i=1:11
  13.     for j=1:11
  14.         XY(11*(i-1)+j,:)=[-0.75+0.15*(j-1),-1+0.2*(i-1)];
  15.     end
  16. end
  17. %所有单元内部的节点排序
  18. EPI=zeros(200,3);
  19. for i=1:10
  20.     for j=1:10
  21.         EPI(20*(i-1)+2*j-1,:)=[11*(i-1)+j 11*i+j+1 11*i+j];
  22.         EPI(20*(i-1)+2*j,:)=[11*(i-1)+j 11*(i-1)+j+1 11*i+j+1];
  23.     end
  24. end
  25. %受力节点的信息:
  26. FPI=zeros(11,3);
  27. for i=1:11
  28.     FPI(i,:)=[11*i 200/11 0];
  29. end
  30. %受约束节点的信息:
  31. BPI=zeros(11,3);
  32. for i=1:11
  33.     BPI(i,:)=[11*(i-1)+1 0 0];
  34. end
  35. %画出网格图形,用星号标记受力点,用叉号标记受约束点
  36. hold on
  37. for i=1:size(EPI,1)
  38.     plot([XY(EPI(i,1),1),XY(EPI(i,2),1),XY(EPI(i,3),1),XY(EPI(i,1),1)], ...
  39.          [XY(EPI(i,1),2),XY(EPI(i,2),2),XY(EPI(i,3),2),XY(EPI(i,1),2)]);
  40. end
  41. for i=1:size(FPI,1);
  42.     plot(XY(FPI(i,1),1),XY(FPI(i,1),2),'r*','MarkerSize',15,'LineWidth',3);
  43. end
  44. for i=1:size(BPI,1);
  45.     plot(XY(BPI(i,1),1),XY(BPI(i,1),2),'kx','MarkerSize',15,'LineWidth',3);
  46. end
  47. axis equal
  48. hold off
  49. %下面整个for循环计算其总刚度矩阵
  50. K=zeros(2*size(XY,1));
  51. for i=1:size(EPI,1)
  52.     x1=XY(EPI(i,1),1);
  53.     y1=XY(EPI(i,1),2);
  54.     x2=XY(EPI(i,2),1);
  55.     y2=XY(EPI(i,2),2);
  56.     x3=XY(EPI(i,3),1);
  57.     y3=XY(EPI(i,3),2);
  58.     A=det([1 x1 y1;1 x2 y2;1 x3 y3])/2;
  59.     syms x y
  60.     N1=((x2*y3-x3*y2)+(y2-y3)*x+(x3-x2)*y)/2/A;
  61.     N2=((x3*y1-x1*y3)+(y3-y1)*x+(x1-x3)*y)/2/A;
  62.     N3=((x1*y2-x2*y1)+(y1-y2)*x+(x2-x1)*y)/2/A;
  63.     %三角形单元的几何矩阵
  64.     B=[];
  65.     B=[diff(N1,'x') 0 diff(N2,'x') 0 diff(N3,'x') 0
  66.        0 diff(N1,'y') 0 diff(N2,'y') 0 diff(N3,'y')
  67.        diff(N1,'y') diff(N1,'x') diff(N2,'y') ...
  68.        diff(N2,'x') diff(N3,'y') diff(N3,'x')];
  69.    B=double(B);
  70.    BB(3*i-2:3*i,1:6)=B;
  71.    ke=zeros(6);
  72.    ke=h*transpose(B)*D*B*A;
  73.    P=[EPI(i,1) EPI(i,2) EPI(i,3)];
  74.    for j=1:3
  75.        for m=1:3
  76.            for o=1:2
  77.                for p=1:2
  78.                    K(2*(P(j)-1)+o,2*(P(m)-1)+p)=K(2*(P(j)-1)+o, ...
  79.                    2*(P(m)-1)+p)+ke(2*(j-1)+o,2*(m-1)+p);
  80.                end
  81.            end
  82.        end
  83.    end
  84. end

  85. t=1;
  86. j=1;
  87. l=1;
  88. a=zeros(1,2*size(XY,1)-2*size(BPI,1));
  89. for i=1:size(XY,1)  
  90.     if i==BPI(j,1)
  91.         j=j+1;
  92.         if j>size(BPI,1)
  93.             j=j-1;
  94.         end
  95.     elseif i==FPI(l,1)
  96.             fp(l)=t;
  97.             for k=1:2
  98.                 a(t)=2*(i-1)+k;
  99.                 t=t+1;
  100.             end
  101.             l=l+1;
  102.             if l>size(FPI,1)
  103.                 l=l-1;
  104.             end
  105.     else
  106.         for k=1:2
  107.             a(t)=2*(i-1)+k;
  108.             t=t+1;
  109.         end
  110.     end
  111. end
  112. k=zeros(length(a));
  113. for i=1:length(a)
  114.     for j=1:length(a)
  115.         k(i,j)=K(a(i),a(j));
  116.     end
  117. end
  118. f=zeros(length(a),1);
  119. for i=1:l
  120.     f(fp(i))=FPI(i,2);
  121.     f(fp(i)+1)=FPI(i,3);
  122. end
  123. u=k\f;
  124. U=zeros(2*size(XY,1),1);
  125. for i=1:length(a)
  126.     U(a(i),1)=u(i);
  127. end
  128. F=K*U;
  129. %计算xy方向的应力和主应力
  130. for i=1:size(EPI,1)
  131.     P=[EPI(i,1),EPI(i,2),EPI(i,3)];
  132.     ue=[U(2*P(1)-1);U(2*P(1))
  133.         U(2*P(2)-1);U(2*P(2))
  134.         U(2*P(3)-1);U(2*P(3))];
  135.     sigma(i,:)=D*BB((3*i-2):3*i,1:6)*ue;
  136.     R=(sigma(i,1)+sigma(i,2))/2;
  137.     Q=((sigma(i,1)-sigma(i,2))/2)^2+sigma(i,3)*sigma(i,3);
  138.     M=2*sigma(i,3)/(sigma(i,1)-sigma(i,2));
  139.     s1=R+sqrt(Q);
  140.     s2=R-sqrt(Q);
  141.     theta=(atan(M)/2)*180/pi;
  142.     yl(i,:)=[s1 s2 theta];
  143. end
  144. %计算最右上角节点的位移:
  145. ys=[U(241) U(242)]
  146. %计算第60个节点(-0.15,0)上的位移和应力:
  147. wy=[U(119) U(120)]
  148. yl=(sigma(87,:)+sigma(88,:)+sigma(89,:)+sigma(108,:)+sigma(109,:)+sigma(110,:))/6
复制代码
希望对LZ有所帮助!!



您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-12 07:44 , Processed in 0.061860 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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