DAMASK_EICMD/src/phase_damage.f90

507 lines
17 KiB
Fortran
Raw Normal View History

2020-07-15 18:05:21 +05:30
!----------------------------------------------------------------------------------------------------
2020-12-22 16:50:00 +05:30
!> @brief internal microstructure state for all damage sources and kinematics constitutive models
2020-07-15 18:05:21 +05:30
!----------------------------------------------------------------------------------------------------
submodule(phase) damagee
2021-01-08 04:20:06 +05:30
enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, &
DAMAGE_ISODUCTILE_ID, &
DAMAGE_ANISOBRITTLE_ID, &
DAMAGE_ANISODUCTILE_ID
end enum
2021-01-21 01:24:31 +05:30
type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
end type tDataContainer
2021-02-13 17:14:47 +05:30
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: &
2021-01-08 04:20:06 +05:30
phase_source !< active sources mechanisms of each phase
2021-02-13 23:11:30 +05:30
integer, dimension(:), allocatable :: &
phase_Nsources
2021-01-21 01:24:31 +05:30
type(tDataContainer), dimension(:), allocatable :: current
interface
2021-02-13 17:14:47 +05:30
module function anisobrittle_init() result(mySources)
logical, dimension(:), allocatable :: mySources
2021-01-26 04:08:32 +05:30
end function anisobrittle_init
2021-02-13 17:14:47 +05:30
module function anisoductile_init() result(mySources)
logical, dimension(:), allocatable :: mySources
2021-01-26 04:08:32 +05:30
end function anisoductile_init
2021-02-13 17:14:47 +05:30
module function isobrittle_init() result(mySources)
logical, dimension(:), allocatable :: mySources
2021-01-26 04:08:32 +05:30
end function isobrittle_init
2021-02-13 17:14:47 +05:30
module function isoductile_init() result(mySources)
logical, dimension(:), allocatable :: mySources
2021-01-26 04:08:32 +05:30
end function isoductile_init
2021-02-13 15:31:08 +05:30
module subroutine isobrittle_deltaState(C, Fe, ph, me)
2021-01-19 14:55:52 +05:30
integer, intent(in) :: ph,me
real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
2021-02-13 15:31:08 +05:30
end subroutine isobrittle_deltaState
2021-01-19 14:55:52 +05:30
2020-07-10 18:29:07 +05:30
2021-02-13 15:31:08 +05:30
module subroutine anisobrittle_dotState(S, ph, me)
integer, intent(in) :: ph,me
2021-01-19 15:00:10 +05:30
real(pReal), intent(in), dimension(3,3) :: &
S
2021-01-26 04:08:32 +05:30
end subroutine anisobrittle_dotState
2021-01-19 15:00:10 +05:30
2021-02-13 15:31:08 +05:30
module subroutine anisoductile_dotState(ph,me)
integer, intent(in) :: ph,me
2021-01-26 04:08:32 +05:30
end subroutine anisoductile_dotState
2021-01-19 15:00:10 +05:30
2021-02-13 15:31:08 +05:30
module subroutine isoductile_dotState(ph,me)
integer, intent(in) :: ph,me
2021-01-26 04:08:32 +05:30
end subroutine isoductile_dotState
2021-02-13 15:31:08 +05:30
module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
integer, intent(in) :: ph,me
2021-01-26 04:08:32 +05:30
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
2021-02-13 15:31:08 +05:30
end subroutine anisobrittle_getRateAndItsTangent
2021-01-26 04:08:32 +05:30
2021-02-13 15:31:08 +05:30
module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: ph,me
2021-01-26 04:08:32 +05:30
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
2021-02-13 15:31:08 +05:30
end subroutine anisoductile_getRateAndItsTangent
2021-01-26 04:08:32 +05:30
2021-02-13 15:31:08 +05:30
module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: ph,me
2021-01-26 04:08:32 +05:30
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
2021-02-13 15:31:08 +05:30
end subroutine isobrittle_getRateAndItsTangent
2021-01-26 04:08:32 +05:30
2021-02-13 15:31:08 +05:30
module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: ph,me
2021-01-26 04:08:32 +05:30
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
2021-02-13 15:31:08 +05:30
end subroutine isoductile_getRateAndItsTangent
2021-01-26 04:08:32 +05:30
module subroutine anisobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine anisobrittle_results
module subroutine anisoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine anisoductile_results
module subroutine isobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine isobrittle_results
module subroutine isoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine isoductile_results
end interface
contains
2020-07-12 20:14:26 +05:30
!----------------------------------------------------------------------------------------------
!< @brief initialize damage sources and kinematics mechanism
!----------------------------------------------------------------------------------------------
module subroutine damage_init
2020-08-15 19:32:10 +05:30
integer :: &
2021-01-21 01:24:31 +05:30
ph, & !< counter in phase loop
2021-03-05 01:46:36 +05:30
Nmembers
2020-08-15 19:32:10 +05:30
class(tNode), pointer :: &
phases, &
phase, &
sources
2020-08-15 19:32:10 +05:30
2021-01-27 15:14:03 +05:30
print'(/,a)', ' <<<+- phase:damage init -+>>>'
phases => config_material%get('phase')
2020-08-15 19:32:10 +05:30
2021-01-21 01:24:31 +05:30
allocate(current(phases%length))
2021-01-08 11:40:38 +05:30
allocate(damageState (phases%length))
2021-01-26 04:08:32 +05:30
allocate(phase_Nsources(phases%length),source = 0)
2020-08-15 19:32:10 +05:30
do ph = 1,phases%length
2021-01-21 01:24:31 +05:30
2021-03-05 01:46:36 +05:30
Nmembers = count(material_phaseAt2 == ph)
2021-01-21 01:24:31 +05:30
2021-03-05 01:46:36 +05:30
allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
2021-01-21 01:24:31 +05:30
2020-08-15 19:32:10 +05:30
phase => phases%get(ph)
sources => phase%get('damage',defaultVal=emptyList)
if (sources%length > 1) error stop
2021-02-13 23:11:30 +05:30
phase_Nsources(ph) = sources%length
2020-08-15 19:32:10 +05:30
enddo
2021-02-13 17:14:47 +05:30
allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID)
2020-08-15 19:32:10 +05:30
! initialize source mechanisms
2020-08-15 19:32:10 +05:30
if(maxval(phase_Nsources) /= 0) then
2021-02-13 17:14:47 +05:30
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(anisoductile_init()) phase_source = DAMAGE_ANISODUCTILE_ID
2020-08-15 19:32:10 +05:30
endif
end subroutine damage_init
2020-07-12 20:14:26 +05:30
!----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force
!----------------------------------------------------------------------------------------------
module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
2020-07-12 20:14:26 +05:30
integer, intent(in) :: ce
2020-07-12 20:14:26 +05:30
real(pReal), intent(in) :: &
2020-12-22 16:50:00 +05:30
phi !< damage parameter
2020-07-12 20:14:26 +05:30
real(pReal), intent(inout) :: &
phiDot, &
dPhiDot_dPhi
real(pReal) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
2021-01-26 04:08:32 +05:30
ph, &
co, &
me
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
2020-12-22 16:50:00 +05:30
2021-03-05 01:46:36 +05:30
do co = 1, homogenization_Nmembers(material_homogenizationAt2(ce))
ph = material_phaseAt2(co,ce)
me = material_phasememberAt2(co,ce)
2021-02-13 23:11:30 +05:30
2021-02-13 17:14:47 +05:30
select case(phase_source(ph))
2021-01-08 04:20:06 +05:30
case (DAMAGE_ISOBRITTLE_ID)
2021-02-13 15:31:08 +05:30
call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
2021-01-08 04:20:06 +05:30
case (DAMAGE_ISODUCTILE_ID)
2021-02-13 15:31:08 +05:30
call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
2021-01-08 04:20:06 +05:30
case (DAMAGE_ANISOBRITTLE_ID)
2021-02-13 15:31:08 +05:30
call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
2021-01-08 04:20:06 +05:30
case (DAMAGE_ANISODUCTILE_ID)
2021-02-13 15:31:08 +05:30
call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
2021-02-09 03:51:53 +05:30
end subroutine phase_damage_getRateAndItsTangents
2021-01-08 04:02:54 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
module function integrateDamageState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: broken
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
size_so
real(pReal) :: &
zeta
2021-02-09 03:51:53 +05:30
real(pReal), dimension(phase_source_maxSizeDotState) :: &
2021-01-08 04:02:54 +05:30
r ! state residuum
real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState
2021-01-08 04:02:54 +05:30
logical :: &
converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
2021-02-13 23:11:30 +05:30
if (damageState(ph)%sizeState == 0) then
2021-02-13 18:45:41 +05:30
broken = .false.
return
endif
2021-01-08 04:02:54 +05:30
converged_ = .true.
2021-02-13 15:31:08 +05:30
broken = phase_damage_collectDotState(ph,me)
2021-01-08 04:02:54 +05:30
if(broken) return
size_so = damageState(ph)%sizeDotState
damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) &
2021-02-13 18:45:41 +05:30
+ damageState(ph)%dotState (1:size_so,me) * dt
source_dotState(1:size_so,2) = 0.0_pReal
2021-01-08 04:02:54 +05:30
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(:,me)
2021-01-08 04:02:54 +05:30
2021-02-13 15:31:08 +05:30
broken = phase_damage_collectDotState(ph,me)
2021-01-08 04:02:54 +05:30
if(broken) exit iteration
zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta &
+ source_dotState(1:size_so,1)* (1.0_pReal - zeta)
r(1:size_so) = damageState(ph)%state (1:size_so,me) &
- damageState(ph)%subState0(1:size_so,me) &
- damageState(ph)%dotState (1:size_so,me) * dt
damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so)
converged_ = converged_ .and. converged(r(1:size_so), &
damageState(ph)%state(1:size_so,me), &
damageState(ph)%atol(1:size_so))
2021-01-08 04:02:54 +05:30
if(converged_) then
2021-02-09 03:51:53 +05:30
broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
2021-01-08 04:02:54 +05:30
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(current,previous,previous2)
real(pReal), dimension(:), intent(in) ::&
current, previous, previous2
real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(current - previous, previous - previous2)
dot_prod22 = dot_product(previous - previous2, previous - previous2)
if ((dot_product(current,previous) < 0.0_pReal .or. 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
2020-07-12 20:14:26 +05:30
!----------------------------------------------------------------------------------------------
2020-08-15 19:32:10 +05:30
!< @brief writes damage sources results to HDF5 output file
2020-07-12 20:14:26 +05:30
!----------------------------------------------------------------------------------------------
2020-12-22 16:50:00 +05:30
module subroutine damage_results(group,ph)
2020-12-22 16:50:00 +05:30
character(len=*), intent(in) :: group
integer, intent(in) :: ph
2020-12-22 16:50:00 +05:30
integer :: so
2020-08-15 19:32:10 +05:30
2020-12-22 16:50:00 +05:30
sourceLoop: do so = 1, phase_Nsources(ph)
2021-02-13 17:14:47 +05:30
if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) &
2020-12-23 02:52:43 +05:30
call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage'
2021-02-13 17:14:47 +05:30
sourceType: select case (phase_source(ph))
2021-01-08 04:20:06 +05:30
case (DAMAGE_ISOBRITTLE_ID) sourceType
2021-01-26 04:08:32 +05:30
call isobrittle_results(ph,group//'sources/')
2020-12-22 16:50:00 +05:30
2021-01-08 04:20:06 +05:30
case (DAMAGE_ISODUCTILE_ID) sourceType
2021-01-26 04:08:32 +05:30
call isoductile_results(ph,group//'sources/')
2020-12-22 16:50:00 +05:30
2021-01-08 04:20:06 +05:30
case (DAMAGE_ANISOBRITTLE_ID) sourceType
2021-01-26 04:08:32 +05:30
call anisobrittle_results(ph,group//'sources/')
2020-12-22 16:50:00 +05:30
2021-01-08 04:20:06 +05:30
case (DAMAGE_ANISODUCTILE_ID) sourceType
2021-01-26 04:08:32 +05:30
call anisoductile_results(ph,group//'sources/')
2020-12-22 16:50:00 +05:30
end select sourceType
enddo SourceLoop
end subroutine damage_results
2021-01-08 04:02:54 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
2021-02-13 15:31:08 +05:30
function phase_damage_collectDotState(ph,me) result(broken)
2021-01-08 04:02:54 +05:30
integer, intent(in) :: &
ph, &
2021-02-13 15:31:08 +05:30
me !< counter in source loop
2021-01-08 04:02:54 +05:30
logical :: broken
broken = .false.
2021-02-13 23:11:30 +05:30
if (damageState(ph)%sizeState > 0) then
2021-01-08 04:02:54 +05:30
2021-02-13 17:14:47 +05:30
sourceType: select case (phase_source(ph))
2021-01-08 04:02:54 +05:30
2021-01-08 04:20:06 +05:30
case (DAMAGE_ISODUCTILE_ID) sourceType
2021-02-13 15:31:08 +05:30
call isoductile_dotState(ph,me)
2021-01-08 04:02:54 +05:30
2021-01-08 04:20:06 +05:30
case (DAMAGE_ANISODUCTILE_ID) sourceType
2021-02-13 15:31:08 +05:30
call anisoductile_dotState(ph,me)
2021-01-08 04:02:54 +05:30
case (DAMAGE_ANISOBRITTLE_ID) sourceType
2021-02-13 15:31:08 +05:30
call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?
2021-01-08 04:02:54 +05:30
end select sourceType
broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me)))
2021-01-08 04:02:54 +05:30
2021-02-13 15:31:08 +05:30
endif
2021-01-08 04:02:54 +05:30
2021-02-09 03:51:53 +05:30
end function phase_damage_collectDotState
2021-01-08 04:02:54 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
function phase_damage_deltaState(Fe, ph, me) result(broken)
2021-01-08 04:02:54 +05:30
integer, intent(in) :: &
ph, &
2021-01-19 14:51:51 +05:30
me
2021-01-08 04:02:54 +05:30
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
integer :: &
myOffset, &
mySize
logical :: &
broken
broken = .false.
2021-02-13 23:11:30 +05:30
if (damageState(ph)%sizeState == 0) return
2021-01-08 04:02:54 +05:30
2021-02-13 17:14:47 +05:30
sourceType: select case (phase_source(ph))
2021-01-08 04:02:54 +05:30
2021-01-08 04:20:06 +05:30
case (DAMAGE_ISOBRITTLE_ID) sourceType
2021-02-13 15:31:08 +05:30
call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
2021-01-08 04:02:54 +05:30
if(.not. broken) then
myOffset = damageState(ph)%offsetDeltaState
mySize = damageState(ph)%sizeDeltaState
damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = &
damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me)
2021-01-08 04:02:54 +05:30
endif
end select sourceType
2021-02-09 03:51:53 +05:30
end function phase_damage_deltaState
2021-01-08 04:02:54 +05:30
2021-01-08 12:07:51 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
2021-02-13 17:14:47 +05:30
function source_active(source_label) result(active_source)
2021-01-08 12:07:51 +05:30
character(len=*), intent(in) :: source_label !< name of source mechanism
2021-02-13 17:14:47 +05:30
logical, dimension(:), allocatable :: active_source
2021-01-08 12:07:51 +05:30
class(tNode), pointer :: &
phases, &
phase, &
sources, &
src
2021-02-13 17:14:47 +05:30
integer :: ph
2021-01-08 12:07:51 +05:30
phases => config_material%get('phase')
2021-02-13 17:14:47 +05:30
allocate(active_source(phases%length))
do ph = 1, phases%length
phase => phases%get(ph)
sources => phase%get('damage',defaultVal=emptyList)
2021-02-13 17:14:47 +05:30
src => sources%get(1)
2021-02-13 23:11:30 +05:30
active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label
2021-01-08 12:07:51 +05:30
enddo
end function source_active
2021-01-21 01:24:31 +05:30
!----------------------------------------------------------------------------------------------
!< @brief Set damage parameter
!----------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module subroutine phase_damage_set_phi(phi,co,ce)
2021-01-21 01:24:31 +05:30
real(pReal), intent(in) :: phi
integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi
2021-02-09 03:51:53 +05:30
end subroutine phase_damage_set_phi
2021-01-21 01:24:31 +05:30
2021-02-09 03:51:53 +05:30
module function phase_damage_get_phi(co,ip,el) result(phi)
2021-01-21 01:24:31 +05:30
integer, intent(in) :: co, ip, el
real(pReal) :: phi
phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
2021-02-09 03:51:53 +05:30
end function phase_damage_get_phi
2021-01-21 01:24:31 +05:30
2021-02-13 11:45:57 +05:30
module function damage_phi(ph,me) result(phi)
integer, intent(in) :: ph, me
real(pReal) :: phi
phi = current(ph)%phi(me)
end function damage_phi
end submodule damagee