《有限单元法基本原理和数值方法》一书的源程序

《有限单元法基本原理和数值方法》一书的源程序--------------------------------------------------------------------------------!*********************************************************************!*!PL----PROGRAMOFPLANEPROBLEM96.1*!*!*********************************************************************CC--------输入数据顺序--------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若NF0则填写C8.MF(2,NF),ZF(NF)2×NF整型,NF实型CMF(1,I)---第I个集中载荷所在的结点号CMF(2,I)=1--给定X方向集中力C=0--给定Y方向集中力CZF(I)-----作用的集中力值C若NP0则填写C9.MP(2,NP),ZP(NP)2×NP整型,NP实型CMP(1,I)----第I个载荷作用边的起始结点号CMP(2,I)----第I个载荷作用边的起始结点号CZP(I)------第I个均布载荷值CC若NX1重复7.-9.(NX-1)次CC最后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实型CC若NWA=1,则输出C3.I,B(7)单元参数NE行,1×NE整型,7×NE实型C每行结构为:'NE='+单元号+Bi+Bj+Bm+Ci+Cj+Cm+ACC若NWE=1,则输出C4.IO,EK(6×6)单元刚度阵NE行,1×NE整型,6×6×NE实型C每行结构为:'NE='+单元号+EK(单元刚度阵)CC若NWK=1,则输出C5.SK(NT,ND)结构刚度矩阵NT×ND=2NG×ND实型CC若NWD=1,则输出C6.I,B结点位移数据NG行,1×NG整型,2×NG实型C每行结构为:单元号+U+VCC7.S1,S2,S3,X1,X2,CTAC单元应力数据6×NE实型C分别代表σx,σy,τxy,σ1,σ2和主应力方向CC若NX1重复6.-7.(NX-1)次C!--------可调数组分配--------CC实型数组C(100000)整型数组IA(100000)CC(1)XY(2,NG)IA(1)IJM(3,NE)CC(N1)ZB(NIA(M1)MB(2,MCC(N2)BCA(7,NE)IA(M2)MF(2,NCC(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结构刚度矩阵主元非正!*********************************************************************CC主程序CDIMENSIONC(500000),IA(50000),EK(36)CHARACTER*12IN,OUTCIN和OUT为输入文件和输出文件的文件名WRITE(*,*)''WRITE(*,*)'PLEASEINPUTTHEINPUT-FILENAME(A12)'WRITE(*,*)''READ(*,5)INC输入输入文件的文件名WRITE(*,*)''WRITE(*,*)'PLEASEINPUTTHEOUTPUT-FILENAME(A12)'WRITE(*,*)''READ(*,5)OUTC输入输出文件的文件名5FORMAT(A12)OPEN(5,FILE=IN,STATUS='OLD')OPEN(6,FILE=OUT,STATUS='UNKNOWN')C打开对应的输入和输出文件10READ(5,*)NGIF(NG.EQ.0)STOPC输入结构结点数;如果结点数为0则停止运行READ(5,*)NE,MC,NX,NB,ND,EO,VO,TC按顺序输入结构单元数,问题类型参数,载荷组数,给定位移个数C结构刚度阵的半带宽,弹性模量,泊松比和结构厚度READ(5,*)NWA,NWE,NWK,NWP,NWDC按顺序输入各输出控制参数NT=2*NGC确定总刚度矩阵阶数NTCC计算变界数组的下限!M1=3*NE+1M2=M1+2*NBN1=2*NG+1N2=N1+NBN3=7*NE+N2N4=N3+NT*NDN5=N4+NTC得到各变界数组在一维大数组中的起始元素编号CC检验实型数组C的下限CNEND=N5IF(NEND.LE.500000)GOTO35WRITE(*,*)'***EXCEEDTHELIMITOFARRAYC(INTHEMIDDLE)!!***'WRITE(*,30)NEND30FORMAT(/,'********NEND=',I6,1X,'80000********')STOP111C若C下限超出500000,则给出错误信息并停止运行CC数据输入C35CALLINPUT(NE,NG,NB,IA(1),C(1),IA(M1),C(N1))C调用INPUT子程输入数据WRITE(*,40)40FORMAT(/10X,'#####INPUTPASSED#####')C显示提示信息IF(MC.EQ.0)GOTO45C检验是否是平面应力问题CC平面应变问题CE=EO/(1.0-VO*VO)V=VO/(1.0-VO)C平面应变问题时,先进行弹性常数替换GOTO50CC平面应力问题C45E=EOV=VO50NX1=NXA1=E/(1.0-V*V)/4.0A2=0.5*(1.0-V)C初始化NX1,A1和A2/NX1为剩余载荷的组数CC计算单元参数CCALLABC(NE,NG,NWA,IA(1),C(1),C(N2))WRITE(*,55)55FORMAT(/10X,'#####ABCPASSED#####')C调用ABC子程计算单元参数并显示提示信息CC集成结构刚度矩阵KCDO60I=N3,N4C(I)=0.060CONTINUEC初始化结构刚度矩阵SKDO65K=1,NEC遍历结构的所有单元IO=KCALLKE(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子程将其集成到结构刚度阵中65CONTINUEWRITE(*,70)70FORMAT(/10X,'#####SUMKPASSED#####')C显示提示信息CALLCHECK(NT,ND,NWK,C(N3))C调用CHECK子程检验结构刚度阵中的主元是否非正WRITE(*,75)75FORMAT(/10X,'#####CHECKPASSED#####')C显示提示信息80READ(5,*)NF,NPC输入集中载荷个数NF和均布载荷个数NPCC再次计算变界数组的下限C并检验实型数组C和整型数组IA的下限CM3=M2+2*NFN6=N5+NFNEND=N6+NP-1MEND=M3+2*NP-1C计算C和IA的下限NM=0IF(NEND.LE.500000)GOTO85WRITE(*,*)'***EXCEEDTHELIMITOFARRAYC(ATTHEEND)!!***'WRITE(*,30)NENDNM=185IF(MEND.LE.50000)GOTO95WRITE(*,*)'***EXCEEDTHELIMITOFARRAYIA(ATTHEEND)!!***'WRITE(*,90)MEND90FORMAT(/,'********MEND=',I6,1X,'500********')STOP22295IF(NM.EQ.1)STOP222C检验两个数组的下限,若下限超出则给出错误信息并停止运行CC集成结构结点载荷列阵PCDO100I=N4,N5C(I)=0.0100CONTINUEC初始化结构结点载荷列阵FIF(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子程集成各均布载荷并显示提示信息CC引入给定位移CCALLDBC(NT,ND,NB,NX,NX1,C(N3),C(N4),IA(M1),C(N1))WRITE(*,115)115FORMAT(/10X,'#####DBCPASSED#####')C调用DBC子程引入给定位移消除系数矩阵的奇异性,并显示提示信息CC求解线性方程组KA=PCCALLGAUSS(NT,ND,NWD,NX,NX1,C(N3),C(N4))WRITE(*,120)120FORMAT(/10X,'#####GAUSSPASSED#####')C调用GAUSS子程求解线性方程组并显示提示信息CC计算单元应力CCALLSTRESS(NE,NT,A1,A2,V,IA(1),C(N2),C(N4))WRITE(*,125)125FORMAT(/10X,'#####STRESSPASSED#####')C调用STRESS子程输出各单元应力并显示提示信息NX1=NX1-1IF(NX1.GT.0)GOTO80C剩余载荷组自减1;若还有载荷剩余则继续计算GOTO10C重新输入ENDC*********************************************************************C1*C子过程名称:INPUT*C子过程功能:按顺序输入并输出计算所需数据*C*C*********************************************************************SUBROUTINEINPUT(NE,NG,NB,IJM,XY,MB,ZC形参说明C输入:CNE整型,结构单元总数CNG整型,结构结点总数CNB整型,给定位移的个数CIJM(3,NE)整型,