有限元程序设计-平面四边形4结点等参有限单元法程序设计

有限元程序设计

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

程序设计

1、程序功能及特点

a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。

b.前处理采用网格自动划分技术，自动生成单元及结点信息。

b.能计算受集中力、自重体力、分布面力和静水压力的作用。

c.计算结点的位移和单元中心点的应力分量及其主应力。

d.后处理采取整体应力磨平求得各个结点的应力分量。

e.算例计算结果与ANSYS计算结果比较，并给出误差分析。

£程序采用VisualFortran5.0编制而成。

2、程序流程及图框

图2-1程序流程图

图2-2子程序框图

其中，各子程序的主要功能为：

INPUT一—输入原始数据

HUAFEN——自动网格划分，形成C00R(2,NP),X,Y的坐标值与单元信息

CBAND一—形成主元素序号指示矩阵MA(*)

SKO-----形成整体刚度矩阵[K]

CONCR——计算集中力引起的等效结点荷载{R}'

BODYR一—计算自重体力引起的等效结点荷载{R}'

FACER一—计算分布面力引起的等效结点荷载{R}'

DECOP一一支配方程LU三角分解

FOBA——LU分解直接解法中的回代过程

OUTDISP一—输出结点位移分量

STRESS——计算单元应力分量

OUTSTRE—-输出单元应力分量

STIF——计算单元刚度矩阵

一dN.取丫

FDNX一—计算形函数对整体坐标的导数一—,i=1,2,3,4。

_dxdy_

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

SFUN——应力磨平-单元下的'K'=NCN'

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

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

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

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

3、输入数据及变量说明

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

在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，

这个文件的名字为INDAT.DATo数据文件包括如下内容：

(1)总控信息。共一条，4个数据。

L,B,NNL,NNB,NM,NR

L——矩形体长度

B—矩形体宽度

NNL—L方向上划分的结点数

NNB—B方向上划分的结点数

NM一单元材料类型数

NR——约束结点总数

(2)结点约束信息。共版条，每条依次输入：

IP,IX,IY

IP——结点号

IX、IY一一分别为IP结点在x,y方向的约束情况，如果约束填0,如果自由填1。

(3)材料信息。共A必条，每条依次输入：

JJ,(/£(/,〃)，1=1,4)

〃一一材料类型号，(4£(乙〃)，/=1,4)――该材料的材料参数，共4个参数,

排列顺序为：弹性模量、泊松比、容重、单元厚度.

(4)荷载信息

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

NCP,IZ

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

IZ一一面力批数

b.若NCP>0,输入

IP,PX,PY

IP——结点号

PX.PY-——分别为/P结点x,y方向的集中力分量。

c.若30,输入面力荷载信息，共/Z批，按批输入：

JS,NSE,(WG(D1=1,4)

JS----面力批号

NSE——第元批面力受到面力作用的单元个数，

(伤(/),/=1,4)一一该面力的特征参数共4个数据，第1个数为面力类型，填

1表示受静水压力作用，填2表示受均布法向压力作用；第2个数为水压密度，如果是均

布压力情况，就填均布压力的集度；第3个数为最高水位的y坐标，如果是均布压力情况，

可以填任意数；第4个数为面力作用的单元面的面号，单元面号的规定见图2-3。

(IEW(M),M=l,NSE),IEW(*)为受面力作用的单元的单元号，共NSE个。

n

图3-1单元结点编码与面号

4、理论基础和求解过程

4.1、构造插值函数：

本有限元计算采取的是四边形八结点等参元进行插值计算的。

直接调用教材115页3..3.21的结果，写出所有插值函数：

^=3+1)(1一叫

»

3=3+1)(1+叫乂彳凶乂=；(1)“+if7VM

/V=1(l-^2)(l-7)1

57M=5(I-91

1

29乂弓(1-丁)(1-4

7V7=-(l-^)(7+l)

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

在有限元方法中单元的位移模式一般采用多项式作为近似函数，多项式的选取应由低

次到高次的完备多项式。位移模式的选取一般为：u=<bB。

u

u=,6为位移模式，B为广义坐标向量。根据方程求解得出广义坐标，可将位

v

88

移函数表示成结点位移的函数，即〃=必V=±NM,写成矩阵的形式为:

1=1i=l

u=[N,N2...N8M=

N称为插值函数矩阵或形函数矩阵，/为单元结点位移列阵。

确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。

单元应变为:

£=%=Lu=LNae=L[N]N2...乂]浦=[4与...优=3/

B称为应变矩阵，L是平面问题的微分算子，其中:

ddN

0i0

dxdx

dNo一dNj

Bj=LN：=00

Sy0Nj_Sy

dd叫dNj

的dxdx

o

d0N20Mo

B=[BtB2

彳[0N]0N?o乂

d

Sy

a

aw

^。

-

&

-

_

4

±

M

o

av

£

-

-一

-一

8

V

^

%

包aA

av

一T

-

-

《&

②

ax

■一

,

应力

单元

求得

可以

方程

物理

根据

e

e

=Sa

DBa

Ds=

(y=

cr=

y

变矩

是应

，B

矩阵

应力

称为

S,s

5

=[$

…司

0H

s=

其中

矩阵

弹性

所以

和V,

取为E

和V0

，E0

力问题

平面应

由于是

阵。

0

1V

E

-

D=-

0

v1

1-U

——

00

2)

1

O

节。

第2.2

教材

考了

容参

分内

这部

变换:

等参元

4.3、

曲的单

形状扭

中几何

体坐标

换成总

单元转

规则的

何形状

标中几

然）坐

部（自

将局

为了

:

转换

坐标

建立

必须

要，

化的需

行离散

解域进

形状求

对一般

以满足

元，

（J*

=g

7）y

式,

的形

函数

插值

示成

中表

上式

是将

方法

便的

最方

,'y

ZN

y=

i

N「X

x=Z

/=0

/=0

:

式为