有限元程序设计.doc

前言
有限单元法是在当今工程分析中获得最广泛应用的数值计算方法。由于它的通用性和有效性,受到工程技术界的高度重视。伴随着计算机科学与技术的快速发展,现已成为计算机辅助设计与计算机辅助制造的重要组成部分。
由于有限元法是通过计算机实现的,因此它的计算机程序的研制和开发是其理论和方法应用于生产和科研实际的前提和基础。同时所研制和开发的计算机程序又是有限元理论的和方法的研究必要平台。本次程序设计是将Fortran程序设计与有限元理论结合。根据有限元理论知识,进行设计程序,从而获得简单平面问题的计算方法。
平面4、8节点有限元公式及计算原理
(1)通过Serendipity四边形单元格式构造插值函数。对于4节点单元,插值函数为:
(i=1,2,3,4)
对于8节点单元,插值函数为:



(i=1,2,3,4)
(2)通过,

计算出D和B矩阵。
(3)运用高斯积分公式求积分公式近似数值积分解:
其中权函数
(4)单元矩阵的变换:
,因此
所以面积微元可表示成。
(5)单元刚度矩阵与载荷列向量可表示为:
,

将其代入整体刚度矩阵与载荷列向量可得整体刚度矩阵K和载荷列向量P。
(6)模型泛函总势能为:。
工程力学系
根据最小势能原理可得:
运用上式即可求得求得各节点位移。
程序结构
(1)有限元程序系统的组成及分析过程前处理
几何模型、材料参数、边界条件、分析类型定义、网格划分等。
核心:各种计算方法的实现。
有限元分析本体程序
以图形、曲线和表格等方式表达、分析计算结果。
后处理
(2)程序框图
输入离散模型数据型数据
计算单元刚度阵
组集结构刚度矩阵
计算单元等效结点载荷
组集结构结点载荷列阵
引入位移边界条件
求解线性方程组
其它辅助计算
输出结果
结束
单元循环
形成 K
形成 P
消除 K 的奇异
求解 Ka=P,得结点位移 a
计算应力、应变等
(3)变量说明
NUMNP:节点数;NUMEL:单元数;NUMMAT:材料种类;NLOAD:载荷数;NDM:每个节点坐标方向数;NDF:节点自由度;NEN:每个单元节点数。
ID(K,J):节点J的第K个自由度约束编号;X(K,J):节点J的第K个坐标。
IX(K,J):J单元中K节点的全局节点编号。
EE:杨氏弹性模量;XNU:泊松比;ITYPE:问题类型。
F(K,J):节点J在K方向上集中载荷大小。
ALFL=TRUE:非对称矩阵集合;ALFL=FALSE:对称矩阵集合:S:单元刚度矩阵;
P:载荷和内力向量;AD:对角元素;AU:上三角元素;AL:下三角元素;
JP:上(下)三角元素每行(列)最后元素;LD:一个单元中各自由度方程数。
JD:列高;ID:边界条件。
(4)数据输入格式
整数:I5格式;小数:。
(5)程序本体
C----------------------------------------
C
C..... FEA2DP---A FINITE ELEMENT ANALYSIS PROGRAM FOR
C 2D ELASTIC PROBLEMS
C
C TANGENT MATRIX IS STORED WITH VARIOUD BAND METHOD
C THIS PROGRAM IS USED TO DEMONSTRTE THE USAGE OF VRIOUS BAND
C STORAGE SCHEM OF SYMMETRIC AND UNSYMMETRIC TANGENT MATRIX
C
C Chen Yang
C AT CHONGQING VNIVERSITY(15/06/2011)
C
C------------------------------------------------------
PROGRAM FEA2DP
C
C A(1)-A(N1-1):X(NDM,NUMMNP); A(N1)-A(N2-1):F(NDF,NUMNP)
C A(N2)-A(N3-1):B(NEQ); A(N3)-A(N4-1):AD(NEQ)
C A(N4)-A(N5-1):AL(NAD); A(N5)-A(N6-1):NU(NAD)
C
C IA(1)-IA(N1-1):IX(NEN1,NUMEL); IA(N1)-IA(N2-1):ID(NDF,NUMNP)
C IA(N2)-IA(N3-1):JD((NDF*NUMNP); IA(N3)-IA(N4-1):IDL(NEN*NUMEL*NDF)
C
IMPLICIT REAL*8(