SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) !PROGRAM UMAT IMPLICIT NONE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! CYCLIC CLAY PLASTICITY IMPLEMENTATION IN ABAQUS !!!! !!! FOR M.TECH PROJECT SUBROUTINE WRTITTEN BY !!!! !!! M.ASHOK KUMAR !!!! !!! UNDER THE GUIDANCE OF Dr. K.RAJAGOPAL !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAME ! INTEGER :: NTENS =6, NSTATV=3 , NPROPS =7 REAL(8) , PARAMETER:: TINY = 1.0E-18 REAL(8), PARAMETER:: PI=3.14159265358979323846264338327950288 ! INTEGER ,PARAMETER :: NTENS =6, NSTATV=3 , NPROPS =7 REAL(8):: STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) REAL(8) :: SSE,SPD,SCD,RPL,DRPLDT,DTIME,NSHR,PNEWDT, 1 CELENT,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,TEMP,DTEMP,E,LAMBDA, 2 KAPPA,P INTEGER , PARAMETER :: MAXNPROPS =50,MAXNASV =20, MAXN=30 INTEGER , PARAMETER :: MAXNY=30 REAL (8) :: DTMIN = 1.0E-17 , PERTURB = 1.0E-5, TOLINT =1.0E-3 INTEGER (8) :: MANINT=10000 INTEGER :: NPARMS, ERROR, NYACT,I,NASV,NFEV,NDI,NTENS,NSTATV, 1 NPROPS REAL(8) Y(MAXNY), Y_N(MAXNY), ASV(MAXNASV),SIG_EF REAL(8) DTSUB, PORE, DEPSV_NP1,NORM_DEPS2,NORM_DEPS,THETA,DEPSV INTEGER MAXNINT REAL(8) :: SIG_N(6),DEPS_NP1(6),DDTAN(6,6),DEPS(NTENS), 1 SIG_INITIAL(NTENS),P_INI,m,FREQ,N_CYC LOGICAL:: ELPRSW = .FALSE. REAL(8) :: PARMS(MAXNPROPS) ! INPUT PARAMTERS DECLARATION REAL(8) :: G, MYIELD, N, KPARAM=0,APARAM=0,SF=1 !!!!!!!!!!!!!!!!!!!!!!! EXECUTION SECTION!!!!!!!!!!!!!!!!!!!!!!!! IF (NPROPS>MAXNPROPS) THEN WRITE(*,*) 'Increase maxnprops in the subroutine' END IF PARMS = 0.D0 NPARMS = NPROPS PARMS(1:NPROPS) = PROPS(1:NPROPS) ! NO OF ADD STATE VARIABLES, INITIALISE THE PARAMETERS CALL INIT_PARMS (PARMS,NPARMS, NASV) IF (NSTATVMAXNASV) THEN WRITE (*,*) 'ERROR : INCREASE MAXNASVDIM IN UMAT' END IF NYACT =6+NASV IF (NYACT > MAXNY) THEN WRITE (*,*) 'ERROR: MAXNY TOO SMALL IN UMAT, PROGRAM TERMINATED' ERROR = 10 END IF Y = 0.D0 Y_N = 0.D0 ! ... Error Management: ! ---------------- ! error = 0 ... no problem in time integration ! error = 1 ... problems in evaluation of the time rate, (e.g. undefined ! stress state), reduce time integration substeps ! error = 3 ... problems in time integration, reduce abaqus load increment (cut-back) ! error = 10 ... severe error, terminate calculation ! ... check problem dimensions IF (NDI /=3) THEN WRITE (*,*) 'ERROR: THIS MODEL CAN BE USED ONLY FOR ELEMENTS' WRITE (*,*) ' WITH 3 DIRECT STRESS/STRAIN COMPONENTS' WRITE (*,*) 'NOEL = ', NOEL WRITE (*,*) 'NDI = ', NDI ERROR = 10 END IF DTSUB = STATEV(NASV+2) PORE = STATEV(NASV+1) ASV (1:NASV) = STATEV(1:NASV) SIG_N=0.D0 DEPS_NP1 = 0.D0 DEPSV_NP1=0.D0 CALL MOVE_SIG(STRESS,NTENS,PORE,SIG_N) CALL MOVE_EPS(DSTRAN,NTENS,DEPS_NP1,DEPSV_NP1) CALL INIY(Y,MAXNY,NASV,NTENS,SIG_N,ASV) Y_N = Y CALL DOT_VECT(2,DEPS_NP1,DEPS_NP1,NORM_DEPS2,NTENS) NORM_DEPS = SQRT(NORM_DEPS2) THETA = -PERTURB*MAX(NORM_DEPS,1.0E-6) MAXNINT = 10000 NFEV = 0 IF ((DTSUBDTIME)) THEN DTSUB = DTIME END IF IF (ERROR /=10) THEN CALL RKF23_UPDATE(Y,NYACT,NASV,DTSUB,TOLINT,MAXNINT,DTMIN, 1 TIME,DEPS_NP1, PARMS,NPARMS,NFEV,ELPRSW,DTIME,ERROR) END IF IF (ERROR == 3) THEN PNEWDT = 0.25 WRITE (*,*) 'SUBROUTINE UMAT: REDUCE STEP SIZE IN ABAQUS' ELSE IF (ERROR ==10) THEN WRITE(*,*) , 'ERROR = 10, UMAT EXITING....' CALL XIT END IF DDTAN = 0.D0 CALL PERTURBATE (Y_N, Y, NYACT, NASV, DTSUB, TOLINT, MAXNINT, 1 DTMIN, DEPS_NP1, PARMS, NPARMS, NFEV, ELPRSW, THETA, 2 NTENS, DDTAN, DTIME, ERROR,TIME) CALL SOLOUT (STRESS,NTENS,ASV,NASV,DDSDDE,Y,MAXNY,PORE, 1 DEPSV_NP1,PARMS,NPARMS,DDTAN,TIME,ERROR) DO I = 1,NASV STATEV(I) = ASV(I) END DO STATEV(NASV+1) = PORE STATEV(NASV+2) = DTSUB STATEV(NASV+3) = NFEV !! STATEV(2) = 1 !SENSTIVITY S = 1 TEMP IT WILL B DELETED RETURN END SUBROUTINE !END PROGRAM UMAT ! END PROGRAM UMAT !!!!!!!!!!!!!! INITIALIZE PARAMETERS~~~~~~~~~~~~~ SUBROUTINE INIT_PARMS(PARMS,NPARMS,NASV) IMPLICIT NONE INTEGER :: NPARMS, NASV,I REAL(8)::PARMS(NPARMS),G,MYIELD,KAPPA,LAMBDA,N, 1 KPARAM,APARAM,SF G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) DO I = 3,NPARMS IF (PARMS(I).LT.0) THEN WRITE(*,*) 'THE PARAMETER ', I ,'HAS NEGATIVE VALUE' CALL XIT END IF END DO NASV = 10 RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!TO FIND EFFECTIVE STRESS TENSOR!!! SUBROUTINE MOVE_SIG(STRESS, NTENS,PORE,SIG) IMPLICIT NONE INTEGER NTENS REAL(8) :: STRESS(NTENS),SIG(6) REAL(8) PORE !WRITE(8,*) 'ENTERING INTO MOVE_SIG' SIG =0.D0 SIG = STRESS SIG(1:3) = STRESS(1:3)+PORE RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! TO FIND VOL STRAIN & TO MOVE DSTRAN TO DEPS SUBROUTINE MOVE_EPS(DSTRAN,NTENS,DEPS,DEPSV) IMPLICIT NONE INTEGER :: NTENS REAL(8) :: DSTRAN(NTENS),DEPS(NTENS) REAL(8) :: DEPSV !WRITE(8,*) 'ENTERING INTO MOVE_EPS' DEPS = DSTRAN DEPSV = DEPS(1)+DEPS(2)+DEPS(3) RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!INITIALIZES THE VECTOR OF STATE VARIABLES !!----------------------------------------!! SUBROUTINE INIY(Y,MAXNY,NASV,NTENS,SIG,QQ2) IMPLICIT NONE !INTEGER , PARAMETER:: MAXNY = 30 INTEGER :: NASV,NTENS,MAXNY REAL(8) :: Y(MAXNY),SIG(NTENS),QQ2(NASV) !WRITE(8,*) 'ENTERING INTO INIY' Y(1:MAXNY) = 0.D0 Y(1:NTENS) = SIG(1:NTENS) Y(7:7+NASV) = QQ2 RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE DOT_VECT(FLAG,A,B,C,N) IMPLICIT NONE INTEGER :: FLAG,I,N !INTEGER, PARAMETER :: N REAL(8) :: A(N),B(N),C,COEFF,DOT_VECT1 !WRITE(8,*) 'ENTERING INTO DOT_VECT' !!! STRESS TENSOR !!!! DOT_VECT1=0.D0 IF (FLAG==1) THEN COEFF=2.D0 ELSE IF (FLAG==2) THEN !!! STAIN TENSOR!!!!! COEFF = 0.5 !!! STANDARD TENSOR!!! ELSE COEFF = 1.D0 END IF DO I=1,N IF (I<=3) THEN DOT_VECT1 = DOT_VECT1 + A(I)*B(I) ELSE DOT_VECT1 = DOT_VECT1 + COEFF * A(I)*B(I) END IF END DO C = DOT_VECT1 !WRITE(8,*) 'EXITING DOT_VECT' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE RKF23_UPDATE(Y,N,NASV,DTSUB,ERR_TOL,MAXNINT,DTMIN, 1 TIME,DEPS_NP1, PARMS, NPARMS, NFEV, ELPRSW, DTIME,ERROR) IMPLICIT NONE !! NUMERICAL SOULTION OF Y' = F(Y) !!EXPLICIT ADAPTIVE RKF23 SCHEME WITH LOCAL TIME STEP EXTRAPOLATION !!!INITIALIZE Y_K VECTOR & OTHER VARIABLES INTEGER , PARAMETER:: MAXN = 30 LOGICAL :: ELPRSW INTEGER :: N,NFEV,MAXNINT,ERROR,NASV,I,NPARMS,KSUBST REAL(8) :: Y(MAXN),Y_K(MAXN),Y_2(MAXN),Y_3(MAXN),Y_TIL(MAXN), 1 Y_HAT(MAXN),KRK_1(MAXN),KRK_2(MAXN),KRK_3(MAXN), 2 DEPS_NP1(6),DTMIN,DTSUB,DTIME,TIME(2) REAL(8) :: T_K, DT_K,NORM_R,S_HULL,PARMS(NPARMS),ERR_TOL IF (N > MAXN) THEN WRITE(*,*) 'INCREASE MAXN IN UMAT SUBROUTINE' RETURN END IF Y_K = 0.D0 Y_2 = 0.D0 Y_3 = 0.D0 Y_TIL = 0.D0 Y_HAT = 0.D0 KRK_1 = 0.D0 KRK_2 = 0.D0 KRK_3 = 0.D0 T_K = 0.D0 DT_K = 1.D0 KSUBST = 0.D0 NFEV = 0 Y_K = Y DO WHILE (T_K<1) KSUBST = KSUBST+1 IF (KSUBST>MAXNINT) THEN WRITE(*,*) 'NO OF SUBSTEPS',KSUBST WRITE(*,*) 'IS TOO BIG STEP IS REJECTED' ERROR = 3 RETURN END IF CALL RHS(Y_K,N,NASV,PARMS,NPARMS,DEPS_NP1,KRK_1,NFEV,ERROR,TIME) IF (ERROR == 10) THEN RETURN END IF !! FIND Y_2 DO I=1,N Y_2(I) = Y_K(I)+DT_K*KRK_1(I)/2 END DO CALL RHS(Y_2,N,NASV,PARMS,NPARMS,DEPS_NP1,KRK_2,NFEV,ERROR,TIME) IF (ERROR == 10) THEN RETURN END IF DO I=1,N Y_3(I) = Y_K(I)-DT_K*KRK_1(I)+2*DT_K*KRK_2(I) END DO CALL RHS(Y_3,N,NASV,PARMS,NPARMS,DEPS_NP1,KRK_3,NFEV,ERROR,TIME) IF (ERROR ==10) THEN RETURN END IF DO I=1,N Y_TIL(I)=Y_K(I)+DT_K*KRK_2(I) Y_HAT(I) = Y_K(I) + DT_K*(KRK_1(I)/6+2*KRK_2(I)/3+KRK_3(I)/6) END DO NORM_R=0.D0 CALL NORM_RES(Y_TIL,Y_HAT,N,NASV,NORM_R) S_HULL=1.D0 IF (NORM_R>0) THEN S_HULL =0.9*DT_K*((ERR_TOL/NORM_R)**(1.D0/3)) END IF IF (NORM_R0) THEN DO I = 1,6 ERR(I) = DEL_SIG(I)/NORM_SIG END DO IF (VOID_HAT /= 0) THEN ERR(7) = DEL_VOID/VOID_HAT END IF END IF CALL DOT_VECT(3,ERR,ERR,NORM_R2,8) NORM_R = SQRT(ABS(NORM_R2)) RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! MAIN SUBROUTINE !!!!!!!!!!!! SUBROUTINE GET_F_SIG_Q(SIG,Q,NASV,PARMS,NPARMS,DEPS,F_SIG,F_Q, 1 ERROR,TIME) IMPLICIT NONE INTEGER , PARAMETER :: MAXNASV =20 INTEGER :: NASV,NPARMS,ERROR,I,MAXITS REAL(8) :: ME(6,6), DEPS_E(6), DEPS_P(6),DFDSIG(6), DGDSIG(6), 1 F_YIELD, TMP , H, LAMBDA, SIG(6),PARMS(NPARMS), 2 TIME(2) ,SIG_T(6),D(6,6) REAL(8) :: DSIG_EL(6),F_SIG(6),Q(MAXNASV),DEPS(6),T,F_Q(NASV) LOGICAL :: PLASTI_LOAD,UNLOADING REAL(8) :: FTOL,SIG1(6),SIG2(6),SIG3(6),a,a0,a1,F_YIELD_NEW, 1 F_YIELD_SAVE, F_YIELD0,F_YIELD1 F_YIELD = 0.D0 TMP = 0.D0 H = 0.D0 LAMBDA = 0.D0 ME = 0.D0 DFDSIG = 0.D0 DGDSIG = 0.D0 DEPS_E = 0.D0 DEPS_P = 0.D0 CALL DM_ELASTIC (SIG,Q,ME,PARMS,NPARMS,UNLOADING,TIME,ERROR) DSIG_EL = 0.D0 DSIG_EL = MATMUL(ME,DEPS) CALL GET_FYIELD(SIG,Q,F_YIELD,TMP,PARMS,NPARMS,UNLOADING,TIME, 1 ERROR) CALL GET_PLASTI_TENS (PARMS,NPARMS,SIG,Q,DFDSIG,DGDSIG,H,TIME) CALL DOT_VECT(1, DSIG_EL, DFDSIG, T, 6) IF ((F_YIELD>=0) .AND. (T >0)) THEN PLASTI_LOAD = .TRUE. ELSE PLASTI_LOAD = .FALSE. END IF IF (PLASTI_LOAD == .FALSE.) THEN IF (T < 0 ) THEN UNLOADING = .TRUE. END IF ELSE !IF (PLASTI_LOAD) THEN CALL INTERSECT (SIG,Q,NASV,PARMS,NPARMS,DEPS,F_SIG,F_Q, 1 ERROR,TIME) CALL GET_FYIELD(SIG,Q,F_YIELD,TMP,PARMS,NPARMS,UNLOADING,TIME, 1 ERROR) !! CALCULATES PLASTIC MULTIPLIER LAMBDA CALL CALC_LAMBDA(DFDSIG, DGDSIG, H, ME, DEPS, LAMBDA) DEPS_E = DEPS-DEPS_P END IF CALL GET_FQ1 (PARMS,NPARMS,DEPS,SIG,Q,F_SIG, F_Q ,TIME,DEPS_P, 1 PLASTI_LOAD,UNLOADING) RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE GET_FQ (PARMS,NPARMS,DEPS, DEPS_E, SIG, Q, F_SIG, F_Q) IMPLICIT NONE INTEGER :: NPARMS REAL(8) :: E, DEPV, DEPS(6), DEPS_P(6), DEPV_P,SENS, EDEV_P(6), 1 F_Q(2), Q(20), DEPS_E(6),DEPSH_P,DDAMAGE,T,APARAM,KPARAM, 2 PARMS(NPARMS),KAPPA,LAMBDA,N, SIG(6),F_SIG(6) INTEGER :: SF !WRITE(8,*) 'ENTERED INTO GET_FQ' KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) E = Q(1) DEPV = DEPS(1) + DEPS(2) + DEPS(3) F_Q(1) = (1+E)*DEPV DEPS_P = DEPS - DEPS_E DEPV_P = DEPS_P(1) + DEPS_P(2) + DEPS_P(3) RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE INTERSECT (SIG,Q,NASV,PARMS,NPARMS,DEPS,F_SIG,F_Q, 1 ERROR,TIME) IMPLICIT NONE INTEGER , PARAMETER :: MAXNASV =20 INTEGER :: NASV,NPARMS,ERROR,I,MAXITS REAL(8) :: ME(6,6), DEPS_E(6), DEPS_P(6),DFDSIG(6), DGDSIG(6), 1 F_YIELD, TMP , H, LAMBDA, SIG(6),PARMS(NPARMS), 2 TIME(2) REAL(8) :: DSIG_EL(6),F_SIG(6),Q(MAXNASV),DEPS(6),T,F_Q(4),D(6,6) LOGICAL :: PLASTI_LOAD,UNLOADING REAL(8) :: FTOL,SIG1(6),SIG2(6),SIG3(6),a,a0,a1,F_YIELD_NEW, 1 F_YIELD_SAVE, F_YIELD0,F_YIELD1,SIG_TR(6),F_YIELD_TR MAXITS = 10 FTOL = 10E-6 CALL DM_ELASTIC (SIG,Q,ME,PARMS,NPARMS,UNLOADING,TIME,ERROR) CALL GET_FYIELD(SIG,Q,F_YIELD,TMP,PARMS,NPARMS,UNLOADING,TIME, 1 ERROR) DSIG_EL = 0.D0 SIG_TR = SIG + DSIG_EL CALL GET_FYIELD(SIG_TR,Q,F_YIELD_TR,TMP,PARMS,NPARMS,UNLOADING, 1 TIME,ERROR) IF (F_YIELD_TR <= FTOL) THEN F_SIG = DSIG_EL RETURN END IF IF ((F_YIELD < -FTOL) .AND. (F_YIELD_TR > FTOL)) CONTINUE !!! CHANGE FROM ELASTIC TO PLASTIC F_YIELD_SAVE = F_YIELD a0 = 0 a1 = 1 SIG1 = SIG + a0 * DSIG_EL CALL GET_FYIELD(SIG1,Q,F_YIELD0,TMP,PARMS,NPARMS,UNLOADING,TIME, 1 ERROR) SIG2 = SIG + a1 * DSIG_EL CALL GET_FYIELD(SIG2,Q,F_YIELD1,TMP,PARMS,NPARMS,UNLOADING,TIME, 1 ERROR) DO I = 1,MAXITS a = a1 - (a1 - a0)*(F_YIELD1/(F_YIELD1 - F_YIELD0)) ! WRITE(*,*) a, F_YIELD1,F_YIELD0 SIG3 = SIG + a * DSIG_EL CALL GET_FYIELD(SIG3,Q,F_YIELD_NEW,TMP,PARMS,NPARMS,UNLOADING, 1 TIME,ERROR) IF (ABS(F_YIELD_NEW) <= FTOL) THEN GO TO 9 END IF IF (((F_YIELD_NEW .LT. 0) .AND. (F_YIELD0 .GT. 0)) .OR. 1 ((F_YIELD_NEW .GT. 0) .AND. (F_YIELD0 .LT. 0))) THEN a1 = a F_YIELD1 = F_YIELD_NEW IF (((F_YIELD_NEW < 0) .AND. (F_YIELD_SAVE < 0)) .OR. 1 ((F_YIELD_NEW > 0) .AND. (F_YIELD_SAVE > 0))) THEN F_YIELD0 = F_YIELD0 /2 END IF ELSE a0 = a F_YIELD0 = F_YIELD_NEW IF (((F_YIELD_NEW < 0) .AND. (F_YIELD_SAVE < 0)) .OR. 1 ((F_YIELD_NEW > 0) .AND. (F_YIELD_SAVE > 0))) THEN F_YIELD1 = F_YIELD1 /2 END IF END IF F_YIELD_SAVE = F_YIELD_NEW END DO a = 0 9 CONTINUE F_SIG = a * DSIG_EL RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE GET_FQ1 (PARMS,NPARMS,DEPS,SIG,Q,F_SIG,F_Q,TIME,DEPS_P, 1 PLASTI_LOAD,UNLOADING) IMPLICIT NONE INTEGER :: NPARMS REAL(8) :: E, DEPV, DEPS(6), DEPS_P(6), DEPV_P,SENS, EDEV_P(6), 1 F_Q(7), Q(20), DEPS_E(6),DEPSH_P,DDAMAGE,T,APARAM,KPARAM, 2 PARMS(NPARMS),KAPPA,LAMBDA,N, SIG(6),F_SIG(6),TIME(2), 3 N_CYC,FREQ,m ,TETA,PC LOGICAL :: PLASTI_LOAD,UNLOADING KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) FREQ = PARMS(1) m = PARMS(9) E = Q(1) DEPV = DEPS(1) + DEPS(2) + DEPS(3) F_Q(1) = (1+E)*DEPV TETA = (1+E)/(LAMBDA-KAPPA) DEPV_P = DEPS_P(1) + DEPS_P(2) + DEPS_P(3) PC = Q(5) IF (PLASTI_LOAD) THEN F_Q(5) = -TETA*PC*DEPV_P END IF F_Q(6) = DEPV_P F_Q(7) = DEPV RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE GET_ELASTIC_MATRIX (SIG,Q,ME,PARMS,NPARMS,ERROR) IMPLICIT NONE INTEGER :: ERROR,NPARMS REAL(8) :: SIG(6),ME(6,6),P, KBULK,POISSON,G,YOUNG,E, 1 MYIELD, KAPPA,LAMBDA,N,PARMS(NPARMS),Q(20) P = -(SIG(1)+SIG(2)+SIG(3))/3 G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) !! G IS ASSUMED TO BE CONSTATNT & THE POISSON RATIO IS CHANGING WITH THE BULK MODULUS !KBULK = P/KAPPA E = Q(2) KBULK = P*(1+E)/KAPPA POISSON = (3*KBULK -2*G)/(2*G+6*KBULK) YOUNG = 2.D0 * G * (1.D0 + POISSON) CALL MAKE_MATRIX_ME(POISSON,YOUNG,ME) RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! FYIELD SUBROUTINE GET_FYIELD (SIG,Q,F_YIELD,G_FLOW,PARMS,NPARMS, 1 UNLOADING,TIME,ERROR) IMPLICIT NONE INTEGER :: NPARMS,ERROR REAL(8) :: GEO_SIGD(9),E,S,P,QDEV,P0,SIG(6),F_YIELD,G_FLOW,Q(20), 1 PARMS (NPARMS),G,MYIELD,KAPPA,LAMBDA,N,TIME(2),ETA,M,PC LOGICAL :: UNLOADING GEO_SIGD = 0.D0 G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) E = Q(1) PC = Q(5) P = -(SIG(1) + SIG(2)+ SIG(3))/3 IF (P < 0) THEN WRITE(*,*) 'ERROR :: P < 0 ' ERROR = 10 RETURN END IF CALL GET_QDEV(SIG,QDEV) !!! E CONSIDERED AS A STATE VARIABLE, PO IS CALCULATED !P0 = EXP((N-KAPPA*LOG(P)-(1+E))/(LAMBDA-KAPPA)) IF (TIME(1) == 0) PC = MAX(PC,P) F_YIELD = QDEV*QDEV-MYIELD*MYIELD*(P*(PC-P)) !WRITE(*,*) F_YIELD G_FLOW = F_YIELD RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE CHECK_UNLOAD (SIG,Q,F_YIELD,G_FLOW,PARMS,NPARMS, 1 UNLOADING,TIME,ERROR,DSIG_EL) IMPLICIT NONE INTEGER :: NPARMS,ERROR REAL(8) :: GEO_SIGD(9),E,S,P,QDEV,P0,SIG(6),F_YIELD,G_FLOW,Q(20), 1 PARMS (NPARMS) , G,MYIELD,KAPPA,LAMBDA,N,TIME(2), 2 DSIG_EL(6),Q2,PY,P_INC,DPY,D_Q LOGICAL :: UNLOADING !WRITE(8,*) 'ENTERED INTO GET_FYIELD' GEO_SIGD = 0.D0 G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) !CALL MAKE_FULLTENS(1,GEO_SIGD,SIG) !GEO_SIGD = -1*GEO_SIGD E = Q(1) P = -(SIG(1) + SIG(2)+ SIG(3))/3 IF (P < 0) THEN WRITE(*,*) 'ERROR :: P < 0 ' ERROR = 10 RETURN END IF CALL GET_QDEV(SIG,Q2) CALL GET_QDEV(DSIG_EL,D_Q) P0 = (FLOOR(1+TIME(1)*PARMS(1)))**(-PARMS(9)) * 1 EXP((N-KAPPA*LOG(P)-(1+E))/(LAMBDA-KAPPA)) P_INC = (DSIG_EL(1)+DSIG_EL(2)+DSIG_EL(3))/3 PY = P+Q2*Q2/(P*MYIELD*MYIELD) DPY = P_INC*(1-(Q2/(MYIELD*P))**2)+(2*Q2/(MYIELD*MYIELD*P))*D_Q IF (DPY < 0.D0) THEN UNLOADING = .TRUE. WRITE(*,*) TIME(1), 'UNLOADING' WRITE(*,*) DPY,PY,P,P_INC,D_Q,Q2 END IF F_YIELD = QDEV*QDEV-MYIELD*MYIELD*(P*(P0-P)) G_FLOW = F_YIELD RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE GET_PLASTI_TENS(PARMS,NPARMS,SIG,Q,DFDSIG,DGDSIG,H, 1 TIME) IMPLICIT NONE REAL(8) :: GEO_SIGD(9),E,S,ZERO(9),DG(9),DF(9),EPPD(9),DEV(9), 1 DFDEP(9), DGDEV(9), KRON_DELTA(3,3),SIG(6),Q(20), 2 G,P,P0,N,KAPPA,LAMBDA,MYIELD,DFDEPPACC,PARMS(NPARMS), 3 KPARAM,APARAM,SF,FDG,TMP1,H,DFDSIG(6),DGDSIG(6), 4 TMP2,TMP3,L,M,K,E0,TERM1,QDE,TIME(2),ETA,Q2 INTEGER :: I,J,NPARMS KRON_DELTA = 0.D0 KRON_DELTA(1,1) = 1.D0 KRON_DELTA(2,2) = 1.D0 KRON_DELTA(3,3) = 1.D0 G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) GEO_SIGD = 0.D0 CALL MAKE_FULLTENS(1,GEO_SIGD, SIG) GEO_SIGD = GEO_SIGD*(-1) E = Q(1) P0 = Q(5) S = 1 SF = 1 P = (GEO_SIGD(1) +GEO_SIGD(5)+ GEO_SIGD(9))/3 CALL GET_QDEV(SIG,Q2) !P0 = EXP((N-KAPPA*LOG(P)-(1+E))/(LAMBDA-KAPPA)) ETA = Q2/P M = MYIELD ZERO = 0.D0 DG = 0.D0 DF = 0.D0 EPPD = 0.D0 DEV = 0.D0 CALL MAKE_DEV (GEO_SIGD, DEV) !!! DFDSIG CALCULATED ANALYTICALLY DO I=0,2 DO J = 1,3 DF(3*I+J) = MYIELD*MYIELD*(2*P-P0)*(KRON_DELTA(I+1,J)/3)+ 1 3*DEV(3*I+J) DF(3*I+J) = -1*DF(3*I+J) END DO END DO DG = DF DFDEP = 0.D0 DGDEV = 0.D0 !DFDEP(1) = MYIELD*MYIELD*P*P0/(LAMBDA-KAPPA) !DFDEP(5) = MYIELD*MYIELD*P*P0/(LAMBDA-KAPPA) !DFDEP(9) = MYIELD*MYIELD*P*P0/(LAMBDA-KAPPA) TMP1 = DF(1)+ DF(5)+ DF(9) TMP2 = MYIELD*MYIELD*P*P0/(LAMBDA-KAPPA) H = -TMP1*TMP2 CALL MAKE_SHORTTENS(1,DF,DFDSIG) CALL MAKE_SHORTTENS(2,DF,DGDSIG) RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!! SUBROUTINE DMCLAY (SIG,Q,D,PARMS,NPARMS,TIME,ERROR) IMPLICIT NONE INTEGER :: NPARMS,J,JJ,ERROR REAL(8) :: GEO_SIGD(9),E,P,QDEV,P0,SIG(6),F_YIELD,G_FLOW,Q(20), 1 PARMS (NPARMS),G,MYIELD,KAPPA,LAMBDA,N, 2 SX,SY,SZ,TXY,TYZ,TZX,Q2,PC,PY,BK,AL,DL,PCS,D(6,6),S(6), 3 A(6),B(6),PB,BB,C,AA,AB,BETA,XI,P1,TIME(2),FREQ,m, 4 DP,DPY,DQ,Q1,SX1,SY1,SZ1,TXY1,TYZ1,TZX1,ETA GEO_SIGD = 0.D0 FREQ = PARMS(1) G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) SX = -SIG(1) SY = -SIG(2) SZ = -SIG(3) TXY = -SIG(4) TYZ = -SIG(5) TZX = -SIG(6) P = (SX+SY+SZ)/3. IF (P .LT. 0.)THEN WRITE(*,*) 'P LESS THAN ZERO USE INITIAL CONDITION' WRITE(*,*) 'STRESS ARRAY =', SIG ERROR = 10 GO TO 99 END IF Q2 = SX*(SX-SY)+SY*(SY-SZ)+SZ*(SZ-SX)+3.*TXY*TXY 1 +3.*TYZ*TYZ+3.*TZX+TZX Q2 = SQRT(Q2) E = Q(1) ETA = Q2/P M = MYIELD ! PC = (FLOOR(1+TIME(1)*PARMS(1)))**(-PARMS(9)) * ! 1 EXP((N-KAPPA*LOG(P)-(1+E))/(LAMBDA-KAPPA)) PC = PARMS(8)*(FLOOR(1+TIME(1)*PARMS(1)))**(-PARMS(9))* 1 (M**2 + ETA**2)/(M**2) PY = P+Q2*Q2/(P*MYIELD*MYIELD) BK = (1.+E)*P/KAPPA !! CALCULATE ELASTIC STRESS - STRAIN MATRIX !!! G = PARMS(3) IF (G .LT. 1.) THEN G = BK*1.5*(1.-2.*PARMS(3))/(1.+PARMS(3)) !WRITE(*,*) TIME(1) , 'G =',G END IF AL = (3.*BK+4.*G)/3. DL = (3.*BK-2.*G)/3. D = 0.D0 D(1,1) = AL D(2,1) = DL D(3,1) = DL D(1,2) = DL D(2,2) = AL D(3,2) = DL D(1,3) = DL D(2,3) = DL D(3,3) = AL D(4,4) = G D(5,5) = G D(6,6) = G IF (PY .LT. 0.9999*PC) GO TO 50 PCS = 0.5*PC PB = P/PCS S(1) = SX - P S(2) = SY - P S(3) = SZ - P S(4) = 2. * TXY S(5) = 2. * TYZ S(6) = 2. * TZX BB = -2.*(1.-PB)/(3.*PCS) C = 3./(PCS*PCS*MYIELD*MYIELD) A(1) = BB+C*S(1) A(2) = BB+C*S(2) A(3) = BB+C*S(3) A(4) = C*S(4) A(5) = C*S(5) A(6) = C*S(6) DO J =1,3 B(J) = 0.D0 DO JJ = 1,3 B(J) = B(J)+D(J,JJ)*A(JJ) END DO END DO B(4) = D(4,4)*A(4) B(5) = D(5,5)*A(5) B(6) = D(6,6)*A(6) XI = (LAMBDA-KAPPA)/(1.+E) AA = -4.*PB*(1.-PB)/(PCS*XI) AB =0.D0 DO J = 1,6 AB = AB+A(J)*B(J) END DO BETA = AA+AB DO J =1,6 DO JJ =1,6 D(JJ,J) = D(JJ,J) - B(JJ)*B(J)/BETA END DO END DO 50 CONTINUE !WRITE(*,*) 'DMCLAY CALL OVER..' 99 IF (P .LT. 0.)THEN WRITE(*,*) 'D MATRIX', D ERROR = 10 END IF RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE DM_ELASTIC (SIG,Q,D,PARMS,NPARMS,UNLOADING,TIME,ERROR) IMPLICIT NONE INTEGER :: NPARMS,J,JJ,ERROR REAL(8) :: GEO_SIGD(9),E,P,QDEV,P0,SIG(6),F_YIELD,G_FLOW,Q(20), 1 PARMS (NPARMS),G,MYIELD,KAPPA,LAMBDA,N, 2 SX,SY,SZ,TXY,TYZ,TZX,Q2,PC,PY,BK,AL,DL,PCS,D(6,6),S(6), 3 A(6),B(6),PB,BB,C,AA,AB,BETA,XI,P1,TIME(2),FREQ,m, 4 DP,DPY,DQ,Q1,SX1,SY1,SZ1,TXY1,TYZ1,TZX1,ETA LOGICAL :: UNLOADING !WRITE(*,*) ' STRESS IN DMCLAY = ', SIG !WRITE(*,*) ' STRESS= ', SIG GEO_SIGD = 0.D0 FREQ = PARMS(1) G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) SX = -SIG(1) SY = -SIG(2) SZ = -SIG(3) TXY = -SIG(4) TYZ = -SIG(5) TZX = -SIG(6) P = (SX+SY+SZ)/3. !IF ((P .LT. 0.) .AND. (TIME(1).NE. 0.D0)) THEN IF (P .LT. 0.)THEN WRITE(*,*) 'P LESS THAN ZERO USE INITIAL CONDITION' WRITE(*,*) 'STRESS ARRAY =', SIG ERROR = 10 END IF Q2 = SX*(SX-SY)+SY*(SY-SZ)+SZ*(SZ-SX)+3.*TXY*TXY 1 +3.*TYZ*TYZ+3.*TZX+TZX Q2 = SQRT(Q2) E = Q(1) PC = Q(5) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PY = P+Q2*Q2/(P*MYIELD*MYIELD) BK = (1.+E)*P/KAPPA G = PARMS(3) IF (G .LT. 1.) THEN G = BK*1.5*(1.-2.*PARMS(3))/(1.+PARMS(3)) END IF AL = (3.*BK+4.*G)/3. DL = (3.*BK-2.*G)/3. D = 0.D0 D(1,1) = AL D(2,1) = DL D(3,1) = DL D(1,2) = DL D(2,2) = AL D(3,2) = DL D(1,3) = DL D(2,3) = DL D(3,3) = AL D(4,4) = G D(5,5) = G D(6,6) = G RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE DM_ELASTIC_UNLOAD (SIG,Q,D,PARMS,NPARMS,UNLOADING,TIME, 1 ERROR) IMPLICIT NONE INTEGER :: NPARMS,J,JJ,ERROR REAL(8) :: GEO_SIGD(9),E,P,QDEV,P0,SIG(6),F_YIELD,G_FLOW,Q(20), 1 PARMS (NPARMS),G,MYIELD,KAPPA,LAMBDA,N, 2 SX,SY,SZ,TXY,TYZ,TZX,Q2,PC,PY,BK,AL,DL,PCS,D(6,6),S(6), 3 A(6),B(6),PB,BB,C,AA,AB,BETA,XI,P1,TIME(2),FREQ,M, 4 DP,DPY,DQ,Q1,SX1,SY1,SZ1,TXY1,TYZ1,TZX1,ETA LOGICAL :: UNLOADING UNLOADING = .TRUE. GEO_SIGD = 0.D0 FREQ = PARMS(1) G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) SX = -SIG(1) SY = -SIG(2) SZ = -SIG(3) TXY = -SIG(4) TYZ = -SIG(5) TZX = -SIG(6) P = (SX+SY+SZ)/3. IF (P .LT. 0.)THEN WRITE(*,*) 'P LESS THAN ZERO USE INITIAL CONDITION' WRITE(*,*) 'STRESS ARRAY =', SIG ERROR = 10 END IF BK = (1.+E)*P/KAPPA !! CALCULATE ELASTIC STRESS - STRAIN MATRIX !!! G = PARMS(3) IF (G .LT. 1.) THEN G = BK*1.5*(1.-2.*PARMS(3))/(1.+PARMS(3)) END IF AL = (3.*BK+4.*G)/3. DL = (3.*BK-2.*G)/3. D = 0.D0 D(1,1) = AL D(2,1) = DL D(3,1) = DL D(1,2) = DL D(2,2) = AL D(3,2) = DL D(1,3) = DL D(2,3) = DL D(3,3) = AL D(4,4) = G D(5,5) = G D(6,6) = G RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE PERTURBATE (Y_N, Y_NP1, N, NASV, DTSUB, ERR_TOL, 1 MAXNINT,DTMIN, DEPS_NP1, PARMS, NPARMS, NFEV, ELPRSW, THETA, 2 NTENS, DD, DTIME, ERROR,TIME) !COMPUTE NUMERICALLY CONSISTENT TANGET STIFFNESS IMPLICIT NONE INTEGER , PARAMETER :: MAXN =30 REAL(8) :: Y_N(MAXN),Y_NP1(MAXN),Y_STAR(MAXN),DEPS_STAR(6), 1 DSIG(6),DD(NTENS,NTENS),DEPS_NP1(NTENS),THETA, 2 DTSUB, ERR_TOL,DTMIN,PARMS(NPARMS),DTIME,TIME(2) INTEGER :: JJ,KK ,ERROR,N,NASV,MAXNINT,NFEV,NTENS,NPARMS LOGICAL :: ELPRSW !WRITE(8,*) 'ENTERED INTO PERTURBATE' Y_STAR = 0.D0 DEPS_STAR =0.D0 DSIG = 0.D0 DD = 0.D0 !! LOOP OVER STRAIN BASIS VECTORS DO JJ = 1,NTENS Y_STAR = Y_N !! PERTURBED STRAIN INCREMENT DEPS_STAR = DEPS_NP1 DEPS_STAR(JJ) = DEPS_STAR(JJ)+THETA !! PERTURBED FINAL STATE, STORED IN Y_STAR IF (ERROR /= 10) THEN CALL RKF23_UPDATE (Y_STAR,N,NASV,DTSUB,ERR_TOL,MAXNINT, 1 DTMIN,TIME,DEPS_STAR,PARMS,NPARMS,NFEV,ELPRSW, 2 DTIME,ERROR) END IF !! STIFFNESS COMPONENTS IN COLUMN JJ !! (KK=ROW INDEX, JJ =COLUMN INDEX) DO KK = 1,NTENS DSIG(KK) = Y_STAR(KK) -Y_NP1(KK) DD(KK,JJ) = DSIG(KK)/THETA END DO END DO !WRITE(8,*) 'EXITING PERTURBATE' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE SOLOUT (STRESS,NTENS,ASV,NASV,DDSDDE,Y,MAXNY,PORE, 1 DEPSV_NP1,PARMS,NPARMS,DD,TIME,ERROR) IMPLICIT NONE REAL(8) :: BULK_W, STRESS(NTENS),ASV(NASV),DDSDDE(6,6),PORE, 1 PARMS(NPARMS),DEPSV_NP1,Y(MAXNY),DD(6,6),TIME(2),ETA,Q2, 2 E,P0,m,FREQ,N_CYC,PC0,P,G,N,MYIELD,LAMBDA,KAPPA,A(NTENS), 3 DEG_TMP,NCYC,T,RAT INTEGER :: I,J,NTENS,NASV,MAXNY,NPARMS,Q,ERROR MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) FREQ = PARMS(1) m = PARMS(9) !! UPDATE EXCESS PORE PRESSURE (IF UNDRAINED CONDITIONS), COMPRESSION POSITIVE BULK_W = PARMS(2) !WRITE(*,*) 'PORE BEFORE =',PORE PORE = PORE - BULK_W*DEPSV_NP1 !WRITE(*,*) 'PORE AFTER =',PORE !!UPDATED TOTAL STRESSES (EFF STRESSES STORED IN Y(1:6) DO I=1,NTENS IF (I<=3) THEN STRESS(I) = Y(I) - PORE ELSE STRESS(I) = Y(I) END IF END DO !! ADDITIONAL STATE VARIABLES (FIRST 6 COMPONENTS ARE INTERGRANULAR STRAINS) DO I=1,NASV ASV(I) = Y(6+I) END DO IF (TIME(1) /= 0) THEN ASV(2) = -(STRESS(1)+STRESS(2)+STRESS(3))/3 ASV(3) = -STRESS(1)+STRESS(2) ASV(4) = ASV(2) - PORE END IF NCYC = FLOOR(1+TIME(1)*PARMS(1)) T = 1./PARMS(1) !DEG_TMP = ASV(5)*((FLOOR(1+TIME(1)*FREQ))**-m) ETA = ASV(3)/ASV(4) RAT = (MYIELD**2+ETA**2)/MYIELD**2 DEG_TMP = PARMS(8)*((FLOOR(1+TIME(1)*FREQ))**-m)*RAT IF ((TIME(1) > (NCYC*T - 0.75*T)) .AND. ASV(9) < NCYC) THEN ASV(5) = DEG_TMP ASV(9) = NCYC END IF ASV(8) = (FLOOR(1+TIME(1)*FREQ)) DDSDDE = DD DO I = 1,3 DO J = 1,3 DDSDDE (I,J) = DDSDDE(I,J) + BULK_W END DO END DO write(*,*) 'solout ending' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE PLASTIC_STR (F_SIG,Q,SIG,DEPS_P,PARMS,NPARMS) IMPLICIT NONE INTEGER :: NPARMS,I REAL(8) :: DSIG(6),DEPS_P(6),DP,DQ,FULL_SIG(9),PARMS(NPARMS), 1 GEO_SIGD(9),G,MYIELD,KAPPA,LAMBDA,N,E0,E,S, 2 SIG(9),P,QDEV,P0,TERM0,TERM1,TERM2,TERM3, 3 KRONDELTA(6),F_SIG(6),Q(20) DP = 0.D0 DQ = 0.D0 CALL MAKE_FULLTENS(1,F_SIG,FULL_SIG) !FULL_SIG = -1*FULL_DSIG !WRITE (8,*) 'FULL_SIG =' ,FULL_SIG DP = (FULL_SIG(1) + FULL_SIG(5)+ FULL_SIG(9))/3 !WRITE (8,*) 'DP =' ,DP CALL GET_QDEV(F_SIG,DQ) !WRITE (8,*) 'DQ =' ,DQ GEO_SIGD = 0.D0 G = PARMS(3) MYIELD = PARMS(4) KAPPA = PARMS(5) LAMBDA = PARMS(6) N = PARMS(7) E0 = PARMS(12) E = Q(1) S = Q(2) CALL MAKE_FULLTENS(1,GEO_SIGD,SIG) GEO_SIGD = -1*GEO_SIGD P = (GEO_SIGD(1) + GEO_SIGD(5)+ GEO_SIGD(9))/3 CALL GET_QDEV(SIG,QDEV) !!! E CONSIDERED AS A STATE VARIABLE, PO IS CALCULATED !P0 = EXP((N-KAPPA*LOG(P)-LOG(1+E))/(LAMBDA-KAPPA)) P0 = PARMS(11) P0 = 206.7*(206.7/P)**(KAPPA/(LAMBDA-KAPPA)) !F_YIELD = QDEV*QDEV-MYIELD*MYIELD*(P*(P0*S-P)) TERM1 = (MYIELD*MYIELD*P*P-QDEV*QDEV)/(MYIELD*MYIELD*P*P 1 +QDEV*QDEV) TERM1 = TERM1/(3*P) TERM2 = 3.D0/((MYIELD*MYIELD*P*P+QDEV*QDEV)) TERM3 = DP + ((2*P*QDEV)*DQ/(MYIELD*MYIELD*P*P-QDEV*QDEV)) TERM0 = (LAMBDA-KAPPA)/(1+E0) KRONDELTA = 0.D0 KRONDELTA (1) = 1.D0 KRONDELTA (2) = 1.D0 KRONDELTA (3) = 1.D0 DO I = 1,6 DEPS_P (I) = TERM0*(TERM1*KRONDELTA(I)+TERM2*(SIG(I) 1 -P*KRONDELTA(I)))*TERM3 END DO !WRITE (8,*) 'DEPS_P =' ,DEPS_P RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE CALC_LAMBDA (DFDSIG, DGDSIG, H, ME, DEPS, LAMBDA) IMPLICIT NONE REAL(8) :: NDOTME(6),CHI,ME(6,6), 1 DFDSIG(6),DGDSIG(6) ,H, LAMBDA_MLP(6),DEPS(6),LAMBDA !INTEGER :: I !WRITE(8,*) 'ENTERED INTO CALC_LAMBDA' NDOTME = 0.D0 !CALL MAT_MUL(ME,DFDSIG,NDOTME,6,6,1) !! ME IS SYMMETRICA NDOTME = MATMUL(ME,DFDSIG) CALL DOT_VECT(2,NDOTME,DGDSIG,CHI,6) CHI = CHI+H IF (CHI == 0) THEN WRITE (*,*) 'PROBLEM IN CALCULATION OF SCALAR MULTIPLIER' STOP END IF CALL ARRAY_MULTIPLY (NDOTME, LAMBDA_MLP,1.D0/CHI,6) CALL DOT_VECT(2,LAMBDA_MLP,DEPS,LAMBDA,6) RETURN END SUBROUTINE !!!!!!!!!ARRAY_MULTIPLY!!!!!! SUBROUTINE ARRAY_MULTIPLY(A,B,C,N) IMPLICIT NONE REAL(8) :: B(N),A(N),C INTEGER :: I,N !WRITE(8,*) 'ENTERED INTO ARRAY_MULTIPLY' DO I = 1,N B(I) = C*A(I) END DO RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE ARRAY_SUBTRACT (A,B,C,N) IMPLICIT NONE REAL(8) :: C(N),A(N),B(N) INTEGER :: I,N !WRITE(8,*) 'ENTERED INTO ARRAY_SUBTRACT' DO I=1,N C(I) = A(I) - B(I) END DO RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE GET_QDEV(SIG,QDEV) IMPLICIT NONE REAL(8) :: QDEV,SIG(6),A1,A2,A3,A4,A5,A6 !WRITE(8,*) 'ENTERED INTO GET_QDEV' A1 = (SIG(1)-SIG(2))**2 A2 = (SIG(2)-SIG(3))**2 A3 = (SIG(1)-SIG(3))**2 A4 = (SIG(4))**2 A5 = (SIG(5))**2 A6 = (SIG(6))**2 QDEV = SQRT(0.5*(A1+A2+A3)+ 3.D0*(A4+A5+A6)) !WRITE(8,*) 'EXITING GET_QDEV' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE MAKE_MATRIX_ME(POISSON,YOUNG,ME) IMPLICIT NONE REAL(8) :: FAC, ME(6,6),POISSON,YOUNG !WRITE(8,*) 'ENTERED INTO MAKE_MATRIX_ME' !WRITE(*,*) 'EVALUATING ELSATIC MATRIX...' ME = 0.D0 FAC = 0.D0 FAC =(YOUNG*(1.D0 -POISSON))/((1.+POISSON)*(1.D0-2.D0*POISSON)) ME(1,1) = FAC ME(2,2) = FAC ME(3,3) = FAC ME(1,2) = FAC*(POISSON/(1.D0-POISSON)) ME(1,3) = FAC*(POISSON/(1.D0-POISSON)) ME(2,3) = FAC*(POISSON/(1.D0-POISSON)) ME(2,1) = FAC*(POISSON/(1.D0-POISSON)) ME(3,1) = FAC*(POISSON/(1.D0-POISSON)) ME(3,2) = FAC*(POISSON/(1.D0-POISSON)) ME(4,4) = FAC*((1.D0-2.D0*POISSON)/(1.D0-POISSON)) ME(5,5) = FAC*((1.D0-2.D0*POISSON)/(1.D0-POISSON)) ME(6,6) = FAC*((1.D0-2.D0*POISSON)/(1.D0-POISSON)) ! WRITE(*,*) 'CHECK @ MAKE_MATRIX_ME, ME =', ME !WRITE(8,*) 'EXITING MAKE_MATRIX_ME' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE MAKE_FULLTENS (FLAG,FULLTENS,SHORTTENS) IMPLICIT NONE REAL(8) :: COEFF , FULLTENS(9),SHORTTENS(6) INTEGER :: FLAG ! WRITE(*,*) , 'USING MAKE_FULLTENS...' !!!STRESS TENSOR COEFF = 1.D0 !!STRAIN TENSOR IF (FLAG == 2) THEN COEFF = 0.5 END IF FULLTENS = 0.D0 FULLTENS(1) = SHORTTENS(1) FULLTENS(5) = SHORTTENS(2) FULLTENS(9) = SHORTTENS(3) FULLTENS(2) = SHORTTENS(4)*COEFF FULLTENS(4) = SHORTTENS(4)*COEFF FULLTENS(3) = SHORTTENS(5)*COEFF FULLTENS(7) = SHORTTENS(5)*COEFF FULLTENS(6) = SHORTTENS(6)*COEFF FULLTENS(8) = SHORTTENS(6)*COEFF !WRITE(*,*) , 'MAKE_FULLTENS OVER...', FULLTENS !WRITE(*,*) 'CHECK @ GET_PLASTI_TENS, FULLTENS',FULLTENS RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 SUBROUTINE MAKE_DEV (TNZ,DEV) IMPLICIT NONE REAL(8) :: TNZM, TNZ(9), DEV(9) INTEGER :: IDIM DEV = 0.D0 TNZM = (TNZ(1)+TNZ(5)+TNZ(9))/3.D0 DEV = TNZ DO IDIM =0,2 DEV(IDIM*3+IDIM+1) = DEV(IDIM*3+IDIM+1) -TNZM END DO RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE ARRAY_INPRODUCT (A,B,RESULT1,N) IMPLICIT NONE REAL(8) :: RESULT1,A(N),B(N) INTEGER :: I,N !WRITE(8,*) 'ENTERED INTO ARRAY_INPRODUCT' RESULT1 = 0.D0 DO I=1,N RESULT1 = RESULT1 + A(I)*B(I) END DO !WRITE(8,*) 'EXITING ARRAY_INPRODUCT' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE MAKE_SHORTTENS (FLAG,FULLTENS,SHORTTENS) IMPLICIT NONE REAL(8):: COEFF =1.D0, SHORTTENS(6), FULLTENS(9) INTEGER :: FLAG !WRITE(8,*) 'ENTERED INTO MAKE_SHORTENS' IF (FLAG==2) THEN COEFF = 0.5 END IF SHORTTENS = 0.D0 SHORTTENS(1) = FULLTENS (1) SHORTTENS(2) = FULLTENS (5) SHORTTENS(3) = FULLTENS (9) SHORTTENS(4) = FULLTENS (2)/COEFF SHORTTENS(5) = FULLTENS (3)/COEFF SHORTTENS(6) = FULLTENS (6)/COEFF !WRITE(8,*) 'EXITING MAKE_SHORTENS' RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE MAKE_DEV_SHRT (TNZ, DEV) IMPLICIT NONE REAL (8) :: DEV(6),TNZ(6),TNZM INTEGER :: I !WRITE(8,*) 'ENTERED INTO MAKE_DEV_SHRT' DEV = 0.D0 TNZM = (TNZ(1) + TNZ(2) + TNZ(3))/3.D0 DEV = TNZ DO I = 1,3 DEV(I) = DEV(I) - TNZM END DO !WRITE(8,*) 'EXITING MAKE_DEV_SHRT' RETURN END SUBROUTINE ! SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT, ! 1 TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER, ! 2 KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO, ! 3 LACCFLA) ! ! INCLUDE 'ABA_PARAM.INC' ! CHARACTER*80 CMNAME,ORNAME ! CHARACTER*3 FLGRAY(15) ! DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3), ! 1 T(3,3),TIME(2) ! DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*), ! 1 COORD(*) ! CALL GETVRM('POR',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP, ! 1 MATLAYO,LACCFLA) ! PORE = ARRAY(1) ! STATEV(11) = PORE ! RETURN ! END SUBROUTINE ! SUBROUTINE SIGINI(SIGMA,COORDS,NTENS,NCRDS,NOEL,NPT,LAYER, ! 1 KSPT,LREBAR,NAMES) ! INCLUDE 'ABA_PARAM.INC' ! DIMENSION SIGMA(NTENS),COORDS(NCRDS) ! CHARACTER NAMES(2)*80 ! SIGMA(1) = -200 ! SIGMA(2) = -200 ! SIGMA(3) = -200 ! SIGMA(4) = 0.D0 ! SIGMA(5) = 0.D0 ! SIGMA(6) = 0.D0 ! RETURN ! END