平面问题有限元应力分析程序

第六章有限元程序设计中的若干问题

基本步骤：

ⅰ.结构离散化,输入或生成

结点信息-结点坐标

单元信息-单元结点编号

ⅱ.计算单元刚度矩阵,形成总体刚度矩阵,包括计算[]B

ⅲ.形成结点载荷向量

ⅳ.引入约束条件

ⅴ.解线性方程组

ⅵ.求出结点位移

ⅶ.计算单元的应力并输出

§6-1 约束条件的处理

1．对称性与反对称性

(1)对称结构承受对称载荷作用时

(2)对称结构承受反对称载荷作用

2. 约束位移的引入

主元置1法

主元赋大值

§6-2 总刚度矩阵的存贮法

1．半带宽存贮法

2．

一维压缩存贮法

考虑到总体刚度矩阵中各行的带宽并不相等，有时由于结构的几何形状的原因，使总体刚度矩阵某些行的带宽特别大。这种情况下如采用半带宽存贮法，就可能把许多零元素也包含了进去，这对节省计算机的存贮量是很不利的。

一维压缩存贮法是将总体刚度矩阵的下三角形中每一行从第一个非零元素开始按行将元素排成一序列，存放于一维数组)(L SK 中。但是为了确定SK 中的元素在[K]中的行列号，还需要将[K]中各行对角线的元素在伊维数组中的序号存放于另一辅助数组KD （N2）中（N2是总刚度矩阵的阶数）。现举例说明这一存贮法： 设有一系数阵

⎥⎥⎥⎥

⎥

⎥⎥⎦

⎤⎢⎢⎢⎢⎢⎢⎢⎣

⎡----1.30

.00.00.06.00.00.00.04.87

.

10.00.00.00

.00.01.50.00.07.10.01.50.00.00.00.00.00.02.100.03.10

.03

.52.03.12.05.4 在一维数组SK （13）中依次存放的是

[]1.3007.16.04.81.52.1003.13.52.05.4--

而辅助数组KD （6）中存放的是

[]1398631

KD （6）其实就是[K]中对角元素在一维数组SK （13）中的地址。

将一结构离散化后，对结点进行编号，就能依据单元号确定出总刚度矩阵[K]各行的带宽，由它依次累加就可得出其对角线元素一维存贮中的序号。

显然，形成了数组Kd ，就确定了[K]中被存贮的元素分布情况以及SK 和[K]中元素的对应关系，例如可求出[K]中第I 行带宽为

)1()(--=I KD I KD L i

也可确定出[K]中第I 行左边第一个非零元素在[K]中的列号

1)1()(+-+-=I KD I KD I M i

此外，也能立即确定出单元刚度矩阵e K ][中的某子矩阵

[]

⎥

⎦

⎤

⎢⎣⎡=⨯)4()3()2()1(2

2I SK I SK I SK I SK K e

ij 组集到一维数组存贮总刚度矩阵SK 中的地址

)(2)12(1J I I KD I -⨯--⨯=，112+=I I 1)(2)12(1--⨯--⨯=J I I KD I ，134+=I I

结构总刚度矩阵的以上两种存贮方法，一般应用于直接解法中。

附录A

平面问题有限元应力分析程序

本程序是用FORTRAN 语言编写的教学使用程序，采用常应变三角形单元，用于解决弹性力学平面应力问题.它极易由学员扩充为求解平面应力和平面应变问题的通用应力分析程序。

程序中总刚度矩阵按一维压缩存贮，线性代数方程组用三角分解直接法求解。

为使教学程序易读易懂，计算时输入计算机的载荷是结点载荷，任何其他载荷都要事先换算成等效的结点载荷.程序所处理的约束仅是X 或Y 轴方向上的零位移约束。

计算结果除输出结点位移外，还输出单元形心处的σx ，σy ， τxy ，主应力σ1,σ2及主方向角θ。

程序的结构便于修改和扩充，易于连接图形库扩充为包含前、后处理程序（网格自动生成及计算结果的图形显示）的更完整的程序系统. 1. 输入数据说明（以输入的先后为序，自由格式）

（1） NN ，NE ，KU ，KV ，KRX ，KRY ，EO ，PO ————6个整型数，2个实型数，

其中

NN 结点总数（≤500）； NE 单元总数（≤700）；

KU x 方向位移受约束的结点数（≤50）； KV y 方向位移受约束的结点数（≤50）；

KRX 在x 方向受结点载荷作用的结点数（≤60）； KRY 在y 方向受结点载荷作用的结点数(≤60)；

EO 材料的弹性模量；

PO 材料的泊桑比；

（2）X（NN）————NN个实型数，为结点的x坐标.

（3）Y（NN）————NN个实型数，为结点的y坐标.

（4）IJM（NE，3）————3*NE个整型数，按行输入，为单元按逆时针向的结点编号（5）JU (KU) ————KU个整型数，为x方向位移受约束的结点号.

（6）JV（KV）————KV个整型数，为y方向位移受约束的结点号.

（7）NRX(KRX) ————KRX个整型数，为在x方向受结点载荷作用的结点号. （8）NRY(KRY) ————KRY个整型数，为在y方向受结点载荷作用的结点号. （9）RX(KRX) ————KRX个实型数，为x方向结点载荷的大小.

（10）RY(KRY) ————KRY个实型数，为y方向结点载荷的大小.

2. 其他标识符说明

NF（=2*NN）方程的总数；

LK 总刚度矩阵下三角形一维存贮的总长；

D1,D2,D3 弹性矩阵中的元素；

B(3),C(3) Bi,Ci(i=1～3)的工作单元；

DEL 三角形单元面积的工作单元；

U（NF）开始存放结点载荷，求解后存放结点位移；

SK(LK) 按一维存贮的总刚度矩阵；

EK(21) 单元刚度矩阵下三角形元素按一维存贮的数组；

BM（3,6）用于存放2Δ〔B〕矩阵的元素；

CM（3,6）用于存放2Δ﹝D﹞〔B〕矩阵的元素；

JD（NF）总刚度矩阵下三角形对角线元素在一维存贮中的序号指示数组；JLL（6）单元刚度矩阵元素和总刚度矩阵元素行列对应关系的工作数组；SGM（3）存放单元形心处的应力σx，σy，τxy；

SGM 1，SGM 2 存放单元形心处的主应力σ1,σ2；

CET A 存放主应力的方向角.

3. 各子程序功能

SKDD 计算总刚度矩阵下三角形对角线元在一维存

贮中的序号

SHAPE（N）用于计算第N个单元的有关常数.

FEK 计算单元刚度矩阵，并存放于EK（21）中.

SKKE 单元刚度矩阵向总刚度矩阵传送.

FIXD 用于处理约束条件.

SOLVE 解线性代数方程组.

STRESS 用于计算单元形心处的应力，主应力和主

向，并输出.

4. 出错信息的处理

本程序中若干数组是半动态数组,如求解问题的规模超过所设

数组下标的界限，则停机，等待处理