空间桁架结构程序设计(Fortran).doc

实用文案
空间桁架静力分析程序及算例
1、变量及数组说明
输入数据
控制数据
NF
单个节点的自由度数
NP
结构离散节点的总数
NE
结构离散单元的总数
NM
结构中单元不同的特征数类的总数
NR
结构受约束节点的总数
NCF
结构受外荷载作用的节点总数
ND
一个单元的节点总数
几何数据
X(NP)
节点X坐标数组
Y(NP)
节点Y坐标数组
Z(NP)
节点Z坐标数组
ME(ND,NE)
单元节点信息存储矩阵
ME(1,NE)存储杆件始端节点号
ME(2,NE)储存杆件末端节点号
RR(2,NR)
结构约束信息矩阵
RR(1NR)存放受有约束的节点号
RR(2,NR)存放节点位移约束情况
单元特征数据
AE(2,IN)
单元特征数类数组
AE(1,IN)单元的弹性模量
AE(2,IN)单元的横截面面积
NAE(NE)
单元特征类信息存储数组
荷载数据
PF(4,NCF)
外荷载信息数组
PF(1,NCF)存放外荷载作用的节点号
PF(2,NCF)存放X方向的外荷载
PF(3,NCF)存放Y方向的外荷载
PF(4,NCF)存放Z方向的外荷载
输出数据
位移
DIST(NPF)
节点位移数组
DIST(NF*I-2)存放I节点X方向的位移
DIST(NF*I-1)存放I节点Y方向的位移
DIST(NF*I) 存放I节点Z方向的位移
力
SG(NE)
单元内力数组
SM(NE)
单元截面应力数组
FL(NF*NR)
支座反力数组
FL(NF*I-2)存放受约束的I节点X方向的反力
FL(NF*I-1)存放受约束的I节点Y方向的反力
FL(NF*I)存放受约束的I节点Z方向的反力

中间变量
NPF=NF*NP
二维总刚度矩阵的最大行数
NDF=ND*NF
一个单元的自由度总数（2*3=6）
IN
单元特征类总数
AKE(2,2)
单元在局部坐标系中的刚度局矩阵
BL
杆件单元长度
T(2,6)
坐标转换矩阵
TAK(6,6)
单元在总体坐标系中的刚度矩阵
IT(NF,NP)
节点联系数组
LMT(NDF,NE)
单元联系数组
MAXA(NPF)
结构二维总刚度矩阵主对角元地址数组
NWK
结构一维总刚度矩阵的总容量
CKK(NWK)
结构一维总刚度矩阵
NN
结构矩阵方程的方程总数（去掉约束）
NNM
NNM=NN+1
V(NN)
已知节点荷载列阵数组，回代完成后为存放结构位移
PP(NPF)
所有节点荷载列阵数组
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
NDF=ND*NF
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
T(1,3)=CZ;T(2,6)=CZ
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
RETURN
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
IJ=MAXA(NJ)+NJ-NI
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).ANDIT(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
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,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
KL=MAXA(N)+1
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(LAB.GT.0.AND.LAB.LE.NN) THEN
DIST(3*(I-1)+J)=FTOOL(LAB)
ENDIF
ENDDO
ENDDO
DO IE=1,NE
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=ITJ,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)
620 FORMAT(//6X,'The Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'σ(MPa)')
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
24
1
6
1
2

第二部分：节点坐标数据The Information of Joints
JOINT
X(m)
Y(m)
Z(m)
1
0
-10
50
2
43.3
-10
25
3
43.3
-10
-25
4
0
-10
-50
5
-43.3
-10
-25
6
-43.3
-10
25
7
-12.5
-2
21.65
8
12.5
-2
21.65
9
25
-2
0
10
12.5
-2
-21.65
11
-12.5
-2
-21.65
12
-25
-2
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
8
4
11
1
9
5
11
1
10
5
12
1
11
6
12
1
12
6
7
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
0.111
2
0.111
3
0.111
4
0.111
5
0.111
6
0.111

第五部分：单元截面信息The Information of Sections
E0(KPa)
A0(m2)
210E6
0.04

第六部分：外荷载信息The Loading at Joints
JOINT
FX(KN)
FY(KN)
FZ(KN)
13
0
-500
0

（2）计算结果
第一部分：节点位移结果
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

第二部分：单元内力及应力结果
MEMBER
FN(kN)
(MPa)
1
-166.66
-4.17
2
-166.66
-4.17
3
-166.66
-4.17
4
-166.65
-4.17
5
-166.66
-4.17
6
-166.66
-4.17
7
-166.66
-4.17
8
-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
第三部分：支反力结果
JOINT
FX(KN)
FY(KN)
FZ(KN)
1
0
83.33
-295.3
2
-255.73
83.33
-147.65
3
-255.73
83.33
147.65
4
0
83.33
295.3
5
255.73
83.33
147.65
6
255.73
83.33
-147.65
（3）输入及输出文件
① 输入文件.txt
3 13 24 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
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
1 0.111
2 0.111
3 0.111
4 0.111
5 0.111
6 0.111
210E6 0.04
13 0