求matlab三维有限元patch test程序（六面体单元）

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

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

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

x

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

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)*[1 mu 0;mu 1 0;0 0 (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,:)=[11*(i-1)+j 11*i+j+1 11*i+j];
EPI(20*(i-1)+2*j,:)=[11*(i-1)+j 11*(i-1)+j+1 11*i+j+1];
end
end
%受力节点的信息：
FPI=zeros(11,3);
for i=1:11
FPI(i,:)=[11*i 200/11 0];
end
%受约束节点的信息：
BPI=zeros(11,3);
for i=1:11
BPI(i,:)=[11*(i-1)+1 0 0];
end
%画出网格图形，用星号标记受力点，用叉号标记受约束点
hold on
for i=1:size(EPI,1)
plot([XY(EPI(i,1),1),XY(EPI(i,2),1),XY(EPI(i,3),1),XY(EPI(i,1),1)], ...
[XY(EPI(i,1),2),XY(EPI(i,2),2),XY(EPI(i,3),2),XY(EPI(i,1),2)]);
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([1 x1 y1;1 x2 y2;1 x3 y3])/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=[EPI(i,1) EPI(i,2) EPI(i,3)];
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=[EPI(i,1),EPI(i,2),EPI(i,3)];
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,:)=[s1 s2 theta];
end
%计算最右上角节点的位移：
ys=[U(241) U(242)]
%计算第60个节点（-0.15，0）上的位移和应力：
wy=[U(119) U(120)]
yl=(sigma(87,:)+sigma(88,:)+sigma(89,:)+sigma(108,:)+sigma(109,:)+sigma(110,:))/6

复制代码

希望对LZ有所帮助！！

账号		自动登录	找回密码
密码			我要加入

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

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