本科生实验报告书四节点等参单元有限元分析的FORTRAN程序目录1．问题概述 (1)2．四节点四边形等参单元介绍 (1)3．单元应力磨平方法介绍 (4)4．程序流程设计 (6)程序设计概述流程图5．程序结构及程序说明 (8)6．程序应用及算例分析 (9)算例概述算例ANSYS求解算例程序数值解算例分析7. 总结 (15)1. 问题概述等参单元是有限元方法中使用最广泛的单元类型。
等参单元的位移模式和坐标变换均采用相同的形函数，这种坐标变换叫做等参变换。
通过等参变换可以将自然（局部）坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中形状扭曲的单元，因而使得单元有较好的适应性。
本问题首先对平面四节点四边形等参单元的形函数、应力矩阵和等效节点力矩阵、应力磨平公式等的推导和计算求解。
并通过设计FORTRAN 求解程序进行编程求解，最后给出算例（受集中荷载的悬臂梁）并进行求解，将解与ANSYS 的解进行比较。
在这个过程中，采用了高斯三点积分和高斯两点积分，这种积分方法的求解效率较高而且精度也较好。
在问题的最后，尝试去分析引起数值解误差的原因，并分析四节点等参单元的若干特性。
2. 四节点四边形等参单元介绍边长为2的正方形单元（如下图所示），在其形心处安置一个局部坐标。
单元角结点i 的坐标（,）分别为,因此单元四条边界的方程可用简单公式和逐一给出。
图2-1 母单元1 234 0图2-2 四边形单元形函数的表达式：位移函数： ∑==41i N ui i u ∑==41i N v i i v坐标变换式：∑==41i x N x i i ，∑==41i y N y i i单元应变矩阵{}[]{}[]{}ee B B B B B x v y u y v x u δδε4321==⎪⎪⎪⎭⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧∂∂+∂∂∂∂∂∂=式中{}[]TTT T T 4321e δδδδδ=——单元节点的位移列阵；[]{}），，，（432100,,,,=⎭⎬⎫⎩⎨⎧=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=i v u N N N N B i i i x i y i y i x i i δ这里记号x N i x ∂∂=/N ,i ，y /N y ,i ∂∂=i N 。
根据复合函数求导规则，有[]⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧y i x i y i x i i i N N J N N y y x x N N ,,,,,,,,,,ηξηξηξ从而有[]⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧-ηξ,,1,,i i y i x i N N J N N因此，单元内的应力可以表示为{}[][][][]{}e e xy y s B D δδτσσσ==⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=x应力矩阵[]S 可以写成分块形式[][]4321S S S S S =对于平面应力情形，[][][]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---==2/)1(2/)1(1,,,,,,2x i y i yi y i x i x i i i N N N N N N E B D S μμμμμ 单元刚度矩阵是[][][][]TeAK B D B hdxdy =⎰⎰其中积分采用三点高斯积分，),(),(111111j i j nip i nipj i f d d f ηξωωζξηξ⎰⎰∑∑--===2nip n =（高斯积分点的总数），i ω和j ω是加权系数，i ξ和j η是单元内的坐标.。
对于三点高斯积分，高斯积分点的位置： 6.01-=η,110.6, 5.09.0ξω=-=，0.02=η220.0,8.0ξω==，6.03=η330.6, 5.0ξω=。
结构刚度矩阵e eK K =∑结构结点荷载列阵e eP P =∑注意，对于上两式中e∑的理解不是简单的叠加而是按照对应的自由度集成。
总刚平衡方程[]{}{}K P δ=从式上式求出：{}[]{}1K P δ-=将{}δ回代入{}[][]{}D B σδ=和{}[]{}e B δε=中，得到{}ε和{}σ。
等效节点力(1) 集中力 将集中荷载作用点取为结点，集中荷载就可以直接转化为等效结点力。
(2) 体积力 等效结点力按下式计算其中积分采用高斯2点积分，),(),(111111j i j nip i nipj i f d d f ηξωωζξηξ⎰⎰∑∑--===2nip n =（高斯积分点的总数），i ω和j ω或i W 是加权系数，i ξ和j η是单元内的坐标。
对于2点高斯积分，高斯积分点的位置：;;;(3) 表面力 单元的ij 边上两个结点的等效结点力按下式计算3. 应力磨平方法介绍为了减少改进应力结果的工作量，可以采用单元应力的局部磨平。
当单元足够小时，磨平可以在各个单元内进行。
对于等参元来说，有了单元内应力局部磨平的方法，可以方便地利用精度较高的高斯积分点(最佳应力点)的应力值来改进等参元结点应力的近似性质。
以二维单元为例，如果是二次等参元，插值函数中的完全多项式是二次，即p＝2。
对于C0型单元m＝1。
这时p - m+1=2,则在2×2个高斯积分点上有限元的应力计算值有较高的精度。
如果进行应力磨平时只要求得到4各角点的改进应力值，即使在计算位移时是8结点二次等参元，但进行应力磨平时，单元结点数可只取4个，即用二维双线性单元进行应力磨平。
磨平式各应力分量分别进行，这时应力磨平插值函数Ni'应采用双线性函数，即（）在此情况下，4个结点上的应力可由高斯点上的应力外推得到，即令在2×2高斯积分点上有。
4个高斯积分点的座标(见图3-1)如下:图3-1 二维等参单元将高斯点坐标代入式得到下面的等式：）式中等号左边的应力列阵是有限元中已求出的4个高斯点相应的应力分量；是磨平后应力的结点值；转换矩阵由式的第二式代入高斯点坐标后的插值函数值构成。
由（）式求逆可得（）其中（）各应力分量均可用（）求解。
这种改进结点应力的方法亦称之为应力插值外推。
求得改进的应力结点值后，如需要求单元内部的应力值仍可按第一式进行计算。
采用单元局部应力磨平的方法，对于同一结点，由不同相邻单元求得的应力改进值通常是不相同的。
可把相关单元求得的改进结点值再取平均作为最后的结点应力值。
4.程序流程设计4.1程序设计概述本程序的设计以P3单元的FORTRAN程序为基础，在其架构之上修改而来。
整体的架构与P3单元相似，但是在应力应变矩阵、主单元刚度矩阵、整体刚度矩阵、等效结点力以及应力磨平的算法方面有着较大的差异。
具体的算法方面，在刚度矩阵的计算用了3点高斯积分，在等效结点力的计算用了2点高斯积分。
应力磨平采用了单元内应力局部磨平的方法，对于同一结点，由不同相邻单元求得的应力改进值通常是不相同的。
把相关单元求得的改进结点值再取平均作为最后的结点应力值。
4.2流程图5.程序结构5.1程序结构5.2子程序功能说明[P4INPUT] 读入单元数据、节点数据、节点约束条件及各类单元属性[DEA] 形成主对角元素地址[P4SSM] 形成总体刚度矩阵[P4STIFB] 计算单元刚度矩阵[P4BMATB] 计算当前单元的应变矩阵[P4MODB] 计算当前单元的弹性矩阵[P4DBE] 计算应力矩阵[P4XJACM] 计算当前单元的形函数和雅可比行列式[BOUNDARY] 边界条件处理[LDLT] LDLT分解结构刚度矩阵[P4LOAK] 计算单元荷载矩阵形成节点总荷载向量用置0置1法处理边界条件时产生的荷载向量修正项，修正节点总荷载向量[SOLV] 处理荷载向量,回代求结点位移[P4STREB] 计算单元应力和主应力[P4GAUSTR] 单元内应力局部磨平[GSTREB] 计算磨平后的单元应力[P4SUMSTRS] 计算磨平后各节点的应力6.程序应用及算例分析6.1算例概述考虑一个右端固定的悬臂梁如下图5-1，其材料参数为：，，；几何参数为：L=8mm, H=2mm, b=1mm。
在左端上方有集中荷载F=-100N，不考虑悬臂梁的自重，求解悬臂梁的挠度以及各截面的应力。
图6-16.2算例ANSYS求解挠度图图6-2 X方向应力图6-3Y方向应力图6-46.3算例程序数值解将悬臂梁进行单元网格划分如下图，划分为8个单元15个结点。
图6-5将划分结果及相关参数输入程序的input文档，启动程序进行求解。
求解结果如下：表6-1 节点应力表6-2 节点位移6.4算例分析为了验证P4求解程序的准确性，我们将悬臂梁最大受拉区（顶边）的位移解和应力解在ANSYS下的和在本程序下得到的进行比较如下图：图6-6 位移解的比较由图可以看出，本程序的位移数值解与ANSYS的解十分接近，最大的节点误差在%左右。
图6-7 受拉区x向应力解比较由图可以看出，本程序的应力数值解与ANSYS的解十分接近，最大的节点误差在%左右。
应力磨平结果分析可见对于四边形等参单元应力磨平前后，局部坐标为（0，0）的点的应力值不变。
7.总结本程序采用的四结点平面等参单元具有较高的求解精度，对于平面的问题也有很好的适应性。
P4等参单元的求解还是和P3三角形常应变单元相比还是比较复杂的。
虽然求解的流程基本一致，但是等参单元的求解需要解决刚度阵和等效节点力求解时的数值积分问题，本程序中采用应用比较广泛的高斯求积法。
关于高斯积分的阶数，在李人宪所著的《有限元法基础》中指出，“用数值积分代替精确积分, 无疑计算时会产生误差. 为尽量减少这一误差,希望采用尽可能高的数值积分阶次”。
文献亦同时指出:鉴于N 很大时计算费时, 建议在二维情况下, 四节点单元的N 取为2 , 八节点单元的N 取为3 较好。
本程序则混合采用二阶和三阶高斯积分以获得较好的精确性和较高的计算效率。
至于在求解本算例和ANSYS求解存在10%左右的误差主要是因为4节点等参单元缺少边中点，不能很好地适应悬臂梁的弯曲变形，加之单元的划分比较少（总共8个单元），网格的解就存在着一定的误差。
其次，在应力磨平方面基于4个结点上的应力可由高斯点上的应力外推得到的原理，可能会存在一定的误差。