平面四边形4结点等参有限单元法

有限元程序设计

平面四边形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}e BODYR――计算自重体力引起的等效结点荷载{R}e FACER――计算分布面力引起的等效结点荷载{R}e DECOP――支配方程LU三角分解

FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量

STRESS――计算单元应力分量

OUTSTRE――输出单元应力分量

STIF――计算单元刚度矩阵

FDNX――计算形函数对整体坐标的导数

T

i

i

y

N

x

N

⎥

⎦

⎤

⎢

⎣

⎡

∂

∂

∂

∂

，=

i1,2,3,4。

FUN8――计算形函数及雅可比矩阵[J]

SFUN ――应力磨平-单元下的‘K’=NCN‘

SCN――应力磨平-单元下的右端项系数‘CN‘

SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘

SUMSTRS――应力磨平-单元下的集成到总体的‘K‘

GAUSTRSS――高斯消元求磨平后的应力

3、输入数据及变量说明

当程序开始运行时，按屏幕提示，键入数据文件的名字。

在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。数据文件包括如下内容：

（1）总控信息。共一条，4个数据。

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）荷载信息

a.荷载控制信息。共一条，3个数据

NCP，IZ

NCP——受集中力作用的结点数

IZ——面力批数

b. 若NCP >0，输入

IP ，PX ，PY IP ——结点号

PX 、PY ——分别为IP 结点x ，y 方向的集中力分量。

c. 若IZ >0，输入面力荷载信息，共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的结果，写出所有插值函数：

158111(1)(1)422N N N ξη=----; 256

111(1)(1)422N N N ξη=+--- 367111(1)(1)422N N N ξη=++-- 478

111(1)(1)422N N N ξη=-+-- 251(1)(1)2N ξη=-- 261(1)(1)

2N ηξ=-+ 271(1)(1)2N ξη=-+ 281(1)(1)

2N ηξ=--

4.2位移插值函数及应变应力求解：

在有限元方法中单元的位移模式一般采用多项式作为近似函数，多项式的选取应由低次到高次的完备多项式。位移模式的选取一般为：u ＝φβ。

u u v ⎡⎤

=⎢⎥⎣⎦

，φ为位移模式，β为广义坐标向量。根据方程求解得出广义坐标，可将位移函数表示成

结点位移的函数，即8

1

i i

i u N u

==

∑ ,8

1

i i

i v N v

==

∑，写成矩阵的形式为：