空间结构内力位移计算

1、 作业四 编制有限元程序求解结构内力位移第一部分 程序设计过程和子程序的说明本例为空间桁架结构有限元分析程序。设计思路为：自然离散桁架结构，确定各节点自由度；为单元和节点编号，输入支撑信息、荷载信息、截面特性；运行程序求得刚度矩阵，继而求得节点位移和杆件内力。程序所需各子程序已给出，在此只需编制主程序然后调用子程序求解即可。1、 主要变量及数组说明程序中要设置许多变量和数组来存放各种数据。在本程序中，变量及数组名称选用习惯中常用的表示方法，同时遵从FORTRAN90语言的隐含规则，即由字母I-N开头的均为整型，否则为实型。程序中的变量及数组说明详见附录源程序变量和数组说明。2、 内力计算及检算

2、程序形成一维存储总刚子程序CONKB、解线性方程组子程序LDLTREBACK、求节点位移及单元内力子程序DISPLS。下面是子程序的介绍。SUBROUTINE FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT)功能：形成IT以及LMT子程序：传入参数：NP,NE ,NR,ME,RR传出参数：IT,LMT,NN,NN1FMAXA(NN1,NE,LMT,MAXA,NWK,NP)功能：形成MAXA数组传入参数：NN1,NE,LMT,NP传出参数：MAXA,NWKSUBROUTINE CONKB(NP,NE,NWK,ME,X,Y,Z,AE,LMT,MAXA,CKK,NN1)功能：根据

3、输入的杆件编号、节点位置、杆件位置信息及截面信息，形成杆件在局部坐标系下刚度矩阵的子程序： SUBROUTINE FKE(IE,NP,NE,X,Y,Z,ME,AE,AKE)传入参数：NF,NP,NE,NR, ME,X,Y,Z,AE传出参数：AKE 由局部坐标系向总体坐标系转换的子程序：SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)传入参数：IE,NP,NE,X,Y,Z,ME传出参数：T矩阵转置子程序：SUBROUTINE MAT(M,N,A,B)传入参数：M,N,A传出参数：B矩阵乘法子程序：SUBROUTINE MUL(A,B,M,N,L,AB)传入参数：A,B,M,N

4、,L传出参数：AB传入参数：NF,NP,NE,NM,NR,ME,X,Y,Z,AE,NAE,NN1,JS传出参数：IT,MAXA,CKK，NWKSUBROUTINE LDLTREBACK(PF,CKK,V,MAXA,NN,NWK,NN1,ISH,IOUT,NCF,NP) 功能：将一维存储的结构刚度矩阵进行LDLT分解子程序 SUBROUTINE LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1) 传入参数：CKK,MAXA,NN,ISH,IOUT,NWK,NN1 传出参数：MAXA,CKK回代求解得节点位移子程序SUBROUTINE REBACK(PF,CKK,V,MAXA,N

5、N,NWK,NN1,NCF,NR)传入参数：PF,CKK,MAXA,NN,NWK,NN1,NCF,NR传出参数：VSUBROUTINE DISPLS(NP,NE,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM) 功能：求单元内力、节点位移及约束反力 传入参数：NP,NE,NM,NN,IT,FTOOL,AE,NAE,X,Y,Z传出参数：DIST,SG,SM,FF第二部分 程序使用说明1、 在使用所述桁架内力计算程序具体步骤如下：（1）绘制结构计算简图，对杆件和节点进行编号；（2）建立整体坐标系，确定各个节点的坐标；（3）填写输入数据表；（4）建立输入数据

6、文件，文件扩展名为“.txt”；（5）运行程序，按屏幕提示语句进行操作；（6）程序运行完毕，查看输出数据文件结果。2、参照附录部分各个变量的含义和说明，填写输入文件（1）输入NP,NE,NR,NM,NCF（2）输入X(M),Y(M),Z(M),M=1,NP（3）输入MM(M),ME(1,M),ME(2,M),M=1,NE（4）输入RR(M),M=1,NR（5）输入PF(1,M),PF(2,M),PF(3,M),PF(4,M),M=1,NCF（6）输入AE(1,1),AE(2,1)第三部分 空间桁架有限元程序1、 源程序PROGRAM MAINOPEN(1,FILE="IN.TXT&q

7、uot;,STATUS='OLD')OPEN(2,FILE="OUT1.TXT",STATUS='UNKNOWN')READ(1,*)NP,NE,NR,NM,NCFCALL LAST(NP,NE,NR,NM,NCF)ENDSUBROUTINE LAST(NP,NE,NR,NM,NCF) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N)DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(2,NE),& RR(2,NR),PF(4,NCF),AE(2,NM),SG(N

8、E),SM(NE),& IT(3,NP),DIST(NP*3),LMT(6,NE),MAXA(1000),V(1000),CKK(1000),& FF(NP*3),PP(NP*3)READ(1,*)(X(M),Y(M),Z(M),M=1,NP)READ(1,*)(MM(I),ME(1,I),ME(2,I),I=1,NE)READ(1,*)(RR(1,M),RR(2,M),M=1,NR)READ(1,*)(PF(1,M),PF(2,M),PF(3,M),PF(4,M),M=1,NCF)READ(1,*)(AE(1,M),AE(2,M),M=1,NM)WRITE(2,'(2

9、0X,A)')'/'WRITE(2,'(20X,A)')'/ program for analysis of space trusses /'WRITE(2,'(20X,A)')'/ school of civil engineering CSU /'WRITE(2,'(20X,A)')'/ 2012.06.23 designed by HuLihua /'WRITE(2,'(20X,A/)')'/'WRITE(2,'(20X,/A/)

10、')'the control data:'WRITE(2,10)'num joints','num members','num supports','num type','num p'10 FORMAT(1X,A,2X,A,2X,A,4X,A,4X,A)WRITE(2,15)NP,NE,NR,NM,NCF15 FORMAT(3X,I4,8X,I4,9X,I4,8X,I5,8X,I5)WRITE(2,'(/A/)')'the information of joints:

11、'WRITE(2,70)'num','x','y','z'70 FORMAT(1X,A,9X,A,9X,A,10X,A)WRITE(2,80)(I,X(I),Y(I),Z(I),I=1,NP)80 FORMAT(1X,I4,3X,E8.2,3X,E8.2,3X,E8.2)WRITE(2,'(/A/)')'the information of members:'WRITE(2,25)'member','start','end'25 FORMAT(

12、1X,A,2X,A,3X,A)WRITE(2,30)(MM(I),ME(1,I),ME(2,I),I=1,NE)!WRITE(*,*)(ME(1,I),I=1,NE)30 FORMAT(1X,I4,3X,I4,3X,I4)WRITE(2,'(/A/)')'the information of supports:'WRITE(2,35)'joint','s'35 FORMAT(1X,A,2X,A)WRITE(2,400)(RR(1,M),RR(2,M),M=1,NR)400 FORMAT(1X,F4.0,2X,F5.3)WRITE(

13、2,'(/A/)')'the information of load:'WRITE(2,900)'member','x','y','z'900 FORMAT(1X,A,2X,A,3X,A,8X,A)WRITE(2,600)(PF(1,M),PF(2,M),PF(3,M),PF(4,M),M=1,NCF)600 FORMAT(1X,F8.0,3X,F8.3,3X,F8.3,3X,F8.3)WRITE(2,'(/A/)')'the information of sections:

14、'WRITE(2,700)'a','e'700 FORMAT(1X,A,2X,A)WRITE(2,800)AE(1,1),AE(2,1)800 FORMAT(1X,1PE8.2,3X,F8.4)CALL FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT)CALL FMAXA(NN1,NE,LMT,MAXA,NWK,NP)!WRITE(*,*)(LMT(1,I),LMT(2,I),LMT(3,I),LMT(4,I),LMT(5,I),LMT(6,I),I=1,26)!WRITE(*,*)(MAXA(I),I=1,NN1)CALL CONK

15、B(NP,NE,NWK,ME,X,Y,Z,AE,LMT,MAXA,CKK,NN1)!WRITE(*,*)(TAK(I,1),TAK(I,2),TAK(I,3),TAK(I,4),TAK(I,5),TAK(I,6),I=1,6)!WRITE(*,*)NWK!WRITE(*,*)(CKK(I),I=1,NWK)!CALL P1(NCF,P,PF,NP)!CALL SUPPORT(NP,NR,JS,IDS,DJS,TK,P)ISH=10CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1)CALL REBACK(CKK,V,MAXA,NN,NWK,NN1,NP,NCF,PF,

16、IT)WRITE(2,'(20X,/A/)')'the result data:'WRITE(2,'(/1X,A/)')'joint displacements'WRITE(2,'(1X,A,6X,A,15X,A,15X,A/)')'joint','U','V','W'!WRITE(2,40)(I,P(3*I-2),P(3*I-1),P(3*I),I=1,NP)!40 FORMAT(I3,4X,E11.4,4X,E11.4,4X,E11.4)!DO I

17、E=1,NECALL DISPLS(NP,NE,NN,IT,V,DIST,AE,X,Y,Z,PP,FF,SG,SM,ME)WRITE(2,40)(I,DIST(3*I-2),DIST(3*I-1),DIST(3*I),I=1,NP)40 FORMAT(I3,4X,E11.4,4X,E11.4,4X,E11.4)!END DOWRITE (2,'(20X,/A/)')'axial force'WRITE(2,'(/1X,A,5X,A/)')'member','N'WRITE(2,50)(M,SG(M),M=1,NE)

18、50 FORMAT(I3,4X,E10.2)WRITE (2,'(20X,/A/)')'stress'WRITE(2,'(/1X,A,5X,A/)')'member',''WRITE(2,60)(M,SM(M),M=1,NE)60 FORMAT(I3,4X,E10.2)RETURNEND!矩阵转置SUBROUTINE MAT(M,N,A,B) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N)DIMENSION A(M,N),B(N,M) DO I=1,M DO J=1

19、,N B(J,I) = A(I,J) END DO END DORETURNEND!形成单元局部坐标系中的刚度矩阵SUBROUTINE FKE(IE,NP,NE,X,Y,Z,ME,AE,AKE) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N)DIMENSION AKE(2,2),X(NP),Y(NP),Z(NP),ME(2,NE),AE(2,1)!WRITE(*,*)ME(1,IE)N1=ME(1,IE)N2=ME(2,IE)X1=X(N1);Y1=Y(N1);Z1=Z(N1)X2=X(N2);Y2=Y(N2);Z2=Z(N2)BL=SQRT

20、(X2-X1)*2+(Y2-Y1)*2+(Z2-Z1)*2)WRITE(*,*)BLE0=AE(1,1)A0=AE(2,1)C=E0*A0/BL!WRITE(*,*)CAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND!局部坐标与总体坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)!WRITE(*,*)NP !WRITE(*,*

21、)ME(1,IE) T=0 N1=ME(1,IE);N2=ME(2,IE) X1=X(N1);Y1=Y(N1);Z1=Z(N1) X2=X(N2);Y2=Y(N2);Z2=Z(N2) !WRITE(*,*)X1,X2 BL=SQRT(X2-X1)*2+(Y2-Y1)*2+(Z2-Z1)*2) CX=(X2-X1)/BL CY=(Y2-Y1)/BL CZ=(Z2-Z1)/BL T(1,1)=CX;T(2,4)=CX T(1,2)=CY;T(2,5)=CY T(1,3)=CZ;T(2,6)=CZ !WRITE(*,*)(T(1,I),T(2,I),I=1,6)RETURNEND!形成IT和LMT数

22、组SUBROUTINE FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N)!This program forms the joint&element numbering matrix IT&LMT DIMENSION IT(3,NP),LMT(6,NE),ME(2,NE),RR(2,NR) NN=0;NN1=0;IT=0;LMT=0 N=0 DO I=1,NP C=0 DO K=1,NR KR=RR(1,K) IF(KR.EQ.I)C=RR(2,K) ENDDO

23、NC=C C=C-NC DO J=1,3 C=C*10.0 L=C+0.1 C=C-L IF(L.EQ.0)THEN N=N+1 IT(J,I)=N ELSE IT(J,I)=0 ENDIF ENDDO ENDDO NN=N NN1=NN+1 DO IE=1,NE DO I=1,2 NI=ME(I,IE) DO J=1,3 LMT(I-1)*3+J,IE)=IT(J,NI) ENDDO ENDDO ENDDO RETURN END !矩阵的乘法 SUBROUTINE MUL(A,B,M,N,L,AB) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I

24、-N) DIMENSION A(M,N),B(N,L),AB(M,L) DO I=1,M DO J=1,L AB(I,J)=0 DO K=1,N AB(I,J)=AB(I,J)+A(I,K)*B(K,J) ENDDO ENDDO ENDDO END !形成MAXA数组 SUBROUTINE FMAXA(NN1,NE,LMT,MAXA,NWK,NP) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N)! This program forms the MDE address matrix MAXA of KDIMENSION MAXA(3*NP),L

25、MT(6,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NN1 IP=I-1 IG=IP DO IE=1,NE DO J=1,6 IF(LMT(J,IE).EQ.IP) THENDO K=1,6 IF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE)ENDDO END IF ENDDO ENDDO MAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK=MAXA(NN1)-1RETURNEND !形成一维存储数组总体刚度矩阵CKK SUBROUTINE CONKB(NP,NE,NWK,ME,X,Y,Z,AE,LMT,

26、MAXA,CKK,NN1) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,1),& LMT(6,NE),ME(2,NE),MAXA(NN1),& AKE(2,2),T(2,6),TT(6,2),TAK(6,6),AK(6,2) !WRITE(*,*)(LMT(1,I),LMT(2,I),LMT(3,I),LMT(4,I),LMT(5,I),LMT(6,I),I=1,26) CKK=0 !WRITE(*,*)(ME(1,IE),IE=1,NE

27、) DO IE=1,NE!WRITE(*,*)NETAK=0K=IE CALL FKE(IE,NP,NE,X,Y,Z,ME,AE,AKE) CALL FT(IE,NP,NE,X,Y,Z,ME,T) !WRITE(*,*)(T(1,I),T(2,I),I=1,6) CALL MAT(2,6,T,TT) CALL MUL(TT,AKE,6,2,2,AK) CALL MUL(AK,T,6,2,6,TAK) !CALL FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT) !CALL FMAXA(NN1,NE,LMT,MAXA,NWK,NP) !WRITE(*,*)LMT(1,IE)

28、DO I=1,6 DO J=1,6 NI=LMT(I,IE) !WRITE(*,*)NI NJ=LMT(J,IE) IF(NJ-NI).GE.0.AND.NI*NJ.GT.0)THEN K=MAXA(NJ)+NJ-NI CKK(K)=CKK(K)+TAK(I,J) ENDIF ENDDOENDDOENDDO RETURN END !解线性方程组得位移值V !SUBROUTINE LDLTREBACK(PF,CKK,V,MAXA,IT,NN,NWK,NN1,ISH,IOUT,NCF,NP) !IMPLICIT REAL*8(A-H,O-Z) !IMPLICIT INTEGER*4(I-N) !D

29、IMENSION CKK(NWK),V(NN,1),MAXA(NN1),PF(4,NCF),IT(3,NP) !CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1) !CALL REBACK(CKK,V,MAXA,NN,NWK,NN1,NP,NCF,PF,IT) !END !将一维存储的结构刚度矩阵进行LDLT分解 SUBROUTINE LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION CKK(NWK),MAXA(NN1) IF(

30、NN.EQ.1) RETURN DO 200 N=1,NN KN=MAXA(N) KL=KN+1 KU=MAXA(N+1)-1 KH=KU-KL IF(KH)304,240,210210 K=N-KH IC=0 KLT=KU DO 260 J=1,KH KLT=KLT-1 IC=IC+1 KI=MAXA(K) ND=MAXA(K+1)-KI-1 IF(ND) 260,260,270270 KK=MIN0(IC,ND) C=0.0 DO 280 L=1,KK280 C=C+CKK(KI+L)*CKK(KLT+L) CKK(KLT)=CKK(KLT)-C260 K=K+1240 K=N B=0.0

31、 DO 300 KK=KL,KU K=K-1 KI=MAXA(K) C=CKK(KK)/CKK(KI) IF(ABS(C).LT.1.0E+07) GOTO 290 WRITE(IOUT,2010) N,C STOP290 B=B+C*CKK(KK)300 CKK(KK)=C CKK(KN)=CKK(KN)-B304 IF(CKK(KN) 310,310,200310 IF(ISH.EQ.0) GOTO 320 IF(CKK(KN).EQ.0.0) CKK(KN)=-1.0E-16 GOTO 200320 WRITE(IOUT,2000) N,CKK(KN) STOP200 CONTINUE

32、RETURN2000FORMAT(/' Stop-stiffness matrix not positive& + definite',/,'nonpositive pivot for equation',& I4,/,' pivot =',E20.10)2010FORMAT(/,' Stop-sturm sequence check failed& + because of multiplier growth for column& + number',I4,/, ' Multiplier

33、 = ',E20.8) END !回代求解子程序SUBROUTINE REBACK(CKK,V,MAXA,NN,NWK,NN1,NP,NCF,PF,IT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION CKK(NWK),V(NN),MAXA(NN1),PF(4,NCF),IT(3,NP) V=0 DO J=1,NCF M=PF(1,J) L=IT(1,M) V(L)=PF(2,J) V(L+1)=PF(3,J) V(L+2)=PF(4,J) ENDDO DO 400 N=1,NN KL=MAXA(N)+1 KU=

34、MAXA(N+1)-1 IF(KU-KL) 400,410,410410 K=N C=0.0 DO 420 KK=KL,KU K=K-1420 C=C+CKK(KK)*V(K) V(N)=V(N)-C400 CONTINUE DO 480 N=1,NN K=MAXA(N)480 V(N)=V(N)/CKK(K) IF(NN.EQ.1)RETURN N=NN DO 500 L=2,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF(KU-KL) 500,510,510510 K=N DO 520 KK=KL,KU K=K-1520 V(K)=V(K)-CKK(KK)*V(N)50

35、0 N=N-1 RETURN END!求单元内力节点位移及约束 SUBROUTINE DISPLS(NP,NE,NN,IT,V,DIST,AE,X,Y,Z,PP,FF,SG,SM,ME)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4(I-N)DIMENSION IT(3,NP),DIST(NP*3),V(NN),T(2,6),X(NP),Y(NP),Z(NP),& TT(6,2),AE(2,1),ME(2,NE),UE(6),U(2),& AKE(2,2),FE1(2),FE(6),FF(NP*3),PP(NP*3),SG(NE),SM(N

36、E) SG=0;SM=0;FF=0 DO I=1,NP DO J=1,3 LAB=IT(J,I) IF(LAB.EQ.0) THEN DIST(3*(I-1)+J)=0.0 ELSEIF(LAB.GT.0.AND.LAB.LE.NN) THEN DIST(3*(I-1)+J)=V(LAB) ENDIF ENDDO ENDDO DO IE=1,NE N1=ME(1,IE);N2=ME(2,IE) DO J=1,3 UE(J)=DIST(3*(N1-1)+J) UE(3+J)=DIST(3*(N2-1)+J) ENDDO CALL FT(IE,NP,NE,X,Y,Z,ME,T) CALL FKE(

37、IE,NP,NE,X,Y,Z,ME,AE,AKE) CALL MUL(T,UE,2,6,1,U) !U=T*UE求局部坐标系中的节点位移 CALL MUL(AKE,U,2,2,1,FE1) !FE1=AKE*U求单元局部坐标系节点力 CALL MAT(2,6,T,TT) CALL MUL(TT,FE1,6,2,1,FE) !FE=TT*FE1 DO J=1,3 FF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J) FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J) ENDDO AO=AE(2,1) SG(IE)=FE1(2) SM(IE)=FE1(2)

38、/AO DO I=1,NP*3 FF(I)=FF(I)-PP(I) ENDDO ENDDO RETURN END第四部分 程序使用算例说明：A=0.002，E=2.1E8输出文件OUT1 / program for analysis of space trusses / school of civil engineering CSU / 2012.06.23 designed by Hulihua /the control data: num joints num members num supports num tyepe num p 12 26 4 1 2the information o

39、f joints: num x y z 1 .00E+00 .00E+00 .00E+00 2 .40E+01 .00E+00 .00E+00 3 .80E+01 .00E+00 .00E+00 4 .00E+00 .40E+01 .00E+00 5 .40E+01 .40E+01 .00E+00 6 .80E+01 .40E+01 .00E+00 7 .00E+00 .00E+00 .40E+01 8 .40E+01 .00E+00 .40E+01 9 .80E+01 .00E+00 .40E+01 10 .00E+00 .40E+01 .40E+01 11 .40E+01 .40E+01

40、.40E+01 12 .80E+01 .40E+01 .40E+01the information of members: member start end 1 1 2 2 2 3 3 1 4 4 2 5 5 3 6 6 4 5 7 5 6 8 7 8 9 8 9 10 7 10 11 8 11 12 9 12 13 10 11 14 11 12 15 1 8 16 3 8 17 4 11 18 6 11 19 1 10 20 3 12 21 1 7 22 2 8 23 3 9 24 4 10 25 5 11 26 6 12the information of supports: joint

41、s 1. .111 3. .111 4. .111 6. .111the information of load: member x y z 10. .000 10.000 .000 12. .000 10.000 .000the information of sections: a e 2.10E+08 .02000the result data: joint displacements joint U V W 1 .0000E+00 .0000E+00 .0000E+00 2 .0000E+00 .0000E+00 .0000E+00 3 .0000E+00 .0000E+00 .0000E+00