Sources Esope | Cast3M

damatc
C DAMATC    SOURCE    PV090527  24/06/12    21:15:02     11940          
      SUBROUTINE DAMATC(IFOUR,XMAT,DDHOOK,LHOOK,SIG0,VAR0,
     >                  DEPST,DSIGT,EPST0,EPIN0, SIGF,VARF,DEFP,
     >                  XMAT0,DDAUX) 
*--------------------------------------------------------------------*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER   IFOUR, LHOOK, NCRIT
      REAL*8    AVALC, AVALT, BVALC, BVALT, BPCOM, DVOLU,
     .          FTULT, FC01D, FCU1D, GVALT, HLENG,
     .          POISS, RHO,   REDUC, RAT45, STRAU, STRAP, STREP, STRA1,
     .          STRE1, STRA2, STRE2, STRU0, STRIN, SIG0(LHOOK),
     .          TANGB, THR0T, YOUNG, VAR0(*), XMAT(*), XMAT0(*),
     .          DDAUX(LHOOK,LHOOK)
C
      INTEGER   IDIM,  ISTEG, ITERA, IRTD, KCAS,  NDIME, NDIM1, NTEMP
      REAL*8    ATINV, CONS0, CONS1, CONS2, CONS3, CLENG,
     .          DAMAT, DAMAC, DEFP(LHOOK), DEPST(LHOOK), DSIGT(LHOOK),
     .          DYOUN, DDHOOK(LHOOK,LHOOK), EPST0(LHOOK), EPIN0(LHOOK),
     .          EQUIC, EXTCT, EXTTS, EQUIT, EBOU1, GTESP, HBOUN,
     .          PROVI, ROOT2, ROO13, ROO23, SIGF(LHOOK),
     .          THR0C, THRUC, THRET, THREC, VARF(*), ZEROM,
     .          DAMAT0, DAMAC0 
C
      CHARACTER*8 STYPE
C
      PARAMETER (ROOT2 = 1.414213562373090 D0,
     &           ROO13 = 0.577350269189626 D0,
     &           ROO23 = 0.816496580927726 D0,
     &           ZEROM = 1.0 D-10)
C
*--------------------------------------------------------------------*
*     INPUT: XMAT=material parameter
*            XMAT( 1)='YOUN'    ! Young Modulus
*            XMAT( 2)='NU  '    ! Poisson coefficient
*            XMAT( 3)='RHO '    ! Density
*            XMAT( 4)='ALPH'    ! Thermal coefficient
*            XMAT( 5)='HLEN'    ! Effective length
*            XMAT( 6)='GVAL'    ! Fracture energy
*            XMAT( 7)='FTUL'    ! Tensile stress
*            XMAT( 8)='REDC'    ! Drop Factor for Peak Tensile Stress
*            XMAT( 9)='FC01'    ! Elastic Limit Compressive Stress
*            XMAT(10)='RT45'    ! Equi-biaxial Compressive Ratio
*            XMAT(11)='FCU1'    ! Compressive Peak Stress
*            XMAT(12)='STRU'    ! Ultimate Limit Strain
*            XMAT(13)='EXTP'    ! Reference Strain for Plastic Parameter
*            XMAT(14)='STRP'    ! Reference Stress for Plastic Parameter
*            XMAT(15)='EXT1'    ! Fitting Point 1 (Strain)
*            XMAT(16)='STR1'    ! Fitting Point 1 (Stress)
*            XMAT(17)='EXT2'    ! Fitting Point 2 (Strain)
*            XMAT(18)='STR2'    ! Fitting Point 2 (Stress)
*            XMAT(19)='NCRI'    ! Tensile Softening Criteria (1=Expon./2=Lin.)
*
*            XMAT(20)='TANB'    ! Drucker-Prager Tangent
*            XMAT(21)='BPLC'    ! Plastic Compressive Parameter
*            XMAT(22)='TR0T'    ! Initial Tension     Damage Threshold
*            XMAT(23)='TR0C'    ! Initial Compression Damage Threshold
*            XMAT(24)='TRUC'    ! Ultimate Compression Damage Threshold
*            XMAT(25)='AVAT'    ! "A" Parameter for Tension
*            XMAT(26)='BVAT'    ! "B" Parameter for Tension
*            XMAT(27)='AVAC'    ! "A" Parameter for Compression
*            XMAT(28)='BVAC'    ! "B" Parameter for Compression
*            XMAT(29)='DVOL'    ! Volume of a finite element
*
*            DDHOOK=Hooke matrix
*
*            VAR0: internal variable
*            VAR0(1)='THRT'     ! Tension theshold
*            VAR0(2)='EQUT'     ! Tension equivalent strain
*            VAR0(3)='DAMT'     ! Tension damage
*            VAR0(4)='THRC'     ! Compre. theshold
*            VAR0(5)='EQUC'     ! Compre. equivalent strain
*            VAR0(6)='DAMC'     ! Compre. damage
*--------------------------------------------------------------------*
C
      IDIM=2
      IF(IFOUR.EQ.2)IDIM=3
      NDIME = IDIM
      NDIM1 = NDIME + 1
C
C     COMPULSORY CARACTERISTICS
C
      YOUNG = XMAT(1)
      POISS = XMAT(2)
      RHO   = XMAT(3)
C     ALPHA = XMAT(4)
C
      HLENG = XMAT(5)
      GVALT = XMAT(6)
      FTULT = XMAT(7)
      REDUC = XMAT(8)
      FC01D = XMAT(9)
      RAT45 = XMAT(10)
      FCU1D = XMAT(11)
      STRAU = XMAT(12)
      STRAP = XMAT(13)
      STREP = XMAT(14)
      STRA1 = XMAT(15)
      STRE1 = XMAT(16)
      STRA2 = XMAT(17)
      STRE2 = XMAT(18)
      NCRIT = nint(XMAT(19))
C
      DVOLU = XMAT(29)
C
C     FACULTATIVE CARACTERISTICS
C
      DYOUN = SQRT (YOUNG)
      THR0T = ABS(FTULT/DYOUN)
C
C=====----------------------------------------------------------------
C  EVALUATING PRELIMINAR QUANTITIES FOR PLASTIC-VISCOUS-DAMAGE MODEL |
C  (INITIAL AND FINAL NONLINEAR THRESHOLDS AND PLASTIC PARAMETER)    |
C=====----------------------------------------------------------------
C
      CONS0 = FC01D / YOUNG
C
      r_z = 1.D0 / (1.0 D0 - 2.0 D0 * RAT45)
      TANGB = ROOT2 * (1.0D0 - RAT45) * r_z
      PROVI = ROO23 * RAT45 * r_z
 
      IF(PROVI*FC01D.LT.0.D0) THEN
        write (*,'(/,''*** PROVI*FC01D is negative ***'')')
        CALL ERREUR(5)
        RETURN
      ENDIF
 
      THR0C = SQRT (PROVI*FC01D)
C
      BPCOM = YOUNG * FCU1D * (STRAP - STREP / YOUNG) /
     .       ((STRAP - CONS0) * (STREP + FCU1D))
C
      CONS1 = YOUNG - BPCOM
      CONS2 = BPCOM * CONS0
      CONS3 = ROO13 * TANGB - ROO23
C
      r_z = CONS1 * STRAU + CONS2
      THRUC = SQRT (ABS (PROVI*r_z))
C
C=====----------------------------------------------------------------
C  EVALUATION OF PARAMETERS AVALC & BVALC & ...                      |
C=====----------------------------------------------------------------
C
      CALL DAM_AB (AVALC, BVALC, CONS1, CONS2, CONS3, EBOU1,
     .             HBOUN, STRA1, STRA2, STRE1, STRE2, THR0C)
C
      EBOU1 = 1.0 D+33
C
C=====================================================================
C   DEFINITION OF L* AND TENSION SOFTENING PARAMETERS.               =
C=====================================================================
C
C=====----------------------------------------------------------------
C  CHARACTERISTIC LENGTH, IF NOT PREVIOUSLY SPECIFIED OR EVALUATED   |
C=====----------------------------------------------------------------
C
      IF (HLENG .EQ. 0.0 D0)  THEN
         IF (DVOLU .EQ. 0.0 D0)  THEN
            write (*,'(/,''*** DVOLU or HLENG is unidentified  ***'')')
            CALL ERREUR(5)
            RETURN
         ENDIF
         CLENG = DVOLU**(1.0D0/FLOAT(NDIME))
      ELSE
         CLENG = HLENG
      ENDIF
C
      GTESP = GVALT / CLENG
C
      IF (NCRIT .EQ. 1 ) THEN
C
         ATINV = GTESP / (THR0T **2) - 0.5 D0
         IF (ATINV .LE. 0.0 D0)  ATINV = ZEROM
         AVALT = 1.0 D0 / ATINV
      ELSE IF (NCRIT .EQ. 2)  THEN
C
         EXTCT = FTULT  / YOUNG
         EXTTS = 1.0 D0 / REDUC * (2.0 D0 * GTESP / FTULT -
     .                            (1.0 D0 - REDUC) * EXTCT  )
         IF (EXTTS .LE. EXTCT)  EXTTS = EXTCT + ZEROM
         BVALT = REDUC * EXTCT / (EXTTS - EXTCT)
         AVALT = BVALT * EXTTS * DYOUN
      ELSE
         write(*,'(/,''*** NCRIT : valeur incorrecte  ***'')')
         CALL ERREUR(5)
      RETURN
C
      END IF
C
C      calcul de la matrice elastique
      STYPE = 'ISOTROPE'
      KCAS=1
      IRTD=1
      CALL DOHMAS (XMAT,STYPE,IFOUR,LHOOK,KCAS,DDHOOK,IRTD)
      IF (IRTD.EQ.0) THEN
        CALL ERREUR(5)
        RETURN
      ENDIF
C
C      calcul de la matrice elastique debut de pas
      IRTD=1
      CALL DOHMAS (XMAT0,STYPE,IFOUR,LHOOK,KCAS,DDAUX,IRTD)
      IF (IRTD.EQ.0) THEN
        CALL ERREUR(5)
        RETURN
      ENDIF
C
       IF (IFOUR.EQ.-2) THEN
        STYPE = 'PSTRES'
       ELSEIF ((IFOUR.EQ.-1).OR.(IFOUR.EQ.-3)) THEN
        STYPE = 'PSTRAI'
       ELSEIF (IFOUR.EQ. 0) THEN
        STYPE = 'AXISYM'
       ENDIF
C
C     ENFORCE NO INITIAL STRESSES DUE TO TEMPERATURE.....
         NTEMP = 0
         ISTEG = 0
         ITERA = 0
         STRIN = 0.D0
C
      THRET = VAR0(1)
      EQUIT = VAR0(2)
      DAMAT = VAR0(3)
      THREC = VAR0(4)
      EQUIC = VAR0(5)
      DAMAC = VAR0(6)
 
      DAMAT0 = DAMAT
      DAMAC0 = DAMAC
C
      IF (THRET.EQ.0.) THRET = THR0T
      IF (THREC.EQ.0.) THREC = THR0C
 
*  AM  16/7/19
 
      IF (THRET.LT.THR0T) THRET = THR0T
      IF (THREC.LT.THR0C) THREC = THR0C
 
 
C
      do IE1=1,LHOOK
        EPST0(IE1)=EPST0(IE1)-EPIN0(IE1)
      enddo
 
 
 
*AM*  CALL MATVE1 (DDHOOK,EPST0,LHOOK,LHOOK,SIG0,1)
      CALL MATVE1 (DDAUX,EPST0,LHOOK,LHOOK,SIG0,1)
 
C
      CALL DAMSTC (ISTEG, ITERA,
     .             NCRIT, NDIME, NDIM1,  LHOOK, LHOOK, NTEMP,
     .             AVALC, AVALT, BPCOM,  BVALC, BVALT,
     .             DAMAC, DAMAT, DDHOOK, DEPST, DSIGT, EBOU1,
     .             EQUIC, EQUIT,
     .             POISS, SIG0 , SIGF ,  STRIN, TANGB,
     .             THREC, THRET, THRUC,  THR0C, THR0T, YOUNG,
     .             STYPE, DAMAT0, DAMAC0)
C
      VARF(1) = THRET
      VARF(2) = EQUIT
      VARF(3) = DAMAT
      VARF(4) = THREC
      VARF(5) = EQUIC
      VARF(6) = DAMAC
C
      IF (IFOUR.EQ.-2)THEN
        CALL INVALM(DDHOOK,LHOOK,2,IRTD,1.d-8)
        DDHOOK(4,4)=1.D0/DDHOOK(4,4)
      ELSE
        CALL INVALM(DDHOOK,LHOOK,LHOOK,IRTD,1.d-8)
      ENDIF
c
      CALL MATVE1(DDHOOK,SIG0,LHOOK,LHOOK,DEFP,1)
c
      do IE1=1,LHOOK
        DEFP(IE1)=EPST0(IE1)+DEPST(IE1)-DEFP(IE1)
      enddo
c
      RETURN
      END