From 0b59519a2a1c54e19a46df13ae84d6f40e8cf430 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Thu, 29 Jan 2015 13:56:09 +0000 Subject: [PATCH] updated damage models: - coupling to plasticity handled within damage module instead of plasticity module - anisotropic models more stable --- code/constitutive.f90 | 130 +++++---------------------- code/damage_anisoBrittle.f90 | 134 +++++++++------------------- code/damage_anisoDuctile.f90 | 157 ++++++++++----------------------- code/damage_isoDuctile.f90 | 55 ++++++------ code/plastic_disloKMC.f90 | 16 +--- code/plastic_disloUCLA.f90 | 16 +--- code/plastic_dislotwin.f90 | 14 ++- code/plastic_j2.f90 | 14 +-- code/plastic_nonlocal.f90 | 19 ++-- code/plastic_phenopowerlaw.f90 | 19 ++-- code/plastic_titanmod.f90 | 8 +- 11 files changed, 165 insertions(+), 417 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 44dd4b6a9..f9d2a2009 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -38,7 +38,6 @@ module constitutive constitutive_getLocalDamage, & constitutive_putLocalDamage, & constitutive_getDamage, & - constitutive_getSlipDamage, & constitutive_getDamageDiffusion33, & constitutive_getAdiabaticTemperature, & constitutive_putAdiabaticTemperature, & @@ -701,23 +700,17 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) use plastic_j2, only: & plastic_j2_LpAndItsTangent use plastic_phenopowerlaw, only: & - plastic_phenopowerlaw_LpAndItsTangent, & - plastic_phenopowerlaw_totalNslip + plastic_phenopowerlaw_LpAndItsTangent use plastic_dislotwin, only: & - plastic_dislotwin_LpAndItsTangent, & - plastic_dislotwin_totalNslip + plastic_dislotwin_LpAndItsTangent use plastic_dislokmc, only: & - plastic_dislokmc_LpAndItsTangent, & - plastic_dislokmc_totalNslip + plastic_dislokmc_LpAndItsTangent use plastic_disloucla, only: & - plastic_disloucla_LpAndItsTangent, & - plastic_disloucla_totalNslip + plastic_disloucla_LpAndItsTangent use plastic_titanmod, only: & - plastic_titanmod_LpAndItsTangent, & - plastic_titanmod_totalNslip + plastic_titanmod_LpAndItsTangent use plastic_nonlocal, only: & - plastic_nonlocal_LpAndItsTangent, & - totalNslip + plastic_nonlocal_LpAndItsTangent implicit none integer(pInt), intent(in) :: & @@ -730,8 +723,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) Lp !< plastic velocity gradient real(pReal), intent(out), dimension(9,9) :: & dLp_dTstar !< derivative of Lp with respect to Tstar (4th-order tensor) - integer(pInt) :: & - nSlip select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -739,44 +730,28 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) Lp = 0.0_pReal dLp_dTstar = 0.0_pReal case (PLASTICITY_J2_ID) - nSlip = 1_pInt - call plastic_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & - ipc,ip,el) + call plastic_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - nSlip = plastic_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) - call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & - ipc,ip,el) + call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) - nSlip = totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) - nSlip = plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) - nSlip = plastic_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) - nSlip = plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_disloucla_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & ipc,ip,el) case (PLASTICITY_TITANMOD_ID) - nSlip = plastic_titanmod_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & ipc,ip,el) end select @@ -1163,23 +1138,17 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su use plastic_j2, only: & plastic_j2_dotState use plastic_phenopowerlaw, only: & - plastic_phenopowerlaw_dotState, & - plastic_phenopowerlaw_totalNslip + plastic_phenopowerlaw_dotState use plastic_dislotwin, only: & - plastic_dislotwin_dotState, & - plastic_dislotwin_totalNslip + plastic_dislotwin_dotState use plastic_dislokmc, only: & - plastic_dislokmc_dotState, & - plastic_dislokmc_totalNslip + plastic_dislokmc_dotState use plastic_disloucla, only: & - plastic_disloucla_dotState, & - plastic_disloucla_totalNslip + plastic_disloucla_dotState use plastic_titanmod, only: & - plastic_titanmod_dotState, & - plastic_titanmod_totalNslip + plastic_titanmod_dotState use plastic_nonlocal, only: & - plastic_nonlocal_dotState, & - totalNslip + plastic_nonlocal_dotState use damage_gurson, only: & damage_gurson_dotState @@ -1203,38 +1172,30 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su tick, tock, & tickrate, & maxticks - integer(pInt) :: & - nSlip if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) & call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_J2_ID) - nSlip = 1_pInt - call plastic_j2_dotState (Tstar_v,nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) + call plastic_j2_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - nSlip = plastic_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) - call plastic_phenopowerlaw_dotState(Tstar_v,nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) + call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) - nSlip = plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_dislotwin_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) + ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) - nSlip = plastic_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_dislokmc_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) + ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) - nSlip = plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) - call plastic_disloucla_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) + call plastic_disloucla_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & + ipc,ip,el) case (PLASTICITY_TITANMOD_ID) call plastic_titanmod_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) - nSlip = totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,constitutive_getTemperature(ipc,ip,el), & - nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),subdt,subfracArray,ip,el) + subdt,subfracArray,ip,el) end select select case (phase_damage(material_phase(ipc,ip,el))) @@ -1497,52 +1458,6 @@ function constitutive_getDamage(ipc, ip, el) end function constitutive_getDamage -!-------------------------------------------------------------------------------------------------- -!> @brief Returns the damage on each slip system -!-------------------------------------------------------------------------------------------------- -function constitutive_getSlipDamage(nSlip, Tstar_v, ipc, ip, el) - use prec, only: & - pReal - use material, only: & - material_phase, & - LOCAL_DAMAGE_isoDuctile_ID, & - LOCAL_DAMAGE_anisoDuctile_ID, & - LOCAL_DAMAGE_gurson_ID, & - phase_damage - use damage_isoDuctile, only: & - damage_isoDuctile_getSlipDamage - use damage_anisoDuctile, only: & - damage_anisoDuctile_getSlipDamage - use damage_gurson, only: & - damage_gurson_getSlipDamage - - implicit none - integer(pInt), intent(in) :: & - nSlip, & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress - real(pReal), dimension(nSlip) :: & - constitutive_getSlipDamage - - select case (phase_damage(material_phase(ipc,ip,el))) - case (LOCAL_DAMAGE_isoDuctile_ID) - constitutive_getSlipDamage = damage_isoDuctile_getSlipDamage(ipc, ip, el) - - case (LOCAL_DAMAGE_anisoDuctile_ID) - constitutive_getSlipDamage = damage_anisoDuctile_getSlipDamage(ipc, ip, el) - - case (LOCAL_DAMAGE_gurson_ID) - constitutive_getSlipDamage = damage_gurson_getSlipDamage(Tstar_v, ipc, ip, el) - - case default - constitutive_getSlipDamage = 1.0_pReal - end select - -end function constitutive_getSlipDamage - !-------------------------------------------------------------------------------------------------- !> @brief returns damage diffusion tensor !-------------------------------------------------------------------------------------------------- @@ -1569,12 +1484,13 @@ function constitutive_getDamageDiffusion33(ipc, ip, el) real(pReal), dimension(3,3) :: & constitutive_getDamageDiffusion33 - constitutive_getDamageDiffusion33 = lattice_DamageDiffusion33(1:3,1:3,material_phase(ipc,ip,el)) select case(phase_damage(material_phase(ipc,ip,el))) case (LOCAL_DAMAGE_isoBrittle_ID) constitutive_getDamageDiffusion33 = damage_isoBrittle_getDamageDiffusion33(ipc, ip, el) case (LOCAL_DAMAGE_phaseField_ID) constitutive_getDamageDiffusion33 = damage_phaseField_getDamageDiffusion33(ipc, ip, el) + case default + constitutive_getDamageDiffusion33 = lattice_DamageDiffusion33(1:3,1:3,material_phase(ipc,ip,el)) end select diff --git a/code/damage_anisoBrittle.f90 b/code/damage_anisoBrittle.f90 index 3b120b6ed..a117a86cb 100644 --- a/code/damage_anisoBrittle.f90 +++ b/code/damage_anisoBrittle.f90 @@ -257,9 +257,8 @@ subroutine damage_anisoBrittle_init(fileUnit) sizeDotState = 0_pInt sizeState = sizeDotState + & 1_pInt + & ! non-local damage - damage_anisoBrittle_totalNcleavage(instance) + & ! opening on each damage system - damage_anisoBrittle_totalNcleavage(instance) + & ! damage on each damage system - 9_pInt ! Fd + 1_pInt + & ! opening on each damage system + 9_pInt ! Fd damageState(phase)%sizeState = sizeState damageState(phase)%sizeDotState = sizeDotState @@ -304,16 +303,12 @@ subroutine damage_anisoBrittle_stateInit(phase,instance) real(pReal), dimension(damageState(phase)%sizeState) :: tempState - tempState(1) = 1.0_pReal - tempState(2 : & - 1 + damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal - tempState(2 + damage_anisoBrittle_totalNcleavage(instance): & - 1 +2*damage_anisoBrittle_totalNcleavage(instance)) = 1.0_pReal - tempState(2 +2*damage_anisoBrittle_totalNcleavage(instance): & - 10+2*damage_anisoBrittle_totalNcleavage(instance)) = reshape(math_I3, shape=[9]) - damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) - damageState(phase)%state0 = damageState(phase)%state - damageState(phase)%partionedState0 = damageState(phase)%state + tempState(1) = 1.0_pReal + tempState(2) = 0.0_pReal + tempState(3:11) = reshape(math_I3, shape=[9]) + + damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:))) + end subroutine damage_anisoBrittle_stateInit !-------------------------------------------------------------------------------------------------- @@ -329,20 +324,20 @@ subroutine damage_anisoBrittle_aTolState(phase,instance) instance ! number specifying the current instance of the damage real(pReal), dimension(damageState(phase)%sizeState) :: tempTol - tempTol(1) = damage_anisoBrittle_aTol_damage(instance) - tempTol(2 : & - 1 + damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_disp (instance) - tempTol(2 + damage_anisoBrittle_totalNcleavage(instance): & - 1 +2*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_damage(instance) - tempTol(2 +2*damage_anisoBrittle_totalNcleavage(instance): & - 10+2*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_damage(instance) + tempTol(1) = damage_anisoBrittle_aTol_damage(instance) + tempTol(2) = damage_anisoBrittle_aTol_disp (instance) + tempTol(3:11) = damage_anisoBrittle_aTol_damage(instance) + damageState(phase)%aTolState = tempTol + end subroutine damage_anisoBrittle_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el) + use numerics, only: & + residualStiffness use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -367,13 +362,9 @@ subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el) phase, & constituent, & instance, & - f, i, index_d, index_o, index_myFamily + f, i, index_myFamily real(pReal) :: & - localDamage, & - resForce, & - deltaState, & - available, & - predicted + localDamage real(pReal) :: & traction_d, traction_t, traction_n, traction_crit @@ -381,55 +372,29 @@ subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el) constituent = mappingConstitutive(1,ipc,ip,el) instance = phase_damageInstance(phase) - localDamage = minval(damageState(phase)%state(2+ damage_anisoBrittle_totalNcleavage(instance): & - 1+2*damage_anisoBrittle_totalNcleavage(instance),constituent)) - - index_o = 2_pInt - index_d = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) + damageState(phase)%state(2,constituent) = damageState(phase)%subState0(2,constituent) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase)) traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) - traction_crit = damage_anisoBrittle_critLoad(f,instance) - - deltaState = subdt*damage_anisoBrittle_sdot_0(instance)* & - ( (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) + & - (abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) + & - (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance)) - if (damageState(phase)%state(index_o,constituent) <= 1.0_pReal) then - damageState(phase)%state(index_o,constituent) = & - damageState(phase)%subState0(index_o,constituent) + & - deltaState - else - if (damageState(phase)%state(index_d,constituent) == localDamage) then - resForce = damage_anisoBrittle_getDamage(ipc, ip, el) - localDamage - else - resForce = 0.0_pReal - endif - available = (damageState(phase)%subState0(index_o,constituent) - resForce)** & - (1.0_pReal - damage_anisoBrittle_N(instance)) - predicted = -(1.0_pReal - damage_anisoBrittle_N(instance))*deltaState - if (available > predicted) then - damageState(phase)%state(index_o,constituent) = resForce + & - (available - predicted)**(1.0_pReal/(1.0_pReal - damage_anisoBrittle_N(instance))) - else - damageState(phase)%state(index_o,constituent) = huge(1.0_pReal) - endif - damageState(phase)%state(index_d,constituent) = & - min(damageState(phase)%state0(index_d,constituent) , & - 1.0_pReal/(max(0.0_pReal,damageState(phase)%state(index_o,constituent) - resForce))) - endif + traction_crit = damage_anisoBrittle_critLoad(f,instance)* & + damage_anisoBrittle_getDamage(ipc, ip, el) + damageState(phase)%state(2,constituent) = & + damageState(phase)%state(2,constituent) + & + subdt*damage_anisoBrittle_sdot_0(instance)* & + ( (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) + & + (abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) + & + (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance))/ & + damage_anisoBrittle_critDisp(f,instance) - index_d = index_d + 1_pInt; index_o = index_o + 1_pInt enddo enddo - localDamage = minval(damageState(phase)%state(2+ damage_anisoBrittle_totalNcleavage(instance): & - 1+2*damage_anisoBrittle_totalNcleavage(instance),constituent)) - + localDamage = min(1.0_pReal,max(residualStiffness,1.0/damageState(phase)%state(2,constituent))) + damageState(phase)%state(1,constituent) = & localDamage + & (damageState(phase)%subState0(1,constituent) - localDamage)* & @@ -441,8 +406,6 @@ end subroutine damage_anisoBrittle_microstructure !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) - use numerics, only: & - residualStiffness use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -472,7 +435,7 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, phase, & constituent, & instance, & - f, i, index, index_myFamily, k, l, m, n + f, i, index_myFamily, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt @@ -483,7 +446,6 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, Ld = 0.0_pReal dLd_dTstar3333 = 0.0_pReal - index = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family @@ -491,7 +453,7 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) traction_crit = damage_anisoBrittle_critLoad(f,instance)* & - (damageState(phase)%state0(index,constituent) + residualStiffness) + damage_anisoBrittle_getDamage(ipc, ip, el) udotd = & sign(1.0_pReal,traction_d)* & damage_anisoBrittle_sdot_0(instance)* & @@ -530,7 +492,6 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, lattice_Scleavage(m,n,3,index_myFamily+i,phase) endif - index = index + 1_pInt enddo enddo dLd_dTstar = math_Plain3333to99(dLd_dTstar3333) @@ -562,11 +523,7 @@ pure function damage_anisoBrittle_getFd(ipc, ip, el) constituent = mappingConstitutive(1,ipc,ip,el) instance = phase_damageInstance(phase) - damage_anisoBrittle_getFd = & - reshape(damageState(phase)% & - state(2*damage_anisoBrittle_totalNcleavage(instance)+2: & - 2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent), & - shape=[3,3]) + damage_anisoBrittle_getFd = reshape(damageState(phase)%state(3:11,constituent),shape=[3,3]) end function damage_anisoBrittle_getFd @@ -574,8 +531,6 @@ end function damage_anisoBrittle_getFd !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el) - use numerics, only: & - residualStiffness use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -605,7 +560,7 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el) phase, & constituent, & instance, & - f, i, index, index_myFamily + f, i, index_myFamily real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, udott, udotn @@ -615,7 +570,6 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el) instance = phase_damageInstance(phase) Ld = 0.0_pReal - index = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family @@ -623,7 +577,7 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el) traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) traction_crit = damage_anisoBrittle_critLoad(f,instance)* & - (damageState(phase)%state0(index,constituent) + residualStiffness) + damage_anisoBrittle_getDamage(ipc, ip, el) udotd = & sign(1.0_pReal,traction_d)* & damage_anisoBrittle_sdot_0(instance)* & @@ -641,11 +595,9 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el) (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance) Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) - index = index + 1_pInt enddo enddo - damageState(phase)%state(2*damage_anisoBrittle_totalNcleavage(instance)+2: & - 2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent) = & + damageState(phase)%state(3:11,constituent) = & reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), & damage_anisoBrittle_getFd0(ipc, ip, el)), shape=[9]) @@ -676,11 +628,7 @@ pure function damage_anisoBrittle_getFd0(ipc, ip, el) constituent = mappingConstitutive(1,ipc,ip,el) instance = phase_damageInstance(phase) - damage_anisoBrittle_getFd0 = & - reshape(damageState(phase)% & - subState0(2*damage_anisoBrittle_totalNcleavage(instance)+2: & - 2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent), & - shape=[3,3]) + damage_anisoBrittle_getFd0 = reshape(damageState(phase)%subState0(3:11,constituent),shape=[3,3]) end function damage_anisoBrittle_getFd0 @@ -710,10 +658,7 @@ pure function damage_anisoBrittle_getPartionedFd0(ipc, ip, el) instance = phase_damageInstance(phase) damage_anisoBrittle_getPartionedFd0 = & - reshape(damageState(phase)% & - partionedState0(2*damage_anisoBrittle_totalNcleavage(instance)+2: & - 2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent), & - shape=[3,3]) + reshape(damageState(phase)%partionedState0(3:11,constituent),shape=[3,3]) end function damage_anisoBrittle_getPartionedFd0 @@ -741,7 +686,7 @@ function damage_anisoBrittle_getDamage(ipc, ip, el) select case(field_damage_type(material_homog(ip,el))) case (FIELD_DAMAGE_LOCAL_ID) damage_anisoBrittle_getDamage = damageState(mappingConstitutive(2,ipc,ip,el))% & - state0(1,mappingConstitutive(1,ipc,ip,el)) + state(1,mappingConstitutive(1,ipc,ip,el)) case (FIELD_DAMAGE_NONLOCAL_ID) damage_anisoBrittle_getDamage = fieldDamage(material_homog(ip,el))% & @@ -799,7 +744,8 @@ end function damage_anisoBrittle_getLocalDamage function damage_anisoBrittle_postResults(ipc,ip,el) use material, only: & mappingConstitutive, & - phase_damageInstance + phase_damageInstance, & + damageState implicit none integer(pInt), intent(in) :: & diff --git a/code/damage_anisoDuctile.f90 b/code/damage_anisoDuctile.f90 index cb4fa0d89..62dbb867a 100644 --- a/code/damage_anisoDuctile.f90 +++ b/code/damage_anisoDuctile.f90 @@ -66,7 +66,6 @@ module damage_anisoDuctile damage_anisoDuctile_getDamage, & damage_anisoDuctile_putLocalDamage, & damage_anisoDuctile_getLocalDamage, & - damage_anisoDuctile_getSlipDamage, & damage_anisoDuctile_postResults contains @@ -254,7 +253,7 @@ subroutine damage_anisoDuctile_init(fileUnit) sizeDotState = 0_pInt sizeState = sizeDotState + & 1_pInt + & ! time regularised damage - damage_anisoDuctile_totalNslip(instance) + & ! slip system damages + 1_pInt + & ! damaged plasticity 9 ! Fd damageState(phase)%sizeState = sizeState damageState(phase)%sizeDotState = sizeDotState @@ -297,16 +296,11 @@ subroutine damage_anisoDuctile_stateInit(phase, instance) real(pReal), dimension(damageState(phase)%sizeState) :: tempState - tempState(1) = 1.0_pReal - tempState(2 : & - 1 + damage_anisoDuctile_totalNslip(instance)) = 1.0_pReal - tempState(damage_anisoDuctile_totalNslip(instance)+2: & - damage_anisoDuctile_totalNslip(instance)+10) = reshape(math_I3, shape=[9]) + tempState(1) = 1.0_pReal + tempState(2) = 0.0_pReal + tempState(3:11) = reshape(math_I3, shape=[9]) - damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) - damageState(phase)%state0 = damageState(phase)%state - damageState(phase)%partionedState0 = damageState(phase)%state - + damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:))) end subroutine damage_anisoDuctile_stateInit @@ -325,20 +319,22 @@ subroutine damage_anisoDuctile_aTolState(phase,instance) tempTol = damage_anisoDuctile_aTol_damage(instance) damageState(phase)%aTolState = tempTol + end subroutine damage_anisoDuctile_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine damage_anisoDuctile_microstructure(subdt, ipc, ip, el) + use numerics, only: & + residualStiffness use material, only: & mappingConstitutive, & phase_damageInstance, & plasticState, & damageState use lattice, only: & - lattice_maxNslipFamily - use lattice, only: & + lattice_maxNslipFamily, & lattice_DamageMobility implicit none @@ -354,36 +350,30 @@ subroutine damage_anisoDuctile_microstructure(subdt, ipc, ip, el) instance, & index, f, i real(pReal) :: & - localDamage, & - drivingForce + localDamage phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) instance = phase_damageInstance(phase) - localDamage = minval(damageState(phase)%state(2:1+damage_anisoDuctile_totalNslip(instance),constituent)) - index = 1_pInt + damageState(phase)%state(2,constituent) = damageState(phase)%subState0(2,constituent) do f = 1_pInt,lattice_maxNslipFamily - do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family - if (localDamage == damageState(phase)%state(index+1,constituent)) then - drivingForce = plasticState(phase)%accumulatedSlip(index,constituent)/damage_anisoDuctile_critPlasticStrain(f,instance) - & - damage_anisoDuctile_getDamage(ipc, ip, el) - damageState(phase)%state(index+1,constituent) = & - min(damageState(phase)%state0(index+1,constituent), & - (sqrt(drivingForce*drivingForce + 4.0_pReal) - drivingForce)/2.0_pReal) - else - drivingForce = plasticState(phase)%accumulatedSlip(index,constituent)/damage_anisoDuctile_critPlasticStrain(f,instance) - damageState(phase)%state(index+1,constituent) = & - min(damageState(phase)%state0(index+1,constituent), & - 1.0_pReal/drivingForce) ! irreversibility - endif + do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family + damageState(phase)%state(2,constituent) = & + damageState(phase)%state(2,constituent) + & + subdt* & + plasticState(phase)%slipRate(index,constituent)/ & + (damage_anisoDuctile_getDamage(ipc, ip, el)**damage_anisoDuctile_N(instance))/ & + damage_anisoDuctile_critPlasticStrain(f,instance) + index = index + 1_pInt enddo enddo - - localDamage = minval(damageState(phase)%state(2:1+damage_anisoDuctile_totalNslip(instance),constituent)) - + + localDamage = & + max(residualStiffness,min(1.0_pReal,1.0_pReal/damageState(phase)%state(2,constituent))) + damageState(phase)%state(1,constituent) = & localDamage + & (damageState(phase)%subState0(1,constituent) - localDamage)* & @@ -395,8 +385,6 @@ end subroutine damage_anisoDuctile_microstructure !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) - use numerics, only: & - residualStiffness use lattice, only: & lattice_maxNslipFamily, & lattice_NslipSystem, & @@ -436,7 +424,7 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, phase, & constituent, & instance, & - f, i, index, index_myFamily, k, l, m, n + f, i, index_myFamily, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt @@ -447,13 +435,9 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, Ld = 0.0_pReal dLd_dTstar3333 = 0.0_pReal - - - index = 2_pInt do f = 1_pInt,lattice_maxNslipFamily index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family - projection_d = math_tensorproduct(lattice_sd(1:3,index_myFamily+i,phase),& lattice_sn(1:3,index_myFamily+i,phase)) projection_t = math_tensorproduct(lattice_st(1:3,index_myFamily+i,phase),& @@ -470,12 +454,13 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, traction_n = dot_product(Tstar_v,projection_n_v(1:6)) traction_crit = damage_anisoDuctile_critLoad(f,instance)* & - (damageState(phase)%state0(index,constituent) + residualStiffness) ! degrading critical load carrying capacity by damage + damage_anisoDuctile_getDamage(ipc, ip, el) ! degrading critical load carrying capacity by damage udotd = & sign(1.0_pReal,traction_d)* & damage_anisoDuctile_sdot_0(instance)* & - (abs(traction_d)/traction_crit)**damage_anisoDuctile_N(instance) + (abs(traction_d)/traction_crit - & + abs(traction_d)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance) if (udotd /= 0.0_pReal) then Ld = Ld + udotd*projection_d dudotd_dt = udotd*damage_anisoDuctile_N(instance)/traction_d @@ -487,7 +472,8 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, udott = & sign(1.0_pReal,traction_t)* & damage_anisoDuctile_sdot_0(instance)* & - (abs(traction_t)/traction_crit)**damage_anisoDuctile_N(instance) + (abs(traction_t)/traction_crit - & + abs(traction_t)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance) if (udott /= 0.0_pReal) then Ld = Ld + udott*projection_t dudott_dt = udott*damage_anisoDuctile_N(instance)/traction_t @@ -497,15 +483,15 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, endif udotn = & damage_anisoDuctile_sdot_0(instance)* & - (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoDuctile_N(instance) + (max(0.0_pReal,traction_n)/traction_crit - & + max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance) if (udotn /= 0.0_pReal) then Ld = Ld + udotn*projection_n(1:3,1:3) dudotn_dt = udotn*damage_anisoDuctile_N(instance)/traction_n forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & dudotn_dt*projection_n(k,l)* projection_n(m,n) - endif - index = index + 1_pInt + endif enddo enddo @@ -539,10 +525,7 @@ pure function damage_anisoDuctile_getFd(ipc, ip, el) instance = phase_damageInstance(phase) damage_anisoDuctile_getFd = & - reshape(damageState(phase)% & - state(damage_anisoDuctile_totalNslip(instance)+2: & - damage_anisoDuctile_totalNslip(instance)+10,constituent), & - shape=[3,3]) + reshape(damageState(phase)%state(3:11,constituent),shape=[3,3]) end function damage_anisoDuctile_getFd @@ -550,8 +533,6 @@ end function damage_anisoDuctile_getFd !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) - use numerics, only: & - residualStiffness use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -559,9 +540,7 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) use math, only: & math_mul33x33, & math_inv33, & - math_Plain3333to99, & math_I3, & - math_identity4th, & math_symmetric33, & math_Mandel33to6, & math_tensorproduct @@ -591,7 +570,7 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) phase, & constituent, & instance, & - f, i, index, index_myFamily + f, i, index_myFamily real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, udott, udotn @@ -599,7 +578,7 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) instance = phase_damageInstance(phase) - index = 2_pInt + Ld = 0.0_pReal do f = 1_pInt,lattice_maxNslipFamily index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family @@ -621,12 +600,13 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) traction_n = dot_product(Tstar_v,projection_n_v(1:6)) traction_crit = damage_anisoDuctile_critLoad(f,instance)* & - (damageState(phase)%state0(index,constituent) + residualStiffness) ! degrading critical load carrying capacity by damage + damage_anisoDuctile_getDamage(ipc, ip, el) ! degrading critical load carrying capacity by damage udotd = & sign(1.0_pReal,traction_d)* & damage_anisoDuctile_sdot_0(instance)* & - (abs(traction_d)/traction_crit)**damage_anisoDuctile_N(instance) + (abs(traction_d)/traction_crit - & + abs(traction_d)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance) if (udotd /= 0.0_pReal) then Ld = Ld + udotd*projection_d endif @@ -634,26 +614,25 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) udott = & sign(1.0_pReal,traction_t)* & damage_anisoDuctile_sdot_0(instance)* & - (abs(traction_t)/traction_crit)**damage_anisoDuctile_N(instance) + (abs(traction_t)/traction_crit) - & + abs(traction_t)/damage_anisoDuctile_critLoad(f,instance)**damage_anisoDuctile_N(instance) if (udott /= 0.0_pReal) then Ld = Ld + udott*projection_t endif udotn = & damage_anisoDuctile_sdot_0(instance)* & - (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoDuctile_N(instance) + (max(0.0_pReal,traction_n)/traction_crit - & + max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance) if (udotn /= 0.0_pReal) then Ld = Ld + udotn*projection_n(1:3,1:3) endif - - index = index + 1_pInt enddo enddo - damageState(phase)%state(damage_anisoDuctile_totalNslip(instance)+2: & - damage_anisoDuctile_totalNslip(instance)+10,constituent) = & - reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), & - damage_anisoDuctile_getFd0(ipc, ip, el)), shape=[9]) + damageState(phase)%state(3:11,constituent) = reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), & + damage_anisoDuctile_getFd0(ipc, ip, el)), & + shape=[9]) end subroutine damage_anisoDuctile_putFd @@ -683,10 +662,7 @@ pure function damage_anisoDuctile_getFd0(ipc, ip, el) instance = phase_damageInstance(phase) damage_anisoDuctile_getFd0 = & - reshape(damageState(phase)% & - subState0(damage_anisoDuctile_totalNslip(instance)+2: & - damage_anisoDuctile_totalNslip(instance)+10,constituent), & - shape=[3,3]) + reshape(damageState(phase)%subState0(3:11,constituent),shape=[3,3]) end function damage_anisoDuctile_getFd0 @@ -716,10 +692,7 @@ pure function damage_anisoDuctile_getPartionedFd0(ipc, ip, el) instance = phase_damageInstance(phase) damage_anisoDuctile_getPartionedFd0 = & - reshape(damageState(phase)% & - partionedState0(damage_anisoDuctile_totalNslip(instance)+2: & - damage_anisoDuctile_totalNslip(instance)+10,constituent), & - shape=[3,3]) + reshape(damageState(phase)%partionedState0(3:11,constituent),shape=[3,3]) end function damage_anisoDuctile_getPartionedFd0 @@ -747,7 +720,7 @@ function damage_anisoDuctile_getDamage(ipc, ip, el) select case(field_damage_type(material_homog(ip,el))) case (FIELD_DAMAGE_LOCAL_ID) damage_anisoDuctile_getDamage = damageState(mappingConstitutive(2,ipc,ip,el))% & - state0(1,mappingConstitutive(1,ipc,ip,el)) + state(1,mappingConstitutive(1,ipc,ip,el)) case (FIELD_DAMAGE_NONLOCAL_ID) damage_anisoDuctile_getDamage = fieldDamage(material_homog(ip,el))% & @@ -799,40 +772,6 @@ function damage_anisoDuctile_getLocalDamage(ipc, ip, el) end function damage_anisoDuctile_getLocalDamage -!-------------------------------------------------------------------------------------------------- -!> @brief returns slip system damage -!-------------------------------------------------------------------------------------------------- -function damage_anisoDuctile_getSlipDamage(ipc, ip, el) - use numerics, only: & - residualStiffness - use material, only: & - mappingConstitutive, & - phase_damageInstance, & - damageState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal) :: & - damage_anisoDuctile_getSlipDamage(damage_anisoDuctile_totalNslip( & - phase_damageInstance(mappingConstitutive(2,ipc,ip,el)))) - integer(pInt) :: & - phase, & - constituent, & - instance - - phase = mappingConstitutive(2,ipc,ip,el) - constituent = mappingConstitutive(1,ipc,ip,el) - instance = phase_damageInstance(phase) - - damage_anisoDuctile_getSlipDamage = & - damageState(phase)%state0(2:1+damage_anisoDuctile_totalNslip(instance),constituent) + & - residualStiffness - -end function damage_anisoDuctile_getSlipDamage - !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- diff --git a/code/damage_isoDuctile.f90 b/code/damage_isoDuctile.f90 index 53fb58cd8..f4f28f1d0 100644 --- a/code/damage_isoDuctile.f90 +++ b/code/damage_isoDuctile.f90 @@ -27,7 +27,8 @@ module damage_isoDuctile real(pReal), dimension(:), allocatable, private :: & damage_isoDuctile_aTol, & - damage_isoDuctile_critPlasticStrain + damage_isoDuctile_critPlasticStrain, & + damage_isoDuctile_N enum, bind(c) enumerator :: undefined_ID, & @@ -44,7 +45,6 @@ module damage_isoDuctile damage_isoDuctile_aTolState, & damage_isoDuctile_microstructure, & damage_isoDuctile_getDamage, & - damage_isoDuctile_getSlipDamage, & damage_isoDuctile_putLocalDamage, & damage_isoDuctile_getLocalDamage, & damage_isoDuctile_getDamagedC66, & @@ -121,6 +121,7 @@ subroutine damage_isoDuctile_init(fileUnit) allocate(damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) allocate(damage_isoDuctile_Noutput(maxNinstance), source=0_pInt) allocate(damage_isoDuctile_critPlasticStrain(maxNinstance), source=0.0_pReal) + allocate(damage_isoDuctile_N(maxNinstance), source=0.0_pReal) allocate(damage_isoDuctile_aTol(maxNinstance), source=0.0_pReal) rewind(fileUnit) @@ -157,6 +158,9 @@ subroutine damage_isoDuctile_init(fileUnit) case ('criticalplasticstrain') damage_isoDuctile_critPlasticStrain(instance) = IO_floatValue(line,positions,2_pInt) + case ('damageratesensitivity') + damage_isoDuctile_N(instance) = IO_floatValue(line,positions,2_pInt) + case ('atol_damage') damage_isoDuctile_aTol(instance) = IO_floatValue(line,positions,2_pInt) @@ -239,10 +243,11 @@ subroutine damage_isoDuctile_stateInit(phase) real(pReal), dimension(damageState(phase)%sizeState) :: tempState - tempState = 1.0_pReal - damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) - damageState(phase)%state0 = damageState(phase)%state - damageState(phase)%partionedState0 = damageState(phase)%state + tempState(1) = 1.0_pReal + tempState(2) = 0.0_pReal + + damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:))) + end subroutine damage_isoDuctile_stateInit !-------------------------------------------------------------------------------------------------- @@ -260,6 +265,7 @@ subroutine damage_isoDuctile_aTolState(phase,instance) tempTol = damage_isoDuctile_aTol(instance) damageState(phase)%aTolState = tempTol + end subroutine damage_isoDuctile_aTolState !-------------------------------------------------------------------------------------------------- @@ -281,24 +287,30 @@ subroutine damage_isoDuctile_microstructure(subdt,ipc, ip, el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & subdt integer(pInt) :: & phase, constituent, instance + real(pReal) :: & + localDamage phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) instance = phase_damageInstance(phase) damageState(phase)%state(2,constituent) = & - max(residualStiffness, & - min(damageState(phase)%state0(2,constituent), & - damage_isoDuctile_critPlasticStrain(instance)/ & - sum(plasticState(phase)%accumulatedSlip(:,constituent)))) + damageState(phase)%subState0(2,constituent) + & + subdt* & + sum(plasticState(phase)%slipRate(:,constituent))/ & + (damage_isoDuctile_getDamage(ipc, ip, el)**damage_isoDuctile_N(instance))/ & + damage_isoDuctile_critPlasticStrain(instance) + + localDamage = & + max(residualStiffness,min(1.0_pReal, 1.0_pReal/damageState(phase)%state(2,constituent))) damageState(phase)%state(1,constituent) = & - damageState(phase)%state(2,constituent) + & - (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* & + localDamage + & + (damageState(phase)%subState0(1,constituent) - localDamage)* & exp(-subdt/lattice_DamageMobility(phase)) end subroutine damage_isoDuctile_microstructure @@ -337,23 +349,6 @@ function damage_isoDuctile_getDamage(ipc, ip, el) end function damage_isoDuctile_getDamage -!-------------------------------------------------------------------------------------------------- -!> @brief returns slip damage -!-------------------------------------------------------------------------------------------------- -function damage_isoDuctile_getSlipDamage(ipc, ip, el) - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal) :: damage_isoDuctile_getSlipDamage, damage - - damage = damage_isoDuctile_getDamage(ipc, ip, el) - damage_isoDuctile_getSlipDamage = damage*damage - -end function damage_isoDuctile_getSlipDamage - !-------------------------------------------------------------------------------------------------- !> @brief puts local damage !-------------------------------------------------------------------------------------------------- diff --git a/code/plastic_disloKMC.f90 b/code/plastic_disloKMC.f90 index 89d342af8..87513436d 100644 --- a/code/plastic_disloKMC.f90 +++ b/code/plastic_disloKMC.f90 @@ -1112,8 +1112,7 @@ end subroutine plastic_disloKMC_microstructure !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage, & - ipc,ip,el) +subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check use math, only: & @@ -1145,10 +1144,9 @@ subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature, LATTICE_fcc_ID implicit none - integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el + integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: Temperature real(pReal), dimension(6), intent(in) :: Tstar_v - real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 @@ -1206,9 +1204,6 @@ subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature, nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_disloKMC_nonSchmidCoeff(k,instance)*& lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) enddo nonSchmidSystems - !* Applying damage to slip system - tau_slip_pos = tau_slip_pos/slipDamage(j) - tau_slip_neg = tau_slip_neg/slipDamage(j) significantPostitiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then !* Stress ratios @@ -1354,7 +1349,7 @@ end subroutine plastic_disloKMC_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check use math, only: & @@ -1385,12 +1380,9 @@ subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage, real(pReal), intent(in) :: & temperature !< temperature at integration point integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage integer(pInt) :: instance,ns,nt,f,i,j,k,index_myFamily,s1,s2, & ph, & @@ -1457,8 +1449,6 @@ subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage, tau_slip_neg = tau_slip_neg + plastic_disloKMC_nonSchmidCoeff(k,instance)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) enddo nonSchmidSystems - tau_slip_pos = tau_slip_pos/slipDamage(j) - tau_slip_neg = tau_slip_pos/slipDamage(j) significantPositiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then !* Stress ratios diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90 index 3529b6b56..77955c14d 100644 --- a/code/plastic_disloUCLA.f90 +++ b/code/plastic_disloUCLA.f90 @@ -1130,8 +1130,7 @@ end subroutine plastic_disloUCLA_microstructure !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage, & - ipc,ip,el) +subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check use math, only: & @@ -1163,10 +1162,9 @@ subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature LATTICE_fcc_ID implicit none - integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el + integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: Temperature real(pReal), dimension(6), intent(in) :: Tstar_v - real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 @@ -1224,9 +1222,6 @@ subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_disloUCLA_nonSchmidCoeff(k,instance)*& lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) enddo nonSchmidSystems - !* Applying damage to slip system - tau_slip_pos = tau_slip_pos/slipDamage(j) - tau_slip_neg = tau_slip_neg/slipDamage(j) significantPostitiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then !* Stress ratio @@ -1437,7 +1432,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check use math, only: & @@ -1468,12 +1463,9 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage real(pReal), intent(in) :: & temperature !< temperature at integration point integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage integer(pInt) :: instance,ns,nt,f,i,j,k,index_myFamily,s1,s2, & ph, & @@ -1539,8 +1531,6 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage tau_slip_neg = tau_slip_neg + plastic_disloUCLA_nonSchmidCoeff(k,instance)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) enddo nonSchmidSystems - tau_slip_pos = tau_slip_pos/slipDamage(j) - tau_slip_neg = tau_slip_pos/slipDamage(j) significantPositiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then !* Stress ratios diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index 25e92f1b7..ccd9438bb 100644 --- a/code/plastic_dislotwin.f90 +++ b/code/plastic_dislotwin.f90 @@ -1363,7 +1363,7 @@ end subroutine plastic_dislotwin_microstructure !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check use math, only: & @@ -1399,10 +1399,9 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature LATTICE_fcc_ID implicit none - integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el + integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: Temperature real(pReal), dimension(6), intent(in) :: Tstar_v - real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 @@ -1461,7 +1460,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature !* Calculation of Lp !* Resolved shear stress on slip system - tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))/slipDamage(j) + tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) if((abs(tau_slip(j))-plasticState(ph)%state(7*ns+4*nt+2*nr+j, of)) > tol_math_check) then !* Stress ratios @@ -1668,7 +1667,7 @@ end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el) use prec, only: & tol_math_check use math, only: & @@ -1706,12 +1705,9 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage real(pReal), intent(in) :: & temperature !< temperature at integration point integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage integer(pInt) :: instance,ns,nt,nr,f,i,j,index_myFamily,s1,s2,a,b,sa,sb,ssa,ssb, & ph, & @@ -1751,7 +1747,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage j = j+1_pInt !* Resolved shear stress on slip system - tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))/slipDamage(j) + tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) if((abs(tau_slip(j))-plasticState(ph)%state(7*ns+4*nt+2*nr+j, of)) > tol_math_check) then !* Stress ratios diff --git a/code/plastic_j2.f90 b/code/plastic_j2.f90 index 00feb3072..a407520a9 100644 --- a/code/plastic_j2.f90 +++ b/code/plastic_j2.f90 @@ -362,7 +362,7 @@ end subroutine plastic_j2_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) use math, only: & math_mul6x6, & math_Mandel6to33, & @@ -388,12 +388,9 @@ subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDa real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage real(pReal), dimension(3,3) :: & Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor @@ -417,7 +414,7 @@ subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDa dLp_dTstar99 = 0.0_pReal else gamma_dot = plastic_j2_gdot0(instance) & - * (sqrt(1.5_pReal) * norm_Tstar_dev / (slipDamage(1)*plastic_j2_fTaylor(instance) * & + * (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * & plasticState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))) & **plastic_j2_n(instance) @@ -441,7 +438,7 @@ end subroutine plastic_j2_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_j2_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el) use math, only: & math_mul6x6 use mesh, only: & @@ -458,12 +455,9 @@ subroutine plastic_j2_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el) real(pReal), dimension(6), intent(in):: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage real(pReal), dimension(6) :: & Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal) :: & @@ -490,7 +484,7 @@ subroutine plastic_j2_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el) ! strain rate gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!----------------------------------------------------------------------------------- - (slipDamage(1)*plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance) + (plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance) !-------------------------------------------------------------------------------------------------- ! hardening coefficient diff --git a/code/plastic_nonlocal.f90 b/code/plastic_nonlocal.f90 index 77d478b44..819fc9e96 100644 --- a/code/plastic_nonlocal.f90 +++ b/code/plastic_nonlocal.f90 @@ -2021,8 +2021,7 @@ end subroutine plastic_nonlocal_kinetics !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, nSlipDamage, slipDamage, & - ipc, ip, el) +subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, ipc, ip, el) use math, only: math_Plain3333to99, & math_mul6x6, & @@ -2048,14 +2047,11 @@ use mesh, only: mesh_ipVolume implicit none !*** input variables -integer(pInt), intent(in) :: nSlipDamage, & - ipc, & +integer(pInt), intent(in) :: ipc, & ip, & !< current integration point el !< current element number real(pReal), intent(in) :: Temperature !< temperature real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation -real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage !*** output variables @@ -2121,7 +2117,7 @@ tauThreshold = plasticState(ph)%state(iTauF(1:ns,instance),of) do s = 1_pInt,ns sLattice = slipSystemLattice(s,instance) - tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))/slipDamage(s) + tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) tauNS(s,1) = tau(s) tauNS(s,2) = tau(s) if (tau(s) > 0.0_pReal) then @@ -2408,7 +2404,7 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature,nSlipDamage,slipDamage, & +subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, & timestep,subfrac, ip,el) use prec, only: DAMASK_NaN @@ -2463,8 +2459,7 @@ implicit none !*** input variables integer(pInt), intent(in) :: ip, & !< current integration point - el, & !< current element number - nSlipDamage + el !< current element number real(pReal), intent(in) :: Temperature, & !< temperature timestep !< substepped crystallite time increment real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation @@ -2473,8 +2468,6 @@ real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient -real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage !*** local variables @@ -2639,7 +2632,7 @@ forall (t = 1_pInt:4_pInt) & do s = 1_pInt,ns ! loop over slip systems sLattice = slipSystemLattice(s,instance) - tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))/slipDamage(s) + tauBack(s) + tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo diff --git a/code/plastic_phenopowerlaw.f90 b/code/plastic_phenopowerlaw.f90 index 85cb0ab01..c1048d1cf 100644 --- a/code/plastic_phenopowerlaw.f90 +++ b/code/plastic_phenopowerlaw.f90 @@ -717,8 +717,7 @@ end subroutine plastic_phenopowerlaw_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDamage, & - ipc,ip,el) +subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) use math, only: & math_Plain3333to99, & math_Mandel6to33 @@ -745,14 +744,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDa dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage integer(pInt) :: & instance, & @@ -807,11 +803,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDa lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) enddo gdot_slip_pos = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* & - ((abs(tau_slip_pos)/(slipDamage(j)*plasticState(ph)%state(j,of))) & + ((abs(tau_slip_pos)/(plasticState(ph)%state(j,of))) & **plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_pos) gdot_slip_neg = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* & - ((abs(tau_slip_neg)/(slipDamage(j)*plasticState(ph)%state(j,of))) & + ((abs(tau_slip_neg)/(plasticState(ph)%state(j,of))) & **plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_neg) Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F @@ -871,7 +867,7 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el) +subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) use lattice, only: & lattice_Sslip_v, & lattice_Stwin_v, & @@ -891,12 +887,9 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip, real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element !< microstructure state - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage integer(pInt) :: & instance,ph, & @@ -968,8 +961,8 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip, dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) enddo nonSchmidSystems gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & - ((abs(tau_slip_pos)/(slipDamage(j)*plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) & - +(abs(tau_slip_neg)/(slipDamage(j)*plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))& + ((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) & + +(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))& *sign(1.0_pReal,tau_slip_pos) enddo slipSystems1 enddo slipFamilies1 diff --git a/code/plastic_titanmod.f90 b/code/plastic_titanmod.f90 index 698976bb3..049340084 100644 --- a/code/plastic_titanmod.f90 +++ b/code/plastic_titanmod.f90 @@ -1332,8 +1332,7 @@ end subroutine plastic_titanmod_microstructure !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,nSlipDamage,slipDamage, & - ipc,ip,el) +subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,ipc,ip,el) use math, only: & math_Plain3333to99, & math_Mandel6to33 @@ -1365,7 +1364,6 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature, dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress integer(pInt), intent(in) :: & - nSlipDamage, & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -1373,8 +1371,6 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature, Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), intent(in) :: & temperature !< temperature at IP - real(pReal), dimension(nSlipDamage), intent(in) :: & - slipDamage integer(pInt) :: & index_myFamily, instance, & ns,nt, & @@ -1432,7 +1428,7 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature, !* Calculation of Lp !* Resolved shear stress on slip system - tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))/slipDamage(j) + tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) if(lattice_structure(ph)==LATTICE_hex_ID) then ! only for prismatic and pyr systems in hex screwvelocity_prefactor=plastic_titanmod_debyefrequency(instance)* & plasticState(ph)%state(4_pInt*ns+nt+j, of)*(plastic_titanmod_burgersPerSlipSys(j,instance)/ &