DAMASK_EICMD/src/constitutive.f90

949 lines
48 KiB
Fortran
Raw Normal View History

!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
2013-03-06 20:11:15 +05:30
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
2020-07-15 18:05:21 +05:30
!> @brief elasticity, plasticity, damage & thermal internal microstructure state
!--------------------------------------------------------------------------------------------------
module constitutive
use prec
2019-06-15 19:10:22 +05:30
use math
use rotations
2019-06-15 19:10:22 +05:30
use IO
use config
use material
use results
use lattice
use discretization
2020-07-14 11:10:42 +05:30
use geometry_plastic_nonlocal, only: &
geometry_plastic_nonlocal_disable
2019-06-15 19:10:22 +05:30
implicit none
private
2020-08-15 19:32:10 +05:30
integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable, protected :: &
phase_elasticity !< elasticity of each phase
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & !ToDo: old intel compiler complains about protected
phase_plasticity !< plasticity of each phase
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & ! ToDo: old intel compiler complains about protected
phase_source, & !< active sources mechanisms of each phase
phase_kinematics, & !< active kinematic mechanisms of each phase
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
integer, dimension(:), allocatable, public :: & ! ToDo: old intel compiler complains about protected
phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_plasticityInstance, & !< instance of particular plasticity of each phase
phase_elasticityInstance !< instance of particular elasticity of each phase
logical, dimension(:), allocatable, public :: & ! ToDo: old intel compiler complains about protected
phase_localPlasticity !< flags phases with local constitutive law
type(tPlasticState), allocatable, dimension(:), public :: &
plasticState
type(tSourceState), allocatable, dimension(:), public :: &
sourceState
2019-06-15 19:10:22 +05:30
integer, public, protected :: &
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
interface
module subroutine plastic_init
end subroutine plastic_init
module subroutine damage_init
end subroutine damage_init
2020-07-03 20:15:11 +05:30
module subroutine thermal_init
end subroutine thermal_init
2020-08-15 19:32:10 +05:30
module function plastic_active(plastic_label) result(active_plastic)
character(len=*), intent(in) :: plastic_label
logical, dimension(:), allocatable :: active_plastic
end function plastic_active
module function source_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label
integer, intent(in) :: src_length
logical, dimension(:,:), allocatable :: active_source
end function source_active
module function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics)
character(len=*), intent(in) :: kinematics_label
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: active_kinematics
end function kinematics_active
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
2020-08-15 19:32:10 +05:30
module subroutine plastic_disloTungsten_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
2020-08-15 19:32:10 +05:30
end subroutine plastic_disloTungsten_dotState
module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
2020-07-15 18:05:21 +05:30
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) :: &
2020-07-12 20:14:26 +05:30
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) :: &
2020-07-12 20:14:26 +05:30
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) :: &
2020-07-12 20:14:26 +05:30
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
2020-07-15 18:05:21 +05:30
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(inout) :: &
phiDot, &
dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
S, & !< current 2nd Piola Kitchoff stress vector
Lp !< plastic velocity gradient
real(pReal), intent(inout) :: &
TDot, &
dTDot_dT
end subroutine constitutive_thermal_getRateAndItsTangents
module function plastic_dislotwin_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_dislotwin_homogenizedC
2020-07-10 21:59:36 +05:30
pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain)
integer, intent(in) :: &
phase, &
homog, &
offset
real(pReal), dimension(3,3) :: &
initialStrain
end function kinematics_thermal_expansion_initialStrain
2020-07-10 21:59:36 +05:30
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
integer, intent(in) :: &
instance, &
i, &
e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
2020-07-02 00:52:05 +05:30
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LiAndItsTangent
module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: &
2020-07-12 20:14:26 +05:30
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
2020-07-12 20:14:26 +05:30
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
2020-07-12 20:14:26 +05:30
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_cleavage_opening_LiAndItsTangent
module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: &
2020-07-12 20:14:26 +05:30
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
2020-07-12 20:14:26 +05:30
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
2020-07-12 20:14:26 +05:30
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_slipplane_opening_LiAndItsTangent
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
integer, intent(in) :: &
2020-07-12 20:14:26 +05:30
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
2020-07-12 20:14:26 +05:30
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
2020-07-12 20:14:26 +05:30
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
end subroutine kinematics_thermal_expansion_LiAndItsTangent
2020-07-02 00:52:05 +05:30
module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_kinehardening_deltaState
2020-07-02 00:52:05 +05:30
module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
2020-03-16 21:59:47 +05:30
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
2020-03-16 21:59:47 +05:30
instance, &
of, &
ip, &
el
end subroutine plastic_nonlocal_deltaState
2019-12-03 02:21:25 +05:30
module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
integer, intent(in) :: &
2020-07-12 20:14:26 +05:30
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
end subroutine source_damage_isoBrittle_deltaState
module subroutine plastic_results
end subroutine plastic_results
module subroutine damage_results
end subroutine damage_results
end interface
2020-07-12 16:57:28 +05:30
interface constitutive_LpAndItsTangents
2020-07-15 18:05:21 +05:30
module subroutine constitutive_plastic_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
end subroutine constitutive_plastic_LpAndItsTangents
2020-07-12 16:57:28 +05:30
end interface constitutive_LpAndItsTangents
2020-07-15 18:05:21 +05:30
2020-07-12 16:57:28 +05:30
interface constitutive_dependentState
2020-07-15 18:05:21 +05:30
module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
F, & !< elastic deformation gradient
Fp !< plastic deformation gradient
end subroutine constitutive_plastic_dependentState
2020-07-15 18:05:21 +05:30
2020-07-12 16:57:28 +05:30
end interface constitutive_dependentState
2020-07-02 02:21:21 +05:30
type :: tDebugOptions
logical :: &
basic, &
extensive, &
selective
integer :: &
element, &
ip, &
grain
end type tDebugOptions
type(tDebugOptions) :: debugConstitutive
2019-06-15 19:10:22 +05:30
public :: &
constitutive_init, &
constitutive_homogenizedC, &
constitutive_LpAndItsTangents, &
constitutive_dependentState, &
2019-06-15 19:10:22 +05:30
constitutive_LiAndItsTangents, &
constitutive_initialFi, &
constitutive_SandItsTangents, &
constitutive_collectDotState, &
2020-04-01 11:32:08 +05:30
constitutive_deltaState, &
2020-07-15 18:05:21 +05:30
constitutive_damage_getRateAndItsTangents, &
constitutive_thermal_getRateAndItsTangents, &
2020-08-15 19:32:10 +05:30
constitutive_results, &
constitutive_allocateState, &
plastic_nonlocal_updateCompatibility, &
plastic_active, &
source_active, &
kinematics_active
contains
!--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics
!--------------------------------------------------------------------------------------------------
2019-04-06 01:56:23 +05:30
subroutine constitutive_init
2019-06-15 19:10:22 +05:30
integer :: &
2020-08-15 19:32:10 +05:30
p, & !< counter in phase loop
s, & !< counter in source loop
stiffDegradationCtr
2020-07-02 02:21:21 +05:30
class (tNode), pointer :: &
2020-08-15 19:32:10 +05:30
debug_constitutive, &
phases, &
phase, &
elastic, &
stiffDegradation
2020-07-02 02:21:21 +05:30
debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic')
debugConstitutive%extensive = debug_constitutive%contains('extensive')
2020-07-03 20:15:11 +05:30
debugConstitutive%selective = debug_constitutive%contains('selective')
debugConstitutive%element = config_debug%get_asInt('element',defaultVal = 1)
debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1)
2020-07-02 02:21:21 +05:30
2020-08-15 19:32:10 +05:30
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
phases => config_material%get('phase')
2020-08-15 19:32:10 +05:30
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
allocate(phase_elasticityInstance(phases%length), source = 0)
allocate(phase_NstiffnessDegradations(phases%length),source=0)
do p = 1, phases%length
phase => phases%get(p)
elastic => phase%get('elasticity')
if(elastic%get_asString('type') == 'hooke') then
phase_elasticity(p) = ELASTICITY_HOOKE_ID
else
call IO_error(200,ext_msg=elastic%get_asString('type'))
endif
stiffDegradation => phase%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms
phase_NstiffnessDegradations(p) = stiffDegradation%length
enddo
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
source=STIFFNESS_DEGRADATION_undefined_ID)
if(maxVal(phase_NstiffnessDegradations)/=0) then
do p = 1, phases%length
phase => phases%get(p)
stiffDegradation => phase%get('stiffness_degradation',defaultVal=emptyList)
do stiffDegradationCtr = 1, stiffDegradation%length
if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID
enddo
enddo
endif
do p = 1, phases%length
phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p))
enddo
2020-07-13 18:18:23 +05:30
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
! initialize constitutive laws
call plastic_init
call damage_init
call thermal_init
write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6)
2019-06-15 19:10:22 +05:30
constitutive_source_maxSizeDotState = 0
2020-08-15 19:32:10 +05:30
PhaseLoop2:do p = 1,phases%length
!--------------------------------------------------------------------------------------------------
2020-06-26 15:14:17 +05:30
! partition and initialize state
2020-08-15 19:32:10 +05:30
plasticState(p)%partionedState0 = plasticState(p)%state0
plasticState(p)%state = plasticState(p)%partionedState0
forall(s = 1:phase_Nsources(p))
sourceState(p)%p(s)%partionedState0 = sourceState(p)%p(s)%state0
sourceState(p)%p(s)%state = sourceState(p)%p(s)%partionedState0
2019-06-15 19:10:22 +05:30
end forall
2020-07-13 18:18:23 +05:30
2020-03-16 23:58:50 +05:30
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, &
2020-08-15 19:32:10 +05:30
maxval(sourceState(p)%p%sizeDotState))
2019-06-15 19:10:22 +05:30
enddo PhaseLoop2
2020-03-16 23:58:50 +05:30
constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
end subroutine constitutive_init
2009-03-06 15:32:36 +05:30
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
module function source_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label !< name of source mechanism
2020-08-15 19:32:10 +05:30
integer, intent(in) :: src_length !< max. number of sources in system
logical, dimension(:,:), allocatable :: active_source
class(tNode), pointer :: &
phases, &
phase, &
sources, &
src
2020-08-15 19:32:10 +05:30
integer :: p,s
phases => config_material%get('phase')
2020-08-15 19:32:10 +05:30
allocate(active_source(src_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
sources => phase%get('source',defaultVal=emptyList)
do s = 1, sources%length
src => sources%get(s)
if(src%get_asString('type') == source_label) active_source(s,p) = .true.
enddo
enddo
end function source_active
!--------------------------------------------------------------------------------------------------
!> @brief checks if a kinematic mechanism is active or not
!--------------------------------------------------------------------------------------------------
module function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics)
character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism
integer, intent(in) :: kinematics_length !< max. number of kinematics in system
logical, dimension(:,:), allocatable :: active_kinematics
class(tNode), pointer :: &
phases, &
phase, &
kinematics, &
kinematics_type
2020-08-15 19:32:10 +05:30
integer :: p,k
phases => config_material%get('phase')
2020-08-15 19:32:10 +05:30
allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
kinematics => phase%get('kinematics',defaultVal=emptyList)
do k = 1, kinematics%length
kinematics_type => kinematics%get(k)
if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true.
enddo
enddo
end function kinematics_active
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
function constitutive_homogenizedC(ipc,ip,el)
2020-07-13 18:18:23 +05:30
real(pReal), dimension(6,6) :: &
constitutive_homogenizedC
integer, intent(in) :: &
2020-07-13 18:18:23 +05:30
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
case default plasticityType
constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el))
end select plasticityType
end function constitutive_homogenizedC
!--------------------------------------------------------------------------------------------------
2015-10-14 00:22:01 +05:30
!> @brief contains the constitutive equation for calculating the velocity gradient
2018-12-30 17:05:26 +05:30
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, ipc, ip, el)
2015-10-14 00:22:01 +05:30
2019-06-15 19:10:22 +05:30
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
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
2019-06-15 19:10:22 +05:30
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
2020-07-03 20:15:11 +05:30
2019-06-15 19:10:22 +05:30
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
2019-06-15 19:10:22 +05:30
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_isotropic_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
2020-07-02 00:52:05 +05:30
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
2019-06-15 19:10:22 +05:30
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
2019-06-15 19:10:22 +05:30
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
2019-06-15 19:10:22 +05:30
KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el))
kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, 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
2019-06-15 19:10:22 +05:30
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)
2019-06-15 19:10:22 +05:30
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)
2019-07-01 10:39:51 +05:30
enddo; enddo
2015-10-14 00:22:01 +05:30
2018-08-28 18:37:39 +05:30
end subroutine constitutive_LiAndItsTangents
!--------------------------------------------------------------------------------------------------
2015-10-14 00:22:01 +05:30
!> @brief collects initial intermediate deformation gradient
!--------------------------------------------------------------------------------------------------
pure function constitutive_initialFi(ipc, ip, el)
2015-10-14 00:22:01 +05:30
2019-06-15 19:10:22 +05:30
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
constitutive_initialFi !< composite initial intermediate deformation gradient
integer :: &
k !< counter in kinematics loop
integer :: &
phase, &
homog, offset
2019-06-15 19:10:22 +05:30
constitutive_initialFi = math_I3
phase = material_phaseAt(ipc,el)
2019-06-15 19:10:22 +05:30
KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption
kinematicsType: select case (phase_kinematics(k,phase))
case (KINEMATICS_thermal_expansion_ID) kinematicsType
homog = material_homogenizationAt(el)
offset = thermalMapping(homog)%p(ip,el)
constitutive_initialFi = &
constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset)
2019-06-15 19:10:22 +05:30
end select kinematicsType
enddo KinematicsLoop
2015-10-14 00:22:01 +05:30
end function constitutive_initialFi
!--------------------------------------------------------------------------------------------------
2015-10-14 00:22:01 +05:30
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic/intermediate deformation gradients depending on the selected elastic law
2018-08-28 18:24:36 +05:30
!! (so far no case switch because only Hooke is implemented)
!--------------------------------------------------------------------------------------------------
2018-08-28 18:24:36 +05:30
subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
2019-06-15 19:10:22 +05:30
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor
real(pReal), intent(out), dimension(3,3,3,3) :: &
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
2015-10-14 00:22:01 +05:30
2019-06-15 19:10:22 +05:30
call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
2015-10-14 00:22:01 +05:30
2018-08-28 18:24:36 +05:30
end subroutine constitutive_SandItsTangents
!--------------------------------------------------------------------------------------------------
2015-10-14 00:22:01 +05:30
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
2020-06-26 15:14:17 +05:30
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ipc, ip, el)
2019-06-15 19:10:22 +05:30
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
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
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer :: &
ho, & !< homogenization
d !< counter in degradation loop
integer :: &
i, j
2019-06-15 19:10:22 +05:30
ho = material_homogenizationAt(el)
C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el))
2019-06-15 19:10:22 +05:30
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el)))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2
end select degradationType
enddo DegradationLoop
2019-06-15 19:10:22 +05:30
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
do i =1, 3;do j=1,3
2019-06-15 19:10:22 +05:30
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo
2015-10-14 00:22:01 +05:30
2018-08-28 18:24:36 +05:30
end subroutine constitutive_hooke_SandItsTangents
!--------------------------------------------------------------------------------------------------
2015-10-14 00:22:01 +05:30
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
2020-04-01 16:25:49 +05:30
function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken)
2019-06-15 19:10:22 +05:30
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
2020-07-12 20:14:26 +05:30
el, & !< element
2020-04-01 16:25:49 +05:30
phase, &
of
2019-06-15 19:10:22 +05:30
real(pReal), intent(in) :: &
subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
2020-03-15 17:39:27 +05:30
FArray, & !< elastic deformation gradient
2019-06-15 19:10:22 +05:30
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)
real(pReal), dimension(3,3) :: &
Mp
integer :: &
2020-07-13 18:18:23 +05:30
ho, & !< homogenization
tme, & !< thermal member position
i, & !< counter in source loop
instance
2020-04-01 11:11:55 +05:30
logical :: broken
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
2020-08-15 19:32:10 +05:30
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
2020-08-15 19:32:10 +05:30
call plastic_dislotwin_dotState(Mp,temperature(ho)%p(tme),instance,of)
2020-08-15 19:32:10 +05:30
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
2020-08-15 19:32:10 +05:30
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
2020-08-15 19:32:10 +05:30
call source_damage_anisoBrittle_dotState(S, ipc, ip, el) ! correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType
2020-08-15 19:32:10 +05:30
call source_damage_isoDuctile_dotState(ipc, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType
2020-08-15 19:32:10 +05:30
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
2020-04-01 11:11:55 +05:30
end function constitutive_collectDotState
2009-03-06 15:32:36 +05:30
2019-06-15 19:10:22 +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
!--------------------------------------------------------------------------------------------------
2020-04-01 16:25:49 +05:30
function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken)
2015-10-14 00:22:01 +05:30
2019-06-15 19:10:22 +05:30
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
2020-04-01 16:25:49 +05:30
el, & !< element
phase, &
of
2019-06-15 19:10:22 +05:30
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: &
Mp
integer :: &
i, &
instance, &
myOffset, &
mySize
2020-04-01 11:32:08 +05:30
logical :: &
broken
2019-06-15 19:10:22 +05:30
Mp = matmul(matmul(transpose(Fi),Fi),S)
2020-04-01 16:25:49 +05:30
instance = phase_plasticityInstance(phase)
2020-04-01 11:32:08 +05:30
plasticityType: select case (phase_plasticity(phase))
2019-06-15 19:10:22 +05:30
case (PLASTICITY_KINEHARDENING_ID) plasticityType
2020-07-02 00:52:05 +05:30
call plastic_kinehardening_deltaState(Mp,instance,of)
2020-04-01 16:25:49 +05:30
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
2019-06-15 19:10:22 +05:30
case (PLASTICITY_NONLOCAL_ID) plasticityType
2020-07-02 00:52:05 +05:30
call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
2020-04-01 16:25:49 +05:30
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
case default
broken = .false.
2019-06-15 19:10:22 +05:30
end select plasticityType
2020-04-01 16:25:49 +05:30
if(.not. broken) then
select case(phase_plasticity(phase))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
myOffset = plasticState(phase)%offsetDeltaState
mySize = plasticState(phase)%sizeDeltaState
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = &
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of)
end select
endif
2020-04-01 11:32:08 +05:30
sourceLoop: do i = 1, phase_Nsources(phase)
2020-04-01 11:32:08 +05:30
sourceType: select case (phase_source(i,phase))
2019-06-15 19:10:22 +05:30
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
ipc, ip, el)
2020-04-01 16:25:49 +05:30
broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of)))
if(.not. broken) then
myOffset = sourceState(phase)%p(i)%offsetDeltaState
mySize = sourceState(phase)%p(i)%sizeDeltaState
sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = &
sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of)
endif
2019-06-15 19:10:22 +05:30
end select sourceType
2019-06-15 19:10:22 +05:30
enddo SourceLoop
2020-04-01 11:32:08 +05:30
end function constitutive_deltaState
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState(state, &
NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: &
state
integer, intent(in) :: &
NipcMyPhase, &
sizeState, &
sizeDotState, &
sizeDeltaState
state%sizeState = sizeState
state%sizeDotState = sizeDotState
state%sizeDeltaState = sizeDeltaState
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal)
end subroutine constitutive_allocateState
2018-12-05 04:25:39 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file
2018-12-05 04:25:39 +05:30
!--------------------------------------------------------------------------------------------------
subroutine constitutive_results
2019-12-19 00:35:51 +05:30
call plastic_results
call damage_results
2019-12-19 00:35:51 +05:30
2018-12-05 04:25:39 +05:30
end subroutine constitutive_results
2020-08-15 19:32:10 +05:30
end module constitutive