From 2d0c9e03c84488dbf3f34a712eff1acbd35f0a37 Mon Sep 17 00:00:00 2001
From: Pratheek Shanthraj
Date: Fri, 15 Aug 2014 09:09:31 +0000
Subject: [PATCH] all constitutive models coupled with damage
---
code/constitutive.f90 | 42 ++++++++++++++++-------------
code/constitutive_j2.f90 | 12 +++------
code/constitutive_phenopowerlaw.f90 | 30 +++++++++------------
3 files changed, 40 insertions(+), 44 deletions(-)
diff --git a/code/constitutive.f90 b/code/constitutive.f90
index bad830b82..16e0dae99 100644
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@@ -397,17 +397,17 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
Lp = 0.0_pReal
dLp_dTstar = 0.0_pReal
case (PLASTICITY_J2_ID)
- call constitutive_j2_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
+ call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v/constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
- call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
+ call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v/constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
- call constitutive_nonlocal_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature, ip,el)
+ call constitutive_nonlocal_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v/constitutive_damageValue(ipc,ip,el),temperature,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
- call constitutive_dislotwin_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
+ call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v/constitutive_damageValue(ipc,ip,el),temperature,ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID)
- call constitutive_dislokmc_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
+ call constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v/constitutive_damageValue(ipc,ip,el),temperature,ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
- call constitutive_titanmod_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
+ call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v/constitutive_damageValue(ipc,ip,el),temperature,ipc,ip,el)
end select
@@ -492,18 +492,18 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
pressure = math_trace33(T)/3.0_pReal
damage = constitutive_damageValue(ipc,ip,el)
T = damage*T
- if (pressure < 0.0_pReal) T = T + (1.0_pReal - damage)*pressure*math_I3
+! if (pressure < 0.0_pReal) T = T + (1.0_pReal - damage)*pressure*math_I3
dT_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) &
dT_dFe(i,j,k,l) = damage*sum(C(i,j,l,1:3)*Fe(k,1:3)) ! dT*_ij/dFe_kl
- if (pressure < 0.0_pReal) then
- do i=1_pInt, 3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do j=1_pInt,3_pInt
- dT_dFe(i,i,k,l) = dT_dFe(i,i,k,l) + &
- (1.0_pReal - damage)*math_trace33(C(1:3,1:3,l,j))*Fe(k,j)/3.0_pReal
- enddo; enddo; enddo; enddo
- endif
+! if (pressure < 0.0_pReal) then
+! do i=1_pInt, 3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do j=1_pInt,3_pInt
+! dT_dFe(i,i,k,l) = dT_dFe(i,i,k,l) + &
+! (1.0_pReal - damage)*math_trace33(C(1:3,1:3,l,j))*Fe(k,j)/3.0_pReal
+! enddo; enddo; enddo; enddo
+! endif
end subroutine constitutive_hooke_TandItsTangent
@@ -721,17 +721,21 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_TITANMOD_ID)
- constitutive_postResults = constitutive_titanmod_postResults (ipc,ip,el)
+ constitutive_postResults = constitutive_titanmod_postResults(ipc,ip,el)
case (PLASTICITY_J2_ID)
- constitutive_postResults= constitutive_j2_postResults (Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
+ constitutive_postResults= constitutive_j2_postResults(Tstar_v/constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
- constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
+ constitutive_postResults = &
+ constitutive_phenopowerlaw_postResults(Tstar_v/constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
- constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
+ constitutive_postResults = &
+ constitutive_dislotwin_postResults(Tstar_v/constitutive_damageValue(ipc,ip,el),Temperature,ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID)
- constitutive_postResults = constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
+ constitutive_postResults = &
+ constitutive_dislokmc_postResults(Tstar_v/constitutive_damageValue(ipc,ip,el),Temperature,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
- constitutive_postResults = constitutive_nonlocal_postResults (Tstar_v,FeArray, ip,el)
+ constitutive_postResults = &
+ constitutive_nonlocal_postResults (Tstar_v/constitutive_damageValue(ipc,ip,el),FeArray,ip,el)
end select
end function constitutive_postResults
diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90
index 23e53c406..04a1ae7af 100644
--- a/code/constitutive_j2.f90
+++ b/code/constitutive_j2.f90
@@ -343,7 +343,7 @@ end subroutine constitutive_j2_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
-subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip,el)
+subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
@@ -368,8 +368,6 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal), intent(in) :: &
- damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@@ -397,7 +395,7 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = constitutive_j2_gdot0(instance) &
- * (sqrt(1.5_pReal) * norm_Tstar_dev / (damage * constitutive_j2_fTaylor(instance) * &
+ * (sqrt(1.5_pReal) * norm_Tstar_dev / (constitutive_j2_fTaylor(instance) * &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))) &
**constitutive_j2_n(instance)
@@ -503,7 +501,7 @@ end subroutine constitutive_j2_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
-function constitutive_j2_postResults(Tstar_v,damage,ipc,ip,el)
+function constitutive_j2_postResults(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use mesh, only: &
@@ -520,8 +518,6 @@ function constitutive_j2_postResults(Tstar_v,damage,ipc,ip,el)
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal), intent(in) :: &
- damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@@ -562,7 +558,7 @@ function constitutive_j2_postResults(Tstar_v,damage,ipc,ip,el)
constitutive_j2_postResults(c+1_pInt) = &
constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!----------------------------------------------------------------------------------
- (damage * constitutive_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** constitutive_j2_n(instance)
+ (constitutive_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** constitutive_j2_n(instance)
c = c + 1_pInt
end select
enddo outputsLoop
diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90
index ce93bfcab..d8c3f93d8 100644
--- a/code/constitutive_phenopowerlaw.f90
+++ b/code/constitutive_phenopowerlaw.f90
@@ -662,7 +662,7 @@ end subroutine constitutive_phenopowerlaw_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
-subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip,el)
+subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use prec, only: &
p_vec
use math, only: &
@@ -696,8 +696,6 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,da
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal), intent(in) :: &
- damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@@ -755,11 +753,11 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,da
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo
gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
- ((abs(tau_slip_pos(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))*&
+ ((abs(tau_slip_pos(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_pos(j))
gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
- ((abs(tau_slip_neg(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))*&
+ ((abs(tau_slip_neg(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_neg(j))
Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
@@ -797,7 +795,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,da
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
- (abs(tau_twin(j))/(damage*plasticState(ph)%state(nSlip+j,of)))**&
+ (abs(tau_twin(j))/plasticState(ph)%state(nSlip+j,of))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
@@ -918,8 +916,8 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
- ((abs(tau_slip_pos(j))/(plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance) &
- +(abs(tau_slip_neg(j))/(plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))&
+ ((abs(tau_slip_pos(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) &
+ +(abs(tau_slip_neg(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos(j))
enddo
enddo slipFamiliesLoop1
@@ -940,7 +938,7 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
- (abs(tau_twin(j))/(plasticState(ph)%state(nslip+j,of)))**&
+ (abs(tau_twin(j))/plasticState(ph)%state(nslip+j,of))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
enddo
enddo twinFamiliesLoop1
@@ -951,7 +949,7 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt
- plasticState(ph)%dotState(j,of) = & ! evolution of slip resistance j
+ plasticState(ph)%dotState(j,of) = & ! evolution of slip resistance j
c_SlipSlip * left_SlipSlip(j) * &
dot_product(constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(j,1:nSlip,instance), &
right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor
@@ -966,7 +964,7 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
j = 0_pInt
twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
- index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
+ index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt
plasticState(ph)%dotState(j+nSlip,of) = & ! evolution of twin resistance j
@@ -990,7 +988,7 @@ end subroutine constitutive_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
-function constitutive_phenopowerlaw_postResults(Tstar_v,damage,ipc,ip,el)
+function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
use prec, only: &
p_vec
use mesh, only: &
@@ -1018,8 +1016,6 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,damage,ipc,ip,el)
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal), intent(in) :: &
- damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@@ -1077,8 +1073,8 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,damage,ipc,ip,el)
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
- ((abs(tau_slip_pos)/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance) &
- +(abs(tau_slip_neg)/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))&
+ ((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) &
+ +(abs(tau_slip_neg)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos)
enddo
@@ -1120,7 +1116,7 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,damage,ipc,ip,el)
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
- (abs(tau)/(damage*plasticState(ph)%state(j+nSlip,of)))**&
+ (abs(tau)/plasticState(ph)%state(j+nSlip,of))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau))
enddo
enddo twinFamiliesLoop1