平面四节点等参单元matlab实现

1、计算力学报告平面四节点等参单元学生姓名： 朱 彬学号： 20100311 一、 问题描述及分析在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：打开ansys，在单元类型中选择solid->Quad 4

2、 node 182单元；建立类半径为0.5外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。如下图： 图1三、有限元程序及求解程序求解使用了matlab语言。具体如下：程序：clcclearE=2e11; %弹性模量NU=0.3; %泊松比t=0.1; %厚度X=xlsread('D:data','nodes'); %读取节点坐标elem=xlsread('D:data','el

3、ements'); %读取单元编号w=1,2,3,4,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28; %有位移约束的节点n=size(X,1); %节点数m=size(elem,1); %单元数K=zeros(2*n); %初始总体刚度矩阵for i=1:m syms Ks Et x y I1 I2 a b B; %定义可能存在的变量 e=1,1;-1,1;-1,-1;1,-1; for j=1:4 %形函数 N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et); end x=0;y=0; f

4、or j=1:4 %标准母单元映射到真实单元 x=x+N(j)*X(elem(i,j),1); y=y+N(j)*X(elem(i,j),2); end J1=jacobian(x;y,Ks;Et); %雅克比矩阵及其转置 J=J1' for j=1:4 I1=diff(N(j),Ks); %形函数分别对Ks和Et的偏导数 I2=diff(N(j),Et); C=(J-1)*I1;I2; a=C(1); %形函数对x，y的偏导数 b=C(2); B(1,2*j-1)=a; %组成B阵 B(1,2*j)=0; B(2,2*j-1)=0; B(2,2*j)=b; B(3,2*j-1)=b;

5、 B(3,2*j)=a; end D=(E/(1-NU*NU)*1,NU,0;NU,1,0;0,0,(1-NU)/2; %D阵 k=zeros(8,8); Kss=-0.906179,-0.538469,0,0.538469,0.906179; %5*5高斯积分点 Ett=-0.906179,-0.538469,0,0.538469,0.906179; H=0.236926,0.478628,0.568888,0.478628,0.236926;%高斯积分权系数 for j=1:5 %高斯积分求单元刚度阵 for l=1:5 Ks=Kss(j);Et=Ett(l); B=subs(B); J=

6、subs(J); k=k+H(j)*H(l)*B'*D*B*det(J); end end G=zeros(8,2*n); %初始总刚变换矩阵 G(1,2*elem(i,1)-1)=1; %总刚变换矩阵 G(2,2*elem(i,1)=1; G(3,2*elem(i,2)-1)=1; G(4,2*elem(i,2)=1; G(5,2*elem(i,3)-1)=1; G(6,2*elem(i,3)=1; G(7,2*elem(i,4)-1)=1; G(8,2*elem(i,4)=1; K=K+G'*k*G; %总体刚度矩阵合成endKK=K;b=size(w,1);for i=1

7、:b K(2*w(i)-1,2*w(i)-1)=1e20; K(2*w(i),2*w(i)=1e20;end %置大数法f=zeros(2*n,1); %初始载荷矩阵f(10)=-10e3; %加载荷10kNU=Kf; %节点位移for i=1:m %将每个单元各个节点位移集合 u(:,i)=U(2*elem(i,1)-1);U(2*elem(i,1);U(2*elem(i,2)-1);U(2*elem(i,2);U(2*elem(i,3)-1);U(2*elem(i,3);U(2*elem(i,4)-1);U(2*elem(i,4);end for i=1:m %求单元应力 syms Ks

8、Et x y I1 I2 a b B; e=1,1;-1,1;-1,-1;1,-1; for j=1:4 N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et); end x=0;y=0; for j=1:4 x=x+N(j)*X(elem(i,j),1); y=y+N(j)*X(elem(i,j),2); end J1=jacobian(x;y,Ks;Et); J=J1' for j=1:4 I1=diff(N(j),Ks); I2=diff(N(j),Et); C=(J-1)*I1;I2; a=C(1); b=C(2); B(1,2*j-1)=a; B(1,2*

9、j)=0; B(2,2*j-1)=0; B(2,2*j)=b; B(3,2*j-1)=b; B(3,2*j)=a; %以上同前面部分为得到B阵 end D=(E/(1-NU*NU)*1,NU,0;NU,1,0;0,0,(1-NU)/2; w=D*B*u(:,i); w1=subs(w,Ks,Et,1,1); %求单元上各节点的应力 sigma1(:,i)=double(w1); w2=subs(w,Ks,Et,-1,1); sigma2(:,i)=double(w2); w3=subs(w,Ks,Et,-1,-1); sigma3(:,i)=double(w3); w4=subs(w,Ks,E

10、t,1,-1); sigma4(:,i)=double(w4); endc=24,29,47,58,78,79,137,149,160,166,186' %如截图选取圆半径方向的节点号d=166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185'%圆周方向选择的节点号e=size(c,1); f=size(d,1);sigmar=zeros(3,e); %r相同角度不同节点应力矩阵sigmat=zeros(3,f); %角度不同r不同节点应力矩阵msigmar=zeros(1,

11、e); %半径方向节点处的mises应力msigmat=zeros(1,f); %圆周方向节点处的mises应力for i=1:e %绕节点平均 g=0; for j=1:m %判断节点在单元的那个位置并加上相应的应力值 if elem(j,1)-c(i)=0 g=g+1; sigmar(:,i)=sigmar(:,i)+sigma1(:,j); end if elem(j,2)-c(i)=0 g=g+1; sigmar(:,i)=sigmar(:,i)+sigma2(:,j); end if elem(j,3)-c(i)=0 g=g+1; sigmar(:,i)=sigmar(:,i)+si

12、gma3(:,j); end if elem(j,4)-c(i)=0 g=g+1; sigmar(:,i)=sigmar(:,i)+sigma4(:,j); end end sigmar(:,i)=sigmar(:,i)/g; %求应力绕节点平均 msigmar(:,i)=(0.5*(sigmar(1,i)-sigmar(2,i)2+sigmar(1,i)2+sigmar(2,i)2+6*(sigmar(3,i)2)0.5; %求节点处的mises应力endmsigmar %mises应力for i=1:f %同上 g=0; for j=1:m if elem(j,1)-d(i)=0 g=g+

13、1; sigmat(:,i)=sigmat(:,i)+sigma1(:,j); end if elem(j,2)-d(i)=0 g=g+1; sigmat(:,i)=sigmat(:,i)+sigma2(:,j); end if elem(j,3)-d(i)=0 g=g+1; sigmat(:,i)=sigmat(:,i)+sigma3(:,j); end if elem(j,4)-d(i)=0 g=g+1; sigmat(:,i)=sigmat(:,i)+sigma4(:,j); end end sigmat(:,i)=sigmat(:,i)/g; msigmat(:,i)=(0.5*(si

14、gmat(1,i)-sigmat(2,i)2+sigmat(1,i)2+sigmat(2,i)2+6*(sigmat(3,i)2)0.5;endmsigmat四、计算结果：1.ANSYS软件计算结果计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。选取半径方向的单元编号为： c=24,29,47,58,78,79,137,149,160,166,186'选取圆周方向上的单元编号为：d=166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185'其在图中的位置为图1中红线标

15、注： 图2在处理ANSYS计算结果时，我导出了节点x，y方向的正应力和剪应力，将其导入excle中并去掉无用项如图2所示，并编写matlab程序将选取的点的mises应力求出来。图2，求节点mises应力的程序以及计算结果如下： 图3求选取点的mises应力程序：clcA=xlsread('D:data','ansys');c=24,29,47,58,78,79,137,149,160,166,186'd=166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184

16、,185'e=size(c,1)f=size(d,1);for i=1:e ansysr(i)=(0.5*(A(c(i),1)-A(c(i),2)2+A(c(i),1)2+A(c(i),2)2+6*(A(c(i),3)2)0.5;endfor i=1:f ansyst(i)=(0.5*(A(d(i),1)-A(d(i),2)2+A(d(i),1)2+A(d(i),2)2+6*(A(d(i),3)2)0.5;endansysransyst计算结果如下：ansysr = 1.0e+004 * Columns 1 through 9 0.0388 2.2423 0.0969 0.0439 0

17、.0521 0.0702 0.1100 0.1444 0.2330 Columns 10 through 11 0.4766 0.5914ansyst = 1.0e+003 * Columns 1 through 9 4.7655 2.3504 0.9826 0.8454 0.6197 2.0322 0.7001 1.0187 1.2914 Columns 10 through 18 1.3306 1.2516 1.0667 1.0775 0.6165 5.4094 4.4263 1.2137 1.3410 Columns 19 through 20 0.5689 0.77002.等参单元编程

18、计算结果msigmar = 1.0e+004 * Columns 1 through 9 0.0365 3.4201 0.0458 0.0415 0.0448 0.0673 0.1396 0.1076 0.1859 Columns 10 through 11 0.4403 3.2183msigmat = 1.0e+004 * Columns 1 through 9 0.4403 0.3125 0.2060 0.1686 0.2256 0.6079 0.2483 0.2184 0.1692 Columns 10 through 18 0.1140 0.2197 0.1035 0.1069 0.0854 1.5036 0.4233 0.1903 0.1553 Columns 19 through 20 0.2950 0.2027五、讨论比较对所选取的节点的mises应力列表做对比：1.沿半径方向（排列顺序为节点号从小到大） 表1节点号242947587879137Ansys解0.03882.24230.09690.04390.05210