程序说明.txt

用于平面四节点单元的静力分析有限元程序

二、程序实现及说明
（一）、总体说明：
1、本程序是平面应力或平面应变问题的矩形单元的静力分析程序，本主要用以较规则的平面结构的有限元分析。
2、本程序要求结构整体坐标水平向右为X轴正向，竖直向上为Y轴正向，矩形单元局部坐标和整体坐标一致。由于结构和矩形单元的规则性，这是很容易做到的。而且要求在编制单元结点号I，J，M，P的规则是从左到右、从下到上，逆时针。这一点也是很容易做到的。
(二)、标识符说明：
NO---算例标号。整型变量，输入参数。
TL(30)—算例标题。字符型数组，输入参数。
NJ---结点总数。整型变量，输入参数。
NE---单元总数。整型变量，输入参数。
NZ---约束结点数。整型变量，输入参数。
NPJ---结点荷载数。整型变量，输入参数。
IPS---问题类型码。整型变量，输入参数。
E---弹性模量。实型变量，输入参数。
PR---泊松比。实型变量，输入参数。
T---弹性体厚度。实型变量，输入参数。
V---容重。实型变量，输入参数。
LND(NE,4)---单元结点编号数组。LND(NE,1)、LND(NE,2)、LND(NE,3)、LND(NE,4)分别为单元结点的单元NE的四个结点编号。整型数组，输入参数。
AA(NE)---单元a值数组。AA(I)为单元号为I的a值。实型数组，输入参数。
BB(NE)---单元b值数组。BB(I)为单元号为I的b值。实型数组，输入参数。
JZ(NZ,3)---支承结点数组。JZ(I,1)为第I个支承结点的结点编号。JZ(I,2)、JZ(I,2)分别是第I个支承结点在位移方向的约束信息。该信息为1时，表示有约束；为0时则无约束。整型数组，输入参数。
PJ(NPJ,3)---结点荷载数组。PJ(I,1)表示第I个结点荷载作用结点的结点编号。PJ(I,2)、PJ(I,3)分别表示该结点沿方向作用的结点荷载数值。实型数组，输入参数。
N---自由度。整型变量，N=NJ×2。
ND---半带宽。整型变量。
NPJ0---存放结点荷载数NPJ信息，以保证当NPJ0=0时，亦可以按规定定义数组。
AK(N,ND)---整体刚度矩阵。
AKE(8,8)---单元刚度矩阵。
B(3,8)--- 位移-应力转换矩阵(几何矩阵)。
D(3,3)---弹性矩阵。
S(3,8)---应力矩阵，S=DB。
P(N)---结点荷载数组，存放结点荷载列向量；解方程后存结点位移。
DE(8)---单元的结点位移。
ST(NE,4)---单元结点应力。ST(I,1)、ST(I,2) 、ST(I,3) 、ST(I,4)分别表示单元I的第I，J，M，P号结点应力。
READ---子程序名，输入数据。
MKE---子程序名，形成单元刚度矩阵。
MD---子程序名，计算弹性矩阵。
MB---子程序名，计算位移-应力转换矩阵。
MF---子程序名，求结点荷载。
RKR---子程序名，引入支承条件，修改刚度方程。
SLOV---子程序名，求结点位移并输出。
MADE---子程序名，计算应力并输出。
(三)、数据输入格式：
算例 
标题 
基本数据 结点总数 单元总数 约束结点数 结点荷载数 问题类型

材料性质 弹性模量 泊松比 单元厚度 容重

单元结点编号 按i，j，m，p结点顺序依次填写

单元a、b值 按单元号及a、b顺序填写

支承结点信息 按约束的结点号、X向信息、Y向信息依次填写

结点荷载信息 按荷载作用结点号、X向荷载值、Y向荷载值依次填写 

注意：该程序与三角形单元程序的对比
1、从所需输入数据看。区别只在于三角形单元需要各结点的坐标值（X(NJ),Y(NJ)），而矩形单元需要的是单元的a、b值（AA(NE),BB(NE)）。
2、从程序的过程看。
（1）、一些数组的阶数有所变化，这是由于三角形是三结点，而四边形是四结点。区别在如下数组：单元结点编号数组、单元刚度矩阵、几何矩阵、应力矩阵S、单元的结点位移、单元(结点)应力，三角形单元分别为LND(NE,3)、AKE(6,6) 、B(3,6) 、S(3,6) 、DE(6) 、ST(3)；矩形单元分别为LND(NE,3)、AKE(8,8)、 B(3,8)、 S(3,8)、 DE(8) 、ST(NE,4)。
（2）、由于矩形单元计算单元面积很方便，所以不需要专门的子程序来计算其单元面积；而三角形单元就不一样，它需要用一个子程序来计算其单元面积（MA子程序）。
（3）、三角形单元是线性单元，从而位移-应力转换矩阵(几何矩阵)B是常量，又应力矩阵S=DB，因D是由材料性质决定，对同一材料而言D是常量，所以S也是常量（均对同一单元各点而言），所以同一单元上的各点应力相同；矩形单元是双线性单元，位移-应力转换矩阵(几何矩阵)B对不同点而言是变量，从而应力矩阵S=DB对不同点也是变量，所以同一单元上各点的应力不同。这就使二者在形成B矩阵、S矩阵的程序方法上不一样，从而导致二者生成单元刚度矩阵K 的程序方式也不一样。