平面弹性力学有限元源程序FORTRANWord文件下载.docx
- 文档编号:19964221
- 上传时间:2023-01-12
- 格式:DOCX
- 页数:20
- 大小:20.62KB
平面弹性力学有限元源程序FORTRANWord文件下载.docx
《平面弹性力学有限元源程序FORTRANWord文件下载.docx》由会员分享,可在线阅读,更多相关《平面弹性力学有限元源程序FORTRANWord文件下载.docx(20页珍藏版)》请在冰豆网上搜索。
K单元K结点编号;
4
结点坐标信息:
((K,XY(1,K),XY(2,K)),K=1,NG)(每行3个整形数,共计NG行)
结点号
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行)
支承位移序号;
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行)
分布荷载序号;
MP(1,K):
第K个分布荷载作用的结点号;
MP(2,K):
第K个分布荷载的坐标方向;
ZP(K):
第K个分布荷载的数值;
(4)若NM≠0,则输入下面数据
若NM≥NE,则表示计所有单元的自重,不需输入计自重的单元号;
若NM<
NE,则需要输入计自重的单元号;
$DEBUG
PROGRAMPLANE
IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N)
ALLOCATABLE:
IJK(:
:
),XY(:
),BCA(:
),SK(:
),STR(:
),MB(:
),ZB(:
),B(:
)
DELD(:
),TOD(:
),DELST(:
),TOST(:
),DELSUP(:
),TOTSUP(:
DIMENSIONEK(6,6)
CHARACTERPN*40,FN*12
WRITE(*,'
(A)'
)'
本程序为计算平面问题的有限元程序'
特点:
(1)采用三结点三角形单元;
'
(2)采用等带宽存贮技术;
(3)采用高斯消元法解线性方程组。
(/A)'
输入计算问题名(PN):
READ(*,'
)PN
CALLFNAME(PN,'
.DAT'
FN)
(2A)'
输入数据文件名为:
FN
OPEN(5,FILE=FN,STATUS='
OLD'
.OUT'
(/2A)'
结果输出数据文件名为:
'
OPEN(6,FILE=FN,STATUS='
UNKNOWN'
.OU1'
参数输出数据文件名为:
OPEN(7,FILE=FN,STATUS='
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)
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
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)
IF(NP.GT.0)THEN
CALLPP(NP,NG,NT,T,XY,B)
IF(NM.GT.0)THEN
CALLPG(NM,NE,NT,DENSITY,T,IJK,BCA,B)
DOI=1,NB
I1=2*(MB(1,I)-1)+2-MB(2,I)
DELSUP(I,IP)=-B(I1)
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,
WRITE(*,'
(2x,A,I4,A)'
IP,"
没有荷载!
"
NX1=NX1-1
IF(NX1.GT.0)GOTO110
120
FORMAT(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)
DOJ=1,ND
IF((I+J-1).LE.NT)THEN
A(I,I+J-1)=SK(I,J)
ENDIF
ENDDO
IF(NWK.EQ.1.AND.NX1.EQ.NX)THEN
WRITE(7,'
结构总刚(列出右上三角元素)"
WRITE(7,100)(A(I,J),J=1,NT)
IF(NWP.EQ.1)THEN
(1X,A,I3,A)'
第"
NX-NX1+1,"
荷载工况的荷载列阵"
(1x,E11.4)'
)B(I)
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
50
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
B(I)=B(I)/A(I,I)
80
FORMAT(1x,4000E10.3)
SUBROUTINECHECK(NT,ND,SK)
DIMENSIONSK(NT,ND)
M=0
IF(SK(I,1).LE.0.0)THEN
M=M+1
IF(M.GT.0)THEN
WRITE(*,*)"
总刚矩阵的对角元素为负值!
!
STOP2222
SUBROUTINESTRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR)
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)
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
CALLMATMUL(3,6,1,S,D1,D2)
STR(I,K)=D2(I)
10
SUBROUTINEOUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP)
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)
WRITE(6,50)
DOJ=1,NE
WRITE(6,60)J,(DELST(L,J,IP),TOST(L,J),L=1,3)
WRITE(6,70)
DOJ=1,NB
IF(MB(2,J).EQ.1)THEN
VE='
x方向'
Y方向'
WRITE(6,80)MB(1,J),VE,DELSUP(J,IP),TOTSUP(J)
荷载工况='
I3,'
的计算结果'
30
FORMAT(/,1X,'
结点号'
5X,'
X位移增量'
X累计位移'
5x,&
Y位移增量'
Y累计位移'
40
FORMAT(1X,I3,6X,F10.4,4x,F10.4,4X,F10.4,4x,F10.4)
单元号'
6x,'
X应力增量'
X累计应力'
6x,&
Y应力增量'
Y累计应力'
剪应力增量'
累计剪应力'
FORMAT(2X,I3,6X,E10.3,5X,E10.3,5X,E10.3,5X,E10.3,6X,E10.3,&
6X,E10.3)
支座结点号'
支反力方向'
支反力增量'
6x,'
支反力累计量'
FORMAT(5X,I3,12X,A,6x,E10.3,8X,E10.3)
END
SUBROUTINEMATMUL(M,N,L,A,B,C)
DIMENSIONA(M,N),B(N,L),C(M,L)
DO100I=1,M
DO100J=1,L
C(I,J)=0.0
DO100K=1,N
C(I,J)=C(I,J)+A(I,K)*B(K,J)
END
SUBROUTINECLEAR(M,N,A)
IMPLICITREAL*8(A-H,O-Z),INTEGER
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 平面 弹性 力学 有限元 源程序 FORTRAN