资源描述：
《平面四节点等参单元matlab实现.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、计算力学报告平面四节点等参单元学生姓名：朱彬学号：一、问题描述及分析在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：打开ansys，在单元类型中选择sol
2、id->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','
3、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=
4、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
5、,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.,-0.,0,0.,0.];%5*5高斯积分点Ett=[-0.,-0.,0,0.,0.];H=[0.,0.,0.,0.,0.];%高斯积分权系数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(
6、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))=
7、1e20;end%置大数法f=zeros(2*n,1);%初始载荷矩阵f(10)=-10e3;%加载荷10kNU=Kf;%节点位移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
8、,-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
