有限元程序的设计大作业

1.不同板宽的孔边应力集中问题姓名：胡宇学号：21201201282.摘要本文采用MATLAB和FOTRAN四节点平面单元，利用有限元数值解法对不同板宽的孔边应力集中问题进行了数值模拟研究。对于不同的板宽系数（半板宽b/孔半径r），得到了不同的应力集中系数1，并且与解析解进行了比较，验证了有限元解的正确性，并且得出了解析解的适用范围。3.引言通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中。MATLAB是由美国MATHWORKS公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。4.MATLAB部分1，计算模型本程序采用MATLAB编程，编制平面四边形四节点等参元程序，用以求解近似平面结构问题。本程序的研究对象为中央开有小孔的长方形板，选取的材料参数为：板厚h=1、材料强度E=1.0e11Pa、泊松比mu=0.3。此外，为方便网格的划分和计算，本文所取板的长度与宽度相等。其孔半径为r=1，板宽为2b待定。由于本程序的目的在于验证有限元解的正确性和确定解析解的适用范围，因此要求网格足够细密，以满足程序的精度要求。同时为了减小计算量，我采取网格径向长度递增的网格划分方法。此种方法特点是，靠近小孔部分的网格细密，在远离小孔的过程中，网格逐渐变得稀疏。以下为网格节点坐标计算和单元节点排序的MATLAB源程序：dr=(b-r)/(m+m^2/2+m^3/6);dfi=(pi/2)/n;gNode=zeros((n+1)*(m+1),2);fori=1:1:mforj=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))];endendfori=1:1:(n/2+1)gNode((n+1)*m+i,:)=[b,b*tan(dfi*(i-1))];endfori=1:1:n/2gNode((n+1)*m+(n/2+1)+i,:)=[b/tan(dfi*(i+n/2)),b];endgElement=zeros(m*n,5);fori=1:mforj=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）为结构的微变形图。图（1）图（2）假定该板所受的外载为p=1.0e9Pa。由于是取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，孔边应力研究无限大板宽的孔边应力集中问题的弹性力学的解析解为：222244222222223112sin2312cos2123112cos212rrprprprrprp在孔边的y轴上有分布:parrp3,90,2321190044220为了验证有限元解的正确性，我取10x20的网格进行了计算，得到了y轴上的x向应力分布曲线，并且与解析解进行了比较。如图（6）。y轴上的x向应力解析解为：payxyryrpxxx3,0,2321104422图（6）同时，小孔边缘的剪切力分布如图（7）所示。剪切力的解析解为2sin4sin21pxy图（7）图中，红色*点为有限元节点解，黑色x点为解析解点。由图可知，有限元解和解析解基本相同，有限元解是有效的。同时，为了研究不同板宽对应力集中系数的影响，我选取了多组数据进行对比。如表（1）所示。表（1）bor（）b-r1m12m22116.028646.11143213.626464.03154323.228973.52485422.909483.31696522.712293.20987622.5751103.14708732.8527113.10699832.7762113.017010932.7102123.0604171642.6707163.0071262552.6589202.9911373662.6545242.9851504972.6526282.9824表中，1m,2m为径向单元的数目。1,2分别为与之对应的应力集中系数（p/max）。为半板宽b与孔半径r的比值，即rb/。解析解中的取值为3.下图（7）为解析解和有限元解的比较。图（7）图中黑色o点为理论应力集中系数，红色x点为较低网格密度下的应力集中系数变化曲线，蓝色*点为较高网格密度下的应力集中系数变化曲线。由图可知，同为有限元解，网格密度越高，精确度就越高。同时，我们还可以知道，解析解并不是适用于所有的情况。当板宽与孔径的比值小于5时，解析解与有限元解的差别变大，即解析解不再适用。5.FORTRAN部分为了进一步验证结果的正确性，我同时运用FORTRAN对相同的网格进行了计算。计算结果表明，有限元法所得的结果是可信的。6.结论有限元解和解析解各有优劣，有限元解的优点是基本上对任意情况都适用，但是一般情况下计算量都很大，结果获取速度慢。解析解则求解迅速，但并不适用于所有情况。本文所研究的孔边应力集中问题就很好的体现了两种方法的特点。解析解并不适用于半板宽和孔半径之比小于5的情况，但是对于板宽远大于孔径的情况具有很高的精确度。而有限元解的精确度依赖于网格的密度和算法，在算法相同的情况下，为增加精度而加密网格会使计算量增加，在精度要求不高的情况下适用。附录一MATLAB源程序functionhuyu2120120128globalrbmnpEmu%定义板的尺寸r=1;%小孔半径b=10;%半板宽formatshorte%定义网格密度m=1;%径向划分的数目n=2;%周向划分的数目（必须为偶数）n=2*m为最优%定义材料性质E=1e11;%弹性模量mu=0.3;%泊松比%定义均布荷载大小p=1.0e9;%单位Pa%FemModel;%定义有限元模型SolveModel;%求解有限元模型DisplayResults;%把计算结果显示出来returnfunctionFemModel%定义有限元模型%全局变量及名称%r-----圆孔半径%b-----板宽%m---------半径方向单元数目%n---------圆周方向单元数目%E---------弹性模量%mu---------泊松比%p----------均布载荷大小%gNode------节点坐标%gElement---单元定义%gMaterial--材料性质%gBC1-------第一类约束条件%gBC3-------斜约束%gDF--------分布力%返回值：%无globalgNodegElementgMaterialgBC1gDFglobalrbmnpEmu%计算结点坐标dr=(b-r)/(m+m^2/2+m^3/6);dfi=(pi/2)/n;gNode=zeros((n+1)*(m+1),2);fori=1:1:mforj=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))];endendfori=1:1:(n/2+1)gNode((n+1)*m+i,:)=[b,b*tan(dfi*(i-1))];endfori=1:1:n/2gNode((n+1)*m+(n/2+1)+i,:)=[b/tan(dfi*(i+n/2)),b];end%单元定义gElement=zeros(m*n,5);fori=1:mforj=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];endendgElement(:,5)=1;%定义材料性质gMaterial=[E,mu,0];%定义约束条件gBC1=zeros((m+1)*2,3);fori=1:1:m+1gBC1(i,1:2)=[(n+1)*(i-1)+1,((n+1)*(i-1)+1)*2];endfori=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);fori=1:1:n/2gDF(i,1:2)=[(n+1)*m+i,(n+1)*m+i+1];endgDF(:,3)=p;gDF(:,4)=p;returnfunctionSolveModel%求解有限元模型%说明：%该函数求解有限元模型，过程如下%1.计算单元刚度矩阵，集成整体刚度矩阵%2.计算单元的等效节点力，集成整体节点力向量%3.处理约束条件，修改整体刚度矩阵和节点力向量%4.求解方程组，得到整体节点位移向量gMaterialglobalgNodegElementgBC1gDFgKgDeltagElementStress%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(gElement);forie=1:1:element_numberfprintf(sprintf('计算单元刚度矩阵并集成，当前单元:%d\n',ie));k=StiffnessMatrix(ie);AssembleStiffnessMatrix(ie,k);end%step3.计算分布压力的等效节点力，并集成到整体节点力向量中[df_number,dummy]=size(gDF);foridf=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)+e