有限元程序设计报告报告

有限元程序设计报告课程名称:有限元程序设计指导教师：张*学校：重庆大学专业：工程力学01班******学号：********2015年7月8日有限元程序设计报告一、前言有限元方法（the Finite Element Method）是起源与上个世纪50、60 年代，基于弹性力学变分原理的一种近似计算方法，也是当今工程分析中获得最广泛应用的数值计算方法。

由于它的通用性和有效性，受到工程技术界的高度重视。

伴随着计算机科学和技术的快速发展，现已成为计算机辅助设计（CAD）和计算机辅助制造（CAM）的重要组成部分。

有限元程序系统通常包括前处理、有限元程序本体和后处理三部分。

前处理包括几何实体模型的建立、材料参数的赋值、位移边界条件的定义、载荷的定义、分析问题类型的定义、单元类型的选择和网格的划分等。

（分析问题类型如静力分析、动力特性分析、动力响应、温度场分析、电磁场分析、流体动力学分析等）有限元程序本体是有限元程序系统的核心部分，其功能是实现各种问题的计算。

后处理则是将计算结果用图形、曲线和表格的形式表达。

（通常包括结构的变形图、应力、应变分布云图等）本课程设计则是针对有限元程序本体，参照教学程序（FEATP），编写简单的有限元程序以计算简单的平面应力、平面应变和轴对称问题，并将其结果与有限元商用软件（ANSYS）的计算结果，以及问题的理论值进行比较，从而验证程序以及问题模型建立的正确性。

1．设计目的1）通过编写简单的有限元程序熟悉用有限元方法解决实际问题的基本步骤和过程，体会这种方法的处理手段。

2）在V isual Fortran 中编写程序，熟悉并巩固F ortran 语言的语法、算法，学习程序的调试方法，并体会其在执行某个具体算例时，文件的输入、输出以及程序的执行过程。

2．设计内容1）以教学程序（FEATP）为参照，编写程序，计算简单的平面应力（Plane Stress），平面应变(Plane Strain)问题，验证程序的正确性。

2）在具体的算例中，对同一问题，在程序和ANSYS 中采用不同的单元和网格划分方式，将其结果与理论值进行对比，体会不同的单元和网格划分对问题解的影响，从而判断模型的正确性和合理性。

3）总结在编写程序和算例中遇到的问题和解决方法，写出自己的心得体会。

二、弹性力学平面问题有限元方法的基本公式⎢ ⎢ ⎥ im 平面问题1．三角形单元(1) 位移⎡ u i ⎤ ⎢ ⎥ ⎢ v i ⎥ ⎡ ⎤ ⎡⎤⎢ ⎥ uN i u = ⎢ ⎥ = ⎢0 N j 0 N m 0 u j⎥⎢ ⎥ ⎣v ⎦ ⎣ 0N iN jN m ⎦⎢ v j ⎥⎢ ⎥⎡ a i ⎤ u m ⎥ ⎢⎣v m ⎥⎦= [IN IN j IN m ]⎢ a j ⎥= [N N ⎢⎣a m ⎥⎦N ]a e = Na e其中，形函数a i ,b i ,c i 是取决于节点坐标的常数应变矩阵 B 的分块子矩阵是⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂==i i i i i i i i b c c b A N N x y y x LN B 002100003 节点单元的应变矩阵是1⎥(2) 形成单元的刚度矩阵和等效节点载荷列阵i j()()m j i y i c x i b A N ,,ia 21i ++=[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==m m j j i i m j i m j i m j i b c b c b c c c c b b b A B B B B 00000021苏世宏《有限元课程设计》⎰B T ⎢ [i ⎥im m j m i j mK e = B T DBdV = B T DBtAV eP e = P e + P e + P e + P efsσ 0 ε0(3) 集成结构的刚度矩阵和等效节点载荷列阵K = ∑ K e = ∑⎰V DBdVeP = P f e+ P s P σ e+ P 0+ P F ∑f S σ 0 ε 0F= (P e + P eP ee+ P e) + P(4) 引入强制边界条件（消除 K 的奇异性） (5) 求解有限元解方程，得到节点位移 aKa = P(6) 计算单元应变和应力⎡ ε x ⎤ε = ⎢ ε ⎢γ ⎥e ⎥ = Lu = LNa = L N i N j⎥ N ]a e⎣ = [B xy ⎦B j B ]a e = Ba e⎡σ x ⎤⎢ ⎢σ y ⎢τ ⎥ = D ε = DBa e = Sa e⎥ ⎣ xy ⎦应力矩阵 S = DB = D [B B B ]= [S S S ]S 的分块子矩阵为⎡ ⎢ b i⎤ν 0 c i ⎥S i = DB i =E 0⎢2⎢ ν 0b ic ⎥（i, j, m ）i⎥2(1 -ν 0 ) A ⎢1 -ν 0 1 -ν 0 ⎥⎣⎢ 2 c i 2b i⎥⎦对于平面应力问题σ = 0 ε yE 0 = E ,ν 0 =ν对于平面应变问题νννν-=-=11020E E3． 四边形单元 用相同的推导方法可以得到四边形单元位移、应变、应力的有限元表达格式，它具有和三角形单元相似的表达形式，这里就不一一列举了。

三、有限元程序设计1．程序功能本程序是在教学程序（FEATP ）的基础上修改，删减而成的，主要应用于解决各项同性 的弹性力学二维问题，平面应力问题，平面应变问题和轴对称问题。

对于原程序中的动态响 应问题（DYNAM ）和动力特性分析问题（EIGENVALUE PROBLEM ）在本程序中将不涉及。

（1）问题类型①平面应力问题（MPROB=1） ②平面应变问题（MPROB=2） ③轴对称问题（MPROB=3） （2）单元类型①3—6 节点三角形单元（NODE=3 或 6） ②4—8 节点四边形单元（NODE=4 或 8） ③9 节点四边形单元（NODE=9）（3）求解类型 静力平衡分析：等带宽三角分解法（MOSLV=1）2．程序框图（1）程序总体框图形成 K的奇异=P ，得结点位移 a（2）程序调用框图3.输入文件变量名说明输入（1） M ND ：计算模型的各类单元中最多的节点数； N UMEL ：计算模型的单元总数； NUMPT ：计算模型的节点总数； MBAND ：半带宽（包括主对角元素）。

输入（2） NFIX ：有位移约束的节点数； NPC ：等效载荷作用的节点数； M PROB ：问题类型； MSOLV ：分析类型。

结束输出结果输入离散模型数据引入位移边界条件求解线性方程组其它辅助计算输入（3）NMATI：材料类型数；G RAV：重力加速度值；若不考虑重力，则输入0.0；M TYPE：输入控制参数；MTYPE＝0 输出全部计算结果；MTYPE＝1 输出除积分点应力以外的全部计算结果；MTYPE＝2 输出除总体质量矩阵、总体刚度矩阵和总体载荷向量以外的全部计算结果；MTYPE＝3 输出除积分点应力以及总体质量矩阵、总体刚度矩阵和总体载荷向量以外的全部计算结果。

输入（4）（4）是节点坐标信息的输入，即（II,(VCOOD(I,J),J=1,2,I=1,NUMPT)。

其中II：模型中的节点号，从1至N UMPT 依次按行输入；VCOOD(I,1)：II 节点处的x 向坐标。

VCOOD(I,2)：II 节点处的y 向坐标。

输入（5）。

其中（5）是单元信息的输入，即（II,(IELEM(I,J),J=1,4+MND),I=1,NUMEL）II：模型中的单元号，从1至N UMEL 依次按行输入。

IELEM(I,1)：II 单元的节点数。

IELEM(I,2)：II 单元的材料类型号。

IELEM(I,3)：II 单元沿x方向的高斯积分点数，对于三角形单元，则是H ammer 积分点数IELEM(I,4)：II 单元沿y方向的高斯积分点数。

对于三角形单元，填1。

IELEM(I,5)：IELEM(I,MND)：依次是II 单元的局部编号所对应的总体编号。

对于IELEM(I,1) 小于I ELEM （I,MND）的情况，在相应的位置上填0。

输入（6）。

（6）是位移约束信息的输入，即（II,(IFIXD(I,J),J=1,NF+1),(VFIXD(I,J),J=1,NF),I=1,NFIX）其中II：约束信息号，从1至N FIX 依次按行输入。

输入（7）7）是等效节点载荷信息的输入，即（II,(ILOAD(I,J),J=1,NF+1),(VLOAD(I,J),J=1,NF,I=1,NPC)。

其中，II：等效节点载荷号。

ILOAD(I,1)：第I I 个等效节点载荷作用的节点号。

ILOAD(I,2)～ILOAD(I,NF+1)：第II 个等效节点载荷所用节点的自由度开关，1 表示有载荷作用，0 表示没有载荷作用。

VLOAD(I,1)～VLOAD(I,NF):第II 个等效节点载荷作用于节点的自由度方向的载荷值大小。

对于平面问题和轴对称问题分别代表X,Y 方向的载荷值；对于M indlin 板分别代表θ x ,θ y和W方向的载荷值。

输入（8）（8）是材料类型和几何信息的输入，即（II，（VMATI(I,J),J=1,4）,I=1,NMATI）I I：材料类型号，从1至N MATI 依次按行输入；VMATI(I,1)：II 号材料的弹性模量（E）。

VMATI(I,2)：II 号材料的泊松比（v）。

VMATI(I,3)：II 号材料的质量密度（dens）。

VMATI(I,4)：II 号材料处板的厚度（th ）。

四、数值算例1． 算例一问题阐述：请采用4节点四边形等参单元对图1所示的无量纲L 型框架结构进行有限元分析。

材料杨氏模量和泊松比分别为7110E =⨯、0.33v =。

(1) 绘制出A 点水平位移u 随均布剪力P 取值变化（0~5000）的关系曲线及结构在3个典型载荷下的变形图；（建议长边采用40个单元，短边采用4个单元对结构进行离散）(2) 结合弹性力学小变形、线弹性假设，谈谈你对有限元分析结果的认识。

解：（1）程序编辑和程序调用以及演算过程演示：（i）整体程序（ii）子程序（stineffness，stress，input，output...）（iii）网格的划分：使用ABAQUS对结构进行离散，即用ABAQUS对结构进行网格划分，获得网格节点和单元的信息，然后形成“.txt”文件，保存到相应的文件目录里。

网格划分步骤如下：(a)安装ABAQUS并且启动ABAQUS（b）按照ABAQUS建模步骤，一步一步进行。

分别为：创建部件创建材料和截面属性定义装配件定义边界条件和载荷划分网格提交分析作业后处理退出ABAQUS/CAE.（c）利用Ultraedit软件，打开job里面的“.inp”文件，获取节点信息，转换成“.txt”问价输出。