有限元程序设计-平面四边形4结点等参有限单元法程序设计

《有限元程序设计-平面四边形4结点等参有限单元法程序设计》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、有限元程序设计-38-平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。b.前处理采用网格自动划分技术，自动生成单元及结点信息。b.能计算受集中力、自重体力、分布面力和静水压力的作用。c.计算结点的位移和单元中心点的应力分量及其主应力。d.后处理采取整体应力磨平求得各个结点的应力分量。e.算例计算结果与ANSYS计算结果比较，并给出误差分析。f.程序采用VisualFortran5.0编制而成。2、程序流程及图框图2-１程序流程图-38-图2

2、-２子程序框图其中，各子程序的主要功能为：INPUT――输入原始数据HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*）SKO――形成整体刚度矩阵[K]CONCR――计算集中力引起的等效结点荷载{R}BODYR――计算自重体力引起的等效结点荷载{R}FACER――计算分布面力引起的等效结点荷载{R}DECOP――支配方程LU三角分解FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量STRESS――计算单元应力分

3、量OUTSTRE――输出单元应力分量STIF――计算单元刚度矩阵FDNX――计算形函数对整体坐标的导数，1,2,3,4。FUN8――计算形函数及雅可比矩阵[J]SFUN――应力磨平-单元下的‘K’=NCN‘SCN――应力磨平-单元下的右端项系数‘CN‘SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘-38-SUMSTRS――应力磨平-单元下的集成到总体的‘K‘GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。在运行程序之前，根据程序中INPUT

4、需要的数据输入建立一个存放原始数据的文件，这个文件的名字为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结点在x，y方向的约束情况，如果约束填0，如果自由填1。（3）材料信息。共NM条，每条依次输入：JJ,（AE（I，JJ），I

5、=1，4）JJ――材料类型号，（AE（I，JJ），I=1，4）――该材料的材料参数，共4个参数，排列顺序为：弹性模量、泊松比、容重、单元厚度。（4）荷载信息a.荷载控制信息。共一条，3个数据NCP，IZNCP——受集中力作用的结点数IZ——面力批数b.若NCP>0，输入IP，PX，PYIP——结点号PX、PY——分别为IP结点x，y方向的集中力分量。c.若IZ>0，输入面力荷载信息，共IZ批，按批输入：JS，NSE，（WG（I）I=1，4）JS——面力批号-38-NSE——第JS批面力受到面力作用的单元个数，（WG（

6、I），I＝1，4）——该面力的特征参数共4个数据，第1个数为面力类型，填1表示受静水压力作用，填2表示受均布法向压力作用；第2个数为水压密度，如果是均布压力情况，就填均布压力的集度；第3个数为最高水位的y坐标，如果是均布压力情况，可以填任意数；第4个数为面力作用的单元面的面号，单元面号的规定见图2-3。（IEW（M），M＝1，NSE），IEW（*）为受面力作用的单元的单元号，共NSE个。图3-1单元结点编码与面号4、理论基础和求解过程4.1、构造插值函数：本有限元计算采取的是四边形八结点等参元进行插值计算的。直接调用

7、教材115页3..3.21的结果，写出所有插值函数：;4.2位移插值函数及应变应力求解：在有限元方法中单元的位移模式一般采用多项式作为近似函数，多项式的选取应由低次到高次的完备多项式。位移模式的选取一般为：u＝φβ。-38-，φ为位移模式，β为广义坐标向量。根据方程求解得出广义坐标，可将位移函数表示成结点位移的函数，即,，写成矩阵的形式为：N称为插值函数矩阵或形函数矩阵，为单元结点位移列阵。确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。单元应变为：B称为应变矩阵，L是平面问题的微分算子，其

8、中：，-38-根据物理方程可以求得单元应力，其中，S称为应力矩阵，B是应变矩阵。由于是平面应力问题，E0和v0取为E和v，所以弹性矩阵。这部分内容参考了教材第2.2节。4.3、等参元变换：为了将局部（自然）坐标中几何形状规则的单元转换成总体坐标中几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立坐标转换：,最方便的方法是将上式