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

计算力学报告平面四节点等参单元学生姓名：朱彬学号：20100311一、问题描述及分析在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：打开ansys，在单元类型中选择solid-Quad4node182单元；建立类半径为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','elements');%读取单元编号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);%初始总体刚度矩阵fori=1:msymsKsEtxyI1I2abB;%定义可能存在的变量e=[1,1;-1,1;-1,-1;1,-1];forj=1:4%形函数N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;forj=1:4%标准母单元映射到真实单元x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(elem(i,j),2);endJ1=jacobian([x;y],[Ks;Et]);%雅克比矩阵及其转置J=J1';forj=1:4I1=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;B(3,2*j)=a;endD=(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];%高斯积分权系数forj=1:5%高斯积分求单元刚度阵forl=1:5Ks=Kss(j);Et=Ett(l);B=subs(B);J=subs(J);k=k+H(j)*H(l)*B'*D*B*det(J);endendG=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);fori=1:bK(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=K\f;%节点位移fori=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))];endfori=1:m%求单元应力symsKsEtxyI1I2abB;e=[1,1;-1,1;-1,-1;1,-1];forj=1:4N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;forj=1:4x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(elem(i,j),2);endJ1=jacobian([x;y],[Ks;Et]);J=J1';forj=1:4I1=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*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阵endD=(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,Et},{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,e);%半径方向节点处的mises应力msigmat=zeros(1,f);%圆周方向节点处的mises应力fori=1:e%绕节点平均g=0;forj=1:m%判断节点在单元的那个位置并加上相应的应力值ifelem(j,1)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma1(:,j);endifelem(j,2)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma2(:,j);endifelem(j,3)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma3(:,j);endifelem(j,4)-c(i)==0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma4(:,j);endendsigmar(:,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应力fori=1:f%同上g=0;forj=1:mifelem(j,1)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma1(:,j);endifelem(j,2)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma2(:,j);endifelem(j,3)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma3(:,j);endifelem(j,4)-d(i)==0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma4(:,j);endendsigmat(:,i)=sigmat(:,i)/g;msigmat(:,i)=(0.5*((sigmat(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中红线标注：图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,185]';e=size(c,1)f=size(d,1);fori=1:eansysr(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;endfori=1:fansyst(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*Columns1through90.03882.24230.09690.04390.05210.07020.11000.14440.2330Columns10through110.47660.5914ansyst=1.0e+003*Columns1through94.76552.35040.98260.84540.61972.03220.70011.01871.2914Columns10through181.33061.25161.06671.07750.61655.40944.42631.21371.3410Columns19through200.56890.77002.等参单元编程计算结果msigmar=1.0e+004*Columns1through90.03653.42010.04580.04150.04480.06730.13960.10760.1859Columns10through110.44033.2183msigmat=1.0e+004*Columns1throu