diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e1e6ae47d..a8070bd52 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -35,6 +35,103 @@ module constitutive module subroutine thermal_init end subroutine thermal_init + + module subroutine plastic_isotropic_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_dotState + + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_dotState + + module subroutine plastic_kinehardening_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_dotState + + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_dotState + + module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dotState + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & + instance,of,ip,el) + real(pReal), dimension(3,3), intent(in) ::& + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + F, & !< deformation gradient + Fp !< plastic deformation gradient + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_dotState + + + module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S + end subroutine source_damage_anisoBrittle_dotState + + module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_anisoDuctile_dotState + + module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_isoDuctile_dotState + + module subroutine source_thermal_externalheat_dotState(phase, of) + + integer, intent(in) :: & + phase, & + of + end subroutine source_thermal_externalheat_dotState + + pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) integer, intent(in) :: & @@ -46,14 +143,14 @@ module constitutive end function kinematics_thermal_expansion_initialStrain - module function plastic_homogenizedC(ipc,ip,el) result(homogenizedC) + module function constitutive_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - end function plastic_homogenizedC + end function constitutive_homogenizedC module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) @@ -109,65 +206,6 @@ module constitutive dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) end subroutine kinematics_thermal_expansion_LiAndItsTangent - module function plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_plastic) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_plastic - end function plastic_dotState - - module function damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_damage) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_damage - end function damage_dotState - - module function thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_thermal) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_thermal - end function thermal_dotState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & @@ -261,7 +299,7 @@ module constitutive Fp !< plastic deformation gradient end subroutine plastic_dependentState - module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -276,7 +314,7 @@ module constitutive dLp_dS, & dLp_dFi !< derivative of Lp with respect to Fi - end subroutine plastic_LpAndItsTangents + end subroutine constitutive_plastic_LpAndItsTangents module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & @@ -331,7 +369,7 @@ module constitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_LpAndItsTangents, & + constitutive_plastic_LpAndItsTangents, & constitutive_dependentState, & constitutive_LiAndItsTangents, & constitutive_initialFi, & @@ -392,25 +430,6 @@ subroutine constitutive_init end subroutine constitutive_init - -subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - - call plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S, Fi, ipc, ip, el) - -end subroutine constitutive_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models @@ -430,23 +449,6 @@ subroutine constitutive_dependentState(F, Fp, ipc, ip, el) end subroutine constitutive_dependentState -!-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenize elasticity matrix -!> ToDo: homogenizedC66 would be more consistent -!-------------------------------------------------------------------------------------------------- -function constitutive_homogenizedC(ipc,ip,el) - - real(pReal), dimension(6,6) :: constitutive_homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - constitutive_homogenizedC = plastic_homogenizedC(ipc,ip,el) - -end function constitutive_homogenizedC - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? @@ -654,19 +656,65 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola Kirchhoff stress (vector notation) + real(pReal), dimension(3,3) :: & + Mp + integer :: & + ho, & + tme, & + i, & + instance logical :: broken - logical :: & - broken_plastic, & - broken_thermal, & - broken_damage + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + instance = phase_plasticityInstance(phase) - broken_plastic = plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) - broken = broken_plastic - broken_damage = damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) - broken = broken .or. broken_damage - broken_thermal = thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) - broken = broken .or. broken_thermal + Mp = matmul(matmul(transpose(Fi),Fi),S) + + plasticityType: select case (phase_plasticity(phase)) + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_dotState (Mp,instance,of) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_dotState(Mp,instance,of) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_dotState(Mp,instance,of) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & + instance,of,ip,el) + end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) + + SourceLoop: do i = 1, phase_Nsources(phase) + + sourceType: select case (phase_source(i,phase)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_dotState ( ipc, ip, el) + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_dotState ( ipc, ip, el) + + case (SOURCE_thermal_externalheat_ID) sourceType + call source_thermal_externalheat_dotState(phase,of) + + end select sourceType + + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + + enddo SourceLoop end function constitutive_collectDotState diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 86ac47aab..b7d1632f8 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -22,32 +22,6 @@ submodule(constitutive) constitutive_damage module subroutine kinematics_slipplane_opening_init end subroutine kinematics_slipplane_opening_init - module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S - end subroutine source_damage_anisoBrittle_dotState - - module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_anisoDuctile_dotState - - module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_isoDuctile_dotState - module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) @@ -134,34 +108,6 @@ module subroutine damage_init end subroutine damage_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -module procedure damage_dotState - - integer :: i - - SourceLoop: do i = 1, phase_Nsources(phase) - - sourceType: select case (phase_source(i,phase)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? - - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState ( ipc, ip, el) - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState ( ipc, ip, el) - - end select sourceType - - enddo sourceLoop - - broken_damage = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) - -end procedure damage_dotState - module procedure damage_source_getRateAndItsTangents diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index f33e7b32f..c91db80ec 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -25,67 +25,6 @@ submodule(constitutive) constitutive_plastic module subroutine plastic_nonlocal_init end subroutine plastic_nonlocal_init - module subroutine plastic_isotropic_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_dotState - - module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_dotState - - module subroutine plastic_kinehardening_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_dotState - - module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_dotState - - module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_dotState - - module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) ::& - Mp !< MandelStress - real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - F, & !< deformation gradient - Fp !< plastic deformation gradient - real(pReal), intent(in) :: & - Temperature, & !< temperature - timestep !< substepped crystallite time increment - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_dotState - module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -236,55 +175,12 @@ module subroutine plastic_init end subroutine plastic_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -module procedure plastic_dotState - - real(pReal), dimension(3,3) :: & - Mp - integer :: & - ho, & !< homogenization - tme, & !< thermal member position - instance - - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - instance = phase_plasticityInstance(phase) - - Mp = matmul(matmul(transpose(Fi),Fi),S) - - plasticityType: select case (phase_plasticity(phase)) - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (Mp,instance,of) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Mp,instance,of) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,of) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & - instance,of,ip,el) - end select plasticityType - broken_plastic = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) - -end procedure plastic_dotState - !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- -module procedure plastic_homogenizedC +module procedure constitutive_homogenizedC plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -293,7 +189,7 @@ module procedure plastic_homogenizedC homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) end select plasticityType -end procedure plastic_homogenizedC +end procedure constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- @@ -327,7 +223,7 @@ end procedure plastic_dependentState ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module procedure plastic_LpAndItsTangents +module procedure constitutive_plastic_LpAndItsTangents real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress @@ -378,7 +274,7 @@ module procedure plastic_LpAndItsTangents dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo; enddo -end procedure plastic_LpAndItsTangents +end procedure constitutive_plastic_LpAndItsTangents end submodule constitutive_plastic diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 4a34ac224..48908b997 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -13,12 +13,6 @@ submodule(constitutive) constitutive_thermal module subroutine kinematics_thermal_expansion_init end subroutine kinematics_thermal_expansion_init - module subroutine source_thermal_externalheat_dotState(phase, of) - integer, intent(in) :: & - phase, & - of - end subroutine source_thermal_externalheat_dotState - module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) @@ -61,28 +55,6 @@ module subroutine thermal_init end subroutine thermal_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -module procedure thermal_dotState - - integer :: i - - SourceLoop: do i = 1, phase_Nsources(phase) - - sourceType: select case (phase_source(i,phase)) - - case (SOURCE_thermal_externalheat_ID) sourceType - call source_thermal_externalheat_dotState(phase,of) - - end select sourceType - - enddo sourceLoop - - broken_thermal = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) - -end procedure thermal_dotState - module procedure thermal_source_getRateAndItsTangents diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fbce3ab47..eaa01270d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -554,7 +554,7 @@ subroutine crystallite_stressTangent dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS endif - call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & crystallite_S (1:3,1:3,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS @@ -915,7 +915,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) - call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ipc, ip, el) !* update current residuum and check for convergence of loop