空间桁架结构程序的设计Fortran.docx
《空间桁架结构程序的设计Fortran.docx》由会员分享，可在线阅读，更多相关《空间桁架结构程序的设计Fortran.docx（27页珍藏版）》请在冰豆网上搜索。
空间桁架结构程序的设计Fortran
空间桁架静力分析程序及算例
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（1,NR）存放受有约束的节点号
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、空间桁架结构有限元分析程序源代码
!
主程序（读入文件，调用总计算程序，输出结果）
CHARACTERIDFUT*20,OUTFUT*20
WRITE（*,*）'InputDataFilename:
'
READ（*,*）IDFUT
OPEN（11,FILE=IDFUT,STATUS='OLD'）
WRITE（*,*）'OutputFilename:
'
READ（*,*）OUTFUT
OPEN（12,FILE=OUTFUT,STATUS='UNKNOWN'）
WRITE（12,*）'*****************************************'
WRITE（12,*）'*ProgramforAnalysisofSpaceTrusses*'
WRITE（12,*）'*SchoolofCivilEngineeringCSU*'
WRITE（12,*）'*2012.6.25DesignedByMuZhaoxiang*'
WRITE（12,*）'*****************************************'
WRITE（12,*）''
WRITE（12,*）'*****************************************'
WRITE（12,*）'*************TheInputData****************'
WRITE（12,*）'*****************************************'
WRITE（12,100）
READ（11,*）NF,NP,NE,NM,NR,NCF,ND
WRITE（12,110）NF,NP,NE,NM,NR,NCF,ND
100FORMAT（6X,'TheGeneralInformation'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',&
5X,'NCF',5X,'ND'）
110FORMAT（2X,I2,6I7）
NPF=NF*NP
NDF=ND*NF
CALLANALYSE（NF,NP,NE,NM,NR,NCF,ND,NPF,NDF）
END
!
********************************************************************
!
总计算程序
SUBROUTINEANALYSE（NF,NP,NE,NM,NR,NCF,ND,NPF,NDF）
DIMENSIONX（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
120FORMAT（/6X,'TheInformationofJoints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z'）
121FORMAT（1X,I4,3F8.1）
130FORMAT（/6X,'TheInformationofMembers'/2x,'Member',2X,'START',4X,'END',6X,'NAE'）
131FORMAT（1X,I4,3I8）
140FORMAT（/6X,'TheInformationofSUPPORTS'/2x,'Joint',5X,'S'）
141FORMAT（1X,I4,F8.3）
150FORMAT（/6X,'TheInformationofSections'/4x,'E0',8X,'A0'）
151FORMAT（1X,1PE8.2,F8.4）
160FORMAT（/6X,'TheLoadingatJoints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ'）
161FORMAT（1X,I4,3F8.2）
CALLFLMT（NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT）
CALLFMAXA（NNM,NE,LMT,MAXA,NWK,NPF,NDF）
CALLLP（V,PP,IT,PF,NN,NCF,NF,NP,NPF）
CALLCONKB（NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM）
ISH=1
CALLLDLT（CKK,MAXA,NN,ISH,IOUT,NWK,NNM）
CALLREBACK（CKK,V,MAXA,NN,NWK,NNM）
CALLDISPLS（NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF）
END
!
********************************************************************
!
矩阵转置子程序
SUBROUTINEMAT（M,N,A,B）
DIMENSIONA（M,N）,B（N,M）
DOI=1,M
DOJ=1,N
B（J,I）=A（I,J）
ENDDO
ENDDO
RETURN
END
!
单元刚度矩阵的形成
SUBROUTINEFKE（NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE）
DIMENSIONX（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
!
单元坐标转换矩阵
SUBROUTINEFT（IE,NP,NE,X,Y,Z,ME,T）
DIMENSIONX（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
SUBROUTINEFLMT（NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT）
DIMENSIONIT（NF,NP）,LMT（NDF,NE）,ME（ND,NE）,RR（2,NR）
NN=0;NNM=0;IT=0;LMT=0
N=0
DOI=1,NP
C=0
DOK=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
DOJ=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
DOIE=1,NE
DOI=1,ND
NI=ME（I,IE）
DOJ=1,NF
LMT（（I-1）*NF+J,IE）=IT（J,NI）
ENDDO
ENDDO
ENDDO
RETURN
END
!
二维总刚中对角线元地址数组
SUBROUTINEFMAXA（NNM,NE,LMT,MAXA,NWK,NPF,NDF）
DIMENSIONMAXA（NPF）,LMT（NDF,NE）
MAXA=0;NWK=0
MAXA
（1）=1
DOI=2,NNM
IP=I-1
IG=IP
DOIE=1,NE
DOJ=1,NDF
IF（LMT（J,IE）.EQ.IP）THEN
DOK=1,NDF
IF（LMT（K,IE）.GT.0.AND.LMT（K,IE）.LE.IG）IG=LMT（K,IE）
ENDDO
ENDIF
ENDDO
ENDDO
MAXA（I）=MAXA（I-1）+IP-IG+1
ENDDO
NWK=MAXA（NNM）-1
RETURN
END
!
生成一维存储结构总刚度矩阵
SUBROUTINECONKB（NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM）
DIMENSIONCKK（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
DO10IE=1,NE
TAK=0
CALLFKE（NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE）
CALLFT（IE,NP,NE,X,Y,Z,ME,T）
CALLMAT（2,6,T,TT）
AK=MATMUL（TT,AKE）
TAK=MATMUL（AK,T）!
总体坐标系下的单元刚度矩阵
DO220I=1,6
DO220J=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
220CONTINUE
10CONTINUE
RETURN
END
!
生成荷载矩阵
SUBROUTINELP（V,PP,IT,PF,NN,NCF,NF,NP,NPF）
DIMENSIONV（NN）,PP（NPF）,IT（NF,NP）,PF（4,NCF）
V=0
PP=0
DOI=1,NF
DOJ=1,NP
DOK=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
DOK=1,NCF
DOI=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）
SUBROUTINELDLT（A,MAXA,NN,ISH,IOUT,NWK,NNM）
DIMENSIONA（NWK）,MAXA（NNM）
IF（NN.EQ.1）RETURN
DO200N=1,NN
KN=MAXA（N）
KL=KN+1
KU=MAXA（N+1）-1
KH=KU-KL
IF（KH）304,240,210
210K=N-KH
IC=0
KLT=KU
DO260J=1,KH
KLT=KLT-1
IC=IC+1
KI=MAXA（K）
ND=MAXA（K+1）-KI-1
IF（ND）260,260,270
270KK=MIN0（IC,ND）
C=0.0
DO280L=1,KK
280C=C+A（KI+L）*A（KLT+L）
A（KLT）=A（KLT）-C
260K=K+1
240K=N
B=0.0
DO300KK=KL,KU
K=K-1
KI=MAXA（K）
C=A（KK）/A（KI）
IF（ABS（C）.LT.1.0E+07）GOTO290
WRITE（IOUT,2010）N,C
STOP
290B=B+C*A（KK）
300A（KK）=C
A（KN）=A（KN）-B
304IF（A（KN））310,310,200
310IF（ISH.EQ.0）GOTO320
IF（A（KN）.EQ.0.0）A（KN）=-1.0E-16
GOTO200
320WRITE（IOUT,2000）N,A（KN）
STOP
200CONTINUE
RETURN
2000FORMAT（//,'Stop-stiffnessmatrixnotpositivedefinite',//,'nopositive&pivotforequation',I4,&
//,'pivot=',E20.10）
2010FORMAT（//,'Stop-strumsequencecheckfailed+becauseofmultiplier&growthforcolumn&
number',I4,//,'Multiplier=',E20.8）
END
!
回代，求得节点位移
SUBROUTINEREBACK（A,V,MAXA,NN,NWK,NNM）
DIMENSIONA（NWK）,V（NN,1）,MAXA（NNM）
NIP=1
DOIP=1,NIP
DO400N=1,NN
KL=MAXA（N）+1
KU=MAXA（N+1）-1
IF（KU-KL）400,410,410
410K=N
C=0.0
DO420KK=KL,KU
K=K-1
420C=C+A（KK）*V（K,IP）
V（N,IP）=V（N,IP）-C
400CONTINUE
DO480N=1,NN
K=MAXA（N）
480V（N,IP）=V（N,IP）/A（K）
IF（NN.EQ.1）RETURN
N=NN
DO500L=2,NN
KL=MAXA（N）+1
KU=MAXA（N+1）-1
IF（KU-KL）500,510,510
510K=N
DO520KK=KL,KU
K=K-1
520V（K,IP）=V（K,IP）-A（KK）*V（N,IP）
500N=N-1
ENDDO
RETURN
END
!
求解杆件力、支反力和位移
SUBROUTINEDISPLS（NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,&
NR,RR,NF）
DIMENSIONIT（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
DOI=1,NP
DOJ=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
DOIE=1,NE
N1=ME（1,IE）;N2=ME（2,IE）
DOJ=1,NF
UE（J）=DIST（3*（N1-1）+J）
UE（3+J）=DIST（3*（N2-1）+J）
ENDDO
CALLFT（IE,NP,NE,X,Y,Z,ME,T）
CALLFKE（NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE）
U=MATMUL（T,UE）
FE1=MATMUL（AKE,U）
CALLMAT（2,6,T,TT）
FE=MATMUL（TT,FE1）
DOJ=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
DOI=1,NPF
FF2（I）=FF（I）-PP（I）
ENDDO
ENDDO
DOI=1,NP
DOJ=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,*）'*********TheResultsofCalculation**********'
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）
600FORMAT（6X,'TheJointDisplacement'/2x,'Joint',6X,'X（mm）',8X,'Y（mm）',6X,'Z（mm）'）
610FORMAT（1X,I4,2X,1P3E12.2）
升级会员 