平面弹性力学有限元源程序(FORTRAN)

1平面弹性力学有限元源程序（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，则表示计所有单元的自重，不需输入计2自重的单元号；若NMNE，则需要输入计自重的单元号；$DEBUGPROGRAMPLANEIMPLICITREAL*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*12WRITE(*,'(A)')'本程序为计算平面问题的有限元程序'WRITE(*,'(A)')'特点：(1)采用三结点三角形单元；'WRITE(*,'(A)')'(2)采用等带宽存贮技术；'WRITE(*,'(A)')'(3)采用高斯消元法解线性方程组。'WRITE(*,'(/A)')'输入计算问题名(PN):'READ(*,'(A)')PNCALLFNAME(PN,'.DAT',FN)WRITE(*,'(2A)')'输入数据文件名为：',FNOPEN(5,FILE=FN,STATUS='OLD')CALLFNAME(PN,'.OUT',FN)WRITE(*,'(/2A)')'结果输出数据文件名为：',FNOPEN(6,FILE=FN,STATUS='UNKNOWN')CALLFNAME(PN,'.OU1',FN)WRITE(*,'(/2A)')'参数输出数据文件名为：',FNOPEN(7,FILE=FN,STATUS='UNKNOWN')READ(5,*)NG,NE,MC,NX,NB,EO,VO,DENSITY,TWRITE(6,120)NG,NE,MC,NX,NBWRITE(6,130)EO,VO,DENSITY,TREAD(5,*)NWA,NWE,NWK,NWPNT=2*NGALLOCATE(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(THENSTOP000ENDIFCALLINPUT(NG,NE,NB,IJK,XY,MB,ZB)CALLCALND(NG,NE,IJK,ND)ALLOCATE(SK(NT,ND))IF(THENE=EOV=VOELSE3E=EO/(1-VO*VO)V=VO/(1-VO)ENDIFIP=0NX1=NXA1=E/(1-V*V)/4.0A2=0.5*(1-V)CALLABC(NG,NE,IJK,XY,BCA,NWA)CALLCLEAR(NT,ND,SK)DO100K=1,NECALLKE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE)CALLSUMK(K,NE,NT,ND,IJK,EK,SK)100CONTINUECALLCHECK(NT,ND,SK)110CONTINUEIP=IP+1WRITE(6,'(///1X,A,I4)')荷载工况=,IPREAD(5,*)NF,NP,NMDOI=1,NTB(I)=0.0ENDDOIF(THENCALLPF(NF,NT,B)ENDIFIF(THENCALLPP(NP,NG,NT,T,XY,B)ENDIFIF(THENCALLPG(NM,NE,NT,DENSITY,T,IJK,BCA,B)ENDIFDOI=1,NBI1=2*(MB(1,I)-1)+2-MB(2,I)DELSUP(I,IP)=-B(I1)ENDDOIF(THENCALLDBC(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)ELSE4WRITE(*,'(2x,A,I4,A)')荷载工况=,IP,没有荷载!WRITE(6,'(2x,A,I4,A)')荷载工况=,IP,没有荷载!ENDIFNX1=NX1-1IF(GOTO110120FORMAT(1X,'结点总数=',I3,1X,'单元总数=',I3,1X,'问题类型=',&&I1,1X,'荷载工况数=',I2,1X,'支承位移数=',I2)130FORMAT(/1X,'弹性模量=',E10.3,1X,'泊松比=',F5.3,1X,'密度=',E10.3,&&1X,'厚度=',F6.3)ENDSUBROUTINEGAUSS(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,NTDOJ=1,NDIF((I+J-1).LE.NT)THENA(I,I+J-1)=SK(I,J)ENDIFENDDOENDDOIF(THENWRITE(7,'(A)')结构总刚(列出右上三角元素)DOI=1,NTWRITE(7,100)(A(I,J),J=1,NT)ENDDOENDIFIF(THENWRITE(7,'(1X,A,I3,A)')第,NX-NX1+1,荷载工况的荷载列阵DOI=1,NTWRITE(7,'(1x,E11.4)')B(I)ENDDOENDIFDO50K=1,NT-1DO60I=K+1,NTC1=A(K,I)/A(K,K)DO70J=I,NTA(I,J)=A(I,J)-C1*A(K,J)70CONTINUEB(I)=B(I)-C1*B(K)60CONTINUE50CONTINUEB(NT)=B(NT)/A(NT,NT)DO80I=NT-1,1,-15DO90J=I+1,NTB(I)=B(I)-A(I,J)*B(J)90CONTINUEB(I)=B(I)/A(I,I)80CONTINUE100FORMAT(1x,4000E10.3)ENDSUBROUTINECHECK(NT,ND,SK)IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N)DIMENSIONSK(NT,ND)M=0DOI=1,NTIF(SK(I,1).THENM=M+1ENDIFENDDOIF(THENWRITE(*,*)总刚矩阵的对角元素为负值!!STOP2222ENDIFENDSUBROUTINESTRESS(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,NEDOI=1,3II=IJK(I,K)D1(2*(I-1)+1)=B(2*(II-1)+1)D1(2*(I-1)+2)=B(2*(II-1)+2)ENDDOCALLCLEAR(3,6,S)C1=2*A1/BCA(7,K)DO20I=1,3S(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)20CONTINUECALLMATMUL(3,6,1,S,D1,D2)DOI=1,3STR(I,K)=D2(I)6ENDDO10CONTINUEENDSUBROUTINEOUTPUT(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*6WRITE(6,20)IPWRITE(6,30)DOI=1,NGWRITE(6,40)I,DELD(1,I,IP),TOD(1,I),DELD(2,I,IP),TOD(2,I)ENDDOWRITE(6,50)DOJ=1,NEWRITE(6,60)J,(DELST(L,J,IP),TOST(L,J),L=1,3)ENDDOWRITE(6,70)DOJ=1,NBIF(MB(2,J).EQ.1)THENVE='x方向'ELSEVE='Y方向'ENDIFWRITE(6,80)MB(1,J),VE,DELSUP(J,IP),TOTSUP(J)ENDDO20FORMAT(/1X,'荷载工况=',I3,'的计算结果')30FORMAT(/,1X,'结点号',5X,'