Fortran语言编写及有限元结构程序.docx
- 文档编号:28301783
- 上传时间:2023-07-10
- 格式:DOCX
- 页数:16
- 大小:57.34KB
Fortran语言编写及有限元结构程序.docx
《Fortran语言编写及有限元结构程序.docx》由会员分享,可在线阅读,更多相关《Fortran语言编写及有限元结构程序.docx(16页珍藏版)》请在冰豆网上搜索。
Fortran语言编写及有限元结构程序
Fortran语言编写及有限元结构程序
算例一计算简图及结果输出
用平面刚架静力计算程序下图结构的内力。
各杆EA,EI相同。
已知:
计算简图如下:
(1)输入原始数据
控制参数3,5,8,7,1,2(NE,NJ,N,NW,NPJ,NPF)
结点坐标集结点未知量编号
单元杆端结点编号及单元EA、EI
结点荷载
非结点荷载
(2)输出结果
NE=3NJ=5N=8NW=7NPJ=1NPF=2
NODEXYXXYYZZ
10.00000.0000000
20.00004.0000123
30.00004.0000124
44.00004.0000567
54.00000.0000008
ELEMENTNODE-INODE-JEAEI
1120.400000E+070.160000E+05
2340.400000E+070.160000E+05
3540.400000E+070.160000E+05
CODEPX-PY-PM
7.-15.0000
ELEMENTINDAQ
1.2.2.0000-18.0000
2.1.4.0000-25.0000
NODEUVCETA
10.000000E+000.000000E+000.000000E+00
2-0.221743E-02-0.464619E-04-0.139404E-02
3-0.221743E-02-0.464619E-040.357876E-02
4-0.222472E-02-0.535381E-04-0.298554E-02
50.000000E+000.000000E+000.658499E-03
ELEMENTNQM
1N1=46.4619Q1=10.7119M1=-6.8477
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
算例二计算简图及结果输出
用平面刚架静力计算程序下图结构的内力。
已知:
桁架单元的抗拉刚度为
,平面刚架单元的抗拉刚度为已知:
,抗弯刚度为
。
计算简图如下:
(1)输入原始数据
控制参数
(NE,NJ,N,NW,NPJ,NPF)
结点坐标集结点未知量编号
单元杆端结点编号及单元EA、EI
非结点荷载
(2)输出结果
NE=5NJ=4N=8NW=7NPJ=0NPF=1
NODEXYXXYYZZ
10.00000.0000001
24.00000.0000234
34.0000-3.0000560
48.00000.0000708
ELEMENTNODE-INODE-JEAEI
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
ELEMENTINDAQ
1.1.4.0000-20.0000
NODEUVCETA
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
ELEMENTNQM
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
N2=38.0173Q2=0.0000M2=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=',I5,2X,'N=',I5,2X,'NW=',I5,2X,
$'NPJ=',I5,2X,'NPF='I5)
20FORMAT(/7X,'NODE',7X,'X',11X,'Y',12X,'XX',8X,'YY',8X,'ZZ'/
$(1X,I10,2F12.4,3I10))
30FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',11X,'EA',13X,'EI'/
$(1X,3I10,2E15.6))
40FORMAT(/7X,'CODE',7X,'PX-PY-PM'/(1X,F10.0,F15.4))
50FORMAT(/4X,'ELEMENT',7X,'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
CALLSCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALLCTM(SI,CO,T)
CALLESM(M,NE,BL,EA,EI,KD)
CALLEJC(M,NE,NJ,JE,JN,JC)
DO110I=1,6
DO115J=1,6
S=0.00
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,'NODE',10X,'U',14X,'V',11X,'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,'ELEMENT',13X,'N',17X,'Q',17X,'M')
DO205M=1,NE
CALLSCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALLESM(M,NE,BL,EA,EI,KD)
CALLCTM(SI,CO,T)
CALLEJC(M,NE,NJ,JE,JN,JC)
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
CALLEFX(I,NPF,BL,PF,FO)
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=',F12.4,3X,'M1=',F12.4
$/14X,'N2=',F12.4,3X,'Q2=',F12.4,3X,'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
END
C(七)求单元常数
SUBROUTINESCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
DIMENSIONJE(2,NE),X(NJ),Y(NJ)
REAL*8BL,SI,CO,DX,DY
J1=JE(1,M)
J2=JE(2,M)
DX=X(J2)-X(J1)
DY=Y(J2)-Y(J1)
BL=DSQRT(DX*DX+DY*DY)
SI=DY/BL
CO=DX/BL
RETURN
END
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
END
C(九)形成单元坐标转换矩阵
SUBROUTINECTM(SI,CO,T)
REAL*8T(6,6),SI,CO
DO10I=1,6
DO10J=1,6
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)
RETURN
END
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
GOTO100
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)
GOTO100
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)
GOTO100
50FO
(1)=-Q*A*(1.00-0.500*C)
FO(4)=-0.500*Q*C*A
GOTO100
60FO
(1)=-Q*B/BL
FO(4)=-Q*C
GOTO100
70S=B/BL
FO
(2)=-Q*G*(3.00*S+C)
FO(5)=-FO
(2)
S=S*B/BL
FO(3)=-Q*S*A
FO(6)=Q*G*B
100RETURN
END
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Fortran 语言 编写 有限元 结构 程序