《有限单元法基本原理和数值方法》一书的源程序剖析.docx
《《有限单元法基本原理和数值方法》一书的源程序剖析.docx》由会员分享，可在线阅读，更多相关《《有限单元法基本原理和数值方法》一书的源程序剖析.docx（33页珍藏版）》请在冰豆网上搜索。
《有限单元法基本原理和数值方法》一书的源程序剖析
《有限单元法基本原理和数值方法》一书的源程序
--------------------------------------------------------------------------------
!
*********************************************************************
!
*
!
PL----PROGRAMOFPLANEPROBLEM96.1*
!
*
!
*********************************************************************
C
C--------输入数据顺序--------
C1.NG1整型
CNG结构结点总数
CNG=0则停止运行
C2.NE,MC,NX,NB,ND,EO,VO,T5整型3实型
CNE结构单元总数
CMC计算控制类型参数
CMC=0平面应力
C=1平面应变
CNX作用载荷组数
CNB给定位移个数
CND结构刚度矩阵的半带宽
CEO弹性模量
CVO泊松比
CT单元（结构）的厚度
C3.NWANWE,NWK,NWP,NWD5整型
C输出控制参数
C=1输出
C=0不输出
CNWA单元参数的输出控制参数
CNWE单元刚度矩阵的输出控制参数
CNWK结构刚度矩阵的输出控制参数
CNWP载荷向量的输出控制参数
CNWD结点位移的输出控制参数
C4.IJM（3,NE）单元结点编码数组3×NE整型
CIJM（1,I）;IJM（2,I）;IJM（3,I）
C第I个三角形单元的结点编号,按结点编号顺序填写
C5.XY（2,NG）结构结点坐标数组2×NG实型
C6.MB（2,N,ZB（N2×NB整型,NB实型
CMB（1,I）---第I个给定位移所在的结点号
CMB（2,I）=1--给定X方向位移
C=0--给定Y方向位移
CZB（I）----给定位移值（以坐标正向为正）
C7.NF,NP2整型
CNF-----作用在结点上的集中载荷（坐标方向）的个数
CNP-----作用均布侧压的单元边数
C若NF>0则填写
C8.MF（2,NF）,ZF（NF）2×NF整型,NF实型
CMF（1,I）---第I个集中载荷所在的结点号
CMF（2,I）=1--给定X方向集中力
C=0--给定Y方向集中力
CZF（I）-----作用的集中力值
C若NP>0则填写
C9.MP（2,NP）,ZP（NP）2×NP整型,NP实型
CMP（1,I）----第I个载荷作用边的起始结点号
CMP（2,I）----第I个载荷作用边的起始结点号
CZP（I）------第I个均布载荷值
C
C若NX>1重复7.-9.（NX-1）次
C
C最后NG=0表示数据结束
C
!
--------输出数据顺序--------
C1.IJM（3,NE）单元结点编码数组3×NE整型
CIJM（1,I）;IJM（2,I）;IJM（3,I）
C第I个三角形单元的结点编号,按结点编号顺序填写
C2.XY（2,NG）结构结点坐标数组2×NG实型
C
C若NWA=1,则输出
C3.I,B（7）单元参数NE行,1×NE整型,7×NE实型
C每行结构为:
'NE='+单元号+Bi+Bj+Bm+Ci+Cj+Cm+A
C
C若NWE=1,则输出
C4.IO,EK（6×6）单元刚度阵NE行,1×NE整型,6×6×NE实型
C每行结构为:
'NE='+单元号+EK（单元刚度阵）
C
C若NWK=1,则输出
C5.SK（NT,ND）结构刚度矩阵NT×ND=2NG×ND实型
C
C若NWD=1,则输出
C6.I,B结点位移数据NG行,1×NG整型,2×NG实型
C每行结构为:
单元号+U+V
C
C7.S1,S2,S3,X1,X2,CTA
C单元应力数据6×NE实型
C分别代表σx,σy,τxy,σ1,σ2和主应力方向
C
C若NX>1重复6.-7.（NX-1）次
C
!
--------可调数组分配--------
C
C实型数组C（100000）整型数组IA（100000）
CC
（1）XY（2,NG）IA
（1）IJM（3,NE）
CC（N1）ZB（NIA（M1）MB（2,M
CC（N2）BCA（7,NE）IA（M2）MF（2,N
CC（N3）SK（NT,ND）IA（M3）MP（2,NP）
CC（N4）F（NT）IA（MEND）下限
CC（N5）ZF（NF）
CC（N6）ZP（NP）
CC（NEND）下限
C
!
--------程序停止代码--------
!
0正常停止
!
111数组C越界
!
222数组C/IA越界
!
333单元面积非正
!
444结构刚度矩阵主元非正
!
*********************************************************************
C
C主程序
C
DIMENSIONC（500000）,IA（50000）,EK（36）
CHARACTER*12IN,OUT
CIN和OUT为输入文件和输出文件的文件名
WRITE（*,*）''
WRITE（*,*）'PLEASEINPUTTHEINPUT-FILENAME（A<12）'
WRITE（*,*）''
READ（*,5）IN
C输入输入文件的文件名
WRITE（*,*）''
WRITE（*,*）'PLEASEINPUTTHEOUTPUT-FILENAME（A<12）'
WRITE（*,*）''
READ（*,5）OUT
C输入输出文件的文件名
5FORMAT（A12）
OPEN（5,FILE=IN,STATUS='OLD'）
OPEN（6,FILE=OUT,STATUS='UNKNOWN'）
C打开对应的输入和输出文件
10READ（5,*）NG
IF（NG.EQ.0）STOP
C输入结构结点数;如果结点数为0则停止运行
READ（5,*）NE,MC,NX,NB,ND,EO,VO,T
C按顺序输入结构单元数,问题类型参数,载荷组数,给定位移个数
C结构刚度阵的半带宽,弹性模量,泊松比和结构厚度
READ（5,*）NWA,NWE,NWK,NWP,NWD
C按顺序输入各输出控制参数
NT=2*NG
C确定总刚度矩阵阶数NT
C
C计算变界数组的下限
!
M1=3*NE+1
M2=M1+2*NB
N1=2*NG+1
N2=N1+NB
N3=7*NE+N2
N4=N3+NT*ND
N5=N4+NT
C得到各变界数组在一维大数组中的起始元素编号
C
C检验实型数组C的下限
C
NEND=N5
IF（NEND.LE.500000）GOTO35
WRITE（*,*）'***EXCEEDTHELIMITOFARRAYC（INTHEMIDDLE）!
!
***'
WRITE（*,30）NEND
30FORMAT（/,'********NEND=',I6,1X,'>80000********'）
STOP111
C若C下限超出500000,则给出错误信息并停止运行
C
C数据输入
C
35CALLINPUT（NE,NG,NB,IA
（1）,C
（1）,IA（M1）,C（N1））
C调用INPUT子程输入数据
WRITE（*,40）
40FORMAT（/10X,'#####INPUTPASSED#####'）
C显示提示信息
IF（MC.EQ.0）GOTO45
C检验是否是平面应力问题
C
C平面应变问题
C
E=EO/（1.0-VO*VO）
V=VO/（1.0-VO）
C平面应变问题时,先进行弹性常数替换
GOTO50
C
C平面应力问题
C
45E=EO
V=VO
50NX1=NX
A1=E/（1.0-V*V）/4.0
A2=0.5*（1.0-V）
C初始化NX1,A1和A2/NX1为剩余载荷的组数
C
C计算单元参数
C
CALLABC（NE,NG,NWA,IA
（1）,C
（1）,C（N2））
WRITE（*,55）
55FORMAT（/10X,'#####ABCPASSED#####'）
C调用ABC子程计算单元参数并显示提示信息
C
C集成结构刚度矩阵K
C
DO60I=N3,N4
C（I）=0.0
60CONTINUE
C初始化结构刚度矩阵SK
DO65K=1,NE
C遍历结构的所有单元
IO=K
CALLKE（IO,NE,NWE,T,A1,A2,V,EK
（1）,C（N2））
CALLSUMK（IO,NE,ND,NT,IA
（1）,C（N3）,EK
（1））
C调用KE子程计算出单元刚度阵并调用SUMK子程将其集成到结构刚度阵中
65CONTINUE
WRITE（*,70）
70FORMAT（/10X,'#####SUMKPASSED#####'）
C显示提示信息
CALLCHECK（NT,ND,NWK,C（N3））
C调用CHECK子程检验结构刚度阵中的主元是否非正
WRITE（*,75）
75FORMAT（/10X,'#####CHECKPASSED#####'）
C显示提示信息
80READ（5,*）NF,NP
C输入集中载荷个数NF和均布载荷个数NP
C
C再次计算变界数组的下限
C并检验实型数组C和整型数组IA的下限
C
M3=M2+2*NF
N6=N5+NF
NEND=N6+NP-1
MEND=M3+2*NP-1
C计算C和IA的下限
NM=0
IF（NEND.LE.500000）GOTO85
WRITE（*,*）'***EXCEEDTHELIMITOFARRAYC（ATTHEEND）!
!
***'
WRITE（*,30）NEND
NM=1
85IF（MEND.LE.50000）GOTO95
WRITE（*,*）'***EXCEEDTHELIMITOFARRAYIA（ATTHEEND）!
!
***'
WRITE（*,90）MEND
90FORMAT（/,'********MEND=',I6,1X,'>500********'）
STOP222
95IF（NM.EQ.1）STOP222
C检验两个数组的下限,若下限超出则给出错误信息并停止运行
C
C集成结构结点载荷列阵P
C
DO100I=N4,N5
C（I）=0.0
100CONTINUE
C初始化结构结点载荷列阵F
IF（NF.GT.0）CALLPF（NF,NP,NT,NWP,C（N4）,IA（M2）,C（N5））
WRITE（*,105）
105FORMAT（/10X,'#####PFPASSED#####'）
C若集中载荷个数>0,调用PF子程集成各结点集中力并显示提示信息
IF（NP.GT.0）CALLPP（NP,NT,NG,NWP,C
（1）,C（N4）,IA（M3）,C（N6））
WRITE（*,110）
110FORMAT（/10X,'#####PPPASSED#####'）
C若均布载荷个数>0,调用PP子程集成各均布载荷并显示提示信息
C
C引入给定位移
C
CALLDBC（NT,ND,NB,NX,NX1,C（N3）,C（N4）,IA（M1）,C（N1））
WRITE（*,115）
115FORMAT（/10X,'#####DBCPASSED#####'）
C调用DBC子程引入给定位移消除系数矩阵的奇异性,并显示提示信息
C
C求解线性方程组KA=P
C
CALLGAUSS（NT,ND,NWD,NX,NX1,C（N3）,C（N4））
WRITE（*,120）
120FORMAT（/10X,'#####GAUSSPASSED#####'）
C调用GAUSS子程求解线性方程组并显示提示信息
C
C计算单元应力
C
CALLSTRESS（NE,NT,A1,A2,V,IA
（1）,C（N2）,C（N4））
WRITE（*,125）
125FORMAT（/10X,'#####STRESSPASSED#####'）
C调用STRESS子程输出各单元应力并显示提示信息
NX1=NX1-1
IF（NX1.GT.0）GOTO80
C剩余载荷组自减1;若还有载荷剩余则继续计算
GOTO10
C重新输入
END
C*********************************************************************
C1*
C子过程名称:
INPUT*
C子过程功能:
按顺序输入并输出计算所需数据*
C*
C*********************************************************************
SUBROUTINEINPUT（NE,NG,NB,IJM,XY,MB,Z
C形参说明
C输入:
CNE整型,结构单元总数
CNG整型,结构结点总数
CNB整型,给定位移的个数
CIJM（3,NE）整型,单元结点编码数组
CXY（2,NG）实型,结构结点坐标数组
CMB（2,N整型,位移约束信息数组
CZB（N实型,位移约束数值数组
DIMENSIONIJM（3,NE）,XY（2,NG）,MB（2,N,ZB（N
C
READ（5,*）（（IJM（I,L）,I=1,3）,L=1,NE）
C输入单元结点编码数组
READ（5,*）（（XY（I,J）,I=1,2）,J=1,NG）
C输入结构结点坐标数组
READ（5,*）（（MB（I,L）,I=1,2）,L=1,N,（ZB（L）,L=1,N
C输入位移约束信息和数值数组
WRITE（6,20）（（IJM（M,I）,M=1,3）,I=1,NE）
20FORMAT（1X,4（3I4,3X）,3I4）
C输出单元结点编码数组
WRITE（6,40）（（XY（M,I）,M=1,2）,I=1,NG）
40FORMAT（1X,6E12.5）
C输出结构结点坐标数组
RETURN
END
C*********************************************************************
C2*
C子过程名称:
ABC*
C子过程功能:
根据各单元的结点坐标,*
C计算并输出所有单元的各参数*
C*
C*********************************************************************
SUBROUTINEABC（NE,NG,NWA,IJM,XY,BCA）
C形参说明
C输入:
CNE整型,结构单元总数
CNG整型,结构结点总数
CNWA整型,单元参数输出控制,0－不输出,1－输出
CIJM（3,NE）整型,单元结点编码数组
CXY（2,NG）实型,结构结点坐标数组
C输出:
CBCA（7,NE）实型,结构单元参数数组
C变量说明
CX（2,5）实型,当前计算单元的结点坐标数组
CB（7）实型,当前计算单元的单元参数数组
DIMENSIONIJM（3,NE）,XY（2,NG）,BCA（7,NE）,X（2,5）,B（7）
C
IF（NWA.EQ.1）WRITE（6,5）
5FORMAT（/10X,'PARAMETERSOFELEMENTSBCA（7,NE）'/）
C若控制打开,则输出单元参数提示信息
DO80I=1,NE
C遍历所有单元
DO10K=1,3
C遍历单元内的3个结点
K1=IJM（K,I）
C取结点在结构中对应的结点号
DO10J=1,2
X（J,K）=XY（J,K1）
C得到当前结点的坐标值
10CONTINUE
DO20J=1,2
X（J,4）=X（J,1）
X（J,5）=X（J,2）
20CONTINUE
C为编程方便,每单元多存两个结点坐标
DO30K=1,3
B（K）=X（2,K+1）-X（2,K+2）
C计算Bm
B（K+3）=X（1,K+2）-X（1,K+1）
C计算Cm
30CONTINUE
B（7）=（B
（1）*B（5）-B（4）*B
（2））*0.5
C计算单元面积A
IF（NWA.GT.0）WRITE（6,40）I,B
40FORMAT（1X,'NE=',I3,/3X,7E10.4）
C若输出控制打开,则输出单元号和对应的参数
IF（B（7）.LE.0.0）GOTO60
C若当前单元面积为负,则出错
DO50J=1,7
BCA（J,I）=B（J）
50CONTINUE
C将当前单元参数数组中的值传送给输出数组
GOTO80
60WRITE（6,70）I,（IJM（J,I）,J=1,3）
70FORMAT（/5X,'ELEMENT',I5,5X,'AREAISNONPOSITIVE',5X,'IJM',3I5）
STOP333
C显示出错信息并停止运行
80CONTINUE
RETURN
END
C*********************************************************************
C3*
C子过程名称:
KE*
C子过程功能:
根据结构单元参数数组,计算出*
C指定单元的单元刚度矩阵*
C*
C*********************************************************************
SUBROUTINEKE（IO,NE,NWE,T,A1,A2,V,EK,BCA）
C形参说明
C输入:
CIO整型,计算单元编号
CNE整型,结构单元总数
CNWE整型,刚度矩阵输出控制,0－不输出,1－输出
CIJM（3,NE）整型,单元结点编码数组
CXY（2,NG）实型,结构结点坐标数组
CT实型,单元厚度
CA1实型,材料系数,A1=E/（4*（1-V**2））
CA2实型,材料系数,A2=（1-V）/2
CV实型,泊松比
CBCA（7,NE）实型,结构单元参数数组
C输出:
CEK（6,6）实型,单元刚度矩阵
DIMENSIONB（7）,BCA（7,NE）,EK（6,6）
DO10I=1,7
B（I）=BCA（I,IO）
10CONTINUE
C得到计算单元的参数
A=A1/B（7）*T
DO20I=1,3
DO20J=I,3
C将单元刚度矩阵分块成3×3个子矩阵
I1=2*I
J1=2*J
EK（I1-1,J1-1）=A*（B（I）*B（J）+A2*B（I+3）*B（J+3））
EK（I1-1,J1）=A*（V*B（I）*B（J+3）+A2*B（I+3）*B（J））
EK（I1,J1-1）=A*（V*B（I+3）*B（J）+A2*B（I）*B（J+3））
EK（I1,J1）=A*（B（I+3）*B（J+3）+A2*B（I）*B（J））
C计算每个子矩阵各元素的值
20CONTINUE
DO30I=3,6
DO30J=1,I
EK（I,J）=EK（J,I）
30CONTINUE
C根据对称性得到左下角矩阵的值
IF（NWE.EQ.0）GOTO60
C若输出控制关闭,则直接结束子过程
WRITE（6,40）IO
40FORMAT（/1X,'EKNE=',I5）
WRITE（6,50）EK
50FORMAT（1X,6E11.4）
C输出单元号和对应的刚度矩阵
60RETURN
END
C*********************************************************************
C4*
C子过程名称:
SUMK*
C子过程功能:
将指定单元的单元刚度矩阵集成到结构刚度*
C矩阵中,结构刚度矩阵以二维等带宽方式存储*
C*
C*********************************************************************
SUBROUTINESUMK（IO,NE,ND,NT,IJM,SK,EK）
C形参说明
C输入:
CIO整型,计算单元编号
CNE整型,结构单元总数
CND整型,结构刚度矩阵的半带宽
CNT整型,结构刚度矩阵的阶数
CIJM（3,NE）实型,单元结点编码数组
CEK（6,6）实型,单元刚度矩阵
C输出:
CSK（NT,ND）实型,结构刚度矩阵
DIMENSIONIJ（3）,SK（NT,ND）,IJM（3,NE）,EK（6,6）
C
DO
升级会员 