Skip to content
Snippets Groups Projects
Commit 83566311 authored by Van Dung NGUYEN's avatar Van Dung NGUYEN
Browse files
parents 2ed3753c 01d47e30
No related branches found
No related tags found
No related merge requests found
......@@ -48,18 +48,27 @@ C INCLUDE 'ABA_PARAM.INC'
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,
& NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1,
& PNEWDT)
return
end
ELSE
call NEOHOOK(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
& NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1)
return
END IF
end
C !--------------------------------------------------------------
......@@ -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,7 +154,8 @@ 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)
......@@ -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
......@@ -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)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment