DAMASK_EICMD/src/phase_mechanical_plastic_no...

1836 lines
103 KiB
Fortran
Raw Normal View History

!--------------------------------------------------------------------------------------------------
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for plasticity including dislocation flux
!--------------------------------------------------------------------------------------------------
submodule(phase:plastic) nonlocal
use geometry_plastic_nonlocal, only: &
2019-06-07 13:50:56 +05:30
nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
2019-06-06 14:38:58 +05:30
IPvolume => geometry_plastic_nonlocal_IPvolume0, &
IParea => geometry_plastic_nonlocal_IParea0, &
2020-12-20 21:41:43 +05:30
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, &
geometry_plastic_nonlocal_disable
type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0
end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom
2019-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-17 21:32:08 +05:30
! storage order of dislocation types
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
sgl = [1,2,3,4,5,6,7,8] !< signed (single)
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
edg = [1,2,5,6,9], & !< edge
scr = [3,4,7,8,10] !< screw
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
mob = [1,2,3,4], & !< mobile
imm = [5,6,7,8] !< immobile (blocked)
2020-03-29 23:37:09 +05:30
integer, dimension(*), parameter :: &
2019-03-17 21:32:08 +05:30
dip = [9,10], & !< dipole
imm_edg = imm(1:2), & !< immobile edge
imm_scr = imm(3:4) !< immobile screw
integer, parameter :: &
mob_edg_pos = 1, & !< mobile edge positive
mob_edg_neg = 2, & !< mobile edge negative
mob_scr_pos = 3, & !< mobile screw positive
mob_scr_neg = 4 !< mobile screw positive
2020-03-16 20:06:34 +05:30
! BEGIN DEPRECATED
2019-12-04 23:30:56 +05:30
integer, dimension(:,:,:), allocatable :: &
2019-03-17 21:32:08 +05:30
iRhoU, & !< state indices for unblocked density
2020-06-26 15:14:17 +05:30
iV, & !< state indices for dislocation velocities
2019-03-17 21:32:08 +05:30
iD !< state indices for stable dipole height
!END DEPRECATED
2019-12-04 23:30:56 +05:30
real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
compatibility !< slip system compatibility between me and my neighbors
2020-03-17 12:05:25 +05:30
type :: tInitialParameters !< container type for internal constitutive parameters
real(pReal) :: &
2020-09-22 21:41:09 +05:30
sigma_rho_u, & !< standard deviation of scatter in initial dislocation density
random_rho_u, &
random_rho_u_binning
2020-03-17 12:05:25 +05:30
real(pReal), dimension(:), allocatable :: &
2020-09-22 21:41:09 +05:30
rho_u_ed_pos_0, & !< initial edge_pos dislocation density
rho_u_ed_neg_0, & !< initial edge_neg dislocation density
rho_u_sc_pos_0, & !< initial screw_pos dislocation density
rho_u_sc_neg_0, & !< initial screw_neg dislocation density
rho_d_ed_0, & !< initial edge dipole dislocation density
rho_d_sc_0 !< initial screw dipole dislocation density
2020-05-31 01:09:35 +05:30
integer, dimension(:), allocatable :: &
2020-03-17 12:05:25 +05:30
N_sl
end type tInitialParameters
2019-12-04 23:30:56 +05:30
type :: tParameters !< container type for internal constitutive parameters
2019-03-17 21:32:08 +05:30
real(pReal) :: &
2020-09-22 21:41:09 +05:30
V_at, & !< atomic volume
D_0, & !< prefactor for self-diffusion coefficient
Q_cl, & !< activation enthalpy for diffusion
2020-03-14 23:41:26 +05:30
atol_rho, & !< absolute tolerance for dislocation density in state integration
2020-09-22 21:41:09 +05:30
rho_significant, & !< density considered significant
2020-09-23 04:25:19 +05:30
rho_min, & !< number of dislocations considered significant
w, & !< width of a doubkle kink in multiples of the Burgers vector length b
2020-09-22 21:41:09 +05:30
Q_sol, & !< activation energy for solid solution in J
2020-09-23 05:36:03 +05:30
f_sol, & !< solid solution obstacle size in multiples of the Burgers vector length
2020-09-22 21:41:09 +05:30
c_sol, & !< concentration of solid solution in atomic parts
2019-03-17 21:32:08 +05:30
p, & !< parameter for kinetic law (Kocks,Argon,Ashby)
q, & !< parameter for kinetic law (Kocks,Argon,Ashby)
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
2019-12-04 23:30:56 +05:30
type :: tNonlocalMicrostructure
2019-03-17 21:32:08 +05:30
real(pReal), allocatable, dimension(:,:) :: &
tau_pass, &
tau_Back
2019-03-17 21:32:08 +05:30
end type tNonlocalMicrostructure
2019-12-04 23:30:56 +05:30
type :: tNonlocalState
2019-03-17 21:32:08 +05:30
real(pReal), pointer, dimension(:,:) :: &
rho, & ! < all dislocations
rhoSgl, &
rhoSglMobile, & ! iRhoU
2019-03-17 21:32:08 +05:30
rho_sgl_mob_edg_pos, &
rho_sgl_mob_edg_neg, &
rho_sgl_mob_scr_pos, &
rho_sgl_mob_scr_neg, &
2020-03-16 20:06:34 +05:30
rhoSglImmobile, &
2019-03-17 21:32:08 +05:30
rho_sgl_imm_edg_pos, &
rho_sgl_imm_edg_neg, &
rho_sgl_imm_scr_pos, &
rho_sgl_imm_scr_neg, &
2020-03-16 20:06:34 +05:30
rhoDip, &
2019-03-17 21:32:08 +05:30
rho_dip_edg, &
rho_dip_scr, &
rho_forest, &
2019-12-01 14:05:44 +05:30
gamma, &
2019-12-01 13:25:24 +05:30
v, &
v_edg_pos, &
v_edg_neg, &
v_scr_pos, &
v_scr_neg
2019-03-17 21:32:08 +05:30
end type tNonlocalState
2020-07-03 20:15:11 +05:30
2020-07-02 02:21:21 +05:30
type(tNonlocalState), allocatable, dimension(:) :: &
2019-03-17 21:32:08 +05:30
deltaState, &
dotState, &
state, &
state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
2019-12-04 23:30:56 +05:30
type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
2019-02-18 14:58:08 +05:30
contains
!--------------------------------------------------------------------------------------------------
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
2020-08-15 19:32:10 +05:30
module function plastic_nonlocal_init() result(myPlasticity)
2019-02-21 04:54:35 +05:30
2020-08-15 19:32:10 +05:30
logical, dimension(:), allocatable :: myPlasticity
2019-03-17 18:05:41 +05:30
integer :: &
Ninstances, &
2021-02-16 20:36:09 +05:30
ph, &
2021-03-05 01:46:36 +05:30
Nmembers, &
2020-03-17 12:05:25 +05:30
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
2019-02-22 13:51:04 +05:30
s1, s2, &
2020-03-17 12:05:25 +05:30
s, t, l
real(pReal), dimension(:), allocatable :: &
a
2019-12-21 16:58:24 +05:30
character(len=pStringLen) :: &
2020-03-16 20:06:34 +05:30
extmsg = ''
2020-03-17 12:05:25 +05:30
type(tInitialParameters) :: &
ini
2020-08-15 19:32:10 +05:30
class(tNode), pointer :: &
phases, &
phase, &
2020-11-03 03:16:46 +05:30
mech, &
2020-08-15 19:32:10 +05:30
pl
2019-03-09 15:32:12 +05:30
2020-08-15 19:32:10 +05:30
myPlasticity = plastic_active('nonlocal')
Ninstances = count(myPlasticity)
if(Ninstances == 0) then
2020-08-15 19:32:10 +05:30
call geometry_plastic_nonlocal_disable
return
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:333348, 2014'
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
phases => config_material%get('phase')
allocate(geom(phases%length))
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-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-08-15 19:32:10 +05:30
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0e4_pReal)
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)
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))
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.
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
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
2021-03-11 22:30:07 +05:30
ini%rho_u_ed_pos_0 = pl%get_as1dFloat('rho_u_ed_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_ed_neg_0 = pl%get_as1dFloat('rho_u_ed_neg_0', requiredSize=size(ini%N_sl))
ini%rho_u_sc_pos_0 = pl%get_as1dFloat('rho_u_sc_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_sc_neg_0 = pl%get_as1dFloat('rho_u_sc_neg_0', requiredSize=size(ini%N_sl))
ini%rho_d_ed_0 = pl%get_as1dFloat('rho_d_ed_0', requiredSize=size(ini%N_sl))
ini%rho_d_sc_0 = pl%get_as1dFloat('rho_d_sc_0', requiredSize=size(ini%N_sl))
2021-03-11 22:30:07 +05:30
prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(ini%N_sl))
prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(ini%N_sl))
2020-09-22 21:41:09 +05:30
prm%i_sl = math_expand(prm%i_sl,ini%N_sl)
prm%b_sl = math_expand(prm%b_sl,ini%N_sl)
2021-03-11 22:30:07 +05:30
prm%d_ed = pl%get_as1dFloat('d_ed', requiredSize=size(ini%N_sl))
prm%d_sc = pl%get_as1dFloat('d_sc', requiredSize=size(ini%N_sl))
2020-09-22 21:41:09 +05:30
prm%d_ed = math_expand(prm%d_ed,ini%N_sl)
prm%d_sc = math_expand(prm%d_sc,ini%N_sl)
2020-03-16 23:36:10 +05:30
allocate(prm%minDipoleHeight(prm%sum_N_sl,2))
2020-09-22 21:41:09 +05:30
prm%minDipoleHeight(:,1) = prm%d_ed
prm%minDipoleHeight(:,2) = prm%d_sc
2021-03-11 22:30:07 +05:30
prm%tau_Peierls_ed = pl%get_as1dFloat('tau_Peierls_ed', requiredSize=size(ini%N_sl))
prm%tau_Peierls_sc = pl%get_as1dFloat('tau_Peierls_sc', requiredSize=size(ini%N_sl))
2020-09-23 04:25:19 +05:30
prm%tau_Peierls_ed = math_expand(prm%tau_Peierls_ed,ini%N_sl)
prm%tau_Peierls_sc = math_expand(prm%tau_Peierls_sc,ini%N_sl)
2020-03-16 23:36:10 +05:30
allocate(prm%peierlsstress(prm%sum_N_sl,2))
2020-09-23 04:25:19 +05:30
prm%peierlsstress(:,1) = prm%tau_Peierls_ed
prm%peierlsstress(:,2) = prm%tau_Peierls_sc
2020-09-22 21:41:09 +05:30
prm%rho_significant = pl%get_asFloat('rho_significant')
2020-09-23 04:25:19 +05:30
prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal)
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-09-22 21:41:09 +05:30
prm%chi_surface = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal)
prm%chi_GB = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal)
prm%f_ed_mult = pl%get_asFloat('f_ed_mult')
2020-08-15 19:32:10 +05:30
prm%shortRangeStressCorrection = pl%get_asBool('short_range_stress_correction', defaultVal = .false.)
2020-04-01 14:39:30 +05:30
2019-02-21 04:54:35 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
2020-09-22 22:56:39 +05:30
if (any(prm%b_sl < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(ini%rho_u_ed_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0'
if (any(ini%rho_u_ed_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0'
if (any(ini%rho_u_sc_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0'
if (any(ini%rho_u_sc_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0'
if (any(ini%rho_d_ed_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0'
if (any(ini%rho_d_sc_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0'
if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' d_ed or d_sc'
if (prm%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'
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-09-22 21:41:09 +05:30
if (prm%Q_sol <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol'
if (prm%f_sol <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol'
if (prm%c_sol <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol'
2019-02-21 04:54:35 +05:30
2020-09-22 21:41:09 +05:30
if (prm%chi_GB > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB'
if (prm%chi_surface < 0.0_pReal .or. prm%chi_surface > 1.0_pReal) &
extmsg = trim(extmsg)//' chi_surface'
2020-09-22 21:41:09 +05:30
if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed_mult'
2019-02-21 04:54:35 +05:30
endif slipActive
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 &
+ size([ 'velocityEdgePos ','velocityEdgeNeg ', &
2019-03-17 18:05:41 +05:30
'velocityScrewPos ','velocityScrewNeg ', &
2020-03-16 23:36:10 +05:30
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
2019-02-21 04:54:35 +05:30
sizeDeltaState = sizeDotState
2021-03-05 01:46:36 +05:30
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
2019-12-21 16:58:24 +05:30
2021-03-05 01:46:36 +05:30
allocate(geom(ph)%V_0(Nmembers))
2021-02-16 20:36:09 +05:30
call storeGeometry(ph)
2021-02-16 20:36:09 +05:30
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
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)
2021-03-05 01:46:36 +05:30
stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
2014-07-08 20:28:23 +05:30
2021-03-05 01:46:36 +05:30
allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
2019-03-17 21:32:08 +05:30
end associate
2021-03-05 01:46:36 +05:30
if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
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-03-16 23:36:10 +05:30
allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,&
discretization_nIPs,discretization_Nelems), source=0.0_pReal)
2019-02-22 14:32:43 +05:30
! BEGIN DEPRECATED----------------------------------------------------------------------------------
allocate(iRhoU(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0)
2019-02-22 14:32:43 +05:30
2021-02-16 20:36:09 +05:30
do ph = 1, phases%length
2020-08-15 19:32:10 +05:30
2021-02-16 20:36:09 +05:30
if(.not. myPlasticity(ph)) cycle
2021-02-16 20:36:09 +05:30
phase => phases%get(ph)
2021-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)'
enddo
2020-08-15 19:32:10 +05:30
end function plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
module subroutine nonlocal_dependentState(ph, me, ip, el)
2020-03-16 21:59:47 +05:30
integer, intent(in) :: &
ph, &
2021-01-26 13:09:17 +05:30
me, &
2020-03-16 21:59:47 +05:30
ip, &
el
2019-03-17 21:32:08 +05:30
integer :: &
no, & !< neighbor offset
neighbor_el, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point
c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index
dir, &
2020-07-02 00:52:05 +05:30
n
2019-03-17 21:32:08 +05:30
real(pReal) :: &
FVsize, &
nRealNeighbors ! number of really existing neighbors
integer, dimension(2) :: &
neighbors
real(pReal), dimension(2) :: &
rhoExcessGradient, &
rhoExcessGradient_over_rho, &
rhoTotal
real(pReal), dimension(3) :: &
rhoExcessDifferences, &
normal_latticeConf
real(pReal), dimension(3,3) :: &
invFe, & !< inverse of elastic deformation gradient
invFp, & !< inverse of plastic deformation gradient
connections, &
invConnections
2019-06-07 13:50:56 +05:30
real(pReal), dimension(3,nIPneighbors) :: &
2019-03-17 21:32:08 +05:30
connection_latticeConf
real(pReal), dimension(2,param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
rhoExcess
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
rho_edg_delta, &
rho_scr_delta
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho, &
rho0, &
rho_neighbor0
real(pReal), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(param(ph)%sum_N_sl,nIPneighbors) :: &
2019-03-17 21:32:08 +05:30
rho_edg_delta_neighbor, &
rho_scr_delta_neighbor
2020-03-16 23:36:10 +05:30
real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: &
2019-03-17 21:32:08 +05:30
neighbor_rhoExcess, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,param(ph)%sum_N_sl,2) :: &
2019-03-17 21:32:08 +05:30
m ! direction of dislocation motion
associate(prm => param(ph),dst => microstructure(ph), stt => state(ph))
2021-02-14 21:59:23 +05:30
rho = getRho(ph,me)
2021-01-26 13:09:17 +05:30
stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+ matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
! coefficients are corrected for the line tension effect
2019-03-17 21:32:08 +05:30
! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
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 &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities
!*** zero for material points of local plasticity
2019-02-20 19:24:26 +05:30
!#################################################################################################
! ToDo: MD: this is most likely only correct for F_i = I
!#################################################################################################
2021-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))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
2020-03-16 13:06:19 +05:30
rhoExcess(1,:) = rho_edg_delta
rhoExcess(2,:) = rho_scr_delta
2021-02-14 21:59:23 +05:30
FVsize = geom(ph)%V_0(me) ** (1.0_pReal/3.0_pReal)
2019-03-17 21:32:08 +05:30
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
2019-03-17 21:32:08 +05:30
nRealNeighbors = 0.0_pReal
neighbor_rhoTotal = 0.0_pReal
2019-06-07 13:50:56 +05:30
do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
2019-06-14 12:47:05 +05:30
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
2019-03-17 21:32:08 +05:30
if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseAt(1,neighbor_el) == ph) then
2019-03-17 21:32:08 +05:30
nRealNeighbors = nRealNeighbors + 1.0_pReal
2021-02-14 21:59:23 +05:30
rho_neighbor0 = getRho0(ph,no)
rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
2019-09-28 02:37:03 +05:30
connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
- discretization_IPcoords(1:3,el+neighbor_ip-1))
2019-06-07 14:03:49 +05:30
normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el))
if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
2019-03-17 21:32:08 +05:30
else
! local neighbor or different lattice structure or different constitution instance -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal
rho_edg_delta_neighbor(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta
endif
else
2019-03-17 21:32:08 +05:30
! free surface -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal
rho_edg_delta_neighbor(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta
endif
2019-03-17 21:32:08 +05:30
enddo
2019-03-17 21:32:08 +05:30
neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor
neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor
2019-03-17 21:32:08 +05:30
!* loop through the slip systems and calculate the dislocation gradient by
!* 1. interpolation of the excess density in the neighorhood
!* 2. interpolation of the dead dislocation density in the central volume
2020-03-16 13:06:19 +05:30
m(1:3,:,1) = prm%slip_direction
m(1:3,:,2) = -prm%slip_transverse
2020-03-16 23:36:10 +05:30
do s = 1,prm%sum_N_sl
2019-02-21 23:48:06 +05:30
! gradient from interpolation of neighboring excess density ...
2019-03-17 21:32:08 +05:30
do c = 1,2
do dir = 1,3
neighbors(1) = 2 * dir - 1
neighbors(2) = 2 * dir
connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) &
- connection_latticeConf(1:3,neighbors(2))
rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) &
- neighbor_rhoExcess(c,s,neighbors(2))
enddo
invConnections = math_inv33(connections)
if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
enddo
2019-03-17 21:32:08 +05:30
! ... plus gradient from deads ...
rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize
rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize
2019-03-17 21:32:08 +05:30
! ... normalized with the total density ...
rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pReal + nRealNeighbors)
rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pReal + nRealNeighbors)
2019-03-17 21:32:08 +05:30
rhoExcessGradient_over_rho = 0.0_pReal
2020-03-16 13:06:19 +05:30
where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
2019-03-17 21:32:08 +05:30
! ... gives the local stress correction when multiplied with a factor
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
#ifdef DEBUG
if (debugConstitutive%extensive &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then
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
endif
#endif
end associate
2019-02-21 23:48:06 +05:30
2021-01-26 16:47:00 +05:30
end subroutine nonlocal_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
2021-01-26 12:03:04 +05:30
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
2021-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) :: &
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
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
2019-03-17 21:32:08 +05:30
rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
2019-03-17 21:32:08 +05:30
v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate
associate(prm => param(ph),dst=>microstructure(ph),&
stt=>state(ph))
2020-03-16 23:36:10 +05:30
ns = prm%sum_N_sl
!*** 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)
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)
! edges
2020-03-16 19:52:44 +05:30
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
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)
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), &
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, ph)
2019-03-17 21:32:08 +05:30
enddo
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)
endif
2019-12-01 13:25:24 +05:30
2021-01-26 13:09:17 +05:30
stt%v(:,me) = pack(v,.true.)
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-09-22 21:41:09 +05:30
gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
2019-03-17 21:32:08 +05:30
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
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) &
* (+ 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)
enddo
2019-03-17 21:32:08 +05:30
end associate
2021-01-26 12:03:04 +05:30
end subroutine nonlocal_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
2021-02-14 21:59:23 +05:30
module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
2020-03-16 21:59:47 +05:30
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
2019-03-17 21:32:08 +05:30
integer, intent(in) :: &
ph, &
2021-02-14 21:59:23 +05:30
me
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
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho ! current dislocation densities
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
2019-03-17 21:32:08 +05:30
v ! dislocation glide velocity
real(pReal), dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
tau ! current resolved shear stress
real(pReal), dimension(param(ph)%sum_N_sl,2) :: &
2019-03-17 21:32:08 +05:30
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
dUpper, & ! current maximum stable dipole distance for edges and screws
dUpperOld, & ! old maximum stable dipole distance for edges and screws
deltaDUpper ! change in maximum stable dipole distance for edges and screws
associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph))
ns = prm%sum_N_sl
!*** shortcut to state variables
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)
2021-02-14 21:59:23 +05:30
rho = getRho(ph,me)
2019-03-17 21:32:08 +05:30
rhoDip = rho(:,dip)
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** dislocation remobilization (bauschinger effect)
where(rho(:,imm) * v < 0.0_pReal)
deltaRhoRemobilization(:,mob) = abs(rho(:,imm))
deltaRhoRemobilization(:,imm) = - rho(:,imm)
rho(:,mob) = rho(:,mob) + abs(rho(:,imm))
rho(:,imm) = 0.0_pReal
elsewhere
deltaRhoRemobilization(:,mob) = 0.0_pReal
deltaRhoRemobilization(:,imm) = 0.0_pReal
endwhere
2019-03-17 21:32:08 +05:30
deltaRhoRemobilization(:,dip) = 0.0_pReal
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** calculate dipole formation and dissociation by stress change
2019-03-17 21:32:08 +05:30
!*** calculate limits for stable dipole height
2020-03-16 23:36:10 +05:30
do s = 1,prm%sum_N_sl
2021-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-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))
2019-03-17 21:32:08 +05:30
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
2020-03-16 13:06:19 +05:30
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
2019-03-17 21:32:08 +05:30
where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) &
2020-03-16 13:06:19 +05:30
dUpper(:,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(:,2))
2019-03-17 21:32:08 +05:30
dUpper = max(dUpper,prm%minDipoleHeight)
deltaDUpper = dUpper - dUpperOld
2019-03-17 21:32:08 +05:30
!*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal
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-03-16 13:06:19 +05:30
forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),me) = dUpper(s,c)
2021-01-26 13:09:17 +05:30
plasticState(ph)%deltaState(:,me) = 0.0_pReal
del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
2019-03-17 21:32:08 +05:30
end associate
end subroutine plastic_nonlocal_deltaState
2019-02-20 05:11:44 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
2021-01-26 16:47:00 +05:30
module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
ph,me,ip,el)
2020-03-16 14:19:59 +05:30
real(pReal), dimension(3,3), intent(in) :: &
2019-03-17 21:32:08 +05:30
Mp !< MandelStress
2020-03-16 14:19:59 +05:30
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
2021-01-26 13:09:17 +05:30
me, &
2020-03-16 14:19:59 +05:30
ip, & !< current integration point
el !< current element number
2019-03-17 21:32:08 +05:30
integer :: &
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
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
2019-03-17 21:32:08 +05:30
rho, &
2020-02-11 10:11:10 +05:30
rho0, & !< dislocation density at beginning of time step
2019-03-17 21:32:08 +05:30
rhoDot, & !< density evolution
rhoDotMultiplication, & !< density evolution by multiplication
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
2019-03-17 21:32:08 +05:30
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
2020-02-11 10:11:10 +05:30
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
2019-03-17 21:32:08 +05:30
v, & !< current dislocation glide velocity
v0, &
2019-03-17 21:32:08 +05:30
gdot !< shear rates
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
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
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
endif
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
2019-03-17 21:32:08 +05:30
tau = 0.0_pReal
gdot = 0.0_pReal
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)
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)
2019-03-17 21:32:08 +05:30
#ifdef DEBUG
if (debugConstitutive%basic &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) &
.or. .not. debugConstitutive%selective)) then
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
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
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))
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)
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
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
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me)
2019-03-17 21:32:08 +05:30
!****************************************************************************
!*** calculate dipole formation and annihilation
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
2019-03-17 21:32:08 +05:30
!*** athermal annihilation
rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) &
2020-09-22 21:41:09 +05:30
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
* ( 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-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
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
rhoDot = rhoDotFlux(timestep, ph,me,ip,el) &
2019-03-17 21:32:08 +05:30
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
+ rhoDotThermalAnnihilation
2020-03-16 13:06:19 +05:30
if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) &
.or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
print'(a)', '<< CONST >> enforcing cutback !!!'
2019-03-17 21:32:08 +05:30
endif
#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
2019-03-17 21:32:08 +05:30
end associate
2021-01-26 16:47:00 +05:30
end subroutine nonlocal_dotState
2019-03-17 21:32:08 +05:30
2020-03-16 19:52:44 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
function rhoDotFlux(timestep,ph,me,ip,el)
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
2021-01-26 13:09:17 +05:30
me, &
ip, & !< current integration point
el !< current element number
2020-07-03 20:15:11 +05:30
integer :: &
neighbor_ph, & !< phase of my neighbor's plasticity
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
n, & !< index of my current neighbor
neighbor_el, & !< element number of my neighbor
neighbor_ip, & !< integration point of my neighbor
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
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
rho, &
rho0, & !< dislocation density at beginning of time step
2020-05-03 14:24:57 +05:30
rhoDotFlux !< density evolution by flux
real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity
v0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates
real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: &
m !< direction of dislocation motion
real(pReal), dimension(3,3) :: &
my_F, & !< my total deformation gradient
neighbor_F, & !< total deformation gradient of my neighbor
my_Fe, & !< my elastic deformation gradient
neighbor_Fe, & !< elastic deformation gradient of my neighbor
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
associate(prm => param(ph), &
dst => microstructure(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
gdot = 0.0_pReal
2021-02-14 21:59:23 +05:30
rho = getRho(ph,me)
rhoSgl = rho(:,sgl)
2021-02-14 21:59:23 +05:30
rho0 = getRho0(ph,me)
my_rhoSgl0 = rho0(:,sgl)
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)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me)
!****************************************************************************
!*** 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 &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal &
2020-09-22 21:41:09 +05:30
.and. prm%f_c * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!'
endif
#endif
2020-05-03 15:34:57 +05:30
rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback
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)))
neighbors: do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
neighbor_n = IPneighborhood(3,n,ip,el)
np = material_phaseAt(1,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_ph = material_phaseAt(1,neighbor_el)
2021-02-09 03:51:53 +05:30
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average
Favg = my_F
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)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal)
endforall
2020-09-23 04:25:19 +05:30
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
2020-09-22 21:41:09 +05:30
.or. neighbor_rhoSgl0 < prm%rho_significant) &
neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
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
2019-03-17 21:32:08 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-17 22:29:01 +05:30
!> @brief Compatibility update
!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip
2019-03-17 21:32:08 +05:30
! plane normals and signed cosine of the angle between the slip directions. Only the largest values
! that sum up to a total of 1 are considered, all others are set to zero.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
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) :: &
ph, &
2019-03-17 21:32:08 +05:30
i, &
e
2019-03-17 21:32:08 +05:30
integer :: &
n, & ! neighbor index
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)
real(pReal), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: &
my_compatibility ! my_compatibility for current element and ip
2019-03-17 21:32:08 +05:30
real(pReal) :: &
my_compatibilitySum, &
thresholdValue, &
nThresholdValues
logical, dimension(param(ph)%sum_N_sl) :: &
2019-03-17 21:32:08 +05:30
belowThreshold
2019-12-02 17:28:23 +05:30
type(rotation) :: mis
2019-03-17 21:32:08 +05:30
associate(prm => param(ph))
2020-03-16 23:36:10 +05:30
ns = prm%sum_N_sl
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-03-16 13:06:19 +05:30
neighbors: do n = 1,nIPneighbors
neighbor_e = IPneighborhood(1,n,i,e)
neighbor_i = IPneighborhood(2,n,i,e)
neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
2020-03-16 14:07:25 +05:30
neighbor_phase = material_phaseAt(1,neighbor_e)
2019-03-17 21:32:08 +05:30
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
2020-03-16 14:07:25 +05:30
!* FREE SURFACE
!* 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)
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))))
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
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
2019-03-17 21:32:08 +05:30
enddo neighbors
2020-03-16 13:06:19 +05:30
compatibility(:,:,:,:,i,e) = my_compatibility
2019-03-17 21:32:08 +05:30
end associate
end subroutine plastic_nonlocal_updateCompatibility
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_results(ph,group)
integer, intent(in) :: ph
character(len=*),intent(in) :: group
2020-02-14 13:56:26 +05:30
integer :: o
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')
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, &
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
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
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)]
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
!--------------------------------------------------------------------------------------------------
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)
ph
2020-03-16 19:52:44 +05:30
real(pReal), intent(in) :: &
Temperature !< temperature
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
2020-03-16 19:52:44 +05:30
tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including non Schmid effects)
tauThreshold !< threshold shear stress
real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
2020-03-16 19:52:44 +05:30
v, & !< velocity
dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
integer :: &
s !< index of my current slip system
real(pReal) :: &
tauRel_P, &
tauRel_S, &
tauEff, & !< effective shear stress
tPeierls, & !< waiting time in front of a peierls barriers
tSolidSolution, & !< waiting time in front of a solid solution obstacle
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
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
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
associate(prm => param(ph))
2020-03-16 12:52:36 +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
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
associate(prm => param(ph))
2020-03-16 12:52:36 +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
subroutine storeGeometry(ph)
integer, intent(in) :: ph
integer :: ce, co
2021-03-01 10:46:16 +05:30
real(pReal), dimension(:), allocatable :: V
2021-03-05 01:46:36 +05:30
2021-03-01 10:46:16 +05:30
V = reshape(IPvolume,[product(shape(IPvolume))])
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)
enddo
enddo
end subroutine
2021-01-26 05:41:32 +05:30
end submodule nonlocal