Skip to content
Snippets Groups Projects
Commit f8cf6170 authored by Mohib's avatar Mohib
Browse files

[Refractor] - Numerical Tangent as an exteral function call

parent 6902956e
No related branches found
No related tags found
1 merge request!7[Refractor] - Numerical Tangent as an exteral function call
Hide whitespace changes
Inline Side-by-side
MaterialModels/FiniteStrain/Finite_VEVP/umat.f
+ 172
− 48
View file @ f8cf6170
...@@ -33,7 +33,7 @@ C !-------------------------------------------------------------- ...@@ -33,7 +33,7 @@ C !--------------------------------------------------------------
C ! UMAT SUBROUTINE C ! UMAT SUBROUTINE
C !-------------------------------------------------------------- C !--------------------------------------------------------------
SUBROUTINE umat_vevp(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT, 1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,
3 DPRED,CMNAME, 3 DPRED,CMNAME,
...@@ -42,24 +42,33 @@ C !-------------------------------------------------------------- ...@@ -42,24 +42,33 @@ C !--------------------------------------------------------------
6 CELENT,DFGRD0,DFGRD1,NOEL,NPT, 6 CELENT,DFGRD0,DFGRD1,NOEL,NPT,
7 LAYER,KSPT,KSTEP,KINC) 7 LAYER,KSPT,KSTEP,KINC)
C INCLUDE 'ABA_PARAM.INC' INCLUDE 'ABA_PARAM.INC'
CHARACTER*80 CMNAME CHARACTER*80 CMNAME
REAL*8 STRESS(NTENS),STATEV(NSTATV), REAL*8 STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 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 SSE, SPD, SCD, RPL, DTIME, TEMP, CELENT, DRPLDT
REAL*8 PNEWDT, DTEMP, R, dR, ddR 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, & NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1,
& PNEWDT) & PNEWDT)
return 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 !-------------------------------------------------------------- C !--------------------------------------------------------------
...@@ -83,7 +92,7 @@ C !-------------------------------------------------------------- ...@@ -83,7 +92,7 @@ C !--------------------------------------------------------------
REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3) REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3)
REAL(prec), INTENT(INOUT) :: statev(nstatv) REAL(prec), INTENT(INOUT) :: statev(nstatv)
REAL(prec), INTENT(OUT) :: stress(ntens) 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 !Declaration of local variables
INTEGER :: ii, jj, O6, order INTEGER :: ii, jj, O6, order
...@@ -101,7 +110,8 @@ C !-------------------------------------------------------------- ...@@ -101,7 +110,8 @@ C !--------------------------------------------------------------
! VE trial variables ! VE trial variables
REAL(prec) :: F_ve_tr(3, 3), C_ve_tr(3, 3), D_ve_tr(3, 3) 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) :: 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) :: BB_tr(8),GGe, KKe, kappa_tr(3, 3), S_tr(3, 3)
REAL(prec) :: temp1(3,3), tau_tr(3, 3) REAL(prec) :: temp1(3,3), tau_tr(3, 3)
! Variables for trial yield ! Variables for trial yield
...@@ -111,7 +121,7 @@ C !-------------------------------------------------------------- ...@@ -111,7 +121,7 @@ C !--------------------------------------------------------------
REAL(prec) :: F_hat(6, 3, 3), J, F_ve_tr_hat(6, 3, 3) 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) :: 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) :: 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) :: AA_tr_hat(6,8,3,3)
REAL(prec) :: BB_tr_hat(6, 8), GGe_hat(6), KKe_hat(6) 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) REAL(prec) :: kappa_tr_hat(6, 3, 3), S_tr_hat(6, 3, 3)
...@@ -124,7 +134,8 @@ C !-------------------------------------------------------------- ...@@ -124,7 +134,8 @@ C !--------------------------------------------------------------
REAL(prec) :: exp_GQ(3, 3), F_vp(3,3), F_vp_inv(3,3) REAL(prec) :: exp_GQ(3, 3), F_vp(3,3), F_vp_inv(3,3)
! VE corrected variables ! VE corrected variables
REAL(prec) :: F_ve(3, 3), C_ve(3, 3), D_ve(3, 3), E_ve(3, 3) 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) :: AA(8,3,3), BB(8), kappa(3, 3), b(3, 3)
REAL(prec) :: S(3, 3), temp2(3, 3), tau(3, 3) REAL(prec) :: S(3, 3), temp2(3, 3), tau(3, 3)
! VP variables perturbated ! VP variables perturbated
...@@ -143,14 +154,15 @@ C !-------------------------------------------------------------- ...@@ -143,14 +154,15 @@ C !--------------------------------------------------------------
REAL(prec) :: F_ve_hat(6,3,3), C_ve_hat(6,3,3) 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) :: D_ve_hat(6,3,3), E_ve_hat(6,3,3)
REAL(prec) :: dlogC_ve_hat(6,3,3,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) :: 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) :: temp2_hat(6,3,3), tau_hat(6,3,3)
REAL(prec) :: tau_hat_v(6,6) REAL(prec) :: tau_hat_v(6,6)
! Tollerances ! Tollerances
REAL(prec), PARAMETER :: TOLL=0.001D0, TOLL_G=0.999D-11 REAL(prec), PARAMETER :: TOLL=0.0001D0, TOLL_G=0.999D-6
INTEGER, PARAMETER :: MAX_i=500 INTEGER, PARAMETER :: MAX_i=100
!Define 2nd order identity tensor !Define 2nd order identity tensor
data I_mat(1,:) /1.D0, 0.D0, 0.D0/ data I_mat(1,:) /1.D0, 0.D0, 0.D0/
...@@ -161,11 +173,11 @@ C !-------------------------------------------------------------- ...@@ -161,11 +173,11 @@ C !--------------------------------------------------------------
! for the first time step (initial condition) in FFTMAD. ! for the first time step (initial condition) in FFTMAD.
! Can be commented out for ABAQUS, If you comment out, then ! Can be commented out for ABAQUS, If you comment out, then
! in abaqus use the inp file to set initial values for F_vp. ! in abaqus use the inp file to set initial values for F_vp.
C IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN
C STATEV(1) = 1.D0 STATEV(1) = 1.D0
C STATEV(2) = 1.D0 STATEV(2) = 1.D0
C STATEV(3) = 1.D0 STATEV(3) = 1.D0
C END IF END IF
! State variables at previous time step ! State variables at previous time step
CALL voit2mat(STATEV(1:9), F_vp_n(:,:)) CALL voit2mat(STATEV(1:9), F_vp_n(:,:))
...@@ -255,7 +267,7 @@ C END IF ...@@ -255,7 +267,7 @@ C END IF
! Approximate VE log Strain (power series) at trial state ! Approximate VE log Strain (power series) at trial state
D_ve_tr(:,:) = C_ve_tr(:, :) - I_mat(:,:) D_ve_tr(:,:) = C_ve_tr(:, :) - I_mat(:,:)
CALL approx_log(D_ve_tr(:,:), order, E_ve_tr(:, :), 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(:, :) E_ve_tr(:, :) = 0.5D0 * E_ve_tr(:, :)
! Calculate the corotated kirchoff stress at trial state along ! Calculate the corotated kirchoff stress at trial state along
...@@ -296,8 +308,7 @@ C END IF ...@@ -296,8 +308,7 @@ C END IF
! Approximate VE log Strain at trial state ! Approximate VE log Strain at trial state
D_ve_tr_hat(O6, :,:) = C_ve_tr_hat(O6, :, :) - I_mat(:,:) D_ve_tr_hat(O6, :,:) = C_ve_tr_hat(O6, :, :) - I_mat(:,:)
CALL approx_log(D_ve_tr_hat(O6, :,:), order, CALL approx_log(D_ve_tr_hat(O6, :,:), order,
1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:), 1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:))
2 ddlogC_ve_tr_hat(O6,:,:,:,:,:,:))
E_ve_tr_hat(O6,:, :) = 0.5D0 * E_ve_tr_hat(O6, :, :) E_ve_tr_hat(O6,:, :) = 0.5D0 * E_ve_tr_hat(O6, :, :)
! Calculate the corotated kirchoff stress at trial state ! Calculate the corotated kirchoff stress at trial state
...@@ -375,7 +386,7 @@ C END IF ...@@ -375,7 +386,7 @@ C END IF
DO ii = 1, 6 DO ii = 1, 6
DO jj = 1, 6 DO jj = 1, 6
DDSDDE(ii, jj) = (1.D0 / (J * TOLL)) 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
END DO END DO
...@@ -400,14 +411,14 @@ C END IF ...@@ -400,14 +411,14 @@ C END IF
GQ(:,:) = GAMMA * Q(:,:) GQ(:,:) = GAMMA * Q(:,:)
CALL approx_exp(GQ(:,:), order, exp_GQ(:,:)) CALL approx_exp(GQ(:,:), order, exp_GQ(:,:))
CALL mat2dot(exp_GQ(:,:), F_vp_n(:,:), F_vp(:,:)) 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 mat2dot(DFGRD1(:,:), F_vp_inv(:,:), F_ve(:, :))
CALL mat2Tdot(F_ve(:,:), F_ve(:,:), C_ve(:,:)) CALL mat2Tdot(F_ve(:,:), F_ve(:,:), C_ve(:,:))
! Approximate VE log Strain (power series) at the corrected state ! Approximate VE log Strain (power series) at the corrected state
D_ve(:,:) = C_ve(:, :) - I_mat(:,:) D_ve(:,:) = C_ve(:, :) - I_mat(:,:)
CALL approx_log(D_ve(:,:), order, E_ve(:, :), CALL approx_log(D_ve(:,:), order, E_ve(:, :),
1 dlogC_ve(:,:,:,:), ddlogC_ve(:,:,:,:,:,:)) 1 dlogC_ve(:,:,:,:))
E_ve(:, :) = 0.5D0 * E_ve(:, :) E_ve(:, :) = 0.5D0 * E_ve(:, :)
! Calculate the corotated kirchoff stress at corrected state along ! Calculate the corotated kirchoff stress at corrected state along
...@@ -513,7 +524,7 @@ C END IF ...@@ -513,7 +524,7 @@ C END IF
CALL approx_exp(GQ_hat(O6,:,:), order, exp_GQ_hat(O6,:,:)) CALL approx_exp(GQ_hat(O6,:,:), order, exp_GQ_hat(O6,:,:))
CALL mat2dot(exp_GQ_hat(O6, :,:), F_vp_n(:,:), CALL mat2dot(exp_GQ_hat(O6, :,:), F_vp_n(:,:),
1 F_vp_hat(O6,:,:)) 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,:,:), CALL mat2dot(F_hat(O6, :,:), F_vp_inv_hat(O6,:,:),
1 F_ve_hat(O6, :, :)) 1 F_ve_hat(O6, :, :))
CALL mat2Tdot(F_ve_hat(O6,:,:), F_ve_hat(O6,:,:), CALL mat2Tdot(F_ve_hat(O6,:,:), F_ve_hat(O6,:,:),
...@@ -523,8 +534,7 @@ C END IF ...@@ -523,8 +534,7 @@ C END IF
! Approximate VE log Strain (power series) at the purturbated state ! Approximate VE log Strain (power series) at the purturbated state
D_ve_hat(O6,:,:) = C_ve_hat(O6,:, :) - I_mat(:,:) D_ve_hat(O6,:,:) = C_ve_hat(O6,:, :) - I_mat(:,:)
CALL approx_log(D_ve_hat(O6,:,:), order, E_ve_hat(O6,:, :), CALL approx_log(D_ve_hat(O6,:,:), order, E_ve_hat(O6,:, :),
1 dlogC_ve_hat(O6,:,:,:,:), 1 dlogC_ve_hat(O6,:,:,:,:))
2 ddlogC_ve_hat(O6,:,:,:,:,:,:))
E_ve_hat(O6,:, :) = 0.5D0 * E_ve_hat(O6,:, :) E_ve_hat(O6,:, :) = 0.5D0 * E_ve_hat(O6,:, :)
! Calculate the corotated kirchoff stress at purturbed state along ! Calculate the corotated kirchoff stress at purturbed state along
...@@ -570,6 +580,119 @@ C END IF ...@@ -570,6 +580,119 @@ C END IF
RETURN RETURN
END SUBROUTINE VEVP 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 ! Helper function to convert a voit array to matrix
!-------------------------------------------------------------- !--------------------------------------------------------------
...@@ -643,7 +766,7 @@ C END IF ...@@ -643,7 +766,7 @@ C END IF
!-------------------------------------------------------------- !--------------------------------------------------------------
! Helper function to compute inverse of 3x3 matrix ! Helper function to compute inverse of 3x3 matrix
!-------------------------------------------------------------- !--------------------------------------------------------------
SUBROUTINE matInv(matrix, inv) SUBROUTINE matInverse(matrix, inv)
IMPLICIT NONE IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0) INTEGER, PARAMETER :: double=kind(1.d0)
...@@ -668,7 +791,7 @@ C END IF ...@@ -668,7 +791,7 @@ C END IF
inv(:, :) = inv(:, :)/ J inv(:, :) = inv(:, :)/ J
RETURN RETURN
END SUBROUTINE matInv END SUBROUTINE matInverse
!------------------------------------------------------------ !------------------------------------------------------------
...@@ -843,7 +966,7 @@ C END IF ...@@ -843,7 +966,7 @@ C END IF
! Compute inverse of int variable F_vp ! Compute inverse of int variable F_vp
CALL matInv(F_vp(:, :), F_vp_inv(:, :)) CALL matInverse(F_vp(:, :), F_vp_inv(:, :))
! Compute F_ve ! Compute F_ve
CALL mat2dot(F(:, :), F_vp_inv(:, :), F_ve(:, :)) CALL mat2dot(F(:, :), F_vp_inv(:, :), F_ve(:, :))
...@@ -929,7 +1052,7 @@ C END IF ...@@ -929,7 +1052,7 @@ C END IF
! Power series approximation ! Power series approximation
!-------------------------------------------------------------- !--------------------------------------------------------------
SUBROUTINE approx_log(AA, order, logAA, dlogAA, ddlogAA) SUBROUTINE approx_log(AA, order, logAA, dlogAA)!, ddlogAA)
IMPLICIT NONE IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0) INTEGER, PARAMETER :: double=kind(1.d0)
INTEGER, PARAMETER :: prec=double INTEGER, PARAMETER :: prec=double
...@@ -937,13 +1060,14 @@ C END IF ...@@ -937,13 +1060,14 @@ C END IF
INTEGER, INTENT(IN) :: order INTEGER, INTENT(IN) :: order
REAL(prec), INTENT(IN) :: AA(3, 3) REAL(prec), INTENT(IN) :: AA(3, 3)
REAL(prec), INTENT(OUT) :: logAA(3, 3), dlogAA(3, 3, 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 ! 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) :: coeffs(order + 1), I_mat(3, 3)
REAL(prec) :: FF(3, 3), dFF(3, 3, 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) REAL(prec) :: temp1(3, 3, 3, 3), temp2(3, 3, 3, 3)
! Define 2nd order identity tensor ! Define 2nd order identity tensor
...@@ -990,12 +1114,12 @@ C END IF ...@@ -990,12 +1114,12 @@ C END IF
END DO END DO
dlogAA(:, :, :, :) = 0.D0 dlogAA(:, :, :, :) = 0.D0
ddlogAA(:, :, :, :, :, :) = 0.D0 !ddlogAA(:, :, :, :, :, :) = 0.D0
DO ii = 1, order DO ii = 1, order
FF(:, :) = 0.D0 FF(:, :) = 0.D0
dFF(:, :, :, :) = 0.D0 dFF(:, :, :, :) = 0.D0
ddFF(:, :, :, :, :, :) = 0.D0 !ddFF(:, :, :, :, :, :) = 0.D0
nn = order + 1 - ii nn = order + 1 - ii
...@@ -1012,14 +1136,14 @@ C END IF ...@@ -1012,14 +1136,14 @@ C END IF
1 + dlogAA(jj, ll, rr, ss) * AA(ll, mm) 1 + dlogAA(jj, ll, rr, ss) * AA(ll, mm)
2 + logAA(jj, ll) * II_mat(ll, mm, rr, ss) 2 + logAA(jj, ll) * II_mat(ll, mm, rr, ss)
DO pp = 1, 3 ! DO pp = 1, 3
DO qq = 1, 3 ! DO qq = 1, 3
ddFF(jj, mm, rr, ss, pp, qq) = ddFF(jj, mm, rr, ss, pp, qq) ! 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) ! 1 + ddlogAA(jj, ll, rr, ss, pp, qq) * AA(ll, mm)
2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq) ! 2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq)
3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss) ! 3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss)
END DO ! END DO
END DO ! END DO
END DO END DO
END DO END DO
...@@ -1030,7 +1154,7 @@ C END IF ...@@ -1030,7 +1154,7 @@ C END IF
END DO END DO
logAA(:, :) = FF(:, :) logAA(:, :) = FF(:, :)
dlogAA(:,:,:,:) = dFF(:,:,:,:) dlogAA(:,:,:,:) = dFF(:,:,:,:)
ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :) !ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :)
END DO END DO
END SUBROUTINE approx_log END SUBROUTINE approx_log
...@@ -1399,8 +1523,8 @@ C END IF ...@@ -1399,8 +1523,8 @@ C END IF
gma_i = gma + Dgma gma_i = gma + Dgma
! Update hardening variables on gma_i ! Update hardening variables on gma_i
CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c, sigma_t, CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c,
1 HHc, HHt, HHb) 1 sigma_t, HHc, HHt, HHb)
! Update Drucker Pragger coefficients along with ecc factor m ! Update Drucker Pragger coefficients along with ecc factor m
CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2) CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
...@@ -1421,7 +1545,7 @@ C END IF ...@@ -1421,7 +1545,7 @@ C END IF
gma = gma_i 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(*,*) "Iter : ", iter_G
! WRITE(*,*) "GAMMA : ", GAMMA ! WRITE(*,*) "GAMMA : ", GAMMA
! WRITE(*,*) "RESI : ", ABS(F) ! WRITE(*,*) "RESI : ", ABS(F)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment