有限元程序设计大作业

1、实用文档1 .不同板宽的孔边应力集中问题姓名：胡宇学号：21201201282 .摘要本文采用MATLAB和FOTRAN四节点平面单元，利用有限元数值解法对不同板宽的孔边应力集中问题进行了数值模拟研究。对于不同的板宽系数（半板宽b/孔半径r）,得到了不同的应力集中系数并且与解析解进行了比较，验证了有限元解的正确性，并且得出了解析解的适用范围。3 .引言通常情况下的有限元分析过程是运用可视化分析软件（如ANSYS、ABAQUS、SAP等）进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。 随着有限元技术的逐

2、步成熟，它被应用在越来越复杂的问题处理中。MATLAB 是由美国 MATHWORKS 公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学 数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易 于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的 众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交 互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。4 . MATLAB部分1 ，计算模型本程序采用MATLAB 编程，编制平面四边形四节点等参元程序，用以求解近似平面结构

3、问题。本程序的研究对象为中央开有小孔的长方形板，选取的材料参数为：板厚h=1、材料强度 E=1.0e11 Pa、泊松比mu=0.3。止匕外，为方便网格的划分和计算，本文所取板的长度与宽度相等。其孔半径为r=1 ，板宽为2b 待定。由于本程序的目的在于验证有限元解的正确性和确定解析解的适用范围，因此要求网格足够细密，以满足程序的精度要求。同时为了减小计算量，我采取网格径向长度递增的网格划分方法。此种方法特点是，靠近小孔部分的网格细密，在远离小孔的过程中，网格逐渐变得稀疏。以下为网格节点坐标计算和单元节点排序的MATLAB 源程序：dr = (b-r)/(m+mA2/2+mA3/6);dfi= (

4、pi/2)/n ;gNode = zeros( (n+1)*(m+1), 2 ) ;for i=1:1:mfor j=1:1:n+1gNode( (i-1)*(n+1)+j, : ) = cos(dfi*(j-1)*(r+.dr*(i-1)+(i-1)A2/2+(i-1)A3/6),sin(dfi*(j-1)*(r+dr*(i-1)+(i-1)A2/2+(i-1)A3/6) ;endendfor i=1:1:(n/2+1)gNode(n+1)*m+i, : )=b,b*tan(dfi*(i-1);endfor i=1:1:n/2gNode(n+1)*m+(n/2+1)+i, : )=b/tan

5、(dfi*(i+n/2),b;endgElement =zeros( m*n, 5 ) ;for i=1:mfor j=1:ngElement( (i-1)*n+j, 1:4) = (i-1)*(n+1)+j,.i*(n+1)+j,.i*(n+1)+j+1,.(i-1)*(n+1)+j+1 ;EndendgElement( :, 5 ) = 1 ;以上源程序中 gNode为节点坐标矩阵，gElement为单元节点排列矩阵Deformed shape024681。图（2）假定该板所受的外载为p=1.0e9 Pa。由于是取1/4板进行计算，需要在分界面上加上合适的边界条件。以小孔圆心为原点，板长方

6、向为x轴，板宽x轴平行的分界面上的节点的y向固定，x向方向为y轴。我所取得边界条件为：与标准文案实用文档可以移动；与 y轴平行的分界面上的节点的x向固定，y向可以移动。划分好网格，约束和载荷加好后。就可以计算了。取b=10, m=5, n=10时的应力图。图（3）为x向应力分布图，图（4）为y向应力分布图，图（5） 为剪切力分布图。2.521.510.50图(3)图（4）Stress *xy distribution* 1(/0246810图(5)2,孔边应力研究2 11-72无限大板宽的孔边应力集中问题的弹性力学的解析解为:p_八. r . 一* cos2 11 1 - 32<P2 人

7、41 3匚。42rp.-上 cos2-：22,：i - -卫 sin 2口221er21 3-7：2标准文案在孔边的y轴上仃日有分布:71rl _ 900=P-I 2 P23d2 ：4一- 900, 二 a)= 3p为了验证有限元解的正确性，我取10x20的网格进行了计算，得到了轴上的x向应力分布曲线，并且与解析解进行了比较。如图（6）。y轴上的x向应力解析解为:221 r0)= p 1 +一I 2V4 、+ 3d2y4Jxx = 0, y = a = 3p乂=0轴上的k向应力变化曲线琏力大小，PaxlO5EE慨盘sw图宅小割些图（6）同时，小孔边缘的剪切力分布如图（7）所示。剪切力的解析解为

8、1.x p sin 4 1-sin 公xy 2卷 A.kK-R&孔边缘翦切力变化曲线 逆时针）0.40.6口,8112141.6角度（:单位为弧度）图（7）图中，红色*点为有限元节点解，黑色 x点为解析解点。由图可知，有限元解 和解析解基本相同，有限元解是有效的。同时，为了研究不同板宽对应力集中系数的影响，我选取了多组数据进行对比。如表（1）所示表（1）b or ( )b-rm1m2K22116.028646.11143213.626464.03154323.228973.52485422.909483.31696522.712293.20987622.5751103.14708732

9、.8527113.10699832.7762113.017010932.7102123.0604171642.6707163.0071262552.6589202.9911373662.6545242.9851504972.6526282.9824表中，m1,m2为径向单元的数目。尸2分别为与之对应的应力集中系数（Knmax/p）。之为半板宽b与孔半径r的比值，即t=b/r。解析解中氐 的取值为3.下图（7）为解析解和有限元解的比较。应力集卬系骸和极宽的关系B -寺5 5 5 4 5 S.4.3.蕊然frM-R囱图(7)图中黑色o 点为理论应力集中系数，红色x 点为较低网格密度下的应力集中系数

10、变化曲线，蓝色 * 点为较高网格密度下的应力集中系数变化曲线。由图可知，同为有限元解，网格密度越高，精确度就越高。同时，我们还可以知道，解析解并不是适用于所有的情况。当板宽与孔径的比值小于5 时，解析解与有限元解的差别变大，即解析解不再适用。5. FORTRA解分为了进一步验证结果的正确性，我同时运用FORTRAN 对相同的网格进行了计算。计算结果表明，有限元法所得的结果是可信的。6. 结论有限元解和解析解各有优劣，有限元解的优点是基本上对任意情况都适用，但是一般情况下计算量都很大，结果获取速度慢。解析解则求解迅速，但并不适用于所有情况。本文所研究的孔边应力集中问题就很好的体现了两种方法的特点

11、。解析解并不适用于半板宽和孔半径之比小于5 的情况，但是对于板宽远大于孔径的情况具有很高的精确度。而有限元解的精确度依赖于网格的密度和算法，在算法相同的情况下，为增加精度而加密网格会使计算量增加，在精度要求不高的情况下适用。MATLAB 源程序实用文档function huyu2120120128global r b m n p E mu% 定义板的尺寸r = 1 ;% 小孔半径b = 10 ;% 半板宽format short e% 定义网格密度m = 1;%径向划分的数目n=2*m 为最优% 定义有限元模型% 求解有限元模型把计算结果显示出来n = 2;%周向划分的数目（必须为偶数）% 定

12、义材料性质E = 1e11 ;%弹性模量mu =0.3;% 泊松比% 定义均布荷载大小p = 1.0e9;%单位Pa%FemModel ;SolveModel ;DisplayResults ;%returnfunction FemModel% 定义有限元模型% 全局变量及名称%r 圆孔半径%b 板宽% m 半径方向单元数目%n 圆周方向单元数目%E 弹性模量%mu 泊松比% p 均布载荷大小%gNode 节点坐标%gElement -单元定义%gMaterial -材料性质%gBC1 第一类约束条件% gBC3 斜约束%gDF 分布力%返回值：%无global gNode gElement

13、gMaterial gBC1 gDFglobal r b m n p E mu% 计算结点坐标dr = (b-r)/(m+mA2/2+mA3/6);dfi= (pi/2)/n ;gNode = zeros( (n+1)*(m+1), 2 ) ;for i=1:1:mfor j=1:1:n+1gNode( (i-1)*(n+1)+j, : ) = cos(dfi*(j-1)*(r+.dr*(i-1)+(i-1)A2/2+(i-1)A3/6),sin(dfi*(j-1)*(r+dr*(i-1)+(i-1)A2/2+(i-1)A3/6); endendfor i=1:1:(n/2+1)gNode(n

14、+1)*m+i, : )=b,b*tan(dfi*(i-1);endfor i=1:1:n/2gNode(n+1)*m+(n/2+1)+i, : )=b/tan(dfi*(i+n/2),b; end% 单元定义gElement =zeros( m*n, 5 ) ;for i=1:mfor j=1:ngElement( (i-1)*n+j, 1:4) = i*(n+1)+j,.i*(n+1)+j+1,.(i-1)*(n+1)+j+1,.(i-1)*(n+1)+j ;end endgElement( :, 5 ) = 1 ;% 定义材料性质gMaterial = E, mu, 0 ;% 定义约束条

15、件gBC1 = zeros( (m+1)*2, 3 ) ;for i=1:1:m+1gBC1(i, 1: 2)=(n+1)*(i-1)+1,(n+1)*(i-1)+1)*2;endfor i=1:1:m+1gBC1(m+1+i,1 :2 )=(n+1)*i,(n+1)*i-1)*2+1;endgBC1( :, 3 )=0;% 定义分布压力 gDF = zeros(n/2,4);标准文案实用文档for i=1:1:n/2gDF(i,1:2)=(n+1)*m+i,(n+1)*m+i+1;endgDF( :, 3 )=p;gDF( :, 4 )=p;returnfunction SolveModel

16、% 求解有限元模型% 说明：%该函数求解有限元模型，过程如下%1.计算单元刚度矩阵，集成整体刚度矩阵%2.计算单元的等效节点力，集成整体节点力向量%3.处理约束条件，修改整体刚度矩阵和节点力向量%4.求解方程组，得到整体节点位移向量gMaterialglobal gNode gElement gBC1 gDF gK gDelta gElementStress% step1. 定义整体刚度矩阵和节点力向量node_number,dummy = size( gNode ) ;gK = sparse( node_number * 2, node_number * 2 ) ;f = sparse( n

17、ode_number * 2, 1 ) ;% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中 element_number,dummy = size( gElement ) ;for ie=1:1:element_numberfprintf( sprintf( '计算单元刚度矩阵并集成，当前单元: %dn', ie ) ) ;k = StiffnessMatrix( ie ) ;AssembleStiffnessMatrix( ie, k ) ;end% step3. 计算分布压力的等效节点力，并集成到整体节点力向量中df_number,dummy = size( gD

18、F ) ;for idf = 1:1:df_numberenf = EquivalentDistPressure( gDF(idf,1), gDF(idf,2),gDF(idf,3),gDF(idf,4) ) ;i = gDF(idf, 1) ;j = gDF(idf, 2) ;f( (i-1)*2+1:(i-1)*2+2, 1 ) = f( (i-1)*2+1:(i-1)*2+2, 1 ) + enf( 1:2 ) ;f( (j-1)*2+1:(j-1)*2+2, 1 ) = f( (j-1)*2+1:(j-1)*2+2, 1 ) + enf( 3:4 ) ;end% step4. 处理约束

19、条件，修改刚度矩阵和节点力向量。采用乘大数法bc1_number,dummy = size( gBC1 ) ;for ibc=1:1:bc1_numberg = gBC1(ibc, 2 ) ;f(g) = gBC1(ibc, 3)* gK(g,g) * 1e15 ;gK(g,g) = gK(g,g) * 1e15 ;end% step5. 求解方程组，得到节点位移向量gDelta = gK f ;% step6. 计算每个单元的结点应力gElementStress = zeros( element_number,5, 3 ) ;delta = zeros( 8, 1 ) ;for ie = 1

20、:element_numberxi = -111-1 ;eta = -1-111 ;for n=1:4B = MatrixB( ie, xi(n), eta(n) ) ;D = MatrixD( ie, 1 ) ;delta(1:2:7) = gDelta( gElement(ie,1:4)*2-1 ) ;delta(2:2:8) = gDelta( gElement(ie,1:4)*2) ;sigma = D*B*delta;gElementStress( ie, n, :) = sigma ;endendreturnfunction k = StiffnessMatrix( ie )% 计

21、算平面应力等参数单元的刚度矩阵% 输入参数：%ie - 单元号% 返回值：%k - 单元刚度矩阵% 说明：% 用高斯积分法求解平面等参数单元的刚度矩阵k = zeros( 8, 8 ) ;D = MatrixD( ie,1) ;% 3 x 3 高斯积分点和权系数%x = -0.774596669241483,0.0,0.774596669241483 ;%w = 0.555555555555556,0.888888888888889,0.555555555555556 ;% 2 x 2 高斯积分点和权系数x = -0.577350269189626, 0.577350269189626 ;w

22、= 1, 1 ;for i=1:1:length(x)for j=1:1:length(x)B = MatrixB( ie, x(i), x(j) ) ;J = Jacobi( ie, x(i), x(j) ) ;k = k + w(i)*w(j)*transpose(B)*D*B*det(J) ;endendreturnfunction D = MatrixD( ie, opt )% 计算单元的弹性矩阵D% 输入参数：% ie 单元号%opt 问题选项%1 => 平面应力%2 => 平面应变% 返回值：% D 弹性矩阵Dglobal gElement gMaterialE = g

23、Material( gElement(ie, 5), 1 ) ;% 弹性模量mu = gMaterial( gElement(ie, 5), 2 ) ;% 泊松比if opt = 1% 平面应力的弹性常数A1 = mu ;A2 = (1-mu)/2 ;A3= E/(1-muA2);else% 平面应变的弹性常数A1 = mu/(1-mu) ;A2 = (1-2*mu)/2/(1-mu) ;A3 = E*(1-mu)/(1+mu)/(1-2*mu) ;endD = A3* 1 A1 0;A110;00 A2 ;returnfunction B = MatrixB( ie, xi, eta )%

24、计算单元的应变矩阵B% 输入参数：%ie 单元号%xi,eta 局部坐标% 返回值：% B 在局部坐标处的应变矩阵BN_x,N_y = N_xy( ie, xi, eta );B = zeros( 3, 8 ) ;for i=1:1:4B(1:3,(2*i-1):2*i) = N_x(i)0;0N_y(i);N_y(i) N_x(i);endreturnfunction N_x, N_y = N_xy( ie, xi, eta )% 计算形函数对整体坐标的导数% 输入参数：%ie 单元号%xi,eta 局部坐标% 返回值：%N_x 在局部坐标处的形函数对x 坐标的导数%N_y 在局部坐标处的形

25、函数对y 坐标的导数J = Jacobi( ie, xi, eta ) ;N_xi,N_eta = N_xieta( ie, xi, eta ) ;A=JN_xi;N_eta ;N_x = A(1,:) ;N_y = A(2,:) ;returnfunction N_xi, N_eta = N_xieta( , xi, eta )% 计算形函数对局部坐标的导数% 输入参数：%ie 单元号%xi,eta 局部坐标% 返回值：% N_xi 在局部坐标处的形函数对xi 坐标的导数% N_eta 在局部 坐标处的形函数对eta 坐标的导数x = -1, 1, 1, -1 ;e = -1, -1, 1,

26、 1 ;N_xi = zeros( 1, 4 ) ;N_eta = zeros( 1, 4 ) ;N_xi(1)= x(1)*(1+e(1)*eta)/4;N_eta(1) = e(1)*(1+x(1)*xi)/4;N_xi(2)= x(2)*(1+e(2)*eta)/4;N_eta(2) = e(2)*(1+x(2)*xi)/4;N_xi(3)= x(3)*(1+e(3)*eta)/4;N_eta(3) = e(3)*(1+x(3)*xi)/4;N_xi(4)= x(4)*(1+e(4)*eta)/4;N_eta(4) = e(4)*(1+x(4)*xi)/4;returnfunction

27、J = Jacobi( ie, xi, eta )% 计算雅克比矩阵% 输入参数：%ie 单元号%xi,eta 局部坐标% 返回值：% J 在局部坐标(xi,eta) 处的雅克比矩阵global gNode gElementx = gNode(gElement(ie,1:4),1) ;y = gNode(gElement(ie,1:4),2) ; N_xi,N_eta = N_xieta( ie, xi, eta ) ;x_xi = N_xi * x ;x_eta = N_eta * x ;y_xi = N_xi * y ;y_eta = N_eta * y ;J = x_xi, y_xi;

28、x_eta, y_eta ;returnfunction AssembleStiffnessMatrix( ie, k )% 把单元刚度矩阵集成到整体刚度矩阵% 输入参数:% ie - 单元号% k - 单元刚度矩阵% 返回值 :% 无global gElement gKfor i=1:1:4for j=1:1:4for p=1:1:2for q=1:1:2s = (i-1)*2+p ;t = (j-1)*2+q ;S = (gElement(ie,i)-1)*2+p ;T = (gElement(ie,j)-1)*2+q ; gK(S,T) = gK(S,T) + k(s,t) ; ende

29、ndendendreturnfunction enf = EquivalentDistPressure( i, j, pi, pj )% 计算线性分布压力的等效节点力% 输入参数:% i,j 结点号% pi,pj 跟结点号对应的压力值% 返回值 :% enf 等效 节点力向量global gNode% 获取结点坐标xi = gNode( i, 1 ) ;xj = gNode( j, 1 ) ;yi = gNode( i, 2 ) ;yj = gNode( j, 2 ) ;sig(1:4,1)=0;% 2 x 2 高斯积分点和权系数x = -0.577350269189626, 0.577350

30、269189626 ;w = 1, 1 ;for i=1:1:length(x)N1=(1-x(i)/2;N2=(1+x(i)/2;dN1=-0.5;dN2=0.5;dx=xi*dN1+xj*dN2;dy=yi*dN1+yj*dN2;p=pi*N1+pj*N2;sig(1)=sig(1)+w(i)*N1*dy*p;sig(2)=sig(2)-w(i)*N1*dx*p;sig(3)=sig(3)+w(i)*N2*dy*p;sig(4)=sig(4)-w(i)*N2*dx*p;endenf = sig ;returnfunction DisplayResults% 显示计算结果，并与解析解结果比较

31、global gNode gElement gMaterial gDelta gElementStress outfile gBC1global gDF n m r p b Stressoutfile=fopen('huyu2120120128.txt','w');Stress = zeros(m+1)*(n+1),3);%集成节点应力矩阵for k=1:1:3for i=1:1:mStress (n+1)*(i-1)+1,k)=gElementStress(n*(i-1)+1,4,k);endfor j=1:1:nStress (n+1)*m+j,k)=gEl

32、ementStress(n*(m-1)+j,1,k);endStress (n+1)*(m+1),k)=gElementStress(n*m,2,k);for i=1:1:mfor j=1:1:nStress (i-1)*(n+1)+j+1,k)=gElementStress(i-1)*n+j,3,k);endendendxx_max=max(max(Stress(:,1);N_xx_max=find(Stress=xx_max);xx_min=min(min(Stress(:,1);N_xx_min=find(Stress=xx_min);yy_max=max(max(Stress(:,2)

33、;N_yy_max=find(Stress=yy_max)-(m+1)*(n+1);yy_min=min(min(Stress(:,2);N_yy_min=find(Stress=yy_min)-(m+1)*(n+1);xy_max=max(max(Stress(:,3);N_xy_max=find(Stress=xy_max)-2*(m+1)*(n+1);xy_min=min(min(Stress(:,3);N_xy_min=find(Stress=xy_min)-2*(m+1)*(n+1);fprintf(outfile, '表 0 节点坐标(mm) n' ) ;fprin

34、tf(outfile, '=n' ) ;fprintf(outfile, ' 结点号X 坐标Y 坐标 n' ) ;fprintf(outfile, 'n' ) ;for i=1:size(gNode(:,1)fprintf(outfile, '%4d , %8.4f , %8.4f n',i,gNode(i,1),gNode(i,2) ; end fprintf( outfile,'=nnn' ) ;fprintf(outfile, '表 1节点位移、应力n' ) ;fprintf(outfile,

35、标准文案实用文档=n');fprintf(outfile,'结点号 x-位移y-位移Stress-xx Stress-yy Stress-xyn'); fprintf(outfile, ' n');for i=1:size(gNode(:,1)fprintf(outfile,'%4d %12.4e %12.4e%12.4e%12.4e%12.4en',i, full(gDelta(2*i-1),full(gDelta(2*i),Stress(i,1),Stress(i,2),Stress(i,3);end fprintf( outfile

36、,'= =nnn' ) ;fprintf(outfile,'表 2单元节点序列n');fprintf(outfile, '=n' ) ;fprintf(outfile,' 单元号 材料组 节点号(逆时针排列)n');fprintf(outfile, '1234n' ) ;fprintf(outfile, 'n' ) ;for i=1:size(gElement(:,1)fprintf(outfile, '%4d , %4d , %4d , %4d , %4d , %4dn',i, g

37、Element(i,5), gElement(i,1:4);endfprintf( outfile,'=nnn');fprintf(outfile,'表 3材料性质n');fprintf(outfile, '=n' ) ;fprintf(outfile,' 材料号 弹性模量 泊松比比重 n');fprintf(outfile, 'n' ) ;for i=1:size(gMaterial(:,1)fprintf(outfile, '%4d%8.3e%8.3e%8.3en',i,gMaterial(i,

38、1),gMaterial(i,2),gMaterial(i,3);endfprintf( outfile,'=nnn');fprintf(outfile,'表 4边界条件n');fprintf(outfile, '=n' ) ;fprintf(outfile,' 序号 结点号 约束自由度 给定值 n');fprintf(outfile, 'n' ) ;for i=1:size(gBC1(:,1)fprintf(outfile, '%4d%4d%4d%4d n',.i,gBC1(i,1),gBC1(i

39、,2),gBC1(i,3) ;endfprintf(outfile, '=nnn');fprintf(outfile,'表 5 分布力n');fprintf(outfile, '=n');fprintf(outfile,' 边界 结点1结点2 均布力1均布力 2n');fprintf(outfile, 'n');for i=1:size(gDF(:,1)fprintf(outfile, '%4d %4d%4d%8.4e%8.4en',i,gDF(i,1),gDF(i,2),gDF(i,3),gDF(

40、i,4);endfprintf( outfile,'=nnn');fprintf( outfile,'表 6 xx 向节点应力极值n');fprintf(outfile, '=n');fprintf(outfile,' 极大点 Stress_xx_max 极小点 Stress_xx_min n');fprintf(outfile, ' n');fprintf(outfile,'%4d%12.4e%4d%12.4en', .N_xx_max,xx_max,N_xx_min,xx_min);fprint

41、f(outfile, '=nnn') ;fprintf( outfile,'表 7 yy 向节点应力 极值n');fprintf(outfile, '=n' ) ;fprintf(outfile,' 极大点 Stress_yy_max 极小点 Stress_yy_min n');fprintf(outfile, ' n');fprintf(outfile,'%4d%12.4e%4d%12.4en', .N_yy_max,yy_max,N_yy_min,yy_min);fprintf(outfile,

42、 '=nnn') ;fprintf( outfile,'表 8 xy 向节点应力极值n');fprintf(outfile, '=n' ) ;fprintf(outfile,' 极大点 Stress_xy_max 极小点 Stress_xy_min n');fprintf(outfile, ' n');fprintf(outfile,'%4d%12.4e%4d%12.4en',.N_xy_max,xy_max,N_xy_min,xy_min) ;fprintf(outfile, '=nnn&

43、#39;) ;fclose(outfile);% (7)画图ord=1,2,3,4,1;for e=1:size(gElement(:,1)标准文案实用文档for i=1:5xpt1(i)=gNode(gElement(e,ord(i),1);ypt1(i)=gNode(gElement(e,ord(i),2);disx1(i)=gDelta(2*gElement(e,ord(i)-1);disy1(i)=gDelta(2*gElement(e,ord(i);spt11(i)=gElementStress(e,ord(i),1);spt21(i)=gElementStress(e,ord(i)

44、,2);spt31(i)=gElementStress(e,ord(i),3);endxpt(1)=xpt1(1);ypt(1)=ypt1(1);disx(1)=disx1(1);disy(1)=disy1(1);spt1(1)=spt11(1);spt2(1)=spt21(1);spt3(1)=spt31(1);xpt(2)=xpt1(2);ypt(2)=ypt1(2);disx(2)=disx1(2);disy(2)=disy1(2);spt1(2)=spt11(2);spt2(2)=spt21(2);spt3(2)=spt31(2);xpt(3)=xpt1(3);ypt(3)=ypt1(

45、3);disx(3)=disx1(3);disy(3)=disy1(3);spt1(3)=spt11(3);spt2(3)=spt21(3);spt3(3)=spt31(3);xpt(4)=xpt1(4);ypt(4)=ypt1(4);disx(4)=disx1(4);disy(4)=disy1(4);spt1(4)=spt11(4);spt2(4)=spt21(4);spt3(4)=spt31(4);xpt(5)=xpt1(1);ypt(5)=ypt1(1);disx(5)=disx1(1);disy(5)=disy1(1);spt1(5)=spt11(1);spt2(5)=spt21(1)

46、;spt3(5)=spt31(1);scale=1e1;% 画变形图subplot(2,3,1)hold onplot(xpt,ypt,'b-')plot(xpt+scale*disx,ypt+scale*disy,'r-') title('Deformed shape')axis equal% axis off% 画应力图subplot(2,3,2)hold onfill(xpt+scale*disx,ypt+scale*disy,spt1) colorbaraxis equalaxis (0 b 0 b)title('Stress-x

47、x distribution')% axis offsubplot(2,3,3)hold onfill(xpt+scale*disx,ypt+scale*disy,spt2) colorbaraxis equalaxis (0 b 0 b)title('Stress-yy distribution')% axis offsubplot(2,3,4)hold onfill(xpt+scale*disx,ypt+scale*disy,spt3)colorbaraxis equalaxis (0 b 0 b)title('Stress-xy distribution&

48、#39;)% axis offend% 画x=0 轴上的 x 向应力变化曲线subplot(2,3,5)hold onA=zeros(1,m+1);B=zeros(1,m+1);C=zeros(1,m+1);D=zeros(1,m+1);for i=1:1:m+1A(i) = gNode(i*(n+1),2);B(i) = Stress(i*(n+1),1);C(i) = p*(1+rA2/2/gNode(i*(n+1),2)A2+3*”4/2/gNode(i*(n+1),2)A4);D(i) = p;endplot(B,A,'r:*');plot(C,A,'k:x&#

49、39;);plot(D,A,'b-');title('x=0 轴上的 x 向应力变化曲线')ylabel(' 点距小孔圆心的距离 (mm)')xlabel(' 应力大小/ Pa')axis fill% 画小孔边缘的剪切力力变化曲线(逆时针)subplot(2,3,6)hold onA=zeros(1,n+1);B=zeros(1,n+1);C=zeros(1,n+1);D=zeros(1,n+1);for i=1:1:n+1A(1,i)=(i-1)*(pi/2/n);B(1,i)=Stress(i,3);C(1,i)=p*(sin

50、(4*A(1,i)-sin(2*A(1,i)/2;D(1,i)=(xy_max+xy_min)/2;endplot(A,B,'r:*');plot(A,C,'k:x');plot(A,D,'b-');title(' 孔边缘剪切力变化曲线(逆时针 )')xlabel(' 角度(单位为弧度) ')ylabel(' 应力大小/ Pa')axis fillreturn10x20 网格计算结果结点号X 坐标Y 坐标x-位移y-位移Stress-xxStress-yyStress-xy11.00000.0000

51、3.0494e-02-6.7903e-19-6.7755e+07-1.0550e+09-4.2692e+0720.99690.07853.0399e-02-8.1180e-04-8.3977e+07-1.0451e+091.0140e+0830.98770.15643.0115e-02-1.6189e-03-1.0225e+08-9.7927e+081.4996e+0840.97240.23342.9645e-02-2.4168e-03-1.1795e+08-8.6325e+081.7814e+0850.95110.30902.8990e-02-3.2006e-03-1.2395e+08-7.0667e+081.