结构有限元分析程序研究生阶段.docx
《结构有限元分析程序研究生阶段.docx》由会员分享，可在线阅读，更多相关《结构有限元分析程序研究生阶段.docx（71页珍藏版）》请在冰豆网上搜索。
结构有限元分析程序研究生阶段
中南大学
CENTRALSOUTHUNIVERSITY
题目有限单元法程序设计
学生姓名陈煜
学号114811114
指导老师周文伟
学院土木工程学院
专业岩土工程
完成时间2011年12月
一、PF程序
1.源程序
C主程序读入并输出控制数据，分配动态数组地址，并通过调用各
C子程序来调用整个计算
CPF.FOR（Aprogramforanalysisofplaneframe）
CVersion4.31994
CMainprogramreadsthecontroldata&organizesthewhole
Ccalculationbycallingsubroutines.
DIMENSIONW（20000）
CHARACTERIDFN*20,TITLE（5）*72
READ（*,'（A12）'）IDFN
OPEN（3,FILE=IDFN,STATUS='OLD'）
READ（3,'（A72）'）（TITLE（M）,M=1,5）
READ（3,*）E,NM,NJ,NS,NLC
L1=1
L2=L1+NM
L3=L2+NM
L4=L3+NM
L11=L4+NM
L12=L11+NJ
L21=L12+NJ
L22=L21+NS
L31=L22+NS
L41=L31+6*NM
CALLIOMJS（TITLE,E,NM,NJ,NS,NLC,W（L1）,W（L2）,W（L3）,
&W（L4）,W（L11）,W（L12）,W（L21）,W（L22））
CALLLCVCT（NM,W（L1）,W（L2）,W（L31）,NJ,N）
CALLLCDIA（NM,N,W（L31）,W（L41）,W（L41）,W（L41）,MAXBDW,NA）
L51=L41+N
L52=L51+36
L53=L52+NA*2
L54=L53
L61=L54+N*2
NW=L61+6*NM-1
WRITE（*,1）NA,NW
1FORMAT（/40X,'（NA=',I6,'）'
&/40X,'（NW=',I6,'）'）
CALLFORMA（E,NM,NJ,N,NA,W（L1）,W（L2）,W（L3）,W（L4）,
&W（L11）,W（L12）,W（L31）,W（L51）,W（L41）,W（L52））
CALLAS（NS,N,NA,W（L21）,W（L41）,W（L52））
CALLLDLT（N,NA,W（L41）,W（L52）,W（L53））
DO100LC=1,NLC
READ（3,*）NLJ
L62=L61+NLJ
L63=L62+NLJ
L64=L63+NLJ
L71=L61
L81=L71+6*NM
CALLB0（LC,N,NLJ,W（L54））
IF（NLJ.NE.0）CALLIOLJB（N,NLJ,W（L61）,W（L62）,
&W（L63）,W（L64）,W（L54））
READ（3,*）NLM
L82=L81+NLM
L83=L82+NLM
L84=L83+NLM
CALLF0（NLM,NM,W（L71））
IF（NLM.NE.0）CALLIOLMFB（NM,NJ,N,NLM,W（L81）,
&W（L82）,W（L83）,W（L84）,W（L1）,W（L2）,W（L11）,
&W（L12）,W（L31）,W（L71）,W（L54））
CALLBS（NS,N,W（L21）,W（L22）,W（L54））
CALLSLVEQ（N,NA,MAXBDW,W（L41）,W（L52）,W（L54））
CALLOJD（NJ,N,W（L54））
CALLCOTF（E,NM,NJ,N,W（L1）,W（L2）,W（L3）,W（L4）,
&W（L11）,W（L12）,W（L31）,W（L54）,W（L71））
NW=L84+NLM-1
WRITE（*,1）NA,NW
100CONTINUE
WRITE（*,'（/）'）
END
SUBROUTINEIOMJS（TITLE,E,NM,NJ,NS,NLC,IST,IEN,
&AR,RI,X,Y,IS,VS）
CReaddataofmembers,joints,supports&printthem
DIMENSIONIST（NM）,IEN（NM）,AR（NM）,RI（NM）,
&X（NJ）,Y（NJ）,IS（NS）,VS（NS）
CHARACTERTITLE（5）*72
WRITE（*,'（/）'）
WRITE（*,1）（TITLE（M）,M=1,5）
1FORMAT（1X,A72）
WRITE（*,2）E,NM,NJ,NS,NLC
2FORMAT（/13X,'TheInputData'
&//10X,'TheGeneralInformation'
&//6X,'E',9X,'NM',5X,'NJ',5X,'NS',5X,'NLC'
&/1X,1PE10.3,4I7）
READ（3,*）（IST（M）,IEN（M）,AR（M）,RI（M）,M=1,NM）
WRITE（*,3）
3FORMAT（/10X,'TheInformationofMembers'
&//1X,'member',2X,'start',2X,'end',9X,'A',15X,'I'）
WRITE（*,4）（M,IST（M）,IEN（M）,AR（M）,RI（M）,M=1,NM）
4FORMAT（1X,I4,I8,I6,1P2E16.6）
READ（3,*）（X（M）,Y（M）,M=1,NJ）
WRITE（*,5）
5FORMAT（/10X,'TheJointCoordinates'
&//1X,'joint',11X,'X',17X,'Y'）
WRITE（*,6）（M,X（M）,Y（M）,M=1,NJ）
6FORMAT（1X,I4,2F18.6）
READ（3,*）（IS（M）,VS（M）,M=1,NS）
WRITE（*,7）
7FORMAT（/10X,'TheInformationofSupports'
&//4X,'IS',9X,'VS'）
WRITE（*,8）（IS（M）,VS（M）,M=1,NS）
8FORMAT（1X,I5,F16.6）
RETURN
END
SUBROUTINELCVCT（NM,IST,IEN,LV,NJ,N）
CDeterminelocationvectorofelement
DIMENSIONIST（NM）,IEN（NM）,LV（6,NM）
DO100M=1,NM
I=IST（M）*3
J=IEN（M）*3
LV（1,M）=I-2
LV（2,M）=I-1
LV（3,M）=I
LV（4,M）=J-2
LV（5,M）=J-1
LV（6,M）=J
100CONTINUE
N=NJ*3
RETURN
END
SUBROUTINELCDIA（NM,N,LV,MIN,IBDW,LD,MAXBDW,NA）
CDeterminelocationofdiagonalelementsofglobalstiffness
CmatrixA
DIMENSIONLV（6,NM）,MIN（N）,IBDW（N）,LD（N）
DO100I=1,N
MIN（I）=I
100CONTINUE
DO400M=1,NM
MINLV=LV（1,M）
DO200I=2,6
IF（LV（I,M）.LT.MINLV）MINLV=LV（I,M）
200CONTINUE
DO300I=1,6
IF（MINLV.LT.MIN（LV（I,M）））MIN（LV（I,M））=MINLV
300CONTINUE
400CONTINUE
MAXBDW=0
DO500I=1,N
IBDW（I）=I-MIN（I）+1
IF（IBDW（I）.GT.MAXBDW）MAXBDW=IBDW（I）
500CONTINUE
LD
（1）=IBDW
（1）
DO600I=2,N
LD（I）=LD（I-1）+IBDW（I）
600CONTINUE
NA=LD（N）
RETURN
END
SUBROUTINERLCS（M,NM,NJ,IST,IEN,X,Y,RL,C,S）
CCalculatelength,cosine&sineofmember
DIMENSIONIST（NM）,IEN（NM）,X（NJ）,Y（NJ）
I=IST（M）
J=IEN（M）
X1=X（J）-X（I）
Y1=Y（J）-Y（I）
RL=SQRT（X1*X1+Y1*Y1）
C=X1/RL
S=Y1/RL
RETURN
END
SUBROUTINEKEBAR（M,E,NM,NJ,IST,IEN,AR,RI,
&X,Y,C,S,E1,E2,E3,E4）
CCalculateelementstiffnessmatrixalonglocalaxes
DIMENSIONIST（NM）,IEN（NM）,X（NJ）,Y（NJ）,AR（NM）,RI（NM）
CALLRLCS（M,NM,NJ,IST,IEN,X,Y,RL,C,S）
E1=E*AR（M）/RL
E2=12.0*E*RI（M）/（RL*RL*RL）
E3=0.5*E2*RL
E4=0.6666667*E3*RL
RETURN
END
SUBROUTINEKE（M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE）
CCalculateelementstiffnessmatrixalongglobalaxes
DIMENSIONIST（NM）,IEN（NM）,AR（NM）,RI（NM）,
&X（NJ）,Y（NJ）,AE（6,6）
CALLKEBAR（M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4）
A1=E1*C*C+E2*S*S
A2=（E1-E2）*C*S
A3=E1*S*S+E2*C*C
A4=E3*S
A5=E3*C
A6=E4
AE（1,1）=A1
AE（2,1）=A2
AE（2,2）=A3
AE（3,1）=-A4
AE（3,2）=A5
AE（3,3）=A6
AE（4,1）=-A1
AE（4,2）=-A2
AE（4,3）=A4
AE（4,4）=A1
AE（5,1）=-A2
AE（5,2）=-A3
AE（5,3）=-A5
AE（5,4）=A2
AE（5,5）=A3
AE（6,1）=-A4
AE（6,2）=A5
AE（6,3）=0.5*A6
AE（6,4）=A4
AE（6,5）=-A5
AE（6,6）=A6
RETURN
END
SUBROUTINEFORMA（E,NM,NJ,N,NA,IST,IEN,AR,RI,
&X,Y,LV,AE,LD,A）
CFormglobalstiffnessmatrixA
DIMENSIONIST（NM）,IEN（NM）,AR（NM）,RI（NM）,X（NJ）,Y（NJ）,
&LV（6,NM）,AE（6,6）,LD（N）
DOUBLEPRECISIONA（NA）
DO300M=1,NM
CALLKE（M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE）
DO200I=1,6
DO100J=1,I
IF（LV（I,M）.GE.LV（J,M））THEN
A（LD（LV（I,M））-LV（I,M）+LV（J,M））
&=A（LD（LV（I,M））-LV（I,M）+LV（J,M））+AE（I,J）
ELSE
A（LD（LV（J,M））-LV（J,M）+LV（I,M））
&=A（LD（LV（J,M））-LV（J,M）+LV（I,M））+AE（I,J）
ENDIF
100CONTINUE
200CONTINUE
300CONTINUE
RETURN
END
SUBROUTINEAS（NS,N,NA,IS,LD,A）
CIntroducesupportconditionsintoglobalstiffnessmatrixA
DIMENSIONIS（NS）,LD（N）
DOUBLEPRECISIONA（NA）
DO100M=1,NS
I=3*（IS（M）/10）-3+MOD（IS（M）,10）
A（LD（I））=1D22
100CONTINUE
RETURN
END
SUBROUTINELDLT（N,NA,LD,A,T）
CSolveequations
（1）-decompositionofmatrixA
DIMENSIONLD（N）
DOUBLEPRECISIONA（NA）,T（N）,SUM
DO300I=2,N
LDI=LD（I）
I1=I-LDI+LD（I-1）+1
DO200J=I1,I-1
LDJ=LD（J）
J1=J-LDJ+LD（J-1）+1
IF（I1.GT.J1）J1=I1
SUM=0.0D0
DO100K=J1,J-1
SUM=SUM+T（K）*A（LDJ-J+K）
100CONTINUE
T（J）=A（LDI-I+J）-SUM
A（LDI-I+J）=T（J）/A（LDJ）
A（LDI）=A（LDI）-T（J）*A（LDI-I+J）
200CONTINUE
300CONTINUE
RETURN
END
SUBROUTINESLVEQ（N,NA,MAXBDW,LD,A,B）
CSolveequations
（2）-forward&backsubstitution
DIMENSIONLD（N）
DOUBLEPRECISIONA（NA）,B（N）
DO200I=2,N
LDI=LD（I）
I1=I-LDI+LD（I-1）+1
DO100J=I1,I-1
B（I）=B（I）-A（LDI-I+J）*B（J）
100CONTINUE
200CONTINUE
DO300I=1,N
B（I）=B（I）/A（LD（I））
300CONTINUE
DO500I=N-1,1,-1
IMIN=I+MAXBDW
IF（IMIN.GT.N）IMIN=N
DO400J=I+1,IMIN
LDJ=LD（J）
J1=J-LDJ+LD（J-1）+1
IF（I.GE.J1）B（I）=B（I）-A（LDJ-J+I）*B（J）
400CONTINUE
500CONTINUE
RETURN
END
SUBROUTINEB0（LC,N,NLJ,B）
CInitializejointloadvectorB
DOUBLEPRECISIONB（N）
WRITE（*,1）LC,NLJ
1FORMAT（/15X,'LoadingCase',I3
&//10X,'TheLoadingsatJoints'
&//17X,'NLJ=',I4）
DO100I=1,N
B（I）=0.0D0
100CONTINUE
RETURN
END
SUBROUTINEIOLJB（N,NLJ,ILJ,PX,PY,PM,B）
CReaddataofloadingatjoint&formjointloadvectorB
DIMENSIONILJ（NLJ）,PX（NLJ）,PY（NLJ）,PM（NLJ）
DOUBLEPRECISIONB（N）
READ（3,*）（ILJ（M）,PX（M）,PY（M）,PM（M）,M=1,NLJ）
WRITE（*,1）
1FORMAT（/2X,'ILJ',11X,'PX',14X,'PY',15X,'PM'）
WRITE（*,2）（ILJ（M）,PX（M）,PY（M）,PM（M）,M=1,NLJ）
2FORMAT（1X,I4,2F16.4,F18.5）
DO100M=1,NLJ
I=ILJ（M）*3
B（I-2）=PX（M）
B（I-1）=PY（M）
B（I）=PM（M）
100CONTINUE
RETURN
END
SUBROUTINEF0（NLM,NM,F）
CInitializeterminalforcesofmembers
DIMENSIONF（6,NM）
WRITE（*,1）NLM
1FORMAT（/10X,'TheLoadingsatMembers'
&//17X,'NLM=',I4）
DO200J=1,NM
DO100I=1,6
F（I,J）=0.0
100CONTINUE
200CONTINUE
RETURN
END
SUBROUTINEIOLMFB（NM,NJ,N,NLM,ILM,ITL,PV,DST,
&IST,IEN,X,Y,LV,F,B）
CReaddataofloadingatmember&calculatefixed-endforces,
CaddequivalentjointloadstovectorB
DIMENSIONILM（NLM）,ITL（NLM）,PV（NLM）,DST（NLM）,IST（NM）,
&IEN（NM）,X（NJ）,Y（NJ）,LV（6,NM）,F（6,NM）
DOUBLEPRECISIONB（N）
READ（3,*）（ILM（M）,ITL（M）,PV（M）,DST（M）,M=1,NLM）
WRITE（*,1）
1FORMAT（/2X,'ILM',2X,'ITL',11X,'PV',12X,'DST'）
WRITE（*,2）（ILM（M）,ITL（M）,PV（M）,DST（M）,M=1,NLM）
2FORMAT（1X,I4,I5,F16.4,F16.6）
DO100M=1,NLM
L=ILM（M）
CALLRLCS（L,NM,NJ,IST,IEN,X,Y,RL,C,S）
D1=DST（M）
D2=RL-D1
IF（ITL（M）.EQ.1.OR.ITL（M）.EQ.3）THEN
P1=PV（M）*C
P2=-PV（M）*S
ENDIF
IF（ITL（M）.EQ.2.OR.ITL（M）.EQ.4）THEN
P1=PV（M）*S
P2=PV（M）*C
ENDIF
IF（ITL（M）.EQ.1.OR.ITL（M）.EQ.2）THEN
F1=-P1*D2/RL
F4=-P1-F1
F2=-P2*D2*D2*（RL+2.0*D1）/（RL*RL*RL）
F5=-P2-F2
F3=-P2*D1*D2*D2/（RL*RL）
F6=P2*D1*D1*D2/（RL*RL）
ENDIF
IF（ITL（M）.EQ.3.OR.ITL（M）.EQ.4）THEN
G=P2*D1*D1/（12.0*RL*RL）
F3=-G*（（6.0*RL-8.0*D1）*RL+3.0*D1*D1）
F6=G*D1*（4.0*RL-3.0*D1）
F5=-6.0*G*D1*（2.0-D1/RL）
F2=-P2*D1-F5
F4=-P1*D1*D1/（2.0*RL）
F1=-P1*D1-F4
ENDIF
F（1,L）=F（1,L）+F1
F（2,L）=F（2,L）+F2
F（3,L）=F（3,L）+F3
F（4,L）=F（4,L）+F4
F（5,L）=F（5,L）+F5
F（6,L）=F（6,L）+F6
B（LV（1,L））=B（LV（1,L））-F1*C+F2*S
B（LV（2,L））=B（LV（2,L））-F1*S-F2*C
B（LV（3,L））=B（LV（3,L））-F3
B（LV（4,L））=B（LV（4,L））-F4*C+F5*S
B（LV（5,L））=B（LV（5,L））-F4*S-F5*C
B（LV（6,L））=B（LV（6,L））-F6
100CONTINUE
RETURN
END
SUBROUTINEBS（NS,N,IS,VS,B）
CIntroducesupportconditionsintojointloadvectorB
DIMENSIONIS（NS）,VS（NS）
DOUBLEPRECISIONB（N）
DO100M=1,NS
I=3*（IS（M）/10）-3+MOD（IS（M）,10）
B（I）=VS（M）*1D22
100CONTINUE
RETURN
END
SUBROUTINEOJD（NJ,N,B）
CPrintjointdisplacements
DOUBLEPRECISIONB（N）
WRITE（*,1）
1FORMAT（/13X,'TheResultsofCalculation'
&//10X,'TheJointDisplacements'
&//1X,'joint',8X,'u',15X,'v',14X,'phi'）
WRITE（*,2）（M,B（3*M-2）,B（3*M-1）,B（3*M）,M=1,NJ）
2FORMAT（1X,I4,1P3E16.6）
RETURN
END
SUBROUTINECOTF（E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F）
CCalculate&printterminalforcesofmembers
DIMENSIONIST（NM）,IEN（NM）,AR（NM）,RI（NM）,X（NJ）,Y（NJ）,
&LV（6,NM）,F（6,NM）
DOUBLEPRECISIONB（N）,U1,U2,U3,U4,U5,U6
WRITE（*,1）
1FORMAT（/10X
升级会员 