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

平面四节点等参单元matlab实现 平面四节点等参单元matlab实现 PAGE 平面四节点等参单元matlab实现 计算力学报告 平面四节点等参单元 学生姓名： 朱 彬 学号：  问题描述及分析 在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。 根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。 二、有限元划分描述 在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下： 打开ansys，在单元类型中选择solid-Quad 4 node 182单元；建立类半径为外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。如下图： 图1 三、有限元程序及求解 程序求解使用了matlab语言。具体如下： 程序： clc clear E=2e11; %弹性模量 NU=; %泊松比 t=; %厚度 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); %初始总体刚度矩阵 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)=*(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); %形函数分别对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; end D=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; %D阵 k=zeros(8,8); Kss=[,,0,,]; %5*5高斯积分点 Ett=[,,0,,]; H=[,,,,];%高斯积分权系数 for j=1:5 %高斯积分求单元刚度阵 for l=1:5 Ks=Kss(j);Et=Ett(l); B=sub