HILL48 +各向同性Voce硬化umat子程序

!=====================================================================

!     HILL48 PLASTICITY UMAT WITH VOCE ISOTROPIC HARDENING

!     IMPLEMENTATION: GENERALIZED RADIAL-RETURN IN EIGENSPACE

!===============================================================================

!

!  AUTHOR: 

!    Mohammad Hasaninia

!    Computational Advanced Manufacturing and Materials Laboratory (CAMML)

!    Department of Mechanical Engineering

!    University of Wyoming

!

!  REFERENCE:

!    Versino, D. and Bennett, K.C. (2018). Generalized radial-return mapping 

!    algorithm for anisotropic von Mises plasticity framed in material 

!    eigenspace. Int. J. Numer. Meth. Engng, 116: 202-222.

!    https://doi.org/10.1002/nme.5921

!

!-------------------------------------------------------------------------------

!  INPUT PROPERTIES (PROPS) DEFINITION:

!    PROPS(1)  : Young's Modulus (E)

!    PROPS(2)  : Poisson's Ratio (nu)

!    PROPS(3)  : Initial Yield Stress (Sigma_Y0)

!    PROPS(4)  : Saturation Stress (Q_inf) - Voce Parameter

!    PROPS(5)  : Hardening Rate (b)        - Voce Parameter

!    PROPS(6)  : Hill F

!    PROPS(7)  : Hill G

!    PROPS(8)  : Hill H

!    PROPS(9)  : Hill L

!    PROPS(10) : Hill M

!    PROPS(11) : Hill N

!

!  STATE VARIABLES (STATEV) OUTPUT:

!    STATEV(1)   : Equivalent Plastic Strain (alpha)

!    STATEV(2)   : Current Yield Stress

!    STATEV(3-8) : Plastic Strain Tensor (11, 22, 33, 12, 13, 23)

!    STATEV(9)   : Convergence Flag (0=Converged, 1=Failed)

!=====================================================================

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,       &

                KSTEP, KINC)

    IMPLICIT NONE

    ! Abaqus UMAT interface variables

    CHARACTER*80 CMNAME

    INTEGER NDI, NSHR, NTENS, NSTATV, NPROPS

    INTEGER NOEL, NPT, LAYER, KSPT, KSTEP, KINC

    REAL(8) STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS,NTENS)

    REAL(8) SSE, SPD, SCD, RPL, DRPLDT, DTIME, TEMP, DTEMP, PNEWDT, CELENT

    REAL(8) STRAN(NTENS), DSTRAN(NTENS), TIME(2), PREDEF(1), DPRED(1)

    REAL(8) PROPS(NPROPS), COORDS(3), DROT(3,3), DFGRD0(3,3), DFGRD1(3,3)

    REAL(8) DDSDDT(NTENS), DRPLDE(NTENS)

    ! Local variables for plasticity computation

    REAL(8) :: alpha_tr, alpha, StressYield, s_y, s_y_D, s_y_tr, f_gamma_tr

    REAL(8) :: Sigma_Y0, Q_inf, b_voce

    REAL(8) :: MatM(6,6), MatD(6,6), MatI(6,6)

    REAL(8) :: MatM_copy(6,6), Mat_inv(6,6)

    REAL(8) :: W(6), Q(6,6), Gamma(6,6), K(6,6), dK_dDelta_gamma_matrix(6,6)

    REAL(8) :: s_tr(6), strain_p(6), Stress_New(6), s_tilde_tr(6)

    REAL(8) :: Vec(6), dev_corrected(6)

    REAL(8) :: E, XNUE, EBULK3, EG2, EG, ELAM

    REAL(8) :: HillF, HillG, HillH, HillL, HillM, HillN

    REAL(8) :: pressure, d_ga, omega, d_al, f_gamma_delta, delta_alpha

    REAL(8) :: Df_gamma_d_d_gamma, d_omega_d_d_gamma, d_sy_d_d_gamma

    REAL(8) :: TOLERANCE, exp_term

    INTEGER :: iteration, K3, K5, K6, K7, K8, INFO

    LOGICAL :: converged

    ! Extract material properties

    E = PROPS(1)              ! Young's modulus

    XNUE = PROPS(2)           ! Poisson's ratio

    Sigma_Y0 = PROPS(3)       ! Initial yield stress

    Q_inf = PROPS(4)          ! Saturation stress (Q∞)

    b_voce = PROPS(5)         ! Voce hardening rate parameter

    HillF = PROPS(6)          ! Hill parameter F

    HillG = PROPS(7)          ! Hill parameter G

    HillH = PROPS(8)          ! Hill parameter H

    HillL = PROPS(9)          ! Hill parameter L

    HillM = PROPS(10)         ! Hill parameter M

    HillN = PROPS(11)         ! Hill parameter N

    ! Initialize state variables

    IF (TIME(2) < 1.0D-6) THEN

        alpha_tr = 0.0D0

        ! Voce: σ_y0 = σ_Y0 + Q∞(1 - e^0) = σ_Y0

        StressYield = Sigma_Y0

        strain_p = 0.0D0

        STATEV(1:9) = 0.0D0

    ELSE

        alpha_tr = STATEV(1)

        StressYield = STATEV(2)

        strain_p(1:6) = STATEV(3:8)

    END IF

    ! Create elasticity matrix

    EBULK3 = E/(1.0D0-2.0D0*XNUE)

    EG2 = E/(1.0D0+XNUE)

    EG = EG2/2.0D0

    ELAM = (EBULK3-EG2)/3.0D0

    MatD = 0.0D0

    MatD(1:3,1:3) = ELAM

    DO K3 = 1, 3

        MatD(K3,K3) = EG2 + ELAM

    END DO

    DO K3 = 4, 6

    MatD(K3, K3) = EG

    END DO

    ! Create Hill anisotropy matrix

    MatM = 0.0D0

    MatM(1,1) = HillG + HillH

    MatM(2,2) = HillF + HillH

    MatM(3,3) = HillF + HillG

    MatM(1,2) = -HillH

    MatM(1,3) = -HillG

    MatM(2,3) = -HillF

    MatM(2,1) = MatM(1,2)

    MatM(3,1) = MatM(1,3)

    MatM(3,2) = MatM(2,3)

    MatM(4,4) = 2.0D0*HillN

    MatM(5,5) = 2.0D0*HillM

    MatM(6,6) = 2.0D0*HillL

    ! Identity matrix

    MatI = 0.0D0

    DO K5 = 1, 6

        MatI(K5,K5) = 1.0D0

    END DO

    ! Compute trial stress (elastic predictor)

    Vec = MATMUL(MatD, DSTRAN)

    Stress_New = Vec + STRESS

    s_tr = Stress_New

    ! Extract deviatoric part

    pressure = (Stress_New(1)+Stress_New(2)+Stress_New(3))/3.0D0

    s_tr(1:3) = s_tr(1:3) - pressure

    ! Eigendecomposition of anisotropy tensor

    MatM_copy = MatM

    CALL compute_eigendecomposition(MatM_copy, W, Q, Gamma, INFO)

    IF (INFO /= 0) THEN

        WRITE(*,*) 'Error in eigendecomposition, INFO = ', INFO

        RETURN

    END IF

    ! Transform trial stress to eigenspace

    s_tilde_tr = MATMUL(TRANSPOSE(Q), s_tr)

    ! Initialize trial values

    alpha_tr = STATEV(1)

    s_y_tr = StressYield

    ! Compute trial yield function

    f_gamma_tr = 0.5D0 * DOT_PRODUCT(s_tilde_tr, MATMUL(Gamma, s_tilde_tr)) - s_y_tr**2

    ! Check yield condition

    IF (f_gamma_tr <= 0.0D0) THEN

        ! Elastic step

        STRESS = Stress_New

        d_ga = 0.0D0

        alpha = alpha_tr

        s_y = StressYield

    ELSE

        ! Plastic step - Newton-Raphson iteration

        tolerance = 5.0D-4

        d_ga = 0.0D0

        alpha = alpha_tr

        converged = .FALSE.

        DO iteration = 1, 250

            ! Build K matrix (diagonal)

            K = 0.0D0

            DO K6 = 1, 6

                K(K6,K6) = Gamma(K6,K6) / (1.0D0 + 2.0D0*EG*d_ga*Gamma(K6,K6))**2

            END DO

            ! Compute auxiliary variable omega

            omega = SQRT(0.5D0 * DOT_PRODUCT(s_tilde_tr, MATMUL(K, s_tilde_tr)))

            IF (omega <= 1.0D-14) THEN

                WRITE(*,*) 'Warning: omega is too s mall'

                EXIT

            END IF

            ! Update alpha

            delta_alpha = 2.0D0 * d_ga * omega

            alpha = alpha_tr + delta_alpha

            ! ========================================================

            ! VOCE HARDENING MODEL

            ! σ_y = σ_Y0 + Q∞(1 - exp(-b·α))

            ! ========================================================

            exp_term = EXP(-b_voce * alpha)

            s_y = Sigma_Y0 + Q_inf * (1.0D0 - exp_term)

            ! Derivative: ∂σ_y/∂α = Q∞·b·exp(-b·α)

            s_y_D = Q_inf * b_voce * exp_term

            ! ========================================================

            ! Compute residual function

            f_gamma_delta = (s_y / omega) - 1.0D0

            ! Check convergence

            IF (ABS(f_gamma_delta) <= tolerance) THEN

                StressYield = s_y

                converged = .TRUE.

                EXIT

            END IF

            ! Compute derivatives for Newton-Raphson

            dK_dDelta_gamma_matrix = 0.0D0

            DO K7 = 1, 6

                dK_dDelta_gamma_matrix(K7,K7) = -4.0D0*EG*Gamma(K7,K7)**2 &

                    / (1.0D0 + 2.0D0*EG*d_ga*Gamma(K7,K7))**3

            END DO

            d_omega_d_d_gamma = (1.0D0 / (4.0D0 * omega)) &

                * DOT_PRODUCT(s_tilde_tr, MATMUL(dK_dDelta_gamma_matrix, s_tilde_tr))

            ! Chain rule: ∂σ_y/∂Δγ = (∂σ_y/∂α)·(∂α/∂Δγ)

            d_sy_d_d_gamma = 2.0D0 * s_y_D * (omega + d_ga * d_omega_d_d_gamma)

            ! Derivative of residual

            Df_gamma_d_d_gamma = (1.0D0 / omega) * d_sy_d_d_gamma &

                - (s_y / (omega**2)) * d_omega_d_d_gamma

            ! Newton-Raphson update

            d_ga = d_ga - f_gamma_delta / Df_gamma_d_d_gamma

        END DO

        ! Update stress after Newton loop

        Mat_inv = 0.0D0

        DO K8 = 1, 6

            Mat_inv(K8,K8) = 1.0D0 / (1.0D0 + 2.0D0*EG*d_ga*Gamma(K8,K8))

        END DO

        dev_corrected = MATMUL(Q, MATMUL(Mat_inv, MATMUL(TRANSPOSE(Q), s_tr)))

        ! Assemble full Cauchy stress

        STRESS(1:3) = dev_corrected(1:3) + pressure

        STRESS(4:6) = dev_corrected(4:6)

        IF (.NOT. converged) THEN

            WRITE(*,*) 'WARNING: Newton-Raphson did not converge!'

            StressYield = s_y

            STATEV(9) = 1.0D0

        END IF

    END IF

    ! Consistent tangent operator

    DDSDDE = MatD

    ! Update state variables

    STATEV(1) = alpha                    ! Equivalent plastic strain

    STATEV(2) = s_y                      ! Current yield stress

    strain_p = strain_p + d_ga * MATMUL(MatM, dev_corrected)

    STATEV(3:8) = strain_p(1:6)         ! Plastic strain components

    RETURN

    CONTAINS

    !**********************************************************************

    SUBROUTINE compute_eigendecomposition(A, W, Q, Gamma, INFO)

        IMPLICIT NONE

        INTEGER, INTENT(OUT) :: INFO

        REAL(8), DIMENSION(6,6), INTENT(INOUT) :: A

        REAL(8), DIMENSION(6), INTENT(OUT) :: W

        REAL(8), DIMENSION(6,6), INTENT(OUT) :: Q

        REAL(8), DIMENSION(6,6), INTENT(OUT) :: Gamma

        INTEGER :: LWORK, i, j

        REAL(8), DIMENSION(:), ALLOCATABLE :: WORK

        LOGICAL :: is_symmetric

        ! Check symmetry

        is_symmetric = .TRUE.

        DO i = 1, 6

            DO j = i + 1, 6

                IF (ABS(A(i,j) - A(j,i)) > 1.0D-12) THEN

                    is_symmetric = .FALSE.

                    EXIT

                END IF

            END DO

            IF (.NOT. is_symmetric) EXIT

        END DO

        IF (.NOT. is_symmetric) THEN

            WRITE(*,*) "Error: Matrix not symmetric"

            INFO = -1

            RETURN

        END IF

        ! LAPACK eigendecomposition

        LWORK = -1

        ALLOCATE(WORK(1))

        CALL DSYEV('V', 'U', 6, A, 6, W, WORK, LWORK, INFO)

        LWORK = INT(WORK(1))

        DEALLOCATE(WORK)

        ALLOCATE(WORK(LWORK))

        CALL DSYEV('V', 'U', 6, A, 6, W, WORK, LWORK, INFO)

        IF (INFO == 0) THEN

            Q = A

            Gamma = 0.0D0

            DO i = 1, 6

                Gamma(i,i) = W(i)

            END DO

        ELSE

            WRITE(*,*) "Eigendecomposition failed, INFO:", INFO

            Gamma = 0.0D0

        END IF

        DEALLOCATE(WORK)

    END SUBROUTINE compute_eigendecomposition

END SUBROUTINE UMAT

!*****************************BOTTOM*****************************