平面弹性力学有限元源程序FORTRAN.docx
《平面弹性力学有限元源程序FORTRAN.docx》由会员分享，可在线阅读，更多相关《平面弹性力学有限元源程序FORTRAN.docx（20页珍藏版）》请在冰豆网上搜索。
平面弹性力学有限元源程序FORTRAN
平面弹性力学有限元源程序（FORTRAN）
说明：
1基本控制参数信息：
NG,NE,MC,NX,NB,EO,VO,DENSITY,T（共计5个整形数，4个实型数）
NG：
结构的结点总数；
NE：
结构的单元总数；
MC：
平面问题的类型，MC=0，为平面应力，MC=1，为平面应变；
NX：
荷载工况数；
NB：
支承位移数；
EO：
材料弹性模量（Pa）；
VO：
材料泊松比； DENSITY：
容重（N/m3） T：
材料厚度（m）；
2打印输出控制参数：
NWA,NEW,NWK,NWP（4个整形数） 等于1时，输出，否则不输出。
3单元结点信息：
（K,（IJK（I,K）,I=1,3）,K=1,NE）（每行4个整形数，共计NE行） K：
单元号； IJK（1,K）：
K单元I结点编号； IJK（2,K）：
K单元J结点编号； IJK（3,K）：
K单元K结点编号；
4 结点坐标信息：
（（K,XY（1,K）,XY（2,K））,K=1,NG）（每行3个整形数，共计NG行） K：
结点号 XY（1,K）：
K结点X坐标； XY（2,K）：
K结点Y坐标；
5支承信息：
（（K,MB（1,K）,MB（2,K）,ZB（K））,K=1,NB）（每行3个整形数，1个实型数，共计NB行） K：
支承位移序号； MB（1,K）：
第K个支承位移所在的结点号； MB（2,K）：
第K个支承位移的坐标方向； ZB（K）：
第K个支承位移的数值；
6按NX荷载工况数输入荷载信息：
每一荷载工况如下 :
（1）NF,NP,NM（3个整型数） NF：
集中荷载个数； NP：
分布荷载个数； NM：
计自重单元数；
（2）若NF≠0，则输入下面数据 K,MF（1,K）,MF（2,K）,ZF（K）（每行3个整形数，1个实型数，共计NF行） K：
集中荷载序号； MF（1,K）：
第K个集中荷载作用的结点号； MF（2,K）：
第K个集中荷载的坐标方向； ZF（K）：
第K个集中荷载的数值；
（3）若NP≠0，则输入下面数据 K,MP（1,K）,MP（2,K）,ZP（K）（每行3个整形数，1个实型数，共计NP行） K：
分布荷载序号； MP（1,K）：
第K个分布荷载作用的结点号； MP（2,K）：
第K个分布荷载的坐标方向； ZP（K）：
第K个分布荷载的数值；
（4）若NM≠0，则输入下面数据 若NM≥NE，则表示计所有单元的自重，不需输入计自重的单元号； 若NM$DEBUG
PROGRAMPLANE
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
ALLOCATABLE:
:
IJK（:
:
）,XY（:
:
）,BCA（:
:
）,SK（:
:
）,STR（:
:
）,MB（:
:
）,ZB（:
）,B（:
）
ALLOCATABLE:
:
DELD（:
:
:
）,TOD（:
:
）,DELST（:
:
:
）,TOST（:
:
）,DELSUP（:
:
）,TOTSUP（:
）
DIMENSIONEK（6,6）
CHARACTERPN*40,FN*12
WRITE（*,'（A）'）'本程序为计算平面问题的有限元程序'
WRITE（*,'（A）'）'特点：
（1）采用三结点三角形单元；'
WRITE（*,'（A）'）'
（2）采用等带宽存贮技术；'
WRITE（*,'（A）'）' （3）采用高斯消元法解线性方程组。
'
WRITE（*,'（/A）'）'输入计算问题名（PN）:
'
READ（*,'（A）'）PN
CALLFNAME（PN,'.DAT',FN）
WRITE（*,'（2A）'）' 输入数据文件名为：
',FN
OPEN（5,FILE=FN,STATUS='OLD'）
CALLFNAME（PN,'.OUT',FN）
WRITE（*,'（/2A）'）'结果输出数据文件名为：
',FN
OPEN（6,FILE=FN,STATUS='UNKNOWN'）
CALLFNAME（PN,'.OU1',FN）
WRITE（*,'（/2A）'）'参数输出数据文件名为：
',FN
OPEN（7,FILE=FN,STATUS='UNKNOWN'）
READ（5,*）NG,NE,MC,NX,NB,EO,VO,DENSITY,T
WRITE（6,120）NG,NE,MC,NX,NB
WRITE（6,130）EO,VO,DENSITY,T
READ（5,*）NWA,NWE,NWK,NWP
NT=2*NG
ALLOCATE（IJK（3,NE）,XY（2,NG）,BCA（7,NE）,STR（3,NE）,MB（2,NB）,ZB（NB）,B（NT））
ALLOCATE（DELD（2,NG,NX）,TOD（2,NG）,DELST（3,NE,NX）,TOST（3,NE）,DELSUP（NB,NX）,TOTSUP（NB））
CALLCLEAR（2,NG,TOD）
CALLCLEAR（3,NE,TOST）
CALLCLEAR1（NB,TOTSUP）
IF（NG.EQ.0）THEN
STOP000
ENDIF
CALLINPUT（NG,NE,NB,IJK,XY,MB,ZB）
CALLCALND（NG,NE,IJK,ND）
ALLOCATE（SK（NT,ND））
IF（MC.EQ.0）THEN
E=EO
V=VO
ELSE
E=EO/（1-VO*VO）
V=VO/（1-VO）
ENDIF
IP=0
NX1=NX
A1=E/（1-V*V）/4.0
A2=0.5*（1-V）
CALLABC（NG,NE,IJK,XY,BCA,NWA）
CALLCLEAR（NT,ND,SK）
DO100K=1,NE
CALLKE（K,NE,T,V,A1,A2,IJK,BCA,EK,NWE）
CALLSUMK（K,NE,NT,ND,IJK,EK,SK）
100 CONTINUE
CALLCHECK（NT,ND,SK）
110 CONTINUE
IP=IP+1
WRITE（6,'（///1X,A,I4）'）"荷载工况=",IP
READ（5,*）NF,NP,NM
DOI=1,NT
B（I）=0.0
ENDDO
IF（NF.GT.0）THEN
CALLPF（NF,NT,B）
ENDIF
IF（NP.GT.0）THEN
CALLPP（NP,NG,NT,T,XY,B）
ENDIF
IF（NM.GT.0）THEN
CALLPG（NM,NE,NT,DENSITY,T,IJK,BCA,B）
ENDIF
DOI=1,NB
I1=2*（MB（1,I）-1）+2-MB（2,I）
DELSUP（I,IP）=-B（I1）
ENDDO
IF（NF.NE.0.OR.NP.NE.0.OR.NM.NE.0）THEN
CALLDBC（NB,ND,NT,NX,NX1,MB,ZB,SK,B）
CALLGAUSS（NT,ND,NX1,NX,SK,B,NWK,NWP）
CALLSTRESS（NE,NT,V,A1,A2,IJK,BCA,B,STR）
CALLTRES（IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST）
CALLSUPT（IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP,
TOTSUP）
CALLOUTPUT（IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,
TOTSUP）
ELSE
WRITE（*,'（2x,A,I4,A）'）"荷载工况=",IP,"没有荷载!
"
WRITE（6,'（2x,A,I4,A）'）"荷载工况=",IP,"没有荷载!
"
ENDIF
NX1=NX1-1
IF（NX1.GT.0）GOTO110
120 FORMAT（1X,'结点总数=',I3,1X,'单元总数=',I3,1X,'问题类型=',&
&I1,1X,'荷载工况数=',I2,1X,'支承位移数=',I2）
130 FORMAT（/1X,'弹性模量=',E10.3,1X,'泊松比=',F5.3,1X,'密度=',E10.3,&
&1X,'厚度=',F6.3）
END
SUBROUTINEGAUSS（NT,ND,NX1,NX,SK,B,NWK,NWP）
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
DIMENSIONSK（NT,ND）,A（NT,NT）,B（NT）
CALLCLEAR（NT,NT,A）
DOI=1,NT
DOJ=1,ND
IF（（I+J-1）.LE.NT）THEN
A（I,I+J-1）=SK（I,J）
ENDIF
ENDDO
ENDDO
IF（NWK.EQ.1.AND.NX1.EQ.NX）THEN
WRITE（7,'（A）'）"结构总刚（列出右上三角元素）"
DOI=1,NT
WRITE（7,100）（A（I,J）,J=1,NT）
ENDDO
ENDIF
IF（NWP.EQ.1）THEN
WRITE（7,'（1X,A,I3,A）'）"第",NX-NX1+1,"荷载工况的荷载列阵"
DOI=1,NT
WRITE（7,'（1x,E11.4）'）B（I）
ENDDO
ENDIF
DO50K=1,NT-1
DO60I=K+1,NT
C1=A（K,I）/A（K,K）
DO70J=I,NT
A（I,J）=A（I,J）-C1*A（K,J）
70 CONTINUE
B（I）=B（I）-C1*B（K）
60 CONTINUE
50 CONTINUE
B（NT）=B（NT）/A（NT,NT）
DO80I=NT-1,1,-1
DO90J=I+1,NT
B（I）=B（I）-A（I,J）*B（J）
90 CONTINUE
B（I）=B（I）/A（I,I）
80 CONTINUE
100 FORMAT（1x,4000E10.3）
END
SUBROUTINECHECK（NT,ND,SK）
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
DIMENSIONSK（NT,ND）
M=0
DOI=1,NT
IF（SK（I,1）.LE.0.0）THEN
M=M+1
ENDIF
ENDDO
IF（M.GT.0）THEN
WRITE（*,*）"总刚矩阵的对角元素为负值!
!
"
STOP2222
ENDIF
END
SUBROUTINESTRESS（NE,NT,V,A1,A2,IJK,BCA,B,STR）
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
DIMENSIOND1（6）,D2（3）,B（NT）,S（3,6）,IJK（3,NE）,BCA（7,NE）,STR（3,NE）
CALLCLEAR（3,NE,STR）
DO10K=1,NE
DOI=1,3
II=IJK（I,K）
D1（2*（I-1）+1）=B（2*（II-1）+1）
D1（2*（I-1）+2）=B（2*（II-1）+2）
ENDDO
CALLCLEAR（3,6,S）
C1=2*A1/BCA（7,K）
DO20I=1,3
S（1,2*（I-1）+1）=C1*BCA（I,K）
S（1,2*（I-1）+2）=C1*V*BCA（I+3,K）
S（2,2*（I-1）+1）=C1*V*BCA（I,K）
S（2,2*（I-1）+2）=C1*BCA（I+3,K）
S（3,2*（I-1）+1）=C1*A2*BCA（I+3,K）
S（3,2*（I-1）+2）=C1*A2*BCA（I,K）
20 CONTINUE
CALLMATMUL（3,6,1,S,D1,D2）
DOI=1,3
STR（I,K）=D2（I）
ENDDO
10 CONTINUE
END
SUBROUTINEOUTPUT（IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP）
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
DIMENSIONMB（2,NB）,DELD（2,NG,NX）,TOD（2,NG）
DIMENSIONDELST（3,NE,NX）,TOST（3,NE）,DELSUP（NB,NX）,TOTSUP（NB）
CHARACTERVE*6
WRITE（6,20）IP
WRITE（6,30）
DOI=1,NG
WRITE（6,40）I,DELD（1,I,IP）,TOD（1,I）,DELD（2,I,IP）,TOD（2,I）
ENDDO
WRITE（6,50）
DOJ=1,NE
WRITE（6,60）J,（DELST（L,J,IP）,TOST（L,J）,L=1,3）
ENDDO
WRITE（6,70）
DOJ=1,NB
IF（MB（2,J）.EQ.1）THEN
VE='x方向'
ELSE
VE='Y方向'
ENDIF
WRITE（6,80）MB（1,J）,VE,DELSUP（J,IP）,TOTSUP（J）
ENDDO
20 FORMAT（/1X,'荷载工况=',I3,'的计算结果'）
30 FORMAT（/,1X,'结点号',5X,'X位移增量',5X,'X累计位移',5x,&
&'Y位移增量',5X,'Y累计位移'）
40 FORMAT（1X,I3,6X,F10.4,4x,F10.4,4X,F10.4,4x,F10.4）
50 FORMAT（/,1X,'单元号',6x,'X应力增量',6x,'X累计应力',6x,&
&'Y应力增量',6x,'Y累计应力',6x,'剪应力增量',6x,&
&'累计剪应力'）
60 FORMAT（2X,I3,6X,E10.3,5X,E10.3,5X,E10.3,5X,E10.3,6X,E10.3,&
&6X,E10.3）
70 FORMAT（/,1X,'支座结点号',6x,'支反力方向',6x,'支反力增量',&
&6x,'支反力累计量'）
80 FORMAT（5X,I3,12X,A,6x,E10.3,8X,E10.3）
END
SUBROUTINEMATMUL（M,N,L,A,B,C）
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
DIMENSIONA（M,N）,B（N,L）,C（M,L）
DO100I=1,M
DO100J=1,L
C（I,J）=0.0
DO100K=1,N
100 C（I,J）=C（I,J）+A（I,K）*B（K,J）
END
SUBROUTINECLEAR（M,N,A）
IMPLICITREAL*8（A-H,O-Z）,INTEGER
升级会员 