有限元程序的设计大作业.doc

《有限元程序的设计大作业.doc》由会员分享，可在线阅读，更多相关《有限元程序的设计大作业.doc（28页珍藏版）》请在三一办公上搜索。

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

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

3、0e11 Pa、泊松比mu=0.3。此外，为方便网格的划分和计算，本文所取板的长度与宽度相等。其孔半径为r=1，板宽为2b待定。由于本程序的目的在于验证有限元解的正确性和确定解析解的适用围，因此要求网格足够细密，以满足程序的精度要求。同时为了减小计算量，我采取网格径向长度递增的网格划分方法。此种方法特点是，靠近小孔部分的网格细密，在远离小孔的过程中，网格逐渐变得稀疏。以下为网格节点坐标计算和单元节点排序的MATLAB源程序：dr = (b-r)/(m+m2/2+m3/6) ;dfi= (pi/2)/n ;gNode = zeros( (n+1)*(m+1), 2 ) ;for i=1:1:mf

4、or j=1:1:n+1 gNode( (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+n/2),b;endgElement =zeros( m*n, 5 ) ;for i=1:m

5、for j=1:n gElement( (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）为结构的微变形图。图（1）图（2）假定该板所受的外载为p=1.0e9 Pa。由于是取1/4板进行计算，需要在分界面上加上合适的边界条件。以小孔圆心为原点，板长方向为x轴，板宽方向为y轴。我所取得边界条件为：与x轴平行的分界面上的节点的

6、y向固定，x向可以移动；与y轴平行的分界面上的节点的x向固定，y向可以移动。划分好网格，约束和载荷加好后。就可以计算了。取b=10，m=5，n=10时的应力图。图（3）为x向应力分布图，图（4）为y向应力分布图，图（5）为剪切力分布图。图（3）图（4）图（5）2，孔边应力研究无限大板宽的孔边应力集中问题的弹性力学的解析解为：在孔边的y轴上有分布:为了验证有限元解的正确性，我取10x20的网格进行了计算，得到了y轴上的x向应力分布曲线，并且与解析解进行了比较。如图（6）。y轴上的x向应力解析解为：图（6）同时，小孔边缘的剪切力分布如图（7）所示。剪切力的解析解为图（7）图中，红色*点为有限元节点

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

8、1504972.6526282.9824表中，,为径向单元的数目。,分别为与之对应的应力集中系数（）。为半板宽b与孔半径r的比值，即。解析解中的取值为3.下图（7）为解析解和有限元解的比较。图（7）图中黑色o点为理论应力集中系数，红色x点为较低网格密度下的应力集中系数变化曲线，蓝色*点为较高网格密度下的应力集中系数变化曲线。由图可知，同为有限元解，网格密度越高，精确度就越高。同时，我们还可以知道，解析解并不是适用于所有的情况。当板宽与孔径的比值小于5时，解析解与有限元解的差别变大，即解析解不再适用。4. FORTRAN部分 为了进一步验证结果的正确性，我同时运用FORTRAN对一样的网格进行了

9、计算。计算结果表明，有限元法所得的结果是可信的。5. 结论有限元解和解析解各有优劣，有限元解的优点是基本上对任意情况都适用，但是一般情况下计算量都很大，结果获取速度慢。解析解则求解迅速，但并不适用于所有情况。本文所研究的孔边应力集中问题就很好的体现了两种方法的特点。解析解并不适用于半板宽和孔半径之比小于5的情况，但是对于板宽远大于孔径的情况具有很高的精确度。而有限元解的精确度依赖于网格的密度和算法，在算法一样的情况下，为增加精度而加密网格会使计算量增加，在精度要求不高的情况下适用。附录一 MATLAB源程序function huyu2120120128 global r b m n p E m

10、u% 定义板的尺寸 r = 1 ; % 小孔半径 b = 10 ; %半板宽 format short e% 定义网格密度 m = 1; % 径向划分的数目 n = 2; % 周向划分的数目（必须为偶数）n=2*m为最优% 定义材料性质 E = 1e11 ; % 弹性模量 mu =0.3; % 泊松比% 定义均布荷载大小 p = 1.0e9; %单位Pa% FemModel ; % 定义有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 把计算结果显示出来 returnfunction FemModel% 定义有限元模型% 全局变量与名称% r-圆孔

11、半径% b-板宽% m - 半径方向单元数目 % n - 圆周方向单元数目% E - 弹性模量% mu - 泊松比% p - 均布载荷大小% gNode - 节点坐标% gElement - 单元定义% gMaterial - 材料性质% gBC1 - 第一类约束条件% gBC3 - 斜约束% gDF - 分布力% 返回值：% 无 global gNode gElement gMaterial gBC1 gDF global r b m n p E mu% 计算结点坐标 dr = (b-r)/(m+m2/2+m3/6) ; dfi= (pi/2)/n ; gNode = zeros( (n+1

12、)*(m+1), 2 ) ; for i=1:1:m for j=1:1:n+1 gNode( (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) ; end end for i=1:1:(n/2+1) gNode(n+1)*m+i, : )=b,b*tan(dfi*(i-1); end for i=1:1:n/2 gNode(n+1)*m+(n/2+1)+i, : )=b/tan(dfi*(i+n/2),b; end%

13、单元定义 gElement =zeros( m*n, 5 ) ; for i=1:m for j=1:n gElement( (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 end gElement( :, 5 ) = 1 ;% 定义材料性质 gMaterial = E, mu, 0 ; % 定义约束条件 gBC1 = zeros( (m+1)*2, 3 ) ; for i=1:1:m+1 gBC1(i, 1: 2)=(n+1)*(i-1)+1,(n+1)*(i-1)+1)*2;

14、 end for 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/2 gDF(i,1:2)=(n+1)*m+i,(n+1)*m+i+1; end gDF( :, 3 )=p; gDF( :, 4 )=p;returnfunction SolveModel% 求解有限元模型% 说明：% 该函数求解有限元模型，过程如下% 1. 计算单元刚度矩阵，集成整体刚度矩阵% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理

15、约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量gMaterial global gNode gElement gBC1 gDF gK gDelta gElementStress % step1. 定义整体刚度矩阵和节点力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number * 2, node_number * 2 ) ; f = sparse( node_number * 2, 1 ) ; % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中 element_number,dumm

16、y = size( gElement ) ; for ie=1:1:element_number fprintf( sprintf( 计算单元刚度矩阵并集成，当前单元: %dn, ie ) ) ; k = StiffnessMatrix( ie ) ; AssembleStiffnessMatrix( ie, k ) ; end% step3. 计算分布压力的等效节点力，并集成到整体节点力向量中 df_number,dummy = size( gDF ) ; for idf = 1:1:df_number enf = EquivalentDistPressure( gDF(idf,1), gD

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

18、1:bc1_number g = 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:element_number xi = -1 1 1 -1 ; eta = -1 -1 1 1 ; for n=

19、1:4 B = 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 ; end endreturnfunction k = StiffnessMatrix( ie ) % 计算平面应力等参数单元的刚度矩阵% 输入参数：% ie - 单元号% 返回值：%

20、 k - 单元刚度矩阵% 说明：% 用高斯积分法求解平面等参数单元的刚度矩阵 k = zeros( 8, 8 ) ; D = MatrixD( ie,1) ; % 3 x 3 高斯积分点和权系数 %x = -0.1483, 0.0,0.1483 ; %w = 0.5556,0.8889,0.5556 ; % 2 x 2 高斯积分点和权系数 x = -0.9626, 0.9626 ; w = 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(

21、j) ) ; k = k + w(i)*w(j)*transpose(B)*D*B*det(J) ; end endreturnfunction D = MatrixD( ie, opt )% 计算单元的弹性矩阵D% 输入参数：% ie - 单元号% opt - 问题选项% 1 = 平面应力% 2 = 平面应变% 返回值：% D - 弹性矩阵D global gElement gMaterial E = gMaterial( gElement(ie, 5), 1 ) ; % 弹性模量 mu = gMaterial( gElement(ie, 5), 2 ) ; % 泊松比 if opt = 1

22、 % 平面应力的弹性常数 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) ; end D = A3* 1 A1 0; A1 1 0; 0 0 A2 ;returnfunction B = MatrixB( ie, xi, eta )% 计算单元的应变矩阵B% 输入参数：% ie - 单元号% xi,eta - 局部坐标 % 返回值：% B - 在局部坐标处的应变矩阵B N_x,N_

23、y = N_xy( ie, xi, eta ); B = zeros( 3, 8 ) ; for i=1:1:4 B(1:3,(2*i-1):2*i) = N_x(i) 0; 0 N_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 - 在局部坐标处的形函数对y坐标的导数 J = Jacobi( ie, xi, eta ) ; N_xi,N

24、_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, 1 ; N_xi = zeros( 1, 4 ) ; N_

25、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 J = Jacobi( ie,

26、xi, eta )% 计算雅克比矩阵% 输入参数：% ie - 单元号% xi,eta - 局部坐标 % 返回值：% J - 在局部坐标(xi,eta)处的雅克比矩阵 global gNode gElement x = 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;

27、x_eta, y_eta ;returnfunction AssembleStiffnessMatrix( ie, k )% 把单元刚度矩阵集成到整体刚度矩阵% 输入参数:% ie - 单元号% k - 单元刚度矩阵% 返回值:% 无 global gElement gK for i=1:1:4 for j=1:1:4 for p=1:1:2 for q=1:1:2 s = (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

28、) ; end end end endreturnfunction 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.9

29、626, 0.9626 ; 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; end enf = sig ;returnfunction DisplayResults%

30、 显示计算结果，并与解析解结果比较 global gNode gElement gMaterial gDelta gElementStress outfile gBC1 global gDF n m r p b Stress outfile=fopen(huyu2120120128.txt,w); Stress = zeros(m+1)*(n+1),3); %集成节点应力矩阵 for k=1:1:3 for i=1:1:m Stress (n+1)*(i-1)+1,k)=gElementStress(n*(i-1)+1,4,k); end for j=1:1:n Stress (n+1)*m+

31、j,k)=gElementStress(n*(m-1)+j,1,k); end Stress (n+1)*(m+1),k)=gElementStress(n*m,2,k); for i=1:1:m for j=1:1:n Stress (i-1)*(n+1)+j+1,k)=gElementStress(i-1)*n+j,3,k); end end end xx_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_ma

32、x=max(max(Stress(:,2); 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 节

33、点坐标（mm) n ) ; fprintf(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, =n ) ; fprintf(outfile, 结点号 x-位移

34、y-位移 Stress-xx Stress-yy Stress-xy n ) ; 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,=nnn ) ; fprintf(outfile, 表2 单元节点序列 n ) ; fprintf(outfile, =n ) ; fprintf(outfile, 单元号 材料组 节点号（逆时针排列） n ) ; fprintf(outfile, 1 2 3 4 n ) ; fprintf(outfile, -n ) ; for i=1:size(