空间桁架结构程序的设计FortranWord格式文档下载.docx
《空间桁架结构程序的设计FortranWord格式文档下载.docx》由会员分享，可在线阅读，更多相关《空间桁架结构程序的设计FortranWord格式文档下载.docx（27页珍藏版）》请在冰豆网上搜索。
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'
）
OutputFilename:
READ（*,*）OUTFUT
OPEN（12,FILE=OUTFUT,STATUS='
UNKNOWN'
WRITE（12,*）'
*****************************************'
*ProgramforAnalysisofSpaceTrusses*'
*SchoolofCivilEngineeringCSU*'
*2012.6.25DesignedByMuZhaoxiang*'
'
*************TheInputData****************'
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'
NE'
NM'
NR'
&
5X,'
NCF'
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'
X'
Y'
Z'
121FORMAT（1X,I4,3F8.1）
130FORMAT（/6X,'
TheInformationofMembers'
Member'
2X,'
START'
4X,'
END'
6X,'
NAE'
131FORMAT（1X,I4,3I8）
140FORMAT（/6X,'
TheInformationofSUPPORTS'
S'
141FORMAT（1X,I4,F8.3）
150FORMAT（/6X,'
TheInformationofSections'
/4x,'
E0'
8X,'
A0'
151FORMAT（1X,1PE8.2,F8.4）
160FORMAT（/6X,'
TheLoadingatJoints'
FX'
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）
矩阵转置子程序
SUBROUTINEMAT（M,N,A,B）
DIMENSIONA（M,N）,B（N,M）
DOI=1,M
DOJ=1,N
B（J,I）=A（I,J）
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
单元坐标转换矩阵
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）
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
生成单元联系数组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
NN=N
NNM=NN+1
DOIE=1,NE
DOI=1,ND
NI=ME（I,IE）
LMT（（I-1）*NF+J,IE）=IT（J,NI）
二维总刚中对角线元地址数组
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
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）
ENDIF
MAXA（I）=MAXA（I-1）+IP-IG+1
NWK=MAXA（NNM）-1
生成一维存储结构总刚度矩阵
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）
220CONTINUE
10CONTINUE
生成荷载矩阵
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）
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）
对一维结构总刚度矩阵进行矩阵分解（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
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）
200CONTINUE
2000FORMAT（//,'
Stop-stiffnessmatrixnotpositivedefinite'
//,'
nopositive&
pivotforequation'
I4,&
//,'
pivot='
E20.10）
2010FORMAT（//,'
Stop-strumsequencecheckfailed+becauseofmultiplier&
growthforcolumn&
number'
I4,//,'
Multiplier='
E20.8）
回代，求得节点位移
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
IF（KU-KL）400,410,410
410K=N
DO420KK=KL,KU
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
IF（KU-KL）500,510,510
510K=N
DO520KK=KL,KU
520V（K,IP）=V（K,IP）-A（KK）*V（N,IP）
500N=N-1
求解杆件力、支反力和位移
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
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）
UE（J）=DIST（3*（N1-1）+J）
UE（3+J）=DIST（3*（N2-1）+J）
U=MATMUL（T,UE）
FE1=MATMUL（AKE,U）
FE=MATMUL（TT,FE1）
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）
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）
IF（LAB.EQ.0）THEN
K=K+1
FL（K）=FF2（3*（I-1）+J）
WRITE（12,*）'
****************************************'
*********TheResultsofCalculation**********'
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'
X（mm）'
Y（mm）'
Z（mm）'
610FORMAT（1X,I4,2X,1P3E12.2）
升级会员 