1、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） 若N
2、F0，则输入下面数据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） 若NP0，则输入下面数据K,MP（1,K）,MP（2,K）,ZP（K）（每行3个整形数，1个实型数，共计NP行）分布荷载序号； MP（1,K）：第K个分布荷载作用的结点号； MP（2,K）：第K个分布荷载的坐标方向； ZP（K）：第K个分布荷载的数值；（4） 若NM0，则输入下面数据若NMNE，则表示计所有单元的自重，不需输入计自重的单
3、元号；若NMNE，则需要输入计自重的单元号；$DEBUGPROGRAM PLANEIMPLICIT REAL*8（A-H,O-Z）,INTEGER（I-N）ALLOCATABLE:IJK（:,:）,XY（:）,BCA（:）,SK（:）,STR（:）,MB（:）,ZB（:）,B（:）DELD（:）,TOD（:）,DELST（:）,TOST（:）,DELSUP（:）,TOTSUP（: DIMENSION EK（6,6）CHARACTER PN*40,FN*12WRITE（*,（A） 本程序为计算平面问题的有限元程序 特点：（1）采用三结点三角形单元； （2）采用等带宽存贮技术； （3）采用高斯消元
4、法解线性方程组。（/A） 输入计算问题名（PN）:READ（*,） PNCALL FNAME（PN,.DAT,FN）（2A）输入数据文件名为：,FNOPEN（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
5、,NWE,NWK,NWP NT=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） CALL CLEAR（2,NG,TOD） CALL CLEAR（3,NE,TOST） CALL CLEAR1（NB,TOTSUP） IF （NG.EQ.0） THENSTOP 000 ENDIF CALL INPUT（NG,NE,NB,IJK,XY,M
6、B,ZB） CALL CALND（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） CALL ABC（NG,NE,IJK,XY,BCA,NWA） CALL CLEAR（NT,ND,SK） DO 100 K=1,NE CALL KE（K,NE,T,V,A1,A2,IJK,BCA,EK,NWE） CALL SUMK（K,NE,NT,ND,IJK,EK,SK）100 CONTINUE
7、 CALL CHECK（NT,ND,SK）110 IP=IP+1 WRITE（6,（/1X,A,I4） 荷载工况=,IP READ（5,*） NF,NP,NM DO I=1,NT B（I）=0.0 ENDDO IF（NF.GT.0） THEN CALL PF（NF,NT,B） IF（NP.GT.0） THEN CALL PP（NP,NG,NT,T,XY,B） IF（NM.GT.0） THEN CALL PG（NM,NE,NT,DENSITY,T,IJK,BCA,B） DO I=1,NB I1=2*（MB（1,I）-1）+2-MB（2,I） DELSUP（I,IP）=-B（I1） IF（NF.N
8、E.0.OR.NP.NE.0.OR.NM.NE.0） THEN CALL DBC（NB,ND,NT,NX,NX1,MB,ZB,SK,B） CALL GAUSS（NT,ND,NX1,NX,SK,B,NWK,NWP） CALL STRESS（NE,NT,V,A1,A2,IJK,BCA,B,STR） CALL TRES（IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST） CALL SUPT（IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DEL SUP, TOTSUP） CALL OUTPUT（IP,NG,NE,NB,NX,MB,DE
9、LD,TOD ,DELST,TOST,DELSUP, WRITE（*,（2x,A,I4,A）,IP,没有荷载! NX1=NX1-1 IF （NX1.GT.0） GOTO 110120 FORMAT（1X, 结点总数=, I3, 1X, 单元总数=问题类型=, & I1, 1X, 荷载工况数=, I2, 1X, 支承位移数=, I2）130 FORMAT（/1X, 弹性模量=, E10.3, 1X, 泊松比=, F5.3, 1X, 密度=, E10.3, & 1X, 厚度=, F6.3） END SUBROUTINE GAUSS（NT,ND,NX1,NX,SK,B,NWK,NWP） IMPLIC
10、IT REAL*8（A-H,O-Z）,INTEGER（I-N） DIMENSION SK（NT,ND）,A（NT,NT）,B（NT） CALL CLEAR（NT,NT,A）DO J=1,NDIF（I+J-1）.LE.NT） THENA（I,I+J-1）=SK（I,J）ENDIFENDDO 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） DO 50
11、K=1,NT-1DO 60 I=K+1,NTC1=A（K,I）/A（K,K）DO 70 J=I,NTA（I,J）=A（I,J）-C1*A（K,J）70CONTINUEB（I）=B（I）-C1*B（K）6050 B（NT）=B（NT）/A（NT,NT） DO 80 I=NT-1,1,-1DO 90 J=I+1,NTB（I）=B（I）-A（I,J）*B（J）90B（I）=B（I）/A（I,I）80 FORMAT（1x, 4000E10.3） SUBROUTINE CHECK（NT,ND,SK） DIMENSION SK（NT,ND） M=0IF（SK（I,1）.LE.0.0） THENM=M+1
12、IF （M.GT.0） THEN WRITE（*,*） 总刚矩阵的对角元素为负值! STOP 2222 SUBROUTINE STRESS（NE,NT,V,A1,A2,IJK,BCA,B,STR） DIMENSION D1（6）,D2（3）,B（NT）,S（3,6）,IJK（3,NE）,BCA（7,NE）,STR（3,NE） CALL CLEAR（3,NE,STR） DO 10 K=1,NEDO I=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）CALL CLEAR（3,6,S）C1=2*A1/BCA（7
13、,K）DO 20 I=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）20CALL MATMUL（3,6,1,S,D1,D2）STR（I,K）=D2（I）10 SUBROUTINE OUTPUT（IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP
14、） DIMENSION MB（2,NB）,DELD（2,NG,NX）,TOD（2,NG） DIMENSION DELST（3,NE,NX）,TOST（3,NE）,DELSUP（NB,NX）,TOTSUP（NB） CHARACTER VE*6 WRITE（6,20） IP WRITE（6,30） DO I=1,NG WRITE（6,40） I,DELD（1,I,IP）,TOD（1,I）,DELD（2,I,IP）,TOD（2,I） WRITE（6,50） DO J=1,NE WRITE（6,60） J,（DELST（L,J,IP）,TOST（L,J）,L=1,3） WRITE（6,70） DO J
15、=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
16、, E10.3, 5X, E10.3, 6X, E10.3, & 6X, E10.3）支座结点号支反力方向支反力增量 6x, 支反力累计量 FORMAT（5X, I3, 12X, A, 6x, E10.3, 8X, E10.3）END SUBROUTINE MATMUL（M,N,L,A,B,C）DIMENSION A（M,N）,B（N,L）,C（M,L）DO 100 I=1,MDO 100 J=1,L C（I,J）=0.0 DO 100 K=1,NC（I,J）=C（I,J）+A（I,K）*B（K,J）END SUBROUTINE CLEAR（M,N,A） IMPLICIT REAL*8（A-H,O-Z）,INTEGER
