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

18-Oct-23-1

This commit is contained in:
achalhp 2023-10-18 18:07:54 +05:30
parent d54484dfa6
commit 5aff0d96b6
1 changed files with 383 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

446
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -1,3 +1,17 @@
! Copyright 2011-2022 Max-Planck-Institut für Eisenforschung GmbH
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU Affero General Public License as
! published by the Free Software Foundation, either version 3 of the
! License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Affero General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License
! along with this program. If not, see <https://www.gnu.org/licenses/>.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
@ -7,132 +21,158 @@
submodule(phase:plastic) phenopowerlaw submodule(phase:plastic) phenopowerlaw
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: & !< declaring real number variables
dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip
dot_gamma_0_tw = 1.0_pReal, & !< reference shear strain rate for twin dot_gamma_0_tw = 1.0_pReal, & !< reference shear strain rate for twin
n_sl = 1.0_pReal, & !< stress exponent for slip n_sl = 1.0_pReal, & !< stress exponent for slip
n_tw = 1.0_pReal, & !< stress exponent for twin n_tw = 1.0_pReal, & !< stress exponent for twin
f_sat_sl_tw = 1.0_pReal, & !< push-up factor for slip saturation due to twinning f_sat_sl_tw = 1.0_pReal, & !< push-up factor for slip saturation due to twinning
c_1 = 1.0_pReal, & c_1 = 1.0_pReal, & !< hardening parameters
c_2 = 1.0_pReal, & c_2 = 1.0_pReal, &
c_3 = 1.0_pReal, & c_3 = 1.0_pReal, &
c_4 = 1.0_pReal, & c_4 = 1.0_pReal, &
h_0_sl_sl = 1.0_pReal, & !< reference hardening slip - slip h_0_sl_sl = 1.0_pReal, & !< reference hardening slip - slip
h_0_tw_sl = 1.0_pReal, & !< reference hardening twin - slip h_0_tw_sl = 1.0_pReal, & !< reference hardening twin - slip
h_0_tw_tw = 1.0_pReal, & !< reference hardening twin - twin h_0_tw_tw = 1.0_pReal, & !< reference hardening twin - twin
a_sl = 1.0_pReal h_0_tw_tw_nuc = 1.0_pReal, & !< Achal twin nucleation
real(pReal), allocatable, dimension(:) :: & h_0_tw_tw_grt = 1.0_pReal, & !< Achal twin growth
a_sl = 1.0_pReal, & !< non-Schmid Coefficient
chkstep_nucl = 1.0_pReal, & !< Achal Monte Carlo sampling frequency (for twin nucleation)
chkstep_grow = 1.0_pReal, & !< Achal Monte Carlo sampling frequency (for twin growth)
chkgrowth_twin = 1.0_pReal, & !< Achal flag for twin growth
prefdecay_slip = 1.0_pReal, & !< Achal pre-factor for the flag_slip decay
twin_inclusion = 1.0_pReal !< Achal flag_slip to introduce specific twin vol fraction to particular phase
real(pReal), allocatable, dimension(:) :: & !< 1D array
xi_inf_sl, & !< maximum critical shear stress for slip xi_inf_sl, & !< maximum critical shear stress for slip
h_int, & !< per family hardening activity (optional) h_int, & !< per family hardening activity (optional)
gamma_char !< characteristic shear for twins gamma_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: & !< 2D array
h_sl_sl, & !< slip resistance from slip activity h_sl_sl, & !< slip resistance from slip activity
h_sl_tw, & !< slip resistance from twin activity h_sl_tw, & !< slip resistance from twin activity
h_tw_sl, & !< twin resistance from slip activity h_tw_sl, & !< twin resistance from slip activity
h_tw_tw !< twin resistance from twin activity !h_tw_tw!, & !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: & h_tw_tw_grow, &
P_sl, & h_tw_tw_nucl
P_tw, & real(pReal), allocatable, dimension(:,:,:) :: & !< 3D array
P_nS_pos, & P_sl, & !< Schmid slip
P_nS_neg P_tw, & !< Schmid twin
P_nS_pos, & !< Non Schmid +ve
P_nS_neg, & !< Non Schmid -ve
Corrs_Matr !< Achal Correspondence Matrix
integer :: & integer :: &
sum_N_sl, & !< total number of active slip system sum_N_sl, & !< total number of active slip system
sum_N_tw !< total number of active twin systems sum_N_tw !< total number of active twin systems
logical :: & logical :: & !< flag if non-Schmid slip is Active
nonSchmidActive = .false. nonSchmidActive = .false.
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: & !< Storing output messages
output output
character(len=:), allocatable, dimension(:) :: & character(len=:), allocatable, dimension(:) :: & !< Labels of active slip and twin systems
systems_sl, & systems_sl, &
systems_tw systems_tw
end type tParameters end type tParameters
type :: tIndexDotState type :: tIndexDotState !< Index to access variables in 2D tPhenopowerlawState
integer, dimension(2) :: & integer, dimension(2) :: &
xi_sl, & xi_sl, &
xi_tw, & !xi_tw, &
gamma_sl, & gamma_sl, &
gamma_tw gamma_tw, &
xi_tw_nucl, &
xi_tw_grow
end type tIndexDotState end type tIndexDotState
type :: tPhenopowerlawState type :: tPhenopowerlawState !< state variables at material points
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
xi_sl, & xi_sl, & !< critical shear stress for slip
xi_tw, & !xi_tw, & !< critical shear stress for twin
gamma_sl, & gamma_sl, & !< shear strain due to slip
gamma_tw gamma_tw, & !< shear strain due to twin
xi_tw_nucl, &
xi_tw_grow, &
f_tw_nucl, &
f_tw_grow, &
fmc_tw_nucl, &
fmc_tw_grwo
!real(pReal), pointer, dimension(:) :: &
! variant_twin, & !< flag used to assign twin variant
! frozen !< 0 to 1
end type tPhenopowerlawState end type tPhenopowerlawState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters, dot state index, and state ! containers for parameters, dot state index, and state
! containers are arrays used to store parameters, state indices and state variables for multiple material points
type(tParameters), allocatable, dimension(:) :: param type(tParameters), allocatable, dimension(:) :: param
type(tIndexDotState), allocatable, dimension(:) :: indexDotState type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tPhenopowerlawState), allocatable, dimension(:) :: state type(tPhenopowerlawState), allocatable, dimension(:) :: state
contains contains
!< "contains" means the above types are used in below functions.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Perform module initialization. !> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_phenopowerlaw_init() result(myPlasticity) module function plastic_phenopowerlaw_init() result(myPlasticity) !< function "plastic_phenopowerlaw_init" returns myPlasticity
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity !< logical array: phenopowerlaw is active or not
integer :: & integer :: & !< integer variables to iterate
ph, i, & ph, i, &
Nmembers, & Nmembers, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: & !< active number of slip and twin systems
N_sl, N_tw N_sl, N_tw
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
xi_0_sl, & !< initial critical shear stress for slip xi_0_sl, & !< initial critical shear stress for slip
xi_0_tw, & !< initial critical shear stress for twin xi_0_tw, & !< initial critical shear stress for twin
a !< non-Schmid coefficients a !< non-Schmid coefficients
character(len=pStringLen) :: & character(len=pStringLen) :: & !< to store error messages
extmsg = '' extmsg = ''
class(tNode), pointer :: & class(tNode), pointer :: & !< pointers variables, tNode objects
phases, & phases, &
phase, & phase, &
mech, & mech, & !< mechanical model info
pl pl !< plastic model info
myPlasticity = plastic_active('phenopowerlaw') myPlasticity = plastic_active('phenopowerlaw') !< function defined in phase_mechanical_plastic
if(count(myPlasticity) == 0) return if(count(myPlasticity) == 0) return !< if phenopowerlaw is active or not
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) !< print for which phases phenopowerlaw is active
phases => config_material%get('phase') phases => config_material%get('phase') !< get "phase" objects from config_material object
allocate(param(phases%length)) allocate(param(phases%length))
allocate(indexDotState(phases%length)) allocate(indexDotState(phases%length))
allocate(state(phases%length)) allocate(state(phases%length))
do ph = 1, phases%length do ph = 1, phases%length !< looping for each phase
if (.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), stt => state(ph), & associate(prm => param(ph), stt => state(ph), & !< parameters, state var, dot state indices of current phase
idx_dot => indexDotState(ph)) idx_dot => indexDotState(ph))
phase => phases%get(ph) phase => phases%get(ph) !< getting objects
mech => phase%get('mechanical') mech => phase%get('mechanical')
pl => mech%get('plastic') pl => mech%get('plastic')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) !< Read and Store no. of slip systems
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl)) !< total no. of slip systems
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then !< if +ve, active slip system exists
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) !<
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then !< cI: BCC
a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray) a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray) !< a: non-Schmid coefficients for projections
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) !< function returns non-schmid projections
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
@ -177,6 +217,14 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'),phase_lattice(ph)) prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'),phase_lattice(ph))
prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
!prm%h_tw_tw_nucl =
!prm%h_tw_tw_grow =
!prm%chkstep_nucl =
!prm%chkstep_grow =
!prm%chkgrowth_twin =
!prm%twin_inclusion =
!prm%prefdecay_slip =
xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw)) xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw))
@ -200,6 +248,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
xi_0_tw = emptyRealArray xi_0_tw = emptyRealArray
allocate(prm%gamma_char,source=emptyRealArray) allocate(prm%gamma_char,source=emptyRealArray)
allocate(prm%h_tw_tw(0,0)) allocate(prm%h_tw_tw(0,0))
!allocate(prm%h_tw_tw_nucl(0,0))
!allocate(prm%h_tw_tw_grow(0,0))
end if twinActive end if twinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -229,8 +279,16 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw !&
sizeState = sizeDotState + size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw & ! Why not size(['xi_tw_nucl','gamma_tw'])?
!+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw &
!+ size(['variant_twin','frozen']) * prm%sum_N_tw &
sizeState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw !&
+ size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw & ! Why not size(['xi_tw_nucl','gamma_tw'])?
!+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw &
!+ size(['fmc_tw_nucl','fmc_tw_grow']) * prm%sum_N_tw &
!+ size(['variant_twin','frozen']) * prm%sum_N_tw &
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
@ -245,12 +303,12 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1 !startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw !endIndex = endIndex + prm%sum_N_tw
idx_dot%xi_tw = [startIndex,endIndex] !idx_dot%xi_tw = [startIndex,endIndex]
stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:) !stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw = spread(xi_0_tw, 2, Nmembers) !stt%xi_tw = spread(xi_0_tw, 2, Nmembers)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
@ -265,6 +323,50 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:) stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
idx_dot%xi_tw_nucl = [startIndex,endIndex]
stt%xi_tw_nucl => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw_nucl = spread(xi_0_tw, 2, Nmembers)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
idx_dot%xi_tw_grow = [startIndex,endIndex]
stt%xi_tw_grow => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw_grow = spread(xi_0_tw, 2, Nmembers)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
!startIndex = endIndex + 1
!endIndex = endIndex + prm%sum_N_tw
!idx_dot%f_twin_nucl = [startIndex,endIndex]
!stt%f_twin_nucl => plasticState(ph)%state(startIndex:endIndex,:)
!plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
!startIndex = endIndex + 1
!endIndex = endIndex + prm%sum_N_tw
!idx_dot%f_twin_grow = [startIndex,endIndex]
!stt%f_twin_grow => plasticState(ph)%state(startIndex:endIndex,:)
!plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
!startIndex = endIndex + 1
!endIndex = endIndex + prm%sum_N_tw
!idx_dot%fmc_twin_nucl = [startIndex,endIndex]
!stt%fmc_twin_nucl => plasticState(ph)%state(startIndex:endIndex,:)
!plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
!startIndex = endIndex + 1
!endIndex = endIndex + prm%sum_N_tw
!idx_dot%fmc_twin_grow = [startIndex,endIndex]
!stt%fmc_twin_grow => plasticState(ph)%state(startIndex:endIndex,:)
!plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
!startIndex = endIndex + 1
!endIndex = endIndex + prm%sum_N_tw
!idx_dot%frozen = [startIndex,endIndex]
!stt%frozen => plasticState(ph)%state(startIndex:endIndex,:)
!plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
end associate end associate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -299,9 +401,9 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl_pos,dot_gamma_sl_neg, & dot_gamma_sl_pos,dot_gamma_sl_neg, &
ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_tw,ddot_gamma_dtau_tw dot_gamma_tw,ddot_gamma_dtau_tw
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
@ -332,7 +434,7 @@ end subroutine phenopowerlaw_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure. !> @brief Calculate the rate of change of microstructure.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) !< Rate of change of state variables
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -348,13 +450,23 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl_pos,dot_gamma_sl_neg, & dot_gamma_sl_pos,dot_gamma_sl_neg, &
right_SlipSlip right_SlipSlip
!real(pReal), dimension(param(ph)%sum_N_tw) :: &
!fdot_twin_nucl, fdot_twin_grow
associate(prm => param(ph), stt => state(ph), & associate(prm => param(ph), stt => state(ph), &
dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), & dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), & !dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), & dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2))) !dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2))!, &
dot_gamma_tw_nucl => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)), &
dot_gamma_tw_grow => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)), &
dot_xi_tw_nucl => dotState(indexDotState(ph)%xi_tw_nucl(1):indexDotState(ph)%xi_tw_nucl(2)), &
dot_xi_tw_grow => dotState(indexDotState(ph)%xi_tw_grow(1):indexDotState(ph)%xi_tw_grow(2))
!sumGamma
!sumF_nucl
!sumF_grow
call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg) call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
@ -368,15 +480,219 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) & * matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot_gamma_tw) + matmul(prm%h_sl_tw,dot_gamma_tw)
!< Rate of change of critical shear stress
!dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
! * matmul(prm%h_tw_sl,dot_gamma_sl) &
! + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw)
dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & dot_xi_tw_nucl = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot_gamma_sl) & * matmul(prm%h_tw_sl,dot_gamma_sl) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw) + prm%h_0_tw_tw_nucl * sumF_nucl**prm%c_4 * matmul(prm%h_tw_tw_nucl,dot_gamma_tw)
dot_xi_tw_grow = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot_gamma_sl) &
+ prm%h_0_tw_tw_grow * sumF_grow**prm%c_4 * matmul(prm%h_tw_tw_grow,dot_gamma_tw)
end associate end associate
end function phenopowerlaw_dotState end function phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates instantaneous incremental change of kinematics and associated jump state
!--------------------------------------------------------------------------------------------------
! module subroutine plastic_kinematic_deltaFp(twin)
! use prec, only: &
! dNeq, &
! dEq0
! #ifdef DEBUG
! use debug, only: &
! debug_level, &
! debug_constitutive,&
! debug_levelExtensive, &
! debug_levelSelective
! #endif
! use mesh, only: &
! mesh_element, &
! mesh_ipNeighborhood, &
! mesh_ipCoordinates, &
! mesh_ipVolume, &
! mesh_ipAreaNormal, &
! mesh_ipArea, &
! FE_NipNeighbors, &
! mesh_maxNipNeighbors, &
! FE_geomtype, &
! FE_celltype
! use lattice
! use math, only: &
! math_I3
! use material, only: &
! phaseAt, phasememberAt, &
! phase_plasticityInstance
! implicit none
! integer(pInt) :: &
! ph, of, instance, &
! neighbor_el, & !< element number of neighboring material point
! neighbor_ip, & !< integration point of neighboring material point
! np, & !< neighbor phase
! no, n !< nieghbor offset and index for loop at neighbor
! integer(pInt), intent(in) :: &
! el, & !< element index
! ip, & !< integration point index
! ipc !< grain index
! real(pReal), dimension(3,3), intent(out) :: &
! deltaFp
! logical , intent(out) :: &
! twinJump
! ! real(pReal), dimension(3,3,param(instance)%totalNslip) :: &
! ! CorrespondanceMatrix
! integer(pInt), dimension(52) :: &
! twin_el_incl
! real(pReal), dimension(6) :: &
! neighbor_stt
! real(pReal) :: &
! random, random1
! integer(pInt) :: &
! i,j,var_growth,var_nucl
! var_growth = 0_pInt
! var_nucl = 0_pInt
! ph = phaseAt(ipc, ip, el)
! of = phasememberAt(ipc, ip, el)
! instance = phase_plasticityInstance(ph)
! associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
! twinJump = .false.
! deltaFp = math_I3
! ! for eshelby circular inclusion
! twin_el_incl = (/ 10913,10914,10915,10916,10917,10918,10919,10920,10921,10922,10923,10924,10925, &
! 10993,10994,10995,10996,10997,10998,10999,11000,11001,11002,11074,11075,11076, &
! 11077,11078,11079,10751,10752,10753,10754,10755,10756,10757,10758,10759,10760, &
! 10670,10671,10672,10673,10674,10675,10676,10677,10678,10679,10680,10681,10682 /)
! ! TwinLooptest: do i=1_pInt, prm%totalNtwin
! ! write(6,*)'CorrespondenceMatrix for system',i, prm%CorrespondanceMatrix(:,:,i)
! ! enddo TwinLooptest
! !Saving the neighbor information in an array
! NeighborLoop1: do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! only 4 neighbors for quasi 2D (1 element in z direction)
! neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
! neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
! np = phaseAt(1,neighbor_ip,neighbor_el)
! no = phasememberAt(1,neighbor_ip,neighbor_el)
! neighbor_stt(n) = state(phase_plasticityInstance(np))%variant_twin(no)
! enddo NeighborLoop1
! !checking if any of my neighbor is twinned if yes recognize the variant and exit
! ! NeighborLoop2: do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! only 4 neighbors for quasi 2D (1 element in z direction)
! ! neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
! ! neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
! ! np = phaseAt(1,neighbor_ip,neighbor_el)
! ! ! if(of == 1) write(6,*)'phaseAt neighbor_ip of neighbor_el', np
! ! no = phasememberAt(1,neighbor_ip,neighbor_el)
! ! ! if(of == 1) write(6,*)'phasememberAt at neighbor_ip of neighbor_el', no
! ! if (state(phase_plasticityInstance(np))%variant_twin(no) > 0_pInt) then
! ! var_growth = state(phase_plasticityInstance(np))%variant_twin(no)
! !
! ! exit NeighborLoop2
! ! endif
! ! enddo NeighborLoop2
! call RANDOM_NUMBER(random)
! call RANDOM_NUMBER(random1)
! ! Sampling: if (var_growth > 0_pInt) then
! ! ! write(6,*)'I am sampling for growth with variant',var_growth
! ! Ability_Growth: if (stt%f_twin_grow(var_growth,of) > stt%fmc_twin_grow(var_growth,of) &
! ! + prm%checkstep_grow) then
! ! stt%fmc_twin_grow(var_growth,of) = stt%fmc_twin_grow(var_growth,of) &
! ! + prm%checkstep_grow
! ! Success_Growth: if (random <= stt%f_twin_grow(var_growth,of) .or. ALL(neighbor_stt > 0_pReal)) then
! ! write(6,*)'growth sampling is successful for elem',el
! ! twinJump = .true.
! ! deltaFp = prm%CorrespondanceMatrix(:,:,var_growth)
! ! dlt%f_twin_grow(:,of) = 0.0_pReal - stt%f_twin_grow(:,of)
! ! dlt%f_twin_nucl(:,of) = 0.0_pReal - stt%f_twin_nucl(:,of)
! ! dlt%fmc_twin_grow(:,of) = 0.0_pReal - stt%fmc_twin_grow(:,of)
! ! dlt%fmc_twin_nucl(:,of) = 0.0_pReal - stt%fmc_twin_nucl(:,of)
! ! dlt%frozen(of) = 1.0_pReal - stt%frozen(of)
! ! dlt%variant_twin(of) = var_growth - stt%variant_twin(of)
! ! endif Success_Growth
! ! endif Ability_Growth
! ! elseif (var_growth == 0_pInt .and. prm%checkgrowth_twin > 0_pReal ) then
! if (var_growth == 0_pInt .and. prm%checkgrowth_twin > 0_pReal ) then
! var_nucl = maxloc(stt%f_twin_nucl(:,of), dim=1)
! ! write(6,*)'I am sampling for nucleation with variant',var_nucl,stt%f_twin_nucl(var_nucl,of)
! Ability_Nucleation: if (stt%f_twin_nucl(var_nucl,of) > stt%fmc_twin_nucl(var_nucl,of) &
! + prm%checkstep_nucl) then
! stt%fmc_twin_nucl(var_nucl,of) = stt%fmc_twin_nucl(var_nucl,of) &
! + prm%checkstep_nucl
! Success_Nucleation: if (random <= stt%f_twin_nucl(var_nucl,of) &
! .and. random1 <= 0.20) then
! write(6,*)'nucleation sampling is successful for elem',el
! twinJump = .true.
! deltaFp = prm%CorrespondanceMatrix(:,:,var_nucl)
! dlt%f_twin_nucl(:,of) = 0.0_pReal - stt%f_twin_nucl(:,of)
! dlt%f_twin_grow(:,of) = 0.0_pReal - stt%f_twin_grow(:,of)
! dlt%fmc_twin_nucl(:,of) = 0.0_pReal - stt%fmc_twin_nucl(:,of)
! dlt%fmc_twin_grow(:,of) = 0.0_pReal - stt%fmc_twin_grow(:,of)
! dlt%frozen(of) = 1.0_pReal - stt%frozen(of)
! dlt%variant_twin(of) = var_nucl - stt%variant_twin(of)
! endif Success_Nucleation
! endif Ability_Nucleation
! endif
! ! endif Sampling
! end associate
! end subroutine plastic_kinematic_deltaFp
!--------------------------------------------------------------------------------------------------
!> @brief calculates (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
!subroutine plastic_phenopowerlaw_deltaState(instance,of)
! use prec, only: &
! dNeq, &
! dEq0
! #ifdef DEBUG
! use debug, only: &
! debug_level, &
! debug_constitutive,&
! debug_levelExtensive, &
! debug_levelSelective
! #endif
! implicit none
! integer(pInt), intent(in) :: &
! instance, &
! of
!
! associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
! #ifdef DEBUG
! if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
! .and. (of == prm%of_debug &
! .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
! write(6,'(a)') '======= phenopowerlaw delta state ======='
! ! write(6,*) sense,state(instance)%sense(:,of)
! endif
! #endif
! !--------------------------------------------------------------------------------------------------
! dlt%f_twin_nucl(:,of) = 0.0_pReal
! dlt%f_twin_grow(:,of) = 0.0_pReal
! dlt%fmc_twin_nucl(:,of) = 0.0_pReal
! dlt%fmc_twin_grow(:,of) = 0.0_pReal
! dlt%frozen(of) = 0.0_pReal
! dlt%variant_twin(of) = 0.0_pInt
!
! end associate
!end subroutine plastic_phenopowerlaw_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file. !> @brief Write results to HDF5 output file.
@ -494,9 +810,10 @@ end subroutine kinetics_sl
!> @details Derivatives are calculated only optionally. !> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
! Mp: Mandel Stress, ph: , en: , dot_gamma_tw: , ddot_gamma_dtau_tw:
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_tw(Mp,ph,en,& pure subroutine kinetics_tw(Mp,ph,en,&
dot_gamma_tw,ddot_gamma_dtau_tw) dot_gamma_tw,ddot_gamma_dtau_tw) !fdot_twin_nucl, fdot_twin_grow)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -505,7 +822,7 @@ pure subroutine kinetics_tw(Mp,ph,en,&
en en
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
dot_gamma_tw dot_gamma_tw !fdot_twin_nucl, fdot_twin_grow
real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
ddot_gamma_dtau_tw ddot_gamma_dtau_tw
@ -518,9 +835,11 @@ pure subroutine kinetics_tw(Mp,ph,en,&
tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)]
where(tau_tw > 0.0_pReal) where(tau_tw > 0.0_pReal) !and stt%frozen(en) < 0.9_pReal)
dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction
* prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw
!fdot_twin_nucl =
!fdot_twin_nucl =
else where else where
dot_gamma_tw = 0.0_pReal dot_gamma_tw = 0.0_pReal
end where end where
@ -538,3 +857,4 @@ pure subroutine kinetics_tw(Mp,ph,en,&
end subroutine kinetics_tw end subroutine kinetics_tw
end submodule phenopowerlaw end submodule phenopowerlaw