DAMASK_EICMD/src/phase_mechanical_plastic_no...

1778 lines
97 KiB
Fortran
Raw Normal View History

!--------------------------------------------------------------------------------------------------
!> @author Christoph Kords, 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
!> @brief material subroutine for plasticity including dislocation flux
!--------------------------------------------------------------------------------------------------
submodule(phase:plastic) nonlocal
use geometry_plastic_nonlocal, only: &
2019-06-07 13:50:56 +05:30
nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
2019-06-06 14:38:58 +05:30
IPvolume => geometry_plastic_nonlocal_IPvolume0, &
IParea => geometry_plastic_nonlocal_IParea0, &
2020-12-20 21:41:43 +05:30
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, &
geometry_plastic_nonlocal_disable
type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0
end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom
2019-03-17 21:32:08 +05:30
! storage order of dislocation types
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
sgl = [1,2,3,4,5,6,7,8] !< signed (single)
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
edg = [1,2,5,6,9], & !< edge
scr = [3,4,7,8,10] !< screw
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
mob = [1,2,3,4], & !< mobile
imm = [5,6,7,8] !< immobile (blocked)
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
dip = [9,10], & !< dipole
imm_edg = imm(1:2), & !< immobile edge
imm_scr = imm(3:4) !< immobile screw
integer, parameter :: &
mob_edg_pos = 1, & !< mobile edge positive
mob_edg_neg = 2, & !< mobile edge negative
mob_scr_pos = 3, & !< mobile screw positive
mob_scr_neg = 4 !< mobile screw positive
2020-03-16 20:06:34 +05:30
! BEGIN DEPRECATED
2019-12-04 23:30:56 +05:30
integer, dimension(:,:,:), allocatable :: &
2019-03-17 21:32:08 +05:30
iRhoU, & !< state indices for unblocked density
2020-06-26 15:14:17 +05:30
iV, & !< state indices for dislocation velocities
2019-03-17 21:32:08 +05:30
iD !< state indices for stable dipole height
!END DEPRECATED
2019-12-04 23:30:56 +05:30
real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
2021-04-25 11:36:52 +05:30
compatibility !< slip system compatibility between en and my neighbors
2020-03-17 12:05:25 +05:30
type :: tInitialParameters !< container type for internal constitutive parameters
real(pReal) :: &
2020-09-22 21:41:09 +05:30
sigma_rho_u, & !< standard deviation of scatter in initial dislocation density
random_rho_u, &
random_rho_u_binning
2020-03-17 12:05:25 +05:30
real(pReal), dimension(:), allocatable :: &
2020-09-22 21:41:09 +05:30
rho_u_ed_pos_0, & !< initial edge_pos dislocation density
rho_u_ed_neg_0, & !< initial edge_neg dislocation density
rho_u_sc_pos_0, & !< initial screw_pos dislocation density
rho_u_sc_neg_0, & !< initial screw_neg dislocation density
rho_d_ed_0, & !< initial edge dipole dislocation density
rho_d_sc_0 !< initial screw dipole dislocation density
2020-05-31 01:09:35 +05:30
integer, dimension(:), allocatable :: &
2020-03-17 12:05:25 +05:30
N_sl
end type tInitialParameters
2019-12-04 23:30:56 +05:30
type :: tParameters !< container type for internal constitutive parameters
2019-03-17 21:32:08 +05:30
real(pReal) :: &
2020-09-22 21:41:09 +05:30
V_at, & !< atomic volume
D_0, & !< prefactor for self-diffusion coefficient
Q_cl, & !< activation enthalpy for diffusion
2020-03-14 23:41:26 +05:30
atol_rho, & !< absolute tolerance for dislocation density in state integration
2020-09-22 21:41:09 +05:30
rho_significant, & !< density considered significant
2020-09-23 04:25:19 +05:30
rho_min, & !< number of dislocations considered significant
w, & !< width of a doubkle kink in multiples of the Burgers vector length b
2020-09-22 21:41:09 +05:30
Q_sol, & !< activation energy for solid solution in J
2020-09-23 05:36:03 +05:30
f_sol, & !< solid solution obstacle size in multiples of the Burgers vector length
2020-09-22 21:41:09 +05:30
c_sol, & !< concentration of solid solution in atomic parts
2019-03-17 21:32:08 +05:30
p, & !< parameter for kinetic law (Kocks,Argon,Ashby)
q, & !< parameter for kinetic law (Kocks,Argon,Ashby)
2021-06-24 20:01:07 +05:30
B, & !< drag coefficient in Pa s
2020-09-22 21:41:09 +05:30
nu_a, & !< attack frequency in Hz
chi_surface, & !< transmissivity at free surface
chi_GB, & !< transmissivity at grain boundary (identified by different texture)
2021-06-24 22:44:19 +05:30
C_CFL, & !< safety factor for CFL flux condition
2020-09-22 21:41:09 +05:30
f_ed_mult, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
f_F, &
f_ed, &
2019-03-17 21:32:08 +05:30
mu, &
nu
2020-03-17 12:05:25 +05:30
real(pReal), dimension(:), allocatable :: &
2020-09-22 21:41:09 +05:30
d_ed, & !< minimum stable edge dipole height
d_sc, & !< minimum stable screw dipole height
2020-09-23 04:25:19 +05:30
tau_Peierls_ed, &
tau_Peierls_sc, &
2020-09-22 21:41:09 +05:30
i_sl, & !< mean free path prefactor for each
2020-09-23 05:36:03 +05:30
b_sl !< absolute length of Burgers vector [m]
2020-03-17 12:05:25 +05:30
real(pReal), dimension(:,:), allocatable :: &
2019-03-17 21:32:08 +05:30
slip_normal, &
slip_direction, &
slip_transverse, &
minDipoleHeight, & ! edge and screw
peierlsstress, & ! edge and screw
2020-09-22 21:41:09 +05:30
h_sl_sl ,& !< coefficients for slip-slip interaction
2019-03-17 21:32:08 +05:30
forestProjection_Edge, & !< matrix of forest projections of edge dislocations
forestProjection_Screw !< matrix of forest projections of screw dislocations
2019-12-04 23:30:56 +05:30
real(pReal), dimension(:,:,:), allocatable :: &
2021-07-22 03:35:49 +05:30
P_sl, & !< Schmid contribution
2021-07-21 19:59:57 +05:30
P_nS_pos, &
P_nS_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
2019-03-17 21:32:08 +05:30
integer :: &
sum_N_sl = 0
2020-03-17 12:05:25 +05:30
integer, dimension(:), allocatable :: &
2019-03-17 21:32:08 +05:30
colinearSystem !< colinear system to the active slip system (only valid for fcc!)
2020-03-17 12:05:25 +05:30
character(len=pStringLen), dimension(:), allocatable :: &
2020-02-14 13:56:26 +05:30
output
2019-12-04 23:30:56 +05:30
logical :: &
2020-03-17 12:05:25 +05:30
shortRangeStressCorrection, & !< use of short range stress correction by excess density gradient term
nonSchmidActive = .false.
character(len=:), allocatable, dimension(:) :: &
systems_sl
2019-03-17 21:32:08 +05:30
end type tParameters
2021-07-27 02:26:40 +05:30
type :: tNonlocalDependentState
2019-03-17 21:32:08 +05:30
real(pReal), allocatable, dimension(:,:) :: &
tau_pass, &
tau_Back
2021-07-27 02:26:40 +05:30
end type tNonlocalDependentState
2019-12-04 23:30:56 +05:30
type :: tNonlocalState
2019-03-17 21:32:08 +05:30
real(pReal), pointer, dimension(:,:) :: &
rho, & ! < all dislocations
rhoSgl, &
rhoSglMobile, & ! iRhoU
2019-03-17 21:32:08 +05:30
rho_sgl_mob_edg_pos, &
rho_sgl_mob_edg_neg, &
rho_sgl_mob_scr_pos, &
rho_sgl_mob_scr_neg, &
2020-03-16 20:06:34 +05:30
rhoSglImmobile, &
2019-03-17 21:32:08 +05:30
rho_sgl_imm_edg_pos, &
rho_sgl_imm_edg_neg, &
rho_sgl_imm_scr_pos, &
rho_sgl_imm_scr_neg, &
2020-03-16 20:06:34 +05:30
rhoDip, &
2019-03-17 21:32:08 +05:30
rho_dip_edg, &
rho_dip_scr, &
rho_forest, &
2019-12-01 14:05:44 +05:30
gamma, &
2019-12-01 13:25:24 +05:30
v, &
v_edg_pos, &
v_edg_neg, &
v_scr_pos, &
v_scr_neg
2019-03-17 21:32:08 +05:30
end type tNonlocalState
2020-07-03 20:15:11 +05:30
2020-07-02 02:21:21 +05:30
type(tNonlocalState), allocatable, dimension(:) :: &
2019-03-17 21:32:08 +05:30
deltaState, &
dotState, &
state, &
state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
2021-07-27 02:26:40 +05:30
type(tNonlocalDependentState), dimension(:), allocatable :: dependentState
2019-02-18 14:58:08 +05:30
contains
!--------------------------------------------------------------------------------------------------
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
2020-08-15 19:32:10 +05:30
module function plastic_nonlocal_init() result(myPlasticity)
2019-02-21 04:54:35 +05:30
2020-08-15 19:32:10 +05:30
logical, dimension(:), allocatable :: myPlasticity
2019-03-17 18:05:41 +05:30
integer :: &
Ninstances, &
2021-02-16 20:36:09 +05:30
ph, &
2021-03-05 01:46:36 +05:30
Nmembers, &
2020-03-17 12:05:25 +05:30
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
2019-02-22 13:51:04 +05:30
s1, s2, &
2020-03-17 12:05:25 +05:30
s, t, l
real(pReal), dimension(:), allocatable :: &
a
2019-12-21 16:58:24 +05:30
character(len=pStringLen) :: &
2020-03-16 20:06:34 +05:30
extmsg = ''
2020-03-17 12:05:25 +05:30
type(tInitialParameters) :: &
ini
2020-08-15 19:32:10 +05:30
class(tNode), pointer :: &
phases, &
phase, &
2020-11-03 03:16:46 +05:30
mech, &
2020-08-15 19:32:10 +05:30
pl
2019-03-09 15:32:12 +05:30
2020-08-15 19:32:10 +05:30
myPlasticity = plastic_active('nonlocal')
Ninstances = count(myPlasticity)
if (Ninstances == 0) then
2020-08-15 19:32:10 +05:30
call geometry_plastic_nonlocal_disable
return
end if
2020-12-20 21:41:43 +05:30
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:nonlocal init -+>>>'
print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
2021-02-13 17:37:35 +05:30
print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333348, 2014'
print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
2020-08-15 19:32:10 +05:30
print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
print'( 1x,a)', 'http://publications.rwth-aachen.de/record/229993'
2020-07-02 00:52:05 +05:30
phases => config_material%get('phase')
2021-06-24 18:38:31 +05:30
allocate(geom(phases%length))
allocate(param(phases%length))
allocate(state(phases%length))
allocate(state0(phases%length))
allocate(dotState(phases%length))
allocate(deltaState(phases%length))
2021-07-27 02:26:40 +05:30
allocate(dependentState(phases%length))
2013-01-22 05:20:28 +05:30
2021-02-16 20:36:09 +05:30
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
2021-02-16 20:36:09 +05:30
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), &
2021-07-27 02:26:40 +05:30
st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
2021-02-16 20:36:09 +05:30
phase => phases%get(ph)
2021-03-25 23:52:59 +05:30
mech => phase%get('mechanical')
pl => mech%get('plastic')
2020-08-15 19:32:10 +05:30
2021-06-25 16:56:08 +05:30
plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.)
2020-08-15 19:32:10 +05:30
#if defined (__GFORTRAN__)
2021-03-11 22:30:07 +05:30
prm%output = output_as1dString(pl)
2020-08-15 19:32:10 +05:30
#else
2021-03-11 22:30:07 +05:30
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
2020-08-15 19:32:10 +05:30
#endif
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
2021-03-11 22:30:07 +05:30
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
2020-03-17 12:05:25 +05:30
prm%sum_N_sl = sum(abs(ini%N_sl))
2020-03-16 23:36:10 +05:30
slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(ini%N_sl,phase_lattice(ph))
2021-07-22 03:35:49 +05:30
prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph))
2021-07-21 19:16:38 +05:30
if (phase_lattice(ph) == 'cI') then
2021-03-11 22:30:07 +05:30
a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
2020-03-17 12:05:25 +05:30
if(size(a) > 0) prm%nonSchmidActive = .true.
2021-07-21 19:59:57 +05:30
prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
2019-02-21 04:54:35 +05:30
else
2021-07-22 03:35:49 +05:30
prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl
end if
2019-02-21 04:54:35 +05:30
2021-07-21 19:16:38 +05:30
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), &
phase_lattice(ph))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph))
prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph))
prm%slip_direction = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_normal = lattice_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
2019-02-21 04:54:35 +05:30
! collinear systems (only for octahedral slip systems in fcc)
2020-03-16 23:36:10 +05:30
allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)
do s1 = 1, prm%sum_N_sl
do s2 = 1, prm%sum_N_sl
2019-02-21 04:54:35 +05:30
if (all(dEq0 (math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. &
2019-02-21 23:48:06 +05:30
any(dNeq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) &
2019-02-21 04:54:35 +05:30
prm%colinearSystem(s1) = s2
end do
end do
2021-03-11 22:30:07 +05:30
ini%rho_u_ed_pos_0 = pl%get_as1dFloat('rho_u_ed_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_ed_neg_0 = pl%get_as1dFloat('rho_u_ed_neg_0', requiredSize=size(ini%N_sl))
ini%rho_u_sc_pos_0 = pl%get_as1dFloat('rho_u_sc_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_sc_neg_0 = pl%get_as1dFloat('rho_u_sc_neg_0', requiredSize=size(ini%N_sl))
ini%rho_d_ed_0 = pl%get_as1dFloat('rho_d_ed_0', requiredSize=size(ini%N_sl))
ini%rho_d_sc_0 = pl%get_as1dFloat('rho_d_sc_0', requiredSize=size(ini%N_sl))
2021-03-11 22:30:07 +05:30
prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(ini%N_sl))
prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(ini%N_sl))
2020-09-22 21:41:09 +05:30
prm%i_sl = math_expand(prm%i_sl,ini%N_sl)
prm%b_sl = math_expand(prm%b_sl,ini%N_sl)
2021-03-11 22:30:07 +05:30
prm%d_ed = pl%get_as1dFloat('d_ed', requiredSize=size(ini%N_sl))
prm%d_sc = pl%get_as1dFloat('d_sc', requiredSize=size(ini%N_sl))
2020-09-22 21:41:09 +05:30
prm%d_ed = math_expand(prm%d_ed,ini%N_sl)
prm%d_sc = math_expand(prm%d_sc,ini%N_sl)
2020-03-16 23:36:10 +05:30
allocate(prm%minDipoleHeight(prm%sum_N_sl,2))
2020-09-22 21:41:09 +05:30
prm%minDipoleHeight(:,1) = prm%d_ed
prm%minDipoleHeight(:,2) = prm%d_sc
2021-03-11 22:30:07 +05:30
prm%tau_Peierls_ed = pl%get_as1dFloat('tau_Peierls_ed', requiredSize=size(ini%N_sl))
prm%tau_Peierls_sc = pl%get_as1dFloat('tau_Peierls_sc', requiredSize=size(ini%N_sl))
2020-09-23 04:25:19 +05:30
prm%tau_Peierls_ed = math_expand(prm%tau_Peierls_ed,ini%N_sl)
prm%tau_Peierls_sc = math_expand(prm%tau_Peierls_sc,ini%N_sl)
2020-03-16 23:36:10 +05:30
allocate(prm%peierlsstress(prm%sum_N_sl,2))
2020-09-23 04:25:19 +05:30
prm%peierlsstress(:,1) = prm%tau_Peierls_ed
prm%peierlsstress(:,2) = prm%tau_Peierls_sc
2020-09-22 21:41:09 +05:30
prm%rho_significant = pl%get_asFloat('rho_significant')
2020-09-23 04:25:19 +05:30
prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal)
2021-06-24 22:44:19 +05:30
prm%C_CFL = pl%get_asFloat('C_CFL',defaultVal=2.0_pReal)
2020-09-22 21:41:09 +05:30
prm%V_at = pl%get_asFloat('V_at')
2020-09-22 22:56:39 +05:30
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
2020-09-22 21:41:09 +05:30
prm%f_F = pl%get_asFloat('f_F')
prm%f_ed = pl%get_asFloat('f_ed')
2020-09-22 21:41:09 +05:30
prm%w = pl%get_asFloat('w')
prm%Q_sol = pl%get_asFloat('Q_sol')
prm%f_sol = pl%get_asFloat('f_sol')
prm%c_sol = pl%get_asFloat('c_sol')
prm%p = pl%get_asFloat('p_sl')
prm%q = pl%get_asFloat('q_sl')
2021-06-24 20:01:07 +05:30
prm%B = pl%get_asFloat('B')
2020-09-22 21:41:09 +05:30
prm%nu_a = pl%get_asFloat('nu_a')
2019-02-21 04:54:35 +05:30
2019-02-21 23:48:06 +05:30
! ToDo: discuss logic
2020-09-22 21:41:09 +05:30
ini%sigma_rho_u = pl%get_asFloat('sigma_rho_u')
ini%random_rho_u = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal)
2020-08-15 19:32:10 +05:30
if (pl%contains('random_rho_u')) &
2020-09-22 21:41:09 +05:30
ini%random_rho_u_binning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default?
2019-02-21 23:48:06 +05:30
! if (rhoSglRandom(instance) < 0.0_pReal) &
! if (rhoSglRandomBinning(instance) <= 0.0_pReal) &
2020-09-22 21:41:09 +05:30
prm%chi_surface = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal)
prm%chi_GB = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal)
prm%f_ed_mult = pl%get_asFloat('f_ed_mult')
2020-08-15 19:32:10 +05:30
prm%shortRangeStressCorrection = pl%get_asBool('short_range_stress_correction', defaultVal = .false.)
2020-04-01 14:39:30 +05:30
2019-02-21 04:54:35 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
2020-09-22 22:56:39 +05:30
if (any(prm%b_sl < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(ini%rho_u_ed_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0'
if (any(ini%rho_u_ed_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0'
if (any(ini%rho_u_sc_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0'
if (any(ini%rho_u_sc_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0'
if (any(ini%rho_d_ed_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0'
if (any(ini%rho_d_sc_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0'
if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' d_ed or d_sc'
if (prm%B < 0.0_pReal) extmsg = trim(extmsg)//' B'
if (prm%Q_cl < 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
2020-09-22 22:56:39 +05:30
if (prm%nu_a <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a'
if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w'
if (prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if (prm%V_at <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in dislotungsten, the atomic volume is given as a factor
2020-09-22 22:56:39 +05:30
2020-09-23 04:25:19 +05:30
if (prm%rho_min < 0.0_pReal) extmsg = trim(extmsg)//' rho_min'
2020-09-22 22:56:39 +05:30
if (prm%rho_significant < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant'
if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
2021-06-24 22:44:19 +05:30
if (prm%C_CFL < 0.0_pReal) extmsg = trim(extmsg)//' C_CFL'
2020-09-22 21:41:09 +05:30
if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p_sl'
if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q_sl'
if (prm%f_F < 0.0_pReal .or. prm%f_F > 1.0_pReal) &
2020-08-15 19:32:10 +05:30
extmsg = trim(extmsg)//' f_F'
2020-09-22 21:41:09 +05:30
if (prm%f_ed < 0.0_pReal .or. prm%f_ed > 1.0_pReal) &
2020-08-15 19:32:10 +05:30
extmsg = trim(extmsg)//' f_ed'
2020-09-22 21:41:09 +05:30
if (prm%Q_sol <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol'
if (prm%f_sol <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol'
if (prm%c_sol <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol'
2019-02-21 04:54:35 +05:30
2020-09-22 21:41:09 +05:30
if (prm%chi_GB > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB'
if (prm%chi_surface < 0.0_pReal .or. prm%chi_surface > 1.0_pReal) &
extmsg = trim(extmsg)//' chi_surface'
2020-09-22 21:41:09 +05:30
if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed_mult'
end if slipActive
2019-02-21 04:54:35 +05:30
!--------------------------------------------------------------------------------------------------
! allocate state arrays
2021-04-06 15:08:44 +05:30
Nmembers = count(material_phaseID == ph)
2020-03-16 22:59:15 +05:30
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
'rhoSglScrewPosImmobile','rhoSglScrewNegImmobile', &
'rhoDipEdge ','rhoDipScrew ', &
2020-03-16 23:36:10 +05:30
'gamma ' ]) * prm%sum_N_sl !< "basic" microstructural state variables that are independent from other state variables
sizeDependentState = size([ 'rhoForest ']) * prm%sum_N_sl !< microstructural state variables that depend on other state variables
2019-02-21 04:54:35 +05:30
sizeState = sizeDotState + sizeDependentState &
+ size([ 'velocityEdgePos ','velocityEdgeNeg ', &
2019-03-17 18:05:41 +05:30
'velocityScrewPos ','velocityScrewNeg ', &
2020-03-16 23:36:10 +05:30
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
2019-02-21 04:54:35 +05:30
sizeDeltaState = sizeDotState
2021-03-05 01:46:36 +05:30
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
2019-12-21 16:58:24 +05:30
2021-03-05 01:46:36 +05:30
allocate(geom(ph)%V_0(Nmembers))
2021-02-16 20:36:09 +05:30
call storeGeometry(ph)
2021-02-16 20:36:09 +05:30
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
2020-07-14 00:43:53 +05:30
call IO_error(212,ext_msg='IPneighborhood does not exist')
2021-02-16 20:36:09 +05:30
plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
2021-02-16 20:36:09 +05:30
st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rhoSgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rhoSglMobile => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rhoSglMobile => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rhoSglMobile => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_pos => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_mob_edg_neg => plasticState(ph)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_neg => plasticState(ph)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_neg => plasticState(ph)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_mob_scr_pos => plasticState(ph)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_pos => plasticState(ph)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_pos => plasticState(ph)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_mob_scr_neg => plasticState(ph)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rhoSglImmobile => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSglImmobile => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSglImmobile => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_pos => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_imm_edg_neg => plasticState(ph)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_neg => plasticState(ph)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_neg => plasticState(ph)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_imm_scr_pos => plasticState(ph)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_pos => plasticState(ph)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_pos => plasticState(ph)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_sgl_imm_scr_neg => plasticState(ph)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rhoDip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rhoDip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rhoDip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
del%rho_dip_edg => plasticState(ph)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-02-16 20:36:09 +05:30
stt%rho_dip_scr => plasticState(ph)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
2020-03-16 23:36:10 +05:30
2021-03-05 01:46:36 +05:30
stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-6_pReal)
2021-02-16 20:36:09 +05:30
if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
2020-03-15 12:58:38 +05:30
extmsg = trim(extmsg)//' atol_gamma'
2021-03-05 01:46:36 +05:30
stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
2014-07-08 20:28:23 +05:30
2021-03-05 01:46:36 +05:30
allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
2019-03-17 21:32:08 +05:30
end associate
2021-03-05 01:46:36 +05:30
if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
2019-02-22 14:32:43 +05:30
2020-03-15 03:23:05 +05:30
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
2020-08-15 19:32:10 +05:30
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(nonlocal)')
2020-03-15 03:23:05 +05:30
end do
2020-03-16 23:36:10 +05:30
allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,&
discretization_nIPs,discretization_Nelems), source=0.0_pReal)
2019-02-22 14:32:43 +05:30
! BEGIN DEPRECATED----------------------------------------------------------------------------------
allocate(iRhoU(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0)
2019-02-22 14:32:43 +05:30
2021-02-16 20:36:09 +05:30
do ph = 1, phases%length
2020-08-15 19:32:10 +05:30
2021-02-16 20:36:09 +05:30
if(.not. myPlasticity(ph)) cycle
2021-02-16 20:36:09 +05:30
phase => phases%get(ph)
2021-04-06 15:08:44 +05:30
Nmembers = count(material_phaseID == ph)
2020-08-15 19:32:10 +05:30
l = 0
do t = 1,4
2021-02-16 20:36:09 +05:30
do s = 1,param(ph)%sum_N_sl
2020-08-15 19:32:10 +05:30
l = l + 1
2021-02-16 20:36:09 +05:30
iRhoU(s,t,ph) = l
end do
end do
2021-02-16 20:36:09 +05:30
l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest
2020-08-15 19:32:10 +05:30
do t = 1,4
2021-02-16 20:36:09 +05:30
do s = 1,param(ph)%sum_N_sl
2020-08-15 19:32:10 +05:30
l = l + 1
2021-02-16 20:36:09 +05:30
iV(s,t,ph) = l
end do
end do
2020-08-15 19:32:10 +05:30
do t = 1,2
2021-02-16 20:36:09 +05:30
do s = 1,param(ph)%sum_N_sl
2020-08-15 19:32:10 +05:30
l = l + 1
2021-02-16 20:36:09 +05:30
iD(s,t,ph) = l
end do
end do
2021-02-16 20:36:09 +05:30
if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) &
error stop 'state indices not properly set (nonlocal)'
end do
2020-08-15 19:32:10 +05:30
end function plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
2021-04-25 11:36:52 +05:30
module subroutine nonlocal_dependentState(ph, en, ip, el)
2020-03-16 21:59:47 +05:30
integer, intent(in) :: &
ph, &
2021-04-25 11:36:52 +05:30
en, &
2020-03-16 21:59:47 +05:30
ip, &
el
2019-03-17 21:32:08 +05:30
integer :: &
no, & !< neighbor offset
neighbor_el, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point
c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index
dir, &
2020-07-02 00:52:05 +05:30
n
2019-03-17 21:32:08 +05:30
real(pReal) :: &
FVsize, &
nRealNeighbors, & ! number of really existing neighbors
mu, &
nu
2019-03-17 21:32:08 +05:30
integer, dimension(2) :: &
neighbors
real(pReal), dimension(2) :: &
rhoExcessGradient, &
rhoExcessGradient_over_rho, &
rhoTotal
real(pReal), dimension(3) :: &
rhoExcessDifferences, &
normal_latticeConf
real(pReal), dimension(3,3) :: &
invFe, & !< inverse of elastic deformation gradient
invFp, & !< inverse of plastic deformation gradient
connections, &
invConnections
2019-06-07 13:50:56 +05:30
real(pReal), dimension(3,nIPneighbors) :: &
2019-03-17 21:32:08 +05:30
connection_latticeConf
real(pReal), dimension(2,param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
rhoExcess
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
rho_edg_delta, &
rho_scr_delta
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho, &
rho0, &
rho_neighbor0
real(pReal), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(param(ph)%sum_N_sl,nIPneighbors) :: &
2019-03-17 21:32:08 +05:30
rho_edg_delta_neighbor, &
rho_scr_delta_neighbor
2020-03-16 23:36:10 +05:30
real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: &
2019-03-17 21:32:08 +05:30
neighbor_rhoExcess, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,param(ph)%sum_N_sl,2) :: &
2019-03-17 21:32:08 +05:30
m ! direction of dislocation motion
2021-07-27 02:26:40 +05:30
associate(prm => param(ph),dst => dependentState(ph), stt => state(ph))
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
2021-04-25 11:36:52 +05:30
rho = getRho(ph,en)
2021-04-25 11:36:52 +05:30
stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+ matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
! coefficients are corrected for the line tension effect
2019-03-17 21:32:08 +05:30
! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
if (any(phase_lattice(ph) == ['cI','cF'])) then
2020-09-22 21:41:09 +05:30
myInteractionMatrix = prm%h_sl_sl &
* spread(( 1.0_pReal - prm%f_F &
+ prm%f_F &
2021-04-25 11:36:52 +05:30
* log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,en),prm%rho_significant))) &
2020-09-22 21:41:09 +05:30
/ log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
2019-03-17 21:32:08 +05:30
else
2020-09-22 21:41:09 +05:30
myInteractionMatrix = prm%h_sl_sl
end if
2019-03-17 21:32:08 +05:30
dst%tau_pass(:,en) = mu * prm%b_sl &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities
!*** zero for material points of local plasticity
2019-02-20 19:24:26 +05:30
!#################################################################################################
! ToDo: MD: this is most likely only correct for F_i = I
!#################################################################################################
2021-04-25 11:36:52 +05:30
rho0 = getRho0(ph,en)
if (plasticState(ph)%nonlocal .and. prm%shortRangeStressCorrection) then
2021-04-25 11:36:52 +05:30
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
2020-03-16 13:06:19 +05:30
rhoExcess(1,:) = rho_edg_delta
rhoExcess(2,:) = rho_scr_delta
2021-04-25 11:36:52 +05:30
FVsize = geom(ph)%V_0(en) ** (1.0_pReal/3.0_pReal)
2019-03-17 21:32:08 +05:30
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
2019-03-17 21:32:08 +05:30
nRealNeighbors = 0.0_pReal
neighbor_rhoTotal = 0.0_pReal
2019-06-07 13:50:56 +05:30
do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
2019-06-14 12:47:05 +05:30
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
2019-03-17 21:32:08 +05:30
if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseAt(1,neighbor_el) == ph) then
2019-03-17 21:32:08 +05:30
nRealNeighbors = nRealNeighbors + 1.0_pReal
2021-02-14 21:59:23 +05:30
rho_neighbor0 = getRho0(ph,no)
rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
2019-09-28 02:37:03 +05:30
connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
- discretization_IPcoords(1:3,el+neighbor_ip-1))
2019-06-07 14:03:49 +05:30
normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el))
if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
2019-03-17 21:32:08 +05:30
else
! local neighbor or different lattice structure or different constitution instance -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal
rho_edg_delta_neighbor(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta
end if
else
2019-03-17 21:32:08 +05:30
! free surface -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal
rho_edg_delta_neighbor(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta
end if
end do
2019-03-17 21:32:08 +05:30
neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor
neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor
2019-03-17 21:32:08 +05:30
!* loop through the slip systems and calculate the dislocation gradient by
!* 1. interpolation of the excess density in the neighorhood
!* 2. interpolation of the dead dislocation density in the central volume
2020-03-16 13:06:19 +05:30
m(1:3,:,1) = prm%slip_direction
m(1:3,:,2) = -prm%slip_transverse
2020-03-16 23:36:10 +05:30
do s = 1,prm%sum_N_sl
2019-02-21 23:48:06 +05:30
! gradient from interpolation of neighboring excess density ...
2019-03-17 21:32:08 +05:30
do c = 1,2
do dir = 1,3
neighbors(1) = 2 * dir - 1
neighbors(2) = 2 * dir
connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) &
- connection_latticeConf(1:3,neighbors(2))
rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) &
- neighbor_rhoExcess(c,s,neighbors(2))
end do
2019-03-17 21:32:08 +05:30
invConnections = math_inv33(connections)
if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
end do
2019-03-17 21:32:08 +05:30
! ... plus gradient from deads ...
rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize
rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize
2019-03-17 21:32:08 +05:30
! ... normalized with the total density ...
rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pReal + nRealNeighbors)
rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pReal + nRealNeighbors)
2019-03-17 21:32:08 +05:30
rhoExcessGradient_over_rho = 0.0_pReal
2020-03-16 13:06:19 +05:30
where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
2019-03-17 21:32:08 +05:30
! ... gives the local stress correction when multiplied with a factor
dst%tau_back(s,en) = - mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - nu) &
2020-03-16 21:59:47 +05:30
+ rhoExcessGradient_over_rho(2))
end do
end if
end associate
2019-02-21 23:48:06 +05:30
2021-01-26 16:47:00 +05:30
end subroutine nonlocal_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
2021-01-26 12:03:04 +05:30
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
2021-04-25 11:36:52 +05:30
Mp,Temperature,ph,en)
2020-03-16 21:59:47 +05:30
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp
2019-03-17 21:32:08 +05:30
integer, intent(in) :: &
ph, &
2021-04-25 11:36:52 +05:30
en
2019-03-17 21:32:08 +05:30
real(pReal), intent(in) :: &
2020-03-16 21:59:47 +05:30
Temperature !< temperature
2019-03-17 21:32:08 +05:30
real(pReal), dimension(3,3), intent(in) :: &
Mp
2020-03-16 21:59:47 +05:30
!< derivative of Lp with respect to Mp
2019-03-17 21:32:08 +05:30
integer :: &
2021-07-22 03:41:16 +05:30
i, j, k, l, &
2019-03-17 21:32:08 +05:30
t, & !< dislocation type
s !< index of my current slip system
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
2019-03-17 21:32:08 +05:30
rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
2019-03-17 21:32:08 +05:30
v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
tau, & !< resolved shear stress including backstress terms
2021-07-21 20:04:33 +05:30
dot_gamma !< shear rate
2021-07-27 02:26:40 +05:30
associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph))
!*** shortcut to state variables
2021-04-25 11:36:52 +05:30
rho = getRho(ph,en)
2019-03-17 21:32:08 +05:30
rhoSgl = rho(:,sgl)
2021-07-22 03:41:16 +05:30
do s = 1,prm%sum_N_sl
2021-07-22 03:35:49 +05:30
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
2019-03-17 21:32:08 +05:30
tauNS(s,1) = tau(s)
tauNS(s,2) = tau(s)
if (tau(s) > 0.0_pReal) then
2021-07-21 19:59:57 +05:30
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
2019-03-17 21:32:08 +05:30
else
2021-07-21 19:59:57 +05:30
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
end if
end do
2021-04-25 11:36:52 +05:30
tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
tau = tau + dst%tau_back(:,en)
! edges
2020-03-16 19:52:44 +05:30
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
2021-04-25 11:36:52 +05:30
tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph)
2020-03-16 13:06:19 +05:30
v(:,2) = v(:,1)
dv_dtau(:,2) = dv_dtau(:,1)
dv_dtauNS(:,2) = dv_dtauNS(:,1)
2019-03-17 21:32:08 +05:30
!screws
2020-03-17 12:05:25 +05:30
if (prm%nonSchmidActive) then
2019-03-17 21:32:08 +05:30
do t = 3,4
2020-03-16 19:52:44 +05:30
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
2021-04-25 11:36:52 +05:30
tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
end do
else
v(:,3:4) = spread(v(:,1),2,2)
dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2)
dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
end if
2019-12-01 13:25:24 +05:30
2021-04-25 11:36:52 +05:30
stt%v(:,en) = pack(v,.true.)
2019-03-17 21:32:08 +05:30
!*** Bauschinger effect
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
2019-03-17 21:32:08 +05:30
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
2021-07-21 20:04:33 +05:30
dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
2019-03-17 21:32:08 +05:30
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
2021-07-22 03:41:16 +05:30
do s = 1,prm%sum_N_sl
2021-07-22 03:35:49 +05:30
Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
2019-03-17 21:32:08 +05:30
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
2021-07-22 03:35:49 +05:30
+ prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
2020-09-22 21:41:09 +05:30
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
2021-07-22 03:35:49 +05:30
+ prm%P_sl(i,j,s) &
2021-07-21 19:59:57 +05:30
* (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
end do
2019-03-17 21:32:08 +05:30
end associate
2021-01-26 12:03:04 +05:30
end subroutine nonlocal_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
2021-04-25 11:36:52 +05:30
module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
2020-03-16 21:59:47 +05:30
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
2019-03-17 21:32:08 +05:30
integer, intent(in) :: &
ph, &
2021-04-25 11:36:52 +05:30
en
2019-03-17 21:32:08 +05:30
integer :: &
c, & ! character of dislocation
t, & ! type of dislocation
2020-07-02 00:52:05 +05:30
s ! index of my current slip system
real(pReal) :: &
mu, &
nu
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho ! current dislocation densities
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
2019-03-17 21:32:08 +05:30
v ! dislocation glide velocity
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
tau ! current resolved shear stress
real(pReal), dimension(param(ph)%sum_N_sl,2) :: &
2019-03-17 21:32:08 +05:30
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
dUpper, & ! current maximum stable dipole distance for edges and screws
dUpperOld, & ! old maximum stable dipole distance for edges and screws
deltaDUpper ! change in maximum stable dipole distance for edges and screws
2021-07-27 02:26:40 +05:30
associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
!*** shortcut to state variables
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
2021-04-25 11:36:52 +05:30
rho = getRho(ph,en)
2019-03-17 21:32:08 +05:30
rhoDip = rho(:,dip)
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** dislocation remobilization (bauschinger effect)
where(rho(:,imm) * v < 0.0_pReal)
deltaRhoRemobilization(:,mob) = abs(rho(:,imm))
deltaRhoRemobilization(:,imm) = - rho(:,imm)
rho(:,mob) = rho(:,mob) + abs(rho(:,imm))
rho(:,imm) = 0.0_pReal
elsewhere
deltaRhoRemobilization(:,mob) = 0.0_pReal
deltaRhoRemobilization(:,imm) = 0.0_pReal
endwhere
2019-03-17 21:32:08 +05:30
deltaRhoRemobilization(:,dip) = 0.0_pReal
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** calculate dipole formation and dissociation by stress change
2019-03-17 21:32:08 +05:30
!*** calculate limits for stable dipole height
2020-03-16 23:36:10 +05:30
do s = 1,prm%sum_N_sl
2021-07-22 03:35:49 +05:30
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en)
2019-03-17 21:32:08 +05:30
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
end do
dUpper(:,1) = mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - nu) * abs(tau))
dUpper(:,2) = mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
2019-03-17 21:32:08 +05:30
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
2020-03-16 13:06:19 +05:30
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
2019-03-17 21:32:08 +05:30
where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) &
2020-03-16 13:06:19 +05:30
dUpper(:,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(:,2))
2019-03-17 21:32:08 +05:30
dUpper = max(dUpper,prm%minDipoleHeight)
deltaDUpper = dUpper - dUpperOld
2019-03-17 21:32:08 +05:30
!*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal
2021-07-22 03:41:16 +05:30
forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pReal .and. &
2019-03-17 21:32:08 +05:30
dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) &
deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) &
2020-03-16 13:06:19 +05:30
/ (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
2020-03-16 13:06:19 +05:30
forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c)
2021-04-25 11:36:52 +05:30
plasticState(ph)%deltaState(:,en) = 0.0_pReal
2021-07-22 03:41:16 +05:30
del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl])
2019-03-17 21:32:08 +05:30
end associate
end subroutine plastic_nonlocal_deltaState
2019-02-20 05:11:44 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
2021-01-26 16:47:00 +05:30
module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
2021-07-22 03:41:16 +05:30
ph,en,ip,el)
2020-03-16 14:19:59 +05:30
real(pReal), dimension(3,3), intent(in) :: &
2019-03-17 21:32:08 +05:30
Mp !< MandelStress
2020-03-16 14:19:59 +05:30
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
2021-04-25 11:36:52 +05:30
en, &
2020-03-16 14:19:59 +05:30
ip, & !< current integration point
el !< current element number
2019-03-17 21:32:08 +05:30
integer :: &
c, & !< character of dislocation
t, & !< type of dislocation
2020-07-02 00:52:05 +05:30
s !< index of my current slip system
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho, &
2020-02-11 10:11:10 +05:30
rho0, & !< dislocation density at beginning of time step
2019-03-17 21:32:08 +05:30
rhoDot, & !< density evolution
rhoDotMultiplication, & !< density evolution by multiplication
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
2019-03-17 21:32:08 +05:30
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
2020-02-11 10:11:10 +05:30
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
2019-03-17 21:32:08 +05:30
v, & !< current dislocation glide velocity
v0, &
2021-07-21 20:04:33 +05:30
dot_gamma !< shear rates
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
tau, & !< current resolved shear stress
2021-06-26 16:23:58 +05:30
v_climb !< climb velocity of edge dipoles
real(pReal), dimension(param(ph)%sum_N_sl,2) :: &
2019-03-17 21:32:08 +05:30
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws
real(pReal) :: &
D_SD, &
mu, &
nu
2019-03-17 21:32:08 +05:30
if (timestep <= 0.0_pReal) then
2020-03-16 14:19:59 +05:30
plasticState(ph)%dotState = 0.0_pReal
2019-03-17 21:32:08 +05:30
return
end if
2021-07-27 02:26:40 +05:30
associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
2019-03-17 21:32:08 +05:30
tau = 0.0_pReal
2021-07-21 20:04:33 +05:30
dot_gamma = 0.0_pReal
2021-04-25 11:36:52 +05:30
rho = getRho(ph,en)
2019-03-17 21:32:08 +05:30
rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip)
2021-04-25 11:36:52 +05:30
rho0 = getRho0(ph,en)
2020-02-11 10:11:10 +05:30
my_rhoSgl0 = rho0(:,sgl)
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
2021-07-21 20:04:33 +05:30
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
2021-06-24 18:55:43 +05:30
! limits for stable dipole height
2021-07-22 03:41:16 +05:30
do s = 1,prm%sum_N_sl
2021-07-22 03:35:49 +05:30
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en)
2019-03-17 21:32:08 +05:30
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
end do
2019-03-17 22:29:01 +05:30
dLower = prm%minDipoleHeight
dUpper(:,1) = mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - nu) * abs(tau))
dUpper(:,2) = mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
2019-03-17 22:29:01 +05:30
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
2020-03-16 13:06:19 +05:30
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
2019-03-17 22:29:01 +05:30
where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) &
2020-03-16 13:06:19 +05:30
dUpper(:,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(:,2))
2019-03-17 22:29:01 +05:30
2019-03-17 21:32:08 +05:30
dUpper = max(dUpper,dLower)
2021-06-24 18:55:43 +05:30
! dislocation multiplication
2019-03-17 21:32:08 +05:30
rhoDotMultiplication = 0.0_pReal
isBCC: if (phase_lattice(ph) == 'cI') then
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pReal)
2021-07-21 20:04:33 +05:30
rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
2021-04-25 11:36:52 +05:30
* sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path
2019-03-17 21:32:08 +05:30
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
2021-07-21 20:04:33 +05:30
rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
2021-04-25 11:36:52 +05:30
* sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path
2019-03-17 21:32:08 +05:30
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall
2019-03-17 21:32:08 +05:30
else isBCC
2020-03-16 13:06:19 +05:30
rhoDotMultiplication(:,1:4) = spread( &
2021-07-21 20:04:33 +05:30
(sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
2021-06-26 16:23:58 +05:30
* sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26
end if isBCC
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
2019-03-17 21:32:08 +05:30
!****************************************************************************
2021-06-24 18:55:43 +05:30
! dipole formation and annihilation
2021-06-24 18:55:43 +05:30
! formation by glide
2019-03-17 21:32:08 +05:30
do c = 1,2
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
2021-07-21 20:04:33 +05:30
* ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile
2020-03-16 13:06:19 +05:30
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
2021-07-21 20:04:33 +05:30
* ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile
2020-03-16 13:06:19 +05:30
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
2021-07-21 20:04:33 +05:30
* rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile
2020-03-16 13:06:19 +05:30
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
2021-07-21 20:04:33 +05:30
* rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile
2020-03-16 13:06:19 +05:30
rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
+ abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
- rhoDotSingle2DipoleGlide(:,2*c-1) &
- rhoDotSingle2DipoleGlide(:,2*c)
end do
2021-06-24 18:55:43 +05:30
! athermal annihilation
2019-03-17 21:32:08 +05:30
rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) &
2020-09-22 21:41:09 +05:30
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
2021-07-21 20:04:33 +05:30
* ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single
+ 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent
2020-03-16 14:19:59 +05:30
! annihilated screw dipoles leave edge jogs behind on the colinear system
if (phase_lattice(ph) == 'cF') &
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, prm%colinearSystem(s) > 0) &
2019-03-17 21:32:08 +05:30
rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
2021-04-25 11:36:52 +05:30
* 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
2021-06-24 18:55:43 +05:30
! thermally activated annihilation of edge dipoles by climb
2019-03-17 21:32:08 +05:30
rhoDotThermalAnnihilation = 0.0_pReal
D_SD = prm%D_0 * exp(-prm%Q_cl / (K_B * Temperature)) ! eq. 3.53
v_climb = D_SD * mu * prm%V_at &
/ (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54
2021-07-22 03:41:16 +05:30
forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
2021-06-26 16:23:58 +05:30
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
2019-03-17 21:32:08 +05:30
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
2020-09-22 21:41:09 +05:30
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
2019-03-17 21:32:08 +05:30
2021-04-25 11:36:52 +05:30
rhoDot = rhoDotFlux(timestep, ph,en,ip,el) &
2019-03-17 21:32:08 +05:30
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
+ rhoDotThermalAnnihilation
2020-03-16 13:06:19 +05:30
if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) &
.or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
print'(a)', '<< CONST >> enforcing cutback !!!'
end if
#endif
2020-03-16 14:19:59 +05:30
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
2019-03-17 21:32:08 +05:30
else
2021-04-25 11:36:52 +05:30
dot%rho(:,en) = pack(rhoDot,.true.)
2021-07-21 20:04:33 +05:30
dot%gamma(:,en) = sum(dot_gamma,2)
end if
2019-03-17 21:32:08 +05:30
end associate
2021-01-26 16:47:00 +05:30
end subroutine nonlocal_dotState
2019-03-17 21:32:08 +05:30
2020-03-16 19:52:44 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
#if __INTEL_COMPILER >= 2020
non_recursive function rhoDotFlux(timestep,ph,en,ip,el)
#else
2021-04-25 11:36:52 +05:30
function rhoDotFlux(timestep,ph,en,ip,el)
#endif
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
2021-04-25 11:36:52 +05:30
en, &
ip, & !< current integration point
el !< current element number
2020-07-03 20:15:11 +05:30
integer :: &
neighbor_ph, & !< phase of my neighbor's plasticity
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
n, & !< index of my current neighbor
neighbor_el, & !< element number of my neighbor
neighbor_ip, & !< integration point of my neighbor
2021-04-25 11:36:52 +05:30
neighbor_n, & !< neighbor index pointing to en when looking from my neighbor
opposite_neighbor, & !< index of my opposite neighbor
opposite_ip, & !< ip of my opposite neighbor
opposite_el, & !< element index of my opposite neighbor
2021-04-25 11:36:52 +05:30
opposite_n, & !< neighbor index pointing to en when looking from my opposite neighbor
t, & !< type of dislocation
no,& !< neighbor offset shortcut
np,& !< neighbor phase shortcut
topp, & !< type of dislocation with opposite sign to t
s !< index of my current slip system
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho, &
rho0, & !< dislocation density at beginning of time step
2020-05-03 14:24:57 +05:30
rhoDotFlux !< density evolution by flux
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity
v0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip
2021-07-21 20:04:33 +05:30
dot_gamma !< shear rates
real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: &
m !< direction of dislocation motion
real(pReal), dimension(3,3) :: &
my_F, & !< my total deformation gradient
neighbor_F, & !< total deformation gradient of my neighbor
my_Fe, & !< my elastic deformation gradient
neighbor_Fe, & !< elastic deformation gradient of my neighbor
2021-04-25 11:36:52 +05:30
Favg !< average total deformation gradient of en and my neighbor
real(pReal), dimension(3) :: &
2021-04-25 11:36:52 +05:30
normal_neighbor2me, & !< interface normal pointing from my neighbor to en in neighbor's lattice configuration
normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to en in shared deformed configuration
normal_me2neighbor, & !< interface normal pointing from en to my neighbor in my lattice configuration
normal_me2neighbor_defConf !< interface normal pointing from en to my neighbor in shared deformed configuration
real(pReal) :: &
area, & !< area of the current interface
transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point
2020-05-03 14:24:57 +05:30
lineLength !< dislocation line length leaving the current interface
associate(prm => param(ph), &
2021-07-27 02:26:40 +05:30
dst => dependentState(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
2021-07-21 20:04:33 +05:30
dot_gamma = 0.0_pReal
2021-04-25 11:36:52 +05:30
rho = getRho(ph,en)
rhoSgl = rho(:,sgl)
2021-04-25 11:36:52 +05:30
rho0 = getRho0(ph,en)
my_rhoSgl0 = rho0(:,sgl)
2021-04-25 11:36:52 +05:30
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
2021-07-21 20:04:33 +05:30
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
2021-04-25 11:36:52 +05:30
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity)
rhoDotFlux = 0.0_pReal
if (plasticState(ph)%nonlocal) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
2021-07-21 20:04:33 +05:30
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
2021-06-24 22:44:19 +05:30
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
2021-07-21 20:04:33 +05:30
maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
2021-06-24 22:44:19 +05:30
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!'
end if
#endif
2020-05-03 15:34:57 +05:30
rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback
return
end if
!*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!!
!*** opposite sign to our t vector in the (s,t,n) triplet !!!
m(1:3,:,1) = prm%slip_direction
m(1:3,:,2) = -prm%slip_direction
m(1:3,:,3) = -prm%slip_transverse
m(1:3,:,4) = prm%slip_transverse
2021-04-25 11:36:52 +05:30
my_F = phase_mechanical_F(ph)%data(1:3,1:3,en)
my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
neighbors: do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
neighbor_n = IPneighborhood(3,n,ip,el)
np = material_phaseAt(1,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_ph = material_phaseAt(1,neighbor_el)
2021-02-09 03:51:53 +05:30
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average
Favg = my_F
end if
neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below
!* FLUX FROM MY NEIGHBOR TO ME
!* This is only considered, if I have a neighbor of nonlocal plasticity
!* (also nonlocal constitutive law with local properties) that is at least a little bit
!* compatible.
!* If it's not at all compatible, no flux is arriving, because everything is dammed in front of
!* my neighbor's interface.
!* The entering flux from my neighbor will be distributed on my slip systems according to the
!* compatibility
if (neighbor_n > 0) then
if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID .and. &
any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then
forall (s = 1:ns, t = 1:4)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal)
endforall
2020-09-23 04:25:19 +05:30
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
2020-09-22 21:41:09 +05:30
.or. neighbor_rhoSgl0 < prm%rho_significant) &
neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
2021-04-25 11:36:52 +05:30
IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => en)
normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
/ math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor
area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2)
2021-04-25 11:36:52 +05:30
if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to en == entering flux for en
.and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) &
rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to equally signed mobile dislocation type
where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) &
rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type
end if
end do
end do
end if; end if
!* FLUX FROM ME TO MY NEIGHBOR
!* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties).
2021-04-25 11:36:52 +05:30
!* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to en.
!* So the net flux in the direction of my neighbor is equal to zero:
!* leaving flux to neighbor == entering flux from opposite neighbor
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
if (opposite_n > 0) then
if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then
normal_me2neighbor_defConf = math_det33(Favg) &
2021-04-25 11:36:52 +05:30
* matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
/ math_det33(my_Fe) ! interface normal in my lattice configuration
area = IParea(n,ip,el) * norm2(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
2021-04-25 11:36:52 +05:30
if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
end if
lineLength = my_rhoSgl0(s,t) * v0(s,t) &
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
end if
end do
end do
end if; end if
end do neighbors
end if
end associate
2020-05-03 16:39:49 +05:30
end function rhoDotFlux
2019-03-17 21:32:08 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-17 22:29:01 +05:30
!> @brief Compatibility update
!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip
2019-03-17 21:32:08 +05:30
! 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.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
2021-05-22 22:12:18 +05:30
type(tRotationContainer), dimension(:), intent(in) :: &
2020-03-16 14:19:59 +05:30
orientation ! crystal orientation
2019-03-17 21:32:08 +05:30
integer, intent(in) :: &
ph, &
2019-03-17 21:32:08 +05:30
i, &
e
2019-03-17 21:32:08 +05:30
integer :: &
n, & ! neighbor index
2021-04-25 11:36:52 +05:30
en, &
2019-03-17 21:32:08 +05:30
neighbor_e, & ! element index of my neighbor
neighbor_i, & ! integration point index of my neighbor
2021-03-05 01:46:36 +05:30
neighbor_me, &
2019-03-17 21:32:08 +05:30
neighbor_phase, &
ns, & ! number of active slip systems
2021-04-25 11:36:52 +05:30
s1, & ! slip system index (en)
2019-03-17 21:32:08 +05:30
s2 ! slip system index (my neighbor)
real(pReal), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: &
my_compatibility ! my_compatibility for current element and ip
2019-03-17 21:32:08 +05:30
real(pReal) :: &
my_compatibilitySum, &
thresholdValue, &
nThresholdValues
logical, dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
belowThreshold
2019-12-02 17:28:23 +05:30
type(rotation) :: mis
2019-03-17 21:32:08 +05:30
associate(prm => param(ph))
2020-03-16 23:36:10 +05:30
ns = prm%sum_N_sl
2021-04-25 11:36:52 +05:30
en = material_phaseMemberAt(1,i,e)
2019-03-17 21:32:08 +05:30
!*** start out fully compatible
my_compatibility = 0.0_pReal
2020-03-16 13:06:19 +05:30
forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
2020-03-16 13:06:19 +05:30
neighbors: do n = 1,nIPneighbors
neighbor_e = IPneighborhood(1,n,i,e)
neighbor_i = IPneighborhood(2,n,i,e)
neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
2020-03-16 14:07:25 +05:30
neighbor_phase = material_phaseAt(1,neighbor_e)
2019-03-17 21:32:08 +05:30
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
2020-03-16 14:07:25 +05:30
!* FREE SURFACE
2020-09-22 21:41:09 +05:30
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
2020-03-16 14:12:58 +05:30
elseif (neighbor_phase /= ph) then
2020-03-16 14:07:25 +05:30
!* PHASE BOUNDARY
if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) &
2020-03-16 13:06:19 +05:30
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
2020-09-22 21:41:09 +05:30
elseif (prm%chi_GB >= 0.0_pReal) then
!* GRAIN BOUNDARY
2021-07-24 18:33:26 +05:30
if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
plasticState(neighbor_phase)%nonlocal) &
2020-09-22 21:41:09 +05:30
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
2019-03-17 21:32:08 +05:30
else
2020-03-16 14:07:25 +05:30
!* GRAIN BOUNDARY ?
!* Compatibility defined by relative orientation of slip systems:
!* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection.
!* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
!* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
!* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
!* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
!* All values below the threshold are set to zero.
2021-05-22 22:12:18 +05:30
mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me))
2019-03-17 21:32:08 +05:30
mySlipSystems: do s1 = 1,ns
neighborSlipSystems: do s2 = 1,ns
2020-03-16 13:06:19 +05:30
my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
2019-12-02 17:28:23 +05:30
mis%rotate(prm%slip_normal(1:3,s2))) &
2020-03-16 13:06:19 +05:30
* abs(math_inner(prm%slip_direction(1:3,s1), &
2019-12-02 17:28:23 +05:30
mis%rotate(prm%slip_direction(1:3,s2))))
my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
2019-12-02 17:28:23 +05:30
mis%rotate(prm%slip_normal(1:3,s2)))) &
2020-03-16 13:06:19 +05:30
* abs(math_inner(prm%slip_direction(1:3,s1), &
2019-12-02 17:28:23 +05:30
mis%rotate(prm%slip_direction(1:3,s2))))
end do neighborSlipSystems
2019-03-17 21:32:08 +05:30
my_compatibilitySum = 0.0_pReal
belowThreshold = .true.
2020-03-16 13:06:19 +05:30
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold))
thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive
nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal)
where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false.
2019-03-17 21:32:08 +05:30
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
2020-03-16 13:06:19 +05:30
where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) &
my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
my_compatibility(:,:,s1,n))
2019-03-17 21:32:08 +05:30
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
end do
2020-03-16 14:12:58 +05:30
2020-03-16 13:06:19 +05:30
where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
2020-03-16 14:12:58 +05:30
end do mySlipSystems
end if
end do neighbors
2020-03-16 13:06:19 +05:30
compatibility(:,:,:,:,i,e) = my_compatibility
2019-03-17 21:32:08 +05:30
end associate
end subroutine plastic_nonlocal_updateCompatibility
!--------------------------------------------------------------------------------------------------
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_results(ph,group)
integer, intent(in) :: ph
character(len=*),intent(in) :: group
2020-02-14 13:56:26 +05:30
integer :: ou
2021-07-27 02:26:40 +05:30
associate(prm => param(ph),dst => dependentState(ph),stt=>state(ph))
do ou = 1,size(prm%output)
select case(trim(prm%output(ou)))
case('rho_u_ed_pos')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_mob_edg_pos,group,trim(prm%output(ou)), &
'positive mobile edge density','1/m²', prm%systems_sl)
case('rho_b_ed_pos')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_imm_edg_pos,group,trim(prm%output(ou)), &
'positive immobile edge density','1/m²', prm%systems_sl)
case('rho_u_ed_neg')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_mob_edg_neg,group,trim(prm%output(ou)), &
'negative mobile edge density','1/m²', prm%systems_sl)
case('rho_b_ed_neg')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_imm_edg_neg,group,trim(prm%output(ou)), &
'negative immobile edge density','1/m²', prm%systems_sl)
case('rho_d_ed')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_dip_edg,group,trim(prm%output(ou)), &
'edge dipole density','1/m²', prm%systems_sl)
case('rho_u_sc_pos')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_mob_scr_pos,group,trim(prm%output(ou)), &
'positive mobile screw density','1/m²', prm%systems_sl)
case('rho_b_sc_pos')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_imm_scr_pos,group,trim(prm%output(ou)), &
'positive immobile screw density','1/m²', prm%systems_sl)
case('rho_u_sc_neg')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_mob_scr_neg,group,trim(prm%output(ou)), &
'negative mobile screw density','1/m²', prm%systems_sl)
case('rho_b_sc_neg')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_sgl_imm_scr_neg,group,trim(prm%output(ou)), &
'negative immobile screw density','1/m²', prm%systems_sl)
case('rho_d_sc')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_dip_scr,group,trim(prm%output(ou)), &
'screw dipole density','1/m²', prm%systems_sl)
case('rho_f')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%rho_forest,group,trim(prm%output(ou)), &
'forest density','1/m²', prm%systems_sl)
case('v_ed_pos')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%v_edg_pos,group,trim(prm%output(ou)), &
'positive edge velocity','m/s', prm%systems_sl)
case('v_ed_neg')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%v_edg_neg,group,trim(prm%output(ou)), &
'negative edge velocity','m/s', prm%systems_sl)
case('v_sc_pos')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%v_scr_pos,group,trim(prm%output(ou)), &
'positive srew velocity','m/s', prm%systems_sl)
case('v_sc_neg')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%v_scr_neg,group,trim(prm%output(ou)), &
'negative screw velocity','m/s', prm%systems_sl)
case('gamma')
2021-08-08 01:32:44 +05:30
call results_writeDataset(stt%gamma,group,trim(prm%output(ou)), &
'plastic shear','1', prm%systems_sl)
case('tau_pass')
2021-08-08 01:32:44 +05:30
call results_writeDataset(dst%tau_pass,group,trim(prm%output(ou)), &
'passing stress for slip','Pa', prm%systems_sl)
end select
end do
end associate
end subroutine plastic_nonlocal_results
2020-03-16 12:52:36 +05:30
2020-03-17 11:47:40 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density
!--------------------------------------------------------------------------------------------------
2021-05-22 15:33:13 +05:30
subroutine stateInit(ini,phase,Nentries)
2020-03-17 11:47:40 +05:30
2020-03-17 12:05:25 +05:30
type(tInitialParameters) :: &
ini
integer,intent(in) :: &
2020-03-17 11:47:40 +05:30
phase, &
2021-05-22 15:33:13 +05:30
Nentries
2020-03-17 11:47:40 +05:30
integer :: &
e, &
2020-03-17 11:47:40 +05:30
f, &
from, &
upto, &
2021-05-22 15:33:13 +05:30
s
2020-03-17 11:47:40 +05:30
real(pReal), dimension(2) :: &
noise, &
rnd
real(pReal) :: &
meanDensity, &
totalVolume, &
densityBinning, &
minimumIpVolume
2020-08-15 19:32:10 +05:30
associate(stt => state(phase))
2020-03-17 11:47:40 +05:30
2021-05-22 15:33:13 +05:30
if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
totalVolume = sum(geom(phase)%V_0)
minimumIPVolume = minval(geom(phase)%V_0)
densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal)
! fill random material points with dislocation segments until the desired overall density is reached
meanDensity = 0.0_pReal
do while(meanDensity < ini%random_rho_u)
call random_number(rnd)
e = nint(rnd(1)*real(Nentries,pReal) + 0.5_pReal)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume
stt%rhoSglMobile(s,e) = densityBinning
end do
2021-05-22 15:33:13 +05:30
else ! homogeneous distribution with noise
do e = 1, Nentries
do f = 1,size(ini%N_sl,1)
from = 1 + sum(ini%N_sl(1:f-1))
upto = sum(ini%N_sl(1:f))
do s = from,upto
noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), &
math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)]
stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1)
stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
end do
2021-05-22 15:33:13 +05:30
stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f)
stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f)
end do
end do
end if
2020-03-17 11:47:40 +05:30
end associate
end subroutine stateInit
2020-03-16 19:52:44 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics
!--------------------------------------------------------------------------------------------------
2021-06-24 18:55:43 +05:30
pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, ph)
2020-03-16 23:36:10 +05:30
2020-03-16 19:52:44 +05:30
integer, intent(in) :: &
c, & !< dislocation character (1:edge, 2:screw)
ph
2020-03-16 19:52:44 +05:30
real(pReal), intent(in) :: &
2021-06-24 18:55:43 +05:30
T !< T
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
2020-03-16 19:52:44 +05:30
tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including non Schmid effects)
tauThreshold !< threshold shear stress
real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
2020-03-16 19:52:44 +05:30
v, & !< velocity
dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
integer :: &
s !< index of my current slip system
real(pReal) :: &
tauRel_P, &
tauRel_S, &
tauEff, & !< effective shear stress
tPeierls, & !< waiting time in front of a peierls barriers
tSolidSolution, & !< waiting time in front of a solid solution obstacle
dtPeierls_dtau, & !< derivative with respect to resolved shear stress
dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress
2021-06-24 18:55:43 +05:30
lambda_S, & !< mean free distance between two solid solution obstacles
lambda_P, & !< mean free distance between two Peierls barriers
2020-03-16 19:52:44 +05:30
activationVolume_P, & !< volume that needs to be activated to overcome barrier
activationVolume_S, & !< volume that needs to be activated to overcome barrier
activationEnergy_P, & !< energy that is needed to overcome barrier
criticalStress_P, & !< maximum obstacle strength
2021-06-24 19:37:54 +05:30
criticalStress_S !< maximum obstacle strength
2020-03-16 19:52:44 +05:30
associate(prm => param(ph))
2020-03-16 19:52:44 +05:30
v = 0.0_pReal
dv_dtau = 0.0_pReal
dv_dtauNS = 0.0_pReal
2021-03-27 00:24:10 +05:30
do s = 1,prm%sum_N_sl
2020-03-16 19:52:44 +05:30
if (abs(tau(s)) > tauThreshold(s)) then
!* Peierls contribution
2021-06-24 18:55:43 +05:30
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
lambda_P = prm%b_sl(s)
activationVolume_P = prm%w *prm%b_sl(s)**3.0_pReal
2020-03-16 19:52:44 +05:30
criticalStress_P = prm%peierlsStress(s,c)
activationEnergy_P = criticalStress_P * activationVolume_P
2021-06-24 18:55:43 +05:30
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
2020-09-22 21:41:09 +05:30
tPeierls = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_P / (K_B * T) &
2020-03-16 19:52:44 +05:30
* (1.0_pReal - tauRel_P**prm%p)**prm%q)
dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
2021-06-24 19:37:54 +05:30
* (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
0.0_pReal, &
tauEff < criticalStress_P)
2020-03-16 19:52:44 +05:30
2021-06-24 18:55:43 +05:30
! Contribution from solid solution strengthening
2020-03-16 19:52:44 +05:30
tauEff = abs(tau(s)) - tauThreshold(s)
2021-06-24 18:55:43 +05:30
lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
activationVolume_S = prm%f_sol * prm%b_sl(s)**3.0_pReal / sqrt(prm%c_sol)
criticalStress_S = prm%Q_sol / activationVolume_S
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
2020-09-22 21:41:09 +05:30
tSolidSolution = 1.0_pReal / prm%nu_a &
* exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
2021-06-24 19:37:54 +05:30
* (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
0.0_pReal, &
tauEff < criticalStress_S)
2020-03-16 19:52:44 +05:30
!* viscous glide velocity
tauEff = abs(tau(s)) - tauThreshold(s)
2021-06-24 19:37:54 +05:30
2020-03-16 19:52:44 +05:30
v(s) = sign(1.0_pReal,tau(s)) &
2021-06-24 20:01:07 +05:30
/ (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal))
2021-06-24 18:55:43 +05:30
dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P
2021-06-24 19:37:54 +05:30
end if
end do
2020-03-16 19:52:44 +05:30
end associate
end subroutine kinetics
2020-03-16 12:52:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho(ph,en) result(rho)
2020-03-16 12:52:36 +05:30
2021-04-25 11:36:52 +05:30
integer, intent(in) :: ph, en
real(pReal), dimension(param(ph)%sum_N_sl,10) :: rho
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
associate(prm => param(ph))
2020-03-16 12:52:36 +05:30
rho = reshape(state(ph)%rho(:,en),[prm%sum_N_sl,10])
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
! ensure positive densities (not for imm, they have a sign)
rho(:,mob) = max(rho(:,mob),0.0_pReal)
rho(:,dip) = max(rho(:,dip),0.0_pReal)
2020-03-16 12:52:36 +05:30
where(abs(rho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
rho = 0.0_pReal
2020-03-16 12:52:36 +05:30
end associate
end function getRho
!--------------------------------------------------------------------------------------------------
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho0(ph,en) result(rho0)
2020-03-16 12:52:36 +05:30
2021-04-25 11:36:52 +05:30
integer, intent(in) :: ph, en
real(pReal), dimension(param(ph)%sum_N_sl,10) :: rho0
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
associate(prm => param(ph))
2020-03-16 12:52:36 +05:30
rho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10])
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
! ensure positive densities (not for imm, they have a sign)
rho0(:,mob) = max(rho0(:,mob),0.0_pReal)
rho0(:,dip) = max(rho0(:,dip),0.0_pReal)
2020-03-16 12:52:36 +05:30
where (abs(rho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
rho0 = 0.0_pReal
2020-03-16 12:52:36 +05:30
end associate
end function getRho0
subroutine storeGeometry(ph)
integer, intent(in) :: ph
integer :: ce, co
2021-03-01 10:46:16 +05:30
real(pReal), dimension(:), allocatable :: V
2021-03-05 01:46:36 +05:30
2021-03-01 10:46:16 +05:30
V = reshape(IPvolume,[product(shape(IPvolume))])
2021-04-06 15:08:44 +05:30
do ce = 1, size(material_homogenizationEntry,1)
2021-03-01 10:46:16 +05:30
do co = 1, homogenization_maxNconstituents
2021-04-06 15:08:44 +05:30
if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
end do
end do
end subroutine
2021-01-26 05:41:32 +05:30
end submodule nonlocal