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

有限元程序设计 平面四边形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——面力批数 若NCP0，输入 IP，PX，PY IP——结点号 PX、PY——分别为IP结点x，y方向的集中力分量。 若IZ0，输入面力荷载信息，共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、等参元变换： 为了将局部（自然）坐标中几何形状规则的单元转换成总体坐标中几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立坐标转换： , 最方便的方法是将上式中表示