平面四节点等参单元分析程序

变分原理与有限元大作业平面四节点等参单元分析程序

*****

学号：SQ***********

完成时间：2011-4-26

一、概述

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

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

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

二、编程思想

本程序采用C语言编程，编制平面四边形四节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。

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

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

三、平面四节点等参单元介绍

四节点等参单元实际单元与基本单元的映射关系如图3-1所示

坐标的映射关系为：

其位移模式和坐标的映射有相同的插值函数，形函数为：

单元应变矩阵为：

{}

x

y

xy

u

x

u

y

u v

y x ε

εε

γ⎧⎫

∂

⎪⎪

∂

⎧⎫⎪⎪

⎪⎪⎪⎪

∂

==

⎨⎬⎨⎬

∂

⎪⎪⎪⎪

⎩⎭⎪⎪

∂∂

+

⎪⎪

∂∂

⎩⎭上式一般简写为：

图3-1

{}[]{}B εδ=

其中[]B 的子块矩阵为

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

∂=⎨

⎬∂⎪⎪⎪⎪

∂∂⎪⎪∂∂⎩⎭

由于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 εεεηηη∂∂∂∂∂⎡⎤⎡⎤

⎡⎤⎡⎤⎢⎥

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

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

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

从而有

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

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

因此，单元应力矩阵为

{}[][]{}D B σδ=

单元刚度矩阵为

[][][][]T

e

A

K B D B hdxdy =⎰⎰

其中积分采用三点高斯积分，

33

11

,11

11

1

(,)()

(,)nip

i

j

i j i

i

i j i f d d f W f ξηξη

ωω

ξηξη--===∑∑∑⎰⎰

2nip n =（高斯积分点的总数），i ω和j ω或i W 是加权系数，i ξ和j η是单元内的坐标.。对

于三点高斯积分，高斯积分点的位置：

11 5.0ξω==，220.0,8.0ξω==

，

33 5.0ξω==。

结构刚度矩阵

e e

K K =∑

结构结点荷载列阵