following best practise from phenopowerlaw

This commit is contained in:
Martin Diehl 2019-01-06 21:18:35 +01:00
parent f354f04f49
commit 15d1789a19
2 changed files with 164 additions and 190 deletions

View File

@ -3,52 +3,41 @@
!> @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
!! and a Voce-type kinematic hardening rule
!--------------------------------------------------------------------------------------------------
module plastic_kinehardening
use prec, only: &
pReal,&
pInt
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_kinehardening_Noutput !< number of outputs per instance
integer(pInt), dimension(:,:), allocatable, private :: &
plastic_kinehardening_Nslip !< active number of slip systems per family (input parameter, per family)
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?)
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, private :: tParameters !< container type for internal constitutive parameters
type, private :: tParameters
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip (input parameter)
n_slip, & !< stress exponent for slip (input parameter)
aTolResistance, &
aTolShear
real(pReal), dimension(:), allocatable, private :: &
crss0, & !< initial critical shear stress for slip (input parameter, per family)
real(pReal), allocatable, dimension(:) :: &
crss0, & !< initial critical shear stress for slip (input parameter, per family)
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 >
@ -56,8 +45,8 @@ module plastic_kinehardening
tau1, &
tau1_b, &
nonSchmidCoeff
real(pReal), dimension(:,:), allocatable, private :: &
interaction_slipslip !< latent hardening matrix
real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid_slip, &
Schmid_twin, &
@ -72,6 +61,8 @@ module plastic_kinehardening
outputID !< ID of each post result output
end type
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type, private :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
@ -80,18 +71,13 @@ module plastic_kinehardening
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
end type
type(tParameters), dimension(:), allocatable, private :: &
param !< containers of constitutive parameters (len Ninstance)
type(tKinehardeningState), allocatable, dimension(:), private :: &
dotState, &
deltaState, &
state
public :: &
plastic_kinehardening_init, &
plastic_kinehardening_LpAndItsTangent, &
@ -105,15 +91,19 @@ module plastic_kinehardening
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use prec, only: &
dEq0
dEq0, &
pStringLen
use debug, only: &
#ifdef DEBUG
debug_e, &
@ -142,68 +132,55 @@ subroutine plastic_kinehardening_init
material_phase, &
plasticState
use config, only: &
config_phase, &
MATERIAL_partPhase
MATERIAL_partPhase, &
config_phase
use lattice
implicit none
integer(kind(undefined_ID)) :: &
output_ID
integer(pInt) :: &
o, i, p, &
instance, &
Ninstance, &
NipcMyPhase, &
outputSize, &
offset_slip, &
startIndex, endIndex, &
sizeDotState, &
sizeState, &
sizeDeltaState
p, i, o, &
NipcMyPhase, outputSize, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID !< ID of each post result output
character(len=65536), dimension(:), allocatable :: &
outputs
character(len=65536) :: &
extmsg = '', &
structure = ''
outputID
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
Ninstance = int(count(phase_plasticity == PLASTICITY_KINEHARDENING_ID),pInt)
if (Ninstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a,1x,i5,/)') '# instances:',Ninstance
allocate(plastic_kinehardening_sizePostResult(maxval(phase_Noutput),Ninstance), &
source=0_pInt)
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(plastic_kinehardening_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(plastic_kinehardening_output(maxval(phase_Noutput),Ninstance))
plastic_kinehardening_output = ''
allocate(plastic_kinehardening_Noutput(Ninstance), source=0_pInt)
allocate(plastic_kinehardening_Nslip(lattice_maxNslipFamily,Ninstance), source=0_pInt)
allocate(param(Ninstance)) ! one container of parameters per instance
plastic_kinehardening_output = ''
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
do p = 1_pInt, size(phase_plasticityInstance)
if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
instance = phase_plasticityInstance(p) ! which instance of my phase
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
delta => deltaState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)))
stt => state(phase_plasticityInstance(p)),&
config => config_phase(p))
#ifdef DEBUG
if (p==material_phase(debug_g,debug_i,debug_e)) then
@ -211,11 +188,12 @@ subroutine plastic_kinehardening_init
endif
#endif
structure = config_phase(p)%getString('lattice_structure')
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config_phase(p)%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
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 '
@ -223,13 +201,13 @@ subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray)
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',&
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
@ -238,16 +216,16 @@ subroutine plastic_kinehardening_init
prm%nonSchmid_neg = prm%Schmid_slip
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config_phase(p)%getFloats('interaction_slipslip'), &
config%getFloats('interaction_slipslip'), &
structure(1:3))
prm%crss0 = config_phase(p)%getFloats('crss0', requiredShape=shape(prm%Nslip))
prm%tau1 = config_phase(p)%getFloats('tau1', requiredShape=shape(prm%Nslip))
prm%tau1_b = config_phase(p)%getFloats('tau1_b', requiredShape=shape(prm%Nslip))
prm%theta0 = config_phase(p)%getFloats('theta0', requiredShape=shape(prm%Nslip))
prm%theta1 = config_phase(p)%getFloats('theta1', requiredShape=shape(prm%Nslip))
prm%theta0_b = config_phase(p)%getFloats('theta0_b', requiredShape=shape(prm%Nslip))
prm%theta1_b = config_phase(p)%getFloats('theta1_b', requiredShape=shape(prm%Nslip))
prm%crss0 = config%getFloats('crss0', requiredShape=shape(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredShape=shape(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredShape=shape(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredShape=shape(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredShape=shape(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredShape=shape(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredShape=shape(prm%Nslip))
! expand: family => system
prm%crss0 = math_expand(prm%crss0, prm%Nslip)
@ -258,15 +236,27 @@ subroutine plastic_kinehardening_init
prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip)
prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip)
prm%gdot0 = config_phase(p)%getFloat('gdot0')
prm%n_slip = config_phase(p)%getFloat('n_slip')
prm%gdot0 = config%getFloat('gdot0')
prm%n_slip = config%getFloat('n_slip')
!--------------------------------------------------------------------------------------------------
! sanity checks
! if (any(plastic_kinehardening_Nslip (1:nSlipFamilies,instance) > 0_pInt &
! .and. param(instance)%crss0 (1:nSlipFamilies) < 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
! if (any(plastic_kinehardening_Nslip (1:nSlipFamilies,instance) > 0_pInt &
! .and. param(instance)%tau1 (1:nSlipFamilies) <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
! if (any(plastic_kinehardening_Nslip (1:nSlipFamilies,instance) > 0_pInt &
! .and. param(instance)%tau1_b(1:nSlipFamilies) < 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
! if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
! if (param(instance)%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
endif slipActive
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray)
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
outputID = undefined_ID
@ -299,7 +289,6 @@ subroutine plastic_kinehardening_init
end select
if (outputID /= undefined_ID) then
plastic_kinehardening_Noutput(instance) = plastic_kinehardening_Noutput(instance) + 1_pInt
plastic_kinehardening_output(i,phase_plasticityInstance(p)) = outputs(i)
plastic_kinehardening_sizePostResult(i,phase_plasticityInstance(p)) = outputSize
prm%outputID = [prm%outputID , outputID]
@ -309,90 +298,71 @@ subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p) ! number of constituents with my phase
NipcMyPhase = count(material_phase == p)
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, &
prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_kinehardening_sizePostResult(:,phase_plasticityInstance(p)))
plasticState(p)%offsetDeltaState = sizeDotState
startIndex = 1_pInt
endIndex = prm%totalNslip
stt%crss => plasticState(p)%state (startIndex:endIndex,1:NipcMyPhase)
dot%crss => plasticState(p)%dotState (startIndex:endIndex,1:NipcMyPhase)
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,1:NipcMyPhase)
dot%crss_back => plasticState(p)%dotState (startIndex:endIndex,1:NipcMyPhase)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%accshear => plasticState(p)%state (startIndex:endIndex,1:NipcMyPhase)
dot%accshear => plasticState(p)%dotState (startIndex:endIndex,1:NipcMyPhase)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
!----------------------------------------------------------------------------------------------
!locally define deltaState alias
o = endIndex
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,1:NipcMyPhase)
delta%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,1:NipcMyPhase)
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,1:NipcMyPhase)
delta%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,1:NipcMyPhase)
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,1:NipcMyPhase)
delta%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,1:NipcMyPhase)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
offset_slip = plasticState(p)%nSlip
plasticState(p)%slipRate => &
plasticState(p)%dotState(offset_slip+1:offset_slip+plasticState(p)%nSlip,1:NipcMyPhase)
plasticState(p)%accumulatedSlip => &
plasticState(p)%state(offset_slip+1:offset_slip+plasticState(p)%nSlip,1:NipcMyPhase)
end associate
end do
end subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
! sanity checks
! if (any(plastic_kinehardening_Nslip (1:nSlipFamilies,instance) > 0_pInt &
! .and. param(instance)%crss0 (1:nSlipFamilies) < 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
! if (any(plastic_kinehardening_Nslip (1:nSlipFamilies,instance) > 0_pInt &
! .and. param(instance)%tau1 (1:nSlipFamilies) <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
! if (any(plastic_kinehardening_Nslip (1:nSlipFamilies,instance) > 0_pInt &
! .and. param(instance)%tau1_b(1:nSlipFamilies) < 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
! if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
! if (param(instance)%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1_pInt
endIndex = prm%totalNslip
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
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_pInt
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_pInt
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
delta%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
delta%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
delta%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)
implicit none
real(pReal), dimension(3,3), intent(out) :: &
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
@ -404,24 +374,24 @@ pure subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
of
integer(pInt) :: &
j,k,l,m,n
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
associate(prm => param(instance), stt => state(instance))
Lp = 0.0_pReal
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance), stt => state(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do j = 1_pInt, prm%totalNslip
Lp = Lp + (gdot_pos(j)+gdot_neg(j))*prm%Schmid_slip(1:3,1:3,j)
do i = 1_pInt, prm%totalNslip
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid_slip(1:3,1:3,i)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtau_pos(j)*prm%Schmid_slip(k,l,j)*prm%nonSchmid_pos(m,n,j) &
+ dgdot_dtau_neg(j)*prm%Schmid_slip(k,l,j)*prm%nonSchmid_neg(m,n,j)
+ dgdot_dtau_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
@ -447,7 +417,7 @@ subroutine plastic_kinehardening_deltaState(Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
instance, &
of
@ -455,7 +425,7 @@ subroutine plastic_kinehardening_deltaState(Mp,instance,of)
gdot_pos,gdot_neg, &
sense
associate( prm => param(instance), stt => state(instance), del => deltaState(instance))
associate(prm => param(instance), stt => state(instance), del => deltaState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,of), & ! keep existing...
@ -496,12 +466,13 @@ subroutine plastic_kinehardening_dotState(Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
integer(pInt), intent(in) :: &
instance, &
of
integer(pInt) :: &
j
i
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
real(pReal) :: &
@ -514,11 +485,11 @@ subroutine plastic_kinehardening_dotState(Mp,instance,of)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
do j = 1_pInt, prm%totalNslip
dot%crss(j,of) = dot_product(prm%interaction_SlipSlip(j,:),dot%accshear(:,of)) &
* ( prm%theta1(j) &
+ (prm%theta0(j) - prm%theta1(j) + prm%theta0(j)*prm%theta1(j)*sumGamma/prm%tau1(j)) &
*exp(-sumGamma*prm%theta0(j)/prm%tau1(j)) &
do i = 1_pInt, prm%totalNslip
dot%crss(i,of) = dot_product(prm%interaction_SlipSlip(i,:),dot%accshear(:,of)) &
* ( prm%theta1(i) &
+ (prm%theta0(i) - prm%theta1(i) + prm%theta0(i)*prm%theta1(i)*sumGamma/prm%tau1(i)) &
*exp(-sumGamma*prm%theta0(i)/prm%tau1(i)) &
)
enddo
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
@ -550,19 +521,18 @@ function plastic_kinehardening_postResults(Mp,instance,of) result(postResults)
real(pReal), dimension(sum(plastic_kinehardening_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: &
o,c,j
o,c,i
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
c = 0_pInt
associate( prm => param(instance), stt => state(instance))
associate(prm => param(instance), stt => state(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
outputsLoop: do o = 1_pInt,plastic_kinehardening_Noutput(instance)
outputsLoop: do o = 1_pInt,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
postResults(c+1_pInt:c+prm%totalNslip) = stt%crss(:,of)
@ -593,13 +563,14 @@ function plastic_kinehardening_postResults(Mp,instance,of) result(postResults)
c = c + prm%totalNslip
case (resolvedstress_ID)
do j = 1_pInt, prm%totalNslip
postResults(c+j) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,j))
do i = 1_pInt, prm%totalNslip
postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i))
enddo
c = c + prm%totalNslip
end select
enddo outputsLoop
end associate
end function plastic_kinehardening_postResults
@ -607,25 +578,27 @@ end function plastic_kinehardening_postResults
!--------------------------------------------------------------------------------------------------
!> @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
!> @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)
pure subroutine kinetics(Mp,instance,of, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
use prec, only: &
dNeq0
use math, only: &
math_mul33xx33
implicit none
real(pReal), dimension(3,3), intent(in) :: &
Mp
integer(pInt), intent(in) :: &
instance, &
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer(pInt), intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNslip), intent(out) :: &
gdot_pos, &
gdot_neg
real(pReal), dimension(param(instance)%totalNslip), optional, intent(out) :: &
real(pReal), dimension(param(instance)%totalNslip), intent(out), optional :: &
dgdot_dtau_pos, &
dgdot_dtau_neg
@ -636,7 +609,7 @@ pure subroutine kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_d
integer(pInt) :: i
logical :: nonSchmidActive
associate( prm => param(instance), stt => state(instance))
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt

View File

@ -77,6 +77,7 @@ module plastic_phenopowerlaw
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type, private :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
@ -392,9 +393,9 @@ end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @details asumme that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume (
! equally (Taylor assumption). Twinning happens only in untwinned volume
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
implicit none
real(pReal), dimension(3,3), intent(out) :: &
@ -411,11 +412,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
integer(pInt) :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
dgdot_dtauslip_pos,dgdot_dtauslip_neg, &
gdot_slip_pos,gdot_slip_neg
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
@ -626,7 +627,7 @@ pure subroutine kinetics_slip(Mp,instance,of, &
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = 0.5_pReal*prm%gdot0_slip &
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