今年最后一贴:平面问题矩形单元有限元分析程序

yp51920

二、程序实现及说明
（一）、总体说明：
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向荷载值依次填写

iwtry

不知道这个跟王勖成老师的书中的程序有什么不同的地方？

yp51920

看说明，这是矩形单元

苍山负雪

具体说说程序的编写技巧如何？ 给个源码作参考就再好不过了

yp51920

个人认为程序的技巧全在算法，语言只是它表达和实现而已，当然也需要一定的硬件基础，懂得程序在机器语言上的实现，即在物理器件上的实现，这主要是对编制大型的程序而言，这样有助于编制高效率的程序。
不过本人硬件不行，向大家讨教

yp51920

to:苍山负雪——给个源码作参考就再好不过了

U盘坏了，过几天传

苍山负雪

讲一讲算法也好嘛，我们大家一起学习学习, 等你的代码和讲解哦 ;)

yp51920

to:苍山负雪

   讲解就不敢当了。我对算法的理解很简单：实际问题—>物理模型—>数学模型—>求解方法（我们通常是有限元了）—>求解—>对比实际状况或实验数据—>OK，问题解决；NO，修正模型重复。求解的具体计算方法我就不敢发话了。
   呵呵，感觉自己说了等于没说。大家一起来讨论吧

  苍山负雪，我QQ3：60309753，邀请加入。不过验证身份哦，呵呵。

苍山负雪

收到，不过不经常上QQ。

fenglinxiaoxue

请问能提供源代码么？

yp51920

2.4  该程序与三角形单元程序的对比
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 的程序方式也不一样。

发帖时请禁用笑脸标记，以防某些字符被转换成笑脸   --- 苍山负雪

yp51920

谢谢提醒，呵呵，又加了２分啊，谢谢咯

fenglinxiaoxue

谢谢！！！

gun3614

好东西，收藏！

96180

怎么不能运行？

lja999

恩 好东西 谢谢了 呵呵 在做有限元
向楼主学习

zhdd0815

正在学习用fortran写的有限元程序，下载下来学学……，谢谢

zlweek197

很好的参考资料，对我们学习有限元变成帮助很大，非常感谢楼主~~~~

zhangjiayong

参考一下，正在编三角形单元的有限元的程序，希望能有启发。

tieqiao

谢谢。。。