空间桁架结构程序的设计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)