2013-11-27 13:34:05 +05:30
!--------------------------------------------------------------------------------------------------
!> @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
!--------------------------------------------------------------------------------------------------
2021-01-27 01:22:48 +05:30
submodule ( phase : plastic ) nonlocal
2019-05-14 15:02:25 +05:30
use geometry_plastic_nonlocal , only : &
2019-06-07 13:50:56 +05:30
nIPneighbors = > geometry_plastic_nonlocal_nIPneighbors , &
2019-05-14 15:02:25 +05:30
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
2020-02-07 16:14:03 +05:30
2021-02-14 19:06:56 +05:30
type :: tGeometry
real ( pReal ) , dimension ( : ) , allocatable :: V_0
end type tGeometry
type ( tGeometry ) , dimension ( : ) , allocatable :: geom
2019-12-04 23:30:56 +05:30
real ( pReal ) , parameter :: &
2020-03-16 03:37:23 +05:30
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
2019-03-16 17:43:48 +05:30
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
2020-02-07 16:14:03 +05:30
2019-12-04 23:30:56 +05:30
real ( pReal ) , dimension ( : , : , : , : , : , : ) , allocatable :: &
compatibility !< slip system compatibility between me and my neighbors
2020-02-07 16:14:03 +05:30
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)
2020-09-22 21:41:09 +05:30
eta , & !< viscosity for dislocation glide in Pa s
nu_a , & !< attack frequency in Hz
chi_surface , & !< transmissivity at free surface
chi_GB , & !< transmissivity at grain boundary (identified by different texture)
f_c , & !< safety factor for CFL flux condition
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 :: &
2019-03-17 21:32:08 +05:30
Schmid , & !< Schmid contribution
nonSchmid_pos , &
nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
integer :: &
2020-03-16 23:36:10 +05:30
sum_N_sl
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 .
2019-03-17 21:32:08 +05:30
end type tParameters
2020-02-07 16:14:03 +05:30
2019-12-04 23:30:56 +05:30
type :: tNonlocalMicrostructure
2019-03-17 21:32:08 +05:30
real ( pReal ) , allocatable , dimension ( : , : ) :: &
2020-02-07 16:14:03 +05:30
tau_pass , &
tau_Back
2019-03-17 21:32:08 +05:30
end type tNonlocalMicrostructure
2020-02-07 16:14:03 +05:30
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 , &
2020-02-07 16:14:03 +05:30
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 , &
2020-02-07 16:14:03 +05:30
state , &
state0
2021-02-14 19:59:10 +05:30
type ( tParameters ) , dimension ( : ) , allocatable :: param !< containers of constitutive parameters
2020-02-07 16:14:03 +05:30
2019-12-04 23:30:56 +05:30
type ( tNonlocalMicrostructure ) , dimension ( : ) , allocatable :: microstructure
2019-02-18 14:58:08 +05:30
2013-11-27 13:34:05 +05:30
contains
2009-08-11 22:01:57 +05:30
2020-09-14 00:30:34 +05:30
2013-10-09 11:42:16 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-14 00:30:34 +05:30
!> @brief Perform module initialization.
2013-10-09 11:42:16 +05:30
!> @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 :: &
2020-10-28 02:03:30 +05:30
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' )
2020-10-28 02:03:30 +05:30
Ninstances = count ( myPlasticity )
if ( Ninstances == 0 ) then
2020-08-15 19:32:10 +05:30
call geometry_plastic_nonlocal_disable
return
endif
2020-12-20 21:41:43 +05:30
2021-02-13 23:27:41 +05:30
print '(/,a)' , ' <<<+- phase:mechanical:plastic:nonlocal init -+>>>'
2021-02-13 17:37:35 +05:30
print '(a,i0)' , ' # phases: ' , Ninstances ; flush ( IO_STDOUT )
2021-03-18 20:06:40 +05:30
print * , 'C. Reuber et al., Acta Materialia 71:333–
2020-09-14 00:30:34 +05:30
print * , 'https://doi.org/10.1016/j.actamat.2014.03.012' / / IO_EOL
2020-08-15 19:32:10 +05:30
2021-03-18 20:06:40 +05:30
print * , 'C. Kords, Dissertation RWTH Aachen, 2014'
2020-09-19 11:50:29 +05:30
print * , 'http://publications.rwth-aachen.de/record/229993'
2020-07-02 00:52:05 +05:30
2021-02-14 19:06:56 +05:30
phases = > config_material % get ( 'phase' )
allocate ( geom ( phases % length ) )
2021-02-14 23:00:57 +05:30
allocate ( param ( phases % length ) )
allocate ( state ( phases % length ) )
allocate ( state0 ( phases % length ) )
allocate ( dotState ( phases % length ) )
allocate ( deltaState ( phases % length ) )
allocate ( microstructure ( 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-14 23:00:57 +05:30
2021-02-16 20:36:09 +05:30
associate ( prm = > param ( ph ) , dot = > dotState ( ph ) , stt = > state ( ph ) , &
st0 = > state0 ( ph ) , del = > deltaState ( ph ) , dst = > microstructure ( ph ) )
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-02-16 20:36:09 +05:30
phase_localPlasticity ( ph ) = . not . pl % contains ( 'nonlocal' )
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
2020-09-14 00:30:34 +05:30
2020-08-15 19:32:10 +05:30
prm % atol_rho = pl % get_asFloat ( 'atol_rho' , defaultVal = 1.0e4_pReal )
2020-02-07 16:14:03 +05:30
2019-02-21 04:54:35 +05:30
! This data is read in already in lattice
2021-02-16 20:36:09 +05:30
prm % mu = lattice_mu ( ph )
prm % nu = lattice_nu ( ph )
2009-08-11 22:01:57 +05:30
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
2020-08-15 19:32:10 +05:30
prm % Schmid = lattice_SchmidMatrix_slip ( ini % N_sl , phase % get_asString ( 'lattice' ) , &
phase % get_asFloat ( 'c/a' , defaultVal = 0.0_pReal ) )
2009-08-11 22:01:57 +05:30
2020-11-29 16:28:39 +05:30
if ( trim ( phase % get_asString ( 'lattice' ) ) == '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-03-26 18:21:09 +05:30
prm % nonSchmid_pos = lattice_nonSchmidMatrix ( ini % N_sl , a , + 1 )
prm % nonSchmid_neg = lattice_nonSchmidMatrix ( ini % N_sl , a , - 1 )
2019-02-21 04:54:35 +05:30
else
2021-03-26 18:21:09 +05:30
prm % nonSchmid_pos = prm % Schmid
prm % nonSchmid_neg = prm % Schmid
2019-02-21 04:54:35 +05:30
endif
2020-09-22 21:41:09 +05:30
prm % h_sl_sl = lattice_interaction_SlipBySlip ( ini % N_sl , &
2021-03-11 22:30:07 +05:30
pl % get_as1dFloat ( 'h_sl_sl' ) , &
2020-09-22 21:41:09 +05:30
phase % get_asString ( 'lattice' ) )
2019-02-21 04:54:35 +05:30
2020-08-15 19:32:10 +05:30
prm % forestProjection_edge = lattice_forestProjection_edge ( ini % N_sl , phase % get_asString ( 'lattice' ) , &
phase % get_asFloat ( 'c/a' , defaultVal = 0.0_pReal ) )
prm % forestProjection_screw = lattice_forestProjection_screw ( ini % N_sl , phase % get_asString ( 'lattice' ) , &
phase % get_asFloat ( 'c/a' , defaultVal = 0.0_pReal ) )
2019-02-21 04:54:35 +05:30
2020-08-15 19:32:10 +05:30
prm % slip_direction = lattice_slip_direction ( ini % N_sl , phase % get_asString ( 'lattice' ) , &
phase % get_asFloat ( 'c/a' , defaultVal = 0.0_pReal ) )
prm % slip_transverse = lattice_slip_transverse ( ini % N_sl , phase % get_asString ( 'lattice' ) , &
phase % get_asFloat ( 'c/a' , defaultVal = 0.0_pReal ) )
prm % slip_normal = lattice_slip_normal ( ini % N_sl , phase % get_asString ( 'lattice' ) , &
phase % get_asFloat ( 'c/a' , defaultVal = 0.0_pReal ) )
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
enddo
enddo
2020-02-07 16:14:03 +05:30
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 ) )
2020-02-07 16:14:03 +05:30
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-02-07 16:14:03 +05:30
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 )
2020-02-07 16:14:03 +05:30
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
2020-02-07 16:14:03 +05:30
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 )
2020-09-22 21:41:09 +05:30
prm % f_c = pl % get_asFloat ( 'f_c' , defaultVal = 2.0_pReal )
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' ) !,'edgejogs'
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' )
prm % eta = pl % get_asFloat ( 'eta' )
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-02-07 16:14:03 +05:30
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-02-07 16:14:03 +05:30
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 % eta < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' eta'
if ( prm % Q_cl < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' Q_cl'
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'
2020-11-29 16:28:39 +05:30
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'
if ( prm % f_c < 0.0_pReal ) extmsg = trim ( extmsg ) / / ' f_c'
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-02-07 16:14:03 +05:30
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-02-07 16:14:03 +05:30
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'
2009-08-11 22:01:57 +05:30
2019-02-21 04:54:35 +05:30
endif slipActive
2019-02-19 15:01:14 +05:30
2019-02-21 04:54:35 +05:30
!--------------------------------------------------------------------------------------------------
! allocate state arrays
2021-03-05 01:46:36 +05:30
Nmembers = count ( material_phaseAt2 == 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 &
2020-02-07 16:14:03 +05:30
+ 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
2020-02-07 16:14:03 +05:30
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-14 19:06:56 +05:30
2021-02-16 20:36:09 +05:30
plasticState ( ph ) % nonlocal = pl % get_asBool ( 'nonlocal' )
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
2020-02-07 16:14:03 +05:30
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 )
2021-02-16 20:36:09 +05:30
plasticState ( ph ) % atol ( 10 * prm % sum_N_sl + 1 : 11 * prm % sum_N_sl ) = pl % get_asFloat ( 'atol_gamma' , defaultVal = 1.0e-2_pReal )
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
plasticState ( ph ) % slipRate = > plasticState ( ph ) % dotState ( 10 * prm % sum_N_sl + 1 : 11 * prm % sum_N_sl , 1 : Nmembers )
2020-02-07 16:14:03 +05:30
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
2019-02-20 18:02:08 +05:30
2021-03-05 01:46:36 +05:30
if ( Nmembers > 0 ) call stateInit ( ini , ph , Nmembers )
2021-02-16 20:36:09 +05:30
plasticState ( ph ) % state0 = plasticState ( ph ) % state
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
2019-02-22 14:32:43 +05:30
enddo
2020-02-07 16:14:03 +05:30
2020-03-16 23:36:10 +05:30
allocate ( compatibility ( 2 , maxval ( param % sum_N_sl ) , maxval ( param % sum_N_sl ) , nIPneighbors , &
2020-10-28 01:57:26 +05:30
discretization_nIPs , discretization_Nelems ) , source = 0.0_pReal )
2020-02-07 16:14:03 +05:30
2019-02-22 14:32:43 +05:30
! BEGIN DEPRECATED----------------------------------------------------------------------------------
2021-02-14 23:00:57 +05:30
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
2020-09-14 00:30:34 +05:30
2021-02-16 20:36:09 +05:30
phase = > phases % get ( ph )
2021-03-05 01:46:36 +05:30
Nmembers = count ( material_phaseAt2 == 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
2019-03-17 21:32:08 +05:30
enddo
2020-08-15 19:32:10 +05:30
enddo
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
2019-03-17 21:32:08 +05:30
enddo
2020-08-15 19:32:10 +05:30
enddo
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
2019-03-17 21:32:08 +05:30
enddo
2020-08-15 19:32:10 +05:30
enddo
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)'
2020-09-14 00:30:34 +05:30
enddo
2009-09-18 21:07:14 +05:30
2020-08-15 19:32:10 +05:30
end function plastic_nonlocal_init
2009-09-18 21:07:14 +05:30
2013-10-09 11:42:16 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
2021-02-14 19:59:10 +05:30
module subroutine nonlocal_dependentState ( ph , me , ip , el )
2020-02-07 16:14:03 +05:30
2020-03-16 21:59:47 +05:30
integer , intent ( in ) :: &
2021-02-14 19:59:10 +05:30
ph , &
2021-01-26 13:09:17 +05:30
me , &
2020-03-16 21:59:47 +05:30
ip , &
el
2020-02-07 16:14:03 +05:30
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
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
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( 2 , param ( ph ) % sum_N_sl ) :: &
2019-03-17 21:32:08 +05:30
rhoExcess
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl ) :: &
2019-03-17 21:32:08 +05:30
rho_edg_delta , &
rho_scr_delta
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 10 ) :: &
2019-03-17 21:32:08 +05:30
rho , &
2020-02-07 17:11:01 +05:30
rho0 , &
rho_neighbor0
2021-02-14 23:00:57 +05:30
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
2021-02-14 23:00:57 +05:30
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
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( 3 , param ( ph ) % sum_N_sl , 2 ) :: &
2019-03-17 21:32:08 +05:30
m ! direction of dislocation motion
2020-02-07 16:14:03 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) , dst = > microstructure ( ph ) , stt = > state ( ph ) )
2020-02-07 16:14:03 +05:30
2021-02-14 21:59:23 +05:30
rho = getRho ( ph , me )
2020-02-07 16:14:03 +05:30
2021-01-26 13:09:17 +05:30
stt % rho_forest ( : , me ) = matmul ( prm % forestProjection_Edge , sum ( abs ( rho ( : , edg ) ) , 2 ) ) &
2020-02-07 16:14:03 +05:30
+ 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)
2020-03-16 13:06:19 +05:30
if ( any ( lattice_structure ( material_phaseAt ( 1 , el ) ) == [ LATTICE_bcc_ID , LATTICE_fcc_ID ] ) ) 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-01-26 13:09:17 +05:30
* log ( 0.35_pReal * prm % b_sl * sqrt ( max ( stt % rho_forest ( : , me ) , 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
2019-03-17 21:32:08 +05:30
endif
2021-01-26 13:09:17 +05:30
dst % tau_pass ( : , me ) = prm % mu * prm % b_sl &
2020-03-16 21:16:16 +05:30
* sqrt ( matmul ( myInteractionMatrix , sum ( abs ( rho ) , 2 ) ) )
2014-06-26 19:23:12 +05:30
2009-08-28 19:20:47 +05:30
!*** calculate the dislocation stress of the neighboring excess dislocation densities
2012-03-12 19:39:37 +05:30
!*** zero for material points of local plasticity
2011-03-29 13:04:33 +05:30
2019-02-20 19:24:26 +05:30
!#################################################################################################
! ToDo: MD: this is most likely only correct for F_i = I
!#################################################################################################
2019-03-16 17:43:48 +05:30
2021-02-14 21:59:23 +05:30
rho0 = getRho0 ( ph , me )
2020-03-16 13:06:19 +05:30
if ( . not . phase_localPlasticity ( material_phaseAt ( 1 , el ) ) . and . prm % shortRangeStressCorrection ) then
2021-02-09 03:51:53 +05:30
invFp = math_inv33 ( phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
invFe = math_inv33 ( phase_mechanical_Fe ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
2020-02-07 16:14:03 +05:30
2020-02-07 17:11:01 +05:30
rho_edg_delta = rho0 ( : , mob_edg_pos ) - rho0 ( : , mob_edg_neg )
rho_scr_delta = rho0 ( : , mob_scr_pos ) - rho0 ( : , mob_scr_neg )
2020-02-07 16:14:03 +05:30
2020-03-16 13:06:19 +05:30
rhoExcess ( 1 , : ) = rho_edg_delta
rhoExcess ( 2 , : ) = rho_scr_delta
2020-02-07 16:14:03 +05:30
2021-02-14 21:59:23 +05:30
FVsize = geom ( ph ) % V_0 ( me ) ** ( 1.0_pReal / 3.0_pReal )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
2020-02-07 16:14:03 +05:30
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
2019-06-06 11:49:29 +05:30
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
2021-02-14 23:00:57 +05:30
if ( material_phaseAt ( 1 , neighbor_el ) == ph ) then
2020-02-07 16:14:03 +05:30
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 )
2020-02-07 16:14:03 +05:30
2020-02-07 17:11:01 +05:30
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 )
2020-02-07 16:14:03 +05:30
2020-02-07 17:11:01 +05:30
neighbor_rhoTotal ( 1 , : , n ) = sum ( abs ( rho_neighbor0 ( : , edg ) ) , 2 )
neighbor_rhoTotal ( 2 , : , n ) = sum ( abs ( rho_neighbor0 ( : , scr ) ) , 2 )
2020-02-07 16:14:03 +05:30
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 ) )
2020-03-16 20:32:14 +05:30
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
endif
2012-01-17 15:56:57 +05:30
else
2019-03-17 21:32:08 +05:30
! free surface -> use central values instead
2012-03-15 15:38:08 +05:30
connection_latticeConf ( 1 : 3 , n ) = 0.0_pReal
2019-03-16 20:59:16 +05:30
rho_edg_delta_neighbor ( : , n ) = rho_edg_delta
2019-03-16 17:43:48 +05:30
rho_scr_delta_neighbor ( : , n ) = rho_scr_delta
2012-01-17 15:56:57 +05:30
endif
2019-03-17 21:32:08 +05:30
enddo
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
neighbor_rhoExcess ( 1 , : , : ) = rho_edg_delta_neighbor
neighbor_rhoExcess ( 2 , : , : ) = rho_scr_delta_neighbor
2020-02-07 16:14:03 +05:30
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-02-07 16:14:03 +05:30
2020-03-16 23:36:10 +05:30
do s = 1 , prm % sum_N_sl
2020-02-07 16:14:03 +05:30
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 ) )
enddo
invConnections = math_inv33 ( connections )
if ( all ( dEq0 ( invConnections ) ) ) call IO_error ( - 1 , ext_msg = 'back stress calculation: inversion error' )
2020-02-07 16:14:03 +05:30
2019-04-03 11:52:04 +05:30
rhoExcessGradient ( c ) = math_inner ( m ( 1 : 3 , s , c ) , matmul ( invConnections , rhoExcessDifferences ) )
2012-01-17 15:56:57 +05:30
enddo
2020-02-07 16:14:03 +05:30
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
2020-02-07 16:14:03 +05:30
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 )
2020-02-07 16:14:03 +05:30
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
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
! ... gives the local stress correction when multiplied with a factor
2021-01-26 13:09:17 +05:30
dst % tau_back ( s , me ) = - prm % mu * prm % b_sl ( s ) / ( 2.0_pReal * PI ) &
2020-03-16 21:59:47 +05:30
* ( rhoExcessGradient_over_rho ( 1 ) / ( 1.0_pReal - prm % nu ) &
+ rhoExcessGradient_over_rho ( 2 ) )
2019-03-17 21:32:08 +05:30
enddo
endif
2011-03-29 13:04:33 +05:30
2017-10-03 18:50:53 +05:30
#ifdef DEBUG
2020-07-02 04:55:24 +05:30
if ( debugConstitutive % extensive &
. and . ( ( debugConstitutive % element == el . and . debugConstitutive % ip == ip ) &
. or . . not . debugConstitutive % selective ) ) then
2020-09-19 14:31:43 +05:30
print '(/,a,i8,1x,i2,1x,i1,/)' , '<< CONST >> nonlocal_microstructure at el ip ' , el , ip
2021-01-26 13:09:17 +05:30
print '(a,/,12x,12(e10.3,1x))' , '<< CONST >> rhoForest' , stt % rho_forest ( : , me )
print '(a,/,12x,12(f10.5,1x))' , '<< CONST >> tauThreshold / MPa' , dst % tau_pass ( : , me ) * 1e-6
print '(a,/,12x,12(f10.5,1x),/)' , '<< CONST >> tauBack / MPa' , dst % tau_back ( : , me ) * 1e-6
2011-03-29 12:57:19 +05:30
endif
#endif
2019-02-17 16:45:46 +05:30
end associate
2019-02-21 23:48:06 +05:30
2021-01-26 16:47:00 +05:30
end subroutine nonlocal_dependentState
2010-02-17 18:51:36 +05:30
2014-06-26 19:23:12 +05:30
2013-10-09 11:42:16 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
2021-01-26 12:03:04 +05:30
module subroutine nonlocal_LpAndItsTangent ( Lp , dLp_dMp , &
2021-02-14 21:59:23 +05:30
Mp , Temperature , ph , me )
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 ) :: &
2021-01-27 03:11:41 +05:30
ph , &
2021-02-14 21:59:23 +05:30
me
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 :: &
ns , & !< short notation for the total number of active slip systems
i , &
j , &
k , &
l , &
t , & !< dislocation type
s !< index of my current slip system
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 8 ) :: &
2019-03-17 21:32:08 +05:30
rhoSgl !< single dislocation densities (including blocked)
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 10 ) :: &
2019-03-17 21:32:08 +05:30
rho
2021-02-14 23:00:57 +05:30
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
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl ) :: &
2019-03-17 21:32:08 +05:30
tau , & !< resolved shear stress including backstress terms
gdotTotal !< shear rate
2020-02-07 16:14:03 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) , dst = > microstructure ( ph ) , &
stt = > state ( ph ) )
2020-03-16 23:36:10 +05:30
ns = prm % sum_N_sl
2020-02-07 16:14:03 +05:30
!*** shortcut to state variables
2021-02-14 21:59:23 +05:30
rho = getRho ( ph , me )
2019-03-17 21:32:08 +05:30
rhoSgl = rho ( : , sgl )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
do s = 1 , ns
2020-03-16 23:13:04 +05:30
tau ( s ) = math_tensordot ( Mp , prm % Schmid ( 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
2020-03-16 23:13:04 +05:30
tauNS ( s , 3 ) = math_tensordot ( Mp , + prm % nonSchmid_pos ( 1 : 3 , 1 : 3 , s ) )
tauNS ( s , 4 ) = math_tensordot ( Mp , - prm % nonSchmid_neg ( 1 : 3 , 1 : 3 , s ) )
2019-03-17 21:32:08 +05:30
else
2020-03-16 23:13:04 +05:30
tauNS ( s , 3 ) = math_tensordot ( Mp , + prm % nonSchmid_neg ( 1 : 3 , 1 : 3 , s ) )
tauNS ( s , 4 ) = math_tensordot ( Mp , - prm % nonSchmid_pos ( 1 : 3 , 1 : 3 , s ) )
2019-03-17 21:32:08 +05:30
endif
enddo
2021-01-26 13:09:17 +05:30
tauNS = tauNS + spread ( dst % tau_back ( : , me ) , 2 , 4 )
tau = tau + dst % tau_back ( : , me )
2020-02-07 16:14:03 +05:30
! edges
2020-03-16 19:52:44 +05:30
call kinetics ( v ( : , 1 ) , dv_dtau ( : , 1 ) , dv_dtauNS ( : , 1 ) , &
2021-02-14 19:59:10 +05:30
tau , tauNS ( : , 1 ) , dst % tau_pass ( : , me ) , 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 )
2020-02-07 16:14:03 +05:30
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-02-14 19:59:10 +05:30
tau , tauNS ( : , t ) , dst % tau_pass ( : , me ) , 2 , Temperature , ph )
2019-03-17 21:32:08 +05:30
enddo
2021-03-26 18:21:09 +05:30
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 )
2014-06-14 02:23:17 +05:30
endif
2019-12-01 13:25:24 +05:30
2021-01-26 13:09:17 +05:30
stt % v ( : , me ) = pack ( v , . true . )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!*** Bauschinger effect
forall ( s = 1 : ns , t = 5 : 8 , rhoSgl ( s , t ) * v ( s , t - 4 ) < 0.0_pReal ) &
rhoSgl ( s , t - 4 ) = rhoSgl ( s , t - 4 ) + abs ( rhoSgl ( s , t ) )
2020-02-07 16:14:03 +05:30
2020-09-22 21:41:09 +05:30
gdotTotal = sum ( rhoSgl ( : , 1 : 4 ) * v , 2 ) * prm % b_sl
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
do s = 1 , ns
Lp = Lp + gdotTotal ( s ) * prm % Schmid ( 1 : 3 , 1 : 3 , s )
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 ) &
+ prm % Schmid ( i , j , s ) * prm % Schmid ( 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 ) &
2019-03-17 21:32:08 +05:30
+ prm % Schmid ( i , j , s ) &
2021-03-26 18:21:09 +05:30
* ( + prm % nonSchmid_pos ( k , l , s ) * rhoSgl ( s , 3 ) * dv_dtauNS ( s , 3 ) &
- prm % nonSchmid_neg ( k , l , s ) * rhoSgl ( s , 4 ) * dv_dtauNS ( s , 4 ) ) * prm % b_sl ( s )
2014-06-14 02:23:17 +05:30
enddo
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
end associate
2020-02-07 16:14:03 +05:30
2021-01-26 12:03:04 +05:30
end subroutine nonlocal_LpAndItsTangent
2014-06-14 02:23:17 +05:30
2014-06-26 19:23:12 +05:30
2014-06-14 02:23:17 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
2021-02-14 21:59:23 +05:30
module subroutine plastic_nonlocal_deltaState ( Mp , ph , me )
2020-02-07 16:14:03 +05:30
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 ) :: &
2021-02-14 19:59:10 +05:30
ph , &
2021-02-14 21:59:23 +05:30
me
2014-07-02 17:57:39 +05:30
2019-03-17 21:32:08 +05:30
integer :: &
ns , & ! short notation for the total number of active slip systems
c , & ! character of dislocation
t , & ! type of dislocation
2020-07-02 00:52:05 +05:30
s ! index of my current slip system
2021-02-14 23:00:57 +05:30
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)
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 10 ) :: &
2019-03-17 21:32:08 +05:30
rho ! current dislocation densities
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 4 ) :: &
2019-03-17 21:32:08 +05:30
v ! dislocation glide velocity
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl ) :: &
2019-03-17 21:32:08 +05:30
tau ! current resolved shear stress
2021-02-14 23:00:57 +05:30
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
2020-02-07 16:14:03 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) , dst = > microstructure ( ph ) , del = > deltaState ( ph ) )
2021-02-14 19:59:10 +05:30
ns = prm % sum_N_sl
2020-02-07 16:14:03 +05:30
!*** shortcut to state variables
2021-02-14 23:00:57 +05:30
forall ( s = 1 : ns , t = 1 : 4 ) v ( s , t ) = plasticState ( ph ) % state ( iV ( s , t , ph ) , me )
forall ( s = 1 : ns , c = 1 : 2 ) dUpperOld ( s , c ) = plasticState ( ph ) % state ( iD ( s , c , ph ) , me )
2020-02-07 16:14:03 +05:30
2021-02-14 21:59:23 +05:30
rho = getRho ( ph , me )
2019-03-17 21:32:08 +05:30
rhoDip = rho ( : , dip )
2020-02-07 16:14:03 +05:30
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
2020-02-07 16:14:03 +05:30
endwhere
2019-03-17 21:32:08 +05:30
deltaRhoRemobilization ( : , dip ) = 0.0_pReal
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** calculate dipole formation and dissociation by stress change
2020-02-07 16:14:03 +05:30
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-01-26 13:09:17 +05:30
tau ( s ) = math_tensordot ( Mp , prm % Schmid ( 1 : 3 , 1 : 3 , s ) ) + dst % tau_back ( s , me )
2019-03-17 21:32:08 +05:30
if ( abs ( tau ( s ) ) < 1.0e-15_pReal ) tau ( s ) = 1.0e-15_pReal
enddo
2020-02-07 16:14:03 +05:30
2020-09-22 21:41:09 +05:30
dUpper ( : , 1 ) = prm % mu * prm % b_sl / ( 8.0_pReal * PI * ( 1.0_pReal - prm % nu ) * abs ( tau ) )
dUpper ( : , 2 ) = prm % mu * prm % b_sl / ( 4.0_pReal * PI * abs ( tau ) )
2020-02-07 16:14:03 +05:30
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 ) )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
dUpper = max ( dUpper , prm % minDipoleHeight )
deltaDUpper = dUpper - dUpperOld
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal
forall ( c = 1 : 2 , s = 1 : ns , deltaDUpper ( s , c ) < 0.0_pReal . and . &
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-02-07 16:14:03 +05:30
2020-03-16 13:06:19 +05:30
forall ( t = 1 : 4 ) deltaRhoDipole2SingleStress ( : , t ) = - 0.5_pReal * deltaRhoDipole2SingleStress ( : , ( t - 1 ) / 2 + 9 )
2021-02-14 23:00:57 +05:30
forall ( s = 1 : ns , c = 1 : 2 ) plasticState ( ph ) % state ( iD ( s , c , ph ) , me ) = dUpper ( s , c )
2020-02-07 16:14:03 +05:30
2021-01-26 13:09:17 +05:30
plasticState ( ph ) % deltaState ( : , me ) = 0.0_pReal
del % rho ( : , me ) = reshape ( deltaRhoRemobilization + deltaRhoDipole2SingleStress , [ 10 * ns ] )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
end associate
2014-06-14 02:23:17 +05:30
2014-12-08 21:25:30 +05:30
end subroutine plastic_nonlocal_deltaState
2014-06-14 02:23:17 +05:30
2019-02-20 05:11:44 +05:30
2014-06-26 19:23:12 +05:30
!---------------------------------------------------------------------------------------------------
2014-06-14 02:23:17 +05:30
!> @brief calculates the rate of change of microstructure
2014-06-26 19:23:12 +05:30
!---------------------------------------------------------------------------------------------------
2021-01-26 16:47:00 +05:30
module subroutine nonlocal_dotState ( Mp , Temperature , timestep , &
2021-01-27 03:11:41 +05:30
ph , me , ip , el )
2020-02-07 16:14:03 +05:30
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 ) :: &
2021-01-27 03:11:41 +05:30
ph , &
2021-01-26 13:09:17 +05:30
me , &
2020-03-16 14:19:59 +05:30
ip , & !< current integration point
el !< current element number
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
integer :: &
ns , & !< short notation for the total number of active slip systems
c , & !< character of dislocation
t , & !< type of dislocation
2020-07-02 00:52:05 +05:30
s !< index of my current slip system
2021-02-14 23:00:57 +05:30
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
2021-02-14 23:00:57 +05:30
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)
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 4 ) :: &
2019-03-17 21:32:08 +05:30
v , & !< current dislocation glide velocity
2020-02-07 16:23:50 +05:30
v0 , &
2019-03-17 21:32:08 +05:30
gdot !< shear rates
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl ) :: &
2019-03-17 21:32:08 +05:30
tau , & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles
2021-02-14 23:00:57 +05:30
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 ) :: &
selfDiffusion !< self diffusion
2020-02-07 16:14:03 +05:30
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
2011-03-29 12:57:19 +05:30
endif
2020-02-07 16:14:03 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) , &
dst = > microstructure ( ph ) , &
dot = > dotState ( ph ) , &
stt = > state ( ph ) )
2020-03-16 23:36:10 +05:30
ns = prm % sum_N_sl
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
tau = 0.0_pReal
gdot = 0.0_pReal
2020-02-07 16:14:03 +05:30
2021-02-14 21:59:23 +05:30
rho = getRho ( ph , me )
2019-03-17 21:32:08 +05:30
rhoSgl = rho ( : , sgl )
rhoDip = rho ( : , dip )
2021-02-14 21:59:23 +05:30
rho0 = getRho0 ( ph , me )
2020-02-11 10:11:10 +05:30
my_rhoSgl0 = rho0 ( : , sgl )
2020-02-07 16:14:03 +05:30
2021-02-14 23:00:57 +05:30
forall ( s = 1 : ns , t = 1 : 4 ) v ( s , t ) = plasticState ( ph ) % state ( iV ( s , t , ph ) , me )
2020-09-22 21:41:09 +05:30
gdot = rhoSgl ( : , 1 : 4 ) * v * spread ( prm % b_sl , 2 , 4 )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
#ifdef DEBUG
2020-07-02 04:55:24 +05:30
if ( debugConstitutive % basic &
. and . ( ( debugConstitutive % element == el . and . debugConstitutive % ip == ip ) &
. or . . not . debugConstitutive % selective ) ) then
2020-09-19 14:31:43 +05:30
print '(a,/,10(12x,12(e12.5,1x),/))' , '<< CONST >> rho / 1/m^2' , rhoSgl , rhoDip
print '(a,/,4(12x,12(e12.5,1x),/))' , '<< CONST >> gdot / 1/s' , gdot
2019-03-17 21:32:08 +05:30
endif
#endif
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!****************************************************************************
2020-03-16 14:19:59 +05:30
!*** limits for stable dipole height
2020-03-16 13:06:19 +05:30
do s = 1 , ns
2021-01-26 13:09:17 +05:30
tau ( s ) = math_tensordot ( Mp , prm % Schmid ( 1 : 3 , 1 : 3 , s ) ) + dst % tau_back ( s , me )
2019-03-17 21:32:08 +05:30
if ( abs ( tau ( s ) ) < 1.0e-15_pReal ) tau ( s ) = 1.0e-15_pReal
enddo
2020-02-07 16:14:03 +05:30
2019-03-17 22:29:01 +05:30
dLower = prm % minDipoleHeight
2020-09-22 21:41:09 +05:30
dUpper ( : , 1 ) = prm % mu * prm % b_sl / ( 8.0_pReal * PI * ( 1.0_pReal - prm % nu ) * abs ( tau ) )
dUpper ( : , 2 ) = prm % mu * prm % b_sl / ( 4.0_pReal * PI * abs ( tau ) )
2020-02-07 16:14:03 +05:30
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 )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!****************************************************************************
2020-03-16 14:19:59 +05:30
!*** dislocation multiplication
2019-03-17 21:32:08 +05:30
rhoDotMultiplication = 0.0_pReal
isBCC : if ( lattice_structure ( ph ) == LATTICE_bcc_ID ) then
forall ( s = 1 : ns , sum ( abs ( v ( s , 1 : 4 ) ) ) > 0.0_pReal )
2020-09-22 21:41:09 +05:30
rhoDotMultiplication ( s , 1 : 2 ) = sum ( abs ( gdot ( s , 3 : 4 ) ) ) / prm % b_sl ( s ) & ! assuming double-cross-slip of screws to be decisive for multiplication
2021-01-26 13:09:17 +05:30
* sqrt ( stt % rho_forest ( s , me ) ) / 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
2020-09-22 21:41:09 +05:30
rhoDotMultiplication ( s , 3 : 4 ) = sum ( abs ( gdot ( s , 3 : 4 ) ) ) / prm % b_sl ( s ) & ! assuming double-cross-slip of screws to be decisive for multiplication
2021-01-26 13:09:17 +05:30
* sqrt ( stt % rho_forest ( s , me ) ) / 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
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
else isBCC
2020-03-16 13:06:19 +05:30
rhoDotMultiplication ( : , 1 : 4 ) = spread ( &
2020-09-22 21:41:09 +05:30
( sum ( abs ( gdot ( : , 1 : 2 ) ) , 2 ) * prm % f_ed_mult + sum ( abs ( gdot ( : , 3 : 4 ) ) , 2 ) ) &
2021-01-26 13:09:17 +05:30
* sqrt ( stt % rho_forest ( : , me ) ) / prm % i_sl / prm % b_sl , 2 , 4 )
2019-03-17 21:32:08 +05:30
endif isBCC
2020-02-07 16:14:03 +05:30
2021-02-14 23:00:57 +05:30
forall ( s = 1 : ns , t = 1 : 4 ) v0 ( s , t ) = plasticState ( ph ) % state0 ( iV ( s , t , ph ) , me )
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** calculate dipole formation and annihilation
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!*** formation by glide
do c = 1 , 2
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide ( : , 2 * c - 1 ) = - 2.0_pReal * dUpper ( : , c ) / prm % b_sl &
2020-03-16 13:06:19 +05:30
* ( rhoSgl ( : , 2 * c - 1 ) * abs ( gdot ( : , 2 * c ) ) & ! negative mobile --> positive mobile
+ rhoSgl ( : , 2 * c ) * abs ( gdot ( : , 2 * c - 1 ) ) & ! positive mobile --> negative mobile
+ abs ( rhoSgl ( : , 2 * c + 4 ) ) * abs ( gdot ( : , 2 * c - 1 ) ) ) ! positive mobile --> negative immobile
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide ( : , 2 * c ) = - 2.0_pReal * dUpper ( : , c ) / prm % b_sl &
2020-03-16 13:06:19 +05:30
* ( rhoSgl ( : , 2 * c - 1 ) * abs ( gdot ( : , 2 * c ) ) & ! negative mobile --> positive mobile
+ rhoSgl ( : , 2 * c ) * abs ( gdot ( : , 2 * c - 1 ) ) & ! positive mobile --> negative mobile
+ abs ( rhoSgl ( : , 2 * c + 3 ) ) * abs ( gdot ( : , 2 * c ) ) ) ! negative mobile --> positive immobile
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide ( : , 2 * c + 3 ) = - 2.0_pReal * dUpper ( : , c ) / prm % b_sl &
2020-03-16 13:06:19 +05:30
* rhoSgl ( : , 2 * c + 3 ) * abs ( gdot ( : , 2 * c ) ) ! negative mobile --> positive immobile
2020-09-22 21:41:09 +05:30
rhoDotSingle2DipoleGlide ( : , 2 * c + 4 ) = - 2.0_pReal * dUpper ( : , c ) / prm % b_sl &
2020-03-16 13:06:19 +05:30
* rhoSgl ( : , 2 * c + 4 ) * abs ( gdot ( : , 2 * c - 1 ) ) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide ( : , c + 8 ) = abs ( rhoDotSingle2DipoleGlide ( : , 2 * c + 3 ) ) &
+ abs ( rhoDotSingle2DipoleGlide ( : , 2 * c + 4 ) ) &
- rhoDotSingle2DipoleGlide ( : , 2 * c - 1 ) &
- rhoDotSingle2DipoleGlide ( : , 2 * c )
2019-03-17 21:32:08 +05:30
enddo
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
!*** athermal annihilation
rhoDotAthermalAnnihilation = 0.0_pReal
2020-02-07 16:14:03 +05:30
forall ( c = 1 : 2 ) &
2020-09-22 21:41:09 +05:30
rhoDotAthermalAnnihilation ( : , c + 8 ) = - 2.0_pReal * dLower ( : , c ) / prm % b_sl &
* ( 2.0_pReal * ( rhoSgl ( : , 2 * c - 1 ) * abs ( gdot ( : , 2 * c ) ) + rhoSgl ( : , 2 * c ) * abs ( gdot ( : , 2 * c - 1 ) ) ) & ! was single hitting single
2020-03-16 13:06:19 +05:30
+ 2.0_pReal * ( abs ( rhoSgl ( : , 2 * c + 3 ) ) * abs ( gdot ( : , 2 * c ) ) + abs ( rhoSgl ( : , 2 * c + 4 ) ) * abs ( gdot ( : , 2 * c - 1 ) ) ) & ! was single hitting immobile single or was immobile single hit by single
2020-09-22 21:41:09 +05:30
+ rhoDip ( : , c ) * ( abs ( gdot ( : , 2 * c - 1 ) ) + abs ( gdot ( : , 2 * c ) ) ) ) ! single knocks dipole constituent
2020-02-07 16:14:03 +05:30
2020-03-16 14:19:59 +05:30
! annihilated screw dipoles leave edge jogs behind on the colinear system
2019-03-17 21:32:08 +05:30
if ( lattice_structure ( ph ) == LATTICE_fcc_ID ) &
forall ( s = 1 : ns , prm % colinearSystem ( s ) > 0 ) &
rhoDotAthermalAnnihilation ( prm % colinearSystem ( s ) , 1 : 2 ) = - rhoDotAthermalAnnihilation ( s , 10 ) &
2021-01-26 13:09:17 +05:30
* 0.25_pReal * sqrt ( stt % rho_forest ( s , me ) ) * ( dUpper ( s , 2 ) + dLower ( s , 2 ) ) * prm % f_ed
2020-02-07 16:14:03 +05:30
2021-01-26 13:09:17 +05:30
!*** thermally activated annihilation me edge dipoles by climb
2019-03-17 21:32:08 +05:30
rhoDotThermalAnnihilation = 0.0_pReal
2020-09-22 21:41:09 +05:30
selfDiffusion = prm % D_0 * exp ( - prm % Q_cl / ( kB * Temperature ) )
vClimb = prm % V_at * selfDiffusion * prm % mu &
2020-03-16 19:52:44 +05:30
/ ( kB * Temperature * PI * ( 1.0_pReal - prm % nu ) * ( dUpper ( : , 1 ) + dLower ( : , 1 ) ) )
2019-03-17 21:32:08 +05:30
forall ( s = 1 : ns , dUpper ( s , 1 ) > dLower ( s , 1 ) ) &
rhoDotThermalAnnihilation ( s , 9 ) = max ( - 4.0_pReal * rhoDip ( s , 1 ) * vClimb ( s ) / ( dUpper ( s , 1 ) - dLower ( s , 1 ) ) , &
- 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-02-14 19:59:10 +05:30
rhoDot = rhoDotFlux ( timestep , ph , me , ip , el ) &
2019-03-17 21:32:08 +05:30
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
2020-02-07 16:14:03 +05:30
+ rhoDotThermalAnnihilation
2012-08-16 16:33:22 +05:30
2012-08-23 11:18:21 +05:30
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
2017-10-03 18:50:53 +05:30
#ifdef DEBUG
2020-07-02 04:55:24 +05:30
if ( debugConstitutive % extensive ) then
2020-09-19 14:31:43 +05:30
print '(a,i5,a,i2)' , '<< CONST >> evolution rate leads to negative density at el ' , el , ' ip ' , ip
print '(a)' , '<< CONST >> enforcing cutback !!!'
2019-03-17 21:32:08 +05:30
endif
2012-08-23 11:18:21 +05:30
#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-01-26 13:09:17 +05:30
dot % rho ( : , me ) = pack ( rhoDot , . true . )
dot % gamma ( : , me ) = sum ( gdot , 2 )
2019-03-17 21:32:08 +05:30
endif
2014-07-02 17:57:39 +05:30
2019-03-17 21:32:08 +05:30
end associate
2014-06-14 02:23:17 +05:30
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
2020-05-03 14:19:36 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
2021-02-14 19:59:10 +05:30
function rhoDotFlux ( timestep , ph , me , ip , el )
2020-05-03 14:19:36 +05:30
real ( pReal ) , intent ( in ) :: &
timestep !< substepped crystallite time increment
integer , intent ( in ) :: &
2021-02-14 19:59:10 +05:30
ph , &
2021-01-26 13:09:17 +05:30
me , &
2020-05-03 14:19:36 +05:30
ip , & !< current integration point
el !< current element number
2020-07-03 20:15:11 +05:30
2020-05-03 14:19:36 +05:30
integer :: &
2021-02-14 23:00:57 +05:30
neighbor_ph , & !< phase of my neighbor's plasticity
2020-05-03 14:19:36 +05:30
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
neighbor_n , & !< neighbor index pointing to me 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
opposite_n , & !< neighbor index pointing to me 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
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 10 ) :: &
2020-05-03 14:19:36 +05:30
rho , &
rho0 , & !< dislocation density at beginning of time step
2020-05-03 14:24:57 +05:30
rhoDotFlux !< density evolution by flux
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 8 ) :: &
2020-05-03 14:19:36 +05:30
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)
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 4 ) :: &
2020-05-03 14:19:36 +05:30
v , & !< current dislocation glide velocity
v0 , &
neighbor_v0 , & !< dislocation glide velocity of enighboring ip
gdot !< shear rates
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( 3 , param ( ph ) % sum_N_sl , 4 ) :: &
2020-05-03 14:19:36 +05:30
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
Favg !< average total deformation gradient of me and my neighbor
real ( pReal ) , dimension ( 3 ) :: &
normal_neighbor2me , & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration
normal_neighbor2me_defConf , & !< interface normal pointing from my neighbor to me in shared deformed configuration
normal_me2neighbor , & !< interface normal pointing from me to my neighbor in my lattice configuration
normal_me2neighbor_defConf !< interface normal pointing from me 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
2020-05-03 14:19:36 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) , &
dst = > microstructure ( ph ) , &
dot = > dotState ( ph ) , &
stt = > state ( ph ) )
2020-05-03 14:19:36 +05:30
ns = prm % sum_N_sl
gdot = 0.0_pReal
2021-02-14 21:59:23 +05:30
rho = getRho ( ph , me )
2020-05-03 14:19:36 +05:30
rhoSgl = rho ( : , sgl )
2021-02-14 21:59:23 +05:30
rho0 = getRho0 ( ph , me )
2020-05-03 14:19:36 +05:30
my_rhoSgl0 = rho0 ( : , sgl )
2021-02-14 23:00:57 +05:30
forall ( s = 1 : ns , t = 1 : 4 ) v ( s , t ) = plasticState ( ph ) % state ( iV ( s , t , ph ) , me ) !ToDo: MD: I think we should use state0 here
2020-09-22 21:41:09 +05:30
gdot = rhoSgl ( : , 1 : 4 ) * v * spread ( prm % b_sl , 2 , 4 )
2020-05-03 14:19:36 +05:30
2021-02-14 23:00:57 +05:30
forall ( s = 1 : ns , t = 1 : 4 ) v0 ( s , t ) = plasticState ( ph ) % state0 ( iV ( s , t , ph ) , me )
2020-05-03 14:19:36 +05:30
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity)
rhoDotFlux = 0.0_pReal
if ( . not . phase_localPlasticity ( material_phaseAt ( 1 , el ) ) ) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if ( any ( abs ( gdot ) > 0.0_pReal & ! any active slip system ...
2020-09-22 21:41:09 +05:30
. and . prm % f_c * abs ( v0 ) * timestep &
2020-05-03 14:19:36 +05:30
> 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
2020-07-02 04:55:24 +05:30
if ( debugConstitutive % extensive ) then
2020-09-19 14:31:43 +05:30
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 ' , &
2020-05-03 14:19:36 +05:30
maxval ( abs ( v0 ) , abs ( gdot ) > 0.0_pReal &
2020-09-22 21:41:09 +05:30
. and . prm % f_c * abs ( v0 ) * timestep &
2020-05-03 14:19:36 +05:30
> IPvolume ( ip , el ) / maxval ( IParea ( : , ip , el ) ) ) , &
' at a timestep of ' , timestep
2020-09-17 22:58:41 +05:30
print * , '<< CONST >> enforcing cutback !!!'
2020-05-03 14:19:36 +05:30
endif
#endif
2020-05-03 15:34:57 +05:30
rhoDotFlux = IEEE_value ( 1.0_pReal , IEEE_quiet_NaN ) ! enforce cutback
2020-05-03 14:19:36 +05:30
return
endif
!*** 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-02-09 03:51:53 +05:30
my_F = phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me )
my_Fe = matmul ( my_F , math_inv33 ( phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ) )
2020-05-03 14:19:36 +05:30
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
2021-02-14 23:00:57 +05:30
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 ) ) )
2020-05-03 14:19:36 +05:30
Favg = 0.5_pReal * ( my_F + neighbor_F )
else ! if no neighbor, take my value as average
Favg = my_F
endif
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 )
2021-02-14 23:00:57 +05:30
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 )
2020-05-03 14:19:36 +05:30
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 ) &
2020-05-03 14:19:36 +05:30
neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33 ( Favg ) * matmul ( math_inv33 ( transpose ( Favg ) ) , &
IPareaNormal ( 1 : 3 , neighbor_n , neighbor_ip , neighbor_el ) ) ! normal of the interface in (average) deformed configuration (pointing neighbor => me)
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 )
if ( neighbor_v0 ( s , t ) * math_inner ( m ( 1 : 3 , s , t ) , normal_neighbor2me ) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me
. 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
endif
enddo
enddo
endif ; endif
!* 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).
!* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
!* 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 ) &
* matmul ( math_inv33 ( transpose ( Favg ) ) , IPareaNormal ( 1 : 3 , n , ip , el ) ) ! normal of the interface in (average) deformed configuration (pointing me => 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
if ( v0 ( s , t ) * math_inner ( m ( 1 : 3 , s , t ) , normal_me2neighbor ) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (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
endif
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
endif
enddo
enddo
endif ; endif
enddo neighbors
endif
end associate
2020-05-03 16:39:49 +05:30
end function rhoDotFlux
2020-05-03 14:19:36 +05:30
2019-03-17 21:32:08 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-17 22:29:01 +05:30
!> @brief Compatibility update
2020-02-07 16:14:03 +05:30
!> @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.
!--------------------------------------------------------------------------------------------------
2021-02-14 19:59:10 +05:30
module subroutine plastic_nonlocal_updateCompatibility ( orientation , ph , i , e )
2020-02-07 16:14:03 +05:30
2020-10-28 01:57:26 +05:30
type ( rotation ) , dimension ( 1 , discretization_nIPs , discretization_Nelems ) , intent ( in ) :: &
2020-03-16 14:19:59 +05:30
orientation ! crystal orientation
2019-03-17 21:32:08 +05:30
integer , intent ( in ) :: &
2021-02-14 19:59:10 +05:30
ph , &
2019-03-17 21:32:08 +05:30
i , &
e
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
integer :: &
2020-02-07 16:14:03 +05:30
n , & ! neighbor index
2021-02-22 23:02:31 +05:30
me , &
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
s1 , & ! slip system index (me)
s2 ! slip system index (my neighbor)
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( 2 , param ( ph ) % sum_N_sl , param ( ph ) % sum_N_sl , nIPneighbors ) :: &
2020-02-07 16:14:03 +05:30
my_compatibility ! my_compatibility for current element and ip
2019-03-17 21:32:08 +05:30
real ( pReal ) :: &
my_compatibilitySum , &
thresholdValue , &
nThresholdValues
2021-02-14 23:00:57 +05:30
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
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) )
2020-03-16 23:36:10 +05:30
ns = prm % sum_N_sl
2020-02-07 16:14:03 +05:30
2021-02-22 23:02:31 +05:30
me = 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-02-07 16:14:03 +05:30
2020-03-16 13:06:19 +05:30
neighbors : do n = 1 , nIPneighbors
2019-06-06 11:49:29 +05:30
neighbor_e = IPneighborhood ( 1 , n , i , e )
neighbor_i = IPneighborhood ( 2 , n , i , e )
2021-02-22 23:02:31 +05:30
neighbor_me = material_phaseMemberAt ( 1 , neighbor_i , neighbor_e )
2020-03-16 14:07:25 +05:30
neighbor_phase = material_phaseAt ( 1 , neighbor_e )
2020-02-07 16:14:03 +05:30
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
!* Set surface transmissivity to the value specified in the material.config
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 we encounter a different nonlocal phase at the neighbor,
!* we consider this to be a real "physical" phase boundary, so completely incompatible.
!* If one of the two phases has a local plasticity law,
!* we do not consider this to be a phase boundary, so completely compatible.
if ( . not . phase_localPlasticity ( neighbor_phase ) . and . . not . phase_localPlasticity ( ph ) ) &
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
2020-03-16 14:07:25 +05:30
!* GRAIN BOUNDARY !
!* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
2021-02-22 23:02:31 +05:30
if ( any ( dNeq ( material_orientation0 ( 1 , ph , me ) % asQuaternion ( ) , &
material_orientation0 ( 1 , neighbor_phase , neighbor_me ) % asQuaternion ( ) ) ) . and . &
2020-03-16 14:12:58 +05:30
( . not . phase_localPlasticity ( neighbor_phase ) ) ) &
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.
2019-12-02 17:28:23 +05:30
mis = orientation ( 1 , i , e ) % misorientation ( orientation ( 1 , neighbor_i , neighbor_e ) )
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 ) ) ) )
2019-04-03 11:52:04 +05:30
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 ) ) ) )
2019-03-17 21:32:08 +05:30
enddo neighborSlipSystems
2020-02-07 16:14:03 +05:30
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
enddo
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
2019-03-17 21:32:08 +05:30
enddo mySlipSystems
endif
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
enddo neighbors
2020-02-07 16:14:03 +05:30
2020-03-16 13:06:19 +05:30
compatibility ( : , : , : , : , i , e ) = my_compatibility
2020-02-07 16:14:03 +05:30
2019-03-17 21:32:08 +05:30
end associate
2020-02-07 16:14:03 +05:30
2014-12-08 21:25:30 +05:30
end subroutine plastic_nonlocal_updateCompatibility
2014-06-14 02:23:17 +05:30
2020-02-07 16:14:03 +05:30
2019-03-10 01:13:31 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
2021-02-14 19:59:10 +05:30
module subroutine plastic_nonlocal_results ( ph , group )
2019-03-10 01:13:31 +05:30
2021-02-14 19:59:10 +05:30
integer , intent ( in ) :: ph
2019-12-05 01:20:46 +05:30
character ( len = * ) , intent ( in ) :: group
2020-02-14 13:56:26 +05:30
2019-03-10 01:13:31 +05:30
integer :: o
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) , dst = > microstructure ( ph ) , stt = > state ( ph ) )
2020-02-14 13:56:26 +05:30
outputsLoop : do o = 1 , size ( prm % output )
select case ( trim ( prm % output ( o ) ) )
2020-08-15 19:32:10 +05:30
case ( 'rho_u_ed_pos' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_mob_edg_pos , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'positive mobile edge density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_b_ed_pos' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_imm_edg_pos , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'positive immobile edge density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_u_ed_neg' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_mob_edg_neg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'negative mobile edge density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_b_ed_neg' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_imm_edg_neg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'negative immobile edge density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_d_ed' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_dip_edg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'edge dipole density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_u_sc_pos' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_mob_scr_pos , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'positive mobile screw density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_b_sc_pos' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_imm_scr_pos , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'positive immobile screw density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_u_sc_neg' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_mob_scr_neg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'negative mobile screw density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_b_sc_neg' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_sgl_imm_scr_neg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'negative immobile screw density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_d_sc' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_dip_scr , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'screw dipole density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'rho_f' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % rho_forest , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'forest density' , '1/m²' )
2020-08-15 19:32:10 +05:30
case ( 'v_ed_pos' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % v_edg_pos , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'positive edge velocity' , 'm/s' )
2020-08-15 19:32:10 +05:30
case ( 'v_ed_neg' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % v_edg_neg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'negative edge velocity' , 'm/s' )
2020-08-15 19:32:10 +05:30
case ( 'v_sc_pos' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % v_scr_pos , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'positive srew velocity' , 'm/s' )
2020-08-15 19:32:10 +05:30
case ( 'v_sc_neg' )
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % v_scr_neg , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'negative screw velocity' , 'm/s' )
2020-02-14 13:56:26 +05:30
case ( 'gamma' )
2020-08-15 19:32:10 +05:30
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , stt % gamma , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'plastic shear' , '1' )
2020-02-14 13:56:26 +05:30
case ( 'tau_pass' )
2020-08-15 19:32:10 +05:30
if ( prm % sum_N_sl > 0 ) call results_writeDataset ( group , dst % tau_pass , trim ( prm % output ( o ) ) , &
2020-03-16 23:36:10 +05:30
'passing stress for slip' , 'Pa' )
2019-03-10 01:13:31 +05:30
end select
enddo outputsLoop
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-03-05 01:46:36 +05:30
subroutine stateInit ( ini , phase , Nmembers )
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-03-05 01:46:36 +05:30
Nmembers
2020-03-17 11:47:40 +05:30
integer :: &
i , &
2021-02-14 23:00:57 +05:30
e , &
2020-03-17 11:47:40 +05:30
f , &
from , &
upto , &
s , &
phasemember
real ( pReal ) , dimension ( 2 ) :: &
noise , &
rnd
real ( pReal ) :: &
meanDensity , &
totalVolume , &
densityBinning , &
minimumIpVolume
2021-03-05 01:46:36 +05:30
real ( pReal ) , dimension ( Nmembers ) :: &
2020-03-17 11:47:40 +05:30
volume
2020-08-15 19:32:10 +05:30
2021-02-14 23:00:57 +05:30
associate ( stt = > state ( phase ) )
2020-03-17 11:47:40 +05:30
2020-09-22 21:41:09 +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
2020-10-28 01:57:26 +05:30
do e = 1 , discretization_Nelems
do i = 1 , discretization_nIPs
2020-03-17 11:47:40 +05:30
if ( material_phaseAt ( 1 , e ) == phase ) volume ( material_phasememberAt ( 1 , i , e ) ) = IPvolume ( i , e )
enddo
enddo
2020-03-17 12:05:25 +05:30
totalVolume = sum ( volume )
2020-03-17 11:47:40 +05:30
minimumIPVolume = minval ( volume )
2020-09-22 21:41:09 +05:30
densityBinning = ini % random_rho_u_binning / minimumIpVolume ** ( 2.0_pReal / 3.0_pReal )
2020-03-17 11:47:40 +05:30
! fill random material points with dislocation segments until the desired overall density is reached
meanDensity = 0.0_pReal
2020-09-22 21:41:09 +05:30
do while ( meanDensity < ini % random_rho_u )
2020-03-17 11:47:40 +05:30
call random_number ( rnd )
2021-03-05 01:46:36 +05:30
phasemember = nint ( rnd ( 1 ) * real ( Nmembers , pReal ) + 0.5_pReal )
2020-03-17 12:05:25 +05:30
s = nint ( rnd ( 2 ) * real ( sum ( ini % N_sl ) , pReal ) * 4.0_pReal + 0.5_pReal )
2020-03-17 11:47:40 +05:30
meanDensity = meanDensity + densityBinning * volume ( phasemember ) / totalVolume
stt % rhoSglMobile ( s , phasemember ) = densityBinning
enddo
else ! homogeneous distribution with noise
2021-03-05 01:46:36 +05:30
do e = 1 , Nmembers
2020-03-17 12:05:25 +05:30
do f = 1 , size ( ini % N_sl , 1 )
from = 1 + sum ( ini % N_sl ( 1 : f - 1 ) )
upto = sum ( ini % N_sl ( 1 : f ) )
2020-03-17 11:47:40 +05:30
do s = from , upto
2020-09-22 21:41:09 +05:30
noise = [ math_sampleGaussVar ( 0.0_pReal , ini % sigma_rho_u ) , &
math_sampleGaussVar ( 0.0_pReal , ini % sigma_rho_u ) ]
2021-02-14 23:00:57 +05:30
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 )
2020-09-22 21:41:09 +05:30
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 )
2020-03-17 11:47:40 +05:30
enddo
2020-09-22 21:41:09 +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 )
2020-03-17 11:47:40 +05:30
enddo
enddo
endif
end associate
end subroutine stateInit
2020-03-16 19:52:44 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics
!--------------------------------------------------------------------------------------------------
2021-02-14 19:59:10 +05:30
pure subroutine kinetics ( v , dv_dtau , dv_dtauNS , tau , tauNS , tauThreshold , c , Temperature , 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)
2021-02-14 19:59:10 +05:30
ph
2020-03-16 19:52:44 +05:30
real ( pReal ) , intent ( in ) :: &
Temperature !< temperature
2021-02-14 23:00:57 +05:30
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
2021-02-14 23:00:57 +05:30
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
vViscous , & !< viscous glide velocity
dtPeierls_dtau , & !< derivative with respect to resolved shear stress
dtSolidSolution_dtau , & !< derivative with respect to resolved shear stress
meanfreepath_S , & !< mean free travel distance for dislocations between two solid solution obstacles
meanfreepath_P , & !< mean free travel distance for dislocations between two Peierls barriers
jumpWidth_P , & !< depth of activated area
jumpWidth_S , & !< depth of activated area
activationLength_P , & !< length of activated dislocation line
activationLength_S , & !< length of activated dislocation line
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
activationEnergy_S , & !< energy that is needed to overcome barrier
criticalStress_P , & !< maximum obstacle strength
criticalStress_S , & !< maximum obstacle strength
mobility !< dislocation mobility
2021-02-14 23:00:57 +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
!* Effective stress includes non Schmid constributions
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = max ( 0.0_pReal , abs ( tauNS ( s ) ) - tauThreshold ( s ) ) ! ensure that the effective stress is positive
2020-09-22 21:41:09 +05:30
meanfreepath_P = prm % b_sl ( s )
jumpWidth_P = prm % b_sl ( s )
activationLength_P = prm % w * prm % b_sl ( s )
activationVolume_P = activationLength_P * jumpWidth_P * prm % b_sl ( s )
2020-03-16 19:52:44 +05:30
criticalStress_P = prm % peierlsStress ( s , c )
activationEnergy_P = criticalStress_P * activationVolume_P
tauRel_P = min ( 1.0_pReal , tauEff / criticalStress_P ) ! ensure that the activation probability cannot become greater than one
2020-09-22 21:41:09 +05:30
tPeierls = 1.0_pReal / prm % nu_a &
2020-03-16 19:52:44 +05:30
* exp ( activationEnergy_P / ( kB * Temperature ) &
* ( 1.0_pReal - tauRel_P ** prm % p ) ** prm % q )
if ( tauEff < criticalStress_P ) then
dtPeierls_dtau = tPeierls * prm % p * prm % q * activationVolume_P / ( kB * Temperature ) &
* ( 1.0_pReal - tauRel_P ** prm % p ) ** ( prm % q - 1.0_pReal ) * tauRel_P ** ( prm % p - 1.0_pReal )
else
dtPeierls_dtau = 0.0_pReal
endif
!* Contribution from solid solution strengthening
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = abs ( tau ( s ) ) - tauThreshold ( s )
2020-09-22 21:41:09 +05:30
meanfreepath_S = prm % b_sl ( s ) / sqrt ( prm % c_sol )
jumpWidth_S = prm % f_sol * prm % b_sl ( s )
activationLength_S = prm % b_sl ( s ) / sqrt ( prm % c_sol )
activationVolume_S = activationLength_S * jumpWidth_S * prm % b_sl ( s )
activationEnergy_S = prm % Q_sol
2020-03-16 19:52:44 +05:30
criticalStress_S = activationEnergy_S / activationVolume_S
tauRel_S = min ( 1.0_pReal , tauEff / criticalStress_S ) ! ensure that the activation probability cannot become greater than one
2020-09-22 21:41:09 +05:30
tSolidSolution = 1.0_pReal / prm % nu_a &
2020-03-16 19:52:44 +05:30
* exp ( activationEnergy_S / ( kB * Temperature ) * ( 1.0_pReal - tauRel_S ** prm % p ) ** prm % q )
if ( tauEff < criticalStress_S ) then
dtSolidSolution_dtau = tSolidSolution * prm % p * prm % q * activationVolume_S / ( kB * Temperature ) &
* ( 1.0_pReal - tauRel_S ** prm % p ) ** ( prm % q - 1.0_pReal ) * tauRel_S ** ( prm % p - 1.0_pReal )
else
dtSolidSolution_dtau = 0.0_pReal
endif
!* viscous glide velocity
tauEff = abs ( tau ( s ) ) - tauThreshold ( s )
2020-09-22 21:41:09 +05:30
mobility = prm % b_sl ( s ) / prm % eta
2020-03-16 19:52:44 +05:30
vViscous = mobility * tauEff
!* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of
!* free flight at glide velocity in between.
!* adopt sign from resolved stress
v ( s ) = sign ( 1.0_pReal , tau ( s ) ) &
/ ( tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous )
dv_dtau ( s ) = v ( s ) ** 2.0_pReal * ( dtSolidSolution_dtau / meanfreepath_S + mobility / vViscous ** 2.0_pReal )
dv_dtauNS ( s ) = v ( s ) ** 2.0_pReal * dtPeierls_dtau / meanfreepath_P
endif
enddo
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
!--------------------------------------------------------------------------------------------------
2021-02-14 21:59:23 +05:30
pure function getRho ( ph , me )
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
integer , intent ( in ) :: ph , me
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 10 ) :: getRho
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) )
2020-03-16 12:52:36 +05:30
2021-02-14 23:00:57 +05:30
getRho = reshape ( state ( ph ) % rho ( : , me ) , [ 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)
getRho ( : , mob ) = max ( getRho ( : , mob ) , 0.0_pReal )
getRho ( : , dip ) = max ( getRho ( : , dip ) , 0.0_pReal )
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
where ( abs ( getRho ) < max ( prm % rho_min / geom ( ph ) % V_0 ( me ) ** ( 2.0_pReal / 3.0_pReal ) , prm % rho_significant ) ) &
getRho = 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
!--------------------------------------------------------------------------------------------------
2021-02-14 21:59:23 +05:30
pure function getRho0 ( ph , me )
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
integer , intent ( in ) :: ph , me
2021-02-14 23:00:57 +05:30
real ( pReal ) , dimension ( param ( ph ) % sum_N_sl , 10 ) :: getRho0
2020-03-16 12:52:36 +05:30
2021-02-14 21:59:23 +05:30
2021-02-14 23:00:57 +05:30
associate ( prm = > param ( ph ) )
2020-03-16 12:52:36 +05:30
2021-02-14 23:00:57 +05:30
getRho0 = reshape ( state0 ( ph ) % rho ( : , me ) , [ 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)
getRho0 ( : , mob ) = max ( getRho0 ( : , mob ) , 0.0_pReal )
getRho0 ( : , dip ) = max ( getRho0 ( : , dip ) , 0.0_pReal )
2020-03-16 12:52:36 +05:30
2021-03-05 01:46:36 +05:30
where ( abs ( getRho0 ) < max ( prm % rho_min / geom ( ph ) % V_0 ( me ) ** ( 2.0_pReal / 3.0_pReal ) , prm % rho_significant ) ) &
2021-02-14 21:59:23 +05:30
getRho0 = 0.0_pReal
2020-03-16 12:52:36 +05:30
end associate
end function getRho0
2021-02-14 19:06:56 +05:30
subroutine storeGeometry ( ph )
integer , intent ( in ) :: ph
2021-03-27 22:30:13 +05:30
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-02-14 19:06:56 +05:30
2021-03-01 10:46:16 +05:30
V = reshape ( IPvolume , [ product ( shape ( IPvolume ) ) ] )
do ce = 1 , size ( material_homogenizationMemberAt2 , 1 )
do co = 1 , homogenization_maxNconstituents
2021-03-05 01:46:36 +05:30
if ( material_phaseAt2 ( co , ce ) == ph ) geom ( ph ) % V_0 ( material_phaseMemberAt2 ( co , ce ) ) = V ( ce )
2021-02-14 19:06:56 +05:30
enddo
enddo
end subroutine
2021-01-26 05:41:32 +05:30
end submodule nonlocal
