《平面四边形4结点等参有限单元法.doc》由会员分享，可在线阅读，更多相关《平面四边形4结点等参有限单元法.doc（38页珍藏版）》请在金锄头文库上搜索。
1、有限元程序设计平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。b.前处理采用网格自动划分技术，自动生成单元及结点信息。b.能计算受集中力、自重体力、分布面力和静水压力的作用。c.计算结点的位移和单元中心点的应力分量及其主应力。d.后处理采取整体应力磨平求得各个结点的应力分量。e.算例计算结果与ANSYS计算结果比较，并给出误差分析。f.程序采用Visual Fortran 5.0编制而成。2、程序流程及图框图2- 程序流程图图2-子程序框图其中，各子程序的主要功能为：INPUT输入原始数据HUAFEN自动网格划分，形
2、成COOR(2,NP),X,Y的坐标值与单元信息CBAND形成主元素序号指示矩阵MA（*）SKO形成整体刚度矩阵KCONCR计算集中力引起的等效结点荷载RBODYR计算自重体力引起的等效结点荷载RFACER计算分布面力引起的等效结点荷载RDECOP支配方程LU三角分解FOBALU分解直接解法中的回代过程OUTDISP输出结点位移分量STRESS计算单元应力分量OUTSTRE输出单元应力分量STIF计算单元刚度矩阵FDNX计算形函数对整体坐标的导数，1,2,3,4。FUN8计算形函数及雅可比矩阵JSFUN 应力磨平-单元下的K=NCNSCN应力磨平-单元下的右端项系数CNSUMSKN应力磨平-单
3、元下的右端项集成到总体的PSUMSTRS应力磨平-单元下的集成到总体的KGAUSTRSS高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。数据文件包括如下内容：（1）总控信息。共一条，4个数据。L，B，NNL，NNB，NM，NRL矩形体长度B-矩形体宽度NNL-L方向上划分的结点数NNB-B方向上划分的结点数NM单元材料类型数NR约束结点总数（2） 结点约束信息。共NR条，每条依次输入：IP，IX，IYIP结点号IX、IY分别为IP结点在
4、x，y方向的约束情况，如果约束填0，如果自由填1。（3）材料信息。共NM条，每条依次输入：JJ,（AE（I，JJ），I=1，4）JJ材料类型号，（AE（I，JJ），I=1，4）该材料的材料参数，共4个参数，排列顺序为：弹性模量、泊松比、容重、单元厚度。（4） 荷载信息a. 荷载控制信息。共一条，3个数据NCP，IZNCP受集中力作用的结点数IZ面力批数b. 若NCP0，输入IP，PX，PYIP结点号PX、PY分别为IP结点x，y方向的集中力分量。c. 若IZ0，输入面力荷载信息，共IZ批，按批输入：JS，NSE，（WG（I）I=1，4）JS面力批号NSE第JS批面力受到面力作用的单元个数，（W
5、G（I），I1，4）该面力的特征参数共4个数据，第1个数为面力类型，填1表示受静水压力作用，填2表示受均布法向压力作用；第2个数为水压密度，如果是均布压力情况，就填均布压力的集度；第3个数为最高水位的y坐标，如果是均布压力情况，可以填任意数；第4个数为面力作用的单元面的面号，单元面号的规定见图2-3。（IEW（M），M1，NSE），IEW（*）为受面力作用的单元的单元号，共NSE个。 图3-1 单元结点编码与面号4、理论基础和求解过程4.1、构造插值函数：本有限元计算采取的是四边形八结点等参元进行插值计算的。直接调用教材115页3.3.21的结果，写出所有插值函数：; 4.2位移插值函数及应变
6、应力求解：在有限元方法中单元的位移模式一般采用多项式作为近似函数，多项式的选取应由低次到高次的完备多项式。位移模式的选取一般为：u。，为位移模式，为广义坐标向量。根据方程求解得出广义坐标，可将位移函数表示成结点位移的函数，即 ,，写成矩阵的形式为： N称为插值函数矩阵或形函数矩阵，为单元结点位移列阵。确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。单元应变为：B称为应变矩阵，L是平面问题的微分算子，其中：，根据物理方程可以求得单元应力，其中，S称为应力矩阵，B是应变矩阵。由于是平面应力问题，E0和v0取为E和v，所以弹性矩阵。这部分内容参考了教材第2.2节。4.3、等
7、参元变换：为了将局部（自然）坐标中几何形状规则的单元转换成总体坐标中几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立坐标转换：,最方便的方法是将上式中表示成插值函数的形式，, 在有限元分析中，为建立求解方程，需要进行各个单元面积内的积分，其一般形式为：而g中包含场函数对于总体坐标x，y的导数。由于场函数是用自然坐标表述的，又因为在自然坐标内的积分限是规格化的，因此希望能在自然坐标内按规格化的数值积分方法进行上述积分。为此需要建立两个坐标系内导数之间的变换关系。设, ，xi，yi，zi是结点在总体坐标内的坐标值，Ni为形状函数，实际上它也是用局部（自然）坐标表示的插值函数。对
8、于等参变换，结点数、结点号和插值函数都不变。函数Ni对的偏导数可以表示成： 集合成矩阵形式就是： 式中的J称为雅克比矩阵，利用, ，J可以表示为自然坐标的函数，即：,Ni对于x,y的偏导数可以用自然坐标显示地表示为:其中是的逆矩阵,它可以按照下式计算得到：,是J的伴随矩阵，它的元素是J的元素的代数余子式。4.4、总体应力磨平根据有限单元法第5章中的（5.3.15），即 其另一种矩阵形式为：令并利用以下表达式，和以及就可以像单元元刚度矩阵到总体刚度矩阵一样，求出，只不过，这里有Ke为12*12的矩阵，而“总刚”K为3*NP阶的矩阵，NP为结点数。5、算例应用本程序计算一个矩形土体受到均布荷载时的
9、位移和应力，如下图所示，三面约束的土体尺寸为40m*10m，取一半为20m*10m，E=1.0104 kN/m2，在右端受均布荷载 kN/m2，不考虑自重体力。给定网格大小为，有限元网格自动划分如图3-2所示，单元总数2000，节点总数2121。由于对称性，右端约束为滑移支座，只限制x方向位移。虽然土体一般不为弹性，但是方法是一样的，只是刚度矩阵改动即可使用弹塑性模型。图5-1 算例模型图5-2 ANSYS计算模型图实际计算结果与ANSYS分析结果的比较（详细计算结果见后面）：（1） 位移比较比较项目X方向最大位移Y方向最大位移本程序计算结果0.0863mm5.890mmANSYS分析结果0.
10、1200mm6.000mm误差28.0%1.83%（2）应力比较比较项目最大主应力本程序计算结果9.996ANSYS分析结果9.940误差0.56%误差分析：本程序计算得到的Y方向最大位移和ansys计算结果很相近，由于x方向上的位移不占主要部分，因此，其误差虽然有些大，但对总体位移影响不大。 由于网格较密，因此，计算得到的单元应力值（未磨平）与ansys结果相近，误差小于1%，不必要应力磨平也可以达到较好的精度。 本程序可以进行总体应力磨平，但是由于网格数较多，因此，磨平矩阵阶数较大，可能计算时间也很长，甚至无法计算。（3）实际计算结果MATLAB图与ANSYS图比较：(a)(b)图5-3
11、X方向应力图，（a)ANSYS （b) MATLAB 图5-4 主应力计算结果的MATLAB 图图5-5 主应力计算结果的ANSYS 图图5-6 剪应力计算结果的MATLAB 图图5-7 剪应力计算结果的ANSYS 图图5-8 ANSYS总位移图图5-9 ANSYS X方向位移图图5-10 ANSYS Y方向位移图图5-11 ANSYS Y方向应力图附录：（1）输入文件数据：平面四边形四结点单元输入数据L B NNL NNB NM NR20.0 10.0 101 21 1 141*约束信息数据结点号 X向约束信息 Y向约束信息1002003004005006007008009001000110
12、012001300140015001600170018001900200021002200430064008500106001270014800169001900021100232002530027400295003160033700358003790040000421004420046300484005050052600547005680058900610006310065200673006940071500736007570077800799008200084100862008830090400925009460096700988001009001030001051001072001093001114001135001156001177001198001219001240001261001282001303001324001345001366001387001408001429001450001471001492001513001534001555001576001597001618001639001660001681001702001723001744001765001786
