223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

sourceState is now damage state

This commit is contained in:
Martin Diehl 2021-01-07 23:32:54 +01:00
parent 27f4e4ce2a
commit 1df409376c
3 changed files with 286 additions and 272 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

283
src/constitutive.f90
View File

@ -235,14 +235,22 @@ module constitutive
! == cleaned:end ===================================================================================
module function integrateThermalState(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
end function
end function integrateThermalState
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
end function integrateDamageState
module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt
@ -395,7 +403,6 @@ module constitutive
public :: &
constitutive_init, &
constitutive_homogenizedC, &
constitutive_LiAndItsTangents, &
constitutive_damage_getRateAndItsTangents, &
constitutive_thermal_getRateAndItsTangents, &
constitutive_results, &
@ -413,7 +420,8 @@ module constitutive
crystallite_push33ToRef, &
constitutive_restartWrite, &
constitutive_restartRead, &
integrateSourceState, &
integrateThermalState, &
integrateDamageState, &
constitutive_thermal_setT, &
constitutive_mech_getP, &
constitutive_mech_setF, &
@ -555,173 +563,6 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki
end function kinematics_active
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Li !< intermediate velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dS, & !< derivative of Li with respect to S
dLi_dFi
real(pReal), dimension(3,3) :: &
my_Li, & !< intermediate velocity gradient
FiInv, &
temp_33
real(pReal), dimension(3,3,3,3) :: &
my_dLi_dS
real(pReal) :: &
detFi
integer :: &
k, i, j, &
instance, of
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_isotropic_ID) plasticityType
of = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el))
kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el)
case default kinematicsType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select kinematicsType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
enddo KinematicsLoop
FiInv = math_inv33(Fi)
detFi = math_det33(Fi)
Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = matmul(FiInv,Li)
do i = 1,3; do j = 1,3
dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
enddo; enddo
end subroutine constitutive_LiAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
integer :: &
so !< counter in source loop
logical :: broken
broken = .false.
SourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),&
co, ip, el) ! correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState(co, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_dotState(co, ip, el)
end select sourceType
broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of)))
enddo SourceLoop
end function constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @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 constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
integer :: &
so, &
myOffset, &
mySize
logical :: &
broken
broken = .false.
sourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, &
co, ip, el)
broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of)))
if(.not. broken) then
myOffset = sourceState(ph)%p(so)%offsetDeltaState
mySize = sourceState(ph)%p(so)%sizeDeltaState
sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = &
sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of)
endif
end select sourceType
enddo SourceLoop
end function constitutive_damage_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
@ -1030,107 +871,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
end function crystallite_push33ToRef
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
function integrateSourceState(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
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
so
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: &
broken, converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) return
do so = 1, phase_Nsources(ph)
size_so(so) = sourceState(ph)%p(so)%sizeDotState
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
+ sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
do so = 1, phase_Nsources(ph)
if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me)
enddo
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration
do so = 1, phase_Nsources(ph)
zeta = damper(sourceState(ph)%p(so)%dotState(:,me), &
source_dotState(1:size_so(so),1,so),&
source_dotState(1:size_so(so),2,so))
sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) &
- sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
- sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
sourceState(ph)%p(so)%state(1:size_so(so),me), &
sourceState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me)
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 integrateSourceState
!--------------------------------------------------------------------------------------------------

193
src/constitutive_damage.f90
View File

@ -214,6 +214,111 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @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, &
so
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: &
converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) return
do so = 1, phase_Nsources(ph)
size_so(so) = sourceState(ph)%p(so)%sizeDotState
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
+ sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
do so = 1, phase_Nsources(ph)
if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me)
enddo
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration
do so = 1, phase_Nsources(ph)
zeta = damper(sourceState(ph)%p(so)%dotState(:,me), &
source_dotState(1:size_so(so),1,so),&
source_dotState(1:size_so(so),2,so))
sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) &
- sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
- sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
sourceState(ph)%p(so)%state(1:size_so(so),me), &
sourceState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me)
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
!----------------------------------------------------------------------------------------------
!< @brief writes damage sources results to HDF5 output file
!----------------------------------------------------------------------------------------------
@ -250,4 +355,92 @@ module subroutine damage_results(group,ph)
end subroutine damage_results
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
integer :: &
so !< counter in source loop
logical :: broken
broken = .false.
SourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),&
co, ip, el) ! correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState(co, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_dotState(co, ip, el)
end select sourceType
broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of)))
enddo SourceLoop
end function constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @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 constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
integer :: &
so, &
myOffset, &
mySize
logical :: &
broken
broken = .false.
sourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, &
co, ip, el)
broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of)))
if(.not. broken) then
myOffset = sourceState(ph)%p(so)%offsetDeltaState
mySize = sourceState(ph)%p(so)%sizeDeltaState
sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = &
sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of)
endif
end select sourceType
enddo SourceLoop
end function constitutive_damage_deltaState
end submodule constitutive_damage

82
src/constitutive_mech.f90
View File

@ -1653,7 +1653,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el)
endif
@ -1938,5 +1938,85 @@ module subroutine constitutive_mech_setF(F,co,ip,el)
end subroutine constitutive_mech_setF
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Li !< intermediate velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dS, & !< derivative of Li with respect to S
dLi_dFi
real(pReal), dimension(3,3) :: &
my_Li, & !< intermediate velocity gradient
FiInv, &
temp_33
real(pReal), dimension(3,3,3,3) :: &
my_dLi_dS
real(pReal) :: &
detFi
integer :: &
k, i, j, &
instance, of
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_isotropic_ID) plasticityType
of = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el))
kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el)
case default kinematicsType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select kinematicsType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
enddo KinematicsLoop
FiInv = math_inv33(Fi)
detFi = math_det33(Fi)
Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = matmul(FiInv,Li)
do i = 1,3; do j = 1,3
dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
enddo; enddo
end subroutine constitutive_LiAndItsTangents
end submodule constitutive_mech