From f8cf61706c8f12edac0bf0e23205c6f7d6cf8953 Mon Sep 17 00:00:00 2001
From: Mohib <mohib.mustafa@gmail.com>
Date: Fri, 19 Apr 2024 09:32:42 +0200
Subject: [PATCH] [Refractor] - Numerical Tangent as an exteral function call
---
.../FiniteStrain/Finite_VEVP/umat.f | 220 ++++++++++++++----
1 file changed, 172 insertions(+), 48 deletions(-)
diff --git a/MaterialModels/FiniteStrain/Finite_VEVP/umat.f b/MaterialModels/FiniteStrain/Finite_VEVP/umat.f
index 8ebe81a..dd24885 100644
--- a/MaterialModels/FiniteStrain/Finite_VEVP/umat.f
+++ b/MaterialModels/FiniteStrain/Finite_VEVP/umat.f
@@ -33,7 +33,7 @@ C !--------------------------------------------------------------
C ! UMAT SUBROUTINE
C !--------------------------------------------------------------
- SUBROUTINE umat_vevp(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
+ SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,
3 DPRED,CMNAME,
@@ -42,24 +42,33 @@ C !--------------------------------------------------------------
6 CELENT,DFGRD0,DFGRD1,NOEL,NPT,
7 LAYER,KSPT,KSTEP,KINC)
-C INCLUDE 'ABA_PARAM.INC'
+ INCLUDE 'ABA_PARAM.INC'
CHARACTER*80 CMNAME
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),DFGRD0(3,3), DFGRD1(3,3)
+ 3 PROPS(NPROPS), COORDS(3),DROT(3,3),DFGRD0(3,3),
+ 4 DFGRD1(3,3)
REAL*8 SSE, SPD, SCD, RPL, DTIME, TEMP, CELENT, DRPLDT
REAL*8 PNEWDT, DTEMP, R, dR, ddR
+
+ IF (PROPS(1) .GT. 2.D0) THEN
- call VEVP(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
+ call VEVP(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
& NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1,
& PNEWDT)
- return
- end
+ return
+
+ ELSE
+ call NEOHOOK(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
+ & NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1)
+ return
+ END IF
+ end
C !--------------------------------------------------------------
@@ -83,7 +92,7 @@ C !--------------------------------------------------------------
REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3)
REAL(prec), INTENT(INOUT) :: statev(nstatv)
REAL(prec), INTENT(OUT) :: stress(ntens)
- REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens), PNEWDT
+ REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens), PNEWDT
!Declaration of local variables
INTEGER :: ii, jj, O6, order
@@ -101,7 +110,8 @@ C !--------------------------------------------------------------
! VE trial variables
REAL(prec) :: F_ve_tr(3, 3), C_ve_tr(3, 3), D_ve_tr(3, 3)
REAL(prec) :: E_ve_tr(3, 3), dlogC_ve_tr(3,3,3,3)
- REAL(prec) :: ddlogC_ve_tr(3,3,3,3,3,3), AA_tr(8,3,3)
+ !REAL(prec) :: ddlogC_ve_tr(3,3,3,3,3,3)
+ REAL(prec) :: AA_tr(8,3,3)
REAL(prec) :: BB_tr(8),GGe, KKe, kappa_tr(3, 3), S_tr(3, 3)
REAL(prec) :: temp1(3,3), tau_tr(3, 3)
! Variables for trial yield
@@ -111,7 +121,7 @@ C !--------------------------------------------------------------
REAL(prec) :: F_hat(6, 3, 3), J, F_ve_tr_hat(6, 3, 3)
REAL(prec) :: C_ve_tr_hat(6, 3, 3), D_ve_tr_hat(6,3,3)
REAL(prec) :: E_ve_tr_hat(6,3,3), dlogC_ve_tr_hat(6,3,3,3,3)
- REAL(prec) :: ddlogC_ve_tr_hat(6,3,3,3,3,3,3)
+ !REAL(prec) :: ddlogC_ve_tr_hat(6,3,3,3,3,3,3)
REAL(prec) :: AA_tr_hat(6,8,3,3)
REAL(prec) :: BB_tr_hat(6, 8), GGe_hat(6), KKe_hat(6)
REAL(prec) :: kappa_tr_hat(6, 3, 3), S_tr_hat(6, 3, 3)
@@ -124,7 +134,8 @@ C !--------------------------------------------------------------
REAL(prec) :: exp_GQ(3, 3), F_vp(3,3), F_vp_inv(3,3)
! VE corrected variables
REAL(prec) :: F_ve(3, 3), C_ve(3, 3), D_ve(3, 3), E_ve(3, 3)
- REAL(prec) :: dlogC_ve(3,3,3,3), ddlogC_ve(3,3,3,3,3,3)
+ REAL(prec) :: dlogC_ve(3,3,3,3)
+ !REAL(prec) :: ddlogC_ve(3,3,3,3,3,3)
REAL(prec) :: AA(8,3,3), BB(8), kappa(3, 3), b(3, 3)
REAL(prec) :: S(3, 3), temp2(3, 3), tau(3, 3)
! VP variables perturbated
@@ -143,14 +154,15 @@ C !--------------------------------------------------------------
REAL(prec) :: F_ve_hat(6,3,3), C_ve_hat(6,3,3)
REAL(prec) :: D_ve_hat(6,3,3), E_ve_hat(6,3,3)
REAL(prec) :: dlogC_ve_hat(6,3,3,3,3)
- REAL(prec) :: ddlogC_ve_hat(6,3,3,3,3,3,3), AA_hat(6,8,3,3)
+ REAL(prec) :: AA_hat(6,8,3,3)
+ !REAL(prec) :: ddlogC_ve_hat(6,3,3,3,3,3,3)
REAL(prec) :: BB_hat(6,8), kappa_hat(6,3,3), S_hat(6,3,3)
REAL(prec) :: temp2_hat(6,3,3), tau_hat(6,3,3)
REAL(prec) :: tau_hat_v(6,6)
! Tollerances
- REAL(prec), PARAMETER :: TOLL=0.001D0, TOLL_G=0.999D-11
- INTEGER, PARAMETER :: MAX_i=500
+ REAL(prec), PARAMETER :: TOLL=0.0001D0, TOLL_G=0.999D-6
+ INTEGER, PARAMETER :: MAX_i=100
!Define 2nd order identity tensor
data I_mat(1,:) /1.D0, 0.D0, 0.D0/
@@ -161,11 +173,11 @@ C !--------------------------------------------------------------
! for the first time step (initial condition) in FFTMAD.
! Can be commented out for ABAQUS, If you comment out, then
! in abaqus use the inp file to set initial values for F_vp.
-C IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN
-C STATEV(1) = 1.D0
-C STATEV(2) = 1.D0
-C STATEV(3) = 1.D0
-C END IF
+ IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN
+ STATEV(1) = 1.D0
+ STATEV(2) = 1.D0
+ STATEV(3) = 1.D0
+ END IF
! State variables at previous time step
CALL voit2mat(STATEV(1:9), F_vp_n(:,:))
@@ -255,7 +267,7 @@ C END IF
! Approximate VE log Strain (power series) at trial state
D_ve_tr(:,:) = C_ve_tr(:, :) - I_mat(:,:)
CALL approx_log(D_ve_tr(:,:), order, E_ve_tr(:, :),
- 1 dlogC_ve_tr(:,:,:,:), ddlogC_ve_tr(:,:,:,:,:,:))
+ 1 dlogC_ve_tr(:,:,:,:))
E_ve_tr(:, :) = 0.5D0 * E_ve_tr(:, :)
! Calculate the corotated kirchoff stress at trial state along
@@ -296,8 +308,7 @@ C END IF
! Approximate VE log Strain at trial state
D_ve_tr_hat(O6, :,:) = C_ve_tr_hat(O6, :, :) - I_mat(:,:)
CALL approx_log(D_ve_tr_hat(O6, :,:), order,
- 1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:),
- 2 ddlogC_ve_tr_hat(O6,:,:,:,:,:,:))
+ 1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:))
E_ve_tr_hat(O6,:, :) = 0.5D0 * E_ve_tr_hat(O6, :, :)
! Calculate the corotated kirchoff stress at trial state
@@ -375,7 +386,7 @@ C END IF
DO ii = 1, 6
DO jj = 1, 6
DDSDDE(ii, jj) = (1.D0 / (J * TOLL))
- 1 * (tau_tr_hat_v(jj, ii) - J*STRESS(ii))
+ 1 * (tau_tr_hat_v(jj, ii) - J*STRESS(ii))
END DO
END DO
@@ -400,14 +411,14 @@ C END IF
GQ(:,:) = GAMMA * Q(:,:)
CALL approx_exp(GQ(:,:), order, exp_GQ(:,:))
CALL mat2dot(exp_GQ(:,:), F_vp_n(:,:), F_vp(:,:))
- CALL matInv(F_vp(:,:), F_vp_inv(:,:))
+ CALL matInverse(F_vp(:,:), F_vp_inv(:,:))
CALL mat2dot(DFGRD1(:,:), F_vp_inv(:,:), F_ve(:, :))
CALL mat2Tdot(F_ve(:,:), F_ve(:,:), C_ve(:,:))
! Approximate VE log Strain (power series) at the corrected state
D_ve(:,:) = C_ve(:, :) - I_mat(:,:)
CALL approx_log(D_ve(:,:), order, E_ve(:, :),
- 1 dlogC_ve(:,:,:,:), ddlogC_ve(:,:,:,:,:,:))
+ 1 dlogC_ve(:,:,:,:))
E_ve(:, :) = 0.5D0 * E_ve(:, :)
! Calculate the corotated kirchoff stress at corrected state along
@@ -513,7 +524,7 @@ C END IF
CALL approx_exp(GQ_hat(O6,:,:), order, exp_GQ_hat(O6,:,:))
CALL mat2dot(exp_GQ_hat(O6, :,:), F_vp_n(:,:),
1 F_vp_hat(O6,:,:))
- CALL matInv(F_vp_hat(O6,:,:), F_vp_inv_hat(O6,:,:))
+ CALL matInverse(F_vp_hat(O6,:,:), F_vp_inv_hat(O6,:,:))
CALL mat2dot(F_hat(O6, :,:), F_vp_inv_hat(O6,:,:),
1 F_ve_hat(O6, :, :))
CALL mat2Tdot(F_ve_hat(O6,:,:), F_ve_hat(O6,:,:),
@@ -523,8 +534,7 @@ C END IF
! Approximate VE log Strain (power series) at the purturbated state
D_ve_hat(O6,:,:) = C_ve_hat(O6,:, :) - I_mat(:,:)
CALL approx_log(D_ve_hat(O6,:,:), order, E_ve_hat(O6,:, :),
- 1 dlogC_ve_hat(O6,:,:,:,:),
- 2 ddlogC_ve_hat(O6,:,:,:,:,:,:))
+ 1 dlogC_ve_hat(O6,:,:,:,:))
E_ve_hat(O6,:, :) = 0.5D0 * E_ve_hat(O6,:, :)
! Calculate the corotated kirchoff stress at purturbed state along
@@ -570,6 +580,119 @@ C END IF
RETURN
END SUBROUTINE VEVP
+C !--------------------------------------------------------------
+C ! UMAT SUBROUTINE FOR ISOTROPIC NEO HOOKIAN MATERIAL MODEL
+C !--------------------------------------------------------------
+
+ SUBROUTINE NEOHOOK(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,
+ 1 PROPS,NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1)
+
+
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: ntens,nprops,nstatv
+ INTEGER, INTENT(IN) :: kinc,kstep,noel
+ REAL(prec), INTENT(IN) :: stran(ntens), dstran(ntens)
+ REAL(prec), INTENT(IN) :: props(nprops), dtime
+ REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3)
+ REAL(prec), INTENT(INOUT) :: statev(nstatv)
+ REAL(prec), INTENT(OUT) :: stress(ntens)
+ REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens)
+
+ !List of internal variables
+ INTEGER :: ii, jj, mm, K1, ll, m2v(3, 3)
+ REAL(prec) :: I(6), B(6), B_bar(6), dev_B_bar(6)
+ REAL(prec) :: F_inv(3, 3)
+ REAL(prec) :: J, I_bar, fac1, fac2, fac3, fac4, cE, cnu
+
+ !Decleration of constants
+ REAL(prec), PARAMETER :: ZERO=0.D0, ONE=1.D0, TWO=2.D0
+ REAL(prec), PARAMETER :: THREE=3.D0, FOUR= 4.D0, SIX=6.D0
+ REAL(prec), PARAMETER :: NINE=9.D0
+
+
+ !Declare matrix to voit mapping
+ data m2v/ 1, 4, 5,
+ 1 4, 2, 6,
+ 2 5, 6, 3/
+
+ !2nd Order Identity
+ data I/ ONE, ONE, ONE, ZERO, ZERO, ZERO/
+
+ !Get material properties
+ cE=PROPS(1)
+ cnu=PROPS(2)
+
+ !Calculate determinant of deformation gradient
+ CALL determinant(DFGRD1(:,:), J)
+
+ !Calculate inverse of deformation gradient
+ CALL matInverse(DFGRD1(:,:), F_inv(:, :))
+
+ !Calculate finger tensor B
+ B(:) = ZERO
+ DO K1 = 1, 3
+ DO ii = 1, 3
+ DO jj = 1, 3
+ IF (ii .EQ. jj) THEN
+ B(m2v(ii, jj)) = B(m2v(ii, jj)) + DFGRD1(ii, K1)
+ 1 * DFGRD1(jj, K1)
+ ELSE
+ B(m2v(ii, jj)) = B(m2v(ii, jj)) + 0.5d0
+ 1 * DFGRD1(ii, K1) * DFGRD1(jj, K1)
+ END IF
+ END DO
+ END DO
+ END DO
+
+ B_bar(:) = B(:) * J ** (-TWO / THREE)
+ I_bar = B_bar(1) + B_bar(2) + B_bar(3)
+ dev_B_bar(:) = B_bar(:) - (ONE / THREE) * I_bar * I(:)
+
+ !Return Cauchy stress to Abaqus
+ STRESS(:) = 2 * cE * (J - ONE) * J * I(:)
+ 1 + 2 * cnu * dev_B_bar(:)
+
+ !Algorithmic constants for the Tangent
+ fac1 = -(FOUR / THREE) * cnu * J ** (-ONE)
+ fac2 = -(TWO * cE * (J - ONE) * J
+ 1 - (TWO/THREE) * cnu * I_bar) * J ** (-ONE)
+ fac3 = (TWO * cE * (J - ONE) * J
+ 1 + (FOUR / NINE) * cnu * I_bar + TWO * cE * J ** TWO)
+ 2 * J ** (-ONE)
+ fac4 = TWO * cnu * J ** (-ONE)
+
+ !Tangent for ABAQUS
+ DDSDDE(:,:) = ZERO
+ DO ii = 1, 3
+ DO jj = 1, 3
+ DO ll = 1, 3
+ DO mm = 1, 3
+ DDSDDE(m2v(ii, jj), m2v(ll, mm)) =
+ 1 DDSDDE(m2v(ii, jj), m2v(ll, mm))
+ 2 + fac1 * (B_bar(m2v(ii, jj)) * I(m2v(ll, mm))
+ 3 + I(m2v(ii, jj)) * B_bar(m2v(ll, mm)))
+ 4 + (fac4 - TWO * cnu * J**(-ONE))
+ 5 * B_bar(m2v(ii, ll)) * I(m2v(jj, mm))
+ 6 + TWO * cnu * (dev_B_bar(m2v(ii, ll))
+ 7 * I(m2v(jj, mm))
+ 8 + I(m2v(ii, ll)) * dev_B_bar(m2v(jj, mm)))
+ 9 + fac2 * I(m2v(ii, mm)) * I(m2v(ll, jj))
+ 1 + fac3 * I(m2v(ii, jj)) * I(m2v(ll, mm))
+ 2 + (fac2 + FOUR * cE * (J - ONE))
+ 3 * I(m2v(ii, ll)) * I(m2v(jj, mm))
+ END DO
+ END DO
+ END DO
+ END DO
+
+
+ RETURN
+ END SUBROUTINE NEOHOOK
+
!--------------------------------------------------------------
! Helper function to convert a voit array to matrix
!--------------------------------------------------------------
@@ -643,7 +766,7 @@ C END IF
!--------------------------------------------------------------
! Helper function to compute inverse of 3x3 matrix
!--------------------------------------------------------------
- SUBROUTINE matInv(matrix, inv)
+ SUBROUTINE matInverse(matrix, inv)
IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0)
@@ -668,7 +791,7 @@ C END IF
inv(:, :) = inv(:, :)/ J
RETURN
- END SUBROUTINE matInv
+ END SUBROUTINE matInverse
!------------------------------------------------------------
@@ -843,7 +966,7 @@ C END IF
! Compute inverse of int variable F_vp
- CALL matInv(F_vp(:, :), F_vp_inv(:, :))
+ CALL matInverse(F_vp(:, :), F_vp_inv(:, :))
! Compute F_ve
CALL mat2dot(F(:, :), F_vp_inv(:, :), F_ve(:, :))
@@ -929,7 +1052,7 @@ C END IF
! Power series approximation
!--------------------------------------------------------------
- SUBROUTINE approx_log(AA, order, logAA, dlogAA, ddlogAA)
+ SUBROUTINE approx_log(AA, order, logAA, dlogAA)!, ddlogAA)
IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0)
INTEGER, PARAMETER :: prec=double
@@ -937,13 +1060,14 @@ C END IF
INTEGER, INTENT(IN) :: order
REAL(prec), INTENT(IN) :: AA(3, 3)
REAL(prec), INTENT(OUT) :: logAA(3, 3), dlogAA(3, 3, 3, 3)
- REAL(prec), INTENT(OUT) :: ddlogAA(3, 3, 3, 3, 3, 3)
+ !REAL(prec), INTENT(OUT) :: ddlogAA(3, 3, 3, 3, 3, 3)
! List of internal variables
- INTEGER :: ii, jj, nn, mm, ll, rr, ss, pp, qq
+ INTEGER :: ii, jj, nn, mm, ll, rr, ss
REAL(prec) :: coeffs(order + 1), I_mat(3, 3)
REAL(prec) :: FF(3, 3), dFF(3, 3, 3, 3)
- REAL(prec) :: ddFF(3, 3, 3, 3, 3, 3), II_mat(3, 3, 3, 3)
+ !REAL(prec) :: ddFF(3, 3, 3, 3, 3, 3)
+ REAL(prec) :: II_mat(3, 3, 3, 3)
REAL(prec) :: temp1(3, 3, 3, 3), temp2(3, 3, 3, 3)
! Define 2nd order identity tensor
@@ -990,12 +1114,12 @@ C END IF
END DO
dlogAA(:, :, :, :) = 0.D0
- ddlogAA(:, :, :, :, :, :) = 0.D0
+ !ddlogAA(:, :, :, :, :, :) = 0.D0
DO ii = 1, order
FF(:, :) = 0.D0
dFF(:, :, :, :) = 0.D0
- ddFF(:, :, :, :, :, :) = 0.D0
+ !ddFF(:, :, :, :, :, :) = 0.D0
nn = order + 1 - ii
@@ -1012,14 +1136,14 @@ C END IF
1 + dlogAA(jj, ll, rr, ss) * AA(ll, mm)
2 + logAA(jj, ll) * II_mat(ll, mm, rr, ss)
- DO pp = 1, 3
- DO qq = 1, 3
- ddFF(jj, mm, rr, ss, pp, qq) = ddFF(jj, mm, rr, ss, pp, qq)
- 1 + ddlogAA(jj, ll, rr, ss, pp, qq) * AA(ll, mm)
- 2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq)
- 3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss)
- END DO
- END DO
+ ! DO pp = 1, 3
+ ! DO qq = 1, 3
+ ! ddFF(jj, mm, rr, ss, pp, qq) = ddFF(jj, mm, rr, ss, pp, qq)
+ ! 1 + ddlogAA(jj, ll, rr, ss, pp, qq) * AA(ll, mm)
+ ! 2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq)
+ ! 3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss)
+ ! END DO
+ ! END DO
END DO
END DO
@@ -1030,7 +1154,7 @@ C END IF
END DO
logAA(:, :) = FF(:, :)
dlogAA(:,:,:,:) = dFF(:,:,:,:)
- ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :)
+ !ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :)
END DO
END SUBROUTINE approx_log
@@ -1399,8 +1523,8 @@ C END IF
gma_i = gma + Dgma
! Update hardening variables on gma_i
- CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c, sigma_t,
- 1 HHc, HHt, HHb)
+ CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c,
+ 1 sigma_t, HHc, HHt, HHb)
! Update Drucker Pragger coefficients along with ecc factor m
CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
@@ -1421,7 +1545,7 @@ C END IF
gma = gma_i
- ! Debuging statements for N.R. Comment out for big sims
+ ! !Debuging statements for N.R. Comment out for big sims
! WRITE(*,*) "Iter : ", iter_G
! WRITE(*,*) "GAMMA : ", GAMMA
! WRITE(*,*) "RESI : ", ABS(F)
--
GitLab