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