diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 2dc684855..05a06b9a7 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -397,9 +397,9 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip Lp = 0.0_pReal dLp_dTstar = math_identity2nd(9) case (PLASTICITY_J2_ID) - call constitutive_j2_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,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,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) case (PLASTICITY_DISLOTWIN_ID) @@ -455,8 +455,8 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) math_mul3333xx33, & math_Mandel66to3333, & math_transpose33, & - MATH_I3, & - math_trace33 + math_trace33, & + math_I3 use material, only: & mappingConstitutive use lattice, only: & @@ -480,42 +480,30 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) dT_dFe !< dT/dFe integer(pInt) :: i, j, k, l - real(pReal), dimension(3,3) :: FeT,strain, CxxDel, CxxDel_undamaged - real(pReal), dimension(3,3,3,3) :: C, C_undamged - real(pReal) :: strain_trace + real(pReal) :: damage, negative_volStrain, negative_volStress, Csum + real(pReal), dimension(3,3) :: strain + real(pReal), dimension(3,3,3,3) :: C - C_undamged = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) - C = constitutive_damageValue(ipc,ip,el)*& - math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) + C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) + damage = constitutive_damageValue(ipc,ip,el) + strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) + negative_volStrain = min(0.0_pReal,math_trace33(strain)/3.0_pReal) + negative_volStress = math_trace33(math_mul3333xx33(C,negative_volStrain*math_I3))/3.0_pReal + T = damage* & + math_mul3333xx33(C,strain - lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))* & + (constitutive_temperature(ipc,ip,el) - & + lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))) + & + (1.0_pReal - damage)*negative_volStress*math_I3 + + Csum = sum(math_I3*math_mul3333xx33(C,math_I3))/9.0_pReal + C = damage*C + forall (i = 1_pInt:3_pInt) & + C(1:3,1:3,i,i) = C(1:3,1:3,i,i) + (1.0_pReal - damage)*Csum*math_I3 - FeT = math_transpose33(Fe) - strain = 0.5_pReal*(math_mul33x33(FeT,Fe)-MATH_I3) - strain_trace = math_trace33(strain) - - - if(strain_trace>=0.0_pReal) then - T = math_mul3333xx33(C,strain) - & - lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))* & - (constitutive_temperature(ipc,ip,el) - & - lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el))) 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) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) ! dT*_ij/dFe_kl - else - strain = strain - (strain_trace/3)*math_I3 ! removing (negative) volumetric part - - T = math_mul3333xx33(C,strain) + & - math_mul3333xx33(C_undamged,((strain_trace/3)*math_I3 )) - & - lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))* & - (constitutive_temperature(ipc,ip,el) - & - lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el))) - CxxDel = math_mul3333xx33(C,math_I3) ! temporary variable C_ijxy * del_xy (damaged part) - CxxDel_undamaged = math_mul3333xx33(C_undamged,math_I3) ! temporary variable C_ijxy * del_xy (Undamaged part) - 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) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) - & - 0.5*CxxDel(i,j)*Fe(k,l) + & - 0.5*CxxDel_undamaged(i,j)*Fe(k,l) ! dT*_ij/dFe_kl = C_ijlt*Fe_kt - 0.5*Cxxdel_ij*Fe_kl - endif + end subroutine constitutive_hooke_TandItsTangent @@ -549,6 +537,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, PLASTICITY_DISLOKMC_ID, & PLASTICITY_TITANMOD_ID, & PLASTICITY_NONLOCAL_ID + use constitutive_damage, only: & + constitutive_damageValue use constitutive_j2, only: & constitutive_j2_dotState use constitutive_phenopowerlaw, only: & @@ -587,9 +577,9 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_J2_ID) - call constitutive_j2_dotState (Tstar_v,ipc,ip,el) + call constitutive_j2_dotState (Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - call constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) + call constitutive_phenopowerlaw_dotState(Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) call constitutive_dislotwin_dotState (Tstar_v,Temperature,ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) @@ -697,6 +687,8 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el) PLASTICITY_DISLOKMC_ID, & PLASTICITY_TITANMOD_ID, & PLASTICITY_NONLOCAL_ID + use constitutive_damage, only: & + constitutive_damageValue use constitutive_j2, only: & #ifdef HDF constitutive_j2_postResults2,& @@ -734,7 +726,7 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el) case (PLASTICITY_J2_ID) constitutive_postResults= constitutive_j2_postResults (Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,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) case (PLASTICITY_DISLOKMC_ID) diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 04a1ae7af..2b73c2d7b 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,ipc,ip,el) +subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip,el) use math, only: & math_mul6x6, & math_Mandel6to33, & @@ -368,6 +368,8 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) 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 @@ -395,7 +397,7 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) dLp_dTstar99 = 0.0_pReal else gamma_dot = constitutive_j2_gdot0(instance) & - * (sqrt(1.5_pReal) * norm_Tstar_dev / (constitutive_j2_fTaylor(instance) * & + * (sqrt(1.5_pReal) * norm_Tstar_dev / (damage * constitutive_j2_fTaylor(instance) * & plasticState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))) & **constitutive_j2_n(instance) @@ -419,7 +421,7 @@ end subroutine constitutive_j2_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el) +subroutine constitutive_j2_dotState(Tstar_v,damage,ipc,ip,el) use math, only: & math_mul6x6 use mesh, only: & @@ -441,6 +443,8 @@ subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el) el !< element real(pReal), dimension(6) :: & Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + damage real(pReal) :: & gamma_dot, & !< strainrate hardening, & !< hardening coefficient @@ -465,7 +469,7 @@ subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el) ! strain rate gamma_dot = constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!----------------------------------------------------------------------------------- - (constitutive_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**constitutive_j2_n(instance) + (damage*constitutive_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**constitutive_j2_n(instance) !-------------------------------------------------------------------------------------------------- ! hardening coefficient diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index d8c3f93d8..ed8340916 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,ipc,ip,el) +subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip,el) use prec, only: & p_vec use math, only: & @@ -696,6 +696,8 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,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 @@ -753,11 +755,11 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ip 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))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*& + ((abs(tau_slip_pos(j))/(damage*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))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*& + ((abs(tau_slip_neg(j))/(damage*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 @@ -795,7 +797,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ip 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))/(damage*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) @@ -819,7 +821,7 @@ end subroutine constitutive_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) +subroutine constitutive_phenopowerlaw_dotState(Tstar_v,damage,ipc,ip,el) use lattice, only: & lattice_Sslip_v, & lattice_Stwin_v, & @@ -842,6 +844,8 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,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 @@ -916,8 +920,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))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance) & + +(abs(tau_slip_neg(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))& *sign(1.0_pReal,tau_slip_pos(j)) enddo enddo slipFamiliesLoop1 @@ -938,7 +942,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))/(damage*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 @@ -988,7 +992,7 @@ end subroutine constitutive_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) +function constitutive_phenopowerlaw_postResults(Tstar_v,damage,ipc,ip,el) use prec, only: & p_vec use mesh, only: & @@ -1016,6 +1020,8 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,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 @@ -1073,8 +1079,8 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,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)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) & - +(abs(tau_slip_neg)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))& + ((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))& *sign(1.0_pReal,tau_slip_pos) enddo @@ -1116,7 +1122,7 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,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)/plasticState(ph)%state(j+nSlip,of))**& + (abs(tau)/(damage*plasticState(ph)%state(j+nSlip,of)))**& constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau)) enddo enddo twinFamiliesLoop1 diff --git a/code/damage_local.f90 b/code/damage_local.f90 index b57ebd3f8..f51283734 100644 --- a/code/damage_local.f90 +++ b/code/damage_local.f90 @@ -255,6 +255,7 @@ end subroutine damage_local_aTolState subroutine damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el) use material, only: & mappingConstitutive, & + phase_damageInstance, & damageState use math, only: & math_Mandel66to3333, & @@ -279,12 +280,13 @@ subroutine damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el) Lp, & Fe integer(pInt) :: & - phase, constituent + phase, constituent, instance real(pReal) :: & trialDamage, strain(3,3), stress(3,3), negative_volStrain phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) ! which instance of my damage is present phase strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) negative_volStrain = min(0.0_pReal,math_trace33(strain)/3.0_pReal)