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

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