Fortran平面钢架有限元分析

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

. 1.2 程序结构及流程 .

. 开 始 标题及数组说明 （读入题目序号NO）

NO是否为零?

形成整体刚度矩阵

计算并打印各杆轴力 解方程并打印杆端位移 形成结点载荷

结 束 子程序READ 子程序MKE

子程序MAKE 子程序MR

子程序MF

是 否 读入数据并打印

子程序MULV6 子程序CALM 子程序MK

子程序MULV 子程序TARN

子程序SOLV 子程序MADE

子程序PE .

. 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 输入输出数据 输入：

输出： .

. .

. 2.4 分析结果 理论分析中，因为力计算应用了力法，所以程序所得结果和理论结果一致。而对位移进行理论分析时忽略了剪力的影响，所以理论位移略小于程序所计算的结果。可以看出软件的的正确性很高，但是此软件只适用计算平面杆系结构，不能解决弹性力学问题 . 结点1位移 u v θ 理论值 1.615e-5 1.611*E-5 0.6617e-3 程序值 1.615e-5 1.640*E-5 0.6617e-3 .

. 3 程序源代码 附上完整的程序源代码。 PROGRAM PFAP C ANALYSIS PROGRAM FOR PLANE FRAME REAL K(200,200),KE(6,6),AKE(6,6),X(100),Y(100),AL(100), #EAI(3,100),PJ(100),PF(2,100),R(6,6),P(100),FF(6), #FE(6),D(100),ADE(6),DE(6),RT(6,6),AFE(6),F(3) INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),M(6),JEAI(100),NO OPEN(6,FILE='ht2.TXT') OPEN(8,FILE='ht.txt',STATUS='NEW') 1 READ(6,*)NO IF(NO.EQ.0)STOP WRITE(8,'(/9X,A5,I3,A1)')'(NO=',NO,')' CALL READ(NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF) DO 5 I=1,N P(I) =0.0 DO 5 J=1,N 5 K(I,J)=0.0 DO 10 IE=1,NEL CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL) CALL MR(R,IE,JE,X,Y) . . CALL MAKE(KE,R,AKE) CALL CALM(M,IE,JN,JE) CALL MK(K,AKE,M) 10 CONTINUE DO 20 IP=1,NPF CALL MR(R,JPF(1,IP),JE,X,Y) CALL TRAN(R,RT) CALL PE(FE,IP,JPF,PF,AL) CALL MULV6(RT,FE,AFE) CALL CALM(M,JPF(1,IP),JN,JE) CALL MF(P,AFE,M) 20 CONTINUE DO 30 I=1,NPJ 30 P(JPJ(I))=P(JPJ(I))+PJ(I) CALL SOLV(K,P,D,N) WRITE(8,'(/2(26(1H*),A))')'RESULTS OF CALCULATION' WRITE(8,'(/28X,A)')'NODEL DISPLACEMENT' WRITE(8,40) 40 FORMAT(9X,'NO.N',4X,'X-DISPLACEMENT',2X, #'Y-DISPLACEMENT',3X,'ANG.ROT.(RAD)') DO 60 KK=1,NJ DO 50 II=1,3