223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

changed state parsing for local models (and for delta_state) such that only the needed part of the state array (for the given material point) is used

This commit is contained in:
Martin Diehl 2014-03-13 06:43:49 +00:00
parent 7c79b31f6c
commit 3f7a389ff7
6 changed files with 316 additions and 316 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

97
code/constitutive.f90
View File

@ -192,8 +192,8 @@ subroutine constitutive_init
select case(phase_plasticity(phase)) ! split per constititution
case (PLASTICITY_NONE_ID)
outputName = PLASTICITY_NONE_label
thisOutput => NULL() ! constitutive_none_output
thisSize => NULL() ! constitutive_none_sizePostResult
thisOutput => null() ! constitutive_none_output
thisSize => null() ! constitutive_none_sizePostResult
case (PLASTICITY_J2_ID)
outputName = PLASTICITY_J2_label
thisOutput => constitutive_j2_output
@ -506,11 +506,11 @@ pure function constitutive_homogenizedC(ipc,ip,el)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID)
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state,ipc,ip,el)
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state(ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state,ipc,ip,el)
constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state(ipc,ip,el), &
ipc,ip,el)
case default
constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase(ipc,ip,el))
@ -550,11 +550,11 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID)
call constitutive_dislotwin_microstructure(temperature,constitutive_state,ipc,ip,el)
call constitutive_dislotwin_microstructure(temperature,constitutive_state(ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
call constitutive_titanmod_microstructure(temperature,constitutive_state,ipc,ip,el)
call constitutive_titanmod_microstructure(temperature,constitutive_state(ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el)
@ -610,20 +610,20 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
dLp_dTstar = math_identity2nd(9)
case (PLASTICITY_J2_ID)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el)
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
temperature,constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el)
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
temperature,constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, constitutive_state(ipc,ip,el), ipc,ip,el)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, &
temperature, constitutive_state(ipc,ip,el), ipc,ip,el)
end select
end subroutine constitutive_LpAndItsTangent
@ -700,7 +700,8 @@ end subroutine constitutive_hooke_TandItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, subfrac, ipc, ip, el)
subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, subdt, subfracArray,&
ipc, ip, el)
use prec, only: &
pLongInt
use debug, only: &
@ -742,10 +743,10 @@ subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, sub
Temperature, &
subdt !< timestep
real(pReal), intent(in), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
subfrac !< subfraction of timestep
subfracArray !< subfraction of timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Fe, & !< elastic deformation gradient
Fp !< plastic deformation gradient
FeArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
integer(pLongInt) :: &
@ -762,20 +763,21 @@ subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, sub
constitutive_dotState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway?
case (PLASTICITY_J2_ID)
constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,&
constitutive_state(ipc,ip,el), ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,&
constitutive_state(ipc,ip,el), ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el), ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el), ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, constitutive_state, &
constitutive_state0, subdt, subfrac, ipc, ip, el)
constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, &
Temperature, constitutive_state, constitutive_state0, subdt, &
subfracArray, ipc, ip, el)
end select
@ -830,7 +832,8 @@ subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONLOCAL_ID)
call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),constitutive_state, Tstar_v,ipc,ip,el)
call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),&
constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el)
case default
constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway?
@ -853,7 +856,7 @@ end subroutine constitutive_collectDeltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results
!--------------------------------------------------------------------------------------------------
function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el)
function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
@ -888,7 +891,7 @@ function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el)
real(pReal), intent(in) :: &
temperature
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Fe !< elastic deformation gradient
FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
@ -899,20 +902,20 @@ function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el)
case (PLASTICITY_NONE_ID)
case (PLASTICITY_TITANMOD_ID)
constitutive_postResults = constitutive_titanmod_postResults(constitutive_state,ipc,ip,el)
constitutive_postResults = constitutive_titanmod_postResults(&
constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_J2_ID)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,constitutive_state,ipc,ip,el)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,&
constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,constitutive_state,ipc,ip,el)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,&
constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, constitutive_state, &
constitutive_dotstate(ipc,ip,el), ipc, ip, el)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, FeArray, &
constitutive_state, constitutive_dotstate(ipc,ip,el), ipc, ip, el)
end select
end function constitutive_postResults

170
code/constitutive_dislotwin.f90
View File

@ -939,7 +939,7 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
integer(pInt) :: instance,ns,nt,i,phase
@ -952,13 +952,13 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el)
nt = constitutive_dislotwin_totalNtwin(instance)
!* Total twin volume fraction
sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
!* Homogenized elasticity matrix
constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*lattice_C66(1:6,1:6,phase)
do i=1_pInt,nt
constitutive_dislotwin_homogenizedC = &
constitutive_dislotwin_homogenizedC + state(ipc,ip,el)%p(3_pInt*ns+i)*lattice_C66(1:6,1:6,phase)
constitutive_dislotwin_homogenizedC + state%p(3_pInt*ns+i)*lattice_C66(1:6,1:6,phase)
enddo
end function constitutive_dislotwin_homogenizedC
@ -992,7 +992,7 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
el !< element
real(pReal), intent(in) :: &
temperature !< temperature at IP
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
type(p_vec), intent(inout) :: &
state !< microstructure state
integer(pInt) :: &
@ -1022,7 +1022,7 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
!* State: 7*ns+5*nt+1 : 7*ns+6*nt twin volume
!* Total twin volume fraction
sumf = sum(state(ipc,ip,el)%p((3*ns+1):(3*ns+nt))) ! safe for nt == 0
sumf = sum(state%p((3*ns+1):(3*ns+nt))) ! safe for nt == 0
!* Stacking fault energy
sfe = constitutive_dislotwin_SFE_0K(instance) + &
@ -1031,61 +1031,61 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
!* rescaled twin volume fraction for topology
forall (t = 1_pInt:nt) &
fOverStacksize(t) = &
state(ipc,ip,el)%p(3_pInt*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,instance)
state%p(3_pInt*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,instance)
!* 1/mean free distance between 2 forest dislocations seen by a moving dislocation
forall (s = 1_pInt:ns) &
state(ipc,ip,el)%p(3_pInt*ns+2_pInt*nt+s) = &
sqrt(dot_product((state(ipc,ip,el)%p(1:ns)+state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
state%p(3_pInt*ns+2_pInt*nt+s) = &
sqrt(dot_product((state%p(1:ns)+state%p(ns+1_pInt:2_pInt*ns)),&
constitutive_dislotwin_forestProjectionEdge(1:ns,s,instance)))/ &
constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,instance)
!* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
!$OMP CRITICAL (evilmatmul)
state(ipc,ip,el)%p((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt)) = 0.0_pReal
state%p((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt)) = 0.0_pReal
if (nt > 0_pInt .and. ns > 0_pInt) &
state(ipc,ip,el)%p((4_pInt*ns+2_pInt*nt+1):(5_pInt*ns+2_pInt*nt)) = &
state%p((4_pInt*ns+2_pInt*nt+1):(5_pInt*ns+2_pInt*nt)) = &
matmul(constitutive_dislotwin_interactionMatrix_SlipTwin(1:ns,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf)
!$OMP END CRITICAL (evilmatmul)
!* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
!$OMP CRITICAL (evilmatmul)
if (nt > 0_pInt) &
state(ipc,ip,el)%p((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt)) = &
state%p((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt)) = &
matmul(constitutive_dislotwin_interactionMatrix_TwinTwin(1:nt,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf)
!$OMP END CRITICAL (evilmatmul)
!* mean free path between 2 obstacles seen by a moving dislocation
do s = 1_pInt,ns
if (nt > 0_pInt) then
state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+s) = &
state%p(5_pInt*ns+3_pInt*nt+s) = &
constitutive_dislotwin_GrainSize(instance)/(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*&
(state(ipc,ip,el)%p(3_pInt*ns+2_pInt*nt+s)+state(ipc,ip,el)%p(4_pInt*ns+2_pInt*nt+s)))
(state%p(3_pInt*ns+2_pInt*nt+s)+state%p(4_pInt*ns+2_pInt*nt+s)))
else
state(ipc,ip,el)%p(5_pInt*ns+s) = &
state%p(5_pInt*ns+s) = &
constitutive_dislotwin_GrainSize(instance)/&
(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*(state(ipc,ip,el)%p(3_pInt*ns+s)))
(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*(state%p(3_pInt*ns+s)))
endif
enddo
!* mean free path between 2 obstacles seen by a growing twin
forall (t = 1_pInt:nt) &
state(ipc,ip,el)%p(6_pInt*ns+3_pInt*nt+t) = &
state%p(6_pInt*ns+3_pInt*nt+t) = &
(constitutive_dislotwin_Cmfptwin(instance)*constitutive_dislotwin_GrainSize(instance))/&
(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*state(ipc,ip,el)%p(5_pInt*ns+2_pInt*nt+t))
(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*state%p(5_pInt*ns+2_pInt*nt+t))
!* threshold stress for dislocation motion
if(lattice_structure(phase) /= LATTICE_BCC_ID) then ! bcc value remains constant
forall (s = 1_pInt:ns) &
state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(instance)+ &
state%p(6_pInt*ns+4_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(instance)+ &
lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(s,instance)*&
sqrt(dot_product((state(ipc,ip,el)%p(1:ns)+state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
sqrt(dot_product((state%p(1:ns)+state%p(ns+1_pInt:2_pInt*ns)),&
constitutive_dislotwin_interactionMatrix_SlipSlip(s,1:ns,instance)))
endif
!* threshold stress for growing twin
forall (t = 1_pInt:nt) &
state(ipc,ip,el)%p(7_pInt*ns+4_pInt*nt+t) = &
state%p(7_pInt*ns+4_pInt*nt+t) = &
constitutive_dislotwin_Cthresholdtwin(instance)*&
(sfe/(3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,instance))+&
3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,instance)*lattice_mu(phase)/&
@ -1093,8 +1093,8 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
!* final twin volume after growth
forall (t = 1_pInt:nt) &
state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+t) = &
(pi/4.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,instance)*state(ipc,ip,el)%p(6*ns+3*nt+t)**(2.0_pReal)
state%p(7_pInt*ns+5_pInt*nt+t) = &
(pi/4.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,instance)*state%p(6*ns+3*nt+t)**(2.0_pReal)
!* equilibrium seperation of partial dislocations
do t = 1_pInt,nt
@ -1145,12 +1145,12 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
LATTICE_fcc_ID
implicit none
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: state
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), intent(inout) :: state
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar
integer(pInt) :: instance,phase,ns,nt,f,i,j,k,l,m,n,index_myFamily,s1,s2
real(pReal) :: sumf,StressRatio_p,StressRatio_pminus1,StressRatio_r,BoltzmannRatio,DotGamma0,Ndot0
@ -1190,7 +1190,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
nt = constitutive_dislotwin_totalNtwin(instance)
!* Total twin volume fraction
sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
@ -1214,16 +1214,16 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))&
StressRatio_p = (abs(tau_slip(j))/state%p(6*ns+4*nt+j))&
**constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))&
StressRatio_pminus1 = (abs(tau_slip(j))/state%p(6*ns+4*nt+j))&
**(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates
DotGamma0 = &
state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*&
state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*&
constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
@ -1234,7 +1234,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
!* Derivatives of shear rates
dgdot_dtauslip(j) = &
((abs(gdot_slip(j))*BoltzmannRatio*constitutive_dislotwin_pPerSlipFamily(f,instance)&
*constitutive_dislotwin_qPerSlipFamily(f,instance))/state(ipc,ip,el)%p(6*ns+4*nt+j))*&
*constitutive_dislotwin_qPerSlipFamily(f,instance))/state%p(6*ns+4*nt+j))*&
StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal)
!* Plastic velocity gradient for dislocation glide
@ -1320,15 +1320,15 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
!* Stress ratios
if (tau_twin(j) > tol_math_check) then
StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance)
StressRatio_r = (state%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance)
!* Shear rates and their derivatives due to twin
select case(lattice_structure(phase))
case (LATTICE_fcc_ID)
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i)
if (tau_twin(j) < constitutive_dislotwin_tau_r(j,instance)) then
Ndot0=(abs(gdot_slip(s1))*(state(ipc,ip,el)%p(s2)+state(ipc,ip,el)%p(ns+s2))+&
abs(gdot_slip(s2))*(state(ipc,ip,el)%p(s1)+state(ipc,ip,el)%p(ns+s1)))/&
Ndot0=(abs(gdot_slip(s1))*(state%p(s2)+state%p(ns+s2))+&
abs(gdot_slip(s2))*(state%p(s1)+state%p(ns+s1)))/&
(constitutive_dislotwin_L0(instance)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))*&
(1.0_pReal-exp(-constitutive_dislotwin_VcrossSlip(instance)/(kB*Temperature)*&
(constitutive_dislotwin_tau_r(j,instance)-tau_twin(j))))
@ -1340,7 +1340,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
end select
gdot_twin(j) = &
(constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*&
state(ipc,ip,el)%p(7*ns+5*nt+j)*Ndot0*exp(-StressRatio_r)
state%p(7*ns+5*nt+j)*Ndot0*exp(-StressRatio_r)
dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_rPerTwinFamily(f,instance))/tau_twin(j))*StressRatio_r
endif
@ -1392,15 +1392,15 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
LATTICE_bcc_ID
implicit none
real(pReal), dimension(6), intent(in):: &
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_dislotwin_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_dislotwin_dotState
@ -1421,7 +1421,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
nt = constitutive_dislotwin_totalNtwin(instance)
!* Total twin volume fraction
sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
constitutive_dislotwin_dotState = 0.0_pReal
@ -1441,16 +1441,16 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
StressRatio_p = (abs(tau_slip(j))/state%p(6_pInt*ns+4_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
StressRatio_pminus1 = (abs(tau_slip(j))/state%p(6_pInt*ns+4_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates
DotGamma0 = &
state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*&
state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*&
constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
@ -1459,7 +1459,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
!* Multiplication
DotRhoMultiplication(j) = abs(gdot_slip(j))/&
(constitutive_dislotwin_burgersPerSlipSystem(j,instance)*state(ipc,ip,el)%p(5*ns+3*nt+j))
(constitutive_dislotwin_burgersPerSlipSystem(j,instance)*state%p(5*ns+3*nt+j))
!* Dipole formation
EdgeDipMinDistance = &
@ -1470,22 +1470,22 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
EdgeDipDistance(j) = &
(3.0_pReal*lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))/&
(16.0_pReal*pi*abs(tau_slip(j)))
if (EdgeDipDistance(j)>state(ipc,ip,el)%p(5*ns+3*nt+j)) EdgeDipDistance(j)=state(ipc,ip,el)%p(5*ns+3*nt+j)
if (EdgeDipDistance(j)>state%p(5*ns+3*nt+j)) EdgeDipDistance(j)=state%p(5*ns+3*nt+j)
if (EdgeDipDistance(j)<EdgeDipMinDistance) EdgeDipDistance(j)=EdgeDipMinDistance
DotRhoDipFormation(j) = &
((2.0_pReal*EdgeDipDistance(j))/constitutive_dislotwin_burgersPerSlipSystem(j,instance))*&
state(ipc,ip,el)%p(j)*abs(gdot_slip(j))
state%p(j)*abs(gdot_slip(j))
endif
!* Spontaneous annihilation of 2 single edge dislocations
DotRhoEdgeEdgeAnnihilation(j) = &
((2.0_pReal*EdgeDipMinDistance)/constitutive_dislotwin_burgersPerSlipSystem(j,instance))*&
state(ipc,ip,el)%p(j)*abs(gdot_slip(j))
state%p(j)*abs(gdot_slip(j))
!* Spontaneous annihilation of a single edge dislocation with a dipole constituent
DotRhoEdgeDipAnnihilation(j) = &
((2.0_pReal*EdgeDipMinDistance)/constitutive_dislotwin_burgersPerSlipSystem(j,instance))*&
state(ipc,ip,el)%p(ns+j)*abs(gdot_slip(j))
state%p(ns+j)*abs(gdot_slip(j))
!* Dislocation dipole climb
AtomicVolume = &
@ -1499,7 +1499,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
((3.0_pReal*lattice_mu(phase)*VacancyDiffusion*AtomicVolume)/(2.0_pReal*pi*kB*Temperature))*&
(1/(EdgeDipDistance(j)+EdgeDipMinDistance))
DotRhoEdgeDipClimb(j) = &
(4.0_pReal*ClimbVelocity(j)*state(ipc,ip,el)%p(ns+j))/(EdgeDipDistance(j)-EdgeDipMinDistance)
(4.0_pReal*ClimbVelocity(j)*state%p(ns+j))/(EdgeDipDistance(j)-EdgeDipMinDistance)
endif
!* Edge dislocation density rate of change
@ -1527,7 +1527,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase))
!* Stress ratios
if (tau_twin(j) > tol_math_check) then
StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance)
StressRatio_r = (state%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance)
!* Shear rates and their derivatives due to twin
select case(lattice_structure(phase))
@ -1535,8 +1535,8 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i)
if (tau_twin(j) < constitutive_dislotwin_tau_r(j,instance)) then
Ndot0=(abs(gdot_slip(s1))*(state(ipc,ip,el)%p(s2)+state(ipc,ip,el)%p(ns+s2))+&
abs(gdot_slip(s2))*(state(ipc,ip,el)%p(s1)+state(ipc,ip,el)%p(ns+s1)))/&
Ndot0=(abs(gdot_slip(s1))*(state%p(s2)+state%p(ns+s2))+&
abs(gdot_slip(s2))*(state%p(s1)+state%p(ns+s1)))/&
(constitutive_dislotwin_L0(instance)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))*&
(1.0_pReal-exp(-constitutive_dislotwin_VcrossSlip(instance)/(kB*Temperature)*&
(constitutive_dislotwin_tau_r(j,instance)-tau_twin(j))))
@ -1548,7 +1548,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
end select
constitutive_dislotwin_dotState(3_pInt*ns+j) = &
(constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*&
state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r)
state%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r)
!* Dotstate for accumulated shear due to twin
constitutive_dislotwin_dotState(3_pInt*ns+nt+j) = constitutive_dislotwin_dotState(3_pInt*ns+j) * &
lattice_sheartwin(index_myfamily+i,phase)
@ -1592,15 +1592,15 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
LATTICE_fcc_ID
implicit none
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_dislotwin_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -1625,7 +1625,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
nt = constitutive_dislotwin_totalNtwin(instance)
!* Total twin volume fraction
sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
!* Required output
c = 0_pInt
@ -1638,10 +1638,10 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
select case(constitutive_dislotwin_outputID(o,instance))
case (edge_density_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(1_pInt:ns)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state%p(1_pInt:ns)
c = c + ns
case (dipole_density_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state%p(ns+1_pInt:2_pInt*ns)
c = c + ns
case (shear_rate_slip_ID)
j = 0_pInt
@ -1657,16 +1657,16 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
StressRatio_p = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
StressRatio_pminus1 = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates
DotGamma0 = &
state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
@ -1677,11 +1677,11 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
c = c + ns
case (accumulated_shear_slip_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = &
state(ipc,ip,el)%p((2_pInt*ns+1_pInt):(3_pInt*ns))
state%p((2_pInt*ns+1_pInt):(3_pInt*ns))
c = c + ns
case (mfp_slip_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) =&
state(ipc,ip,el)%p((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt))
state%p((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt))
c = c + ns
case (resolved_stress_slip_ID)
j = 0_pInt
@ -1695,7 +1695,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
c = c + ns
case (threshold_stress_slip_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = &
state(ipc,ip,el)%p((6_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+4_pInt*nt))
state%p((6_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+4_pInt*nt))
c = c + ns
case (edge_dipole_distance_ID)
j = 0_pInt
@ -1706,8 +1706,8 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
constitutive_dislotwin_postResults(c+j) = &
(3.0_pReal*lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))/&
(16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))))
constitutive_dislotwin_postResults(c+j)=min(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(5*ns+3*nt+j))
! constitutive_dislotwin_postResults(c+j)=max(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(4*ns+2*nt+j))
constitutive_dislotwin_postResults(c+j)=min(constitutive_dislotwin_postResults(c+j),state%p(5*ns+3*nt+j))
! constitutive_dislotwin_postResults(c+j)=max(constitutive_dislotwin_postResults(c+j),state%p(4*ns+2*nt+j))
enddo; enddo
c = c + ns
case (resolved_stress_shearband_ID)
@ -1740,7 +1740,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
enddo
c = c + 6_pInt
case (twin_fraction_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))
c = c + nt
case (shear_rate_twin_ID)
if (nt > 0_pInt) then
@ -1758,16 +1758,16 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
StressRatio_p = (abs(tau)/state%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
StressRatio_pminus1 = (abs(tau)/state%p(5_pInt*ns+3_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates
DotGamma0 = &
state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
@ -1784,7 +1784,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Resolved shear stress on twin system
tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase))
!* Stress ratios
StressRatio_r = (state(ipc,ip,el)%p(7_pInt*ns+4_pInt*nt+j)/tau)**constitutive_dislotwin_rPerTwinFamily(f,instance)
StressRatio_r = (state%p(7_pInt*ns+4_pInt*nt+j)/tau)**constitutive_dislotwin_rPerTwinFamily(f,instance)
!* Shear rates due to twin
if ( tau > 0.0_pReal ) then
@ -1793,8 +1793,8 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i)
if (tau < constitutive_dislotwin_tau_r(j,instance)) then
Ndot0=(abs(gdot_slip(s1))*(state(ipc,ip,el)%p(s2)+state(ipc,ip,el)%p(ns+s2))+&
abs(gdot_slip(s2))*(state(ipc,ip,el)%p(s1)+state(ipc,ip,el)%p(ns+s1)))/&
Ndot0=(abs(gdot_slip(s1))*(state%p(s2)+state%p(ns+s2))+&
abs(gdot_slip(s2))*(state%p(s1)+state%p(ns+s1)))/&
(constitutive_dislotwin_L0(instance)*&
constitutive_dislotwin_burgersPerSlipSystem(j,instance))*&
(1.0_pReal-exp(-constitutive_dislotwin_VcrossSlip(instance)/(kB*Temperature)*&
@ -1807,17 +1807,17 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
end select
constitutive_dislotwin_postResults(c+j) = &
(constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*&
state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r)
state%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r)
endif
enddo ; enddo
endif
c = c + nt
case (accumulated_shear_twin_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((3_pInt*ns+nt+1_pInt):(3_pInt*ns+2_pInt*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((3_pInt*ns+nt+1_pInt):(3_pInt*ns+2_pInt*nt))
c = c + nt
case (mfp_twin_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt))
c = c + nt
case (resolved_stress_twin_ID)
if (nt > 0_pInt) then
@ -1831,7 +1831,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
endif
c = c + nt
case (threshold_stress_twin_ID)
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((7_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+5_pInt*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((7_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+5_pInt*nt))
c = c + nt
case (stress_exponent_ID)
j = 0_pInt
@ -1846,16 +1846,16 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
StressRatio_p = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
StressRatio_pminus1 = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates
DotGamma0 = &
state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
@ -1866,7 +1866,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
dgdot_dtauslip = &
((abs(gdot_slip(j))*BoltzmannRatio*&
constitutive_dislotwin_pPerSlipFamily(f,instance)*constitutive_dislotwin_qPerSlipFamily(f,instance))/&
state(ipc,ip,el)%p(6*ns+4*nt+j))*StressRatio_pminus1*(1_pInt-StressRatio_p)**&
state%p(6*ns+4*nt+j))*StressRatio_pminus1*(1_pInt-StressRatio_p)**&
(constitutive_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal)
!* Stress exponent

34
code/constitutive_j2.f90
View File

@ -319,18 +319,18 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(3,3) :: &
@ -355,7 +355,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = constitutive_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / constitutive_j2_fTaylor(instance) / state(ipc,ip,el)%p(1)) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / constitutive_j2_fTaylor(instance) / state%p(1)) &
**constitutive_j2_n(instance)
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(instance)
@ -395,13 +395,13 @@ pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el)
implicit none
real(pReal), dimension(1) :: &
constitutive_j2_dotState
real(pReal), dimension(6), intent(in):: &
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(6) :: &
@ -425,7 +425,7 @@ pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el)
! strain rate
gamma_dot = constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!-----------------------------------------------------------------------------------
(constitutive_j2_fTaylor(instance) * state(ipc,ip,el)%p(1)) ) ** constitutive_j2_n(instance)
(constitutive_j2_fTaylor(instance) * state%p(1)) ) ** constitutive_j2_n(instance)
!--------------------------------------------------------------------------------------------------
! hardening coefficient
@ -447,8 +447,8 @@ pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el)
)
endif
hardening = ( constitutive_j2_h0(instance) + constitutive_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
* abs( 1.0_pReal - state(ipc,ip,el)%p(1)/saturation )**constitutive_j2_a(instance) &
* sign(1.0_pReal, 1.0_pReal - state(ipc,ip,el)%p(1)/saturation)
* abs( 1.0_pReal - state%p(1)/saturation )**constitutive_j2_a(instance) &
* sign(1.0_pReal, 1.0_pReal - state%p(1)/saturation)
else
hardening = 0.0_pReal
endif
@ -476,13 +476,13 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
phase_Noutput
implicit none
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -511,13 +511,13 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_j2_outputID(o,instance))
case (flowstress_ID)
constitutive_j2_postResults(c+1_pInt) = state(ipc,ip,el)%p(1)
constitutive_j2_postResults(c+1_pInt) = state%p(1)
c = c + 1_pInt
case (strainrate_ID)
constitutive_j2_postResults(c+1_pInt) = &
constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!----------------------------------------------------------------------------------
(constitutive_j2_fTaylor(instance) * state(ipc,ip,el)%p(1)) ) ** constitutive_j2_n(instance)
(constitutive_j2_fTaylor(instance) * state%p(1)) ) ** constitutive_j2_n(instance)
c = c + 1_pInt
end select
enddo outputsLoop

72
code/constitutive_nonlocal.f90
View File

@ -1297,7 +1297,7 @@ end function constitutive_nonlocal_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, gr, ip, el)
subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, ipc, ip, el)
use IO, only: &
IO_error
@ -1347,7 +1347,7 @@ use lattice, only: &
implicit none
!*** input variables
integer(pInt), intent(in) :: gr, & ! current grain ID
integer(pInt), intent(in) :: ipc, & ! current grain ID
ip, & ! current integration point
el ! current element
real(pReal), dimension(3,3), intent(in) :: &
@ -1387,7 +1387,7 @@ real(pReal), dimension(2) :: rhoExcessGradient, &
rhoTotal
real(pReal), dimension(3) :: rhoExcessDifferences, &
normal_latticeConf
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
rhoForest, & ! forest dislocation density
tauBack, & ! back stress from pileup on same slip system
tauThreshold ! threshold shear stress
@ -1397,24 +1397,24 @@ real(pReal), dimension(3,3) :: invFe, & ! inverse of elast
invConnections
real(pReal), dimension(3,mesh_maxNipNeighbors) :: &
connection_latticeConf
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
rhoExcess
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: &
rhoDip ! dipole dislocation density (edge, screw)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),8) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: &
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))), &
totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))), &
totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: &
neighbor_rhoExcess, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: &
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: &
m ! direction of dislocation motion
logical inversionError
phase = material_phase(gr,ip,el)
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase)
ns = totalNslip(instance)
@ -1422,11 +1422,11 @@ ns = totalNslip(instance)
!*** get basic states
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(gr,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(gr,ip,el)%p(iRhoB(s,t,instance))
rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance))
endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
rhoDip(s,c) = max(state(gr,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities
rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) &
.or. abs(rhoSgl) < significantRho(instance)) &
rhoSgl = 0.0_pReal
@ -1488,7 +1488,7 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
neighbor_phase = material_phase(gr,neighbor_ip,neighbor_el)
neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el)
neighbor_instance = phase_plasticityInstance(neighbor_phase)
neighbor_ns = totalNslip(neighbor_instance)
if (.not. phase_localPlasticity(neighbor_phase) &
@ -1497,14 +1497,14 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
nRealNeighbors = nRealNeighbors + 1_pInt
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
neighbor_rhoExcess(c,s,n) = &
max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles
- max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) ! negative mobiles
max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles
- max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) ! negative mobiles
neighbor_rhoTotal(c,s,n) = &
max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles
+ max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) & ! negative mobiles
+ abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads
+ abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) & ! negative deads
+ max(state(gr,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_instance)), 0.0_pReal) ! dipoles
max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles
+ max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) & ! negative mobiles
+ abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads
+ abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) & ! negative deads
+ max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_instance)), 0.0_pReal) ! dipoles
endforall
connection_latticeConf(1:3,n) = &
math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
@ -1586,17 +1586,17 @@ endif
!*** set dependent states
state(gr,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest
state(gr,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold
state(gr,ip,el)%p(iTauB(1:ns,instance)) = tauBack
state(ipc,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest
state(ipc,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold
state(ipc,ip,el)%p(iTauB(1:ns,instance)) = tauBack
#ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == gr)&
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,gr
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,ipc
write(6,*)
write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest
write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
@ -2004,7 +2004,7 @@ integer(pInt), intent(in) :: ipc, & ! current
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
!*** input/output variables
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
type(p_vec), intent(inout) :: &
state ! current microstructural state
!*** output variables
@ -2056,15 +2056,15 @@ ns = totalNslip(instance)
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance))
v(s,t) = state(ipc,ip,el)%p(iV(s,t,instance))
rhoSgl(s,t) = max(state%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state%p(iRhoB(s,t,instance))
v(s,t) = state%p(iV(s,t,instance))
endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities
dUpperOld(s,c) = state(ipc,ip,el)%p(iD(s,c,instance))
rhoDip(s,c) = max(state%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities
dUpperOld(s,c) = state%p(iD(s,c,instance))
endforall
tauBack = state(ipc,ip,el)%p(iTauB(1:ns,instance))
tauBack = state%p(iTauB(1:ns,instance))
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) &
.or. abs(rhoSgl) < significantRho(instance)) &
@ -2168,7 +2168,7 @@ forall (t=1_pInt:4_pInt) &
!*** store new maximum dipole height in state
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
state(ipc,ip,el)%p(iD(s,c,instance)) = dUpper(s,c)
state%p(iD(s,c,instance)) = dUpper(s,c)
@ -3362,14 +3362,14 @@ pure function constitutive_nonlocal_postResults(Tstar_v,Fe,state,dotState,ipc,ip
lattice_nu
implicit none
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe !< elastic deformation gradient
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element

78
code/constitutive_phenopowerlaw.f90
View File

@ -630,18 +630,18 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
integer(pInt) :: &
@ -693,12 +693,12 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,phase)
enddo
gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
((abs(tau_slip_pos(j))/state%p(j))**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))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
((abs(tau_slip_neg(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_neg(j))
Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
Lp = Lp + (1.0_pReal-state%p(index_F))*& ! 1-F
(gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase)
!--------------------------------------------------------------------------------------------------
@ -723,16 +723,16 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
j = 0_pInt
twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! 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
!--------------------------------------------------------------------------------------------------
! Calculation of Lp
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
gdot_twin(j) = (1.0_pReal-state%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
(abs(tau_twin(j))/state%p(nSlip+j))**&
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,phase)
@ -777,13 +777,13 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -820,17 +820,17 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(instance)*&
(1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*state(ipc,ip,el)%p(index_F)**&
(1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*state%p(index_F)**&
constitutive_phenopowerlaw_twinB(instance))
c_SlipTwin = 0.0_pReal
c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(instance)*&
state(ipc,ip,el)%p(index_Gamma)**constitutive_phenopowerlaw_twinE(instance)
state%p(index_Gamma)**constitutive_phenopowerlaw_twinE(instance)
c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(instance)*&
state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(instance)
state%p(index_F)**constitutive_phenopowerlaw_twinD(instance)
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors and calculate dot gammas
ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state(ipc,ip,el)%p(index_F))
ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state%p(index_F))
j = 0_pInt
slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
@ -838,10 +838,10 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
j = j+1_pInt
left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part
left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part
right_SlipSlip(j) = abs(1.0_pReal-state(ipc,ip,el)%p(j) / &
right_SlipSlip(j) = abs(1.0_pReal-state%p(j) / &
(constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) &
**constitutive_phenopowerlaw_a_slip(instance)&
*sign(1.0_pReal,1.0_pReal-state(ipc,ip,el)%p(j) / &
*sign(1.0_pReal,1.0_pReal-state%p(j) / &
(constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset))
right_TwinSlip(j) = 1.0_pReal ! no system-dependent part
@ -856,8 +856,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase))
enddo
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))&
((abs(tau_slip_pos(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos(j))
enddo
enddo slipFamiliesLoop1
@ -875,9 +875,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
gdot_twin(j) = (1.0_pReal-state%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
(abs(tau_twin(j))/state%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
enddo
enddo twinFamiliesLoop1
@ -913,7 +913,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
c_TwinTwin * left_TwinTwin(j) * &
dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(j,1:nTwin,instance), &
right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor
if (state(ipc,ip,el)%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0
if (state%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0
constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + &
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,phase)
constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j))
@ -950,13 +950,13 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
mesh_maxNips
implicit none
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -987,11 +987,11 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_phenopowerlaw_outputID(o,instance))
case (resistance_slip_ID)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(1:nSlip)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state%p(1:nSlip)
c = c + nSlip
case (accumulatedshear_slip_ID)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(index_accshear_slip:&
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state%p(index_accshear_slip:&
index_accshear_slip+nSlip)
c = c + nSlip
@ -1010,8 +1010,8 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase))
enddo
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))&
((abs(tau_slip_pos)/state%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos)
enddo
enddo slipFamiliesLoop1
@ -1031,17 +1031,17 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (totalshear_ID)
constitutive_phenopowerlaw_postResults(c+1_pInt) = &
state(ipc,ip,el)%p(index_Gamma)
state%p(index_Gamma)
c = c + 1_pInt
case (resistance_twin_ID)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state(ipc,ip,el)%p(1_pInt+nSlip:nTwin+nSlip)
state%p(1_pInt+nSlip:nTwin+nSlip)
c = c + nTwin
case (accumulatedshear_twin_ID)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state(ipc,ip,el)%p(index_accshear_twin:index_accshear_twin+nTwin)
state%p(index_accshear_twin:index_accshear_twin+nTwin)
c = c + nTwin
case (shearrate_twin_ID)
@ -1051,9 +1051,9 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau)/state(ipc,ip,el)%p(j+nSlip))**&
(abs(tau)/state%p(j+nSlip))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau))
enddo
enddo twinFamiliesLoop1
@ -1072,7 +1072,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
c = c + nTwin
case (totalvolfrac_ID)
constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_F)
constitutive_phenopowerlaw_postResults(c+1_pInt) = state%p(index_F)
c = c + 1_pInt
end select

181
code/constitutive_titanmod.f90
View File

@ -1108,7 +1108,7 @@ implicit none
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
volumefraction_PerTwinSys
@ -1130,7 +1130,7 @@ real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance
!--------------------------------------------------------------------------------------------------
! total twin volume fraction
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
sumf = sum(abs(volumefraction_PerTwinSys(1:nt))) ! safe for nt == 0
@ -1170,7 +1170,7 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
el !< element
real(pReal), intent(in) :: &
temperature !< temperature at IP
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
type(p_vec), intent(inout) :: &
state !< microstructure state
integer(pInt) :: &
@ -1192,13 +1192,10 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! total twin volume fraction
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ &
forall (i = 1_pInt:nt) &
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
sumf = sum(abs(volumefraction_PerTwinSys(1:nt))) ! safe for nt == 0
sfe = 0.0002_pReal*Temperature-0.0396_pReal
@ -1206,43 +1203,43 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! average segment length for edge dislocations in matrix
forall (s = 1_pInt:ns) &
state(ipc,ip,el)%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(state(ipc,ip,el)%p(1:ns), &
state%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(state%p(1:ns), &
constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ &
dot_product(state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns), &
dot_product(state%p(ns+1_pInt:2_pInt*ns), &
constitutive_titanmod_forestProjectionScrew(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! average segment length for screw dislocations in matrix
forall (s = 1_pInt:ns) &
state(ipc,ip,el)%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(state(ipc,ip,el)%p(1:ns), &
state%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(state%p(1:ns), &
constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ &
dot_product(state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns), &
dot_product(state%p(ns+1_pInt:2_pInt*ns), &
constitutive_titanmod_forestProjectionScrew(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! threshold stress or slip resistance for edge dislocation motion
forall (s = 1_pInt:ns) &
state(ipc,ip,el)%p(5_pInt*ns+nt+s) = &
state%p(5_pInt*ns+nt+s) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)*&
sqrt(dot_product((state(ipc,ip,el)%p(1:ns)),&
sqrt(dot_product((state%p(1:ns)),&
constitutive_titanmod_interactionMatrix_ee(1:ns,s,instance))+ &
dot_product((state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
dot_product((state%p(ns+1_pInt:2_pInt*ns)),&
constitutive_titanmod_interactionMatrix_es(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! threshold stress or slip resistance for screw dislocation motion
forall (s = 1_pInt:ns) &
state(ipc,ip,el)%p(6_pInt*ns+nt+s) = &
state%p(6_pInt*ns+nt+s) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)*&
sqrt(dot_product((state(ipc,ip,el)%p(1:ns)),&
sqrt(dot_product((state%p(1:ns)),&
constitutive_titanmod_interactionMatrix_es(1:ns,s,instance))+ &
dot_product((state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
dot_product((state%p(ns+1_pInt:2_pInt*ns)),&
constitutive_titanmod_interactionMatrix_ss(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! threshold stress or slip resistance for dislocation motion in twin
forall (t = 1_pInt:nt) &
state(ipc,ip,el)%p(7_pInt*ns+nt+t) = &
state%p(7_pInt*ns+nt+t) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerTwinSys(t,instance)*&
(dot_product((abs(state(ipc,ip,el)%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),&
(dot_product((abs(state%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),&
constitutive_titanmod_interactionMatrixTwinTwin(1:nt,t,instance)))
end subroutine constitutive_titanmod_microstructure
@ -1278,20 +1275,20 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
temperature !< temperature at IP
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
type(p_vec), intent(inout) :: &
state !< microstructure state
integer(pInt) :: &
index_myFamily, instance,phase, &
@ -1318,7 +1315,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
nt = constitutive_titanmod_totalNtwin(instance)
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
@ -1346,7 +1343,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
if(lattice_structure(phase)==LATTICE_hex_ID) then ! only for prismatic and pyr <a> systems in hex
screwvelocity_prefactor=constitutive_titanmod_debyefrequency(instance)* &
state(ipc,ip,el)%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSys(j,instance)/ &
state%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSys(j,instance)/ &
constitutive_titanmod_kinkcriticallength_PerSlipSys(j,instance))**2
!* Stress ratio for screw ! No slip resistance for screw dislocations, only Peierls stress
@ -1371,7 +1368,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
else ! if the structure is not hex or the slip family is basal
screwvelocity_prefactor=constitutive_titanmod_v0s_PerSlipSys(j,instance)
bottomstress_screw=constitutive_titanmod_tau0s_PerSlipSys(j,instance)+state(ipc,ip,el)%p(6*ns+nt+j)
bottomstress_screw=constitutive_titanmod_tau0s_PerSlipSys(j,instance)+state%p(6*ns+nt+j)
StressRatio_screw_p = ((abs(tau_slip(j)))/( bottomstress_screw ))**constitutive_titanmod_ps_PerSlipSys(j,instance)
if((1.0_pReal-StressRatio_screw_p)>0.001_pReal) then
@ -1389,7 +1386,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
endif
!* Stress ratio for edge
bottomstress_edge=constitutive_titanmod_tau0e_PerSlipSys(j,instance)+state(ipc,ip,el)%p(5*ns+nt+j)
bottomstress_edge=constitutive_titanmod_tau0e_PerSlipSys(j,instance)+state%p(5*ns+nt+j)
StressRatio_edge_p = ((abs(tau_slip(j)))/ &
( bottomstress_edge) &
)**constitutive_titanmod_pe_PerSlipSys(j,instance)
@ -1415,29 +1412,29 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
constitutive_titanmod_qe_PerSlipSys(j,instance))
!* Shear rates due to edge slip
gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state(ipc,ip,el)%p(j)* &
gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state%p(j)* &
edge_velocity(j))* sign(1.0_pReal,tau_slip(j))
!* Shear rates due to screw slip
gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state(ipc,ip,el)%p(ns+j) * &
gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state%p(ns+j) * &
screw_velocity(j))* sign(1.0_pReal,tau_slip(j))
!Total shear rate
gdot_slip(j) = gdot_slip_edge(j) + gdot_slip_screw(j)
state(ipc,ip,el)%p(7*ns+2*nt+j)=edge_velocity(j)
state(ipc,ip,el)%p(8*ns+2*nt+j)=screw_velocity(j)
state(ipc,ip,el)%p(9*ns+2*nt+j)=tau_slip(j)
state(ipc,ip,el)%p(10*ns+2*nt+j)=gdot_slip_edge(j)
state(ipc,ip,el)%p(11*ns+2*nt+j)=gdot_slip_screw(j)
state(ipc,ip,el)%p(12*ns+2*nt+j)=StressRatio_edge_p
state(ipc,ip,el)%p(13*ns+2*nt+j)=StressRatio_screw_p
state%p(7*ns+2*nt+j)=edge_velocity(j)
state%p(8*ns+2*nt+j)=screw_velocity(j)
state%p(9*ns+2*nt+j)=tau_slip(j)
state%p(10*ns+2*nt+j)=gdot_slip_edge(j)
state%p(11*ns+2*nt+j)=gdot_slip_screw(j)
state%p(12*ns+2*nt+j)=StressRatio_edge_p
state%p(13*ns+2*nt+j)=StressRatio_screw_p
!* Derivatives of shear rates
dgdot_dtauslip(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(( &
( &
( &
( &
(edge_velocity(j)*state(ipc,ip,el)%p(j))) * &
(edge_velocity(j)*state%p(j))) * &
BoltzmannRatioedge*&
constitutive_titanmod_pe_PerSlipSys(j,instance)* &
constitutive_titanmod_qe_PerSlipSys(j,instance) &
@ -1450,7 +1447,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
( &
( &
( &
(state(ipc,ip,el)%p(ns+j) * screw_velocity(j)) * &
(state%p(ns+j) * screw_velocity(j)) * &
BoltzmannRatioscrew* &
constitutive_titanmod_ps_PerSlipSys(j,instance)* &
constitutive_titanmod_qs_PerSlipSys(j,instance) &
@ -1492,20 +1489,20 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
!**************************************************************************************
!* Stress ratios
! StressRatio_r = (state(ipc,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_titanmod_r(instance)
! StressRatio_r = (state%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_titanmod_r(instance)
!* Shear rates and their derivatives due to twin
! if ( tau_twin(j) > 0.0_pReal ) !then
! gdot_twin(j) = 0.0_pReal!&
! (constitutive_titanmod_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*&
! state(ipc,ip,el)%p(6*ns+4*nt+j)*constitutive_titanmod_Ndot0PerTwinSys(f,instance)*exp(-StressRatio_r)
! state%p(6*ns+4*nt+j)*constitutive_titanmod_Ndot0PerTwinSys(f,instance)*exp(-StressRatio_r)
! dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_titanmod_r(instance))/tau_twin(j))*StressRatio_r
! endif
!**************************************************************************************
!* Stress ratio for edge
twinStressRatio_p = ((abs(tau_twin(j)))/ &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state(ipc,ip,el)%p(7*ns+nt+j)) &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) &
)**constitutive_titanmod_twinp_PerTwinSys(j,instance)
if((1.0_pReal-twinStressRatio_p)>0.001_pReal) then
@ -1515,7 +1512,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
endif
twinStressRatio_pminus1 = ((abs(tau_twin(j)))/ &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state(ipc,ip,el)%p(7*ns+nt+j)) &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) &
)**(constitutive_titanmod_twinp_PerTwinSys(j,instance)-1.0_pReal)
!* Boltzmann ratio
@ -1584,15 +1581,15 @@ function constitutive_titanmod_dotState(Tstar_v,temperature,state,ipc,ip,el)
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in):: &
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_titanmod_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_titanmod_dotState
@ -1622,7 +1619,7 @@ implicit none
nt = constitutive_titanmod_totalNtwin(instance)
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
@ -1638,13 +1635,13 @@ implicit none
j = j+1_pInt
DotRhoEdgeGeneration(j) = & ! multiplication of edge dislocations
state(ipc,ip,el)%p(ns+j)*state(ipc,ip,el)%p(8*ns+2*nt+j)/state(ipc,ip,el)%p(4*ns+nt+j)
state%p(ns+j)*state%p(8*ns+2*nt+j)/state%p(4*ns+nt+j)
DotRhoScrewGeneration(j) = & ! multiplication of screw dislocations
state(ipc,ip,el)%p(j)*state(ipc,ip,el)%p(7*ns+2*nt+j)/state(ipc,ip,el)%p(3*ns+nt+j)
DotRhoEdgeAnnihilation(j) = -((state(ipc,ip,el)%p(j))**2)* & ! annihilation of edge dislocations
constitutive_titanmod_capre_PerSlipSys(j,instance)*state(ipc,ip,el)%p(7*ns+2*nt+j)*0.5_pReal
DotRhoScrewAnnihilation(j) = -((state(ipc,ip,el)%p(ns+j))**2)* & ! annihilation of screw dislocations
constitutive_titanmod_caprs_PerSlipSys(j,instance)*state(ipc,ip,el)%p(8*ns+2*nt+j)*0.5_pReal
state%p(j)*state%p(7*ns+2*nt+j)/state%p(3*ns+nt+j)
DotRhoEdgeAnnihilation(j) = -((state%p(j))**2)* & ! annihilation of edge dislocations
constitutive_titanmod_capre_PerSlipSys(j,instance)*state%p(7*ns+2*nt+j)*0.5_pReal
DotRhoScrewAnnihilation(j) = -((state%p(ns+j))**2)* & ! annihilation of screw dislocations
constitutive_titanmod_caprs_PerSlipSys(j,instance)*state%p(8*ns+2*nt+j)*0.5_pReal
constitutive_titanmod_dotState(j) = & ! edge dislocation density rate of change
DotRhoEdgeGeneration(j)+DotRhoEdgeAnnihilation(j)
@ -1652,7 +1649,7 @@ implicit none
DotRhoScrewGeneration(j)+DotRhoScrewAnnihilation(j)
constitutive_titanmod_dotState(2*ns+j) = & ! sum of shear due to edge and screw
state(ipc,ip,el)%p(10*ns+2*nt+j)+state(ipc,ip,el)%p(11*ns+2*nt+j)
state%p(10*ns+2*nt+j)+state%p(11*ns+2*nt+j)
enddo
enddo slipFamiliesLoop
@ -1668,7 +1665,7 @@ implicit none
!* Stress ratio for edge
twinStressRatio_p = ((abs(tau_twin(j)))/ &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state(ipc,ip,el)%p(7*ns+nt+j)) &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) &
)**(constitutive_titanmod_twinp_PerTwinSys(j,instance))
@ -1708,11 +1705,11 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
phase_Noutput
implicit none
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_titanmod_postResults
@ -1734,7 +1731,7 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
nt = constitutive_titanmod_totalNtwin(instance)
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
@ -1749,91 +1746,91 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_titanmod_outputID(o,instance))
case (rhoedge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(1_pInt:ns)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(1_pInt:ns)
c = c + ns
case (rhoscrew_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(ns+1_pInt:2_pInt*ns)
c = c + ns
case (segment_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt)
c = c + ns
case (segment_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt)
c = c + ns
case (resistance_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt)
c = c + ns
case (resistance_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt)
c = c + ns
case (velocity_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(7*ns+2*nt+1:8*ns+2*nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(7*ns+2*nt+1:8*ns+2*nt)
c = c + ns
case (velocity_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(8*ns+2*nt+1:9*ns+2*nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(8*ns+2*nt+1:9*ns+2*nt)
c = c + ns
case (tau_slip_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(9*ns+2*nt+1:10*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(9*ns+2*nt+1:10*ns+2*nt))
c = c + ns
case (gdot_slip_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(10*ns+2*nt+1:11*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(10*ns+2*nt+1:11*ns+2*nt))
c = c + ns
case (gdot_slip_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(11*ns+2*nt+1:12*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(11*ns+2*nt+1:12*ns+2*nt))
c = c + ns
case (gdot_slip_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(10*ns+2*nt+1:11*ns+2*nt)) + &
abs(state(ipc,ip,el)%p(11*ns+2*nt+1:12*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(10*ns+2*nt+1:11*ns+2*nt)) + &
abs(state%p(11*ns+2*nt+1:12*ns+2*nt))
c = c + ns
case (stressratio_edge_p_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(12*ns+2*nt+1:13*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(12*ns+2*nt+1:13*ns+2*nt))
c = c + ns
case (stressratio_screw_p_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(13*ns+2*nt+1:14*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(13*ns+2*nt+1:14*ns+2*nt))
c = c + ns
case (shear_system_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(2*ns+1:3*ns))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(2*ns+1:3*ns))
c = c + ns
case (shear_basal_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+1:2*ns+3)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+1:2*ns+3)))
c = c + 1_pInt
case (shear_prism_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+4:2*ns+6)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+4:2*ns+6)))
c = c + 1_pInt
case (shear_pyra_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+7:2*ns+12)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+7:2*ns+12)))
c = c + 1_pInt
case (shear_pyrca_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+13:2*ns+24)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+13:2*ns+24)))
c = c + 1_pInt
case (rhoedge_basal_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(1:3))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(1:3))
c = c + 1_pInt
case (rhoedge_prism_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(4:6))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(4:6))
c = c + 1_pInt
case (rhoedge_pyra_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(7:12))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(7:12))
c = c + 1_pInt
case (rhoedge_pyrca_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(13:24))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(13:24))
c = c + 1_pInt
case (rhoscrew_basal_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+1:ns+3))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+1:ns+3))
c = c + 1_pInt
case (rhoscrew_prism_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+4:ns+6))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+4:ns+6))
c = c + 1_pInt
case (rhoscrew_pyra_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+7:ns+12))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+7:ns+12))
c = c + 1_pInt
case (rhoscrew_pyrca_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+13:ns+24))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+13:ns+24))
c = c + 1_pInt
case (shear_total_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+1:3*ns)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+1:3*ns)))
c = c + 1_pInt
case (twin_fraction_ID)
constitutive_titanmod_postResults(c+1_pInt:c+nt) = abs(volumefraction_PerTwinSys(1:nt))