diff --git a/PRIVATE b/PRIVATE index dab246aea..c0716e46b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit dab246aea2cc3f29d57e0043edfebd3aedded124 +Subproject commit c0716e46b72439a49212c5e9803a5886c03b4739 diff --git a/examples/config/phase/damage/anisobrittle_cubic.yaml b/examples/config/phase/damage/anisobrittle_cubic.yaml index 410bcf2c3..3539a24d2 100644 --- a/examples/config/phase/damage/anisobrittle_cubic.yaml +++ b/examples/config/phase/damage/anisobrittle_cubic.yaml @@ -1,8 +1,11 @@ -N_cl: [3] -dot_o: 1e-3 -g_crit: [0.50e7] -q: 20 -s_crit: [0.006666] type: anisobrittle +N_cl: [3] +g_crit: [0.50e7] +s_crit: [0.006666] +dot_o: 1e-3 +q: 20 + +output: [f_phi] + D_11: 1.0 M: 0.001 diff --git a/examples/config/phase/damage/isobrittle_generic.yaml b/examples/config/phase/damage/isobrittle_generic.yaml new file mode 100644 index 000000000..b4416f995 --- /dev/null +++ b/examples/config/phase/damage/isobrittle_generic.yaml @@ -0,0 +1,9 @@ +type: isobrittle +W_crit: 1400000.0 +m: 1.0 +isoBrittle_atol: 0.01 + +output: [f_phi] + +D_11: 1.0 +M: 0.001 diff --git a/python/damask/_table.py b/python/damask/_table.py index dd6a981b7..a65971f39 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -384,7 +384,7 @@ class Table: Returns ------- - udated : damask.Table + updated : damask.Table Updated table. """ @@ -412,7 +412,7 @@ class Table: Returns ------- - udated : damask.Table + updated : damask.Table Updated table. """ @@ -435,7 +435,7 @@ class Table: Returns ------- - udated : damask.Table + updated : damask.Table Updated table. """ @@ -461,7 +461,7 @@ class Table: Returns ------- - udated : damask.Table + updated : damask.Table Updated table. """ @@ -494,7 +494,7 @@ class Table: Returns ------- - udated : damask.Table + updated : damask.Table Updated table. """ @@ -519,7 +519,7 @@ class Table: Returns ------- - udated : damask.Table + updated : damask.Table Updated table. """ diff --git a/src/phase.f90 b/src/phase.f90 index 927ba6b2b..4c9f170ee 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -67,10 +67,6 @@ module phase damageState - integer, public, protected :: & - phase_plasticity_maxSizeDotState, & - phase_source_maxSizeDotState - interface ! == cleaned:begin ================================================================================= @@ -235,11 +231,16 @@ module phase logical :: converged_ end function crystallite_stress + !ToDo: Try to merge the all stiffness functions module function phase_homogenizedC(ph,en) result(C) integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C end function phase_homogenizedC - + module function phase_damage_C(C_homogenized,ph,en) result(C) + real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized + integer, intent(in) :: ph,en + real(pReal), dimension(3,3,3,3) :: C + end function phase_damage_C module function phase_f_phi(phi,co,ce) result(f) integer, intent(in) :: ce,co @@ -372,19 +373,6 @@ subroutine phase_init call damage_init call thermal_init(phases) - - phase_source_maxSizeDotState = 0 - PhaseLoop2:do ph = 1,phases%length -!-------------------------------------------------------------------------------------------------- -! partition and initialize state - plasticState(ph)%state = plasticState(ph)%state0 - if(damageState(ph)%sizeState > 0) & - damageState(ph)%state = damageState(ph)%state0 - enddo PhaseLoop2 - - phase_source_maxSizeDotState = maxval(damageState%sizeDotState) - phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) - end subroutine phase_init @@ -540,8 +528,7 @@ subroutine crystallite_init() phases => config_material%get('phase') do ph = 1, phases%length - if (damageState(ph)%sizeState > 0) & - allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack + if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -555,7 +542,7 @@ subroutine crystallite_init() do ip = 1, size(material_phaseMemberAt,2) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) call crystallite_orientations(co,ip,el) - call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states + call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -571,9 +558,9 @@ end subroutine crystallite_init subroutine crystallite_orientations(co,ip,el) integer, intent(in) :: & - co, & !< counter in integration point component loop - ip, & !< counter in integration point loop - el !< counter in element loop + co, & !< counter in integration point component loop + ip, & !< counter in integration point loop + el !< counter in element loop call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5a1d5970b..45898f54d 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -15,13 +15,15 @@ submodule(phase) damage DAMAGE_ANISOBRITTLE_ID end enum + integer :: phase_damage_maxSizeDotState + type :: tDataContainer real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi end type tDataContainer integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & - phase_source !< active sources mechanisms of each phase + phase_damage !< active sources mechanisms of each phase type(tDataContainer), dimension(:), allocatable :: current @@ -126,17 +128,38 @@ module subroutine damage_init enddo - allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID) + allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) if (damage_active) then - where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID - where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID - where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID + where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID + where(isoductile_init() ) phase_damage = DAMAGE_ISODUCTILE_ID + where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID endif + phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) + end subroutine damage_init +!-------------------------------------------------------------------------------------------------- +!> @brief returns the degraded/modified elasticity matrix +!-------------------------------------------------------------------------------------------------- +module function phase_damage_C(C_homogenized,ph,en) result(C) + + real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized + integer, intent(in) :: ph,en + real(pReal), dimension(3,3,3,3) :: C + + damageType: select case (phase_damage(ph)) + case (DAMAGE_ISOBRITTLE_ID) damageType + C = C_homogenized * damage_phi(ph,en)**2 + case default damageType + C = C_homogenized + end select damageType + +end function phase_damage_C + + !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- @@ -155,7 +178,7 @@ module function phase_f_phi(phi,co,ce) result(f) ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) - select case(phase_source(ph)) + select case(phase_damage(ph)) case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID) f = 1.0_pReal & - phi*damageState(ph)%state(1,en) @@ -187,9 +210,9 @@ module function integrateDamageState(dt,co,ip,el) result(broken) size_so real(pReal) :: & zeta - real(pReal), dimension(phase_source_maxSizeDotState) :: & + real(pReal), dimension(phase_damage_maxSizeDotState) :: & r ! state residuum - real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState + real(pReal), dimension(phase_damage_maxSizeDotState,2) :: source_dotState logical :: & converged_ @@ -275,10 +298,10 @@ module subroutine damage_results(group,ph) integer, intent(in) :: ph - if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) & + if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'damage')) - sourceType: select case (phase_source(ph)) + sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_results(ph,group//'damage/') @@ -309,7 +332,7 @@ function phase_damage_collectDotState(ph,me) result(broken) if (damageState(ph)%sizeState > 0) then - sourceType: select case (phase_source(ph)) + sourceType: select case (phase_damage(ph)) case (DAMAGE_ISODUCTILE_ID) sourceType call isoductile_dotState(ph,me) @@ -376,7 +399,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) if (damageState(ph)%sizeState == 0) return - sourceType: select case (phase_source(ph)) + sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index e88a87352..64b49b1bc 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -66,14 +66,14 @@ module function isobrittle_init() result(mySources) Nmembers = count(material_phaseID==ph) call phase_allocateState(damageState(ph),Nmembers,1,1,1) - damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + damageState(ph)%atol = src%get_asFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoBrittle)') + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') endif enddo diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 997b948fe..42b7db434 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -70,14 +70,14 @@ module function isoductile_init() result(mySources) Nmembers=count(material_phaseID==ph) call phase_allocateState(damageState(ph),Nmembers,1,1,0) - damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + damageState(ph)%atol = src%get_asFloat('isoductile_atol',defaultVal=1.0e-3_pReal) if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoDuctile)') + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoductile)') endif enddo diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 29e9e8ada..be0bca892 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -41,6 +41,7 @@ submodule(phase) mechanical integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & phase_plasticity !< plasticity of each phase + integer :: phase_plasticity_maxSizeDotState interface @@ -68,7 +69,7 @@ submodule(phase) mechanical dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient end subroutine phase_hooke_SandItsTangents - + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -294,6 +295,11 @@ module subroutine mechanical_init(materials,phases) call plastic_init() + do ph = 1,phases%length + plasticState(ph)%state0 = plasticState(ph)%state + enddo + phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) + num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 9c18bfa7f..e25e5e8d2 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -2,18 +2,11 @@ submodule(phase:mechanical) elastic enum, bind(c); enumerator :: & ELASTICITY_UNDEFINED_ID, & - ELASTICITY_HOOKE_ID, & - STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID + ELASTICITY_HOOKE_ID end enum - - integer, dimension(:), allocatable :: & - phase_NstiffnessDegradations integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase - integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & - phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase contains @@ -24,18 +17,15 @@ module subroutine elastic_init(phases) phases integer :: & - ph, & - stiffDegradationCtr + ph class(tNode), pointer :: & phase, & mech, & - elastic, & - stiffDegradation + elastic print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>' allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) - allocate(phase_NstiffnessDegradations(phases%length),source=0) do ph = 1, phases%length phase => phases%get(ph) @@ -46,25 +36,8 @@ module subroutine elastic_init(phases) else call IO_error(200,ext_msg=elastic%get_asString('type')) endif - stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms - phase_NstiffnessDegradations(ph) = stiffDegradation%length enddo - allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & - source=STIFFNESS_DEGRADATION_undefined_ID) - - if(maxVal(phase_NstiffnessDegradations)/=0) then - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanical') - stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) - do stiffDegradationCtr = 1, stiffDegradation%length - if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & - phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID - enddo - enddo - endif - end subroutine elastic_init @@ -94,13 +67,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j C = math_66toSym3333(phase_homogenizedC(ph,en)) - - DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) - degradationType: select case(phase_stiffnessDegradation(d,ph)) - case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage_phi(ph,en)**2 - end select degradationType - enddo DegradationLoop + C = phase_damage_C(C,ph,en) E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 4a5c6c6d8..0372d7f63 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -254,8 +254,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) - plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally - end associate !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 151250f14..fdfd630c6 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -467,8 +467,6 @@ module function plastic_dislotwin_init() result(myPlasticity) allocate(dst%tau_r_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) - plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally - end associate !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 5ab73895f..0e7920f1d 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -140,8 +140,6 @@ module function plastic_isotropic_init() result(myPlasticity) ! global alias plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:) - plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally - end associate !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 18f467617..2e9d713f2 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -213,8 +213,6 @@ module function plastic_kinehardening_init() result(myPlasticity) stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:) dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) - plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally - end associate !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index b2db9e346..13335df4c 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -265,8 +265,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' - plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally - end associate !--------------------------------------------------------------------------------------------------