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

有限元程序的设计大作业 有限元程序的设计大作业 PAGE 有限元程序的设计大作业 不同板宽的孔边应力集中问题 姓名：胡宇 学号：28 摘要 本文采用MATLAB和FOTRAN四节点平面单元，利用有限元数值解法对不同板宽的孔边应力集中问题进行了数值模拟研究。对于不同的板宽系数（半板宽b/孔半径r），得到了不同的应力集中系数，并且与解析解进行了比较，验证了有限元解的正确性，并且得出了解析解的适用范围。 引言 通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中。 MATLAB是由美国MATHWORKS公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB部分 1，计算模型 本程序采用MATLAB编程，编制平面四边形四节点等参元程序，用以求解近似平面结构问题。 本程序的研究对象为中央开有小孔的长方形板，选取的材料参数为：板厚h=1、材料强度E= Pa、泊松比mu=。此外，为方便网格的划分和计算，本文所取板的长度与宽度相等。其孔半径为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 ) ; 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 gElement =zeros( m*n, 5 ) ; for i=1:m 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] ; End end gElement( :, 5 ) = 1 ; 以上源程序中gNode为节点坐标矩阵，gElement为单元节点排列矩阵。 图（1）为5x10的网格划分图，图（2）为结构的微变形图。 图（1） 图（2） 假定该板所受的外载为p= 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，孔边应力研究 无限大板宽的孔边应力集中问题的弹性力学的解析解为： 在孔边的y轴上有分布: 为了验证有限元解的正确性，我取10x20