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

平面四节点等参 单元 初级会计实务单元训练题天津单元检测卷六年级下册数学单元教学设计框架单元教学设计的基本步骤主题单元教学设计 matlab实现 计算力学 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 平面四节点等参单元 学生姓名: 朱 彬 学号: 20100311 一、 问题描述及分析 在无限大平面内有一个小圆孔。孔内有一集中力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 % 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 母单元映射到真实单元 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.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=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; %总体刚度矩阵合成 end KK=K; b=size(w,1); for i=1: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; %加载荷10kN U=K\f; %节点位移 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*ele m(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 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*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,Et},{1,-1}); sigma4(:,i)=double(w4); end 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); sigmar=zeros(3,e); %r相同角度不同节点应力矩阵 sigmat=zeros(3,f); %角度不同r不同节点应力矩阵 msigmar=zeros(1,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)+sigma3(:,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*(si gmar(3,i))^2))^0.5; %求节点处的mises应力 end msigmar %mises应力 for i=1:f %同上 g=0; for j=1:m if elem(j,1)-d(i)==0 g=g+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*((sigmat(1,i)-sigmat(2,i))^2+sigmat(1,i)^2+sigmat(2,i)^2+6*(si gmat(3,i))^2))^0.5; end msigmat 四、计算结果: 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应力程序: clc A=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,18 3,184,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; end for 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; end ansysr ansyst 计算结果如下: ansysr = 1.0e+004 * Columns 1 through 9 0.0388 2.2423 0.0969 0.0439 0.0521 0.0702 0.1100 0.1444 0.2330 Columns 10 through 11 0.4766 0.5914 ansyst = 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.7700 2.等参单元编程计算结果 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.2183 msigmat = 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 节点号 24 29 47 58 78 79 137 Ansys解 0.0388 2.2423 0.0969 0.0439 0.0521 0.0702 0.1100 程序解 0.0365 3.4201 0.0458 0.0415 0.0448 0.0673 0.1396 节点号 149 160 166 186 Ansys解 0.1444 0.2330 0.4766 0.5914 程序解 0.1076 0.1859 0.4403 3.2183 (注:单位为1.0e4Pa) 2沿圆周方向(排列顺序按单元从大到小) 表2 节点号 166 167 168 169 170 171 Ansys解 0.4765 0.2350 0.0983 0.0885 0.0620 0.2032 程序解 0.4403 0.3125 0.2060 0.1686 0.2256 0.6079 节点号 172 173 174 175 176 177 Ansys解 0.0700 0.1019 0.1291 0.1331 0.1252 0.1067 程序解 0.2483 0.2184 0.1692 0.1140 0.2197 0.1035 节点号 178 179 180 181 182 183 Ansys解 0.1078 0.0617 0.5409 0.4426 0.1214 0.1341 程序解 0.1069 0.0854 1.5036 0.4233 0.1903 0.1553 节点号 184 185 Ansys解 0.0569 0.0770 程序解 0.2950 0.2027 (注:单位为1.0e4Pa) 对比分析:在沿半径方向数据吻合的很好，误差基本上是很小的;在圆周方向数据基本吻合，但是在一些应力值很小的点存在较大的差别，这可能与ansys处理四节点单元节点应力与程序处理方式差异有关。