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

Merge branch 'development' into MiscImprovements

This commit is contained in:
Martin Diehl 2020-02-07 17:11:01 +01:00
parent 1ceba73d31 177babecb2
commit e932b386b9
18 changed files with 2704 additions and 2566 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 2077cffcf72afd273751c85dfb77bf7c7f372575 Subproject commit 64432754ce3c590c882cf4987695539cee524ee8

2
VERSION
View File

@ -1 +1 @@
v2.0.3-1531-g3d6ec695 v2.0.3-1624-g47109b90

14
src/commercialFEM_fileList.f90
View File

@ -36,14 +36,14 @@
#include "kinematics_cleavage_opening.f90" #include "kinematics_cleavage_opening.f90"
#include "kinematics_slipplane_opening.f90" #include "kinematics_slipplane_opening.f90"
#include "kinematics_thermal_expansion.f90" #include "kinematics_thermal_expansion.f90"
#include "plastic_none.f90"
#include "plastic_isotropic.f90"
#include "plastic_phenopowerlaw.f90"
#include "plastic_kinematichardening.f90"
#include "plastic_dislotwin.f90"
#include "plastic_disloUCLA.f90"
#include "plastic_nonlocal.f90"
#include "constitutive.f90" #include "constitutive.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_disloUCLA.f90"
#include "constitutive_plastic_nonlocal.f90"
#include "crystallite.f90" #include "crystallite.f90"
#include "thermal_isothermal.f90" #include "thermal_isothermal.f90"
#include "thermal_adiabatic.f90" #include "thermal_adiabatic.f90"

289
src/constitutive.f90
View File

@ -4,7 +4,9 @@
!> @brief elasticity, plasticity, internal microstructure state !> @brief elasticity, plasticity, internal microstructure state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module constitutive module constitutive
use prec
use math use math
use rotations
use debug use debug
use numerics use numerics
use IO use IO
@ -13,13 +15,6 @@ module constitutive
use results use results
use lattice use lattice
use discretization use discretization
use plastic_none
use plastic_isotropic
use plastic_phenopowerlaw
use plastic_kinehardening
use plastic_dislotwin
use plastic_disloucla
use plastic_nonlocal
use geometry_plastic_nonlocal use geometry_plastic_nonlocal
use source_thermal_dissipation use source_thermal_dissipation
use source_thermal_externalheat use source_thermal_externalheat
@ -38,7 +33,284 @@ module constitutive
constitutive_plasticity_maxSizeDotState, & constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState constitutive_source_maxSizeDotState
interface
module subroutine plastic_none_init
end subroutine plastic_none_init
module subroutine plastic_isotropic_init
end subroutine plastic_isotropic_init
module subroutine plastic_phenopowerlaw_init
end subroutine plastic_phenopowerlaw_init
module subroutine plastic_kinehardening_init
end subroutine plastic_kinehardening_init
module subroutine plastic_dislotwin_init
end subroutine plastic_dislotwin_init
module subroutine plastic_disloUCLA_init
end subroutine plastic_disloUCLA_init
module subroutine plastic_nonlocal_init
end subroutine plastic_nonlocal_init
module subroutine plastic_isotropic_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_isotropic_LpAndItsTangent
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
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_disloUCLA_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_disloUCLA_LpAndItsTangent
module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, &
Mp, Temperature, volume, 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, &
volume
integer, intent(in) :: &
ip, & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_LpAndItsTangent
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 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_disloUCLA_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_disloUCLA_dotState
module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
timestep,ip,el)
integer, intent(in) :: &
ip, & !< current integration point
el !< current element number
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
real(pReal), dimension(3,3), intent(in) ::&
Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
Fe, & !< elastic deformation gradient
Fp !< plastic deformation gradient
end subroutine plastic_nonlocal_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_disloUCLA_dependentState(instance,of)
integer, intent(in) :: &
instance, &
of
end subroutine plastic_disloUCLA_dependentState
module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
integer, intent(in) :: &
ip, &
el
real(pReal), dimension(3,3), intent(in) :: &
Fe, &
Fp
end subroutine plastic_nonlocal_dependentState
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 plastic_nonlocal_deltaState(Mp,ip,el)
integer, intent(in) :: &
ip, &
el
real(pReal), dimension(3,3), intent(in) :: &
Mp
end subroutine plastic_nonlocal_deltaState
module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
real(pReal), dimension(6,6) :: &
homogenizedC
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
end function plastic_dislotwin_homogenizedC
module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
integer, intent(in) :: &
i, &
e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_isotropic_results
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_phenopowerlaw_results
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_kinehardening_results
module subroutine plastic_dislotwin_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_dislotwin_results
module subroutine plastic_disloUCLA_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_disloUCLA_results
module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results
end interface
public :: & public :: &
plastic_nonlocal_updateCompatibility, &
constitutive_init, & constitutive_init, &
constitutive_homogenizedC, & constitutive_homogenizedC, &
constitutive_microstructure, & constitutive_microstructure, &
@ -179,8 +451,7 @@ end subroutine constitutive_microstructure
! Mp in, dLp_dMp out ! Mp in, dLp_dMp out
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el) S, Fi, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point

49
src/plastic_disloUCLA.f90 → src/constitutive_plastic_disloUCLA.f90
View File

@ -5,21 +5,9 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief crystal plasticity model for bcc metals, especially Tungsten !> @brief crystal plasticity model for bcc metals, especially Tungsten
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_disloUCLA submodule(constitutive) plastic_disloUCLA
use prec
use debug
use math
use IO
use material
use config
use lattice
use discretization
use results
implicit none real(pReal), parameter :: &
private
real(pReal), parameter, private :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
enum, bind(c) enum, bind(c)
@ -33,7 +21,7 @@ module plastic_disloUCLA
tau_pass_ID tau_pass_ID
end enum end enum
type, private :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
aTol_rho, & aTol_rho, &
D, & !< grain size D, & !< grain size
@ -75,14 +63,14 @@ module plastic_disloUCLA
dipoleFormation !< flag indicating consideration of dipole formation dipoleFormation !< flag indicating consideration of dipole formation
end type !< container type for internal constitutive parameters end type !< container type for internal constitutive parameters
type, private :: tDisloUCLAState type :: tDisloUCLAState
real(pReal), dimension(:,:), pointer :: & real(pReal), dimension(:,:), pointer :: &
rho_mob, & rho_mob, &
rho_dip, & rho_dip, &
gamma_sl gamma_sl
end type tDisloUCLAState end type tDisloUCLAState
type, private :: tDisloUCLAdependentState type :: tDisloUCLAdependentState
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
Lambda_sl, & Lambda_sl, &
threshold_stress threshold_stress
@ -90,20 +78,11 @@ module plastic_disloUCLA
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters and state ! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param type(tParameters), allocatable, dimension(:) :: param
type(tDisloUCLAState), allocatable, dimension(:), private :: & type(tDisloUCLAState), allocatable, dimension(:) :: &
dotState, & dotState, &
state state
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState type(tDisloUCLAdependentState), allocatable, dimension(:) :: dependentState
public :: &
plastic_disloUCLA_init, &
plastic_disloUCLA_dependentState, &
plastic_disloUCLA_LpAndItsTangent, &
plastic_disloUCLA_dotState, &
plastic_disloUCLA_results
private :: &
kinetics
contains contains
@ -112,7 +91,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_init() module subroutine plastic_disloUCLA_init
integer :: & integer :: &
Ninstance, & Ninstance, &
@ -328,7 +307,7 @@ end subroutine plastic_disloUCLA_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
Mp,T,instance,of) Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
@ -371,7 +350,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -431,7 +410,7 @@ end subroutine plastic_disloUCLA_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_dependentState(instance,of) module subroutine plastic_disloUCLA_dependentState(instance,of)
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
@ -457,7 +436,7 @@ end subroutine plastic_disloUCLA_dependentState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_results(instance,group) module subroutine plastic_disloUCLA_results(instance,group)
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -604,4 +583,4 @@ pure subroutine kinetics(Mp,T,instance,of, &
end subroutine kinetics end subroutine kinetics
end module plastic_disloUCLA end submodule plastic_disloUCLA

1156
src/constitutive_plastic_dislotwin.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

32
src/plastic_isotropic.f90 → src/constitutive_plastic_isotropic.f90
View File

@ -7,18 +7,7 @@
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal !! untextured polycrystal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_isotropic submodule(constitutive) plastic_isotropic
use prec
use debug
use math
use IO
use material
use config
use discretization
use results
implicit none
private
enum, bind(c) enum, bind(c)
enumerator :: & enumerator :: &
@ -64,20 +53,13 @@ module plastic_isotropic
dotState, & dotState, &
state state
public :: &
plastic_isotropic_init, &
plastic_isotropic_LpAndItsTangent, &
plastic_isotropic_LiAndItsTangent, &
plastic_isotropic_dotState, &
plastic_isotropic_results
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_init module subroutine plastic_isotropic_init
integer :: & integer :: &
Ninstance, & Ninstance, &
@ -207,7 +189,7 @@ end subroutine plastic_isotropic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
@ -268,7 +250,7 @@ end subroutine plastic_isotropic_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient Li !< inleastic velocity gradient
@ -320,7 +302,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_dotState(Mp,instance,of) module subroutine plastic_isotropic_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -369,7 +351,7 @@ end subroutine plastic_isotropic_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_results(instance,group) module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -388,4 +370,4 @@ subroutine plastic_isotropic_results(instance,group)
end subroutine plastic_isotropic_results end subroutine plastic_isotropic_results
end module plastic_isotropic end submodule plastic_isotropic

522
src/constitutive_plastic_kinehardening.f90 Normal file
View File

@ -0,0 +1,522 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Michigan State University
!> @author Zhuowen Zhao, Michigan State University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates
!! and a Voce-type kinematic hardening rule
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_kinehardening
enum, bind(c)
enumerator :: &
undefined_ID, &
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, &
shearrate_ID, &
resolvedstress_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip
n, & !< stress exponent for slip
aTolResistance, &
aTolShear
real(pReal), allocatable, dimension(:) :: &
crss0, & !< initial critical shear stress for slip
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip
theta0_b, & !< initial hardening rate of back stress for each slip
theta1_b, & !< asymptotic hardening rate of back stress for each slip
tau1, &
tau1_b, &
nonSchmidCoeff
real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
totalNslip, & !< total number of active slip system
of_debug = 0
integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type tParameters
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress
sense, & !< sense of acting shear stress (-1 or +1)
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
end type tKinehardeningState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tKinehardeningState), allocatable, dimension(:) :: &
dotState, &
deltaState, &
state
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_init
integer :: &
Ninstance, &
p, i, o, &
NipcMyPhase, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
do p = 1, size(phase_plasticityInstance)
if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
dlt => deltaState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)),&
config => config_phase(p))
#ifdef DEBUG
if (p==material_phaseAt(debug_g,debug_e)) then
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
endif
#endif
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
! expand: family => system
prm%crss0 = math_expand(prm%crss0, prm%Nslip)
prm%tau1 = math_expand(prm%tau1, prm%Nslip)
prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip)
prm%theta0 = math_expand(prm%theta0, prm%Nslip)
prm%theta1 = math_expand(prm%theta1, prm%Nslip)
prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip)
prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip)
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
!ToDo: Any sensible checks for theta?
endif slipActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance')
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0)
case ('accumulatedshear')
outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0)
case ('shearrate')
outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0)
case ('resolvedstress')
outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0)
case ('backstress')
outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0)
case ('sense')
outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0)
case ('chi0')
outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0)
case ('gamma0')
outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0)
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID , outputID]
endif
enddo
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip
sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
endIndex = prm%totalNslip
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
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
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%totalNslip
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
end subroutine plastic_kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
sumGamma
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
* exp(-sumGamma*prm%theta0/prm%tau1) &
)
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + &
(prm%theta0_b - prm%theta1_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
)
end associate
end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
sense
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...
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 (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
write(6,'(a)') '======= kinehardening delta state ======='
write(6,*) sense,state(instance)%sense(:,of)
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
else where
dlt%sense (:,of) = 0.0_pReal
dlt%chi0 (:,of) = 0.0_pReal
dlt%gamma0(:,of) = 0.0_pReal
end where
end associate
end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
case(crss_back_ID)
call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case (sense_ID)
call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
case (chi0_ID)
call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
case (gamma0_ID)
call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
case (accshear_ID)
call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
end select
enddo outputsLoop
end associate
end subroutine plastic_kinehardening_results
!--------------------------------------------------------------------------------------------------
!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
!> @details: Shear rates are calculated only optionally.
! 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, &
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
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_pos, &
gdot_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_pos, &
dgdot_dtau_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_pos, &
tau_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_pos))
gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
else where
gdot_pos = 0.0_pReal
end where
where(dNeq0(tau_neg))
gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
else where
gdot_neg = 0.0_pReal
end where
if (present(dgdot_dtau_pos)) then
where(dNeq0(gdot_pos))
dgdot_dtau_pos = gdot_pos*prm%n/tau_pos
else where
dgdot_dtau_pos = 0.0_pReal
end where
endif
if (present(dgdot_dtau_neg)) then
where(dNeq0(gdot_neg))
dgdot_dtau_neg = gdot_neg*prm%n/tau_neg
else where
dgdot_dtau_neg = 0.0_pReal
end where
endif
end associate
end subroutine kinetics
end submodule plastic_kinehardening

51
src/plastic_none.f90 → src/constitutive_plastic_none.f90
View File

@ -4,16 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Dummy plasticity for purely elastic material !> @brief Dummy plasticity for purely elastic material
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_none submodule(constitutive) plastic_none
use material
use discretization
use debug
implicit none
private
public :: &
plastic_none_init
contains contains
@ -21,27 +12,27 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_none_init module subroutine plastic_none_init
integer :: & integer :: &
Ninstance, & Ninstance, &
p, & p, &
NipcMyPhase NipcMyPhase
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do p = 1, size(phase_plasticity) do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
call material_allocatePlasticState(p,NipcMyPhase,0,0,0) call material_allocatePlasticState(p,NipcMyPhase,0,0,0)
enddo enddo
end subroutine plastic_none_init end subroutine plastic_none_init
end module plastic_none end submodule plastic_none

56
src/plastic_nonlocal.f90 → src/constitutive_plastic_nonlocal.f90
View File

@ -4,18 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for plasticity including dislocation flux !> @brief material subroutine for plasticity including dislocation flux
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_nonlocal submodule(constitutive) plastic_nonlocal
use prec
use IO
use math
use debug
use material
use lattice
use rotations
use results
use config
use lattice
use discretization
use geometry_plastic_nonlocal, only: & use geometry_plastic_nonlocal, only: &
nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
@ -23,11 +12,9 @@ module plastic_nonlocal
IParea => geometry_plastic_nonlocal_IParea0, & IParea => geometry_plastic_nonlocal_IParea0, &
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0
implicit none
private
real(pReal), parameter :: & real(pReal), parameter :: &
KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin
! storage order of dislocation types ! storage order of dislocation types
integer, dimension(8), parameter :: & integer, dimension(8), parameter :: &
sgl = [1,2,3,4,5,6,7,8] !< signed (single) sgl = [1,2,3,4,5,6,7,8] !< signed (single)
@ -199,25 +186,13 @@ module plastic_nonlocal
type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
public :: &
plastic_nonlocal_init, &
plastic_nonlocal_dependentState, &
plastic_nonlocal_LpAndItsTangent, &
plastic_nonlocal_dotState, &
plastic_nonlocal_deltaState, &
plastic_nonlocal_updateCompatibility, &
plastic_nonlocal_results
private :: &
plastic_nonlocal_kinetics
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_init module subroutine plastic_nonlocal_init
integer :: & integer :: &
sizeState, sizeDotState,sizeDependentState, sizeDeltaState, & sizeState, sizeDotState,sizeDependentState, sizeDeltaState, &
@ -254,7 +229,6 @@ subroutine plastic_nonlocal_init
allocate(dotState(maxNinstances)) allocate(dotState(maxNinstances))
allocate(deltaState(maxNinstances)) allocate(deltaState(maxNinstances))
allocate(microstructure(maxNinstances)) allocate(microstructure(maxNinstances))
allocate(totalNslip(maxNinstances), source=0) allocate(totalNslip(maxNinstances), source=0)
@ -738,7 +712,7 @@ end subroutine plastic_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure !> @brief calculates quantities characterizing the microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & ip, &
@ -1100,8 +1074,8 @@ end subroutine plastic_nonlocal_kinetics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, &
Mp, Temperature, volume, ip, el) Mp, Temperature, volume, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< current integration point ip, & !< current integration point
@ -1230,7 +1204,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure !> @brief (instantaneous) incremental change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_deltaState(Mp,ip,el) module subroutine plastic_nonlocal_deltaState(Mp,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & ip, &
@ -1346,8 +1320,8 @@ end subroutine plastic_nonlocal_deltaState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
timestep,ip,el) timestep,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< current integration point ip, & !< current integration point
@ -1795,13 +1769,13 @@ end subroutine plastic_nonlocal_dotState
! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! plane normals and signed cosine of the angle between the slip directions. Only the largest values
! that sum up to a total of 1 are considered, all others are set to zero. ! that sum up to a total of 1 are considered, all others are set to zero.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
integer, intent(in) :: & integer, intent(in) :: &
i, & i, &
e e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
orientation ! crystal orientation in quaternions orientation ! crystal orientation
integer :: & integer :: &
Nneighbors, & ! number of neighbors Nneighbors, & ! number of neighbors
@ -1960,10 +1934,10 @@ end function getRho
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_results(instance,group) module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*) :: group character(len=*),intent(in) :: group
integer :: o integer :: o
associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance))
@ -2026,4 +2000,4 @@ subroutine plastic_nonlocal_results(instance,group)
end subroutine plastic_nonlocal_results end subroutine plastic_nonlocal_results
end module plastic_nonlocal end submodule plastic_nonlocal

596
src/constitutive_plastic_phenopowerlaw.f90 Normal file
View File

@ -0,0 +1,596 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_phenopowerlaw
enum, bind(c)
enumerator :: &
undefined_ID, &
resistance_slip_ID, &
accumulatedshear_slip_ID, &
shearrate_slip_ID, &
resolvedstress_slip_ID, &
resistance_twin_ID, &
accumulatedshear_twin_ID, &
shearrate_twin_ID, &
resolvedstress_twin_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip
gdot0_twin, & !< reference shear strain rate for twin
n_slip, & !< stress exponent for slip
n_twin, & !< stress exponent for twin
spr, & !< push-up factor for slip saturation due to twinning
c_1, &
c_2, &
c_3, &
c_4, &
h0_SlipSlip, & !< reference hardening slip - slip
h0_TwinSlip, & !< reference hardening twin - slip
h0_TwinTwin, & !< reference hardening twin - twin
a_slip, &
aTolResistance, & !< absolute tolerance for integration of xi
aTolShear, & !< absolute tolerance for integration of gamma
aTolTwinfrac !< absolute tolerance for integration of f
real(pReal), allocatable, dimension(:) :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
xi_slip_sat, & !< maximum critical shear stress for slip
nonSchmidCoeff, &
H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity
interaction_SlipTwin, & !< slip resistance from twin activity
interaction_TwinSlip, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid_slip, &
Schmid_twin, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
totalNslip, & !< total number of active slip system
totalNtwin !< total number of active twin systems
integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type tParameters
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
xi_twin, &
gamma_slip, &
gamma_twin
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, &
state
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_init
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, outputSize, &
sizeState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal)
prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal)
prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal)
prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal)
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else
allocate(prm%nonSchmidCoeff(0))
prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip)
prm%H_int = math_expand(prm%H_int, prm%Nslip)
! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip'
if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0'
if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat'
else slipActive
allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0))
endif slipActive
!--------------------------------------------------------------------------------------------------
! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
prm%gdot0_twin = config%getFloat('gdot0_twin')
prm%n_twin = config%getFloat('n_twin')
prm%spr = config%getFloat('s_pr')
prm%h0_TwinTwin = config%getFloat('h0_twintwin')
! expand: family => system
prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin)
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin'
else twinActive
allocate(prm%interaction_TwinTwin(0,0))
allocate(prm%xi_twin_0(0))
allocate(prm%gamma_twin_char(0))
endif twinActive
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
prm%h0_TwinSlip = config%getFloat('h0_twinslip')
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), &
config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), &
config%getString('lattice_structure'))
else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance_slip')
outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('accumulatedshear_slip')
outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('shearrate_slip')
outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resolvedstress_slip')
outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resistance_twin')
outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('accumulatedshear_twin')
outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('shearrate_twin')
outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('resolvedstress_twin')
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip &
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
endIndex = prm%totalNslip
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @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)
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
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(instance)%totalNtwin) :: &
gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%totalNslip
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%totalNtwin
Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
enddo twinSystems
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
sumGamma,sumF
real(pReal), dimension(param(instance)%totalNslip) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int
xi_slip_sat_offset = prm%spr*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+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))
!--------------------------------------------------------------------------------------------------
! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
end associate
end subroutine plastic_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (resistance_slip_ID)
call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case (accumulatedshear_slip_ID)
call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case (resistance_twin_ID)
call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case (accumulatedshear_twin_ID)
call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
end select
enddo outputsLoop
end associate
end subroutine plastic_phenopowerlaw_results
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
!> @details Derivatives are calculated only optionally.
! 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, &
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
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos, &
gdot_slip_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_slip_pos, &
dgdot_dtau_slip_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_slip_pos, &
tau_slip_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
else where
gdot_slip_pos = 0.0_pReal
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg)
else where
gdot_slip_neg = 0.0_pReal
end where
if (present(dgdot_dtau_slip_pos)) then
where(dNeq0(gdot_slip_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
else where
dgdot_dtau_slip_pos = 0.0_pReal
end where
endif
if (present(dgdot_dtau_slip_neg)) then
where(dNeq0(gdot_slip_neg))
dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg
else where
dgdot_dtau_slip_neg = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
! twinning is assumed to take place only in untwinned volume.
!> @details Derivates are calculated only optionally.
! 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,&
gdot_twin,dgdot_dtau_twin)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
gdot_twin
real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
dgdot_dtau_twin
real(pReal), dimension(param(instance)%totalNtwin) :: &
tau_twin
integer :: i
associate(prm => param(instance), stt => state(instance))
do i = 1, prm%totalNtwin
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
else where
gdot_twin = 0.0_pReal
end where
if (present(dgdot_dtau_twin)) then
where(dNeq0(gdot_twin))
dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
else where
dgdot_dtau_twin = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_twin
end submodule plastic_phenopowerlaw

1
src/crystallite.f90
View File

@ -21,7 +21,6 @@ module crystallite
use constitutive use constitutive
use discretization use discretization
use lattice use lattice
use plastic_nonlocal
use results use results
implicit none implicit none

141
src/element.f90
View File

@ -42,38 +42,39 @@ module element
integer, dimension(NELEMTYPE), parameter :: NNODE = & integer, dimension(NELEMTYPE), parameter :: NNODE = &
[ & [ &
3, & ! 2D 3node 1ip 3, & ! 2D, 1 IP
6, & ! 2D 6node 3ip 6, & ! 2D, 3 IP
4, & ! 2D 4node 4ip 4, & ! 2D, 4 IP
8, & ! 2D 8node 9ip 8, & ! 2D, 9 IP
8, & ! 2D 8node 4ip 8, & ! 2D, 4 IP
!-------------------- !----------------------
4, & ! 3D 4node 1ip 4, & ! 3D, 1 IP
5, & ! 3D 5node 4ip 5, & ! 3D, 4 IP
10, & ! 3D 10node 4ip 10, & ! 3D, 4 IP
6, & ! 3D 6node 6ip 6, & ! 3D, 6 IP
8, & ! 3D 8node 1ip 8, & ! 3D, 1 IP
8, & ! 3D 8node 8ip 8, & ! 3D, 8 IP
20, & ! 3D 20node 8ip 20, & ! 3D, 8 IP
20 & ! 3D 20node 27ip 20 & ! 3D, 27 IP
] !< number of nodes that constitute a specific type of element ] !< number of nodes that constitute a specific type of element
integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = & integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = &
[ & [ &
1, & 1, & ! 1 triangle
2, & 2, & ! 3 quadrilaterals
3, & 3, & ! 4 quadrilaterals
4, & 4, & ! 9 quadrilaterals
3, & 3, & ! 4 quadrilaterals
5, & !----------------------
6, & 5, & ! 1 tetrahedron
6, & 6, & ! 4 hexahedrons
7, & 6, & ! 4 hexahedrons
8, & 7, & ! 6 hexahedrons
9, & 8, & ! 1 hexahedron
9, & 9, & ! 8 hexahedrons
10 & 9, & ! 8 hexahedrons
] !< geometry type of particular element type 10 & ! 27 hexahedrons
] !< geometry type (same number of cell nodes and IPs)
integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = & integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = &
[ & [ &
@ -87,7 +88,7 @@ module element
8, & 8, &
27, & 27, &
64 & 64 &
] !< number of cell nodes in a specific geometry type ] !< number of cell nodes
integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = & integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = &
[ & [ &
@ -101,45 +102,45 @@ module element
1, & 1, &
8, & 8, &
27 & 27 &
] !< number of IPs in a specific geometry type ] !< number of IPs
integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = & integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = &
[ & [ &
1, & ! 2D 3node 1, & ! 2D, 3 node (Triangle)
2, & ! 2D 4node 2, & ! 2D, 4 node (Quadrilateral)
2, & ! 2D 4node 2, & ! - " -
2, & ! 2D 4node 2, & ! - " -
3, & ! 3D 4node 3, & ! 3D, 4 node (Tetrahedron)
4, & ! 3D 8node 4, & ! 3D, 4 node (Hexahedron)
4, & ! 3D 8node 4, & ! - " -
4, & ! 3D 8node 4, & ! - " -
4, & ! 3D 8node 4, & ! - " -
4 & ! 3D 8node 4 & ! - " -
] !< cell type that is used by each geometry type ] !< cell type
integer, dimension(maxval(CELLTYPE)), parameter :: NIPNEIGHBOR = & integer, dimension(maxval(CELLTYPE)), parameter :: NIPNEIGHBOR = &
[ & [ &
3, & ! 2D 3node 3, &
4, & ! 2D 4node 4, &
4, & ! 3D 4node 4, &
6 & ! 3D 8node 6 &
] !< number of ip neighbors / cell faces in a specific cell type ] !< number of ip neighbors / cell faces
integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELLFACE = & integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELLFACE = &
[ & [ &
2, & ! 2D 3node 2, &
2, & ! 2D 4node 2, &
3, & ! 3D 4node 3, &
4 & ! 3D 8node 4 &
] !< number of cell nodes in a specific cell type ] !< number of cell nodes per face
integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELL = & integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELL = &
[ & [ &
3, & ! 2D 3node 3, &
4, & ! 2D 4node 4, &
4, & ! 3D 4node 4, &
8 & ! 3D 8node 8 &
] !< number of cell nodes in a specific cell type ] !< number of total cell nodes
! *** IPneighbor *** ! *** IPneighbor ***
! list of the neighborhood of each IP. ! list of the neighborhood of each IP.
@ -493,14 +494,14 @@ module element
0, 1, 1, 0, 0, 0, 0, 0, & ! 10 0, 1, 1, 0, 0, 0, 0, 0, & ! 10
0, 0, 1, 1, 0, 0, 0, 0, & ! 0, 0, 1, 1, 0, 0, 0, 0, & !
1, 0, 0, 1, 0, 0, 0, 0, & ! 1, 0, 0, 1, 0, 0, 0, 0, & !
1, 0, 0, 0, 1, 0, 0, 0, & !
0, 1, 0, 0, 0, 1, 0, 0, & !
0, 0, 1, 0, 0, 0, 1, 0, & ! 15
0, 0, 0, 1, 0, 0, 0, 1, & !
0, 0, 0, 0, 1, 1, 0, 0, & ! 0, 0, 0, 0, 1, 1, 0, 0, & !
0, 0, 0, 0, 0, 1, 1, 0, & ! 0, 0, 0, 0, 0, 1, 1, 0, & !
0, 0, 0, 0, 0, 0, 1, 1, & ! 0, 0, 0, 0, 0, 0, 1, 1, & ! 15
0, 0, 0, 0, 1, 0, 0, 1, & ! 20 0, 0, 0, 0, 1, 0, 0, 1, & !
1, 0, 0, 0, 1, 0, 0, 0, & !
0, 1, 0, 0, 0, 1, 0, 0, & !
0, 0, 1, 0, 0, 0, 1, 0, & !
0, 0, 0, 1, 0, 0, 0, 1, & ! 20
1, 1, 1, 1, 0, 0, 0, 0, & ! 1, 1, 1, 1, 0, 0, 0, 0, & !
1, 1, 0, 0, 1, 1, 0, 0, & ! 1, 1, 0, 0, 1, 1, 0, 0, & !
0, 1, 1, 0, 0, 1, 1, 0, & ! 0, 1, 1, 0, 0, 1, 1, 0, & !
@ -717,14 +718,14 @@ module element
integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = & integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = &
reshape([& reshape([&
1, 9,21,12,13,22,27,25, & 1, 9,21,12,17,22,27,25, &
9, 2,10,21,22,14,23,27, & 9, 2,10,21,22,18,23,27, &
12,21,11, 4,25,27,24,16, & 12,21,11, 4,25,27,24,20, &
21,10, 3,11,27,23,15,24, & 21,10, 3,11,27,23,19,24, &
13,22,27,25, 5,17,26,20, & 17,22,27,25, 5,13,26,16, &
22,14,23,27,17, 6,18,26, & 22,18,23,27,13, 6,14,26, &
25,27,24,16,20,26,19, 8, & 25,27,24,20,16,26,15, 8, &
27,23,15,24,26,18, 7,19 & 27,23,19,24,26,14, 7,15 &
#if !defined(__GFORTRAN__) #if !defined(__GFORTRAN__)
],shape(CELL9)) ],shape(CELL9))
#else #else

2
src/homogenization.f90
View File

@ -160,7 +160,7 @@ subroutine homogenization_init
allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.) allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.)
allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.) allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.)
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6)
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)

27
src/mesh_marc.f90
View File

@ -585,26 +585,23 @@ subroutine inputRead_elemType(elem, &
select case (IO_lc(what)) select case (IO_lc(what))
case ( '6') case ( '6')
mapElemtype = 1 ! Two-dimensional Plane Strain Triangle mapElemtype = 1 ! Two-dimensional Plane Strain Triangle
case ( '155', & case ( '125') ! 155, 128 (need test)
'125', &
'128')
mapElemtype = 2 ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric) mapElemtype = 2 ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric)
case ( '11') !case ( '11') ! need test
mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain ! mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain
case ( '27') case ( '27')
mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral
case ( '54') case ( '54')
mapElemtype = 5 ! Plane Strain, Eight-node Distorted Quadrilateral with reduced integration mapElemtype = 5 ! Plane Strain, Eight-node Distorted Quadrilateral with reduced integration
case ( '134') !case ( '134') ! need test
mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron ! mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron
case ( '157') !case ( '157') ! need test
mapElemtype = 7 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations ! mapElemtype = 7 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations
case ( '127') !case ( '127') ! need test
mapElemtype = 8 ! Three-dimensional Ten-node Tetrahedron ! mapElemtype = 8 ! Three-dimensional Ten-node Tetrahedron
case ( '136') !case ( '136') ! need test
mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral ! mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral
case ( '117', & case ( '117') ! 123 (need test)
'123')
mapElemtype = 10 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration mapElemtype = 10 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration
case ( '7') case ( '7')
mapElemtype = 11 ! Three-dimensional Arbitrarily Distorted Brick mapElemtype = 11 ! Three-dimensional Arbitrarily Distorted Brick

1175
src/plastic_dislotwin.f90
View File

File diff suppressed because it is too large Load Diff

541
src/plastic_kinematichardening.f90
View File

@ -1,541 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Michigan State University
!> @author Zhuowen Zhao, Michigan State University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates
!! and a Voce-type kinematic hardening rule
!--------------------------------------------------------------------------------------------------
module plastic_kinehardening
use prec
use debug
use math
use IO
use material
use config
use lattice
use discretization
use results
implicit none
private
enum, bind(c)
enumerator :: &
undefined_ID, &
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, &
shearrate_ID, &
resolvedstress_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip
n, & !< stress exponent for slip
aTolResistance, &
aTolShear
real(pReal), allocatable, dimension(:) :: &
crss0, & !< initial critical shear stress for slip
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip
theta0_b, & !< initial hardening rate of back stress for each slip
theta1_b, & !< asymptotic hardening rate of back stress for each slip
tau1, &
tau1_b, &
nonSchmidCoeff
real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
totalNslip, & !< total number of active slip system
of_debug = 0
integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type tParameters
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress
sense, & !< sense of acting shear stress (-1 or +1)
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
end type tKinehardeningState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tKinehardeningState), allocatable, dimension(:) :: &
dotState, &
deltaState, &
state
public :: &
plastic_kinehardening_init, &
plastic_kinehardening_LpAndItsTangent, &
plastic_kinehardening_dotState, &
plastic_kinehardening_deltaState, &
plastic_kinehardening_results
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_init
integer :: &
Ninstance, &
p, i, o, &
NipcMyPhase, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
do p = 1, size(phase_plasticityInstance)
if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
dlt => deltaState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)),&
config => config_phase(p))
#ifdef DEBUG
if (p==material_phaseAt(debug_g,debug_e)) then
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
endif
#endif
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
! expand: family => system
prm%crss0 = math_expand(prm%crss0, prm%Nslip)
prm%tau1 = math_expand(prm%tau1, prm%Nslip)
prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip)
prm%theta0 = math_expand(prm%theta0, prm%Nslip)
prm%theta1 = math_expand(prm%theta1, prm%Nslip)
prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip)
prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip)
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
!ToDo: Any sensible checks for theta?
endif slipActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance')
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0)
case ('accumulatedshear')
outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0)
case ('shearrate')
outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0)
case ('resolvedstress')
outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0)
case ('backstress')
outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0)
case ('sense')
outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0)
case ('chi0')
outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0)
case ('gamma0')
outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0)
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID , outputID]
endif
enddo
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip
sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
endIndex = prm%totalNslip
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
pure 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
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%totalNslip
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
end subroutine plastic_kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
sumGamma
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
* exp(-sumGamma*prm%theta0/prm%tau1) &
)
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + &
(prm%theta0_b - prm%theta1_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
)
end associate
end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_deltaState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
sense
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...
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 (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
write(6,'(a)') '======= kinehardening delta state ======='
write(6,*) sense,state(instance)%sense(:,of)
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
else where
dlt%sense (:,of) = 0.0_pReal
dlt%chi0 (:,of) = 0.0_pReal
dlt%gamma0(:,of) = 0.0_pReal
end where
end associate
end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
case(crss_back_ID)
call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case (sense_ID)
call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
case (chi0_ID)
call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
case (gamma0_ID)
call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
case (accshear_ID)
call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
end select
enddo outputsLoop
end associate
end subroutine plastic_kinehardening_results
!--------------------------------------------------------------------------------------------------
!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
!> @details: Shear rates are calculated only optionally.
! 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, &
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
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_pos, &
gdot_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_pos, &
dgdot_dtau_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_pos, &
tau_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_pos))
gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
else where
gdot_pos = 0.0_pReal
end where
where(dNeq0(tau_neg))
gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
else where
gdot_neg = 0.0_pReal
end where
if (present(dgdot_dtau_pos)) then
where(dNeq0(gdot_pos))
dgdot_dtau_pos = gdot_pos*prm%n/tau_pos
else where
dgdot_dtau_pos = 0.0_pReal
end where
endif
if (present(dgdot_dtau_neg)) then
where(dNeq0(gdot_neg))
dgdot_dtau_neg = gdot_neg*prm%n/tau_neg
else where
dgdot_dtau_neg = 0.0_pReal
end where
endif
end associate
end subroutine kinetics
end module plastic_kinehardening

614
src/plastic_phenopowerlaw.f90
View File

@ -1,614 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!--------------------------------------------------------------------------------------------------
module plastic_phenopowerlaw
use prec
use debug
use math
use IO
use material
use config
use lattice
use discretization
use results
implicit none
private
enum, bind(c)
enumerator :: &
undefined_ID, &
resistance_slip_ID, &
accumulatedshear_slip_ID, &
shearrate_slip_ID, &
resolvedstress_slip_ID, &
resistance_twin_ID, &
accumulatedshear_twin_ID, &
shearrate_twin_ID, &
resolvedstress_twin_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip
gdot0_twin, & !< reference shear strain rate for twin
n_slip, & !< stress exponent for slip
n_twin, & !< stress exponent for twin
spr, & !< push-up factor for slip saturation due to twinning
twinB, &
twinC, &
twinD, &
twinE, &
h0_SlipSlip, & !< reference hardening slip - slip
h0_TwinSlip, & !< reference hardening twin - slip
h0_TwinTwin, & !< reference hardening twin - twin
a_slip, &
aTolResistance, & !< absolute tolerance for integration of xi
aTolShear, & !< absolute tolerance for integration of gamma
aTolTwinfrac !< absolute tolerance for integration of f
real(pReal), allocatable, dimension(:) :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
xi_slip_sat, & !< maximum critical shear stress for slip
nonSchmidCoeff, &
H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity
interaction_SlipTwin, & !< slip resistance from twin activity
interaction_TwinSlip, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid_slip, &
Schmid_twin, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
totalNslip, & !< total number of active slip system
totalNtwin !< total number of active twin systems
integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type tParameters
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
xi_twin, &
gamma_slip, &
gamma_twin
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, &
state
public :: &
plastic_phenopowerlaw_init, &
plastic_phenopowerlaw_LpAndItsTangent, &
plastic_phenopowerlaw_dotState, &
plastic_phenopowerlaw_results
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_init
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, outputSize, &
sizeState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal)
prm%twinC = config%getFloat('twin_c',defaultVal=0.0_pReal)
prm%twinD = config%getFloat('twin_d',defaultVal=0.0_pReal)
prm%twinE = config%getFloat('twin_e',defaultVal=0.0_pReal)
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1)
else
allocate(prm%nonSchmidCoeff(0))
prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip)
prm%H_int = math_expand(prm%H_int, prm%Nslip)
! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip'
if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0'
if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat'
else slipActive
allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0))
endif slipActive
!--------------------------------------------------------------------------------------------------
! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
prm%gdot0_twin = config%getFloat('gdot0_twin')
prm%n_twin = config%getFloat('n_twin')
prm%spr = config%getFloat('s_pr')
prm%h0_TwinTwin = config%getFloat('h0_twintwin')
! expand: family => system
prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin)
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin'
else twinActive
allocate(prm%interaction_TwinTwin(0,0))
allocate(prm%xi_twin_0(0))
allocate(prm%gamma_twin_char(0))
endif twinActive
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
prm%h0_TwinSlip = config%getFloat('h0_twinslip')
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), &
config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), &
config%getString('lattice_structure'))
else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance_slip')
outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('accumulatedshear_slip')
outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('shearrate_slip')
outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resolvedstress_slip')
outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resistance_twin')
outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('accumulatedshear_twin')
outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('shearrate_twin')
outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('resolvedstress_twin')
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip &
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
endIndex = prm%totalNslip
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume
!--------------------------------------------------------------------------------------------------
pure 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
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(instance)%totalNtwin) :: &
gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%totalNslip
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%totalNtwin
Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
enddo twinSystems
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
sumGamma,sumF
real(pReal), dimension(param(instance)%totalNslip) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int
xi_slip_sat_offset = prm%spr*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+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))
!--------------------------------------------------------------------------------------------------
! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
end associate
end subroutine plastic_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (resistance_slip_ID)
call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case (accumulatedshear_slip_ID)
call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case (resistance_twin_ID)
call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case (accumulatedshear_twin_ID)
call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
end select
enddo outputsLoop
end associate
end subroutine plastic_phenopowerlaw_results
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
!> @details Derivatives are calculated only optionally.
! 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, &
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
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos, &
gdot_slip_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_slip_pos, &
dgdot_dtau_slip_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_slip_pos, &
tau_slip_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
else where
gdot_slip_pos = 0.0_pReal
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg)
else where
gdot_slip_neg = 0.0_pReal
end where
if (present(dgdot_dtau_slip_pos)) then
where(dNeq0(gdot_slip_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
else where
dgdot_dtau_slip_pos = 0.0_pReal
end where
endif
if (present(dgdot_dtau_slip_neg)) then
where(dNeq0(gdot_slip_neg))
dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg
else where
dgdot_dtau_slip_neg = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
! twinning is assumed to take place only in untwinned volume.
!> @details Derivates are calculated only optionally.
! 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,&
gdot_twin,dgdot_dtau_twin)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
gdot_twin
real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
dgdot_dtau_twin
real(pReal), dimension(param(instance)%totalNtwin) :: &
tau_twin
integer :: i
associate(prm => param(instance), stt => state(instance))
do i = 1, prm%totalNtwin
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
else where
gdot_twin = 0.0_pReal
end where
if (present(dgdot_dtau_twin)) then
where(dNeq0(gdot_twin))
dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
else where
dgdot_dtau_twin = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_twin
end module plastic_phenopowerlaw