Fortran平面钢架有限元分析.doc

结构力学刚架问题的有限元分析1 有限元分析软件的开发1.1 程序功能该程序为平面刚架静力分析程序，能针对平面刚架间问题进行有限元计算，计算杆端位移及杆端力大小。程序从磁盘文件中读取单元编号、节点编号及坐标、材料属性、荷载、边界条件等信息；将杆端位移，杆端力等计算结果以磁盘文件的形式输出，采用等带宽二维数组存储整体刚度矩阵并使用高斯消去法进行求解。251.2 程序结构及流程开 始标题及数组说明（读入题目序号NO）NO是否为零?形成整体刚度矩阵计算并打印各杆轴力解方程并打印杆端位移形成结点载荷结 束子程序READ子程序MKE子程序MAKE子程序MR子程序MF是否读入数据并打印子程序MULV6子程序CALM子程序MK子程序MULV子程序TARN子程序SOLV子程序MADE子程序PE1.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号类型。第七行：分别为，弹性模量为E2108 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 算例说明已知图示刚架，各杆的材料及截面均相同，弹性模量E2108 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*107*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位移uv理论值1.615e-51.611*E-50.6617e-3程序值1.615e-51.640*E-50.6617e-33 程序源代码附上完整的程序源代码。 PROGRAM PFAPC ANALYSIS PROGRAM FOR PLANE FRAMEREAL 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),NOOPEN(6,FILE=ht2.TXT)OPEN(8,FILE=ht.txt,STATUS=NEW)1READ(6,*)NOIF(NO.EQ.0)STOPWRITE(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,NP(I) =0.0DO 5 J=1,N5K(I,J)=0.0DO 10 IE=1,NELCALL 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 CONTINUEDO 20 IP=1,NPFCALL 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)20CONTINUEDO 30 I=1,NPJ30P(JPJ(I)=P(JPJ(I)+PJ(I)CALL SOLV(K,P,D,N)WRITE(8,(/2(26(1H*),A)RESULTS OF CALCULATIONWRITE(8,(/28X,A)NODEL DISPLACEMENTWRITE(8,40)40FORMAT(9X,NO.N,4X,X-DISPLACEMENT,2X,#Y-DISPLACEMENT,3X,ANG.ROT.(RAD)DO 60 KK=1,NJDO 50 II=1,3F(II)=0.0I1=JN(II,KK)50IF(I1.GT.0)F(II)=D(I1)60WRITE(8,70)KK,F(1),F(2),F(3)70FORMAT(4X,I8,2X,3G16.5)WRITE(8,(/30X,A)ELEMANT FORCESWRITE(8,80)80FORMAT(2X,NO.E,4X,N(1),9X,Q(1),9X,M(1),#9X,N(2),9X,Q(2),9X,M(2)DO 120 IE=1,NELCALL MADE(IE,JN,JE,D,ADE)CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALL MR(R,IE,JE,X,Y)CALL MULV6(R,ADE,DE)CALL MULV6(KE,DE,FF)DO 100 IP=1,NPFIF(JPF(1,IP).EQ.IE)THENCALL PE(FE,IP,JPF,PF,AL)DO 90 I=1,690FF(I)=FF(I)-FE(I) ENDIF100CONTINUEWRITE(8,110)IE,(FF(I),I=1,6)110FORMAT(I5,2X,6G13.7)120CONTINUEGOTO 1ENDSUBROUTINE READ(NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI, #JPJ,PJ,JPF,PF)REAL X(100),Y(100),EAI(3,100),PJ(100),PF(2,100)INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),JEAI(100),#TITLE(20)READ(6,(20A4)(TITLE(I),I=1,20)WRITE(8,(/9X,20A4)TITLEREAD(6,*)NJ,N,NEL,NM,NPJ,NPFWRITE(8,(/3(8X,A4,1H:,I2)NJ=,NJ,#N=,N,NE=,NEL,NM=,NM,NPJ=,NPJ,NPF=,NPFWRITE(8,10)10FORMAT(/8X,NO.N (1) (2) (3),10X,X,9X,Y) READ(6,*)(JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)DO 20 I=1,NJ20WRITE(8,(8X,1H(,I2,1H),3I6,4X,2F10.3)I,JN(1,I),JN(2,I),#JN(3,I),X(I),Y(I)READ(6,*)(JE(1,I),JE(2,I),JEAI(I),I=1,NEL)WRITE(8,25)25FORMAT(/7X,NO.E (1) (2) NO.MAT NO.E (1) (2) NO.MAT) N2=(NEL+1)/2DO 30 I=1,N2-130WRITE(8,40)I,(JE(J,I),J=1,2),JEAI(I), #I+N2,(JE(J,I+N2),J=1,2),JEAI(I+N2)IF(N2*2.NE.NEL)WRITE(8,40)N2,JE(1,N2),JE(2,N2),JEAI(N2)IF(N2*2.EQ.NEL)WRITE(8,40)N2,JE(1,N2),JE(2,N2),JEAI(N2),#NEL,JE(1,NEL),JE(2,NEL),JEAI(NEL)40FORMAT(4X,4I6,2X,4I6) READ(6,*)(EAI(I,J),I=1,3),J=1,NM)WRITE(8,50)(J,(EAI(I,J),I=1,3),J=1,NM)50 FORMAT(/3X,NO.MAT,6X,ELASTIC MODULUS,5X, #AREA,7X,MOMENT OF INERTIA/(I6,5X,3G16.4)IF(NPJ.EQ.0)GOTO 70WRITE(8,(/20X,12H NODEL LOADS)WRITE(8,(16XA)NO.DISP. VALUEREAD(6,*)(JPJ(I),PJ(I),I=1,NPJ)DO 60 I=1,NPJ60WRITE(8,(14X,I7,F16.3)JPJ(I),PJ(I)70 CONTINUEIF(NPF.EQ.0)GOTO 100WRITE(8,(/20X,16H NON-NODEL LOADS)WRITE(8,(7X,A,8X,A,9X,A)NO.E NO.LOAD.MODEL,A,CREAD(6,*)(JPF(1,I),JPF(2,I),PF(1,I),PF(2,I),I=1,NPF)DO 80 I=1,NPF80WRITE(8,90)(JPF(J,I),J=1,2),PF(1,I),PF(2,I)90 FORMAT(6X,I3,8X,I4,8X,2F10.3)100 CONTINUE RETURNENDSUBROUTINE MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)REAL KE(6,6),X(100),Y(100),EAI(3,100),AL(100),LINTEGER JE(2,100),JEAI(100)II=JE(1,IE)JJ=JE(2,IE)MT=JEAI(IE)L=SQRT(X(JJ)-X(II)*2+(Y(JJ)-Y(II)*2)AL(IE)=LA1=EAI(1,MT)*EAI(2,MT)/LA2=EAI(1,MT)*EAI(3,MT)/L*3A3=EAI(1,MT)*EAI(3,MT)/L*2A4=EAI(1,MT)*EAI(3,MT)/LKE(1,1)=A1KE(1,4)=-A1KE(2,2)=12*A2KE(2,3)=6*A3KE(2,5)=-12*A2KE(2,6)=6*A3KE(3,3)=4*A4KE(3,5)=-6*A3KE(3,6)=2*A4 KE(4,4)=A1 KE(5,5)=12*A2 KE(5,6)=-6*A3KE(6,6)=4*A4DO 10 I=1,6DO 10 K=1,610KE(K,I)=KE(I,K) RETURNENDSUBROUTINE MR(R,IE,JE,X,Y)REAL R(6,6),X(100),Y(100),L,CX,CYINTEGER JE(2,100)I=JE(1,IE)J=JE(2,IE)L=SQRT(X(J)-X(I)*2+(Y(J)-Y(I)*2)CX=(X(J)-X(I)/LCY=(Y(J)-Y(I)/LDO 10 J=1,6DO 10 I=1,610R(I,J)=0.0DO 20 I=1,4,3R(I,I)=CXR(I,I+1)=CYR(I+1,I)=-CYR(I+1,I+1)=CX20R(I+2,I+2)=1 RETURNENDSUBROUTINE MAKE(KE,R,AKE)REAL KE(6,6),R(6,6),RT(6,6),TMP(6,6),AKE(6,6)CALL TRAN(R,RT)CALL MULV(RT,KE,TMP)CALL MULV(TMP,R,AKE)RETURNENDSUBROUTINE CALM(M,IE,JN,JE)INTEGER M(6),JN(3,100),JE(2,100),IEDO 10 I=1,3M(I)=JN(I,JE(1,IE)10M(I+3)=JN(I,JE(2,IE)RETURNENDSUBROUTINE MK(K,AKE,M)REAL K(200,200),AKE(6,6)INTEGER M(6)DO 10 I=1,6DO 10 J=1,6IF(M(I).NE.0.AND.M(J).NE.0)#K(M(I),M(J)=K(M(I),M(J)+AKE(I,J)10CONTINUERETURNENDSUBROUTINE PE(FE,IP,JPF,PF,AL)REAL FE(6),PF(2,100),AL(100),LINTEGER JPF(2,100)A=PF(1,IP)C=PF(2,IP)L=AL(JPF(1,IP)IND=JPF(2,IP)DO 5 I=1,65FE(I)=0.0GOTO(10,20,30,40,50,60),IND10FE(2)=(7.*A/20.+3.*C/20.)*LFE(3)=(A/20.+C/30.)*L*2FE(5)=(3.*A/20.+7*C/20.)*LFE(6)=-(A/30.+C/20.)*L*2RETURN20FE(5)=A*C*3*(2.*L-C)/2./L*3FE(2)=A*C-FE(5)FE(3)=A*C*2*(6.*L*L-8.*C*L+3.*C*C)/12./L/LFE(6)=-A*C*3*(4.*L-3.*C)/12./L/LRETURN30FE(2)=A*(L-C)*2*(L+2.*C)/L*3FE(3)=A*C*(L-C)*2/L*2FE(5)=A-FE(2)FE(6)=-A*C*2*(L-C)/L*2RETURN40FE(2)=-6.*A*C*(L-C)/L*3FE(3)=A*(L-C)*(L-3.*C)/L*2FE(5)=-FE(2)FE(6)=A*C*(3.*C-2.*L)/L*2RETURN50FE(1)=A*(1.-C/L)FE(4)=A*C/LRETURN60FE(1)=C*L/2.FE(4)=FE(1)RETURNENDSUBROUTINE MULV6(A,B,C)REAL C(6),A(6,6),B(6)DO 10 I=1,6C(I)=0.0DO 10 J=1,610C(I)=C(I)+A(I,