plastic_beam_fencheng_module_内部文件.for

欧文的程序

		ESTIF(2,3)=-VALU1
		ESTIF(2,4)=-VALU3+VALU4
		ESTIF(3,3)=VALU2
		ESTIF(3,4)=-VALU1
		ESTIF(4,4)=VALU3+VALU4
		DO ISTIF=1,4
			DO JSTIF=ISTIF,4
				ESTIF(JSTIF,ISTIF)=ESTIF(ISTIF,JSTIF)
			ENDDO
		ENDDO
		WRITE(1)ESTIF
!		DO ISTIF=1,4
!			WRITE(1,100)(ESTIF(ISTIF,JSTIF),JSTIF=1,4)
!		ENDDO
!100		FORMAT(4F25.5)
	ENDDO
	RETURN
	END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     求每层的EI和GA,并求和
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE LAYER(IELEM,EIVAL,SVALU)
	USE CONFIG_ARRAY
	IMPLICIT NONE
	DOUBLEPRECISION::EIVAL,SVALU,SHEAR,HARDS,THKTO,ZMIDL
	DOUBLEPRECISION::YOUNG,BRDTH,THICK
	INTEGER::LPROP,KLAYR,KOUNT,ILAYR,IELEM
	EIVAL=0.0
	SVALU=0.0
	LPROP=MATNO(IELEM)
	KLAYR=(IELEM-1)*NLAYR
	SHEAR=PROPS(LPROP,2)
	HARDS=PROPS(LPROP,4)
	THKTO=PROPS(LPROP,5)
	ZMIDL=-THKTO/2.0
	KOUNT=5
	DO ILAYR=1,NLAYR
		KLAYR=KLAYR+1
		YOUNG=PROPS(LPROP,1)
		IF(PLAST(KLAYR).NE.0) YOUNG=YOUNG*(1.0-YOUNG/(YOUNG+HARDS))
!		IF(YOUNG.EQ.0)YOUNG=1.0E-10
		KOUNT=KOUNT+1
		BRDTH=PROPS(LPROP,KOUNT)
		KOUNT=KOUNT+1
		THICK=PROPS(LPROP,KOUNT)
		ZMIDL=ZMIDL+THICK/2.0
		EIVAL=EIVAL+YOUNG*BRDTH*THICK*ZMIDL*ZMIDL
		SVALU=SVALU+SHEAR*BRDTH*THICK
		ZMIDL=ZMIDL+THICK/2.0
	ENDDO
	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE GREDUC
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::KOUNT,NEQNS,IEQNS,IEQN1,ICOLS,IROWS
	DOUBLEPRECISION::PIVOT,FACTR 
	KOUNT=0
	NEQNS=NSVAB
	DO  IEQNS=1,NEQNS
		IF(IFPRE(IEQNS).EQ.1) GOTO 40
		PIVOT=ASTIF(IEQNS,IEQNS)
		IF(ABS(PIVOT).LT.1.0E-10)GOTO 60
		IF(IEQNS.EQ.NEQNS)GOTO 70
		IEQN1=IEQNS+1
		DO IROWS=IEQN1,NEQNS
			KOUNT=KOUNT+1
			FACTR=ASTIF(IROWS,IEQNS)/PIVOT
			FRESV(KOUNT)=FACTR
			IF(FACTR.EQ.0.0) GOTO 30
			DO  ICOLS=IEQNS,NEQNS
				ASTIF(IROWS,ICOLS)=ASTIF(IROWS,ICOLS)-
     $						FACTR*ASTIF(IEQNS,ICOLS)
			ENDDO
			ASLOD(IROWS)=ASLOD(IROWS)-FACTR*ASLOD(IEQNS)
30		ENDDO
		GOTO 70
40		DO IROWS=IEQNS,NEQNS
			ASLOD(IROWS)=ASLOD(IROWS)-ASTIF(IROWS,IEQNS)*FIXED(IEQNS)
		ENDDO
		GOTO 70
60		WRITE(6,900)
900		FORMAT(5X,'INCORRECT PIVOT')
		STOP
70	ENDDO
	RETURN
	END		
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE BAKSUB
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::NEQNS,IEQNS,NEQN1,NBACK,NBAC1,IPOIN,IDOFN,KOUNT,IELEM
	INTEGER::ICOLS,INODE,NLOCA
	INTEGER::NPOSN,IEVAB
	DOUBLEPRECISION::PIVOT,RESID 
	NEQNS=NSVAB
	DO IEQNS=1,NEQNS
		REACT(IEQNS)=0.0
	ENDDO
	NEQN1=NEQNS+1
	DO IEQNS=1,NEQNS
		NBACK=NEQN1-IEQNS
		PIVOT=ASTIF(NBACK,NBACK)
		RESID=ASLOD(NBACK)
		IF(NBACK.EQ.NEQNS) GOTO 30
		NBAC1=NBACK+1
		DO ICOLS=NBAC1,NEQNS
			RESID=RESID-ASTIF(NBACK,ICOLS)*XDISP(ICOLS)
		ENDDO
30		IF(IFPRE(NBACK).EQ.0) XDISP(NBACK)=RESID/PIVOT
		IF(IFPRE(NBACK).EQ.1) XDISP(NBACK)=FIXED(NBACK)
		IF(IFPRE(NBACK).EQ.1) REACT(NBACK)=-RESID
	ENDDO
	KOUNT=0
	DO IPOIN=1,NPOIN
		DO IDOFN=1,NDOFN
			KOUNT=KOUNT+1
			TDISP(IPOIN,IDOFN)=TDISP(IPOIN,IDOFN)+XDISP(KOUNT)
			TREAC(IPOIN,IDOFN)=TREAC(IPOIN,IDOFN)+REACT(KOUNT)
		ENDDO
	ENDDO
	DO IPOIN=1,NPOIN
		DO IELEM=1,NELEM
			DO INODE=1,NNODE
				NLOCA=LNODS(IELEM,INODE)
				IF(IPOIN.EQ.NLOCA) GOTO 70
			ENDDO
		ENDDO
70		DO IDOFN=1,NDOFN
			NPOSN=(IPOIN-1)*NDOFN+IDOFN
			IEVAB=(INODE-1)*NDOFN+IDOFN
			TLOAD(IELEM,IEVAB)=TLOAD(IELEM,IEVAB)+REACT(NPOSN)
		ENDDO
	ENDDO
	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE RESOLV
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::KOUNT,NEQNS,IEQNS,IEQN1,IROWS
	DOUBLEPRECISION::FACTR

	KOUNT=0
	NEQNS=NSVAB
	DO IEQNS=1,NEQNS
		IF(IFPRE(IEQNS).EQ.1)GOTO 20
		IF(IEQNS.EQ.NEQNS)GOTO 40
		IEQN1=IEQNS+1
		DO IROWS=IEQN1,NEQNS
			KOUNT=KOUNT+1
			FACTR=FRESV(KOUNT)
			IF(FACTR.EQ.0)GOTO 10
			ASLOD(IROWS)=ASLOD(IROWS)-FACTR*ASLOD(IEQNS)
10		ENDDO
		GOTO 40
20		DO IROWS=IEQNS,NEQNS
			ASLOD(IROWS)=ASLOD(IROWS)-ASTIF(IROWS,IEQNS)*FIXED(IEQNS)
30		ENDDO
40	ENDDO
	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE REFORBL
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::IELEM,IEVAB,IDOFN,KLAYR,LPROP,NODE1,NODE2,ILAYR,KOUNT
	DOUBLEPRECISION::YOUNG,SHEAR,YIELD,HARDS,THKTO,ELENG,WNOD1,WNOD2
	DOUBLEPRECISION::THTA1,THTA2,THICK
	DOUBLEPRECISION::ZMIDL,STLIN,STCUR,PREYS,ESCUR,RFACT,REDUC,BRDTH 
	DOUBLEPRECISION::STRAN(2)
	DO IELEM=1,NELEM
		DO IEVAB=1,NEVAB
			ELOAD(IELEM,IEVAB)=0.0
		ENDDO
		DO IDOFN=1,NDOFN
			STRES(IELEM,IDOFN)=0.0
		ENDDO
	ENDDO
	KLAYR=0
	DO IELEM=1,NELEM

		LPROP=MATNO(IELEM)
		YOUNG=PROPS(LPROP,1)
		SHEAR=PROPS(LPROP,2)
		YIELD=PROPS(LPROP,3)
		HARDS=PROPS(LPROP,4)
		THKTO=PROPS(LPROP,5)

		NODE1=LNODS(IELEM,1)
		NODE2=LNODS(IELEM,2)
	
		ELENG=ABS(COORD(NODE2)-COORD(NODE1))
		WNOD1=XDISP(NODE1*NDOFN-1)
		WNOD2=XDISP(NODE2*NDOFN-1)
		THTA1=XDISP(NODE1*NDOFN)
		THTA2=XDISP(NODE2*NDOFN)

		STRAN(1)=(THTA1-THTA2)/ELENG
		STRAN(2)=(WNOD2-WNOD1)/ELENG-0.5*(THTA1+THTA2)

		ZMIDL=-THKTO/2.0
		KOUNT=5
		DO ILAYR=1,NLAYR
			KLAYR=KLAYR+1
			KOUNT=KOUNT+1
			BRDTH=PROPS(LPROP,KOUNT)
			KOUNT=KOUNT+1
			THICK=PROPS(LPROP,KOUNT)
			ZMIDL=ZMIDL+THICK/2.0
			STLIN=YOUNG*STRAN(1)*ZMIDL
			STCUR=STRSL(KLAYR,1)+STLIN
			PREYS=YIELD+HARDS*ABS(PLAST(KLAYR))

			IF(ABS(STRSL(KLAYR,1)).GE.PREYS)GOTO 20
			ESCUR=ABS(STCUR)-PREYS
			IF(ESCUR.LE.0.0)GOTO 40
			RFACT=ESCUR/ABS(STLIN)
			GOTO 30
20			IF(STRSL(KLAYR,1).GT.0.0.AND.STLIN.LE.0.0)GOTO 40
			IF(STRSL(KLAYR,1).LT.0.0.AND.STLIN.GE.0.0)GOTO 40
			RFACT=1.0
30			REDUC=1.0-RFACT
			STRSL(KLAYR,1)=STRSL(KLAYR,1)+REDUC*STLIN+
     $			 RFACT*YOUNG*(1.0-YOUNG/(YOUNG+HARDS))*STRAN(1)*ZMIDL
			PLAST(KLAYR)=PLAST(KLAYR)+RFACT*STRAN(1)*YOUNG/(YOUNG+
     $			HARDS)*ZMIDL
			GOTO 45
40			STRSL(KLAYR,1)=STRSL(KLAYR,1)+STLIN	
C			ZF=1
45			STRSL(KLAYR,2)=STRSL(KLAYR,2)+STRAN(2)*SHEAR
			STRES(IELEM,1)=STRES(IELEM,1)+STRSL(KLAYR,1)*
     $						BRDTH*THICK*ZMIDL
			STRES(IELEM,2)=STRES(IELEM,2)+STRSL(KLAYR,2)*
     $						BRDTH*THICK 
			ZMIDL=ZMIDL+THICK/2.0
		ENDDO  
		ELOAD(IELEM,1)=ELOAD(IELEM,1)-STRES(IELEM,2)
		ELOAD(IELEM,2)=ELOAD(IELEM,2)+STRES(IELEM,1)
     $	-0.5*ELENG*STRES(IELEM,2)
		ELOAD(IELEM,3)=ELOAD(IELEM,3)+STRES(IELEM,2)
		ELOAD(IELEM,4)=ELOAD(IELEM,4)-STRES(IELEM,1)
     $	-0.5*ELENG*STRES(IELEM,2)
	ENDDO
	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	

	SUBROUTINE CONUND
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::ISVAB,IELEM,IEVAB,INODE,NODNO,IDOFN,NPOSN
	DOUBLEPRECISION::RESID,RETOT,RATIO,REFOR 
	DOUBLEPRECISION::STFOR(52),TOFOR(52) 

	NCHEK=0
	RESID=0.0
	RETOT=0.0
	DO ISVAB=1,NSVAB
		STFOR(ISVAB)=0.0
		TOFOR(ISVAB)=0.0
	ENDDO
	DO IELEM=1,NELEM
		IEVAB=0
		DO INODE=1,NNODE
			NODNO=IABS(LNODS(IELEM,INODE))
			DO IDOFN=1,NDOFN
				IEVAB=IEVAB+1
				NPOSN=(NODNO-1)*NDOFN+IDOFN
				STFOR(NPOSN)=STFOR(NPOSN)+ELOAD(IELEM,IEVAB)
				TOFOR(NPOSN)=TOFOR(NPOSN)+TLOAD(IELEM,IEVAB)
			ENDDO
		ENDDO
	ENDDO
 	DO ISVAB=1,NSVAB
 		REFOR=TOFOR(ISVAB)-STFOR(ISVAB)
 		RESID=RESID+REFOR*REFOR
 		RETOT=RETOT+TOFOR(ISVAB)*TOFOR(ISVAB)
 	ENDDO

	DO IELEM=1,NELEM
		DO IEVAB=1,NEVAB
			ELOAD(IELEM,IEVAB)=TLOAD(IELEM,IEVAB)-
     $			ELOAD(IELEM,IEVAB)
		ENDDO
	ENDDO
	RATIO=100.0*SQRT(RESID/RETOT)

	IF(RATIO.GT.TOLER)NCHEK=1
	IF(IITER.EQ.1)GOTO 50
	IF(RATIO.GT.PVALU)NCHEK=999
50	PVALU=RATIO
	WRITE(6,900)IITER,NCHEK,RATIO
900	FORMAT(6X,'ITERATION NUMBER=',I5/ 
     $	   6X,'CONVERGENCE CODE=',I4,3X,
     $	   'NORM OF RESIDUAL SUM RATIO=',E14.6)
	RETURN
	END

!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE RESULT
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::IPOIN,IELEM,IDOFN
	IF (NDOFN.EQ.1) WRITE(6,900)
900	FORMAT( 5X,'NODE',4X,'DISPL.',12X,'REACTINOS')
	IF(NDOFN.EQ.2)WRITE(6,910)
910	FORMAT( 5X,'NODE',4X,'DISPL.',12X,'REACTION',7X
     $		,'DISPL.',12X,'REACTION')
	DO IPOIN =1,NPOIN
		WRITE(6,920)IPOIN,(TDISP(IPOIN,IDOFN),TREAC(IPOIN,IDOFN),
     $		IDOFN=1,NDOFN)
		WRITE(*,920)IPOIN,(TDISP(IPOIN,IDOFN),TREAC(IPOIN,IDOFN),
     $		IDOFN=1,NDOFN)
	ENDDO
920	FORMAT(I10,2(E14.6,5X,E14.6))

C
C  ***荷载挠度曲线数值
C
	DO IPOIN =6,6
		WRITE(7,90)(TDISP(IPOIN,IDOFN),IDOFN=1,1)
		WRITE(*,90)(TDISP(IPOIN,IDOFN),IDOFN=1,1)
	ENDDO
90	FORMAT(F20.6)


	IF(NDOFN.EQ.2)WRITE(6,930) 
930	FORMAT(3X,'ELEMENT',12X,'STRESSES',12X,'PL.STRAIN')
	IF(NDOFN.EQ.1)WRITE(6,940)
940	FORMAT(3X,'ELEMENT',5X,'STRESSES',5X,'PL.STRAIN')
	DO IELEM=1,NELEM
		WRITE(6,950)IELEM,(STRES(IELEM,IDOFN),IDOFN=1,NDOFN),
     $				PLAST(IELEM)
		WRITE(*,950)IELEM,(STRES(IELEM,IDOFN),IDOFN=1,NDOFN),
     $				PLAST(IELEM)
	ENDDO
950	FORMAT(I10,3E14.6)
	RETURN
	END

!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE FNAME(PN,FN2,FN)
          
		IMPLICIT REAL*8(A-H,O-Z)

		CHARACTER PN*40,FN2*4,FN*12
!	   去掉PN中前面的空格
		DO 10 I=1,40
			IF(PN(I:I).EQ.' ') GOTO 10
			IP=I
			GOTO 20
10		CONTINUE

20		CONTINUE

		FN(1:8)=PN(IP:IP+7)
!	   去掉FN中后面的空格
		DO 30 I=8,1,-1
			IF(FN(I:I).EQ.' ') GOTO 30
			IL=I
			GOTO 40
30		CONTINUE

40		CONTINUE
!         生成文件名FN=PN+FN2
		FN(IL+1:IL+4)=FN2(1:4)
	END