八节点等参单元平面有限元

八节点等参单元平面有限元

分析程序

土木工程学院

2011.2

目录

1.概述 (1)

2.编程思想 (2)

2.1.八节点矩形单元介绍 (2)

2.2.有限元分析的模块组织 (5)

3.程序变量及函数说明 (6)

3.1.主要变量说明： (6)

3.2.主要函数说明 (7)

4.程序流程图 (8)

5.程序应用与ANSYS分析的比较 (9)

5.1.问题说明 (9)

5.2.ANSYS分析结果 (10)

5.3.自编程序分析结果 (12)

5.4.结果比较分析 (12)

参考文献 (14)

附录程序源代码 (15)

《计算力学》课程大作业

1.概述

通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。

随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。

C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。

2. 编程思想

本程序采用C 语言编程，编制平面四边形八节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移。

在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）和节点温度（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），单元类型信息（桁架、梁、板、壳等）（整型） ，单元特性信息（厚度、内力矩等）（实型），材料信息（弹性模量，泊松比等）（实型）等。

在FORTRAN 程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如需要进行单元网络的自适应划分，节点及单元的修改将非常困难。在进行C 语言编译过程中，采用结构struct 使每个节点信息存储在一个结构体数组中，提高程序的可读性，使数据结构更趋于合理。

2.1. 八节点矩形单元介绍

八节点矩形单元编号如图 2-1所示

八节点矩形单元的位移函数为：

222212345678u x y x xy y x y xy αααααααα=+++++++ (2.1)

2222910111213141516v x y x xy y x y xy αααααααα=+++++++ (2.2)

其形函数为

1122334455667788u N u N u N u N u N u N u N u N u =+++++++ (2.3) 1122334455667788v N v N v N v N v N v N v N v N v =+++++++ (2.4)

式(2.3)和式(2.4)中(,)i i N N εη=，并且采用等参变换，则单元的坐标变换式可取为

1122334455667788

1122334455667788X N x N x N x N x N x N x N x N x Y N y N y N y N y N y N y N y N y =+++++++⎧⎨

=+++++++⎩

(2.5)

图 2-1

单元应变矩阵为

{}x y xy u x u y u v y x εεεγ⎧⎫∂⎪⎪∂⎧⎫⎪⎪

⎪⎪⎪⎪∂==⎨⎬⎨⎬∂⎪⎪⎪⎪⎩⎭⎪⎪

∂∂+⎪⎪∂∂⎩⎭

(2.6)

式(2.6)一般简写为

{}[]{}B εδ= (2.7)

其中[]B 的子块矩阵为

[]i i i i i N x N B y N N y x ⎧⎫∂⎪⎪∂⎪⎪⎪

⎪

∂=⎨

⎬∂⎪⎪⎪⎪

∂∂⎪⎪∂∂⎩⎭

(2.8) 由于i N 是ε、η的函数，在[]i B 中的x 、y 要按照复合函数来求导，即

[]i i i i i i N N N x y x x J N N N x y y y εεεηηη∂∂∂∂∂⎡⎤⎡⎤

⎡⎤⎡⎤⎢⎥

⎢⎥⎢⎥⎢⎥∂∂∂∂∂⎢

⎥⎢⎥⎢⎥⎢⎥==∂∂∂∂∂⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢

⎥∂∂∂∂∂⎣⎦⎣⎦⎣⎦⎣⎦

(2.9) 从而有

[]1

i i i i N N x J N N y εη-∂∂⎡⎤⎡⎤

⎢⎥⎢⎥∂∂⎢⎥⎢⎥=∂∂⎢⎥⎢⎥⎢⎥⎢⎥∂∂⎣⎦⎣⎦

(2.10) 因此，单元应力矩阵为

{}[][]{}D B σδ= (2.11) 单元刚度矩阵为

[][][][]T

e

A

K B D B hdxdy =⎰⎰ (2.12)