a05.for

ADINA84有限元编程学习的好例子

      SUBROUTINE ELT2D3
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1            ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIMEL/ N101,N102,N103,N104,N105,N106,N107,N108,N109,N110,
     1               N111,N112,N113,N114,N120,N121,N122,N123,N124,N125
      COMMON /DPR/ ITWO
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON A(1)
C
      REAL A
C
      DIMENSION IA(1)
C
      EQUIVALENCE (NPAR(7),MXNODS),(A(1),IA(1))
      EQUIVALENCE (NPAR(17),NCON)
C
C
C
C
C
      NN=N110+(NEL-1)*MXNODS
      IF(IND.NE.0) GO TO 100
C
C
      CALL ITHEL2(A(NN))
      GO TO 200
C
C
  100 MATP=IA(N107+NEL-1)
C
C
      NM=N109+(MATP-1)*NCON*ITWO
C
C
      ND5DIM=MXNODS-4
      LL=N111+(NEL-1)*ND5DIM
C
C
      CALL THEL2(A(NM),A(NN),A(N6B + ITWO),A(LL))
C
  200 CONTINUE
C
      RETURN
      END
      SUBROUTINE ITHEL2(IWA)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EM2D/ S(378),XM(27),B(4,18),RE(27),EDIS(27),EDISI(27),
     1              XX(27),NOD(9),NODM(9),NOD5M(5)
      COMMON /DPR/ ITWO
C
      DIMENSION IWA(*)
C
C
C
      II=0
      DO 15 K=1,9
      IF(NODM(K).EQ.0) GO TO 15
      II=II+1
      IWA(II)=NODM(K)
   15 CONTINUE
C
      RETURN
C
      END
      SUBROUTINE THEL2(PROP,NDS,TEMPV2,NOD5M)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SOLPM2/ TREF,XINTP,XCON1,XCON2,XPARM1,XPARM2,TOLIL,TOLPC,
     1               TOL1,TOL2,TOL3,TOL4,TOL5,TOL6,TCHK,CRPCON(8),DTMN,
     2               SUBDD,RNGL,RNGU,DTT,TOL7,
     3               KCRP,NITE,NALG,IINTP,NPTS,ITCHK,IST,ISR
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,
     1            ISTAT,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /MATMOD/ STRESS(4),STRAIN(4),C(4,4),IPT,NEL,IPS,ISVE
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /TODIM/ BET,THIC,DE,IEL,NND5,ISOCOR
      COMMON /PORT/ INPORT,JNPORT,NPUTSV,IBPORT,IFPORT,JDC,JVC,JAC
C
      DIMENSION PROP2(3),PROP(16,*),EPST(4),NDS(*),TEMPV2(*),NOD5M(*),
     1          H(9),XDM1(2,9),XDM2(2,2),XDM3(2,1)
C
      EQUIVALENCE (NPAR(5),ITYP2D),(NPAR(3),INDNL)
C
      DATA NGLAST/1000/
C
C
C
      IF(IPT.GT.1) GO TO 5
      TOLMT=1.D-2
      NPTS=INT(PROP(1,5))
      TREF=PROP(2,5)
C
      IINTP=1
      IST=4
      IF (ITYP2D .GE. 2) IST=3
      ISR=3
      IF(ITYP2D.EQ.0) ISR=4
C
      TOLL=TOLMT*ABS(PROP(1,1))
      IF(TOLL.EQ.0.D0) TOLL=TOLMT
      TOLU=TOLMT*ABS(PROP(NPTS,1))
      IF(TOLU.EQ.0.D0) TOLU=TOLMT
C
      RNGL=PROP(1,1) - TOLL
      RNGU=PROP(NPTS,1) + TOLU
C
C
    5 CALL FUNCT2(E1,E2,H,XDM1,NOD5M,XDM2,XDUM,XDM3,NEL,IINTP)
      TEMP2=0.D0
      DO 10 K=1,IEL
      KK=NDS(K)
   10 TEMP2=TEMP2+H(K)*TEMPV2(KK)
C
C
      IF(TEMP2.GE.RNGL) GO TO 35
      write(66,2008)
      STOP
C
   35 L=0
      DO 40 K=2,NPTS
      L=L+1
      DUM=PROP(K,1)
      IF(K.EQ.NPTS) DUM=RNGU
      IF(TEMP2.GT.DUM) GO TO 40
      GO TO 45
   40 CONTINUE
      write(66,2008)
      STOP
C
   45 RATIO=(TEMP2-PROP(L,1))/(PROP(L+1,1)-PROP(L,1))
C
C
      IF(RATIO.GT.1.D0) RATIO=1.D0
      IF(RATIO.LT.0.D0) RATIO=0.D0
C
      DO 60 M=2,4
   60 PROP2(M-1)=PROP(L,M)+RATIO*(PROP(L+1,M)-PROP(L,M))
C
C
      YM2=PROP2(1)
      PV2=PROP2(2)
      A22=YM2/(1.+PV2)
      C12=0.5*A22
      A22=A22/(1.-2.*PV2)
      B22=A22*PV2
      A22=A22-B22
      E12=YM2/(1.-2.*PV2)
      IF (ITYP2D .GE. 2) GO TO 75
C
C
      A12=A22
      B12=B22
      GO TO 80
C
C     PLANE STRESS **
C
   75 A12=YM2/(1.-PV2*PV2)
      B12=A12*PV2
      D12=PV2/(PV2-1.)
      E12=YM2/(1.-PV2)
C
C
   80 EPST(1)=PROP2(3)*(TEMP2-TREF)
      EPST(2)=EPST(1)
      EPST(3)=0.D0
      EPST(4)=EPST(1)
C
C
C     PLANE STRESS **
C
      STRESS(1)=A12*STRAIN(1)+B12*STRAIN(2)-E12*EPST(1)
      STRESS(2)=B12*STRAIN(1)+A12*STRAIN(2)-E12*EPST(2)
      STRESS(3)=C12*STRAIN(3)
      STRESS(4)=0.D0
      IF (ITYP2D .GE. 2) GO TO 90
C
C     PLANE STRAIN **
C
      STRESS(4)=B12*(STRAIN(1)+STRAIN(2))-E12*EPST(4)
      IF(ITYP2D.EQ.1) GO TO 90
C
C     AXISYMMETRIC **
C
      STRESS(1)=STRESS(1)+B12*STRAIN(4)
      STRESS(2)=STRESS(2)+B12*STRAIN(4)
      STRESS(4)=STRESS(4)+A12*STRAIN(4)
C
C
   90 IF(KPRI.EQ.0) GO TO 190
C
C
      IF(ICOUNT.EQ.3) RETURN
C
C           END OF STEP
C
      DO 115 I=1,4
      DO 115 J=1,4
  115 C(I,J)=0.D0
C
C     PLANE STRAIN **
C
      C(1,1)=A12
      C(1,2)=B12
      C(2,1)=B12
      C(2,2)=A12
      C(3,3)=C12
C
      IF(ITYP2D-1) 125,120,130
  120 RETURN
C
C     AXISYMMETRIC **
C
  125 C(1,4)=B12
      C(2,4)=B12
      C(4,1)=B12
      C(4,2)=B12
      C(4,4)=A12
      RETURN
C
C     PLANE STRESS **
C
  130 C(4,1)=B22
      C(4,2)=B22
      C(4,3)=0.D0
      C(4,4)=A22
      RETURN
C
C
C
  190 IF(INDNL.NE.2) GO TO 200
C
C
      CALL CAUCHY
C
  200 IF (IPRI.NE.0 .OR. IPS.EQ.0) GO TO 250
      IF (NG.NE.NGLAST) GO TO 202
      IF(NEL.GT.NELAST) GO TO 206
      IF(IPT-1) 210,208,210
C
  202 NGLAST=NG
  208 IF(ITYP2D-1) 205,205,203
  203 write(66,2002)
      GO TO 206
  205 write(66,2003)
  206 NELAST=NEL
      write(66,2004) NEL
  210 CALL MAXMIN(STRESS,SX,SY,SM)
      IF(ITYP2D-1) 215,215,213
  213 write(66,2005) IPT,(STRESS(I),I=1,3),SX,SY,SM,TEMP2
      GO TO 250
C
  215 write(66,2007) IPT,STRESS(4),(STRESS(I),I=1,3),SX,SY,SM,TEMP2
C
C
  250 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 260
      IF (ISVE.EQ.0) GO TO 260
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-2',NEL,IPT,(STRESS(I),I=1,4),
     2                       (STRAIN(I),I=1,4),TEMP2
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'OUTPUT-2',NEL,IPT,(STRESS(I),I=1,4),
     2                       (STRAIN(I),I=1,4),TEMP2
C
 9000 FORMAT ( A,/,2I10,/,(4E20.13) )
C
C
  260 CONTINUE
      RETURN
C
 2002 FORMAT (1X,7HELEMENT,5X,9HSTRESS-YY,5X,9HSTRESS-ZZ,5X,
     1        9HSTRESS-YZ,4X,10HMAX STRESS,4X,10HMIN STRESS,4X,
     2        5HANGLE,3X,11HTEMPERATURE,/,1X,7HNUM/IPT)
 2003 FORMAT (1X,7HELEMENT,5X,9HSTRESS-XX,5X,9HSTRESS-YY,5X,
     1        9HSTRESS-ZZ,5X,9HSTRESS-YZ,4X,10HMAX STRESS,4X,
     2        10HMIN STRESS,4X,5HANGLE,3X,11HTEMPERATURE,/,1X,
     3        7HNUM/IPT)
 2004 FORMAT(I4,/)
 2005 FORMAT (5X,I3,3E14.6,2E14.6,3X,F6.2,E14.6)
 2007 FORMAT (5X,I3,4E14.6,2E14.6,3X,F6.2,E14.6)
 2008 FORMAT(92H    ERROR   TEMPERATURE OUTSIDE RANGE OF MATERIAL PROPER
     1TY TEMPERATURES   (SUBROUTINE THEL2))
 2009 FORMAT(6X,14HTEMPERATURE = ,E14.6,/)
C
C*FILE END
      END