!DEC$ FREEFORM ! !This code is distributed under CC BY 4.0 license ! !Copyright 2017 Andreas Seupel, Geralf Hütter, Meinhard Kuna ! !Institute of Mechanics and Fluid Dynamics, TU Bergakademie Freiberg ! !ABAQUS-Version 6.14-2 !USED FORTRAN COMPILER: Fortran Intel(R) 64, Version 11.0 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! MODULE Abaqus_Interface INCLUDE 'ABA_PARAM.INC' PRIVATE INTEGER,PARAMETER,PUBLIC::realkind=KIND(r),intkind=KIND(i) END MODULE ! ! *************************************************************************** ! SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& & RPL,DDSDDT,DRPLDE,DRPLDT,& & STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,& & NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,& & CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC) ! -MATERIAL ROUTINE FOR LOCAL OR NONLOCAL (GRADIENT-ENHANCED) DAMAGE ! -DAMAGE DRIVEN BY EQUIVALENT PLASTIC STRAIN ! -DECREASE OF STRAIN HARDENING ! -THE UMAT CAN BE USED IN 2D PLANE STRAIN, 2D AXISYMMETRIC AND 3D CONTINUUM SIMULATIONS !!!!! MATERIAL PARAMETERS TO BE DEFINED PROPS (NPROPS=15) ! ELASTICITY ! YOUNGS MODULUS ! E =PROPS(1) ! POISSONS RATIO ! NU =PROPS(2) ! ! STRAIN HARDENING: R =TAU_0+H*RR**N+RINF*(ONE_R-EXP(-A*RR)); RR=HARDENING VARIABLE ! INITIAL YIELD STRESS ! TAU_0 =PROPS(3) ! PREFACTOR HOLLOMON ! H =PROPS(4) ! EXPONENT HOLLOMON ! N =PROPS(5) ! SATURATION VOCE ! RINF =PROPS(6) ! PARAMETER VOCE ! A =PROPS(7) ! ! JOHNSON-COOK DAMAGE LOCUS: E_CI =C_2*exp(-C_3*ETA)+C_1; ETA=STRESS TRIAXIALITY ! C_1 =PROPS(8) ! C_2 =PROPS(9) ! C_3 =PROPS(10) ! C_4 =PROPS(11) ! DAMAGE PARAMETER ! S_1 =PROPS(12) ! E_C =E_1 IN NONLOCAL CASE ! E_1 =PROPS(13) ! DECISION PARAMETER, WHETER COMPUTATION IS INTERRUPTED AT INCIPIENT CRACK INITIATION (SET DDS>0.5) ! DDS =PROPS(14) ! LOCAL (DAMAGE_CASE<0.5) OR NONLOCAL (DAMAGE_CASE<0.5) DAMAGE CASE ! DAMAGE_CASE =PROPS(15) ! !!!!! STATEV (INTERNAL AND STATE VARIABLES), NSTATEV=12 ! STATEV(1) =EQUIVALENT PLASTIC STRAIN ! STATEV(2) =HARDENING VARIABLE=EQUIVALENT PLASTIC STRAIN ! STATEV(3) =DAMAGE PARAMETER ! STATEV(4) =INDICATOR FOR DAMAGE SOURCE ACTIVATION, I. E. EQUIVALENT PLASTIC STRAIN>E_I ! STATEV(5) =E_I^FIX, STRAIN AT SOURCE ACTIVATION ! STATEV(6) =EQUIVALENT DEVIATORIC STRAIN FROM WHOLE STRAIN INCREMENT (FOR POSTPROCESSING) ! STATEV(7) =MISES EQUIVALENT STRESS ! STATEV(8) =KAPPA ! STATEV(9) =STRESS_TRIAXILITY ! STATEV(10) =HEAT SOURCE/SOURCE TERM R ! STATEV(11) =DERIVATIVE OF THE HEAT SOURCE/SOURCE TERM R WITH RESPECT TO TEMPERATURE ! STATEV(12) =VARIABLE ACTING AS DELETE INDICATOR: 1.0-NOTHING HAPPENS, 0.0-ELEMENT CAN BE DELETED USE Abaqus_Interface IMPLICIT NONE INTEGER(KIND=intkind),INTENT(IN) ::NTENS,NSTATV,NDI,NSHR,NOEL,NPT,LAYER INTEGER(KIND=intkind),INTENT(IN) ::KSPT,KINC INTEGER(KIND=intkind),DIMENSION(4),INTENT(IN) ::JSTEP ! REAL(KIND=realkind),INTENT(IN) ::PREDEF,DPRED,DTIME,TEMP REAL(KIND=realkind),INTENT(IN) ::DTEMP,CELENT ! REAL(KIND=realkind),INTENT(INOUT) ::PNEWDT REAL(KIND=realkind),INTENT(INOUT) ::SSE,SPD,SCD ! REAL(KIND=realkind) ::RPL,DRPLDT ! REAL(KIND=realkind),DIMENSION(2),INTENT(IN) ::TIME REAL(KIND=realkind),DIMENSION(3),INTENT(IN) ::COORDS REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::STRAN,DSTRAN REAL(KIND=realkind),DIMENSION(3,3),INTENT(IN) ::DROT,DFGRD0,DFGRD1 ! REAL(KIND=realkind),DIMENSION(NSTATV),INTENT(INOUT) ::STATEV REAL(KIND=realkind),DIMENSION(NTENS),INTENT(INOUT) ::STRESS ! REAL(KIND=realkind),DIMENSION(NTENS,NTENS),INTENT(OUT) ::DDSDDE REAL(KIND=realkind),DIMENSION(NTENS),INTENT(OUT) ::DDSDDT,DRPLDE CHARACTER(80)::CMNAME !!! MATERIAL PARAMETERS ! ELASTICITY REAL(KIND=realkind) ::E,NU,K,G ! STRAIN HARDENING REAL(KIND=realkind) ::TAU_0,H,A,RINF,N ! DAMAGE REAL(KIND=realkind) ::C_1,C_2,C_3,C_4,S_1,E_1 ! SIZE PARAMETERS INTEGER(KIND=intkind),INTENT(IN) ::NPROPS REAL(KIND=realkind),DIMENSION(NPROPS),INTENT(IN) ::PROPS ! DAMAGE CASE: LOCAL - NON LOCAL REAL(KIND=realkind) ::DAMAGE_CASE ! FINISH/ABORT CALCULATION AT CRACK INITIATION REAL(KIND=realkind) ::DDS ! INTEGER(KIND=intkind) ::PNEW,SUBSTEPCONTROL ! FLOAT NUMBERS REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! REAL(KIND=realkind),PARAMETER ::PI=3.141592653589793_realkind !!! DAMAGE RELATED VARIABLES REAL(KIND=realkind) ::E_EQ_WHOLE,DELTA_E_EQ_WHOLE,E_C,A_D,B_D,C_D ! ARRAYS REAL(KIND=realkind),DIMENSION(NTENS,1) ::DEVIATOR_EPS_MATRIX,DJAC_DEPS,STRESS_MATRIX REAL(KIND=realkind),DIMENSION(NTENS) ::DEVIATOR_EPS INTEGER(KIND=intkind) ::COUNTER ! ITERATION VARIABLES INTEGER(KIND=intkind) ::i,j,LOCAL_NLOCAL ! SCALARS REAL(KIND=realkind) ::JAC,EPSILON_MISES,RR,Z,DAMAGE,DAMAGE_SWITCH,T,THETA,DELTATHETA,DELTATHETAMIN,I1T,J2T,J3T,COS3THETA,TAU_MISES_T,TAU_HYD_T,TAU_MISES,TAU_HYD,DG_TR_DI1,NULLTEST,E_CI REAL(KIND=realkind) ::TAUM_0,R,HH,YIELD_FUNCTION,DELTA_LAMBDA_PL,DELTA_RR,DELTA_E_D_BAR,DTIME_INK, EPS_D_BAR,E_I,DDAMAGE_D_DLAMBDA_PL,DRR_DDELTA_LAMBDA_PL ! MATRICES AND VECTORS REAL(KIND=realkind),DIMENSION(NTENS) ::EPS,KRONECKER,TAU_ROT,EPSGES,DEVIATORT,N_TR,TAU_T REAL(KIND=realkind),DIMENSION(3) ::STARTDELTA REAL(KIND=realkind),DIMENSION(NTENS,NTENS) ::CMATRIX,P_SCALING_MATRIX,IMATRIX,IDEV,DDSDDE2 REAL(KIND=realkind),DIMENSION(NTENS,NTENS) ::DTAU_ROT_DEPS !!!!! MATERIAL PARAMETERS FROM PROPS ! ELASTICITY ! YOUNGS MODULUS E =PROPS(1) ! POISSONS RATIO NU =PROPS(2) ! STRAIN HARDENING ! INITIAL YIELD STRESS TAU_0 =PROPS(3) ! PREFACTOR HOLLOMON H =PROPS(4) ! EXPONENT HOLLOMON N =PROPS(5) ! SATURATION VOCE RINF =PROPS(6) ! PARAMETER VOCE A =PROPS(7) ! JOHNSON-COOK DAMAGE LOCUS C_1 =PROPS(8) C_2 =PROPS(9) C_3 =PROPS(10) C_4 =PROPS(11) ! DAMAGE PARAMETER S_1 =PROPS(12) ! E_C =E_1 IN NONLOCAL CASE E_1 =PROPS(13) DDS =PROPS(14) DAMAGE_CASE =PROPS(15) ! DECISION BETWEEN LOCAL - NONLOCAL IF (DAMAGE_CASE<0.5) THEN LOCAL_NLOCAL =0 ELSE LOCAL_NLOCAL =1 END IF ! COMPRESSION AND SHEAR MODULUS ! K = E/(THREE_R*(ONE_R-TWO_R*NU)) G = E/(TWO_R*(ONE_R+NU)) PNEW =0 ! TEST NUMBER OF MATERIAL PARAMETERS IF (NPROPS==15) THEN ELSE PNEW =1 IF (NOEL==1 .AND. NPT==1 .AND. KINC==1) THEN WRITE(6,*) "***FATAL ERROR: INPUTERROR! NUMBER OF MATERIAL PARAMETERS IS WRONG! CHECK THE PROPERTY MODUL INPUT!" WRITE(6,*) "NUMBER ALLOWED: 15" END IF END IF ! TEST NUMBER OF STATEV VARIABLES IF (NTENS==6) THEN IF (NSTATV==12) THEN ELSE PNEW =1 IF (NOEL==1 .AND. NPT==1 .AND. KINC==1) THEN WRITE(6,*) "***FATAL ERROR: INPUTERROR! NUMBER OF SDV WRONG! CHECK THE PROPERTY MODUL INPUT!" WRITE(6,*) "NUMBER ALLOWED FOR 3D ANALYSIS: 12" END IF END IF ELSEIF (NTENS==4) THEN IF (NSTATV==12) THEN ELSE PNEW =1 IF (NOEL==1 .AND. NPT==1 .AND. KINC==1) THEN WRITE(6,*) "***FATAL ERROR: INPUTERROR! NUMBER OF SDV WRONG! CHECK THE PROPERTY MODUL INPUT!" WRITE(6,*) "NUMBER ALLOWED FOR 2D/AXISYMMETRIC ANALYSIS: 12" END IF END IF END IF IF (NTENS>3) THEN ELSE PNEW =1 IF (NOEL==1 .AND. NPT==1 .AND. KINC==1) THEN WRITE(6,*) "***FATAL ERROR: Plane Stress State is not supported by the UMAT-subroutine!" END IF END IF ! INITALIZE STATEV IF (KINC==1 .AND. JSTEP(1)==1) THEN DO i=1,NSTATV,1 STATEV(i)=ZERO_R END DO SPD =ZERO_R STATEV(12) =ONE_R STATEV(5) =10000.0 END IF ! DECLARE INTERNAL VARIABLES ! EQUIVALENT PLASTIC STRAIN EPSILON_MISES =STATEV(1) ! HARDENING VARIABLE = QUIVALENT PLASTIC STRAIN RR =STATEV(2) ! DAMAGE VARIABLE DAMAGE =STATEV(3) DAMAGE_SWITCH =STATEV(4) ! EQUIVALENT STRAIN CALCULATED FROM WHOLE STRAIN INCREMENT E_EQ_WHOLE =STATEV(6) ! KRONECKER-DELTA UNITY TENSOR AS VECTOR KRONECKER =ZERO_R DO i=1,NDI,1 KRONECKER(i) =ONE_R END DO ! ! ENTRY OF MATERIAL TANGENT CAPITAL K JAC =EXP(DOT_PRODUCT(STRAN+DSTRAN,KRONECKER)) DJAC_DEPS(1:NTENS,1) =JAC*KRONECKER(1:NTENS) ! REWRITE STRAIN IN SCIENTIFIC FORMAT EPS =DSTRAN DO i=NDI+1,NTENS,1 EPS(i) =EPS(i)/TWO_R END DO DTAU_ROT_DEPS =ZERO_R ! ! ! STANDARD VECTORS/MATRICES ! ISOTROPIC ELASTICITY TENSOR AS MATRIX CMATRIX CMATRIX =ZERO_R IF (NTENS==3) THEN DO i=1,NTENS,1 DO j=1,NTENS,1 IF (i<=NDI .AND. j<=NDI .AND. i/=j) THEN CMATRIX(i,j) =(E/(ONE_R-NU**2)*(NU)) ELSE CMATRIX(i,j) =ZERO_R END IF END DO IF (i<=NDI) THEN CMATRIX(i,i) =E/(ONE_R-NU**2) ELSE CMATRIX(i,i) =TWO_R*G END IF END DO ELSE DO i=1,NTENS,1 DO j=1,NTENS,1 IF (i<=NDI .AND. j<=NDI .AND. i/=j) THEN CMATRIX(i,j) =(E/(ONE_R+NU)*(NU/(ONE_R-TWO_R*NU))) ELSE CMATRIX(i,j) =ZERO_R END IF END DO IF (i<=NDI) THEN CMATRIX(i,i) =(E/(ONE_R+NU)*(ONE_R+NU/(ONE_R-TWO_R*NU))) ELSE CMATRIX(i,i) =TWO_R*G END IF END DO END IF ! SCALING MATRIX P_SCALING_MATRIX =ZERO_R IF (NTENS==3) THEN P_SCALING_MATRIX(1,1) =TWO_R P_SCALING_MATRIX(2,2) =TWO_R P_SCALING_MATRIX(1,2) =ONE_R P_SCALING_MATRIX(2,1) =ONE_R P_SCALING_MATRIX(3,3) =TWO_R ELSE DO i=1,NTENS,1 DO j=1,NTENS,1 P_SCALING_MATRIX(i,j) =ZERO_R END DO END DO ! DO i=1,NTENS,1 IF (i<=NDI) THEN P_SCALING_MATRIX(i,i) =ONE_R ELSE P_SCALING_MATRIX(i,i) =TWO_R END IF END DO END IF ! IDENTITY MATRIX IMATRIX =ZERO_R DO i=1,NTENS,1 IMATRIX(i,i) =ONE_R END DO ! DEVIATORIC OPERATOR CALL VECTORPRODUCT (NDI,NTENS,KRONECKER,KRONECKER,IDEV) IDEV =IMATRIX-ONE_R/THREE_R*IDEV !!!!! ELASTIC TRIAL STEP ! CALL STRAIN_HARDENING (PROPS, NPROPS,RR,1.0,1.0,R,HH,E_CI) TAU_ROT =STRESS NULLTEST=ZERO_R DO i=1,NTENS,1 NULLTEST =NULLTEST+EPS(i) END DO ! TEST ZERO STRAIN CASE IF (ABS(NULLTEST)<=EPSILON(ZERO_R) .AND. KINC>1) THEN STRESS =TAU_ROT DDSDDE =MATMUL(IMATRIX,(DTAU_ROT_DEPS+CMATRIX)) ! REWRITE ABAQUS-FORMAT DO i=1,NTENS,1 DO j=4,NTENS,1 DDSDDE(i,j) =DDSDDE(i,j)*ONE_R/TWO_R END DO END DO DDSDDT =ZERO_R DRPLDE =ZERO_R DRPLDT =ZERO_R RPL =ZERO_R ! ELSE ! INTERNAL SUBSTEPPING (DEACTIVATED BY DEFAULT) THETA =ZERO_R DELTATHETA =ONE_R DELTATHETAMIN =ONE_R SUBSTEPCONTROL =1 EPSGES =EPS STARTDELTA =ZERO_R ! LOOP FOR SUBSTEPPING DO WHILE (SUBSTEPCONTROL==1 .AND. PNEW==0) ! STRAIN INCREMENT, TIME INCREMENT, NONLOCAL VARIABLE EPS =EPSGES*(THETA+DELTATHETA) DTIME_INK =DTIME*(THETA+DELTATHETA) DELTA_E_D_BAR =DTEMP*(THETA+DELTATHETA) IF (LOCAL_NLOCAL==1) THEN !NONLOCAL CASE EPS_D_BAR =TEMP+DELTA_E_D_BAR E_I =ZERO_R E_C =E_I+E_1 ! CALCULATE DAMAGE: MEDIAVILLA DAMAGE LAW IF (EPS_D_BAR>=E_C) THEN DAMAGE =0.99_realkind DDAMAGE_D_DLAMBDA_PL =ZERO_R DRR_DDELTA_LAMBDA_PL =ZERO_R ELSEIF (EPS_D_BAR<=ZERO_R .OR. EPS_D_BAR2.5) THEN DAMAGE =STATEV(3) DDAMAGE_D_DLAMBDA_PL =ZERO_R DRR_DDELTA_LAMBDA_PL =ZERO_R ELSE A_D =THREE_R/TANH(THREE_R)/(E_C-E_I) B_D =TANH(THREE_R) C_D =TWO_R*A_D*B_D*E_I+atanh(B_D) C_D =C_D/TWO_R/B_D DAMAGE =0.99_realkind*ONE_R/TWO_R*(ONE_R+TANH(TWO_R*(A_D*B_D*(EPS_D_BAR)-B_D*C_D))/B_D) DDAMAGE_D_DLAMBDA_PL =0.99_realkind*A_D*(ONE_R/COSH(TWO_R*(A_D*B_D*(EPS_D_BAR)-B_D*C_D)))**2 DRR_DDELTA_LAMBDA_PL =-DELTA_LAMBDA_PL*DDAMAGE_D_DLAMBDA_PL END IF END IF TAU_T =TAU_ROT+MATMUL(CMATRIX,EPS) ! INVARIANTS I1, J2; DEVIATOR OF STRESS AND FLOW DIRECTION (N_TR) FROM TRIAL STRESS I1T =DOT_PRODUCT(TAU_T,KRONECKER) DEVIATORT =TAU_T-ONE_R/THREE_R*KRONECKER*I1T J2T =ONE_R/TWO_R*DOT_PRODUCT(DEVIATORT,MATMUL(P_SCALING_MATRIX,DEVIATORT)) CALL J3CALC (NTENS,DEVIATORT,DEVIATORT,DEVIATORT,J3T) ! COS3THETA IS THE LODE PARAMETER IF (J2T==ZERO_R) THEN COS3THETA =ONE_R TAU_MISES_T =sqrt(THREE_R*J2T) N_TR =ZERO_R ELSE COS3THETA =THREE_R*sqrt(THREE_R)/TWO_R*J3T/J2T**(THREE_R/TWO_R) TAU_MISES_T =sqrt(THREE_R*J2T) N_TR =THREE_R/TWO_R/TAU_MISES_T*DEVIATORT END IF ! HYDROSTATIC STRESS TAU_HYD_T =I1T/THREE_R ! DEVIATOR FROM WHOLE STRAIN INCREMENT (JUST FOR POSTPROCESSING) DEVIATOR_EPS=MATMUL(IDEV,EPS) DEVIATOR_EPS_MATRIX(1:NTENS,1) =DEVIATOR_EPS(1:NTENS) DELTA_E_EQ_WHOLE =SQRT(FOUR_R/THREE_R*(ONE_R/TWO_R*DOT_PRODUCT(DEVIATOR_EPS,MATMUL(P_SCALING_MATRIX,DEVIATOR_EPS)))) !!! ELASTIC-PLASTIC OPERATOR SPLIT CALL STRAIN_HARDENING (PROPS, NPROPS,RR,1.0,1.0,R,HH,E_CI) ! EVALUATE YIELD FUNCTION YIELD_FUNCTION =TAU_MISES_T-R*(ONE_R-DAMAGE) !!! CASE DIFFERENTIATION !!! ! TEST ELASTICITY IF (SIGN(ONE_R,YIELD_FUNCTION)0.5) THEN DDSDDE =MATMUL(DDSDDE,(DTAU_ROT_DEPS+CMATRIX)) DDSDDT =ZERO_R ELSE DDSDDE =MATMUL(DDSDDE,(DTAU_ROT_DEPS+CMATRIX)) DDSDDT =ZERO_R END IF STARTDELTA(1) =DELTA_LAMBDA_PL DRPLDE =ZERO_R DRPLDT =ZERO_R ELSE DDSDDE =CMATRIX DDSDDT =ZERO_R STRESS =TAU_ROT DRPLDE =ZERO_R DRPLDT =ZERO_R END IF ELSE !NONLOCAL CASE CALL PLASTICITY_MODUL_NL(DELTA_E_D_BAR,TEMP+DELTA_E_D_BAR,DELTA_E_EQ_WHOLE,NTENS,NPROPS,NSTATV,STARTDELTA, TAU_T,DEVIATORT,N_TR,CMATRIX,KRONECKER,IMATRIX,IDEV,P_SCALING_MATRIX,PROPS,STATEV,TAU_MISES_T,& & STRESS,DDSDDE,DDSDDT,DRPLDE,DRPLDT,DELTA_LAMBDA_PL,DELTA_RR,DAMAGE,PNEW) IF (PNEW==0) THEN DDSDDE =MATMUL(DDSDDE,(DTAU_ROT_DEPS+CMATRIX)) DRPLDE =MATMUL((DTAU_ROT_DEPS+CMATRIX),DRPLDE) DDSDDT =DDSDDT STARTDELTA(1) =DELTA_LAMBDA_PL ELSE DDSDDE =CMATRIX DDSDDT =ZERO_R DRPLDE =ZERO_R DRPLDT =ZERO_R STRESS =TAU_ROT END IF END IF ! END CASE LOCAL NONLOCAL END IF ! END ELASTIC PLASTIC !! CONVERGENCE SUB-STEPPING IF (THETA+DELTATHETA-ONE_R==ZERO_R .AND. PNEW==0) THEN SUBSTEPCONTROL =0 ELSE IF (PNEW==1) THEN DELTATHETA =DELTATHETA/TWO_R IF (DELTATHETA=ZERO_R) THEN DELTATHETA =ONE_R-THETA END IF END IF END IF ! END CONVERGENCE SUB-STEPPING END DO ! END SUBSTEPPING ! TEST NAN IF (PNEW==0) THEN DO i=1,NTENS,1 IF (STRESS(i).NE. STRESS(i)) THEN IF (PNEW==0) THEN PNEW=1 write(6,*) "NAN" write(6,*) "STRESS" !write(6,*) STRESS !write(6,*) SWITCH END IF END IF END DO END IF IF (PNEW==0) THEN IF (DELTA_LAMBDA_PL .NE. DELTA_LAMBDA_PL) THEN PNEW=1 write(6,*) "NAN" write(6,*) "DELTA_LAMBDA_PL" !write(6,*) DELTA_LAMBDA_PL !write(6,*) SWITCH END IF END IF ! ACTIVATE END AT CRACK INITATION IF (STATEV(3)>0.989_realkind .AND. DDS>0.5) THEN PNEW=1 write(6,*) "*** SIMULATION INTERRUPTED BY UMAT: INITIATION OF AN INCIPIENT CRACK" write(7,*) "*** SIMULATION INTERRUPTED BY UMAT: INITIATION OF AN INCIPIENT CRACK" END IF IF (PNEW==0) THEN ! GLOBAL SUBSTEPPING CONTROL OF ABAQUS PNEWDT =1.25_realkind ! UPDATE STATEV STATEV(1) =EPSILON_MISES+DELTA_LAMBDA_PL STATEV(2) =RR+DELTA_RR CALL STRAIN_HARDENING (PROPS, NPROPS,RR+DELTA_RR,1.0,1.0,R,HH,E_CI) TAU_MISES =TAU_MISES_T-THREE_R*G*DELTA_LAMBDA_PL TAU_MISES =abs(TAU_MISES) ! UPDATE KAPPA = STATEV(8) IF (STATEV(8)>TEMP+DTEMP) THEN ELSE STATEV(8) =TEMP+DTEMP END IF STATEV(3) =DAMAGE TAU_HYD =TAU_HYD_T T =TAU_HYD/TAU_MISES IF (TAU_MISES > EPSILON(ZERO_R))THEN ELSE T =1.0E8 END IF STATEV(7) =TAU_MISES STATEV(9) =T ! UPDATE NONLOCAL SOURCE TERM IF (LOCAL_NLOCAL==1) THEN IF (DAMAGE_SWITCH>0.5) THEN IF (DAMAGE_SWITCH>2.5) THEN DRPLDT =ZERO_R DRPLDE =ZERO_R ELSE STATEV(10) =STATEV(10)+DELTA_LAMBDA_PL END IF STATEV(6) =STATEV(6)+DELTA_E_EQ_WHOLE END IF ELSE IF (SIGN(ONE_R,YIELD_FUNCTION)>ZERO_R) THEN STATEV(10) =ZERO_R STATEV(6) =STATEV(6)+DELTA_E_EQ_WHOLE END IF END IF ! TEST DAMAGE INITIATION CALL STRAIN_HARDENING (PROPS, NPROPS,RR+DELTA_RR,T,COS3THETA,R,HH,E_CI) IF (LOCAL_NLOCAL==1) THEN IF (STATEV(1)>E_CI .AND. DAMAGE_SWITCH<0.5) THEN DAMAGE_SWITCH =ONE_R STATEV(5) =STATEV(1) END IF STATEV(4) =DAMAGE_SWITCH ELSE IF (STATEV(1)>E_CI .AND. DAMAGE_SWITCH<0.5) THEN !E_CI DAMAGE_SWITCH =ONE_R STATEV(5) =STATEV(1) END IF STATEV(4) =DAMAGE_SWITCH END IF ! CALCULATE STRESS AND STRESS TANGENT STRESS_MATRIX(1:NTENS,1) =STRESS DDSDDE =(JAC*DDSDDE+MATMUL(STRESS_MATRIX,TRANSPOSE(DJAC_DEPS)))/JAC DDSDDT =(DDSDDT) ! RESET TANGENT ENTRIES STATEV(11) =DRPLDT DRPLDT =ZERO_R RPL =STATEV(10) ! SIGNAL FOR ELEMENT EROSION OR TOTAL FAILURE IF (DAMAGE>=0.989_REALKIND .AND. DAMAGE_SWITCH<2.5) THEN DAMAGE_SWITCH =THREE_R STATEV(4) =DAMAGE_SWITCH END IF STATEV(6) =TAU_HYD ! REWRITE INTO ABAQUS FORMAT IF (NTENS==3) THEN DO i=1,NTENS,1 DO j=3,NTENS,1 DDSDDE(i,j) =DDSDDE(i,j)*ONE_R/TWO_R END DO END DO DO j=3,NTENS,1 DRPLDE(j) =DRPLDE(j)*ONE_R/TWO_R END DO ELSE DO i=1,NTENS,1 DO j=4,NTENS,1 DDSDDE(i,j) =DDSDDE(i,j)*ONE_R/TWO_R END DO END DO DO j=4,NTENS,1 DRPLDE(j) =DRPLDE(j)*ONE_R/TWO_R END DO END IF ELSE ! CASE OF STEPSIZE DECREASE PNEWDT =ONE_R/TWO_R DDSDDE =CMATRIX DDSDDT =ZERO_R DRPLDE =ZERO_R DRPLDT =ZERO_R RPL =ZERO_R ! REWRITE INTO ABAQUS FORMAT IF (NTENS==3) THEN DO i=1,NTENS,1 DO j=3,NTENS,1 DDSDDE(i,j) =DDSDDE(i,j)*ONE_R/TWO_R END DO END DO ELSE DO i=1,NTENS,1 DO j=4,NTENS,1 DDSDDE(i,j) =DDSDDE(i,j)*ONE_R/TWO_R END DO END DO END IF STRESS =ZERO_R write(6,*) "Problem in Element:" write(6,*) NOEL,KINC write(6,*) "kinc" write(6,*) kinc END IF END IF !END NULLTEST RETURN END SUBROUTINE UMAT ! ! ! SUBROUTINE VECTORPRODUCT (NDI,NTENS,VECTOR1,VECTOR2,MATRIX) ! ! CALCULATES VECTOR1*(VECTOR2)T USE Abaqus_Interface IMPLICIT NONE ! ! RESULT MATRIX REAL(KIND=realkind),DIMENSION(NTENS,NTENS),INTENT(OUT) ::MATRIX ! INPUT VECTORS REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::VECTOR1,VECTOR2 ! INTEGER(KIND=intkind) ::i,j INTEGER(KIND=intkind),INTENT(IN) ::NTENS,NDI REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! ! MATRIX=ZERO_R ! DO i=1,NTENS,1 DO j=1,NTENS,1 MATRIX(i,j)=VECTOR1(i)*VECTOR2(j) END DO END DO ! RETURN ! END SUBROUTINE VECTORPRODUCT ! ! SUBROUTINE J3CALC (NTENS,VECTOR1,VECTOR2,VECTOR3,J3) ! ! CALCULATES THIRD INVARIANT OF STRESS DEVIATOR USE Abaqus_Interface IMPLICIT NONE REAL(KIND=realkind),INTENT(OUT) ::J3 ! INPUTVECTORS REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::VECTOR1,VECTOR2,VECTOR3 REAL(KIND=realkind),DIMENSION(NTENS) ::ZEVECTOR1,ZEVECTOR2 INTEGER(KIND=intkind) ::i,j INTEGER(KIND=intkind),INTENT(IN) ::NTENS REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! ! IF (NTENS==3) THEN J3 =ONE_R/THREE_R*(-THREE_R*VECTOR1(2)*VECTOR1(1)**2+THREE_R*VECTOR1(2)*VECTOR1(3)**2+THREE_R*VECTOR1(1)*(VECTOR1(3)**2-VECTOR1(2)**2)) ELSE CALL SCALARPRODUCT2 (NTENS,VECTOR1,VECTOR2,ZEVECTOR1) CALL SCALARPRODUCT2 (NTENS,ZEVECTOR1,VECTOR3,ZEVECTOR2) J3 =ONE_R/THREE_R*(ZEVECTOR2(1)+ZEVECTOR2(2)+ZEVECTOR2(3)) END IF RETURN ! END SUBROUTINE J3CALC ! SUBROUTINE SCALARPRODUCT2 (NTENS,VECTOR1,VECTOR2,RESULTVECTOR) ! ! SINGLE CONTRACTION OF SYMMETRIC SECOND ORDER TENSORS -> OUTPUT IS VECTOR-FORMAT USE Abaqus_Interface IMPLICIT NONE ! REAL(KIND=realkind),DIMENSION(NTENS),INTENT(OUT) ::RESULTVECTOR REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::VECTOR1,VECTOR2 REAL(KIND=realkind),DIMENSION(3,3) ::MATRIX1,MATRIX2 INTEGER(KIND=intkind),INTENT(IN) ::NTENS REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! INTRINSIC MATMUL ! MATRIX1 =ZERO_R MATRIX2 =ZERO_R RESULTVECTOR =ZERO_R ! MATRIX1(1,1) =VECTOR1(1) MATRIX1(2,2) =VECTOR1(2) MATRIX1(1,2) =VECTOR1(4) MATRIX1(2,1) =VECTOR1(4) MATRIX1(3,3) =VECTOR1(3) ! IF (NTENS==6) THEN MATRIX1(1,3) =VECTOR1(5) MATRIX1(3,1) =VECTOR1(5) MATRIX1(2,3) =VECTOR1(6) MATRIX1(3,2) =VECTOR1(6) END IF ! MATRIX2(1,1) =VECTOR2(1) MATRIX2(2,2) =VECTOR2(2) MATRIX2(1,2) =VECTOR2(4) MATRIX2(2,1) =VECTOR2(4) MATRIX2(3,3) =VECTOR2(3) ! IF (NTENS==6) THEN MATRIX2(1,3) =VECTOR2(5) MATRIX2(3,1) =VECTOR2(5) MATRIX2(2,3) =VECTOR2(6) MATRIX2(3,2) =VECTOR2(6) END IF ! ! CALCULATE MATRIX1 =MATMUL(MATRIX1,MATRIX2) ! ! REWRITE RESULTVECTOR(1) =MATRIX1(1,1) RESULTVECTOR(2) =MATRIX1(2,2) RESULTVECTOR(3) =MATRIX1(3,3) RESULTVECTOR(4) =MATRIX1(1,2) ! IF (NTENS==6) THEN RESULTVECTOR(5) =MATRIX1(1,3) RESULTVECTOR(6) =MATRIX1(2,3) END IF ! RETURN ! END SUBROUTINE SCALARPRODUCT2 ! SUBROUTINE STRAIN_HARDENING (PROPS, NPROPS,RRIN,ETA,LODE,R,HH,E_CI) ! ! CALCULATES STRAIN HARDENING AND DAMAGE INITIATION USE Abaqus_Interface IMPLICIT NONE ! ! REAL(KIND=realkind) ::E,NU,K,G REAL(KIND=realkind) ::TAU_0,H,A,RINF,N REAL(KIND=realkind) ::C_1,C_2,C_3,C_4,S_1,E_1 REAL(KIND=realkind) ::DDS,DAMAGE_CASE INTEGER(KIND=intkind),INTENT(IN) ::NPROPS REAL(KIND=realkind),DIMENSION(NPROPS),INTENT(IN) ::PROPS ! INPUT REAL(KIND=realkind),INTENT(IN) ::RRIN,ETA,LODE ! OUTPUT REAL(KIND=realkind),INTENT(OUT) ::R,HH,E_CI REAL(KIND=realkind) ::RR,ERROR,TOL,RES,ABL,THETA REAL(KIND=realkind),PARAMETER ::PI=3.141592653589793_realkind INTEGER(KIND=intkind) ::COUNTER,MAXIMUM REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! INTRINSIC EXP,LOG,ERF,SQRT,EPSILON,HUGE,ABS ! !!!!! MATERIAL PARAMETERS FROM PROPS ! ELASTICITY ! YOUNGS MODULUS E =PROPS(1) ! POISSONS RATIO NU =PROPS(2) ! STRAIN HARDENING ! INITIAL YIELD STRESS TAU_0 =PROPS(3) ! PREFACTOR HOLLOMON H =PROPS(4) ! EXPONENT HOLLOMON N =PROPS(5) ! SATURATION VOCE RINF =PROPS(6) ! PARAMETER VOCE A =PROPS(7) ! JOHNSON-COOK DAMAGE LOCUS C_1 =PROPS(8) C_2 =PROPS(9) C_3 =PROPS(10) C_4 =PROPS(11) ! DAMAGE PARAMETER S_1 =PROPS(12) ! E_C =E_1 IN NONLOCAL CASE E_1 =PROPS(13) DDS =PROPS(14) DAMAGE_CASE =PROPS(15) IF (RRIN=TOL .AND. PNEW==0) CALL DERIVATIVES (NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,0,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE) DET_K =KMATRIX L1 =-RMATRIX L1 =L1/DET_K IF (L1 .NE. L1) THEN L1 =DELTA_LAMBDA_PL*0.1 END IF DELTA_LAMBDA_PL =DELTA_LAMBDA_PL+L1 ERROR =sqrt(RMATRIX**2+L1**2) IF (COUNTER>=MAX_ITER) THEN PNEW =1 write(6,*) "MAXIMUM NUMBER OF NEWTON ATTEMPTS" END IF COUNTER=COUNTER+1 END DO ! END NEWTON ! STRESS UPDATE IF (PNEW==0) THEN CALL DERIVATIVES (NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,1,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE) STRESS =TAU_T-MATMUL(CMATRIX,(DELTA_LAMBDA_PL*N_TR)) ELSE STRESS =TAU_T END IF !!! MATERIAL TANGENT IF (PNEW==0) THEN DET_K =KMATRIX A1(1,1:NTENS) =-BMATRIX(1,1:NTENS) A1 =A1/DET_K DEVIATORTMATRIX(1:NTENS,1) =DEVIATORT(1:NTENS) N_TRMATRIX(1:NTENS,1) =N_TR(1:NTENS) PHIMATRIX =(IDEV-(THREE_R/TWO_R/TAU_MISES_T**2)*MATMUL(DEVIATORTMATRIX,MATMUL(TRANSPOSE(DEVIATORTMATRIX),P_SCALING_MATRIX)))*THREE_R/TWO_R/TAU_MISES_T DRPLDE =ZERO_R DRPLDE =MATMUL(CMATRIX,DRPLDE) STRESS_MATRIX(1:NTENS,1) =STRESS(1:NTENS) DDSDDE =IMATRIX-MATMUL(CMATRIX,((DELTA_LAMBDA_PL)*PHIMATRIX+MATMUL(N_TRMATRIX,A1))) DDSDDE2 =-RMATRIX*MATMUL(STRESS_MATRIX,A1) CALL DERIVATIVES (NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,2,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE) DET_K =KMATRIX L1 =ZERO_R DRPLDT =ZERO_R DDSDDT =ZERO_R ELSE DDSDDE =IMATRIX DDSDDE2 =ZERO_R DDSDDT =ZERO_R END IF RETURN ! END SUBROUTINE PLASTICITY_MODUL SUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX,PREDEF,DPRED) ! SUBROUTINE TO DEFINE THE HEAT SOURCE AND MATERIAL TANGENT ENTRY USE Abaqus_Interface IMPLICIT NONE ! ! CHARACTERS (MATERIALNAME) CHARACTER(80)::CMNAME REAL(KIND=realkind),DIMENSION(2),INTENT(IN) ::TEMP REAL(KIND=realkind),DIMENSION(2),INTENT(IN) ::TIME REAL(KIND=realkind),INTENT(IN) ::DTIME REAL(KIND=realkind),DIMENSION(2),INTENT(OUT) ::FLUX ! CHANGE THE DIMENSION OF STATEV TO THE CONSIDERED CASE REAL(KIND=realkind),dimension(12) ::STATEV REAL(KIND=realkind),dimension(:), allocatable ::PREDEF REAL(KIND=realkind),dimension(:), allocatable ::DPRED ! PLUG IN THE STATEV VARIABLE NUMBER OF YOUR STORED LOCAL DAMAGE PARAMETER FLUX(1) =-(TEMP(1)) FLUX(2) =-1.0_realkind+STATEV(11) RETURN END SUBROUTINE HETVAL SUBROUTINE PLASTICITY_MODUL_NL(DELTA_E_D_BAR,EPS_D_BAR,DELTA_E_EQ_WHOLE,NTENS,NPROPS,NSTATV,STARTDELTA, TAU_T,DEVIATORT,N_TR,CMATRIX,KRONECKER,IMATRIX,IDEV,P_SCALING_MATRIX,PROPS,STATEV,TAU_MISES_T,& & STRESS,DDSDDE,DDSDDT,DRPLDE,DRPLDT,DELTA_LAMBDA_PL,DELTA_RR,DAMAGE,PNEW) ! ! COMPUTES STRESS STATE (RADIAL RETURN) AND MATERIAL TANGENT ENTRIES (NONLOCAL VERSION) ! CRAMERS RULE TO SOLVE LINEAR SYSTEM OF EQUATIONS USE Abaqus_Interface IMPLICIT NONE ! ! INPUT REAL(KIND=realkind),DIMENSION(NTENS,NTENS),INTENT(IN) ::CMATRIX,IMATRIX,IDEV,P_SCALING_MATRIX REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::TAU_T,DEVIATORT,N_TR,KRONECKER INTEGER(KIND=intkind),INTENT(IN) ::NTENS,NPROPS,NSTATV REAL(KIND=realkind),DIMENSION(3),INTENT(IN) ::STARTDELTA REAL(KIND=realkind),INTENT(IN) ::TAU_MISES_T,EPS_D_BAR,DELTA_E_EQ_WHOLE,DELTA_E_D_BAR REAL(KIND=realkind),DIMENSION(NPROPS),INTENT(IN) ::PROPS REAL(KIND=realkind),DIMENSION(NSTATV),INTENT(IN) ::STATEV ! OUTPUT REAL(KIND=realkind),DIMENSION(NTENS),INTENT(OUT) ::STRESS,DDSDDT REAL(KIND=realkind),DIMENSION(NTENS,NTENS),INTENT(OUT) ::DDSDDE REAL(KIND=realkind),INTENT(OUT) ::DRPLDT REAL(KIND=realkind),DIMENSION(NTENS),INTENT(OUT) ::DRPLDE REAL(KIND=realkind),INTENT(OUT) ::DELTA_LAMBDA_PL,DAMAGE,DELTA_RR INTEGER(KIND=intkind),INTENT(INOUT) ::PNEW REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! REAL(KIND=realkind) ::E,NU,K,G REAL(KIND=realkind) ::TAU_0,H,A,RINF,N REAL(KIND=realkind) ::C_1,C_2,C_3,C_4,S_1,E_1 REAL(KIND=realkind) ::DDS,DAMAGE_CASE INTEGER(KIND=intkind) ::MAX_ITER,COUNTER REAL(KIND=realkind) ::TOL,ERROR,TAU_MISES REAL(KIND=realkind) ::KMATRIX REAL(KIND=realkind) ::RMATRIX REAL(KIND=realkind),DIMENSION(1,NTENS) ::BMATRIX REAL(KIND=realkind),DIMENSION(NTENS,1) ::STRESS_MATRIX REAL(KIND=realkind) ::DET_K, L1,EPSILON_MISES,RR,DAMAGE_SWITCH,DDAMAGE_D_DLAMBDA_PL REAL(KIND=realkind),DIMENSION(1,NTENS) ::A1 REAL(KIND=realkind),DIMENSION(NTENS,NTENS) ::PHIMATRIX REAL(KIND=realkind),DIMENSION(NTENS,1) ::DEVIATORTMATRIX,N_TRMATRIX INTRINSIC MATMUL,TRANSPOSE !!!! MATERIAL PARAMETERS FROM PROPS ! ELASTICITY ! YOUNGS MODULUS E =PROPS(1) ! POISSONS RATIO NU =PROPS(2) ! STRAIN HARDENING ! INITIAL YIELD STRESS TAU_0 =PROPS(3) ! PREFACTOR HOLLOMON H =PROPS(4) ! EXPONENT HOLLOMON N =PROPS(5) ! SATURATION VOCE RINF =PROPS(6) ! PARAMETER VOCE A =PROPS(7) ! JOHNSON-COOK DAMAGE LOCUS C_1 =PROPS(8) C_2 =PROPS(9) C_3 =PROPS(10) C_4 =PROPS(11) ! DAMAGE PARAMETER S_1 =PROPS(12) ! E_C =E_1 IN NONLOCAL CASE E_1 =PROPS(13) DDS =PROPS(14) DAMAGE_CASE =PROPS(15) K = E/(THREE_R*(ONE_R-TWO_R*NU)) G = E/(TWO_R*(ONE_R+NU)) EPSILON_MISES =STATEV(1) RR =STATEV(2) DAMAGE_SWITCH =STATEV(4) ! PREPARE LOCAL NEWTON ITERATION DELTA_LAMBDA_PL =STARTDELTA(1) IF (DELTA_LAMBDA_PL<=ZERO_R) THEN DELTA_LAMBDA_PL =DELTA_E_EQ_WHOLE*0.1 ELSE DELTA_LAMBDA_PL =STARTDELTA(1) END IF DELTA_RR =ZERO_R TOL =1E-10_REALKIND MAX_ITER =30 COUNTER =0 ERROR =TWO_R*TOL PNEW =0 !!! NEWTON-LOOP DO WHILE (ERROR>=TOL .AND. PNEW==0) CALL DERIVATIVES_NL (DELTA_E_D_BAR,EPS_D_BAR,NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,0,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE,DDAMAGE_D_DLAMBDA_PL) DET_K =KMATRIX L1 =-RMATRIX L1 =L1/DET_K IF (L1 .NE. L1) THEN L1 =DELTA_LAMBDA_PL*0.1 END IF DELTA_LAMBDA_PL =DELTA_LAMBDA_PL+L1 ERROR =sqrt(RMATRIX**2+L1**2) IF (COUNTER>=MAX_ITER) THEN PNEW =1 write(6,*) "MAXIMUM NUMBER OF NEWTON ATTEMPTS" END IF COUNTER=COUNTER+1 END DO ! END NEWTON ! STRESS UPDATE IF (PNEW==0) THEN CALL DERIVATIVES_NL (DELTA_E_D_BAR,EPS_D_BAR,NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,1,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE,DDAMAGE_D_DLAMBDA_PL) STRESS =TAU_T-MATMUL(CMATRIX,(DELTA_LAMBDA_PL*N_TR)) ELSE STRESS =TAU_T END IF !!! MATERIAL TANGENT IF (PNEW==0) THEN DET_K =KMATRIX A1(1,1:NTENS) =-BMATRIX(1,1:NTENS) A1 =A1/DET_K DEVIATORTMATRIX(1:NTENS,1) =DEVIATORT(1:NTENS) N_TRMATRIX(1:NTENS,1) =N_TR(1:NTENS) PHIMATRIX =(IDEV-(THREE_R/TWO_R/TAU_MISES_T**2)*MATMUL(DEVIATORTMATRIX,MATMUL(TRANSPOSE(DEVIATORTMATRIX),P_SCALING_MATRIX)))*THREE_R/TWO_R/TAU_MISES_T DRPLDE(1:NTENS) =A1(1,1:NTENS) STRESS_MATRIX(1:NTENS,1) =STRESS(1:NTENS) DDSDDE =IMATRIX-MATMUL(CMATRIX,((DELTA_LAMBDA_PL)*PHIMATRIX+MATMUL(N_TRMATRIX,A1))) CALL DERIVATIVES_NL (DELTA_E_D_BAR,EPS_D_BAR,NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,2,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE,DDAMAGE_D_DLAMBDA_PL) DET_K =KMATRIX L1 =-RMATRIX L1 =L1/DET_K IF (DAMAGE_SWITCH>0.5) THEN DRPLDT =L1 ELSE DRPLDT =ZERO_R DRPLDE =ZERO_R END IF DDSDDT =-MATMUL(CMATRIX,(N_TR*L1)) ELSE DDSDDE =IMATRIX DDSDDT =ZERO_R END IF RETURN ! END SUBROUTINE PLASTICITY_MODUL_NL SUBROUTINE DERIVATIVES (NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,FALL,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE) ! ! COMPUTES RESIDUALS AND DERIVATIVE OF RESIDUALS FOR LOCAL CASE USE Abaqus_Interface IMPLICIT NONE ! ! INPUT INTEGER(KIND=intkind),INTENT(IN) ::NTENS,NPROPS,NSTATV,FALL REAL(KIND=realkind),INTENT(IN) ::DELTA_LAMBDA_PL,TAU_MISES_T REAL(KIND=realkind),DIMENSION(NSTATV),INTENT(IN)::STATEV REAL(KIND=realkind),DIMENSION(NPROPS),INTENT(IN)::PROPS REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::DEVIATORT,KRONECKER REAL(KIND=realkind),DIMENSION(NTENS,NTENS),INTENT(IN)::P_SCALING_MATRIX,IDEV ! OUTPUT REAL(KIND=realkind),INTENT(OUT) ::KMATRIX REAL(KIND=realkind),INTENT(OUT) ::RMATRIX REAL(KIND=realkind),DIMENSION(1,NTENS),INTENT(OUT) ::BMATRIX REAL(KIND=realkind),INTENT(OUT) ::DELTA_RR,DAMAGE ! REAL(KIND=realkind),DIMENSION(NTENS,1) ::DEVIATORTMATRIX,DTAU_MISES_DTAU_T REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! REAL(KIND=realkind),PARAMETER ::PI=3.141592653589793_realkind REAL(KIND=realkind) ::E,NU,K,G REAL(KIND=realkind) ::TAU_0,H,A,RINF,N REAL(KIND=realkind) ::C_1,C_2,C_3,C_4,S_1,E_1 REAL(KIND=realkind) ::DDS,DAMAGE_CASE REAL(KIND=realkind) ::EPSILON_MISES,RR,DAMAGE_SWITCH REAL(KIND=realkind) ::R,HH REAL(KIND=realkind) ::TAU_MISES,DDAMAGE_D_DLAMBDA_PL,DRR_DDELTA_LAMBDA_PL,DR1_DTAU_MISES,DTAU_MISES_DLAMBDA_PL,DR1_DR,DR_DRR REAL(KIND=realkind) ::E_I,E_CI,A_D,B_D,C_D,E_C REAL(KIND=realkind),DIMENSION(NTENS,1) ::B1 REAL(KIND=realkind),DIMENSION(1,NTENS) ::B1TR INTRINSIC MATMUL,ATAN,EPSILON,EXP,LOG,HUGE,cosh,SINH ! !!! MATERIAL PARAMETERS FROM PROPS ! ELASTICITY ! YOUNGS MODULUS E =PROPS(1) ! POISSONS RATIO NU =PROPS(2) ! STRAIN HARDENING ! INITIAL YIELD STRESS TAU_0 =PROPS(3) ! PREFACTOR HOLLOMON H =PROPS(4) ! EXPONENT HOLLOMON N =PROPS(5) ! SATURATION VOCE RINF =PROPS(6) ! PARAMETER VOCE A =PROPS(7) ! JOHNSON-COOK DAMAGE LOCUS C_1 =PROPS(8) C_2 =PROPS(9) C_3 =PROPS(10) C_4 =PROPS(11) ! DAMAGE PARAMETER S_1 =PROPS(12) ! E_C =E_1 IN NONLOCAL CASE E_1 =PROPS(13) DDS =PROPS(14) DAMAGE_CASE =PROPS(15) K = E/(THREE_R*(ONE_R-TWO_R*NU)) G = E/(TWO_R*(ONE_R+NU)) EPSILON_MISES =STATEV(1) RR =STATEV(2) DAMAGE_SWITCH =STATEV(4) E_I =STATEV(5) E_C =E_I+E_1 ! COMPUTE DAMAGE IF (EPSILON_MISES=E_C) THEN DAMAGE =0.999 A_D =THREE_R/TANH(THREE_R)/(E_C-E_I) B_D =TANH(THREE_R) C_D =TWO_R*A_D*B_D*E_I+atanh(B_D) C_D =C_D/TWO_R/B_D DDAMAGE_D_DLAMBDA_PL =ZERO_R DRR_DDELTA_LAMBDA_PL =ONE_R ELSE A_D =THREE_R/TANH(THREE_R)/(E_C-E_I) B_D =TANH(THREE_R) C_D =TWO_R*A_D*B_D*E_I+atanh(B_D) C_D =C_D/TWO_R/B_D DAMAGE =0.999*ONE_R/TWO_R*(ONE_R+TANH(TWO_R*(A_D*B_D*(EPSILON_MISES+DELTA_LAMBDA_PL)-B_D*C_D))/B_D) DDAMAGE_D_DLAMBDA_PL =0.999*A_D*(ONE_R/COSH(TWO_R*(A_D*B_D*(EPSILON_MISES+DELTA_LAMBDA_PL)-B_D*C_D)))**2 DRR_DDELTA_LAMBDA_PL =ONE_R END IF TAU_MISES =TAU_MISES_T-THREE_R*G*DELTA_LAMBDA_PL DELTA_RR =DELTA_LAMBDA_PL CALL STRAIN_HARDENING (PROPS, NPROPS,RR+DELTA_RR,1.0,1.0,R,HH,E_CI) DR1_DTAU_MISES =ONE_R DTAU_MISES_DLAMBDA_PL =-THREE_R*G DR1_DR =-ONE_R*(ONE_R-DAMAGE) DR_DRR =HH KMATRIX =ZERO_R KMATRIX =DR1_DTAU_MISES*DTAU_MISES_DLAMBDA_PL+DR1_DR*DR_DRR*DRR_DDELTA_LAMBDA_PL+R*DDAMAGE_D_DLAMBDA_PL !!! SELECT, WHETHER NEWTON LOOP OR MATERIAL TANGENT SELECT CASE (FALL) CASE (0) ! RESIDUAL RMATRIX =ZERO_R RMATRIX =TAU_MISES-R*(ONE_R-DAMAGE) BMATRIX =ZERO_R CASE (1) ! MATERIAL TANGENT DEVIATORTMATRIX(1:NTENS,1) =DEVIATORT(1:NTENS) DTAU_MISES_DTAU_T(1:NTENS,1) =THREE_R/TWO_R/TAU_MISES_T*MATMUL(P_SCALING_MATRIX,DEVIATORT) BMATRIX =ZERO_R B1 =DR1_DTAU_MISES*DTAU_MISES_DTAU_T B1TR =TRANSPOSE(B1) BMATRIX(1,1:NTENS) =B1TR(1,1:NTENS) RMATRIX =DDAMAGE_D_DLAMBDA_PL CASE (2) RMATRIX =ZERO_R BMATRIX =ZERO_R END SELECT RETURN ! END SUBROUTINE DERIVATIVES SUBROUTINE DERIVATIVES_NL (DELTA_E_D_BAR,EPS_D_BAR,NTENS,NPROPS,NSTATV,DELTA_LAMBDA_PL,DEVIATORT,KRONECKER,P_SCALING_MATRIX,IDEV,PROPS,STATEV,TAU_MISES_T,FALL,& & KMATRIX,RMATRIX,BMATRIX,DELTA_RR,DAMAGE,DDAMAGE_D_DLAMBDA_PL) ! ! COMPUTES RESIDUALS AND DERIVATIVE OF RESIDUALS FOR NONLOCAL CASE USE Abaqus_Interface IMPLICIT NONE ! ! INPUT INTEGER(KIND=intkind),INTENT(IN) ::NTENS,NPROPS,NSTATV,FALL REAL(KIND=realkind),INTENT(IN) ::DELTA_LAMBDA_PL,TAU_MISES_T,EPS_D_BAR,DELTA_E_D_BAR REAL(KIND=realkind),DIMENSION(NSTATV),INTENT(IN)::STATEV REAL(KIND=realkind),DIMENSION(NPROPS),INTENT(IN)::PROPS REAL(KIND=realkind),DIMENSION(NTENS),INTENT(IN) ::DEVIATORT,KRONECKER REAL(KIND=realkind),DIMENSION(NTENS,NTENS),INTENT(IN)::P_SCALING_MATRIX,IDEV ! OUTPUT REAL(KIND=realkind),INTENT(OUT) ::KMATRIX REAL(KIND=realkind),INTENT(OUT) ::RMATRIX REAL(KIND=realkind),DIMENSION(1,NTENS),INTENT(OUT) ::BMATRIX REAL(KIND=realkind),INTENT(OUT) ::DELTA_RR,DAMAGE,DDAMAGE_D_DLAMBDA_PL ! REAL(KIND=realkind),DIMENSION(NTENS,1) ::DEVIATORTMATRIX,DTAU_MISES_DTAU_T REAL(KIND=realkind),PARAMETER ::ZERO_R=0.0_realkind,ONE_R=1.0_realkind,TWO_R=2.0_realkind,THREE_R=3.0_realkind,FOUR_R=4.0_realkind,FIVE_R=5.0_realkind,SIX_R=6.0_realkind,SEVEN_R=7.0_realkind,EIGHT_R=8.0_realkind,NINE_R=9.0_realkind,TEN_R=10.0_realkind ! REAL(KIND=realkind),PARAMETER ::PI=3.141592653589793_realkind REAL(KIND=realkind) ::E,NU,K,G REAL(KIND=realkind) ::TAU_0,H,A,RINF,N REAL(KIND=realkind) ::C_1,C_2,C_3,C_4,S_1,E_1 REAL(KIND=realkind) ::DDS,DAMAGE_CASE REAL(KIND=realkind) ::EPSILON_MISES,RR,DAMAGE_SWITCH REAL(KIND=realkind) ::R,HH REAL(KIND=realkind) ::TAU_MISES,DRR_DDELTA_LAMBDA_PL,DR1_DTAU_MISES,DTAU_MISES_DLAMBDA_PL,DR1_DR,DR_DRR REAL(KIND=realkind) ::E_I,E_CI,A_D,B_D,C_D,E_C REAL(KIND=realkind),DIMENSION(NTENS,1) ::B1 REAL(KIND=realkind),DIMENSION(1,NTENS) ::B1TR INTRINSIC MATMUL,ATAN,EPSILON,EXP,LOG,HUGE,cosh,SINH ! !!! MATERIAL PARAMETERS FROM PROPS ! ELASTICITY ! YOUNGS MODULUS E =PROPS(1) ! POISSONS RATIO NU =PROPS(2) ! STRAIN HARDENING ! INITIAL YIELD STRESS TAU_0 =PROPS(3) ! PREFACTOR HOLLOMON H =PROPS(4) ! EXPONENT HOLLOMON N =PROPS(5) ! SATURATION VOCE RINF =PROPS(6) ! PARAMETER VOCE A =PROPS(7) ! JOHNSON-COOK DAMAGE LOCUS C_1 =PROPS(8) C_2 =PROPS(9) C_3 =PROPS(10) C_4 =PROPS(11) ! DAMAGE PARAMETER S_1 =PROPS(12) ! E_C =E_1 IN NONLOCAL CASE E_1 =PROPS(13) DDS =PROPS(14) DAMAGE_CASE =PROPS(15) K = E/(THREE_R*(ONE_R-TWO_R*NU)) G = E/(TWO_R*(ONE_R+NU)) EPSILON_MISES =STATEV(1) RR =STATEV(2) DAMAGE_SWITCH =STATEV(4) E_I =ZERO_R E_C =E_I+E_1 ! COMPUTE DAMAGE IF (EPS_D_BAR>=E_C ) THEN DAMAGE =0.99_realkind A_D =THREE_R/TANH(THREE_R)/(E_C-E_I) B_D =TANH(THREE_R) C_D =TWO_R*A_D*B_D*E_I+atanh(B_D) C_D =C_D/TWO_R/B_D DDAMAGE_D_DLAMBDA_PL =ZERO_R DRR_DDELTA_LAMBDA_PL =ZERO_R ELSEIF (EPS_D_BAR<=ZERO_R .OR. EPS_D_BAR2.5) THEN DAMAGE =STATEV(3) DDAMAGE_D_DLAMBDA_PL =ZERO_R DRR_DDELTA_LAMBDA_PL =ZERO_R ELSE A_D =THREE_R/TANH(THREE_R)/(E_C-E_I) B_D =TANH(THREE_R) C_D =TWO_R*A_D*B_D*E_I+atanh(B_D) C_D =C_D/TWO_R/B_D DAMAGE =0.99_realkind*ONE_R/TWO_R*(ONE_R+TANH(TWO_R*(A_D*B_D*(EPS_D_BAR)-B_D*C_D))/B_D) DDAMAGE_D_DLAMBDA_PL =0.99_realkind*A_D*(ONE_R/COSH(TWO_R*(A_D*B_D*(EPS_D_BAR)-B_D*C_D)))**2 DRR_DDELTA_LAMBDA_PL =ZERO_R END IF TAU_MISES =TAU_MISES_T-THREE_R*G*DELTA_LAMBDA_PL DELTA_RR =DELTA_LAMBDA_PL CALL STRAIN_HARDENING (PROPS, NPROPS,RR+DELTA_RR,1.0,1.0,R,HH,E_CI) DR1_DTAU_MISES =ONE_R DTAU_MISES_DLAMBDA_PL =-THREE_R*G DR1_DR =-ONE_R*(ONE_R-DAMAGE) DR_DRR =HH KMATRIX =ZERO_R KMATRIX =DR1_DTAU_MISES*DTAU_MISES_DLAMBDA_PL+DR1_DR*DR_DRR !*(ONE_R-DAMAGE) !!! SELECT, WHETHER NEWTON LOOP OR MATERIAL TANGENT SELECT CASE (FALL) CASE (0) ! RESIDUAL RMATRIX =ZERO_R RMATRIX =TAU_MISES-R*(ONE_R-DAMAGE) BMATRIX =ZERO_R CASE (1) ! MATERIAL TANGENT DEVIATORTMATRIX(1:NTENS,1) =DEVIATORT(1:NTENS) DTAU_MISES_DTAU_T(1:NTENS,1) =THREE_R/TWO_R/TAU_MISES_T*MATMUL(P_SCALING_MATRIX,DEVIATORT) BMATRIX =ZERO_R B1 =DR1_DTAU_MISES*DTAU_MISES_DTAU_T B1TR =TRANSPOSE(B1) BMATRIX(1,1:NTENS) =B1TR(1,1:NTENS) RMATRIX =(ONE_R-DAMAGE) CASE (2) RMATRIX =R*DDAMAGE_D_DLAMBDA_PL BMATRIX =ZERO_R END SELECT RETURN ! END SUBROUTINE DERIVATIVES_NL