平面四边形4结等参有限单元法.docVIP

有限元程序设计 平面四边形4结点等参有限单元法 程序设计 1、程序功能及特点 a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。 b.前处理采用网格自动划分技术，自动生成单元及结点信息。 b.能计算受集中力、自重体力、分布面力和静水压力的作用。 c.计算结点的位移和单元中心点的应力分量及其主应力。 d.后处理采取整体应力磨平求得各个结点的应力分量。 e.算例计算结果与ANSYS计算结果比较，并给出误差分析。 f.程序采用Visual Fortran 5.0编制而成。 2、程序流程及图框 图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――计算单元应力分量 OUTSTRE――输出单元应力分量 STIF――计算单元刚度矩阵 FDNX――计算形函数对整体坐标的导数，1,2,3,4。 FUN8――计算形函数及雅可比矩阵[J] SFUN ――应力磨平-单元下的‘K’ NCN‘ SCN――应力磨平-单元下的右端项系数‘CN‘ SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘ SUMSTRS――应力磨平-单元下的集成到总体的‘K‘ GAUSTRSS――高斯消元求磨平后的应力 3、输入数据及变量说明L，B，NNL，NNB，NM，NR L——矩形体长度 B--矩形体宽度 NNL--L方向上划分的结点数 NNB--B方向上划分的结点数 NM—单元材料类型数 NR——约束结点总数 （2） 结点约束信息。共NR条，每条依次输入： IP，IX，IY IP——结点号 IX、IY——分别为IP结点在x，y方向的约束情况，如果约束填0，如果自由填1。 （3）材料信息。共NM条，每条依次输入： JJ,（AE（I，JJ），I 1，4） JJ――材料类型号，（AE（I，JJ），I 1，4）――该材料的材料参数，共4个参数，排列顺序为：弹性模量、泊松比、容重、单元厚度。 （4） 荷载信息 荷载控制信息。共一条，3个数据 NCP，IZ NCP——受集中力作用的结点数 IZ——面力批数 若NCP 0，输入 IP，PX，PY IP——结点号 PX、PY——分别为IP结点x，y方向的集中力分量。 若IZ 0，输入面力荷载信息，共IZ批，按批输入： JS，NSE，（WG（I）I 1，4） JS——面力批号 NSE——第JS批面力受到面力作用的单元个数， （WG（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、构造插值函数： 本有限元计算采取的是四边形八结点等参元进行插值计算的。 直接调用教材115页3..3.21的结果，写出所有插值函数： ; 4.2位移插值函数及应变应力求解： 在有限元方法中单元的位移模式一般采用多项式作为近似函数，多项式的选取应由低次到高次的完备多项式。位移模式的选取一般为：u＝φβ。 ，φ为位移模式，β为广义坐标向量。根据方程求解得出广义坐标，可将位移函数表示成结点位移的函数，即 ,，写成矩阵的形式为： N称为插值函数矩阵或形函数矩阵，为单元结点位移列阵。 确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。单元应变为： B称为应变矩阵，L是平面问题的微分算子，其中： ， 根据物理方程可以求得单元应力， 其中，S称为应力矩阵，B是应变矩阵。由于是平面应力问题，E0和v0取为E和v，所以弹性矩阵 。 这部分内容参考了教材第2.2节。 4.3、等参元变换： 为了将局部（自然）坐标中几何形状规则的单元转换成总体坐标中几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立坐标转换： , 最方便的方法是将上式中表示成插值函数的形式， , 在有限元分析中，为建立求解方程，需要进行各个单元面积内的积分，其一般形式