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