八节点四边形等参单元有限元程序

CC-----------PARAMETRIC FEM FOR PLANE PROBLEMS-----------------

C

COMMON/CONTR/NPOIN,NELEM,NNODE,NDOFN,NDIME,NGAUS,NPROP,NMATS COMMON/CONTR/NVFIX,NEVAB,NSTRE,NTYPE,NTOTV

COMMON/LGDAT/COORD(40,2),PROPS(9,4),PRESC(40,2),ASDIS(80)

COMMON/LGDAT/ELOAD(9,16),NOFIX(40),IFPRE(40,2),LNODS(9,8),MATNO(9) COMMON/WORKS/ELCOD(2,8),SHAPE(8),DERIV(2,8),CARTD(2,8)

COMMON/WORKS/POSGP(3),WEIGP(3),GPCOD(2,9),BMATX(3,16),DMATX(3,3) COMMON/WORKS/SMATX(3,16,9),DBMAT(3,16)

COMMON/GENEL/ASLOD(80),ASTIF(80,80)

CHARACTER *12 NAME

C P231-232

WRITE(*,210)

210 FORMAT(/3X,'PLEASE ENTER DATAFILE NAME:')

READ(*,220) NAME

220 FORMAT(A12)

REWIND 1

REWIND 3

OPEN(5,FILE=NAME,status='old')

OPEN(1,FILE='STIF.dat',form='unformatted')

OPEN(3,FILE='SMAT.dat',form='unformatted')

OPEN(6,FILE='OUT1.txt')

C

NNODE=8

NDOFN=2

NEVAB=NNODE*NDOFN

NDIME=2

NSTRE=3

NPROP=4

READ(5,*)NPOIN,NELEM,NVFIX,NMATS,NGAUS,NTYPE

NGASP=NGAUS*NGAUS

NTOTV=NPOIN*NDOFN

WRITE(6,100)NPOIN,NELEM,NVFIX,NMATS,NGAUS,NTYPE

WRITE(6,110)

DO 10 LELEM=1,NELEM

READ(5,*) IELEM,MATNO(IELEM),(LNODS(IELEM,INODE),INODE=1,NNODE) WRITE(6,115) IELEM,MATNO(IELEM),(LNODS(IELEM,INODE),INODE=1,NNODE) 10 CONTINUE

WRITE(6,120)

WRITE(6,125)

DO 20 IPOIN=1,NPOIN

20 READ(5,*) JPOIN,(COORD(JPOIN,IDIME),IDIME=1,NDIME)

c

WRITE(6,135)(I,(COORD(I,IDIME),IDIME=1,NDIME),I=1,NPOIN)

WRITE(6,140)

WRITE(6,145)

DO 80 IVFIX=1,NVFIX

READ(5,*) NOFIX(IVFIX),(IFPRE(IVFIX,IDOFN),IDOFN=1,NDOFN),

1(PRESC(IVFIX,IDOFN),IDOFN=1,NDOFN)

WRITE(6,155) NOFIX(IVFIX),(IFPRE(IVFIX,IDOFN),IDOFN=1,NDOFN), 1(PRESC(IVFIX,IDOFN),IDOFN=1,NDOFN)

80 CONTINUE

WRITE(6,160)

WRITE(6,165)

DO 90 IMATS=1,NMATS READ(5,*) NUMAT,(PROPS(NUMAT,IPROP),IPROP=1,NPROP)

WRITE(6,170) NUMAT,(PROPS(NUMAT,IPROP),IPROP=1,NPROP)

90 CONTINUE

100 FORMAT(1X,'NPOIN=',I3,1X,'NELEM=',I3,1X,'NVFIX=',I3,1X,

1 'NMATS='I3,1X,'NGAUS='I2,1X,'NTYPE=',I2)

110 FORMAT(/1X,'ELEMENT',3X,'PROPERTY',6X,'NODE NUMBER')

115 FORMAT(1X,I5,I9,6X,8I4)

120 FORMAT(/24H NODAL POINT COORDINATES)

125 FORMAT(/2(7H NODE ,7X,1HX,9X,1HY,5X))

135 FORMAT(2(1X,I4,2X,2F10.3,2X))

140 FORMAT(/16HRESTRAINED NODES)

145 FORMAT(/1X,'NODE',4X,'CODE',6X,'FIXED VALUES')

155 FORMAT(1X,I4,5X,2I1,2F10.5)

160 FORMAT(/1X,'MATERAL PROPERTIES')

165 FORMAT(/1X,'NUMBER',7X,'PROPERTIES')

170 FORMAT(1X,I4,4X,4E14.4)

C

CALL GAUSSQ

CALL STIFPS

CALL LOADPS

CALL ASSEMB

C

CALL GAUSS(ASTIF,ASLOD,NTOTV)

DO 175 I=1,NTOTV

175 ASDIS(I)=ASLOD(I)

WRITE(6,180)

180 FORMAT(/1X,'NODE',5X,'X-DISP',5X,'Y-DISP')

WRITE(6,185) (I,ASDIS(2*I-1),ASDIS(2*I),I=1,NPOIN)

185 FORMAT(1X,I4,2E12.4,4X,I2,2E12.4)

C

CALL STREPS

STOP

END

C

SUBROUTINE GAUSSQ

C----------------------------------------------------------------CC1 P65 C

C GIVES THE WEIGHTING FACTORS FOR GAUSS POINTS CC----------------------------------------------------------------C COMMON/CONTR/NPOIN,NELEM,NNODE,NDOFN,NDIME,NGAUS,NPROP,NMATS COMMON/CONTR/NVFIX,NEVAB,NSTRE,NTYPE,NTOTV

COMMON/WORKS/ELCOD(2,8),SHAPE(8),DERIV(2,8),CARTD(2,8)

COMMON/WORKS/POSGP(3),WEIGP(3),GPCOD(2,9),BMATX(3,16),DMATX(3,3) COMMON/WORKS/SMATX(3,16,9),DBMAT(3,16)

IF(NGAUS.GT.2) GO TO 10

POSGP(1)=-0.577350269189626

WEIGP(1)=1.0

GO TO 20

10 POSGP(1)=-0.774596669241483

POSGP(2)=0.0

WEIGP(1)=0.5555555555555556

WEIGP(2)=0.8888888888888889

20 KGAUS=NGAUS/2

DO 30 IGASH=1,KGAUS

JGASH=NGAUS+1-IGASH POSGP(JGASH)=-POSGP(IGASH)

WEIGP(JGASH)=WEIGP(IGASH)

30 CONTINUE

RETURN

END

C

SUBROUTINE STIFPS

C-----------------------------------------------------------------CC2 P114-116 C

C CALCUATES THE STIFFNESS MATRIX OF ELEMENT CC-----------------------------------------------------------------C COMMON/CONTR/NPOIN,NELEM,NNODE,NDOFN,NDIME,NGAUS,NPROP,NMATS COMMON/CONTR/NVFIX,NEVAB,NSTRE,NTYPE,NTOTV

COMMON/LGDAT/COORD(40,2),PROPS(9,4),PRESC(40,2),ASDIS(80)

COMMON/LGDAT/ELOAD(9,16),NOFIX(40),IFPRE(40,2),LNODS(9,8),MATNO(9) COMMON/WORKS/ELCOD(2,8),SHAPE(8),DERIV(2,8),CARTD(2,8)

COMMON/WORKS/POSGP(3),WEIGP(3),GPCOD(2,9),BMATX(3,16),DMATX(3,3) COMMON/WORKS/SMATX(3,16,9),DBMAT(3,16)

COMMON/GENEL/ASLOD(80),ASTIF(80,80)

DIMENSION ESTIF(16,16)

C

DO 70 IELEM=1,NELEM

LPROP=MATNO(IELEM)

C

C EVALUATE THE COORDINATES

C

DO 10 INODE=1,NNODE

LNODE=LNODS(IELEM,INODE)

DO 10 IDIME=1,NDIME

10 ELCOD(IDIME,INODE)=COORD(LNODE,IDIME)

C

CALL MODPS(LPROP)

THICK=PROPS(LPROP,3)

C

DO 20 IEVAB=1,NEVAB

DO 20 JEVAB=1,NEVAB

20 ESTIF(IEVAB,JEVAB)=0.0

KGASP=0

C

C ENTER NUMERICAL INTEGRATION

C

DO 50 IGAUS=1,NGAUS

DO 50 JGAUS=1,NGAUS

KGASP=KGASP+1

EXISP=POSGP(IGAUS)

ETASP=POSGP(JGAUS)

C

CALL SFR2(EXISP,ETASP)

CALL JACOB2(IELEM,DJACB,KGASP)

DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)

IF (THICK.NE.0.0) DVOLU=DVOLU*THICK

CALL BMATPS

CALL DBE

C

C CALCULATE THE ELEMENT STIFFNESSESC

DO 30 IEVAB=1,NEVAB

DO 30 JEVAB=IEVAB,NEVAB

DO 30 ISTRE=1,NSTRE

30 ESTIF(IEVAB,JEVAB)=ESTIF(IEVAB,JEVAB)+

1 BMATX(ISTRE,IEVAB)*DBMAT(ISTRE,JEVAB)*DVOLU

C

C STORE DB MATRIX

C

DO 40 ISTRE=1,NSTRE

DO …… 此处隐藏：15259字，全部文档内容请下载后查看。喜欢就下载吧 ……