Fortran语言编写的有限元结构程序Word下载.docx
《Fortran语言编写的有限元结构程序Word下载.docx》由会员分享,可在线阅读,更多相关《Fortran语言编写的有限元结构程序Word下载.docx(14页珍藏版)》请在冰豆网上搜索。
N2=-46.4619Q2=7.2881M2=0.0000
2N1=7.2881Q1=46.4619M1=0.0000
N2=-7.2881Q2=53.5381M2=14.1523
3N1=53.5381Q1=7.2881M1=0.0000
N2=-53.5381Q2=-7.2881M2=-29.1523
算例二计算简图及结果输出
桁架单元的抗拉刚度为
,平面刚架单元的抗拉刚度为已知:
,抗弯刚度为
。
控制参数
(NE,NJ,N,NW,NPJ,NPF)
NE=5NJ=4N=8NW=7NPJ=0NPF=1
NODEXYXXYYZZ
10.00000.0000001
24.00000.0000234
34.0000-3.0000560
48.00000.0000708
1120.600000E+070.184000E+06
2240.600000E+070.184000E+06
3310.200000E+070.000000E+00
4320.200000E+070.000000E+00
5340.200000E+070.000000E+00
1.1.4.0000-20.0000
10.000000E+000.000000E+000.312593E-03
2-0.202759E-04-0.253871E-03-0.144928E-03
3-0.202759E-04-0.185440E-030.000000E+00
4-0.405518E-040.000000E+00-0.227378E-04
1N1=30.4138Q1=37.1896M1=0.0000
N2=-30.4138Q2=42.8104M2=11.2415
2N1=30.4138Q1=2.8104M1=-11.2415
N2=-30.4138Q2=-2.8104M2=0.0000
3N1=-38.0173Q1=0.0000M1=0.0000
N2=38.0173Q2=0.0000M2=0.0000
4N1=45.6207Q1=0.0000M1=0.0000
N2=-45.6207Q2=0.0000M2=0.0000
5N1=-38.0173Q1=0.0000M1=0.0000
C主程序
C
(一)输入原始数据
DIMENSIONJE(2,100),JN(3,100),JC(6),EA(100),EI(100),X(100),
$Y(100),PJ(2,50),PF(4,100)
REAL*8KE(6,6),KD(6,6),T(6,6),P(300),KB(200,20),F(6),FO(6),
$D(6),BL,SI,CO,S,C
OPEN(5,FILE='
RPF1.TXT'
)
open(6,file='
jieguo1.dat'
status='
new'
READ(5,*)NE,NJ,N,NW,NPJ,NPF
READ(5,*)(X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)
READ(5,*)((JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)
IF(NPJ.NE.0)READ(5,*)((PJ(I,J),I=1,2),J=1,NPJ)
IF(NPF.NE.0)READ(5,*)((PF(I,J),I=1,4),J=1,NPF)
WRITE(6,10)NE,NJ,N,NW,NPJ,NPF
WRITE(6,20)(J,X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)
WRITE(6,30)(J,(JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)
IF(NPJ.NE.0)WRITE(6,40)((PJ(I,J),I=1,2),J=1,NPJ)
IF(NPF.NE.0)WRITE(6,50)((PF(I,J),I=1,4),J=1,NPF)
10FORMAT(/6X,'
NE='
I5,2X,'
NJ='
N='
NW='
I5,2X,
$'
NPJ='
NPF='
I5)
20FORMAT(/7X,'
NODE'
7X,'
X'
11X,'
Y'
12X,'
XX'
8X,'
YY'
ZZ'
/
$(1X,I10,2F12.4,3I10))
30FORMAT(/4X,'
ELEMENT'
4X,'
NODE-I'
NODE-J'
EA'
13X,'
EI'
$(1X,3I10,2E15.6))
40FORMAT(/7X,'
CODE'
PX-PY-PM'
/(1X,F10.0,F15.4))
50FORMAT(/4X,'
IND'
10X,'
A'
14X,'
Q'
/
$(1X,2F10.0,2F15.4))
C
(二)形成总结点荷载向量
DO55I=1,N
55P(I)=0.00
IF(NPJ.EQ.0)GOTO65
DO60I=1,NPJ
L=PJ(1,I)
60P(L)=PJ(2,I)
65IF(NPF.EQ.0)GOTO90
DO70I=1,NPF
M=PF(1,I)
CALLSCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALLEFX(I,NPF,BL,PF,FO)
CALLCTM(SI,CO,T)
CALLEJC(M,NE,NJ,JE,JN,JC)
DO75L=1,6
S=0.00
DO80K=1,6
80S=S-T(K,L)*FO(K)
F(L)=S
75CONTINUE
DO85J=1,6
L=JC(J)
IF(L.EQ.0)GOTO85
P(L)=P(L)+F(J)
85CONTINUE
70CONTINUE
C(三)形成整体刚度矩阵
90DO95I=1,N
DO100J=1,NW
100KB(I,J)=0.00
95CONTINUE
DO105M=1,NE
CALLESM(M,NE,BL,EA,EI,KD)
DO110I=1,6
DO115J=1,6
DO120L=1,6
DO125K=1,6
125S=S+T(L,I)*KD(L,K)*T(K,J)
120CONTINUE
KE(I,J)=S
115CONTINUE
110CONTINUE
DO130L=1,6
I=JC(L)
IF(I.EQ.0)GOTO130
DO135K=1,6
J=JC(K)
IF(J.EQ.0.OR.J.LT.I)GOTO135
JJ=J-I+1
KB(I,JJ)=KB(I,JJ)+KE(L,K)
135CONTINUE
130CONTINUE
105CONTINUE
C(四)解线性方程组
N1=N-1
DO140K=1,N1
IM=K+NW-1
IF(N.LT.IM)IM=N
I1=K+1
DO145I=I1,IM
L=I-K+1
C=KB(K,L)/KB(K,1)
JM=NW-L+1
DO150J=1,JM
JJ=J+I-K
150KB(I,J)=KB(I,J)-C*KB(K,JJ)
145P(I)=P(I)-C*P(K)
140CONTINUE
P(N)=P(N)/KB(N,1)
DO155K=1,N1
I=N-K
JM=K+1
IF(NW.LT.JM)JM=NW
DO160J=2,JM
L=J+I-1
160P(I)=P(I)-KB(I,J)*P(L)
155P(I)=P(I)/KB(I,1)
WRITE(6,165)
165FORMAT(/7X,'
U'
V'
CETA'
DO170I=1,NJ
DO175J=1,3
D(J)=0.00
L=JN(J,I)
IF(L.EQ.0)GOTO175
D(J)=P(L)
175CONTINUE
WRITE(6,180)I,D
(1),D
(2),D(3)
180FORMAT(1X,I10,3E15.6)
170CONTINUE
C(五)求单元杆端内力
WRITE(6,200)
200FORMAT(/4X,'
N'
17X,'
M'
)
DO205M=1,NE
DO210I=1,6
L=JC(I)
D(I)=0.00
IF(L.EQ.0)GOTO210
D(I)=P(L)
210CONTINUE
DO220I=1,6
F(I)=0.00
DO230J=1,6
DO240K=1,6
240F(I)=F(I)+KD(I,J)*T(J,K)*D(K)
230CONTINUE
220CONTINUE
IF(NPF.EQ.0)GOTO270
DO250I=1,NPF
L=PF(1,I)
IF(M.NE.L)GOTO250
DO260J=1,6
260F(J)=F(J)+FO(J)
250CONTINUE
270WRITE(6,280)M,(F(I),I=1,6)
280FORMAT(/1X,I10,3X,'
N1='
F12.4,3X,'
Q1='
M1='
F12.4
$/14X,'
N2='
Q2='
M2='
F12.4)
205CONTINUE
CLOSE(5)
STOP
END
C子程序
C(六)形成单元定位向量
SUBROUTINEEJC(M,NE,NJ,JE,JN,JC)
DIMENSIONJE(2,NE),JN(3,NJ),JC(6)
J1=JE(1,M)
J2=JE(2,M)
DO10I=1,3
JC(I)=JN(I,J1)
10JC(I+3)=JN(I,J2)
RETURN
C(七)求单元常数
SUBROUTINESCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
DIMENSIONJE(2,NE),X(NJ),Y(NJ)
REAL*8BL,SI,CO,DX,DY
DX=X(J2)-X(J1)
DY=Y(J2)-Y(J1)
BL=DSQRT(DX*DX+DY*DY)
SI=DY/BL
CO=DX/BL
C(八)形成单元刚度矩阵
SUBROUTINEESM(M,NE,BL,EA,EI,KD)
DIMENSIONEA(NE),EI(NE)
REAL*8KD(6,6),BL,S,G,G1,G2,G3
G=EA(M)/BL
G1=2.00*EI(M)/BL
G2=3.00*G1/BL
G3=2.00*G2/BL
DO10I=1,6
DO10J=1,6
10KD(I,J)=0.00
KD(1,1)=G
KD(1,4)=-G
KD(4,4)=G
KD(2,2)=G3
KD(5,5)=G3
KD(2,5)=-G3
KD(2,3)=-G2
KD(2,6)=-G2
KD(3,5)=G2
KD(5,6)=G2
KD(3,3)=2.00*G1
KD(6,6)=2.00*G1
KD(3,6)=G1
DO20I=1,5
I1=I+1
DO30J=I1,6
30KD(J,I)=KD(I,J)
20CONTINUE
RETURN
C(九)形成单元坐标转换矩阵
SUBROUTINECTM(SI,CO,T)
REAL*8T(6,6),SI,CO
10T(I,J)=0.00
T(1,1)=CO
T(1,2)=SI
T(2,1)=-SI
T(2,2)=CO
T(3,3)=1.00
DO20I=1,3
DO20J=1,3
20T(I+3,J+3)=T(I,J)
C(十)形成单元固端力
SUBROUTINEEFX(I,NPF,BL,PF,FO)
DIMENSIONPF(4,NPF)
REAL*8FO(6),A,B,C,G,Q,S,BL
IND=PF(2,I)
A=PF(3,I)
Q=PF(4,I)
C=A/BL
G=C*C
B=BL-A
DO5J=1,6
5FO(J)=0.00
GOTO(10,20,30,40,50,60,70),IND
10S=Q*A*0.50
FO
(2)=-S*(2.00-2.00*G+C*G)
FO(5)=-S*G*(2.00-C)
S=S*A/6.00
FO(3)=S*(6.00-8.00*C+3.00*G)
FO(6)=-S*C*(4.00-3.00*C)
GOTO100
20S=B/BL
FO
(2)=-Q*S*S*(1.00+2.00*C)
FO(5)=-Q*G*(1.00+2.00*S)
FO(3)=Q*S*S*A
FO(6)=-Q*B*G
30S=B/BL
FO
(2)=-6.00*Q*C*S/BL
FO(5)=-FO
(2)
FO(3)=Q*S*(2.00-3.00*S)
FO(6)=Q*C*(2.00-3.00*C)
40S=Q*A*0.2500
FO
(2)=-S*(2.00-3.00*G+1.60*G*C)
FO(5)=-S*G*(3.00-1.600*C)
S=S*A
FO(3)=S*(2.00-3.00*C+1.200*G)/1.500
FO(6)=-S*C*(1.00-0.800*C)
50FO
(1)=-Q*A*(1.00-0.500*C)
FO(4)=-0.500*Q*C*A
60FO
(1)=-Q*B/BL
FO(4)=-Q*C
70S=B/BL
FO
(2)=-Q*G*(3.00*S+C)
S=S*B/BL
FO(3)=-Q*S*A
FO(6)=Q*G*B
100RETURN