平面三角形3节点有限元程序

平面三角形３结点有限元程序

1、程序名：，

2、程序功能

本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。

程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。

3、程序流程及框图

图１-1 程序流程图

图１-2 程序框图

其中，各子程序的功能如下：

INPUT——输入结点坐标、单元信息和材料参数；

MR——形成结点自由度序号矩阵；

FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；

LOAD——形成整体结点荷载列阵F；

OUTPUT——输出结点位移或结点荷载；

TREAT——由于有非零已知位移，对K和F进行处理；

DECOMP——整体劲度矩阵的分解运算；

FOBA——前代、回代求出未知结点位移δ；

ERFAC——计算约束结点的支座反力；

KRS——计算单元劲度矩阵中的子块K rs。

4、输入数据及变量说明

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

在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。数据文件包括如下内容：

⑴总控信息，共一条，9个数据

NP，NE，NM，NR，NI，NL，NG，ND，NC

NP——结点总数；

NE——单元总数；

NM——材料类型总数；

NR——约束结点总数；

NI ——问题类型标识，0为平面应力问题，1为平面应变问题；

NL ——受荷载作用的结点的数目；

NG ——考虑自重作用为1，不计自重为0；

ND ——非零已知位移结点的数目；

NC ——要计算支座约束反力的结点数目。

⑵ 材料信息，共NM 条，每条依次输入

EO ，VO ，W ，t

EO ——弹性模量（kN/m 2）；

VO ——泊松比；

W ——材料容重（kN/m 3）；

t ——单元厚度（m ）。

这些信息都存放在数组AE （4，NM ）中。

⑶ 坐标信息，共NP 条，每条依次输入

IP ，X ，Y

IP ——结点号；

X ，Y ——分别为结点的x 坐标和y 坐标。

坐标信息存放在数组X(2，NP)中。

⑷ 单元信息，共NE 条，每条依次输入

JE ，L ，Io ，Jo ，Mo

JE ——单元号；

L ——为该单元的材料类型号。

Io ，Jo ，Mo ——分别为该单元i ，j ，m 的整体编码。

单元信息存放在数组MEO(4,NE)中。

⑸ 约束信息共NR 条，每条依次输入一个数 ⨯⨯⨯ ×× IP x I y I

IP ——结点号；

I x ，I y ——分别为该结点的约束情况，如果方向受约束时填0，如果自由则填1。

⑹ 荷载信息，共NL 条，每条依次输入y

IP ，F x ，F y

IP ——结点号；

F x ，F y ——分别为该结点的x ，y 方向的荷载分量（kN ）。

结点号存放在数组NF （NL ）中，结点荷载分量存放在数组FV （2，NL ）中。

⑺ 若ND ＞0，输入非零已知位移信息，共ND 条，每条依次输入

IP ，u x ，u y

IP ——结点号；

u x ，u y ——分别为该结点x ，y 方向已知位移分量（m ），若其中某方向为自由，则其相应分量为0。

结点号存放在数NDI （ND ）中，已知位移分量存放在数组DV(2，ND)中。

⑻ 支座反力信息，共NC 条，每条依次输入

IP ，M1，M2，M3，M4

IP ——支座结点号；

M1，M2，M3，M4 ——为与该支座结点相关的单元号，若不足4个，则用0补充。支座结点号存放在数组NCI(NC)中，相关单元号存放在数组NCE(4，NC)中。

以上数据须按如上顺序存放在数据文件中。除此之外，程序中还用到其他一些主要变量

和数组，说明如下：

N ——结构自由度总数；

NH ——按一维存贮的整体劲度矩阵的总容量；

NX ——最大半带宽；

SK(10000)——维存贮的劲度矩阵；

R(1000)——开始存放等效结点荷载，求解方程以后，用来存放结点位移；

B(6)——存放单元应力σx ，σy ，τxy ，σ1，σ2，α；

MA(1000)——主元素序号指标矩阵；

JR(2，500)——结点自由度序号矩阵；

ME(3)——存放单元结点i ，j ，m 的整体编码；

NN(6)——单元结点自由度序号；

BI(3)，CI(3)——单元劲度矩阵计算公式中的b i ，b j ，b m 和c i ，c j ，c m ；

S ——三角形单元的面积；

H 11，H 12，H 21，H 22——单元劲度矩阵中子块K rs 的4个元素。

5、算例

一个正方形弹性体，厚度为1m ，四边受单位均布法向力作用，由于对称性，取其1/4

进行计算，其有限元网格如图2-3所示，设2

7/1042m kn E ⨯⋅=，167.0=μ，不考虑自重。该问题的精确解应力为x σ＝12kN /m ，y σ＝12kN /m ，xy τ＝0。