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

Merge branch 'rename-and-restructure' into 'development' 

Rename and restructure

See merge request damask/DAMASK!329
This commit is contained in:
Franz Roters 2021-01-28 16:41:17 +01:00
parent 8ec385140c 0d3b9b9f2a
commit 85cfa0bab0
35 changed files with 1584 additions and 1544 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CPFEM.f90
View File

@ -19,7 +19,7 @@ module CPFEM
use HDF5_utilities
use results
use lattice
use constitutive
use phase
implicit none
private

2
src/CPFEM2.f90
View File

@ -19,7 +19,7 @@ module CPFEM2
use discretization
use HDF5_utilities
use homogenization
use constitutive
use phase
#if defined(Mesh)
use FEM_quadrature
use discretization_mesh

53
src/commercialFEM_fileList.f90
View File

@ -7,7 +7,6 @@
#include "IO.f90"
#include "YAML_types.f90"
#include "YAML_parse.f90"
#include "future.f90"
#include "config.f90"
#include "LAPACK_interface.f90"
#include "math.f90"
@ -17,38 +16,38 @@
#include "results.f90"
#include "geometry_plastic_nonlocal.f90"
#include "discretization.f90"
#ifdef Marc4DAMASK
#include "marc/discretization_marc.f90"
#endif
#include "material.f90"
#include "lattice.f90"
#include "constitutive.f90"
#include "constitutive_mech.f90"
#include "constitutive_plastic_none.f90"
#include "constitutive_plastic_isotropic.f90"
#include "constitutive_plastic_phenopowerlaw.f90"
#include "constitutive_plastic_kinehardening.f90"
#include "constitutive_plastic_dislotwin.f90"
#include "constitutive_plastic_disloTungsten.f90"
#include "constitutive_plastic_nonlocal.f90"
#include "constitutive_thermal.f90"
#include "constitutive_thermal_dissipation.f90"
#include "constitutive_thermal_externalheat.f90"
#include "kinematics_thermal_expansion.f90"
#include "constitutive_damage.f90"
#include "source_damage_isoBrittle.f90"
#include "source_damage_isoDuctile.f90"
#include "source_damage_anisoBrittle.f90"
#include "source_damage_anisoDuctile.f90"
#include "kinematics_cleavage_opening.f90"
#include "kinematics_slipplane_opening.f90"
#include "phase.f90"
#include "phase_mechanics.f90"
#include "phase_mechanics_plastic.f90"
#include "phase_mechanics_plastic_none.f90"
#include "phase_mechanics_plastic_isotropic.f90"
#include "phase_mechanics_plastic_phenopowerlaw.f90"
#include "phase_mechanics_plastic_kinehardening.f90"
#include "phase_mechanics_plastic_dislotwin.f90"
#include "phase_mechanics_plastic_dislotungsten.f90"
#include "phase_mechanics_plastic_nonlocal.f90"
#include "phase_mechanics_eigendeformation.f90"
#include "phase_mechanics_eigendeformation_cleavageopening.f90"
#include "phase_mechanics_eigendeformation_slipplaneopening.f90"
#include "phase_mechanics_eigendeformation_thermalexpansion.f90"
#include "phase_thermal.f90"
#include "phase_thermal_dissipation.f90"
#include "phase_thermal_externalheat.f90"
#include "phase_damage.f90"
#include "phase_damage_isobrittle.f90"
#include "phase_damage_isoductile.f90"
#include "phase_damage_anisobrittle.f90"
#include "phase_damage_anisoductile.f90"
#include "damage_none.f90"
#include "damage_nonlocal.f90"
#include "homogenization.f90"
#include "homogenization_mech.f90"
#include "homogenization_mech_none.f90"
#include "homogenization_mech_isostrain.f90"
#include "homogenization_mech_RGC.f90"
#include "homogenization_mechanics.f90"
#include "homogenization_mechanics_none.f90"
#include "homogenization_mechanics_isostrain.f90"
#include "homogenization_mechanics_RGC.f90"
#include "homogenization_thermal.f90"
#include "homogenization_damage.f90"
#include "CPFEM.f90"

137
src/constitutive_thermal_externalheat.f90
View File

@ -1,137 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Michigan State University
!> @brief material subroutine for variable heat source
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_thermal) source_externalheat
integer, dimension(:), allocatable :: &
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: &
t_n, &
f_T
integer :: &
nIntervals
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function source_thermal_externalheat_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: &
phases, &
phase, &
sources, thermal, &
src
integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- thermal_externalheat init -+>>>'
mySources = thermal_active('externalheat',source_length)
Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return
phases => config_material%get('phase')
allocate(param(Ninstances))
allocate(source_thermal_externalheat_offset (phases%length), source=0)
allocate(source_thermal_externalheat_instance(phases%length), source=0)
do p = 1, phases%length
phase => phases%get(p)
if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p))
if(count(mySources(:,p)) == 0) cycle
thermal => phase%get('thermal')
sources => thermal%get('source')
do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then
source_thermal_externalheat_offset(p) = sourceOffset
associate(prm => param(source_thermal_externalheat_instance(p)))
src => sources%get(sourceOffset)
prm%t_n = src%get_asFloats('t_n')
prm%nIntervals = size(prm%t_n) - 1
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0)
end associate
endif
enddo
enddo
end function source_thermal_externalheat_init
!--------------------------------------------------------------------------------------------------
!> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers
!--------------------------------------------------------------------------------------------------
module subroutine source_thermal_externalheat_dotState(ph, me)
integer, intent(in) :: &
ph, &
me
integer :: &
sourceOffset
sourceOffset = source_thermal_externalheat_offset(ph)
thermalState(ph)%p(sourceOffset)%dotState(1,me) = 1.0_pReal ! state is current time
end subroutine source_thermal_externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate
!--------------------------------------------------------------------------------------------------
module subroutine thermal_externalheat_getRate(TDot, ph, me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(out) :: &
TDot
integer :: &
sourceOffset, interval
real(pReal) :: &
frac_time
sourceOffset = source_thermal_externalheat_offset(ph)
associate(prm => param(source_thermal_externalheat_instance(ph)))
do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (thermalState(ph)%p(sourceOffset)%state(1,me) - prm%t_n(interval)) &
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside me bounds
enddo
end associate
end subroutine thermal_externalheat_getRate
end submodule source_externalheat

2
src/damage_nonlocal.f90
View File

@ -8,7 +8,7 @@ module damage_nonlocal
use config
use YAML_types
use lattice
use constitutive
use phase
use results
implicit none

35
src/future.f90
View File

@ -1,35 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief New fortran functions for compiler versions that do not support them
!--------------------------------------------------------------------------------------------------
module future
use prec
implicit none
public
contains
#if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800
!--------------------------------------------------------------------------------------------------
!> @brief substitute for the findloc intrinsic (only for integer, dimension(:) at the moment)
!--------------------------------------------------------------------------------------------------
function findloc(a,v)
integer, intent(in), dimension(:) :: a
integer, intent(in) :: v
integer :: i,j
integer, allocatable, dimension(:) :: findloc
allocate(findloc(count(a==v)))
j = 1
do i = 1, size(a)
if (a(i)==v) then
findloc(j) = i
j = j + 1
endif
enddo
end function findloc
#endif
end module future

4
src/homogenization.f90
View File

@ -10,7 +10,7 @@ module homogenization
use config
use math
use material
use constitutive
use phase
use discretization
use damage_none
use damage_nonlocal
@ -145,8 +145,6 @@ module homogenization
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co
real(pReal) :: M
end function damage_nonlocal_getMobility

4
src/homogenization_damage.f90
View File

@ -34,8 +34,8 @@ module subroutine damage_init()
integer :: ho
print'(/,a)', ' <<<+- homogenization_damage init -+>>>'
print'(/,a)', ' <<<+- homogenization:damage init -+>>>'
print'(/,a)', ' <<<+- homogenization:damage:isodamage init -+>>>'
configHomogenizations => config_material%get('homogenization')
allocate(param(configHomogenizations%length))

6
src/homogenization_mech.f90 → src/homogenization_mechanics.f90
View File

@ -2,7 +2,7 @@
!> @author Martin Diehl, KU Leuven
!> @brief Partition F and homogenize P/dPdF
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech
submodule(homogenization) mechanics
interface
@ -86,7 +86,7 @@ module subroutine mech_init(num_homog)
class(tNode), pointer :: &
num_homogMech
print'(/,a)', ' <<<+- homogenization_mech init -+>>>'
print'(/,a)', ' <<<+- homogenization:mechanics init -+>>>'
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity
@ -253,4 +253,4 @@ module subroutine mech_results(group_base,h)
end subroutine mech_results
end submodule homogenization_mech
end submodule mechanics

6
src/homogenization_mech_RGC.f90 → src/homogenization_mechanics_RGC.f90
View File

@ -6,7 +6,7 @@
!> @brief Relaxed grain cluster (RGC) homogenization scheme
!> N_constituents is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
submodule(homogenization:homogenization_mech) homogenization_mech_RGC
submodule(homogenization:mechanics) RGC
use rotations
use lattice
@ -88,7 +88,7 @@ module subroutine mech_RGC_init(num_homogMech)
homog, &
homogMech
print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>'
print'(/,a)', ' <<<+- homogenization:mechanics:RGC init -+>>>'
Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
@ -947,4 +947,4 @@ pure function interface1to4(iFace1D, nGDim)
end function interface1to4
end submodule homogenization_mech_RGC
end submodule RGC

6
src/homogenization_mech_isostrain.f90 → src/homogenization_mechanics_isostrain.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Isostrain (full constraint Taylor assuption) homogenization scheme
!--------------------------------------------------------------------------------------------------
submodule(homogenization:homogenization_mech) homogenization_mech_isostrain
submodule(homogenization:mechanics) isostrain
enum, bind(c); enumerator :: &
parallel_ID, &
@ -37,7 +37,7 @@ module subroutine mech_isostrain_init
homog, &
homogMech
print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>'
print'(/,a)', ' <<<+- homogenization:mechanics:isostrain init -+>>>'
Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
@ -114,4 +114,4 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
end subroutine mech_isostrain_averageStressAndItsTangent
end submodule homogenization_mech_isostrain
end submodule isostrain

6
src/homogenization_mech_none.f90 → src/homogenization_mechanics_none.f90
View File

@ -4,7 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief dummy homogenization homogenization scheme for 1 constituent per material point
!--------------------------------------------------------------------------------------------------
submodule(homogenization:homogenization_mech) homogenization_mech_none
submodule(homogenization:mechanics) none
contains
@ -18,7 +18,7 @@ module subroutine mech_none_init
h, &
Nmaterialpoints
print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>'
print'(/,a)', ' <<<+- homogenization:mechanics:none init -+>>>'
Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
@ -38,4 +38,4 @@ module subroutine mech_none_init
end subroutine mech_none_init
end submodule homogenization_mech_none
end submodule none

24
src/homogenization_thermal.f90
View File

@ -1,7 +1,7 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_thermal
submodule(homogenization) thermal
use lattice
@ -34,7 +34,9 @@ module subroutine thermal_init()
integer :: ho
print'(/,a)', ' <<<+- homogenization_thermal init -+>>>'
print'(/,a)', ' <<<+- homogenization:thermal init -+>>>'
print'(/,a)', ' <<<+- homogenization:thermal:isotemperature init -+>>>'
configHomogenizations => config_material%get('homogenization')
@ -225,15 +227,21 @@ module subroutine thermal_conduction_getSource(Tdot, ip,el)
real(pReal), intent(out) :: &
Tdot
integer :: &
homog
integer :: co, ho,ph,me
real(pReal) :: dot_T_temp
homog = material_homogenizationAt(el)
call constitutive_thermal_getRate(TDot, ip,el)
ho = material_homogenizationAt(el)
Tdot = 0.0_pReal
do co = 1, homogenization_Nconstituents(ho)
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
call constitutive_thermal_getRate(dot_T_temp, ph,me)
Tdot = Tdot + dot_T_temp
enddo
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
Tdot = Tdot/real(homogenization_Nconstituents(ho),pReal)
end subroutine thermal_conduction_getSource
end submodule homogenization_thermal
end submodule thermal

142
src/constitutive.f90 → src/phase.f90
View File

@ -3,7 +3,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief elasticity, plasticity, damage & thermal internal microstructure state
!--------------------------------------------------------------------------------------------------
module constitutive
module phase
use prec
use math
use rotations
@ -15,18 +15,10 @@ module constitutive
use discretization
use parallelization
use HDF5_utilities
use results
implicit none
private
enum, bind(c); enumerator :: &
KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, &
KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID
end enum
type(rotation), dimension(:,:,:), allocatable :: &
crystallite_orientation !< current orientation
@ -65,10 +57,6 @@ module constitutive
type(tDebugOptions) :: debugCrystallite
integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_kinematics !< active kinematic mechanisms of each phase
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
thermal_Nsources, &
phase_Nsources, & !< number of source mechanisms active in each phase
@ -188,6 +176,11 @@ module constitutive
real(pReal) :: T
end function thermal_T
module function thermal_dot_T(ph,me) result(dot_T)
integer, intent(in) :: ph,me
real(pReal) :: dot_T
end function thermal_dot_T
module subroutine constitutive_mech_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F
@ -246,16 +239,12 @@ module constitutive
dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRate(TDot, ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
module subroutine constitutive_thermal_getRate(TDot, ph,me)
integer, intent(in) :: ph, me
real(pReal), intent(out) :: &
TDot
end subroutine constitutive_thermal_getRate
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
integer, intent(in) :: &
instance, &
@ -265,66 +254,31 @@ module constitutive
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
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, co, ip, el)
integer, intent(in) :: &
co, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
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, co, ip, el)
integer, intent(in) :: &
co, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
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, ph,me)
integer, intent(in) :: ph, me
!< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
end subroutine kinematics_thermal_expansion_LiAndItsTangent
module subroutine constitutive_plastic_dependentState(co,ip,el)
module subroutine plastic_dependentState(co,ip,el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine constitutive_plastic_dependentState
end subroutine plastic_dependentState
end interface
type(tDebugOptions) :: debugConstitutive
#if __INTEL_COMPILER >= 1900
public :: &
prec, &
math, &
rotations, &
IO, &
config, &
material, &
results, &
lattice, &
discretization, &
HDF5_utilities
#endif
public :: &
constitutive_init, &
@ -336,7 +290,6 @@ module constitutive
constitutive_forward, &
constitutive_restore, &
plastic_nonlocal_updateCompatibility, &
kinematics_active, &
converged, &
crystallite_init, &
crystallite_stress, &
@ -353,15 +306,10 @@ module constitutive
constitutive_mech_getP, &
constitutive_mech_setF, &
constitutive_mech_getF, &
constitutive_windForward, &
KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, &
KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID
constitutive_windForward
contains
!--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics
!--------------------------------------------------------------------------------------------------
@ -375,7 +323,7 @@ subroutine constitutive_init
phases
print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT)
print'(/,a)', ' <<<+- phase init -+>>>'; flush(IO_STDOUT)
debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic')
@ -410,38 +358,6 @@ subroutine constitutive_init
end subroutine constitutive_init
!--------------------------------------------------------------------------------------------------
!> @brief checks if a kinematic mechanism is active or not
!--------------------------------------------------------------------------------------------------
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
integer :: p,k
phases => config_material%get('phase')
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
!--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
@ -562,9 +478,7 @@ subroutine crystallite_init()
class(tNode), pointer :: &
num_crystallite, &
debug_crystallite, & ! pointer to debug options for crystallite
phases, &
phase, &
mech
phases
print'(/,a)', ' <<<+- crystallite init -+>>>'
@ -630,7 +544,7 @@ subroutine crystallite_init()
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call crystallite_orientations(co,ip,el)
call constitutive_plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
enddo
enddo
enddo
@ -788,4 +702,4 @@ subroutine constitutive_restartRead(fileHandle)
end subroutine constitutive_restartRead
end module constitutive
end module phase

243
src/constitutive_damage.f90 → src/phase_damage.f90
View File

@ -1,7 +1,7 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all damage sources and kinematics constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_damage
submodule(phase) damagee
enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, &
@ -22,35 +22,26 @@ submodule(constitutive) constitutive_damage
interface
module function source_damage_anisoBrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function source_damage_anisoBrittle_init
module function anisobrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function anisobrittle_init
module function source_damage_anisoDuctile_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function source_damage_anisoDuctile_init
module function anisoductile_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function anisoductile_init
module function source_damage_isoBrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function source_damage_isoBrittle_init
module function isobrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function isobrittle_init
module function source_damage_isoDuctile_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function source_damage_isoDuctile_init
module function isoductile_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function isoductile_init
module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_cleavage_opening_init
module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_slipplane_opening_init
module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me)
integer, intent(in) :: ph,me
@ -61,93 +52,93 @@ submodule(constitutive) constitutive_damage
end subroutine source_damage_isoBrittle_deltaState
module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
module subroutine anisobrittle_dotState(S, 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
end subroutine source_damage_anisoBrittle_dotState
end subroutine anisobrittle_dotState
module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
module subroutine anisoductile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_anisoDuctile_dotState
end subroutine anisoductile_dotState
module subroutine source_damage_isoDuctile_dotState(co, ip, el)
module subroutine isoductile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_isoDuctile_dotState
end subroutine isoductile_dotState
module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_anisoBrittle_getRateAndItsTangent
module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_anisoBrittle_getRateAndItsTangent
module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_anisoDuctile_getRateAndItsTangent
module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_anisoDuctile_getRateAndItsTangent
module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_isoBrittle_getRateAndItsTangent
module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_isoBrittle_getRateAndItsTangent
module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_isoDuctile_getRateAndItsTangent
module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_isoDuctile_getRateAndItsTangent
module subroutine source_damage_anisoBrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_anisoBrittle_results
module subroutine anisobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine anisobrittle_results
module subroutine source_damage_anisoDuctile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_anisoDuctile_results
module subroutine anisoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine anisoductile_results
module subroutine source_damage_isoBrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_isoBrittle_results
module subroutine isobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine isobrittle_results
module subroutine source_damage_isoDuctile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_isoDuctile_results
module subroutine isoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine isoductile_results
end interface
@ -164,19 +155,21 @@ module subroutine damage_init
class(tNode), pointer :: &
phases, &
phase, &
sources, &
kinematics
sources
print'(/,a)', ' <<<+- phase:damage init -+>>>'
phases => config_material%get('phase')
allocate(current(phases%length))
allocate(damageState (phases%length))
allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics
allocate(phase_Nsources(phases%length),source = 0)
do ph = 1,phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == ph)
allocate(current(ph)%phi(Nconstituents),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal)
@ -191,26 +184,10 @@ module subroutine damage_init
! initialize source mechanisms
if(maxval(phase_Nsources) /= 0) then
where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID
where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID
where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID
where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID
endif
!--------------------------------------------------------------------------------------------------
! initialize kinematic mechanisms
allocate(phase_Nkinematics(phases%length),source = 0)
do ph = 1,phases%length
phase => phases%get(ph)
kinematics => phase%get('kinematics',defaultVal=emptyList)
phase_Nkinematics(ph) = kinematics%length
enddo
allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID)
if(maxval(phase_Nkinematics) /= 0) then
where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID
where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
where(isobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID
where(isoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID
where(anisobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID
where(anisoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID
endif
end subroutine damage_init
@ -234,30 +211,30 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
localphiDot, &
dLocalphiDot_dPhi
integer :: &
phase, &
grain, &
source, &
constituent
ph, &
co, &
so, &
me
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
do so = 1, phase_Nsources(ph)
select case(phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case (DAMAGE_ISODUCTILE_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case (DAMAGE_ANISOBRITTLE_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case (DAMAGE_ANISODUCTILE_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case default
localphiDot = 0.0_pReal
@ -395,16 +372,16 @@ module subroutine damage_results(group,ph)
sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_results(ph,group//'sources/')
call isobrittle_results(ph,group//'sources/')
case (DAMAGE_ISODUCTILE_ID) sourceType
call source_damage_isoDuctile_results(ph,group//'sources/')
call isoductile_results(ph,group//'sources/')
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call source_damage_anisoBrittle_results(ph,group//'sources/')
call anisobrittle_results(ph,group//'sources/')
case (DAMAGE_ANISODUCTILE_ID) sourceType
call source_damage_anisoDuctile_results(ph,group//'sources/')
call anisoductile_results(ph,group//'sources/')
end select sourceType
@ -436,13 +413,13 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken)
sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISODUCTILE_ID) sourceType
call source_damage_isoDuctile_dotState(co, ip, el)
call isoductile_dotState(co, ip, el)
case (DAMAGE_ANISODUCTILE_ID) sourceType
call source_damage_anisoDuctile_dotState(co, ip, el)
call anisoductile_dotState(co, ip, el)
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call source_damage_anisoBrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress?
call anisobrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress?
end select sourceType
@ -551,4 +528,4 @@ module function constitutive_damage_get_phi(co,ip,el) result(phi)
end function constitutive_damage_get_phi
end submodule constitutive_damage
end submodule damagee

18
src/source_damage_anisoBrittle.f90 → src/phase_damage_anisobrittle.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating anisotropic brittle damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
submodule (phase:damagee) anisobrittle
integer, dimension(:), allocatable :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
@ -35,7 +35,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function source_damage_anisoBrittle_init(source_length) result(mySources)
module function anisobrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
@ -49,7 +49,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
integer, dimension(:), allocatable :: N_cl
character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_anisoBrittle init -+>>>'
print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>'
mySources = source_active('damage_anisoBrittle',source_length)
Ninstances = count(mySources)
@ -114,13 +114,13 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
enddo
enddo
end function source_damage_anisoBrittle_init
end function anisobrittle_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
module subroutine anisobrittle_dotState(S, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -163,7 +163,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
enddo
end associate
end subroutine source_damage_anisoBrittle_dotState
end subroutine anisobrittle_dotState
!--------------------------------------------------------------------------------------------------
@ -196,7 +196,7 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_anisoBrittle_results(phase,group)
module subroutine anisobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
@ -213,6 +213,6 @@ module subroutine source_damage_anisoBrittle_results(phase,group)
enddo outputsLoop
end associate
end subroutine source_damage_anisoBrittle_results
end subroutine anisobrittle_results
end submodule source_damage_anisoBrittle
end submodule anisobrittle

18
src/source_damage_anisoDuctile.f90 → src/phase_damage_anisoductile.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating anisotropic ductile damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
submodule(phase:damagee) anisoductile
integer, dimension(:), allocatable :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
@ -28,7 +28,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function source_damage_anisoDuctile_init(source_length) result(mySources)
module function anisoductile_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
@ -44,7 +44,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_anisoDuctile init -+>>>'
print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>'
mySources = source_active('damage_anisoDuctile',source_length)
Ninstances = count(mySources)
@ -101,13 +101,13 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
enddo
end function source_damage_anisoDuctile_init
end function anisoductile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
module subroutine anisoductile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -132,7 +132,7 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
= sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate
end subroutine source_damage_anisoDuctile_dotState
end subroutine anisoductile_dotState
!--------------------------------------------------------------------------------------------------
@ -165,7 +165,7 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_anisoDuctile_results(phase,group)
module subroutine anisoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
@ -182,6 +182,6 @@ module subroutine source_damage_anisoDuctile_results(phase,group)
enddo outputsLoop
end associate
end subroutine source_damage_anisoDuctile_results
end subroutine anisoductile_results
end submodule source_damage_anisoDuctile
end submodule anisoductile

14
src/source_damage_isoBrittle.f90 → src/phase_damage_isobrittle.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine incoprorating isotropic brittle damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_damage) source_damage_isoBrittle
submodule(phase:damagee) isobrittle
integer, dimension(:), allocatable :: &
source_damage_isoBrittle_offset, &
@ -26,7 +26,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function source_damage_isoBrittle_init(source_length) result(mySources)
module function isobrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
@ -39,7 +39,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
integer :: Ninstances,sourceOffset,Nconstituents,p
character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_isoBrittle init -+>>>'
print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>'
mySources = source_active('damage_isoBrittle',source_length)
Ninstances = count(mySources)
@ -88,7 +88,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
enddo
end function source_damage_isoBrittle_init
end function isobrittle_init
!--------------------------------------------------------------------------------------------------
@ -161,7 +161,7 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoBrittle_results(phase,group)
module subroutine isobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
@ -178,6 +178,6 @@ module subroutine source_damage_isoBrittle_results(phase,group)
enddo outputsLoop
end associate
end subroutine source_damage_isoBrittle_results
end subroutine isobrittle_results
end submodule source_damage_isoBrittle
end submodule isobrittle

18
src/source_damage_isoDuctile.f90 → src/phase_damage_isoductile.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating isotropic ductile damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule (constitutive:constitutive_damage) source_damage_isoDuctile
submodule(phase:damagee) isoductile
integer, dimension(:), allocatable :: &
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
@ -28,7 +28,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function source_damage_isoDuctile_init(source_length) result(mySources)
module function isoductile_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
@ -41,7 +41,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
integer :: Ninstances,sourceOffset,Nconstituents,p
character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_isoDuctile init -+>>>'
print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>'
mySources = source_active('damage_isoDuctile',source_length)
Ninstances = count(mySources)
@ -92,13 +92,13 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
enddo
end function source_damage_isoDuctile_init
end function isoductile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoDuctile_dotState(co, ip, el)
module subroutine isoductile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -123,7 +123,7 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el)
sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate
end subroutine source_damage_isoDuctile_dotState
end subroutine isoductile_dotState
!--------------------------------------------------------------------------------------------------
@ -156,7 +156,7 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoDuctile_results(phase,group)
module subroutine isoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
@ -173,6 +173,6 @@ module subroutine source_damage_isoDuctile_results(phase,group)
enddo outputsLoop
end associate
end subroutine source_damage_isoDuctile_results
end subroutine isoductile_results
end submodule source_damage_isoDuctile
end submodule isoductile

645
src/constitutive_mech.f90 → src/phase_mechanics.f90
View File

@ -1,7 +1,7 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all plasticity constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_mech
submodule(phase) mechanics
enum, bind(c); enumerator :: &
ELASTICITY_UNDEFINED_ID, &
@ -15,9 +15,15 @@ submodule(constitutive) constitutive_mech
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOTUNGSTEN_ID, &
PLASTICITY_NONLOCAL_ID
PLASTICITY_NONLOCAL_ID, &
KINEMATICS_UNDEFINED_ID, &
KINEMATICS_CLEAVAGE_OPENING_ID, &
KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID
end enum
integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_kinematics
integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: &
@ -48,225 +54,85 @@ submodule(constitutive) constitutive_mech
interface
module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_none_init
module subroutine eigendeformation_init(phases)
class(tNode), pointer :: phases
end subroutine eigendeformation_init
module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_isotropic_init
module subroutine plastic_init
end subroutine plastic_init
module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_phenopowerlaw_init
module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_kinehardening_init
module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_dislotwin_init
module function plastic_dislotungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_dislotungsten_init
module function plastic_nonlocal_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_nonlocal_init
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
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) :: &
Mp !< Mandel stress
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LpAndItsTangent
me
end subroutine plastic_isotropic_LiAndItsTangent
pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_phenopowerlaw_LpAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken)
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_kinehardening_LpAndItsTangent
module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
instance, &
of
end subroutine plastic_dislotwin_LpAndItsTangent
pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
instance, &
of
end subroutine plastic_dislotungsten_LpAndItsTangent
module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,instance,of,ip,el)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
Temperature
integer, intent(in) :: &
instance, &
of, &
ip, & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_LpAndItsTangent
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
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
end subroutine plastic_disloTungsten_dotState
module subroutine plastic_nonlocal_dotState(Mp,Temperature,timestep,instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
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
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
real(pReal), intent(in) :: &
subdt !< timestep
logical :: broken
end function plastic_dotState
module subroutine plastic_dislotwin_dependentState(T,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal), intent(in) :: &
T
end subroutine plastic_dislotwin_dependentState
module subroutine plastic_dislotungsten_dependentState(instance,of)
integer, intent(in) :: &
instance, &
of
end subroutine plastic_dislotungsten_dependentState
module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
module function plastic_deltaState(co, ip, el, ph, me) result(broken)
integer, intent(in) :: &
instance, &
of, &
ip, & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_dependentState
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
logical :: &
broken
end function plastic_deltaState
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
module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el)
module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
instance, &
of, &
ip, &
el
end subroutine plastic_nonlocal_deltaState
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
end subroutine constitutive_LiAndItsTangents
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_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
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 plastic_LpAndItsTangents
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
@ -340,7 +206,7 @@ module subroutine mech_init(phases)
elastic, &
stiffDegradation
print'(/,a)', ' <<<+- constitutive_mech init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics init -+>>>'
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
@ -444,13 +310,7 @@ module subroutine mech_init(phases)
allocate(phase_plasticityInstance(phases%length),source = 0)
allocate(phase_localPlasticity(phases%length), source=.true.)
where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID
where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID
where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID
where(plastic_kinehardening_init()) phase_plasticity = PLASTICITY_KINEHARDENING_ID
where(plastic_dislotwin_init()) phase_plasticity = PLASTICITY_DISLOTWIN_ID
where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID
where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID
call plastic_init()
do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
@ -481,36 +341,13 @@ module subroutine mech_init(phases)
end select
call eigendeformation_init(phases)
end subroutine mech_init
!--------------------------------------------------------------------------------------------------
!> @brief checks if a plastic module is active or not
!--------------------------------------------------------------------------------------------------
function plastic_active(plastic_label) result(active_plastic)
character(len=*), intent(in) :: plastic_label !< type of plasticity model
logical, dimension(:), allocatable :: active_plastic
class(tNode), pointer :: &
phases, &
phase, &
mech, &
pl
integer :: ph
phases => config_material%get('phase')
allocate(active_plastic(phases%length), source = .false. )
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanics')
pl => mech%get('plasticity')
if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true.
enddo
end function plastic_active
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
@ -559,216 +396,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
end subroutine constitutive_hooke_SandItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different plasticity constitutive models
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_plastic_dependentState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
ph, &
instance, me
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dependentState(thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_dislotungsten_dependentState(instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dependentState(instance,me,ip,el)
end select plasticityType
end subroutine constitutive_plastic_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_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
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
real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp
integer :: &
i, j, instance, me, ph
Mp = matmul(matmul(transpose(Fi),Fi),S)
me = material_phasememberAt(co,ip,el)
ph = material_phaseAt(co,el)
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me)
end select plasticityType
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo; enddo
end subroutine constitutive_plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function mech_collectDotState(subdt,co,ip,el,ph,me) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
real(pReal), intent(in) :: &
subdt !< timestep
real(pReal), dimension(3,3) :: &
Mp
integer :: &
instance
logical :: broken
instance = phase_plasticityInstance(ph)
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me))
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState(Mp,instance,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Mp,instance,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Mp,instance,me)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState(Mp,thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_disloTungsten_dotState(Mp,thermal_T(ph,me),instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState(Mp,thermal_T(ph,me),subdt,instance,me,ip,el)
end select plasticityType
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me)))
end function mech_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_deltaState(co, ip, el, ph, of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
logical :: &
broken
real(pReal), dimension(3,3) :: &
Mp
integer :: &
instance, &
myOffset, &
mySize
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,of)),constitutive_mech_S(ph)%data(1:3,1:3,of))
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Mp,instance,of)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of)))
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of)))
case default
broken = .false.
end select plasticityType
if(.not. broken) then
select case(phase_plasticity(ph))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) + plasticState(ph)%deltaState(1:mySize,of)
end select
endif
end function constitutive_deltaState
module subroutine mech_results(group,ph)
character(len=*), intent(in) :: group
@ -862,7 +489,6 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
ierr, & ! error indicator for LAPACK
o, &
p, &
m, &
ph, &
me, &
jacoCounterLp, &
@ -875,7 +501,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call constitutive_plastic_dependentState(co,ip,el)
call plastic_dependentState(co,ip,el)
Lpguess = constitutive_mech_Lp(ph)%data(1:3,1:3,me) ! take as first guess
Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me) ! take as first guess
@ -915,7 +541,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, co, ip, el)
call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
S, Fi_new, co, ip, el)
!* update current residuum and check for convergence of loop
@ -1051,7 +677,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(ph)%sizeDotState
@ -1067,7 +693,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) exit iteration
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me)
if(broken) exit iteration
zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),&
@ -1080,7 +706,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) &
- r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then
broken = constitutive_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(co,ip,el,ph,me)
exit iteration
endif
@ -1136,14 +762,14 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = subState0 &
+ plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t
broken = constitutive_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(co,ip,el,ph,me)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -1176,7 +802,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(ph)%sizeDotState
@ -1185,13 +811,13 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
plasticState(ph)%state(1:sizeDotState,me) = subState0 &
+ plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t
broken = constitutive_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(co,ip,el,ph,me)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) return
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me)
if(broken) return
broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
@ -1294,7 +920,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(Delta_t,co,ip,el,ph,me)
broken = plastic_dotState(Delta_t,co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(ph)%sizeDotState
@ -1315,7 +941,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el)
if(broken) exit
broken = mech_collectDotState(Delta_t*C(stage),co,ip,el,ph,me)
broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,me)
if(broken) exit
enddo
@ -1334,7 +960,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
if(broken) return
broken = constitutive_deltaState(co,ip,el,ph,me)
broken = plastic_deltaState(co,ip,el,ph,me)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -1351,7 +977,6 @@ subroutine crystallite_results(group,ph)
integer, intent(in) :: ph
integer :: ou
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
real(pReal), allocatable, dimension(:,:) :: selected_rotations
character(len=:), allocatable :: structureLabel
@ -1708,7 +1333,7 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
constitutive_mech_S(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
@ -1880,86 +1505,4 @@ 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, me, ph
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
me = material_phaseMemberAt(co,ip,el)
ph = material_phaseAt(co,el)
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me)
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
end submodule mechanics

207
src/phase_mechanics_eigendeformation.f90 Normal file
View File

@ -0,0 +1,207 @@
submodule(phase:mechanics) eigendeformation
interface
module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_cleavage_opening_init
module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_slipplane_opening_init
module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_thermal_expansion_init
module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
integer, intent(in) :: &
co, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
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, co, ip, el)
integer, intent(in) :: &
co, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_slipplane_opening_LiAndItsTangent
module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
integer, intent(in) :: ph, me
!< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
end subroutine thermalexpansion_LiAndItsTangent
end interface
contains
module subroutine eigendeformation_init(phases)
class(tNode), pointer :: &
phases
integer :: &
ph
class(tNode), pointer :: &
phase, &
kinematics
print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>'
!--------------------------------------------------------------------------------------------------
! initialize kinematic mechanisms
allocate(phase_Nkinematics(phases%length),source = 0)
do ph = 1,phases%length
phase => phases%get(ph)
kinematics => phase%get('kinematics',defaultVal=emptyList)
phase_Nkinematics(ph) = kinematics%length
enddo
allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID)
if(maxval(phase_Nkinematics) /= 0) then
where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID
where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID
endif
end subroutine eigendeformation_init
!--------------------------------------------------------------------------------------------------
!> @brief checks if a kinematic mechanism is active or not
!--------------------------------------------------------------------------------------------------
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
integer :: p,k
phases => config_material%get('phase')
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
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
module 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, me, ph
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
me = material_phaseMemberAt(co,ip,el)
ph = material_phaseAt(co,el)
call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me)
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 eigendeformation

6
src/kinematics_cleavage_opening.f90 → src/phase_mechanics_eigendeformation_cleavageopening.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_damage) kinematics_cleavage_opening
submodule(phase:eigendeformation) cleavageopening
integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance
@ -44,7 +44,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
kinematics, &
kinematic_type
print'(/,a)', ' <<<+- kinematics_cleavage_opening init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>'
myKinematics = kinematics_active('cleavage_opening',kinematics_length)
Ninstances = count(myKinematics)
@ -162,4 +162,4 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
end subroutine kinematics_cleavage_opening_LiAndItsTangent
end submodule kinematics_cleavage_opening
end submodule cleavageopening

6
src/kinematics_slipplane_opening.f90 → src/phase_mechanics_eigendeformation_slipplaneopening.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating kinematics resulting from opening of slip planes
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_damage) kinematics_slipplane_opening
submodule(phase:eigendeformation) slipplaneopening
integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance
@ -49,7 +49,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
kinematics, &
kinematic_type
print'(/,a)', ' <<<+- kinematics_slipplane init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>'
myKinematics = kinematics_active('slipplane_opening',kinematics_length)
Ninstances = count(myKinematics)
@ -193,4 +193,4 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
end subroutine kinematics_slipplane_opening_LiAndItsTangent
end submodule kinematics_slipplane_opening
end submodule slipplaneopening

40
src/kinematics_thermal_expansion.f90 → src/phase_mechanics_eigendeformation_thermalexpansion.f90
View File

@ -3,7 +3,7 @@
!> @brief material subroutine incorporating kinematics resulting from thermal expansion
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion
submodule(phase:eigendeformation) thermalexpansion
integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
@ -16,7 +16,6 @@ submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion
type(tParameters), dimension(:), allocatable :: param
contains
@ -37,7 +36,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
kinematics, &
kinematic_type
print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:thermalexpansion init -+>>>'
myKinematics = kinematics_active('thermal_expansion',kinematics_length)
Ninstances = count(myKinematics)
@ -84,7 +83,7 @@ end function kinematics_thermal_expansion_init
!--------------------------------------------------------------------------------------------------
!> @brief constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
integer, intent(in) :: ph, me
real(pReal), intent(out), dimension(3,3) :: &
@ -92,28 +91,25 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, p
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer :: &
phase, &
homog
real(pReal) :: T, dot_T
T = current(ph)%T(me)
dot_T = current(ph)%dot_T(me)
T = thermal_T(ph,me)
dot_T = thermal_dot_T(ph,me)
associate(prm => param(kinematics_thermal_expansion_instance(ph)))
Li = dot_T * ( &
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
)
end associate
Li = dot_T * ( &
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
)
end associate
dLi_dTstar = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent
end subroutine thermalexpansion_LiAndItsTangent
end submodule kinematics_thermal_expansion
end submodule thermalexpansion

472
src/phase_mechanics_plastic.f90 Normal file
View File

@ -0,0 +1,472 @@
submodule(phase:mechanics) plastic
interface
module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_none_init
module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_isotropic_init
module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_phenopowerlaw_init
module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_kinehardening_init
module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_dislotwin_init
module function plastic_dislotungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_dislotungsten_init
module function plastic_nonlocal_init() result(myPlasticity)
logical, dimension(:), allocatable :: &
myPlasticity
end function plastic_nonlocal_init
module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
ph, &
me
end subroutine isotropic_LpAndItsTangent
pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
ph, &
me
end subroutine phenopowerlaw_LpAndItsTangent
pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
ph, &
me
end subroutine kinehardening_LpAndItsTangent
module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
real(pReal), dimension(3,3), intent(in) :: &
Mp
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
me
end subroutine dislotwin_LpAndItsTangent
pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
real(pReal), dimension(3,3), intent(in) :: &
Mp
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
me
end subroutine dislotungsten_LpAndItsTangent
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,ph,me,ip,el)
real(pReal), dimension(3,3), intent(out) :: &
Lp
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
Temperature
integer, intent(in) :: &
ph, &
me, &
ip, & !< current integration point
el !< current element number
end subroutine nonlocal_LpAndItsTangent
module subroutine isotropic_dotState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
me
end subroutine isotropic_dotState
module subroutine phenopowerlaw_dotState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
me
end subroutine phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
me
end subroutine plastic_kinehardening_dotState
module subroutine dislotwin_dotState(Mp,T,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
me
end subroutine dislotwin_dotState
module subroutine dislotungsten_dotState(Mp,T,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
me
end subroutine dislotungsten_dotState
module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,me,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
me, &
ip, & !< current integration point
el !< current element number
end subroutine nonlocal_dotState
module subroutine dislotwin_dependentState(T,instance,me)
integer, intent(in) :: &
instance, &
me
real(pReal), intent(in) :: &
T
end subroutine dislotwin_dependentState
module subroutine dislotungsten_dependentState(instance,me)
integer, intent(in) :: &
instance, &
me
end subroutine dislotungsten_dependentState
module subroutine nonlocal_dependentState(instance, me, ip, el)
integer, intent(in) :: &
instance, &
me, &
ip, & !< current integration point
el !< current element number
end subroutine nonlocal_dependentState
module subroutine plastic_kinehardening_deltaState(Mp,instance,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
me
end subroutine plastic_kinehardening_deltaState
module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer, intent(in) :: &
instance, &
me, &
ip, &
el
end subroutine plastic_nonlocal_deltaState
end interface
contains
module subroutine plastic_init
print'(/,a)', ' <<<+- phase:mechanics:plastic init -+>>>'
where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID
where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID
where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID
where(plastic_kinehardening_init()) phase_plasticity = PLASTICITY_KINEHARDENING_ID
where(plastic_dislotwin_init()) phase_plasticity = PLASTICITY_DISLOTWIN_ID
where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID
where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID
end subroutine plastic_init
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_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
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 me Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp
integer :: &
i, j, me, ph
Mp = matmul(matmul(transpose(Fi),Fi),S)
me = material_phasememberAt(co,ip,el)
ph = material_phaseAt(co,el)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
end select plasticityType
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo; enddo
end subroutine plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
real(pReal), intent(in) :: &
subdt !< timestep
real(pReal), dimension(3,3) :: &
Mp
logical :: broken
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me))
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call isotropic_dotState(Mp,ph,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call phenopowerlaw_dotState(Mp,ph,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Mp,ph,me)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el)
end select plasticityType
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me)))
end function plastic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different plasticity constitutive models
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dependentState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
ph, &
instance, me
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call dislotwin_dependentState(thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call dislotungsten_dependentState(instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call nonlocal_dependentState(instance,me,ip,el)
end select plasticityType
end subroutine plastic_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
module function plastic_deltaState(co, ip, el, ph, me) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
me
logical :: &
broken
real(pReal), dimension(3,3) :: &
Mp
integer :: &
instance, &
myOffset, &
mySize
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me))
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Mp,instance,me)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case default
broken = .false.
end select plasticityType
if(.not. broken) then
select case(phase_plasticity(ph))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) + plasticState(ph)%deltaState(1:mySize,me)
end select
endif
end function plastic_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief checks if a plastic module is active or not
!--------------------------------------------------------------------------------------------------
function plastic_active(plastic_label) result(active_plastic)
character(len=*), intent(in) :: plastic_label !< type of plasticity model
logical, dimension(:), allocatable :: active_plastic
class(tNode), pointer :: &
phases, &
phase, &
mech, &
pl
integer :: ph
phases => config_material%get('phase')
allocate(active_plastic(phases%length), source = .false. )
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanics')
pl => mech%get('plasticity')
if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true.
enddo
end function plastic_active
end submodule plastic

91
src/constitutive_plastic_disloTungsten.f90 → src/phase_mechanics_plastic_dislotungsten.f90
View File

@ -5,7 +5,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief crystal plasticity model for bcc metals, especially Tungsten
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_dislotungsten
submodule(phase:plastic) dislotungsten
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
@ -17,7 +17,7 @@ submodule(constitutive:constitutive_mech) plastic_dislotungsten
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: &
b_sl, & !< magnitude of Burgers vector [m]
b_sl, & !< magnitude me Burgers vector [m]
D_a, &
i_sl, & !< Adj. parameter for distance between 2 forest dislocations
f_at, & !< factor to calculate atomic volume
@ -97,7 +97,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotungsten init -+>>>'
myPlasticity = plastic_active('dislotungsten')
Ninstances = count(myPlasticity)
@ -222,7 +222,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt == p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState
@ -272,8 +272,8 @@ end function plastic_dislotungsten_init
!--------------------------------------------------------------------------------------------------
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
Mp,T,instance,of)
pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
Mp,T,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
@ -284,21 +284,21 @@ pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
instance, &
of
ph, &
me
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
dot_gamma_pos,dot_gamma_neg, &
ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
associate(prm => param(phase_plasticityInstance(ph)))
call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
call kinetics(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl
Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -309,25 +309,25 @@ pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
end associate
end subroutine plastic_dislotungsten_LpAndItsTangent
end subroutine dislotungsten_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotungsten_dotState(Mp,T,instance,of)
module subroutine dislotungsten_dotState(Mp,T,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
instance, &
of
ph, &
me
real(pReal) :: &
VacancyDiffusion
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
gdot_pos, gdot_neg,&
tau_pos,&
tau_neg, &
@ -336,13 +336,14 @@ module subroutine plastic_dislotungsten_dotState(Mp,T,instance,of)
dot_rho_dip_climb, &
dip_distance
associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),&
dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph)))
call kinetics(Mp,T,instance,of,&
call kinetics(Mp,T,phase_plasticityInstance(ph),me,&
gdot_pos,gdot_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs
dot%gamma_sl(:,me) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs
VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg
@ -351,50 +352,50 @@ module subroutine plastic_dislotungsten_dotState(Mp,T,instance,of)
else where
dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), &
prm%D_a, & ! lower limit
dst%Lambda_sl(:,of)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
dst%Lambda_sl(:,me)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) &
* (1.0_pReal/(dip_distance+prm%D_a))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency?
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,me))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency?
end where
dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication
dot%rho_mob(:,me) = abs(dot%gamma_sl(:,me))/(prm%b_sl*dst%Lambda_sl(:,me)) & ! multiplication
- dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations
dot%rho_dip(:,of) = dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me)) ! Spontaneous annihilation of 2 single edge dislocations
dot%rho_dip(:,me) = dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,me)*abs(dot%gamma_sl(:,me)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
- dot_rho_dip_climb
end associate
end subroutine plastic_dislotungsten_dotState
end subroutine dislotungsten_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Calculate derived quantities from state.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotungsten_dependentState(instance,of)
module subroutine dislotungsten_dependentState(instance,me)
integer, intent(in) :: &
instance, &
of
me
real(pReal), dimension(param(instance)%sum_N_sl) :: &
dislocationSpacing
associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dst%threshold_stress(:,of) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
dst%threshold_stress(:,me) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
dst%Lambda_sl(:,me) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
end associate
end subroutine plastic_dislotungsten_dependentState
end subroutine dislotungsten_dependentState
!--------------------------------------------------------------------------------------------------
@ -439,7 +440,7 @@ end subroutine plastic_dislotungsten_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,T,instance,of, &
pure subroutine kinetics(Mp,T,instance,me, &
dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out)
real(pReal), dimension(3,3), intent(in) :: &
@ -448,7 +449,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
T !< temperature
integer, intent(in) :: &
instance, &
of
me
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
dot_gamma_pos, &
@ -479,11 +480,11 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(tau_neg_out)) tau_neg_out = tau_neg
associate(BoltzmannRatio => prm%Q_s/(kB*T), &
dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v_0, &
effectiveLength => dst%Lambda_sl(:,of) - prm%w)
dot_gamma_0 => stt%rho_mob(:,me)*prm%b_sl*prm%v_0, &
effectiveLength => dst%Lambda_sl(:,me) - prm%w)
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_Peierls
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,me))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@ -499,7 +500,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
end where significantPositiveTau
if (present(ddot_gamma_dtau_pos)) then
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_pos
@ -512,8 +513,8 @@ pure subroutine kinetics(Mp,T,instance,of, &
end where significantPositiveTau2
endif
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_Peierls
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,me))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@ -529,7 +530,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
end where significantNegativeTau
if (present(ddot_gamma_dtau_neg)) then
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_neg
@ -547,4 +548,4 @@ pure subroutine kinetics(Mp,T,instance,of, &
end subroutine kinetics
end submodule plastic_dislotungsten
end submodule dislotungsten

158
src/constitutive_plastic_dislotwin.f90 → src/phase_mechanics_plastic_dislotwin.f90
View File

@ -7,7 +7,7 @@
!> @brief material subroutine incoprorating dislocation and twinning physics
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_dislotwin
submodule(phase:plastic) dislotwin
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
@ -144,7 +144,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotwin init -+>>>'
myPlasticity = plastic_active('dislotwin')
Ninstances = count(myPlasticity)
@ -408,7 +408,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt == p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr
@ -521,12 +521,12 @@ end function plastic_dislotwin_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp
real(pReal), dimension(3,3), intent(in) :: Mp
integer, intent(in) :: instance,of
integer, intent(in) :: ph,me
real(pReal), intent(in) :: T
integer :: i,k,l,m,n
@ -535,11 +535,11 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
BoltzmannRatio, &
ddot_gamma_dtau, &
tau
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip
real(pReal), dimension(param(instance)%sum_N_tw) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: &
dot_gamma_twin,ddot_gamma_dtau_twin
real(pReal), dimension(param(instance)%sum_N_tr) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: &
dot_gamma_tr,ddot_gamma_dtau_trans
real(pReal):: dot_gamma_sb
real(pReal), dimension(3,3) :: eigVectors, P_sb
@ -564,16 +564,16 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
0, 1, 1 &
],pReal),[ 3,6])
associate(prm => param(instance), stt => state(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) &
- sum(stt%f_tr(1:prm%sum_N_tr,of))
- sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me))
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
call kinetics_slip(Mp,T,instance,of,dot_gamma_sl,ddot_gamma_dtau_slip)
call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip)
slipContribution: do i = 1, prm%sum_N_sl
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -581,7 +581,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin)
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin)
twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -589,7 +589,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans)
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans)
transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -628,21 +628,21 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
end associate
end subroutine plastic_dislotwin_LpAndItsTangent
end subroutine dislotwin_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
module subroutine dislotwin_dotState(Mp,T,ph,me)
real(pReal), dimension(3,3), intent(in):: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T !< temperature at integration point
integer, intent(in) :: &
instance, &
of
ph, &
me
integer :: i
real(pReal) :: &
@ -653,25 +653,25 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
tau, &
sigma_cl, & !< climb stress
b_d !< ratio of Burgers vector to stacking fault width
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
dot_rho_dip_formation, &
dot_rho_dip_climb, &
rho_dip_distance_min, &
dot_gamma_sl
real(pReal), dimension(param(instance)%sum_N_tw) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: &
dot_gamma_twin
real(pReal), dimension(param(instance)%sum_N_tr) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: &
dot_gamma_tr
associate(prm => param(instance), stt => state(instance), &
dot => dotState(instance), dst => dependentState(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph)))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) &
- sum(stt%f_tr(1:prm%sum_N_tr,of))
- sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me))
call kinetics_slip(Mp,T,instance,of,dot_gamma_sl)
dot%gamma_sl(:,of) = abs(dot_gamma_sl)
call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl)
dot%gamma_sl(:,me) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%D_a*prm%b_sl
@ -683,12 +683,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress
rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,of))
rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me))
rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i))
if (prm%dipoleFormation) then
dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
* stt%rho_mob(i,of)*abs(dot_gamma_sl(i))
* stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
else
dot_rho_dip_formation(i) = 0.0_pReal
endif
@ -707,39 +707,39 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) &
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,me) &
/ (rho_dip_distance-rho_dip_distance_min(i))
endif
endif significantSlipStress
enddo slipState
dot%rho_mob(:,of) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,of)) &
dot%rho_mob(:,me) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,me)) &
- dot_rho_dip_formation &
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,of)*abs(dot_gamma_sl)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,me)*abs(dot_gamma_sl)
dot%rho_dip(:,of) = dot_rho_dip_formation &
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,of)*abs(dot_gamma_sl) &
dot%rho_dip(:,me) = dot_rho_dip_formation &
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
- dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin)
dot%f_tw(:,of) = f_unrotated*dot_gamma_twin/prm%gamma_char
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin)
dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr)
dot%f_tr(:,of) = f_unrotated*dot_gamma_tr
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr)
dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
end associate
end subroutine plastic_dislotwin_dotState
end subroutine dislotwin_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Calculate derived quantities from state.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dislotwin_dependentState(T,instance,of)
module subroutine dislotwin_dependentState(T,instance,me)
integer, intent(in) :: &
instance, &
of
me
real(pReal), intent(in) :: &
T
@ -763,18 +763,18 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
stt => state(instance),&
dst => dependentState(instance))
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of))
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me))
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
!* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ...
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_trans/prm%t_tr ! but this not
! ToDo ...Physically correct, but naming could be adjusted
inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, &
stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%i_sl
stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
@ -787,41 +787,41 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here
dst%Lambda_sl(:,of) = prm%D &
dst%Lambda_sl(:,me) = prm%D &
/ (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr))
else
dst%Lambda_sl(:,of) = prm%D &
dst%Lambda_sl(:,me) = prm%D &
/ (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct?
endif
dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
dst%Lambda_tr(:,of) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion
dst%tau_pass(:,of) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dst%tau_pass(:,me) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
!* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) &
dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) &
dst%tau_hat_tw(:,me) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
if(prm%sum_N_tr == prm%sum_N_sl) &
dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) &
dst%tau_hat_tr(:,me) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal
dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal
dst%V_tw(:,me) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,me)**2.0_pReal
dst%V_tr(:,me) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,me)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
dst%tau_r_tw(:,me) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
dst%tau_r_tr(:,me) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
end associate
end subroutine plastic_dislotwin_dependentState
end subroutine dislotwin_dependentState
!--------------------------------------------------------------------------------------------------
@ -882,7 +882,7 @@ end subroutine plastic_dislotwin_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,T,instance,of, &
pure subroutine kinetics_slip(Mp,T,instance,me, &
dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip)
real(pReal), dimension(3,3), intent(in) :: &
@ -891,7 +891,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
T !< temperature
integer, intent(in) :: &
instance, &
of
me
real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: &
dot_gamma_sl
@ -920,7 +920,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
enddo
tau_eff = abs(tau)-dst%tau_pass(:,of)
tau_eff = abs(tau)-dst%tau_pass(:,me)
significantStress: where(tau_eff > tol_math_check)
stressRatio = tau_eff/prm%tau_0
@ -929,7 +929,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q)
v_run_inverse = prm%B/(tau_eff*prm%b_sl)
dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
dot_gamma_sl = sign(stt%rho_mob(:,me)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio &
* (stressRatio**(prm%p-1.0_pReal)) &
@ -938,7 +938,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff
dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
/ (v_wait_inverse+v_run_inverse)**2.0_pReal
ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,of)*prm%b_sl
ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,me)*prm%b_sl
else where significantStress
dot_gamma_sl = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal
@ -959,7 +959,7 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,&
dot_gamma_twin,ddot_gamma_dtau_twin)
real(pReal), dimension(3,3), intent(in) :: &
@ -968,7 +968,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
T !< temperature
integer, intent(in) :: &
instance, &
of
me
real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: &
dot_gamma_sl
@ -992,11 +992,11 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tw(i,of)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state
if (tau(i) < dst%tau_r_tw(i,me)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/& ! ToDo: MD: it would be more consistent to use shearrates from state
(prm%L_tw*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,of)-tau(i)))) ! P_ncs
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,me)-tau(i)))) ! P_ncs
else
Ndot0=0.0_pReal
end if
@ -1006,8 +1006,8 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
enddo
significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r
dot_gamma_twin = prm%gamma_char * dst%V_tw(:,of) * Ndot0*exp(-StressRatio_r)
StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r
dot_gamma_twin = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r)
ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r
else where significantStress
dot_gamma_twin = 0.0_pReal
@ -1028,7 +1028,7 @@ end subroutine kinetics_twin
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,me,&
dot_gamma_tr,ddot_gamma_dtau_trans)
real(pReal), dimension(3,3), intent(in) :: &
@ -1037,7 +1037,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
T !< temperature
integer, intent(in) :: &
instance, &
of
me
real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: &
dot_gamma_sl
@ -1060,11 +1060,11 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tr(i,of)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state
if (tau(i) < dst%tau_r_tr(i,me)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/& ! ToDo: MD: it would be more consistent to use shearrates from state
(prm%L_tr*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,of)-tau(i)))) ! P_ncs
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,me)-tau(i)))) ! P_ncs
else
Ndot0=0.0_pReal
end if
@ -1074,8 +1074,8 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
enddo
significantStress: where(tau > tol_math_check)
StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s
dot_gamma_tr = dst%V_tr(:,of) * Ndot0*exp(-StressRatio_s)
StressRatio_s = (dst%tau_hat_tr(:,me)/tau)**prm%s
dot_gamma_tr = dst%V_tr(:,me) * Ndot0*exp(-StressRatio_s)
ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
else where significantStress
dot_gamma_tr = 0.0_pReal
@ -1088,4 +1088,4 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
end subroutine kinetics_trans
end submodule plastic_dislotwin
end submodule dislotwin

68
src/constitutive_plastic_isotropic.f90 → src/phase_mechanics_plastic_isotropic.f90
View File

@ -7,7 +7,7 @@
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_isotropic
submodule(phase:plastic) isotropic
type :: tParameters
real(pReal) :: &
@ -68,7 +68,7 @@ module function plastic_isotropic_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_isotropic init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:isotropic init -+>>>'
myPlasticity = plastic_active('isotropic')
Ninstances = count(myPlasticity)
@ -131,7 +131,7 @@ module function plastic_isotropic_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt == p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState
@ -168,7 +168,7 @@ end function plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
@ -178,8 +178,8 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
ph, &
me
real(pReal), dimension(3,3) :: &
Mp_dev !< deviatoric part of the Mandel stress
@ -190,24 +190,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)))
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(me))) **prm%n
Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev
#ifdef DEBUG
if (debugConstitutive%extensive .and. (of == prm%of_debug .or. .not. debugConstitutive%selective)) then
print'(/,a,/,3(12x,3(f12.4,1x)/))', '<< CONST isotropic >> Tstar (dev) / MPa', &
transpose(Mp_dev)*1.0e-6_pReal
print'(/,a,/,f12.5)', '<< CONST isotropic >> norm Tstar / MPa', norm_Mp_dev*1.0e-6_pReal
print'(/,a,/,f12.5)', '<< CONST isotropic >> gdot', dot_gamma
end if
#endif
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
forall (k=1:3,l=1:3) &
@ -222,13 +214,13 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
end associate
end subroutine plastic_isotropic_LpAndItsTangent
end subroutine isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate inelastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
@ -239,7 +231,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
me
real(pReal) :: &
tr !< trace of spherical part of Mandel stress (= 3 x pressure)
@ -252,19 +244,10 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
#ifdef DEBUG
if (debugConstitutive%extensive .and. (of == prm%of_debug .or. .not. debugConstitutive%selective)) then
print'(/,a,/,f12.5)', '<< CONST isotropic >> pressure / MPa', tr/3.0_pReal*1.0e-6_pReal
print'(/,a,/,f12.5)', '<< CONST isotropic >> gdot', prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
end if
#endif
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
else
Li = 0.0_pReal
@ -279,20 +262,21 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_dotState(Mp,instance,of)
module subroutine isotropic_dotState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
ph, &
me
real(pReal) :: &
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)))
if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp))
@ -300,7 +284,7 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(me))) **prm%n
if (dot_gamma > 1e-12_pReal) then
if (dEq0(prm%c_1)) then
@ -310,19 +294,19 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
endif
dot%xi(of) = dot_gamma &
dot%xi(me) = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a &
* sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star)
* abs( 1.0_pReal - stt%xi(me)/xi_inf_star )**prm%a &
* sign(1.0_pReal, 1.0_pReal - stt%xi(me)/xi_inf_star)
else
dot%xi(of) = 0.0_pReal
dot%xi(me) = 0.0_pReal
endif
dot%gamma(of) = dot_gamma ! ToDo: not really used
dot%gamma(me) = dot_gamma ! ToDo: not really used
end associate
end subroutine plastic_isotropic_dotState
end subroutine isotropic_dotState
!--------------------------------------------------------------------------------------------------
@ -348,4 +332,4 @@ module subroutine plastic_isotropic_results(instance,group)
end subroutine plastic_isotropic_results
end submodule plastic_isotropic
end submodule isotropic

87
src/constitutive_plastic_kinehardening.f90 → src/phase_mechanics_plastic_kinehardening.f90
View File

@ -5,7 +5,7 @@
!> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates
!! and a Voce-type kinematic hardening rule
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_kinehardening
submodule(phase:plastic) kinehardening
type :: tParameters
real(pReal) :: &
@ -80,7 +80,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:kinehardening init -+>>>'
myPlasticity = plastic_active('kinehardening')
Ninstances = count(myPlasticity)
@ -175,7 +175,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt == p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
sizeState = sizeDotState + sizeDeltaState
@ -240,7 +240,7 @@ end function plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
@ -250,21 +250,21 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
ph, &
me
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
associate(prm => param(phase_plasticityInstance(ph)))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
@ -276,44 +276,45 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
end associate
end subroutine plastic_kinehardening_LpAndItsTangent
end subroutine kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,instance,of)
module subroutine plastic_kinehardening_dotState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
ph, &
me
real(pReal) :: &
sumGamma
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
gdot_pos,gdot_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),&
dot => dotState(phase_plasticityInstance(ph)))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg)
dot%accshear(:,me) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,me))
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
dot%crss(:,me) = matmul(prm%interaction_SlipSlip,dot%accshear(:,me)) &
* ( prm%h_inf_f &
+ (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
)
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
dot%crss_back(:,me) = stt%sense(:,me)*dot%accshear(:,me) * &
( prm%h_inf_b + &
(prm%h_0_b - prm%h_inf_b &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,of))) &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,me))*(stt%accshear(:,me)-stt%gamma0(:,me))&
) *exp(-(stt%accshear(:,me)-stt%gamma0(:,me)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,me))) &
)
end associate
@ -324,13 +325,13 @@ end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief Calculate (instantaneous) incremental change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
module subroutine plastic_kinehardening_deltaState(Mp,instance,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
me
real(pReal), dimension(param(instance)%sum_N_sl) :: &
gdot_pos,gdot_neg, &
@ -338,29 +339,29 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,of), & ! keep existing...
call kinetics(Mp,instance,me,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,me), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
if (debugConstitutive%extensive &
.and. (of == prm%of_debug .or. .not. debugConstitutive%selective)) then
.and. (me == prm%of_debug .or. .not. debugConstitutive%selective)) then
print*, '======= kinehardening delta state ======='
print*, sense,state(instance)%sense(:,of)
print*, sense,state(instance)%sense(:,me)
endif
#endif
!--------------------------------------------------------------------------------------------------
! switch in sense of shear?
where(dNeq(sense,stt%sense(:,of),0.1_pReal))
dlt%sense (:,of) = sense - stt%sense(:,of) ! switch sense
dlt%chi0 (:,of) = abs(stt%crss_back(:,of)) - stt%chi0(:,of) ! remember current backstress magnitude
dlt%gamma0(:,of) = stt%accshear(:,of) - stt%gamma0(:,of) ! remember current accumulated shear
! switch in sense me shear?
where(dNeq(sense,stt%sense(:,me),0.1_pReal))
dlt%sense (:,me) = sense - stt%sense(:,me) ! switch sense
dlt%chi0 (:,me) = abs(stt%crss_back(:,me)) - stt%chi0(:,me) ! remember current backstress magnitude
dlt%gamma0(:,me) = stt%accshear(:,me) - stt%gamma0(:,me) ! remember current accumulated shear
else where
dlt%sense (:,of) = 0.0_pReal
dlt%chi0 (:,of) = 0.0_pReal
dlt%gamma0(:,of) = 0.0_pReal
dlt%sense (:,me) = 0.0_pReal
dlt%chi0 (:,me) = 0.0_pReal
dlt%gamma0(:,me) = 0.0_pReal
end where
end associate
@ -413,14 +414,14 @@ end subroutine plastic_kinehardening_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,instance,of, &
pure subroutine kinetics(Mp,instance,me, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
me
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
gdot_pos, &
@ -437,21 +438,21 @@ pure subroutine kinetics(Mp,instance,of, &
associate(prm => param(instance), stt => state(instance))
do i = 1, prm%sum_N_sl
tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,me)
tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,me), &
0.0_pReal, prm%nonSchmidActive)
enddo
where(dNeq0(tau_pos))
gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
* sign(abs(tau_pos/stt%crss(:,me))**prm%n, tau_pos)
else where
gdot_pos = 0.0_pReal
end where
where(dNeq0(tau_neg))
gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
* sign(abs(tau_neg/stt%crss(:,me))**prm%n, tau_neg)
else where
gdot_neg = 0.0_pReal
end where
@ -474,4 +475,4 @@ pure subroutine kinetics(Mp,instance,of, &
end subroutine kinetics
end submodule plastic_kinehardening
end submodule kinehardening

8
src/constitutive_plastic_none.f90 → src/phase_mechanics_plastic_none.f90
View File

@ -4,7 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Dummy plasticity for purely elastic material
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_none
submodule(phase:plastic) none
contains
@ -25,7 +25,7 @@ module function plastic_none_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_none init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:none init -+>>>'
phases => config_material%get('phase')
allocate(myPlasticity(phases%length), source = .false.)
@ -43,11 +43,11 @@ module function plastic_none_init() result(myPlasticity)
do p = 1, phases%length
phase => phases%get(p)
if(.not. myPlasticity(p)) cycle
Nconstituents = count(material_phaseAt == p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0)
enddo
end function plastic_none_init
end submodule plastic_none
end submodule none

169
src/constitutive_plastic_nonlocal.f90 → src/phase_mechanics_plastic_nonlocal.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for plasticity including dislocation flux
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_nonlocal
submodule(phase:plastic) nonlocal
use geometry_plastic_nonlocal, only: &
nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
@ -187,7 +187,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:nonlocal init -+>>>'
myPlasticity = plastic_active('nonlocal')
Ninstances = count(myPlasticity)
@ -393,7 +393,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -552,17 +552,16 @@ end function plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
module subroutine nonlocal_dependentState(instance, me, ip, el)
integer, intent(in) :: &
instance, &
of, &
me, &
ip, &
el
integer :: &
ph, &
me, &
no, & !< neighbor offset
neighbor_el, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point
@ -612,9 +611,9 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
associate(prm => param(instance),dst => microstructure(instance), stt => state(instance))
rho = getRho(instance,of,ip,el)
rho = getRho(instance,me,ip,el)
stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+ matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
@ -624,13 +623,13 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
myInteractionMatrix = prm%h_sl_sl &
* spread(( 1.0_pReal - prm%f_F &
+ prm%f_F &
* log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,of),prm%rho_significant))) &
* log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,me),prm%rho_significant))) &
/ log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
else
myInteractionMatrix = prm%h_sl_sl
endif
dst%tau_pass(:,of) = prm%mu * prm%b_sl &
dst%tau_pass(:,me) = prm%mu * prm%b_sl &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities
@ -640,10 +639,9 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
! ToDo: MD: this is most likely only correct for F_i = I
!#################################################################################################
rho0 = getRho0(instance,of,ip,el)
rho0 = getRho0(instance,me,ip,el)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
ph = material_phaseAt(1,el)
me = material_phaseMemberAt(1,ip,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me))
@ -734,7 +732,7 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
! ... gives the local stress correction when multiplied with a factor
dst%tau_back(s,of) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
dst%tau_back(s,me) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
+ rhoExcessGradient_over_rho(2))
enddo
@ -745,29 +743,29 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then
print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip
print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,of)
print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6
print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6
print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,me)
print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,me)*1e-6
print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,me)*1e-6
endif
#endif
end associate
end subroutine plastic_nonlocal_dependentState
end subroutine nonlocal_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,instance,of,ip,el)
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,ph,me,ip,el)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
integer, intent(in) :: &
instance, &
of, &
ph, &
me, &
ip, & !< current integration point
el !< current element number
real(pReal), intent(in) :: &
@ -784,24 +782,25 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
l, &
t, & !< dislocation type
s !< index of my current slip system
real(pReal), dimension(param(instance)%sum_N_sl,8) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: &
rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: &
rho
real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: &
v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate
associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance))
associate(prm => param(phase_plasticityInstance(ph)),dst=>microstructure(phase_plasticityInstance(ph)),&
stt=>state(phase_plasticityInstance(ph)))
ns = prm%sum_N_sl
!*** shortcut to state variables
rho = getRho(instance,of,ip,el)
rho = getRho(phase_plasticityInstance(ph),me,ip,el)
rhoSgl = rho(:,sgl)
do s = 1,ns
@ -816,12 +815,12 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s))
endif
enddo
tauNS = tauNS + spread(dst%tau_back(:,of),2,4)
tau = tau + dst%tau_back(:,of)
tauNS = tauNS + spread(dst%tau_back(:,me),2,4)
tau = tau + dst%tau_back(:,me)
! edges
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
tau, tauNS(:,1), dst%tau_pass(:,of),1,Temperature, instance)
tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticityInstance(ph))
v(:,2) = v(:,1)
dv_dtau(:,2) = dv_dtau(:,1)
dv_dtauNS(:,2) = dv_dtauNS(:,1)
@ -834,11 +833,11 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
else
do t = 3,4
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
tau, tauNS(:,t), dst%tau_pass(:,of),2,Temperature, instance)
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticityInstance(ph))
enddo
endif
stt%v(:,of) = pack(v,.true.)
stt%v(:,me) = pack(v,.true.)
!*** Bauschinger effect
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
@ -861,19 +860,19 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
end associate
end subroutine plastic_nonlocal_LpAndItsTangent
end subroutine nonlocal_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
integer, intent(in) :: &
instance, & ! current instance of this plasticity
of, & !< offset
me, & !< offset
ip, &
el
@ -904,10 +903,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
ns = prm%sum_N_sl
!*** shortcut to state variables
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of)
forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),of)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me)
forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),me)
rho = getRho(instance,of,ip,el)
rho = getRho(instance,me,ip,el)
rhoDip = rho(:,dip)
!****************************************************************************
@ -928,7 +927,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
!*** calculate limits for stable dipole height
do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,of)
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,me)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
@ -952,10 +951,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
/ (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),me) = dUpper(s,c)
plasticState(ph)%deltaState(:,of) = 0.0_pReal
del%rho(:,of) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
plasticState(ph)%deltaState(:,me) = 0.0_pReal
del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
#ifdef DEBUG
if (debugConstitutive%extensive &
@ -974,8 +973,8 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
instance,of,ip,el)
module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
ph,me,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
@ -983,18 +982,17 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
instance, &
of, &
ph, &
me, &
ip, & !< current integration point
el !< current element number
integer :: &
ph, &
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
t, & !< type of dislocation
s !< index of my current slip system
real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: &
rho, &
rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution
@ -1002,45 +1000,44 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(param(instance)%sum_N_sl,8) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity
v0, &
gdot !< shear rates
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
tau, & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles
real(pReal), dimension(param(instance)%sum_N_sl,2) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws
real(pReal) :: &
selfDiffusion !< self diffusion
ph = material_phaseAt(1,el)
if (timestep <= 0.0_pReal) then
plasticState(ph)%dotState = 0.0_pReal
return
endif
associate(prm => param(instance), &
dst => microstructure(instance), &
dot => dotState(instance), &
stt => state(instance))
associate(prm => param(phase_plasticityInstance(ph)), &
dst => microstructure(phase_plasticityInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)), &
stt => state(phase_plasticityInstance(ph)))
ns = prm%sum_N_sl
tau = 0.0_pReal
gdot = 0.0_pReal
rho = getRho(instance,of,ip,el)
rho = getRho(phase_plasticityInstance(ph),me,ip,el)
rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip)
rho0 = getRho0(instance,of,ip,el)
rho0 = getRho0(phase_plasticityInstance(ph),me,ip,el)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticityInstance(ph)),me)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG
@ -1055,7 +1052,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
!****************************************************************************
!*** limits for stable dipole height
do s = 1,ns
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,of)
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,me)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
@ -1076,20 +1073,20 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path
* sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path
* sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall
else isBCC
rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) &
* sqrt(stt%rho_forest(:,of)) / prm%i_sl / prm%b_sl, 2, 4)
* sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticityInstance(ph)),me)
!****************************************************************************
@ -1132,10 +1129,10 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
if (lattice_structure(ph) == LATTICE_fcc_ID) &
forall (s = 1:ns, prm%colinearSystem(s) > 0) &
rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
* 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
* 0.25_pReal * sqrt(stt%rho_forest(s,me)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
!*** thermally activated annihilation of edge dipoles by climb
!*** thermally activated annihilation me edge dipoles by climb
rhoDotThermalAnnihilation = 0.0_pReal
selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))
vClimb = prm%V_at * selfDiffusion * prm%mu &
@ -1145,7 +1142,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, instance,of,ip,el) &
rhoDot = rhoDotFlux(timestep, phase_plasticityInstance(ph),me,ip,el) &
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
@ -1162,25 +1159,25 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
#endif
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
else
dot%rho(:,of) = pack(rhoDot,.true.)
dot%gamma(:,of) = sum(gdot,2)
dot%rho(:,me) = pack(rhoDot,.true.)
dot%gamma(:,me) = sum(gdot,2)
endif
end associate
end subroutine plastic_nonlocal_dotState
end subroutine nonlocal_dotState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
function rhoDotFlux(timestep,instance,of,ip,el)
function rhoDotFlux(timestep,instance,me,ip,el)
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
instance, &
of, &
me, &
ip, & !< current integration point
el !< current element number
@ -1243,16 +1240,16 @@ function rhoDotFlux(timestep,instance,of,ip,el)
gdot = 0.0_pReal
rho = getRho(instance,of,ip,el)
rho = getRho(instance,me,ip,el)
rhoSgl = rho(:,sgl)
rho0 = getRho0(instance,of,ip,el)
rho0 = getRho0(instance,me,ip,el)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !ToDo: MD: I think we should use state0 here
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),me)
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity)
@ -1287,8 +1284,8 @@ function rhoDotFlux(timestep,instance,of,ip,el)
m(1:3,:,3) = -prm%slip_transverse
m(1:3,:,4) = prm%slip_transverse
my_F = constitutive_mech_F(ph)%data(1:3,1:3,of)
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
my_F = constitutive_mech_F(ph)%data(1:3,1:3,me)
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)))
neighbors: do n = 1,nIPneighbors
@ -1789,14 +1786,14 @@ end subroutine kinetics
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho(instance,of,ip,el)
pure function getRho(instance,me,ip,el)
integer, intent(in) :: instance, of,ip,el
integer, intent(in) :: instance, me,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho
associate(prm => param(instance))
getRho = reshape(state(instance)%rho(:,of),[prm%sum_N_sl,10])
getRho = reshape(state(instance)%rho(:,me),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign)
getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
@ -1814,14 +1811,14 @@ end function getRho
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho0(instance,of,ip,el)
pure function getRho0(instance,me,ip,el)
integer, intent(in) :: instance, of,ip,el
integer, intent(in) :: instance, me,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0
associate(prm => param(instance))
getRho0 = reshape(state0(instance)%rho(:,of),[prm%sum_N_sl,10])
getRho0 = reshape(state0(instance)%rho(:,me),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
@ -1834,4 +1831,4 @@ pure function getRho0(instance,of,ip,el)
end function getRho0
end submodule plastic_nonlocal
end submodule nonlocal

79
src/constitutive_plastic_phenopowerlaw.f90 → src/phase_mechanics_plastic_phenopowerlaw.f90
View File

@ -4,7 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_mech) plastic_phenopowerlaw
submodule(phase:plastic) phenopowerlaw
type :: tParameters
real(pReal) :: &
@ -89,7 +89,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
mech, &
pl
print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>'
print'(/,a)', ' <<<+- phase:mechanics:plastic:phenopowerlaw init -+>>>'
myPlasticity = plastic_active('phenopowerlaw')
Ninstances = count(myPlasticity)
@ -225,7 +225,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt == p) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
@ -285,7 +285,7 @@ end function plastic_phenopowerlaw_init
!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume
!--------------------------------------------------------------------------------------------------
pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
@ -295,23 +295,23 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
ph, &
me
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(instance)%sum_N_tw) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: &
gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
associate(prm => param(phase_plasticityInstance(ph)))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -320,7 +320,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
+ dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
call kinetics_twin(Mp,phase_plasticityInstance(ph),me,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%sum_N_tw
Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -330,32 +330,33 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
end subroutine phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
module subroutine phenopowerlaw_dotState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
ph, &
me
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
sumGamma,sumF
real(pReal), dimension(param(instance)%sum_N_sl) :: &
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)))
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)
sumGamma = sum(stt%gamma_slip(:,me))
sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char)
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
@ -367,26 +368,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%h_int
xi_slip_sat_offset = prm%f_sl_sat_tw*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset))
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset))
!--------------------------------------------------------------------------------------------------
! shear rates
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of))
call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,phase_plasticityInstance(ph),me,dot%gamma_twin(:,me))
!--------------------------------------------------------------------------------------------------
! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%h_sl_sl,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_twin(:,of))
dot%xi_slip(:,me) = c_SlipSlip * left_SlipSlip * &
matmul(prm%h_sl_sl,dot%gamma_slip(:,me)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_twin(:,me))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,of))
dot%xi_twin(:,me) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,me)) &
+ c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,me))
end associate
end subroutine plastic_phenopowerlaw_dotState
end subroutine phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
@ -431,14 +432,14 @@ end subroutine plastic_phenopowerlaw_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,of, &
pure subroutine kinetics_slip(Mp,instance,me, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
me
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
gdot_slip_pos, &
@ -462,14 +463,14 @@ pure subroutine kinetics_slip(Mp,instance,of, &
where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_sl, tau_slip_pos)
* sign(abs(tau_slip_pos/stt%xi_slip(:,me))**prm%n_sl, tau_slip_pos)
else where
gdot_slip_pos = 0.0_pReal
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_sl, tau_slip_neg)
* sign(abs(tau_slip_neg/stt%xi_slip(:,me))**prm%n_sl, tau_slip_neg)
else where
gdot_slip_neg = 0.0_pReal
end where
@ -500,14 +501,14 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,instance,of,&
pure subroutine kinetics_twin(Mp,instance,me,&
gdot_twin,dgdot_dtau_twin)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
me
real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: &
gdot_twin
@ -525,8 +526,8 @@ pure subroutine kinetics_twin(Mp,instance,of,&
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction
* prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_tw
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,me)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction
* prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,me))**prm%n_tw
else where
gdot_twin = 0.0_pReal
end where
@ -543,4 +544,4 @@ pure subroutine kinetics_twin(Mp,instance,of,&
end subroutine kinetics_twin
end submodule plastic_phenopowerlaw
end submodule phenopowerlaw

137
src/constitutive_thermal.f90 → src/phase_thermal.f90
View File

@ -1,7 +1,7 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all thermal sources and kinematics constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_thermal
submodule(phase) thermal
enum, bind(c); enumerator :: &
THERMAL_UNDEFINED_ID ,&
@ -18,49 +18,44 @@ submodule(constitutive) constitutive_thermal
type(tDataContainer), dimension(:), allocatable :: current
integer :: thermal_source_maxSizeDotState
interface
module function source_thermal_dissipation_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function source_thermal_dissipation_init
module function dissipation_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function dissipation_init
module function source_thermal_externalheat_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function source_thermal_externalheat_init
module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_thermal_expansion_init
module function externalheat_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
end function externalheat_init
module subroutine source_thermal_externalheat_dotState(ph, me)
module subroutine externalheat_dotState(ph, me)
integer, intent(in) :: &
ph, &
me
end subroutine source_thermal_externalheat_dotState
end subroutine externalheat_dotState
module subroutine dissipation_getRate(TDot, ph,me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(out) :: &
TDot
end subroutine dissipation_getRate
module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase)
integer, intent(in) :: &
phase !< phase ID of element
real(pReal), intent(in), dimension(3,3) :: &
Tstar !< 2nd Piola Kirchhoff stress tensor for a given element
real(pReal), intent(in), dimension(3,3) :: &
Lp !< plastic velocuty gradient for a given element
real(pReal), intent(out) :: &
TDot
end subroutine thermal_dissipation_getRate
module subroutine thermal_externalheat_getRate(TDot, ph,me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(out) :: &
TDot
end subroutine thermal_externalheat_getRate
module subroutine externalheat_getRate(TDot, ph,me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(out) :: &
TDot
end subroutine externalheat_getRate
end interface
@ -82,7 +77,7 @@ module subroutine thermal_init(phases)
Nconstituents
print'(/,a)', ' <<<+- constitutive_thermal init -+>>>'
print'(/,a)', ' <<<+- phase:thermal init -+>>>'
allocate(current(phases%length))
@ -91,7 +86,7 @@ module subroutine thermal_init(phases)
do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
Nconstituents = count(material_phaseAt2 == ph)
allocate(current(ph)%T(Nconstituents),source=300.0_pReal)
allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal)
@ -108,8 +103,8 @@ module subroutine thermal_init(phases)
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
if(maxval(thermal_Nsources) /= 0) then
where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
endif
thermal_source_maxSizeDotState = 0
@ -123,54 +118,40 @@ module subroutine thermal_init(phases)
maxval(thermalState(ph)%p%sizeDotState))
enddo PhaseLoop2
!--------------------------------------------------------------------------------------------------
!initialize kinematic mechanisms
if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) &
phase_kinematics = KINEMATICS_thermal_expansion_ID
end subroutine thermal_init
!----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate
!----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRate(TDot, ip, el)
module subroutine constitutive_thermal_getRate(TDot, ph,me)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer, intent(in) :: ph, me
real(pReal), intent(out) :: &
TDot
real(pReal) :: &
my_Tdot
integer :: &
ph, &
homog, &
me, &
so, &
co
so
homog = material_homogenizationAt(el)
TDot = 0.0_pReal
do co = 1, homogenization_Nconstituents(homog)
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
do so = 1, thermal_Nsources(ph)
select case(thermal_source(so,ph))
case (THERMAL_DISSIPATION_ID)
call thermal_dissipation_getRate(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph)
case (THERMAL_EXTERNALHEAT_ID)
call thermal_externalheat_getRate(my_Tdot, ph,me)
do so = 1, thermal_Nsources(ph)
select case(thermal_source(so,ph))
case (THERMAL_DISSIPATION_ID)
call dissipation_getRate(my_Tdot, ph,me)
case (THERMAL_EXTERNALHEAT_ID)
call externalheat_getRate(my_Tdot, ph,me)
case default
my_Tdot = 0.0_pReal
end select
Tdot = Tdot + my_Tdot
enddo
case default
my_Tdot = 0.0_pReal
end select
Tdot = Tdot + my_Tdot
enddo
enddo
end subroutine constitutive_thermal_getRate
@ -191,7 +172,7 @@ function constitutive_thermal_collectDotState(ph,me) result(broken)
SourceLoop: do i = 1, thermal_Nsources(ph)
if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) &
call source_thermal_externalheat_dotState(ph,me)
call externalheat_dotState(ph,me)
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me)))
@ -206,8 +187,6 @@ module function thermal_stress(Delta_t,ph,me) result(converged_)
integer, intent(in) :: ph, me
logical :: converged_
integer :: so
converged_ = .not. integrateThermalState(Delta_t,ph,me)
@ -268,6 +247,20 @@ module function thermal_T(ph,me) result(T)
end function thermal_T
!----------------------------------------------------------------------------------------------
!< @brief Get rate of temperature (for use by non-thermal physics)
!----------------------------------------------------------------------------------------------
module function thermal_dot_T(ph,me) result(dot_T)
integer, intent(in) :: ph, me
real(pReal) :: dot_T
dot_T = current(ph)%dot_T(me)
end function thermal_dot_T
!----------------------------------------------------------------------------------------------
!< @brief Set temperature
!----------------------------------------------------------------------------------------------
@ -318,4 +311,4 @@ function thermal_active(source_label,src_length) result(active_source)
end function thermal_active
end submodule constitutive_thermal
end submodule thermal

59
src/constitutive_thermal_dissipation.f90 → src/phase_thermal_dissipation.f90
View File

@ -4,11 +4,7 @@
!> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_thermal) source_dissipation
integer, dimension(:), allocatable :: &
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism
submodule(phase:thermal) dissipation
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
@ -25,7 +21,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function source_thermal_dissipation_init(source_length) result(mySources)
module function dissipation_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
@ -35,9 +31,9 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
phase, &
sources, thermal, &
src
integer :: Ninstances,sourceOffset,Nconstituents,p
integer :: Ninstances,so,Nconstituents,ph
print'(/,a)', ' <<<+- thermal_dissipation init -+>>>'
print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>'
mySources = thermal_active('dissipation',source_length)
@ -46,25 +42,21 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
if(Ninstances == 0) return
phases => config_material%get('phase')
allocate(param(Ninstances))
allocate(source_thermal_dissipation_offset (phases%length), source=0)
allocate(source_thermal_dissipation_instance(phases%length), source=0)
allocate(param(phases%length))
do p = 1, phases%length
phase => phases%get(p)
if(any(mySources(:,p))) source_thermal_dissipation_instance(p) = count(mySources(:,1:p))
if(count(mySources(:,p)) == 0) cycle
do ph = 1, phases%length
phase => phases%get(ph)
if(count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1
thermal => phase%get('thermal')
sources => thermal%get('source')
do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then
source_thermal_dissipation_offset(p) = sourceOffset
associate(prm => param(source_thermal_dissipation_instance(p)))
src => sources%get(sourceOffset)
do so = 1, sources%length
if(mySources(so,ph)) then
associate(prm => param(ph))
src => sources%get(so)
prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0)
prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt2 == ph)
call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0)
end associate
endif
@ -72,28 +64,23 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
enddo
end function source_thermal_dissipation_init
end function dissipation_init
!--------------------------------------------------------------------------------------------------
!> @brief Ninstancess dissipation rate
!--------------------------------------------------------------------------------------------------
module subroutine thermal_dissipation_getRate(TDot, Tstar, Lp, phase)
integer, intent(in) :: &
phase
real(pReal), intent(in), dimension(3,3) :: &
Tstar
real(pReal), intent(in), dimension(3,3) :: &
Lp
module subroutine dissipation_getRate(TDot, ph,me)
integer, intent(in) :: ph, me
real(pReal), intent(out) :: &
TDot
associate(prm => param(source_thermal_dissipation_instance(phase)))
TDot = prm%kappa*sum(abs(Tstar*Lp))
associate(prm => param(ph))
TDot = prm%kappa*sum(abs(mech_S(ph,me)*mech_L_p(ph,me)))
end associate
end subroutine thermal_dissipation_getRate
end subroutine dissipation_getRate
end submodule source_dissipation
end submodule dissipation

134
src/phase_thermal_externalheat.f90 Normal file
View File

@ -0,0 +1,134 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Michigan State University
!> @brief material subroutine for variable heat source
!--------------------------------------------------------------------------------------------------
submodule(phase:thermal) externalheat
integer, dimension(:), allocatable :: &
source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism?
type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: &
t_n, &
f_T
integer :: &
nIntervals
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function externalheat_init(source_length) result(mySources)
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: &
phases, &
phase, &
sources, thermal, &
src
integer :: Ninstances,so,Nconstituents,ph
print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>'
mySources = thermal_active('externalheat',source_length)
Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(source_thermal_externalheat_offset (phases%length), source=0)
do ph = 1, phases%length
phase => phases%get(ph)
if(count(mySources(:,ph)) == 0) cycle
thermal => phase%get('thermal')
sources => thermal%get('source')
do so = 1, sources%length
if(mySources(so,ph)) then
source_thermal_externalheat_offset(ph) = so
associate(prm => param(ph))
src => sources%get(so)
prm%t_n = src%get_asFloats('t_n')
prm%nIntervals = size(prm%t_n) - 1
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt2 == ph)
call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0)
end associate
endif
enddo
enddo
end function externalheat_init
!--------------------------------------------------------------------------------------------------
!> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers
!--------------------------------------------------------------------------------------------------
module subroutine externalheat_dotState(ph, me)
integer, intent(in) :: &
ph, &
me
integer :: &
so
so = source_thermal_externalheat_offset(ph)
thermalState(ph)%p(so)%dotState(1,me) = 1.0_pReal ! state is current time
end subroutine externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate
!--------------------------------------------------------------------------------------------------
module subroutine externalheat_getRate(TDot, ph, me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(out) :: &
TDot
integer :: &
so, interval
real(pReal) :: &
frac_time
so = source_thermal_externalheat_offset(ph)
associate(prm => param(ph))
do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (thermalState(ph)%p(so)%state(1,me) - prm%t_n(interval)) &
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside me bounds
enddo
end associate
end subroutine externalheat_getRate
end submodule externalheat