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

《平面四节点等参单元matlab实现》由会员分享，可在线阅读，更多相关《平面四节点等参单元matlab实现(15页珍藏版)》请在人人文库网上搜索。

1、计算力学报告平面四节点等参单元学生姓名： 朱 彬学号： 一、 问题描述及分析在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：打开ansys，在单元类型中选择solid-Quad 4 node 182单元；建。

2、立类半径为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,。

3、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:msyms 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);endx=0;y=0;for j=1:4 %标准母单元映射到真实单元x=x+N(j)*X(elem(i,j),1);y=y+N(j)*X(ele。

4、m(i,j),2);endJ1=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;endD=(E/(1-NU*NU)*1,NU,0;NU,1,0;0,0,(1-NU)/2; %D阵k=zeros(8,8); Kss=-。

5、0.,-0.,0,0.,0.; %5*5高斯积分点Ett=-0.,-0.,0,0.,0.;H=0.,0.,0.,0.,0.;%高斯积分权系数for j=1:5 %高斯积分求单元刚度阵for l=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)-。

6、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);for i=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=Kf; %节点位移for i=1:m %将每个单元各个节点位移集合u(:,i)=U(2*elem(i,1)-1);U(2*elem(i,1);U(2*elem(i,2)-1)。

7、;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);endfor 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:4N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);endx=0;y=0;for j=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;for j=1:4。

8、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阵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)=doub。

9、le(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相同角。

10、度不同节点应力矩阵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)=0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma1(:,j);endif elem(j,2)-c(i)=0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma2(:,j);endif ele。

11、m(j,3)-c(i)=0g=g+1;sigmar(:,i)=sigmar(:,i)+sigma3(:,j);endif elem(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应力for i=1:f %同上g=0。

12、;for j=1:mif elem(j,1)-d(i)=0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma1(:,j);endif elem(j,2)-d(i)=0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma2(:,j);endif elem(j,3)-d(i)=0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma3(:,j);endif elem(j,4)-d(i)=0g=g+1;sigmat(:,i)=sigmat(:,i)+sigma4(:,j);endendsigmat(:,i)=sigmat(:,i)/g; msig。

13、mat(:,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;其在图中的。

14、位置为图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。

15、=size(d,1);for i=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;endfor i=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 *Columns 1 through 9 0.0388 2.2423 0.0969 0.0439 0.0521 0.0702 0.1100 0.144。

16、4 0.2330Columns 10 through 11 0.4766 0.5914ansyst =1.0e+003 *Columns 1 through 9 4.7655 2.3504 0.9826 0.8454 0.6197 2.0322 0.7001 1.0187 1.2914Columns 10 through 18 1.3306 1.2516 1.0667 1.0775 0.6165 5.4094 4.4263 1.2137 1.3410Columns 19 through 20 0.5689 0.77002.等参单元编程计算结果msigmar =1.0e+004 *Columns。

17、 1 through 9 0.0365 3.4201 0.0458 0.0415 0.0448 0.0673 0.1396 0.1076 0.1859Columns 10 through 11 0.4403 3.2183msigmat =1.0e+004 *Columns 1 through 9 0.4403 0.3125 0.2060 0.1686 0.2256 0.6079 0.2483 0.2184 0.1692Columns 10 through 18 0.1140 0.2197 0.1035 0.1069 0.0854 1.5036 0.4233 0.1903 0.1553Colum。

18、ns 19 through 20 0.2950 0.2027五、讨论比较对所选取的节点的mises应力列表做对比：1.沿半径方向(排列顺序为节点号从小到大)表1节点号242947587879137Ansys解0.03882.24230.09690.04390.05210.07020.1100程序解0.03653.42010.04580.04150.04480.06730.1396节点号149160166186Ansys解0.14440.23300.47660.5914程序解0.10760.18590.44033.2183(注：单位为1.0e4Pa)2沿圆周方向(排列顺序按单元从大到小)表2节点。

19、号166167168169170171Ansys解0.47650.23500.09830.08850.06200.2032程序解0.44030.31250.20600.16860.22560.6079节点号172173174175176177Ansys解0.07000.10190.12910.13310.12520.1067程序解0.24830.21840.16920.11400.21970.1035节点号178179180181182183Ansys解0.10780.06170.54090.44260.12140.1341程序解0.10690.08541.50360.42330.19030.1553节点号184185Ansys解0.05690.0770程序解0.29500.2027(注：单位为1.0e4Pa)对比分析：在沿半径方向数据吻合的很好，误差基本上是很小的；在圆周方向数据基本吻合，但是在一些应力值很小的点存在较大的差别，这可能与ansys处理四节点单元节点应力与程序处理方式差异有关。

相关资源：欧德克连杆仿真设计软件Linkage_linkage软件-其它工具类资源-CSDN…

来源：二池一