Fortran平面钢架有限元分析

1 有限元分析软件的开发

1.1 程序功能

该程序为平面刚架静力分析程序，能针对平面刚架间问题进行有限元计算，计算杆端位移及杆端力大小。程序从磁盘文件中读取单元编号、节点编号及坐标、材料属性、荷载、边界条件等信息；将杆端位移，杆端力等计算结果以磁盘文件的形式输出，采用等带宽二维数组存储整体刚度矩阵并使用高斯消去法进行求解。

1.2 程序结构及流程

1.3 程序的输入与输出

详细介绍输入输出数据的格式。如：数据文件分几个部分，各有几行，分别包含哪些容及其类型、先后次序，等等。

输入，共有九行。第一行：7,13,5,1,2,2。分别为，7个结点，13个自由度，5个单元，1个类型，2个结点荷载，2个非结点荷载。

第二行：1,2,3,0.0,0.0,0,0,,6.0,0.0。分别为：一号结点的位移序号，x方向为1，y方向为2，转角为3，坐标为(0.0,0.0)，因为二号结点固结在地面，所以二号结点的位移序号，x方向为0，y方向为0，转角为0，坐标为(6.0,0.0)。

第三行：4,5,6,0.0,6.0,4,5,7,0.0,6.0。分别为：三号结点的位移序号，x方向为4，y方向为5，转角为6，坐标为(0.0,6.0), 四号结点位移序号x方向和y相同，转角为7，坐标为(,0.0,6.0)。

第四行：8,9,10,6.0,6.0,0,0,11,0.0,12.0.五号结点位移序号，x方向为8，y方向为9，转角为10，坐标为（6.0,6.0）。因为六号结点铰接在地面，所以六号结点的位移序号，x方向和y方向为0，转角为11，坐标为（0.0,12.0）。

第五行：12,0,13,6.0,12.0. 因为七号结点与地面用滑动支座固定，所以七号结点的位移序号，x方向为12，y方向为0，转角为13，坐标（6.0,12.0）.

第六行：1,2,1,1,3,1,4,5,1,3,6,1,5,7,1,分别为，1号和2号结点组成的单元为1号类型。1号和3号结点组成的单元为1号类型，4号和5号结点组成的为1号类型，3号和6号结点组成的单元为1号类型，5号和7号结点组成的单元为1号类型。

第七行：分别为，弹性模量为E＝2×108 kN/m2，截面面积A=0.16m2,惯性矩I=0.002m4。

第八行：1号结点转角方向的集中力偶为-20.0kN，3号结点集中力为10.0KN。

第九行：1号单元，受集中力（集中力型号为3），大小为15.0kN，到始端的距离为3.0。5号单元，受均布力（均布力型号为1），大小为5.0kN，到端点的距离为5.0。

第十行：0为计算终止符。

输出：第一部分为输入的数据。RESULTS OF CALCULATION以下为输出结果，第二部分的第一段为4个结点的x,y方向的位移和转角。第二段为1,2,3号单元的轴力，剪力和弯矩。

1.4 程序求解中遇到的问题

1对实例进行计算时，坐标原点选用不同的点，会导致整个题目的坐标值发生改变，输入的容会有所不同，最后的结果也不相同

2对结点荷载和非结点荷载的正负判断不同，结点荷载的方向和整体坐标有关，非结点荷载方向判断和局部坐标有关。

3在非结点荷载中，均布荷载和集中力到始端的距离判断不同。

2 有限元分析算例

2.1 算例说明

已知图示刚架，各杆的材料及截面均相同，弹性模量E＝2×108 kN/m2,A=0.16m2,惯性矩I=0.002m4，q=5kN/m.，一号单元集中力为15KN,一号结点集中力偶为20KN*M,三号结点集

中力为10KN. 试求刚架的力。节点编号如图

2.2 理论分析

对所选取的力学问题进行理论分析，要有详细的推导过程和计算结果。

1 力计算

对结构进行分析，可以看出1,2,4单元组成的是二次超静定结构，3,5单元是静定结构。因此先对3,5单元组成的结构进行分析。

如上图所示，可以根据x，y方向力平衡，对结点七力矩平衡算得支座反力。再画出其弯矩，剪力轴力图。然后对1,2,4单元组成的结构分析。

用力法解超静定，将结点六的约束解除，加上支座反力x1=1，x2=1.画出M1，M2,MP图。

MP图M1图M2图

然后画出其弯矩，剪力，轴力图

弯矩图剪力图轴力图

2 位移计算

计算结点1位移，x方向加单位力1

其剪力与弯矩图都为零，轴力图为

根据公式:

轴力图图乘Δ1x=（12.92*6*1）/EA=（12.92*6*1）/3*10^7*0.16=1.615e-5

y方向加单位力1，忽略剪力的影响，弯矩图图乘

Δ1y=-1/EI(3*15.05*3/2+1/2*3*2.82*2/3*3+1/2*1.27*(3+1.27/3)*17.87)+

1/EI(1/2*1.73*24.32*(4.27+2*1.73/3))=1.611*E-5（略小于程序结果）

加单位力偶，剪力与轴力图为零，弯矩图为

θ1=1/EI(3*15.05*1+1/2*3*2.82*1+1/2*1.27*1*17.87)-

1/EI(1/2*1.73*24.32*1)=0.6617e-3

计算结点3的位移，x方向加单位力1.

弯矩 轴力 Δ

1X=-1/EI(3*15.05*6+1/2*3*2.82*6+1/2*1.27*6*17.87)+1/EI(1/2*1.73*24.32*6)-1/EI(1/2*2.17*35.05*(3.29+2*2.71/3) )+1/EI(1/2*3.29*42.47*2.71/3)+（12.92*6*1）/EA =-0.675E-2.其它位移同理可得。

2.3 输入输出数据

输入：

输出：