空间桁架静力分析程序及算例1、变量及数组说明
2、空间桁架结构有限元分析程序源代码
!主程序（读入文件，调用总计算程序，输出结果）
CHARACTER IDFUT*20,OUTFUT*20
WRITE(*,*) 'Input Data File name:'
READ (*,*)IDFUT
OPEN (11,FILE=IDFUT,STATUS='OLD')
WRITE(*,*) 'Output File name:'
READ (*,*)OUTFUT
OPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')
WRITE(12,*)'*****************************************'
WRITE(12,*)'* Program for Analysis of Space Trusses *'
WRITE(12,*)'* School of Civil Engineering CSU *'
WRITE(12,*)'* 2012.6.25 Designed By MuZhaoxiang *'
WRITE(12,*)'*****************************************'
WRITE(12,*)' '
WRITE(12,*)'*****************************************'
WRITE(12,*)'*************The Input Data****************'
WRITE(12,*)'*****************************************'
WRITE(12,100)
READ(11,*)NF,NP,NE,NM,NR,NCF,ND
WRITE(12,110)NF,NP,NE,NM,NR,NCF,ND
100 FORMAT(6X,'The General Information'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',& 5X,'NCF',5X,'ND')
110 FORMAT(2X,I2,6I7)
NPF=NF*NP
CALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)
END
!********************************************************************
!总计算程序
SUBROUTINE ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)
DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR), NAE(NE),& AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),&
PP(NPF),FF(NPF),SG(NE),SM(NE)
READ(11,*)(X(I),Y(I),Z(I),I=1,NP)
READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)
READ(11,*)(RR(1,J),RR(2,J),J=1,NR)
READ(11,*)(AE(1,J),J=1,2)
WRITE(12,120)
WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)
WRITE(12,130)
WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)
WRITE(12,140)
WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)
WRITE(12,150)
WRITE(12,151)(AE(1,J),J=1,2)
IF(NCF/=0)THEN
READ(11,*)((PF(I,J),I=1,4),J=1,NCF)
WRITE(12,160)
WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)
ENDIF
120 FORMAT(/6X,'The Information of Joints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')
121 FORMAT(1X,I4,3F8.1)
130 FORMAT(/6X,'The Information of Members'/2x,'Member',2X,'START',4X,'END',6X,'NAE')
131 FORMAT(1X,I4,3I8)
140 FORMAT(/6X,'The Information of SUPPORTS'/2x,'Joint',5X,'S')
141 FORMAT(1X,I4,F8.3)
150 FORMAT(/6X,'The Information of Sections'/4x,'E0',8X,'A0')
151 FORMAT(1X,1PE8.2,F8.4)
160 FORMAT(/6X,'The Loading at Joints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')
161 FORMAT(1X,I4,3F8.2)
CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)
CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)
CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)
CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)
ISH=1
CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)
CALL REBACK(CKK,V,MAXA,NN,NWK,NNM)
CALL DISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF)
END
!矩阵转置子程序
SUBROUTINE MAT(M,N,A,B)
DIMENSION A(M,N),B(N,M)
DO I=1,M
DO J=1,N
B(J,I)=A(I,J)
END DO
END DO
RETURN
END
!单元刚度矩阵的形成
SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)
DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),NAE(NE),AE(2,NM) ,AKE(2,2) 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((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2)
NMI=NAE(IE)
E0=AE(1,NMI);A0=AE(2,NMI)
C=E0*A0/BL
AKE(1,1)=C
AKE(1,2)=-C
AKE(2,1)=-C
AKE(2,2)=C
RETURN
END
!单元坐标转换矩阵
SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)
DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)
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)
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
RETURN
END
!生成单元联系数组LMT
SUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT) DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)
NN=0;NNM=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
NC=C !NC=0,提取了整数部分
C=C-NC !C=0.***,例如C=0.111
DO J=1,NF
C=C*10.0 !例如C=1.21
L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位
!上的数字,这里"+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
NNM=NN+1
DO IE=1,NE
DO I=1,ND
NI=ME(I,IE)
DO J=1,NF
LMT((I-1)*NF+J,IE)=IT(J,NI)
ENDDO
ENDDO
ENDDO
END
!二维总刚中对角线元地址数组
SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)
DIMENSION MAXA(NPF),LMT(NDF,NE)
MAXA=0;NWK=0
MAXA(1)=1
DO I=2,NNM
IP=I-1
IG=IP
DO IE=1,NE
DO J=1,NDF
IF(LMT(J,IE).EQ.IP) THEN
DO K=1,NDF
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+1
ENDDO
NWK= MAXA(NNM)-1
RETURN
END
!生成一维存储结构总刚度矩阵
SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM) DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),& MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)
CKK=0
DO 10 IE=1,NE
TAK=0
CALL FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)
CALL FT(IE,NP,NE,X,Y,Z,ME,T)
CALL MAT(2,6,T,TT)
AK=MATMUL(TT,AKE)
TAK=MATMUL(AK,T) !总体坐标系下的单元刚度矩阵
DO 220 I=1,6
DO 220 J=1,6
NI=LMT(I,IE)
NJ=LMT(J,IE)
IF((NJ-NI).GE.0.AND.NI*NJ.GT.0) THEN
CKK(IJ)=CKK(IJ)+TAK(I,J)
ENDIF
220 CONTINUE
10 CONTINUE
RETURN
END
!生成荷载矩阵
SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF) DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF)
V=0
PP=0
DO I=1,NF
DO J=1,NP
DO K=1,NCF
IF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THEN
V(IT(I,J))=PF(I+1,K)
ENDIF
ENDDO
ENDDO
ENDDO
DO K=1,NCF
DO I=1,NP
IF(I.EQ.PF(1,K))THEN
PP(NF*(I-1)+1)=PF(2,K)
PP(NF*(I-1)+2)=PF(3,K)
PP(NF*(I-1)+3)=PF(4,K)
ENDIF
ENDDO
ENDDO
RETURN
END
!对一维结构总刚度矩阵进行矩阵分解（LDLT）
SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM)
DIMENSION A(NWK),MAXA(NNM)
IF(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,210
210 K=N-KH
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,270
270 KK=MIN0(IC,ND)
C=0.0
DO 280 L=1,KK
280 C=C+A(KI+L)*A(KLT+L)
A(KLT)=A(KLT)-C
260 K=K+1
240 K=N
B=0.0
DO 300 KK=KL,KU
K=K-1
KI=MAXA(K)
C=A(KK)/A(KI)
IF(ABS(C).LT.1.0E+07) GOTO 290
WRITE(IOUT,2010) N,C
STOP
290 B=B+C*A(KK)
300 A(KK)=C
A(KN)=A(KN)-B
304 IF(A(KN)) 310,310,200
310 IF(ISH.EQ.0) GOTO 320
IF(A(KN).EQ.0.0) A(KN)=-1.0E-16
GOTO 200
320 WRITE(IOUT,2000) N,A(KN)
STOP
200 CONTINUE
RETURN
2000 FORMAT(//,' Stop-stiffness matrix not positive definite',//,'no positive& pivot for equation',I4,& //,' pivot =',E20.10)
2010 FORMAT(//,' Stop-strum sequence check failed + because of multiplier& growth for column & number',I4,//,' Multiplier = ',E20.8)
END
!回代，求得节点位移
SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM)
DIMENSION A(NWK),V(NN,1),MAXA(NNM)
NIP=1
DO IP=1,NIP
DO 400 N=1,NN
KU=MAXA(N+1)-1
IF(KU-KL) 400,410,410
410 K=N
C=0.0
DO 420 KK=KL,KU
K=K-1
420 C=C+A(KK)*V(K,IP)
V(N,IP)=V(N,IP)-C
400 CONTINUE
DO 480 N=1,NN
K=MAXA(N)
480 V(N,IP)=V(N,IP)/A(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,510
510 K=N
DO 520 KK=KL,KU
K=K-1
520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)
500 N=N-1
ENDDO
RETURN
END
!求解杆件力、支反力和位移
SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y, Z,PP,FF,SG,SM,ME,& NR,RR,NF)
DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),T(2,6),TT(6,2), AE(2,NM),& ME(2,NE),NAE(NE),UE(6),U(2),AKE(2,2),FE1(2),FE(6),FF(NPF),PP(NPF),&
SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)
SG=0;SM=0;FF=0;FF2=0
DO I=1,NP
DO J=1,NF
LAB=IT(J,I)
IF(LAB.EQ.0) THEN
DIST(3*(I-1)+J)=0.0
ELSEIF(https://www.doczj.com/doc/861204965.html,B.LE.NN) THEN
DIST(3*(I-1)+J)=FTOOL(LAB)
ENDIF
ENDDO
ENDDO
N1=ME(1,IE);N2=ME(2,IE)
DO J=1,NF
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(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)
U=MATMUL(T,UE)
FE1=MATMUL(AKE,U)
CALL MAT(2,6,T,TT)
FE=MATMUL(TT,FE1)
DO J=1,NF
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
ISW=NAE(IE)
AO=AE(2,ISW)
SG(IE)=FE1(2)
SM(IE)=FE1(2)/AO
DO I=1,NPF
FF2(I)=FF(I)-PP(I)
ENDDO
ENDDO
DO I=1,NP
DO J=1,NF
LAB=IT(J,I)
IF(LAB.EQ.0)THEN
K=K+1
FL(K)=FF2(3*(I-1)+J)
ENDIF
ENDDO
ENDDO
WRITE(12,*)' '
WRITE(12,*)'****************************************'
WRITE(12,*)'*********The Results of Calculation**********'
WRITE(12,*)'****************************************'
WRITE(12,600)
WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&
DIST(3*I)*1000, I=1,NP)
WRITE(12,620)
WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)
WRITE(12,640)
WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)
600 FORMAT(6X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)') 610 FORMAT(1X,I4,2X,1P3E12.2)
630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)
640 FORMAT(//6X,'The Bearing Force'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')
650 FORMAT(2X,I4,2X,3F10.2)
RETURN
END
3、算例
以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（b）所示。
（a）空间桁架立面图（尺寸：m，荷载：kN）
（b）空间桁架平面图（结点和杆件编号）
（1）输入数据说明
第一部分：控制数据The General Information
NF NP NE NM NR NCF ND
3 13 2
4 1 6 1 2
第二部分：节点坐标数据The Information of Joints
JOINT X(m) Y(m) Z(m)
1 0 -10 50
第三部分：单元信息数据The Information of Members
第四部分：约束信息The Information of Supports
第五部分：单元截面信息The Information of Sections
第六部分：外荷载信息The Loading at Joints
（2）计算结果
第一部分：节点位移结果
第二部分：单元力及应力结果
第三部分：支反力结果
（3）输入及输出文件
①输入文件.txt
3 13 2
4 1 6 1 2
0 -10 50
43.3 -10 25
43.3 -10 -25
0 -10 -50
-43.3 -10 -25
-43.3 -10 25
-12.5 -2 21.65
12.5 -2 21.65
25 -2 0
12.5 -2 -21.65
-12.5 -2 -21.65
-25 -2 0
0 0 0
1 1 7 1
2 1 8 1
3 2 8 1
4 2 9 1
5 3 9 1
6 3 10 1
7 4 10 1
8 4 11 1
9 5 11 1
10 5 12 1
11 6 12 1
12 6 7 1
13 7 8 1
15 9 10 1
16 10 11 1
17 11 12 1
18 7 12 1
19 7 13 1
20 8 13 1
21 9 13 1
22 10 13 1
23 11 13 1
24 12 13 1
1 0.111
2 0.111
3 0.111
4 0.111
5 0.111
6 0.111
210E6 0.04
13 0 -500 0
②输出文件.txt
*****************************************
* Program for Analysis of Space Trusses *
* School of Civil Engineering CSU *
* 2012.6.25 Designed By MuZhaoxiang *
*****************************************
*****************************************
*************The Input Data****************
*****************************************
The General Information
NF NP NE NM NR NCF ND
3 13 2
4 1 6 1 2
The Information of Joints
Joint X Y Z
1 .0 -10.0 50.0
2 43.
3 -10.0 25.0
3 43.3 -10.0 -25.0
4 .0 -10.0 -50.0
5 -43.3 -10.0 -25.0
6 -43.3 -10.0 25.0
7 -12.5 -2.0 21.6
8 12.5 -2.0 21.6
9 25.0 -2.0 .0
10 12.5 -2.0 -21.6
11 -12.5 -2.0 -21.6
12 -25.0 -2.0 .0
13 .0 .0 .0
The Information of Members
Member START END NAE
1 1 7 1
2 1 8 1
3 2 8 1
4 2 9 1
5 3 9 1
6 3 10 1
7 4 10 1
13 7 8 1
14 8 9 1
15 9 10 1
16 10 11 1
17 11 12 1
18 7 12 1
19 7 13 1
20 8 13 1
21 9 13 1
22 10 13 1
23 11 13 1
24 12 13 1
The Information of SUPPORTS
Joint S
1 .111
2 .111
3 .111
4 .111
5 .111
6 .111
The Information of Sections
E0 A0
2.10E+08 .0400
The Loading at Joints
Joint FX FY FZ
13 .00 -500.00 .00
**************************************** *********The Results of Calculation********** **************************************** The Joint Displacement
Joint X(mm) Y(mm) Z(mm)
1 0.00E+00 0.00E+00 0.00E+00
2 0.00E+00 0.00E+00 0.00E+00
3 0.00E+00 0.00E+00 0.00E+00
4 0.00E+00 0.00E+00 0.00E+00
5 0.00E+00 0.00E+00 0.00E+00
6 0.00E+00 0.00E+00 0.00E+00
7 -1.27E+00 3.25E+00 2.19E+00
8 1.27E+00 3.25E+00 2.19E+00
9 2.53E+00 3.25E+00 -5.25E-08
10 1.27E+00 3.25E+00 -2.19E+00
11 -1.27E+00 3.25E+00 -2.19E+00
12 -2.53E+00 3.25E+00 1.93E-07
13 -3.38E-07 -6.75E+01 5.95E-07
The Terminal Forces
Member FN(kN) σ(MPa)
1 -166.66 -4.17
2 -166.66 -4.17
3 -166.66 -4.17
9 -166.66 -4.17
10 -166.65 -4.17
11 -166.66 -4.17
12 -166.66 -4.17
13 851.03 21.28
14 851.01 21.28
15 851.01 21.28
16 851.03 21.28
17 851.01 21.28
18 851.01 21.28
19 -1044.98 -26.12
20 -1044.98 -26.12
21 -1044.98 -26.12
22 -1044.98 -26.12
23 -1044.98 -26.12
24 -1044.98 -26.12
The Bearing Force
Joint X Y Z
1 .00 83.33 -295.30
2 -255.7
3 83.33 -147.65
3 -255.73 83.33 147.65
4 .00 83.33 295.30
5 255.73 83.33 147.65
6 255.73 83.33 -147.65