有限元编程大作业报告

1、本科生实验报告书四节点等参单元有限元分析的FORTRAN程序目录1 问题概述 (1)2 四节点四边形等参单元介绍 (1)3 单元应力磨平方法介绍 (4)4 程序流程设计 (6)4.1 程序设计概述 4.2 流程图5 程序结构及程序说明 (8)6 程序应用及算例分析 (9)6.1 算例概述6.2 算例ANSYS求解 6.3 算例程序数值解6.4 算例分析7. 总结 (15)7. 附录 ()1. 问题概述等参单元是有限元方法中使用最广泛的单元类型。等参单元的位移模式和坐标变换均采用相同的形函数，这种坐标变换叫做等参变换。通过等参变换可以将自然（局部）坐标中几何形状规则的单元转换成总体（笛卡尔）坐标

2、中形状扭曲的单元，因而使得单元有较好的适应性。本问题首先对平面四节点四边形等参单元的形函数、应力矩阵和等效节点力矩阵、应力磨平公式等的推导和计算求解。并通过设计FORTRAN求解程序进行编程求解，最后给出算例（受集中荷载的悬臂梁）并进行求解，将解与ANSYS的解进行比较。在这个过程中，采用了高斯三点积分和高斯两点积分，这种积分方法的求解效率较高而且精度也较好。在问题的最后，尝试去分析引起数值解误差的原因，并分析四节点等参单元的若干特性。2. 四节点四边形等参单元介绍边长为2的正方形单元（如下图所示），在其形心处安置一个局部坐标。单元角结点i的坐标（i,i）分别为±1,因此单元四条边界

3、的方程可用简单公式±1=0和±1=0逐一给出。12340 图2-1 母单元 图2-2 四边形单元形函数Ni的表达式：Ni=1+i(1+i)/4位移函数： 坐标变换式： ，单元应变矩阵 式中单元节点的位移列阵； 这里记号，。根据复合函数求导规则，有 从而有 因此，单元内的应力可以表示为 应力矩阵可以写成分块形式对于平面应力情形， 单元刚度矩阵是其中积分采用三点高斯积分，（高斯积分点的总数），和是加权系数，和是单元内的坐标.。对于三点高斯积分，高斯积分点的位置： ,，。结构刚度矩阵结构结点荷载列阵 注意，对于上两式中的理解不是简单的叠加而是按照对应的自由度集成。总刚平衡方程从式

4、上式求出：将回代入和中，得到和。等效节点力(1) 集中力 将集中荷载作用点取为结点，集中荷载就可以直接转化为等效结点力。(2) 体积力 等效结点力按下式计算Fvie=FvxiFvyie=-11-11NiPvxPvyhJdd (i=1,2,3,4)其中积分采用高斯2点积分，（高斯积分点的总数），和或是加权系数，和是单元内的坐标。对于2点高斯积分，高斯积分点的位置： 1=1=-1/3 ; 1=1 ; 2=2=1/3 ; 2=1(3) 表面力 单元的ij边上两个结点的等效结点力按下式计算 Fsie=-11Nixi,+yi,yi,-xi,hd (i=3,4)3. 应力磨平方法介绍 为了减少改进应力结果

5、的工作量，可以采用单元应力的局部磨平。当单元足够小时，磨平可以在各个单元内进行。对于等参元来说，有了单元内应力局部磨平的方法，可以方便地利用精度较高的高斯积分点(最佳应力点)的应力值来改进等参元结点应力的近似性质。以二维单元为例，如果是二次等参元，插值函数中的完全多项式是二次，即p2。对于C0型单元m1。这时p - m+1=2,则在2×2个高斯积分点上有限元的应力计算值 有较高的精度。如果进行应力磨平时只要求得到4各角点的改进应力值，即使在计算位移时是8结点二次等参元，但进行应力磨平时，单元结点数可只取4个，即用二维双线性单元进行应力磨平。磨平式各应力分量分别进行，这时应力磨平插值函

6、数Ni'应采用双线性函数，即    （3.1）在此情况下，4个结点上的应力可由高斯点上的应力外推得到，即令在2×2高斯积分点上有。4个高斯积分点的座标(见图3-1)如下: 图3-1 二维等参单元将高斯点坐标代入(3.1)式得到下面的等式： (3.2） 式中等号左边的应力列阵是有限元中已求出的4个高斯点相应的应力分量；是磨平后应力的结点值；转换矩阵由(3.22)式的第二式代入高斯点坐标后的插值函数值构成。由（3.2）式求逆可得 （3.3）其中 （3.4）各应力分量均可用（3.3）求解。这种改进结点应力的方法亦称之为应力插值外推。求得改进的应力结点值后，如需要

7、求单元内部的应力值仍可按(3.1)第一式进行计算。 采用单元局部应力磨平的方法，对于同一结点，由不同相邻单元求得的应力改进值通常是不相同的。可把相关单元求得的改进结点值再取平均作为最后的结点应力值。4. 程序流程设计4.1 程序设计概述本程序的设计以P3单元的FORTRAN程序为基础，在其架构之上修改而来。整体的架构与P3单元相似，但是在应力应变矩阵、主单元刚度矩阵、整体刚度矩阵、等效结点力以及应力磨平的算法方面有着较大的差异。具体的算法方面，在刚度矩阵的计算用了3点高斯积分，在等效结点力的计算用了2点高斯积分。应力磨平采用了单元内应力局部磨平的方法，对于同一结点，由不同相邻单元求得的应力改进

8、值通常是不相同的。把相关单元求得的改进结点值再取平均作为最后的结点应力值。4.2 流程图5. 程序结构5.1 程序结构5.2 子程序功能说明P4INPUT 读入单元数据、节点数据、节点约束条件及各类单元属性DEA 形成主对角元素地址P4SSM 形成总体刚度矩阵P4STIFB 计算单元刚度矩阵P4BMATB 计算当前单元的应变矩阵P4MODB 计算当前单元的弹性矩阵P4DBE 计算应力矩阵P4XJACM 计算当前单元的形函数和雅可比行列式BOUNDARY 边界条件处理LDLT LDLT分解结构刚度矩阵P4LOAK 计算单元荷载矩阵形成节点总荷载向量用置0置1法处理边界条件时产生的荷载向量修正项，

9、修正节点总荷载向量SOLV 处理荷载向量,回代求结点位移P4STREB 计算单元应力和主应力P4GAUSTR 单元内应力局部磨平GSTREB 计算磨平后的单元应力P4SUMSTRS 计算磨平后各节点的应力6. 程序应用及算例分析6.1 算例概述考虑一个右端固定的悬臂梁如下图5-1，其材料参数为：=0.24kg/mm3 ，E=2e5MPa，=0.2；几何参数为：L=8mm, H=2mm, b=1mm。在左端上方有集中荷载F=-100N，不考虑悬臂梁的自重，求解悬臂梁的挠度以及各截面的应力。 图6-16.2 算例ANSYS求解挠度图图6-2X方向应力图6-3Y方向应力图6-46.3 算例程序数值解

10、将悬臂梁进行单元网格划分如下图，划分为8个单元15个结点。 图6-5将划分结果及相关参数输入程序的input文档，启动程序进行求解。求解结果如下： 表6-1 节点应力表6-2 节点位移6.4 算例分析为了验证P4求解程序的准确性，我们将悬臂梁最大受拉区（顶边）的位移解和应力解在ANSYS下的和在本程序下得到的进行比较如下图：图6-6 位移解的比较由图可以看出，本程序的位移数值解与ANSYS的解十分接近，最大的节点误差在11.417%左右。图6-7 受拉区x向应力解比较由图可以看出，本程序的应力数值解与ANSYS的解十分接近，最大的节点误差在15.39%左右。应力磨平结果分析可见对于四边形等参单

11、元应力磨平前后，局部坐标为（0，0）的点的应力值不变。7. 总结本程序采用的四结点平面等参单元具有较高的求解精度，对于平面的问题也有很好的适应性。P4等参单元的求解还是和P3三角形常应变单元相比还是比较复杂的。虽然求解的流程基本一致，但是等参单元的求解需要解决刚度阵和等效节点力求解时的数值积分问题，本程序中采用应用比较广泛的高斯求积法。关于高斯积分的阶数，在李人宪所著的有限元法基础中指出，“用数值积分代替精确积分, 无疑计算时会产生误差. 为尽量减少这一误差,希望采用尽可能高的数值积分阶次” 。文献亦同时指出:鉴于N 很大时计算费时, 建议在二维情况下, 四节点单元的N 取为2 , 八节点单元的N 取为3 较好。本程序则混合采用