给楼主一个关于“三节点平面应力单元有限元分析,计算出了其节点力、位移和应力”的程序。程序代码:- %清空数据
- 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有所帮助!!
|