!---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(phase) damage type :: tDamageParameters real(pReal) :: mu = 0.0_pReal !< viscosity real(pReal), dimension(3,3) :: D = 0.0_pReal !< conductivity/diffusivity end type tDamageParameters enum, bind(c); enumerator :: & DAMAGE_UNDEFINED_ID, & DAMAGE_ISOBRITTLE_ID, & DAMAGE_ANISOBRITTLE_ID end enum integer :: phase_damage_maxSizeDotState type :: tDataContainer real(pReal), dimension(:), allocatable :: phi end type tDataContainer integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & phase_damage !< active sources mechanisms of each phase type(tDataContainer), dimension(:), allocatable :: current type(tDamageParameters), dimension(:), allocatable :: param interface module function anisobrittle_init() result(mySources) logical, dimension(:), allocatable :: mySources end function anisobrittle_init module function isobrittle_init() result(mySources) logical, dimension(:), allocatable :: mySources end function isobrittle_init module subroutine isobrittle_deltaState(C, Fe, ph, en) integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & C end subroutine isobrittle_deltaState module subroutine anisobrittle_dotState(S, ph, en) integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & S end subroutine anisobrittle_dotState module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group end subroutine anisobrittle_results module subroutine isobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group end subroutine isobrittle_results end interface contains !---------------------------------------------------------------------------------------------- !< @brief initialize damage sources and kinematics mechanism !---------------------------------------------------------------------------------------------- module subroutine damage_init integer :: & ph, & Nmembers class(tNode), pointer :: & phases, & phase, & sources, & source logical:: damage_active print'(/,a)', ' <<<+- phase:damage init -+>>>' phases => config_material%get('phase') allocate(current(phases%length)) allocate(damageState (phases%length)) allocate(param(phases%length)) damage_active = .false. do ph = 1,phases%length Nmembers = count(material_phaseID == ph) allocate(current(ph)%phi(Nmembers),source=1.0_pReal) phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) if (sources%length > 1) error stop if (sources%length == 1) then damage_active = .true. source => sources%get(1) param(ph)%mu = source%get_asFloat('mu') param(ph)%D(1,1) = source%get_asFloat('D_11') if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33') param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph)) endif enddo allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) if (damage_active) then where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID endif phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & ip, & el logical :: converged_ integer :: & ph, en ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) converged_ = .not. integrateDamageState(Delta_t,ph,en) end function phase_damage_constitutive !-------------------------------------------------------------------------------------------------- !> @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 Restore data after homog cutback. !-------------------------------------------------------------------------------------------------- module subroutine damage_restore(ce) integer, intent(in) :: ce integer :: & co do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) if (damageState(material_phaseID(co,ce))%sizeState > 0) & damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) enddo end subroutine damage_restore !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- module function phase_f_phi(phi,co,ce) result(f) integer, intent(in) :: ce,co real(pReal), intent(in) :: & phi !< damage parameter real(pReal) :: & f integer :: & ph, & en ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) select case(phase_damage(ph)) case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) f = 1.0_pReal & - phi*damageState(ph)%state(1,en) case default f = 0.0_pReal end select end function phase_f_phi !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- function integrateDamageState(Delta_t,ph,en) result(broken) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & ph, & en logical :: broken integer :: & NiterationState, & !< number of iterations in state loop size_so real(pReal) :: & zeta real(pReal), dimension(phase_damage_maxSizeDotState) :: & r ! state residuum real(pReal), dimension(phase_damage_maxSizeDotState,2) :: source_dotState logical :: & converged_ if (damageState(ph)%sizeState == 0) then broken = .false. return endif converged_ = .true. broken = phase_damage_collectDotState(ph,en) if(broken) return size_so = damageState(ph)%sizeDotState damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & + damageState(ph)%dotState(1:size_so,en) * Delta_t source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) broken = phase_damage_collectDotState(ph,en) if(broken) exit iteration zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) damageState(ph)%dotState(:,en) = damageState(ph)%dotState(:,en) * zeta & + source_dotState(1:size_so,1)* (1.0_pReal - zeta) r(1:size_so) = damageState(ph)%state (1:size_so,en) & - damageState(ph)%State0 (1:size_so,en) & - damageState(ph)%dotState(1:size_so,en) * Delta_t damageState(ph)%state(1:size_so,en) = damageState(ph)%state(1:size_so,en) - r(1:size_so) converged_ = converged_ .and. converged(r(1:size_so), & damageState(ph)%state(1:size_so,en), & damageState(ph)%atol(1:size_so)) if(converged_) then broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) exit iteration endif enddo iteration broken = broken .or. .not. converged_ contains !-------------------------------------------------------------------------------------------------- !> @brief calculate the damping for correction of state and dot state !-------------------------------------------------------------------------------------------------- real(pReal) pure function damper(omega_0,omega_1,omega_2) real(pReal), dimension(:), intent(in) :: & omega_0, omega_1, omega_2 real(pReal) :: dot_prod12, dot_prod22 dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2) dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2) if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) else damper = 1.0_pReal endif end function damper end function integrateDamageState !---------------------------------------------------------------------------------------------- !< @brief writes damage sources results to HDF5 output file !---------------------------------------------------------------------------------------------- module subroutine damage_results(group,ph) character(len=*), intent(in) :: group integer, intent(in) :: ph if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'damage')) sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_results(ph,group//'damage/') case (DAMAGE_ANISOBRITTLE_ID) sourceType call anisobrittle_results(ph,group//'damage/') end select sourceType end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- function phase_damage_collectDotState(ph,en) result(broken) integer, intent(in) :: & ph, & en !< counter in source loop logical :: broken broken = .false. if (damageState(ph)%sizeState > 0) then sourceType: select case (phase_damage(ph)) case (DAMAGE_ANISOBRITTLE_ID) sourceType call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! correct stress? end select sourceType broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) endif end function phase_damage_collectDotState !-------------------------------------------------------------------------------------------------- !> @brief Damage viscosity. !-------------------------------------------------------------------------------------------------- module function phase_mu_phi(co,ce) result(mu) integer, intent(in) :: co, ce real(pReal) :: mu mu = param(material_phaseID(co,ce))%mu end function phase_mu_phi !-------------------------------------------------------------------------------------------------- !> @brief Damage conductivity/diffusivity in reference configuration. !-------------------------------------------------------------------------------------------------- module function phase_K_phi(co,ce) result(K) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: K real(pReal), parameter :: l = 1.0_pReal K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) & * l**2.0_pReal end function phase_K_phi !-------------------------------------------------------------------------------------------------- !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- function phase_damage_deltaState(Fe, ph, en) result(broken) integer, intent(in) :: & ph, & en real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & myOffset, & mySize logical :: & broken broken = .false. if (damageState(ph)%sizeState == 0) return sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if(.not. broken) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en) endif end select sourceType end function phase_damage_deltaState !-------------------------------------------------------------------------------------------------- !> @brief checks if a source mechanism is active or not !-------------------------------------------------------------------------------------------------- function source_active(source_label) result(active_source) character(len=*), intent(in) :: source_label !< name of source mechanism logical, dimension(:), allocatable :: active_source class(tNode), pointer :: & phases, & phase, & sources, & src integer :: ph phases => config_material%get('phase') allocate(active_source(phases%length)) do ph = 1, phases%length phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) src => sources%get(1) active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label enddo end function source_active !---------------------------------------------------------------------------------------------- !< @brief Set damage parameter !---------------------------------------------------------------------------------------------- module subroutine phase_set_phi(phi,co,ce) real(pReal), intent(in) :: phi integer, intent(in) :: ce, co current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi end subroutine phase_set_phi module function damage_phi(ph,en) result(phi) integer, intent(in) :: ph, en real(pReal) :: phi phi = current(ph)%phi(en) end function damage_phi !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. !-------------------------------------------------------------------------------------------------- module subroutine damage_forward() integer :: ph do ph = 1, size(damageState) if (damageState(ph)%sizeState > 0) & damageState(ph)%state0 = damageState(ph)%state enddo end subroutine damage_forward end submodule damage