有限元程序的设计大作业

1、1. 不同板宽的孔边应力集中问题姓名：胡宇学号： 21201201282. 摘要本文采用 MATLAB 和 FOTRAN 四节点平面单元， 利用有限元数值解法 对不同板宽的孔边应力集中问题进行了数值模拟研究。对于不同的板宽系数 （半板宽 b/ 孔半径 r ），得到了不同的应力集中系数1 ，并且与解析解进行了比较，验证了有限元解的正确性，并且得出了解析解的适用范围。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 待 定。由于本程序的目的在于验证有限元解的正确性和确定解析解的适用范 围，因此要求网格足够细密，以满足程序的精度要求。同时为了减小计算量， 我采取网格径向长度递增的网格划分方法。此种方法特点是，靠近小孔部分 的网格细密，在远离小孔的过程中，网格逐渐变得稀疏。 以下为网格节点坐标计算和单元节点排序的

4、MATLAB 源程序： dr = (b-r)/(m+m2/2+m3/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)2/2+(i-1)3/6),sin(dfi*(j-1)*(r+dr*(i-1)+(i-1)2/2+(i-1)3/6) ;endendfor i=1:1:(n/2+1)gNode(n+1)*m+i, : )=b,b*tan(dfi*(i-1);end for i

5、=1:1:n/2 gNode(n+1)*m+(n/2+1)+i, : )=b/tan(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 为单元节点排列矩阵 图（ 1 ）为 5x10 的网格划分图，图（ 2 ）为结构的微变形图。图（2）

6、假定该板所受的外载为 p=1.0e9 Pa 。由于是取 1/4 板进行计算，需要 在分界面上加上合适的边界条件。以小孔圆心为原点，板长方向为 x 轴，板 宽方向为 y 轴。我所取得边界条件为：与 x 轴平行的分界面上的节点的 y 向固定， x 向 可以移动；与 y 轴平行的分界面上的节点的 x 向固定， y 向可以移动。划分好网格， 约束和载荷加好后。 就可以计算了。 取 b=10 ，m=5 ，n=10 时的应力图。图（ 3）为 x 向应力分布图，图（ 4 ）为 y 向应力分布图，图（ 5） 为剪切力分布图。图（ 3 ）图（ 4 ）图（ 5 ）2 ，孔边应力研究无限大板宽的孔边应力集中问题的弹

7、性力学的解析解为：pcos22r 2 r 2r2 1 3r 22pcos243r44psin22r22r22在孔边的 y 轴上 有分布 :900 p 1r2223 r424900,a 3p的网格进行了计算， 得到了 y 如图（ 6 ）。为了验证有限元解的正确性， 我取 10x20轴上的 x 向应力分布曲线，并且与解析解进行了比较y 轴上的 x 向应力解析解为：x 0,y a 3p1 r2 3r4x x 0 p 1 2 4 ,x2 y2 2 y4图（6）同时，小孔边缘的剪切力分布如图（1 xy 27 ）所示。剪切力的解析解为p sin 4 sin 2图（ 7 ）图中，红色 * 点为有限元节点解，

8、黑色 x 点为解析解点。由图可知，有限元解 和解析解基本相同，有限元解是有效的。同时，为了研究不同板宽对应力集中系数的影响，我选取了多组数据进 行对比。如表（ 1 ）所示。表（ 1 ）b or （ ）b-rm11m222116.028646.11143213.626464.03154323.228973.52485422.909483.31696522.712293.20987622.5751103.14708732.8527113.10699832.7762113.017010932.7102123.0604171642.6707163.0071262552.6589202.99113736

9、62.6545242.9851504972.6526282.982412分别为与之对应的应力集中系数表中， m1 , m2为径向单元的数目。max / p）。 为半板宽 b 与孔半径 r 的比值， 即 b/r 。解析解中的取值为 3. 下图（ 7 ）为解析解和有限元解的比较。图（ 7 ）图中黑色 o 点为理论应力集中系数，红色 x 点为较低网格密度下的应力集中系数变化曲线， 蓝色 * 点为较高网格密度下的应力集中系数变化曲线。由图可知，同为有限元解，网格密度越高，精确度就越高。同时，我们还可以知道， 解析解并不是适用于所有的情况。当板宽与孔径的比值小于 5 时，解析解与 有限元解的差别变大，即

10、解析解不再适用。5. FORTRAN 部分为了进一步验证结果的正确性， 我同时运用 FORTRAN 对相同的网格进行了计算。计算结果表明，有限元法所得的结果是可信的6. 结论有限元解和解析解各有优劣，有限元解的优点是基本上对任意情况都适 用，但是一般情况下计算量都很大，结果获取速度慢。解析解则求解迅速， 但并不适用于所有情况。本文所研究的孔边应力集中问题就很好的体现了两种方法的特点。解析 解并不适用于半板宽和孔半径之比小于 5 的情况，但是对于板宽远大于孔径 的情况具有很高的精确度。而有限元解的精确度依赖于网格的密度和算法，在算法相同的情况下，为增加精度而加密网格会使计算量增加，在精度要求不高

11、的情况下适用。附录一 MATLAB 源程序function huyu2120120128global r b m n p E mu % 定义板的尺寸r = 1 ;% 小孔半径b = 10 ;% 半板宽format short e % 定义网格密度m = 1;% 径向划分 的数目n = 2; % 周向划分 的数目（必须为偶数） n=2*m 为最优 % 定义材料性质E = 1e11 ; % 弹性模 量mu =0.3; % 泊松比% 定义均布荷载大小p = 1.0e9;% 单位 PaFemModel ;% 定义有 限元模型SolveModel ;% 求解有限 元模型DisplayResults ;%

12、 把计算结果 显示出来returnfunction FemModel%定义有限元模型%全局变量及名称%r 圆孔半径%b 板宽%m 半径 方向单元数目%n 圆周 方向单元数目%E 弹性 模量%mu 泊松 比%p 均布 载荷大小%gNode 节点坐标%gElement -单元 定义%gMaterial -材料性 质%gBC1 第一类约束条件%gBC3 斜约束%gDF 分布力%返回值：%无global gNode gElement gMaterial gBC1 gDF globalr b m n p E mu% 计算结点坐标dr = (b-r)/(m+m2/2+m3/6) ;dfi= (pi/2)/

13、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)2/2+(i-1)3/6),sin(dfi*(j-1)*(r+dr*(i-1)+(i-1)2/2+(i-1)3/6) ; endendfor i=1:1:(n/2+1)gNode(n+1)*m+i, : )=b,b*tan(dfi*(i-1);endfor i=1:1:n/2 gNode(n+1)*m+(n/2+1)+i, : )=b/tan(dfi*(i+

14、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 ; endend gElement( :, 5 ) = 1 ;% 定义材料性质 gMaterial = E, mu, 0 ;% 定义约束条件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)*

15、2;endfor i=1:1:m+1 gBC1(m+1+i,1 :2 )=(n+1)*i,(n+1)*i-1)*2+1;end gBC1( :, 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% 求解有限元模型% 说明：% 该函数求解有限元 模型，过程如下%1.计算单元刚度矩阵，集成整体刚度矩阵%2.计算单元的等效节点力，集成整体节点力向量%3.处理约束条件，修改整体刚度矩

16、阵和节点力向量%4.求解方程组，得到整体节点位移向量 gMaterialglobal gNode gElement gBC1 gDF gKgDelta gElementStress% step1. 定义整体刚度矩阵 和节点力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number * 2, node_number * 2 ) ; f = sparse( node_number * 2, 1 ) ;% step2. 计算单元刚度矩阵 ，并集成到整体刚度矩阵中 element_number,dummy = size( gElem

17、ent ) ;for ie=1:1:element_numberfprintf( sprintf( 计算单元刚度矩阵并集成，当前单元: %dn, iek = StiffnessMatrix( ie ) ;AssembleStiffnessMatrix( ie, k ) ;end% step3. 计算分布压力的等 效节点力，并集成到整体节点力向量中 df_number,dummy = size( gDF ) ;for idf = 1:1:df_numberenf = EquivalentDistPressure( gDF(idf,1), gDF(idf,2),gDF(idf,3),gDF(idf

18、,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. 处理约束条件，修 改刚度矩阵和节点力向量。采用乘大数法 bc1_number,dummy = size( gBC1 ) ;for ibc=1:1:bc1_numberg = gBC1(ibc, 2 )

19、;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:element_numberxi = -1 1 1 -1 ;eta = -1 -1 1 1 ;for n=1:4B = MatrixB( ie, xi(n), eta(n) )

20、 ;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 )% 计算平面应力等参数单元 的刚度矩阵% 输入参数：%ie - 单元号% 返回值：%k - 单元刚度矩阵% 说明：等参数单元的刚度矩阵%用高斯积分法求解平面k = zer

21、os( 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 = 1, 1 ;for i=1:1:length(x)for j=1:1:length(x)B = MatrixB( ie, x(i), x(j) ) ;J =

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

23、应 力的弹性常数A1 = mu ;A2 = (1-mu)/2 ;A3 = E/(1-mu2) ;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;A1 1 0;0 0 A2returnfunction B = MatrixB( ie, xi, eta )% 计算单元的应变矩阵 B% 输入参数：% ie 单元号%xi,eta 局部坐标% 返回值：%B 在局部 坐标处的应变矩阵N_x,N_y = N_xy( ie, xi, eta );B

24、 = zeros( 3, 8 ) ;for i=1:1:4B(1:3,(2*i-1):2*i) = N_x(i)0% 弹性模 量% 泊松比0;N_y(i);N_x(i);N_y(i) end return function N_x, N_y = N_xy( ie, xi, eta )% 计算形函数对整体坐标的 导数% 输入参数：%ie 单元号%xi,eta 局部坐标%返回值：%N_x - 在局部坐标处的形函数对x 坐标的导数%N_y - 在局部 坐标处的形函数对y 坐标的导数J = Jacobi( ie, xi, eta ) ;N_xi,N_eta = N_xieta( ie, xi, eta

25、 ) ; 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, 1 ;N_xi = zeros( 1, 4 ) ;N_eta = zeros( 1, 4 ) ;N_xi(1) = x(1)*(1

26、+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 J = Jacobi( ie, xi, eta )%计算雅克比矩阵%输入参数：%ie 单元号%xi,eta 局部坐标

27、%返回值：%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; x_eta, y_eta ;returnfunction AssembleStiffnessMatrix( ie, k

28、)% 把单元刚度矩阵集成到整 体刚度矩阵% 输入参数 :%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) ;end end endendreturnfunction enf = EquivalentDistPressure( i, j, pi, p

29、j )% 计算线性分布压力的等效节 点力% 输入参数 :%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.577350269189626 ;w = 1, 1 ;for i=1:1:length(x) N1=(1-x(i)/2;N2=(1+x(

30、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;end enf = sig ;returnfunction DisplayResults% 显示计算结果，并与解析 解结果比较global gNode gElement gMaterial gDelta gElementStress o

31、utfile gBC1 global gDFn 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)=gElementStress(n*(m-1)+j,1,k);endStress (n+1)*(m+1),k)=gElementStress(n*m,

32、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);N_yy_max=find(Stress=yy_max)-(m+1)*(n+1); yy_min=min(min(Stress(:,2);N

33、_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 ) ;fprintf(outfile, =n ) ; fprintf(outfile, 结点号 X 坐标 Y 坐标 n ) ; fprintf(outfile, n ) ;f

34、or i=1:size(gNode(:,1)fprintf(outfile, %4d , %8.4f , %8.4f n,i,gNode(i,1),gNode(i,2) ; endfprintf( outfile,=nnn ) ;fprintf(outfile, 表 1 节点位移、应力 n ) ;fprintf(outfile,=n ) ;fprintf(outfile, 结点号x- 位移y- 位移Stress-xxStress-yyStress-xy n ) ;fprintf(outfile, n ) ;for i=1:size(gNode(:,1)fprintf(outfile, %4d

35、%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) ; endfprintf( outfile,= =nnn ) ;fprintf(outfile, 表 2 单元节点序列 n ) ;fprintf(outfile, =n ) ; fprintf(outfile, 单元号 材料组 节点号(逆时针 排列) n ) ;fprintf(outfile, 1 2 3 4 n ) ;fprintf(outfile, n ) ;for i=1

36、:size(gElement(:,1) fprintf(outfile, %4d , %4d , %4d , %4d , %4d , %4dn,.i, gElement(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.

37、3e %8.3e%8.3en,.i,gMaterial(i,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%4dn,.i,gBC1(i,1),gBC1(i,2),gBC1(i,3) ;endf

38、printf(outfile, =nnn ) ;fprintf(outfile, 表 5 分布力 n ) ; fprintf(outfile, =n ) ;fprintf(outfile, 边界 结 点 1 结点 2 均布力 1 均布力 2 n ) ; 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(i,4) ;endnnn ) ;fprintf( outfile,fprintf( outfile, 表 6

39、 xx 向节点应力极值n ) ;fprintf(outfile, =n ) ;fprintf(outfile, 极大点 Stress_xx_max 极小 点 Stress_xx_minn ) ;fprintf(outfile, n ) ;fprintf(outfile,%4d %12.4e %4d %12.4en, .N_xx_max,xx_max,N_xx_min,xx_min) ; fprintf(outfile, =nnn) ;fprintf( outfile, 表 7 yy 向节点应 力极值n ) ;fprintf(outfile, =n ) ; fprintf(outfile, 极大

40、点 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, =nnn) ;fprintf( outfile,表 8 xy 向节点 应力极值 n ) ;fprintf(outfile, =n ) ; fprintf(outfile, 极大点 Stress_xy_max 极小 点 Stress_xy_min n ) ;fprintf(outfile,

41、 n ) ;fprintf(outfile,%4d %12.4e %4d%12.4en,.N_xy_max,xy_max,N_xy_min,xy_min) ; fprintf(outfile, =nnn) ; 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)=gDel

42、ta(2*gElement(e,ord(i); spt11(i)=gElementStress(e,ord(i),1); spt21(i)=gElementStress(e,ord(i),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

43、);disy(2)=disy1(2); spt1(2)=spt11(2);spt2(2)=spt21(2);spt3(2)=spt31(2);xpt(3)=xpt1(3);ypt(3)=ypt1(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)=spt3

44、1(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);spt3(5)=spt31(1);scale=1e1;% 画变形图 subplot(2,3,1) hold on plot(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+sc

45、ale*disx,ypt+scale*disy,spt1)colorbaraxis equalaxis (0 b 0 b)title(Stress-xx distribution)% axis off subplot(2,3,3) hold onfill(xpt+scale*disx,ypt+scale*disy,spt2)colorbaraxis equalaxis (0 b 0 b)title(Stress-yy distribution)% axis off subplot(2,3,4) hold onfill(xpt+scale*disx,ypt+scale*disy,spt3)col

46、orbaraxis equalaxis (0 b 0 b)title(Stress-xy distribution)% 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+r2/2/gNode(i*(n+1),2)2+3*r4/2/gNode(i*(n+1),2)4); D

47、(i) = p;end plot(B,A,r:*); plot(C,A,k:x);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(4*A(1,i)-sin(2*A(1,i)/2; D(1,i)=(xy_max+xy_min)/2;end plot(A,B,r:*); plot(A,C,k:x); plot(A,D,b-);title( 孔边缘剪切力变化曲线( 逆时针 ) xlabel( 角度(单位为弧度) ) yla