1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
19 real(pReal),
parameter :: &
23 integer,
dimension(8),
parameter :: &
24 sgl = [1,2,3,4,5,6,7,8]
25 integer,
dimension(5),
parameter :: &
26 edg = [1,2,5,6,9], & !< edge
28 integer,
dimension(4),
parameter :: &
29 mob = [1,2,3,4], & !< mobile
31 integer,
dimension(2),
parameter :: &
32 dip = [9,10], & !< dipole
35 integer,
parameter :: &
36 mob_edg_pos = 1, & !< mobile edge positive
42 integer,
dimension(:,:,:),
allocatable :: &
43 iRhoU, & !< state indices for unblocked density
44 iV, & !< state indices for dislcation velocities
48 real(pReal),
dimension(:,:,:,:,:,:),
allocatable :: &
51 type :: tinitialparameters
53 rhoSglScatter, & !< standard deviation of scatter in initial dislocation density
56 real(pReal),
dimension(:),
allocatable :: &
57 rhoSglEdgePos0, & !< initial edge_pos dislocation density
58 rhoSglEdgeNeg0, & !< initial edge_neg dislocation density
59 rhoSglScrewPos0, & !< initial screw_pos dislocation density
60 rhoSglScrewNeg0, & !< initial screw_neg dislocation density
61 rhoDipEdge0, & !< initial edge dipole dislocation density
63 integer,
dimension(:) ,
allocatable :: &
65 end type tinitialparameters
69 atomicVolume, & !< atomic volume
70 Dsd0, & !< prefactor for self-diffusion coefficient
71 selfDiffusionEnergy, & !< activation enthalpy for diffusion
72 atol_rho, & !< absolute tolerance for dislocation density in state integration
73 significantRho, & !< density considered significant
74 significantN, & !< number of dislocations considered significant
75 doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b
76 solidSolutionEnergy, & !< activation energy for solid solution in J
77 solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length
78 solidSolutionConcentration, & !< concentration of solid solution in atomic parts
79 p, & !< parameter for kinetic law (Kocks,Argon,Ashby)
80 q, & !< parameter for kinetic law (Kocks,Argon,Ashby)
81 viscosity, & !< viscosity for dislocation glide in Pa s
82 fattack, & !< attack frequency in Hz
83 surfaceTransmissivity, & !< transmissivity at free surface
84 grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture)
85 CFLfactor, & !< safety factor for CFL flux condition
86 fEdgeMultiplication, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
91 real(pReal),
dimension(:),
allocatable :: &
92 minDipoleHeight_edge, & !< minimum stable edge dipole height
93 minDipoleHeight_screw, & !< minimum stable screw dipole height
95 peierlsstress_screw, &
96 lambda0, & !< mean free path prefactor for each
98 real(pReal),
dimension(:,:),
allocatable :: &
104 interactionSlipSlip ,& !< coefficients for slip-slip interaction
105 forestProjection_Edge, & !< matrix of forest projections of edge dislocations
106 forestProjection_Screw
107 real(pReal),
dimension(:,:,:),
allocatable :: &
108 Schmid, & !< Schmid contribution
113 integer,
dimension(:),
allocatable :: &
115 character(len=pStringLen),
dimension(:),
allocatable :: &
118 shortRangeStressCorrection, & !< use of short range stress correction by excess density gradient term
119 nonSchmidActive = .false.
122 type :: tnonlocalmicrostructure
123 real(pReal),
allocatable,
dimension(:,:) :: &
126 end type tnonlocalmicrostructure
128 type :: tnonlocalstate
129 real(pReal),
pointer,
dimension(:,:) :: &
133 rho_sgl_mob_edg_pos, &
134 rho_sgl_mob_edg_neg, &
135 rho_sgl_mob_scr_pos, &
136 rho_sgl_mob_scr_neg, &
138 rho_sgl_imm_edg_pos, &
139 rho_sgl_imm_edg_neg, &
140 rho_sgl_imm_scr_pos, &
141 rho_sgl_imm_scr_neg, &
152 end type tnonlocalstate
154 type(tNonlocalState),
allocatable,
dimension(:) :: &
160 type(tParameters),
dimension(:),
allocatable :: param
162 type(tNonlocalMicrostructure),
dimension(:),
allocatable :: microstructure
170 module subroutine plastic_nonlocal_init
176 sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
179 real(pReal),
dimension(:),
allocatable :: &
181 character(len=pStringLen) :: &
183 type(tInitialParameters) :: &
186 write(6,
'(/,a)')
' <<<+- constitutive_'//plasticity_nonlocal_label//
' init -+>>>';
flush(6)
188 write(6,
'(/,a)')
' Reuber et al., Acta Materialia 71:333–348, 2014'
189 write(6,
'(a)')
' https://doi.org/10.1016/j.actamat.2014.03.012'
191 write(6,
'(/,a)')
' Kords, Dissertation RWTH Aachen, 2014'
192 write(6,
'(a)')
' http://publications.rwth-aachen.de/record/229993'
194 ninstance = count(phase_plasticity == plasticity_nonlocal_id)
195 if (iand(debug_level(debug_constitutive),debug_levelbasic) /= 0) &
196 write(6,
'(a16,1x,i5,/)')
'# instances:',ninstance
198 allocate(param(ninstance))
199 allocate(state(ninstance))
200 allocate(state0(ninstance))
201 allocate(dotstate(ninstance))
202 allocate(deltastate(ninstance))
203 allocate(microstructure(ninstance))
205 do p=1,
size(config_phase)
206 if (phase_plasticity(p) /= plasticity_nonlocal_id) cycle
208 associate(prm => param(phase_plasticityinstance(p)), &
209 dot => dotstate(phase_plasticityinstance(p)), &
210 stt => state(phase_plasticityinstance(p)), &
211 st0 => state0(phase_plasticityinstance(p)), &
212 del => deltastate(phase_plasticityinstance(p)), &
213 dst => microstructure(phase_plasticityinstance(p)), &
214 config => config_phase(p))
216 prm%output =
config%getStrings(
'(output)',defaultval=emptystringarray)
218 prm%atol_rho =
config%getFloat(
'atol_rho',defaultval=1.0e4_preal)
221 prm%mu = lattice_mu(p)
222 prm%nu = lattice_nu(p)
224 ini%N_sl =
config%getInts(
'nslip',defaultval=emptyintarray)
225 prm%sum_N_sl = sum(abs(ini%N_sl))
226 slipactive:
if (prm%sum_N_sl > 0)
then
227 prm%Schmid = lattice_schmidmatrix_slip(ini%N_sl,
config%getString(
'lattice_structure'),&
228 config%getFloat(
'c/a',defaultval=0.0_preal))
230 if(trim(
config%getString(
'lattice_structure')) ==
'bcc')
then
231 a =
config%getFloats(
'nonschmid_coefficients',defaultval = emptyrealarray)
232 if(
size(a) > 0) prm%nonSchmidActive = .true.
233 prm%nonSchmid_pos = lattice_nonschmidmatrix(ini%N_sl,a,+1)
234 prm%nonSchmid_neg = lattice_nonschmidmatrix(ini%N_sl,a,-1)
236 prm%nonSchmid_pos = prm%Schmid
237 prm%nonSchmid_neg = prm%Schmid
240 prm%interactionSlipSlip = lattice_interaction_slipbyslip(ini%N_sl, &
241 config%getFloats(
'interaction_slipslip'), &
242 config%getString(
'lattice_structure'))
244 prm%forestProjection_edge = lattice_forestprojection_edge(ini%N_sl,
config%getString(
'lattice_structure'),&
245 config%getFloat(
'c/a',defaultval=0.0_preal))
246 prm%forestProjection_screw = lattice_forestprojection_screw(ini%N_sl,
config%getString(
'lattice_structure'),&
247 config%getFloat(
'c/a',defaultval=0.0_preal))
249 prm%slip_direction = lattice_slip_direction(ini%N_sl,
config%getString(
'lattice_structure'),&
250 config%getFloat(
'c/a',defaultval=0.0_preal))
251 prm%slip_transverse = lattice_slip_transverse(ini%N_sl,
config%getString(
'lattice_structure'),&
252 config%getFloat(
'c/a',defaultval=0.0_preal))
253 prm%slip_normal = lattice_slip_normal(ini%N_sl,
config%getString(
'lattice_structure'),&
254 config%getFloat(
'c/a',defaultval=0.0_preal))
257 allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)
258 do s1 = 1, prm%sum_N_sl
259 do s2 = 1, prm%sum_N_sl
260 if (all(deq0(math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. &
261 any(dneq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) &
262 prm%colinearSystem(s1) = s2
266 ini%rhoSglEdgePos0 =
config%getFloats(
'rhosgledgepos0', requiredsize=
size(ini%N_sl))
267 ini%rhoSglEdgeNeg0 =
config%getFloats(
'rhosgledgeneg0', requiredsize=
size(ini%N_sl))
268 ini%rhoSglScrewPos0 =
config%getFloats(
'rhosglscrewpos0', requiredsize=
size(ini%N_sl))
269 ini%rhoSglScrewNeg0 =
config%getFloats(
'rhosglscrewneg0', requiredsize=
size(ini%N_sl))
270 ini%rhoDipEdge0 =
config%getFloats(
'rhodipedge0', requiredsize=
size(ini%N_sl))
271 ini%rhoDipScrew0 =
config%getFloats(
'rhodipscrew0', requiredsize=
size(ini%N_sl))
273 prm%lambda0 =
config%getFloats(
'lambda0', requiredsize=
size(ini%N_sl))
274 prm%burgers =
config%getFloats(
'burgers', requiredsize=
size(ini%N_sl))
276 prm%lambda0 = math_expand(prm%lambda0,ini%N_sl)
277 prm%burgers = math_expand(prm%burgers,ini%N_sl)
279 prm%minDipoleHeight_edge =
config%getFloats(
'minimumdipoleheightedge', requiredsize=
size(ini%N_sl))
280 prm%minDipoleHeight_screw =
config%getFloats(
'minimumdipoleheightscrew', requiredsize=
size(ini%N_sl))
281 prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge, ini%N_sl)
282 prm%minDipoleHeight_screw = math_expand(prm%minDipoleHeight_screw,ini%N_sl)
283 allocate(prm%minDipoleHeight(prm%sum_N_sl,2))
284 prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge
285 prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw
287 prm%peierlsstress_edge =
config%getFloats(
'peierlsstressedge', requiredsize=
size(ini%N_sl))
288 prm%peierlsstress_screw =
config%getFloats(
'peierlsstressscrew', requiredsize=
size(ini%N_sl))
289 prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge, ini%N_sl)
290 prm%peierlsstress_screw = math_expand(prm%peierlsstress_screw,ini%N_sl)
291 allocate(prm%peierlsstress(prm%sum_N_sl,2))
292 prm%peierlsstress(:,1) = prm%peierlsstress_edge
293 prm%peierlsstress(:,2) = prm%peierlsstress_screw
295 prm%significantRho =
config%getFloat(
'significantrho')
296 prm%significantN =
config%getFloat(
'significantn', 0.0_preal)
297 prm%CFLfactor =
config%getFloat(
'cflfactor',defaultval=2.0_preal)
299 prm%atomicVolume =
config%getFloat(
'atomicvolume')
300 prm%Dsd0 =
config%getFloat(
'selfdiffusionprefactor')
301 prm%selfDiffusionEnergy =
config%getFloat(
'selfdiffusionenergy')
302 prm%linetensionEffect =
config%getFloat(
'linetension')
303 prm%edgeJogFactor =
config%getFloat(
'edgejog')
304 prm%doublekinkwidth =
config%getFloat(
'doublekinkwidth')
305 prm%solidSolutionEnergy =
config%getFloat(
'solidsolutionenergy')
306 prm%solidSolutionSize =
config%getFloat(
'solidsolutionsize')
307 prm%solidSolutionConcentration =
config%getFloat(
'solidsolutionconcentration')
309 prm%p =
config%getFloat(
'p')
310 prm%q =
config%getFloat(
'q')
311 prm%viscosity =
config%getFloat(
'viscosity')
312 prm%fattack =
config%getFloat(
'attackfrequency')
315 ini%rhoSglScatter =
config%getFloat(
'rhosglscatter')
316 ini%rhoSglRandom =
config%getFloat(
'rhosglrandom',0.0_preal)
317 if (
config%keyExists(
'/rhosglrandom/')) &
318 ini%rhoSglRandomBinning =
config%getFloat(
'rhosglrandombinning',0.0_preal)
322 prm%surfaceTransmissivity =
config%getFloat(
'surfacetransmissivity',defaultval=1.0_preal)
323 prm%grainboundaryTransmissivity =
config%getFloat(
'grainboundarytransmissivity',defaultval=-1.0_preal)
324 prm%fEdgeMultiplication =
config%getFloat(
'edgemultiplication')
325 prm%shortRangeStressCorrection =
config%keyExists(
'/shortrangestresscorrection/')
329 if (any(prm%burgers < 0.0_preal)) extmsg = trim(extmsg)//
' burgers'
330 if (any(prm%lambda0 <= 0.0_preal)) extmsg = trim(extmsg)//
' lambda0'
332 if (any(ini%rhoSglEdgePos0 < 0.0_preal)) extmsg = trim(extmsg)//
' rhoSglEdgePos0'
333 if (any(ini%rhoSglEdgeNeg0 < 0.0_preal)) extmsg = trim(extmsg)//
' rhoSglEdgeNeg0'
334 if (any(ini%rhoSglScrewPos0 < 0.0_preal)) extmsg = trim(extmsg)//
' rhoSglScrewPos0'
335 if (any(ini%rhoSglScrewNeg0 < 0.0_preal)) extmsg = trim(extmsg)//
' rhoSglScrewNeg0'
336 if (any(ini%rhoDipEdge0 < 0.0_preal)) extmsg = trim(extmsg)//
' rhoDipEdge0'
337 if (any(ini%rhoDipScrew0 < 0.0_preal)) extmsg = trim(extmsg)//
' rhoDipScrew0'
339 if (any(prm%peierlsstress < 0.0_preal)) extmsg = trim(extmsg)//
' peierlsstress'
340 if (any(prm%minDipoleHeight < 0.0_preal)) extmsg = trim(extmsg)//
' minDipoleHeight'
342 if (prm%viscosity <= 0.0_preal) extmsg = trim(extmsg)//
' viscosity'
343 if (prm%selfDiffusionEnergy <= 0.0_preal) extmsg = trim(extmsg)//
' selfDiffusionEnergy'
344 if (prm%fattack <= 0.0_preal) extmsg = trim(extmsg)//
' fattack'
345 if (prm%doublekinkwidth <= 0.0_preal) extmsg = trim(extmsg)//
' doublekinkwidth'
346 if (prm%Dsd0 < 0.0_preal) extmsg = trim(extmsg)//
' Dsd0'
347 if (prm%atomicVolume <= 0.0_preal) extmsg = trim(extmsg)//
' atomicVolume'
349 if (prm%significantN < 0.0_preal) extmsg = trim(extmsg)//
' significantN'
350 if (prm%significantrho < 0.0_preal) extmsg = trim(extmsg)//
' significantrho'
351 if (prm%atol_rho < 0.0_preal) extmsg = trim(extmsg)//
' atol_rho'
352 if (prm%CFLfactor < 0.0_preal) extmsg = trim(extmsg)//
' CFLfactor'
354 if (prm%p <= 0.0_preal .or. prm%p > 1.0_preal) extmsg = trim(extmsg)//
' p'
355 if (prm%q < 1.0_preal .or. prm%q > 2.0_preal) extmsg = trim(extmsg)//
' q'
357 if (prm%linetensionEffect < 0.0_preal .or. prm%linetensionEffect > 1.0_preal) &
358 extmsg = trim(extmsg)//
' linetensionEffect'
359 if (prm%edgeJogFactor < 0.0_preal .or. prm%edgeJogFactor > 1.0_preal) &
360 extmsg = trim(extmsg)//
' edgeJogFactor'
362 if (prm%solidSolutionEnergy <= 0.0_preal) extmsg = trim(extmsg)//
' solidSolutionEnergy'
363 if (prm%solidSolutionSize <= 0.0_preal) extmsg = trim(extmsg)//
' solidSolutionSize'
364 if (prm%solidSolutionConcentration <= 0.0_preal) extmsg = trim(extmsg)//
' solidSolutionConcentration'
366 if (prm%grainboundaryTransmissivity > 1.0_preal) extmsg = trim(extmsg)//
' grainboundaryTransmissivity'
367 if (prm%surfaceTransmissivity < 0.0_preal .or. prm%surfaceTransmissivity > 1.0_preal) &
368 extmsg = trim(extmsg)//
' surfaceTransmissivity'
370 if (prm%fEdgeMultiplication < 0.0_preal .or. prm%fEdgeMultiplication > 1.0_preal) &
371 extmsg = trim(extmsg)//
' fEdgeMultiplication'
377 nipcmyphase = count(material_phaseat==p) * discretization_nip
378 sizedotstate =
size([
'rhoSglEdgePosMobile ',
'rhoSglEdgeNegMobile ', &
379 'rhoSglScrewPosMobile ',
'rhoSglScrewNegMobile ', &
380 'rhoSglEdgePosImmobile ',
'rhoSglEdgeNegImmobile ', &
381 'rhoSglScrewPosImmobile',
'rhoSglScrewNegImmobile', &
382 'rhoDipEdge ',
'rhoDipScrew ', &
383 'gamma ' ]) * prm%sum_N_sl
384 sizedependentstate =
size([
'rhoForest ']) * prm%sum_N_sl
385 sizestate = sizedotstate + sizedependentstate &
386 +
size([
'velocityEdgePos ',
'velocityEdgeNeg ', &
387 'velocityScrewPos ',
'velocityScrewNeg ', &
388 'maxDipoleHeightEdge ',
'maxDipoleHeightScrew' ]) * prm%sum_N_sl
389 sizedeltastate = sizedotstate
391 call material_allocateplasticstate(p,nipcmyphase,sizestate,sizedotstate,sizedeltastate)
393 plasticstate(p)%nonlocal = .true.
394 plasticstate(p)%offsetDeltaState = 0
396 st0%rho => plasticstate(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
397 stt%rho => plasticstate(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
398 dot%rho => plasticstate(p)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
399 del%rho => plasticstate(p)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
400 plasticstate(p)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
402 stt%rhoSgl => plasticstate(p)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
403 dot%rhoSgl => plasticstate(p)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
404 del%rhoSgl => plasticstate(p)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
406 stt%rhoSglMobile => plasticstate(p)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
407 dot%rhoSglMobile => plasticstate(p)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
408 del%rhoSglMobile => plasticstate(p)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
410 stt%rho_sgl_mob_edg_pos => plasticstate(p)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
411 dot%rho_sgl_mob_edg_pos => plasticstate(p)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
412 del%rho_sgl_mob_edg_pos => plasticstate(p)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
414 stt%rho_sgl_mob_edg_neg => plasticstate(p)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
415 dot%rho_sgl_mob_edg_neg => plasticstate(p)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
416 del%rho_sgl_mob_edg_neg => plasticstate(p)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
418 stt%rho_sgl_mob_scr_pos => plasticstate(p)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
419 dot%rho_sgl_mob_scr_pos => plasticstate(p)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
420 del%rho_sgl_mob_scr_pos => plasticstate(p)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
422 stt%rho_sgl_mob_scr_neg => plasticstate(p)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
423 dot%rho_sgl_mob_scr_neg => plasticstate(p)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
424 del%rho_sgl_mob_scr_neg => plasticstate(p)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
426 stt%rhoSglImmobile => plasticstate(p)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
427 dot%rhoSglImmobile => plasticstate(p)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
428 del%rhoSglImmobile => plasticstate(p)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
430 stt%rho_sgl_imm_edg_pos => plasticstate(p)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
431 dot%rho_sgl_imm_edg_pos => plasticstate(p)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
432 del%rho_sgl_imm_edg_pos => plasticstate(p)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
434 stt%rho_sgl_imm_edg_neg => plasticstate(p)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
435 dot%rho_sgl_imm_edg_neg => plasticstate(p)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
436 del%rho_sgl_imm_edg_neg => plasticstate(p)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
438 stt%rho_sgl_imm_scr_pos => plasticstate(p)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
439 dot%rho_sgl_imm_scr_pos => plasticstate(p)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
440 del%rho_sgl_imm_scr_pos => plasticstate(p)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
442 stt%rho_sgl_imm_scr_neg => plasticstate(p)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
443 dot%rho_sgl_imm_scr_neg => plasticstate(p)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
444 del%rho_sgl_imm_scr_neg => plasticstate(p)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
446 stt%rhoDip => plasticstate(p)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
447 dot%rhoDip => plasticstate(p)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
448 del%rhoDip => plasticstate(p)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
450 stt%rho_dip_edg => plasticstate(p)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
451 dot%rho_dip_edg => plasticstate(p)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
452 del%rho_dip_edg => plasticstate(p)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
454 stt%rho_dip_scr => plasticstate(p)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
455 dot%rho_dip_scr => plasticstate(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
456 del%rho_dip_scr => plasticstate(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
458 stt%gamma => plasticstate(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
459 dot%gamma => plasticstate(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
460 del%gamma => plasticstate(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
461 plasticstate(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) =
config%getFloat(
'atol_gamma', defaultval = 1.0e-2_preal)
462 if(any(plasticstate(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_preal)) &
463 extmsg = trim(extmsg)//
' atol_gamma'
464 plasticstate(p)%slipRate => plasticstate(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
466 stt%rho_forest => plasticstate(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:nipcmyphase)
467 stt%v => plasticstate(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:nipcmyphase)
468 stt%v_edg_pos => plasticstate(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:nipcmyphase)
469 stt%v_edg_neg => plasticstate(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:nipcmyphase)
470 stt%v_scr_pos => plasticstate(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:nipcmyphase)
471 stt%v_scr_neg => plasticstate(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:nipcmyphase)
473 allocate(dst%tau_pass(prm%sum_N_sl,nipcmyphase),source=0.0_preal)
474 allocate(dst%tau_back(prm%sum_N_sl,nipcmyphase),source=0.0_preal)
477 if (nipcmyphase > 0)
call stateinit(ini,p,nipcmyphase)
478 plasticstate(p)%state0 = plasticstate(p)%state
482 if (extmsg /=
'')
call io_error(211,ext_msg=trim(extmsg)//
'('//plasticity_nonlocal_label//
')')
486 allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nipneighbors,&
487 discretization_nip,discretization_nelem), source=0.0_preal)
490 allocate(irhou(maxval(param%sum_N_sl),4,ninstance), source=0)
491 allocate(iv(maxval(param%sum_N_sl),4,ninstance), source=0)
492 allocate(id(maxval(param%sum_N_sl),2,ninstance), source=0)
494 initializeinstances:
do p = 1,
size(phase_plasticity)
495 nipcmyphase = count(material_phaseat==p) * discretization_nip
496 myphase2:
if (phase_plasticity(p) == plasticity_nonlocal_id)
then
499 do s = 1,param(phase_plasticityinstance(p))%sum_N_sl
501 irhou(s,t,phase_plasticityinstance(p)) = l
504 l = l + (4+2+1+1)*param(phase_plasticityinstance(p))%sum_N_sl
506 do s = 1,param(phase_plasticityinstance(p))%sum_N_sl
508 iv(s,t,phase_plasticityinstance(p)) = l
512 do s = 1,param(phase_plasticityinstance(p))%sum_N_sl
514 id(s,t,phase_plasticityinstance(p)) = l
517 if (id(param(phase_plasticityinstance(p))%sum_N_sl,2,phase_plasticityinstance(p)) /= plasticstate(p)%sizeState) &
518 call io_error(0, ext_msg =
'state indices not properly set ('//plasticity_nonlocal_label//
')')
520 enddo initializeinstances
522 end subroutine plastic_nonlocal_init
528 module subroutine plastic_nonlocal_dependentstate(f, fp, instance, of, ip, el)
530 real(pReal),
dimension(3,3),
intent(in) :: &
533 integer,
intent(in) :: &
540 no, & !< neighbor offset
551 integer,
dimension(2) :: &
553 real(pReal),
dimension(2) :: &
555 rhoExcessGradient_over_rho, &
557 real(pReal),
dimension(3) :: &
558 rhoExcessDifferences, &
560 real(pReal),
dimension(3,3) :: &
561 invFe, & !< inverse of elastic deformation gradient
562 invFp, & !< inverse of plastic deformation gradient
565 real(pReal),
dimension(3,nIPneighbors) :: &
566 connection_latticeConf
567 real(pReal),
dimension(2,param(instance)%sum_N_sl) :: &
569 real(pReal),
dimension(param(instance)%sum_N_sl) :: &
572 real(pReal),
dimension(param(instance)%sum_N_sl,10) :: &
576 real(pReal),
dimension(param(instance)%sum_N_sl,param(instance)%sum_N_sl) :: &
578 real(pReal),
dimension(param(instance)%sum_N_sl,nIPneighbors) :: &
579 rho_edg_delta_neighbor, &
580 rho_scr_delta_neighbor
581 real(pReal),
dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: &
582 neighbor_rhoExcess, &
584 real(pReal),
dimension(3,param(instance)%sum_N_sl,2) :: &
587 associate(prm => param(instance),dst => microstructure(instance), stt => state(instance))
589 rho = getrho(instance,of,ip,el)
591 stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
592 + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
597 if (any(lattice_structure(material_phaseat(1,el)) == [lattice_bcc_id,lattice_fcc_id]))
then
598 myinteractionmatrix = prm%interactionSlipSlip &
599 * spread(( 1.0_preal - prm%linetensionEffect &
600 + prm%linetensionEffect &
601 * log(0.35_preal * prm%burgers * sqrt(max(stt%rho_forest(:,of),prm%significantRho))) &
602 / log(0.35_preal * prm%burgers * 1e6_preal))** 2.0_preal,2,prm%sum_N_sl)
604 myinteractionmatrix = prm%interactionSlipSlip
607 dst%tau_pass(:,of) = prm%mu * prm%burgers &
608 * sqrt(matmul(myinteractionmatrix,sum(abs(rho),2)))
617 rho0 = getrho0(instance,of,ip,el)
618 if (.not. phase_localplasticity(material_phaseat(1,el)) .and. prm%shortRangeStressCorrection)
then
619 invfp = math_inv33(fp)
620 invfe = matmul(fp,math_inv33(f))
622 rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
623 rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
625 rhoexcess(1,:) = rho_edg_delta
626 rhoexcess(2,:) = rho_scr_delta
628 fvsize = ipvolume(ip,el) ** (1.0_preal/3.0_preal)
632 nrealneighbors = 0.0_preal
633 neighbor_rhototal = 0.0_preal
634 do n = 1,nipneighbors
635 neighbor_el = ipneighborhood(1,n,ip,el)
636 neighbor_ip = ipneighborhood(2,n,ip,el)
637 no = material_phasememberat(1,neighbor_ip,neighbor_el)
638 if (neighbor_el > 0 .and. neighbor_ip > 0)
then
639 neighbor_instance = phase_plasticityinstance(material_phaseat(1,neighbor_el))
640 if (neighbor_instance == instance)
then
642 nrealneighbors = nrealneighbors + 1.0_preal
643 rho_neighbor0 = getrho0(instance,no,neighbor_ip,neighbor_el)
645 rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
646 rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
648 neighbor_rhototal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
649 neighbor_rhototal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
651 connection_latticeconf(1:3,n) = matmul(invfe, discretization_ipcoords(1:3,neighbor_el+neighbor_ip-1) &
652 - discretization_ipcoords(1:3,el+neighbor_ip-1))
653 normal_latticeconf = matmul(transpose(invfp), ipareanormal(1:3,n,ip,el))
654 if (math_inner(normal_latticeconf,connection_latticeconf(1:3,n)) < 0.0_preal) &
655 connection_latticeconf(1:3,n) = normal_latticeconf * ipvolume(ip,el)/iparea(n,ip,el)
658 connection_latticeconf(1:3,n) = 0.0_preal
659 rho_edg_delta_neighbor(:,n) = rho_edg_delta
660 rho_scr_delta_neighbor(:,n) = rho_scr_delta
664 connection_latticeconf(1:3,n) = 0.0_preal
665 rho_edg_delta_neighbor(:,n) = rho_edg_delta
666 rho_scr_delta_neighbor(:,n) = rho_scr_delta
670 neighbor_rhoexcess(1,:,:) = rho_edg_delta_neighbor
671 neighbor_rhoexcess(2,:,:) = rho_scr_delta_neighbor
676 m(1:3,:,1) = prm%slip_direction
677 m(1:3,:,2) = -prm%slip_transverse
679 do s = 1,prm%sum_N_sl
684 neighbors(1) = 2 * dir - 1
685 neighbors(2) = 2 * dir
686 connections(dir,1:3) = connection_latticeconf(1:3,neighbors(1)) &
687 - connection_latticeconf(1:3,neighbors(2))
688 rhoexcessdifferences(dir) = neighbor_rhoexcess(c,s,neighbors(1)) &
689 - neighbor_rhoexcess(c,s,neighbors(2))
691 invconnections = math_inv33(connections)
692 if (all(deq0(invconnections)))
call io_error(-1,ext_msg=
'back stress calculation: inversion error')
694 rhoexcessgradient(c) = math_inner(m(1:3,s,c), matmul(invconnections,rhoexcessdifferences))
698 rhoexcessgradient(1) = rhoexcessgradient(1) + sum(rho(s,imm_edg)) / fvsize
699 rhoexcessgradient(2) = rhoexcessgradient(2) + sum(rho(s,imm_scr)) / fvsize
702 rhototal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhototal(1,s,:))) / (1.0_preal + nrealneighbors)
703 rhototal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhototal(2,s,:))) / (1.0_preal + nrealneighbors)
705 rhoexcessgradient_over_rho = 0.0_preal
706 where(rhototal > 0.0_preal) rhoexcessgradient_over_rho = rhoexcessgradient / rhototal
709 dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_preal * pi) &
710 * ( rhoexcessgradient_over_rho(1) / (1.0_preal - prm%nu) &
711 + rhoexcessgradient_over_rho(2))
715 # 721 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
719 end subroutine plastic_nonlocal_dependentstate
725 module subroutine plastic_nonlocal_lpanditstangent(lp,dlp_dmp, &
726 mp,temperature,instance,of,ip,el)
727 real(pReal),
dimension(3,3),
intent(out) :: &
729 real(pReal),
dimension(3,3,3,3),
intent(out) :: &
731 integer,
intent(in) :: &
734 ip, & !< current integration point
736 real(pReal),
intent(in) :: &
739 real(pReal),
dimension(3,3),
intent(in) :: &
743 ns, & !< short notation for the total number of active slip systems
748 t, & !< dislocation type
750 real(pReal),
dimension(param(instance)%sum_N_sl,8) :: &
752 real(pReal),
dimension(param(instance)%sum_N_sl,10) :: &
754 real(pReal),
dimension(param(instance)%sum_N_sl,4) :: &
756 tauNS, & !< resolved shear stress including non Schmid and backstress terms
757 dv_dtau, & !< velocity derivative with respect to the shear stress
759 real(pReal),
dimension(param(instance)%sum_N_sl) :: &
760 tau, & !< resolved shear stress including backstress terms
763 associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance))
767 rho = getrho(instance,of,ip,el)
771 tau(s) = math_tensordot(mp, prm%Schmid(1:3,1:3,s))
774 if (tau(s) > 0.0_preal)
then
775 tauns(s,3) = math_tensordot(mp, +prm%nonSchmid_pos(1:3,1:3,s))
776 tauns(s,4) = math_tensordot(mp, -prm%nonSchmid_neg(1:3,1:3,s))
778 tauns(s,3) = math_tensordot(mp, +prm%nonSchmid_neg(1:3,1:3,s))
779 tauns(s,4) = math_tensordot(mp, -prm%nonSchmid_pos(1:3,1:3,s))
782 tauns = tauns + spread(dst%tau_back(:,of),2,4)
783 tau = tau + dst%tau_back(:,of)
786 call kinetics(v(:,1), dv_dtau(:,1), dv_dtauns(:,1), &
787 tau, tauns(:,1), dst%tau_pass(:,of),1,temperature, instance)
789 dv_dtau(:,2) = dv_dtau(:,1)
790 dv_dtauns(:,2) = dv_dtauns(:,1)
793 if (prm%nonSchmidActive)
then
794 v(:,3:4) = spread(v(:,1),2,2)
795 dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2)
796 dv_dtauns(:,3:4) = spread(dv_dtauns(:,1),2,2)
799 call kinetics(v(:,t), dv_dtau(:,t), dv_dtauns(:,t), &
800 tau, tauns(:,t), dst%tau_pass(:,of),2,temperature, instance)
804 stt%v(:,of) = pack(v,.true.)
807 forall (s = 1:ns, t = 5:8, rhosgl(s,t) * v(s,t-4) < 0.0_preal) &
808 rhosgl(s,t-4) = rhosgl(s,t-4) + abs(rhosgl(s,t))
810 gdottotal = sum(rhosgl(:,1:4) * v, 2) * prm%burgers
815 lp = lp + gdottotal(s) * prm%Schmid(1:3,1:3,s)
816 forall (i=1:3,j=1:3,k=1:3,l=1:3) &
817 dlp_dmp(i,j,k,l) = dlp_dmp(i,j,k,l) &
818 + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
819 * sum(rhosgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) &
820 + prm%Schmid(i,j,s) &
821 * ( prm%nonSchmid_pos(k,l,s) * rhosgl(s,3) * dv_dtauns(s,3) &
822 - prm%nonSchmid_neg(k,l,s) * rhosgl(s,4) * dv_dtauns(s,4)) * prm%burgers(s)
827 end subroutine plastic_nonlocal_lpanditstangent
833 module subroutine plastic_nonlocal_deltastate(mp,instance,of,ip,el)
835 real(pReal),
dimension(3,3),
intent(in) :: &
837 integer,
intent(in) :: &
849 real(pReal),
dimension(param(instance)%sum_N_sl,10) :: &
850 deltaRhoRemobilization, &
851 deltaRhoDipole2SingleStress
852 real(pReal),
dimension(param(instance)%sum_N_sl,10) :: &
854 real(pReal),
dimension(param(instance)%sum_N_sl,4) :: &
856 real(pReal),
dimension(param(instance)%sum_N_sl) :: &
858 real(pReal),
dimension(param(instance)%sum_N_sl,2) :: &
864 ph = material_phaseat(1,el)
866 associate(prm => param(instance),dst => microstructure(instance),del => deltastate(instance))
870 forall (s = 1:ns, t = 1:4) v(s,t) = plasticstate(ph)%state(iv(s,t,instance),of)
871 forall (s = 1:ns, c = 1:2) dupperold(s,c) = plasticstate(ph)%state(id(s,c,instance),of)
873 rho = getrho(instance,of,ip,el)
878 where(rho(:,imm) * v < 0.0_preal)
879 deltarhoremobilization(:,mob) = abs(rho(:,imm))
880 deltarhoremobilization(:,imm) = - rho(:,imm)
881 rho(:,mob) = rho(:,mob) + abs(rho(:,imm))
882 rho(:,imm) = 0.0_preal
884 deltarhoremobilization(:,mob) = 0.0_preal
885 deltarhoremobilization(:,imm) = 0.0_preal
887 deltarhoremobilization(:,dip) = 0.0_preal
893 do s = 1,prm%sum_N_sl
894 tau(s) = math_tensordot(mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,of)
895 if (abs(tau(s)) < 1.0e-15_preal) tau(s) = 1.0e-15_preal
898 dupper(:,1) = prm%mu * prm%burgers/(8.0_preal * pi * (1.0_preal - prm%nu) * abs(tau))
899 dupper(:,2) = prm%mu * prm%burgers/(4.0_preal * pi * abs(tau))
901 where(dneq0(sqrt(sum(abs(rho(:,edg)),2)))) &
902 dupper(:,1) = min(1.0_preal/sqrt(sum(abs(rho(:,edg)),2)),dupper(:,1))
903 where(dneq0(sqrt(sum(abs(rho(:,scr)),2)))) &
904 dupper(:,2) = min(1.0_preal/sqrt(sum(abs(rho(:,scr)),2)),dupper(:,2))
906 dupper = max(dupper,prm%minDipoleHeight)
907 deltadupper = dupper - dupperold
911 deltarhodipole2singlestress = 0.0_preal
912 forall (c=1:2, s=1:ns, deltadupper(s,c) < 0.0_preal .and. &
913 dneq0(dupperold(s,c) - prm%minDipoleHeight(s,c))) &
914 deltarhodipole2singlestress(s,8+c) = rhodip(s,c) * deltadupper(s,c) &
915 / (dupperold(s,c) - prm%minDipoleHeight(s,c))
917 forall (t=1:4) deltarhodipole2singlestress(:,t) = -0.5_preal * deltarhodipole2singlestress(:,(t-1)/2+9)
918 forall (s = 1:ns, c = 1:2) plasticstate(ph)%state(id(s,c,instance),of) = dupper(s,c)
920 plasticstate(ph)%deltaState(:,of) = 0.0_preal
921 del%rho(:,of) = reshape(deltarhoremobilization + deltarhodipole2singlestress, [10*ns])
923 # 936 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
927 end subroutine plastic_nonlocal_deltastate
933 module subroutine plastic_nonlocal_dotstate(mp, f, fp, temperature,timestep, &
936 real(pReal),
dimension(3,3),
intent(in) :: &
938 real(pReal),
dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem),
intent(in) :: &
939 F, & !< elastic deformation gradient
941 real(pReal),
intent(in) :: &
942 Temperature, & !< temperature
944 integer,
intent(in) :: &
947 ip, & !< current integration point
952 neighbor_instance, & !< instance of my neighbor
's plasticity
953 ns, & !< short notation for the total number of active slip systems
954 c, & !< character of dislocation
955 n, & !< index of my current neighbor
956 neighbor_el, & !< element number of my neighbor
957 neighbor_ip, & !< integration point of my neighbor
958 neighbor_n, & !< neighbor index pointing to me when looking from my neighbor
959 opposite_neighbor, & !< index of my opposite neighbor
960 opposite_ip, & !< ip of my opposite neighbor
961 opposite_el, & !< element index of my opposite neighbor
962 opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor
963 t, & !< type of dislocation
964 no,& !< neighbor offset shortcut
965 np,& !< neighbor phase shortcut
966 topp, & !< type of dislocation with opposite sign to t
967 s !< index of my current slip system
968 real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
970 rho0, & !< dislocation density at beginning of time step
971 rhoDot, & !< density evolution
972 rhoDotMultiplication, & !< density evolution by multiplication
973 rhoDotFlux, & !< density evolution by flux
974 rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
975 rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
976 rhoDotThermalAnnihilation !< density evolution by thermal annihilation
977 real(pReal), dimension(param(instance)%sum_N_sl,8) :: &
978 rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
979 neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
980 my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
981 real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
982 v, & !< current dislocation glide velocity
984 neighbor_v0, & !< dislocation glide velocity of enighboring ip
986 real(pReal), dimension(param(instance)%sum_N_sl) :: &
987 tau, & !< current resolved shear stress
988 vClimb !< climb velocity of edge dipoles
989 real(pReal), dimension(param(instance)%sum_N_sl,2) :: &
990 rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
991 dLower, & !< minimum stable dipole distance for edges and screws
992 dUpper !< current maximum stable dipole distance for edges and screws
993 real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: &
994 m !< direction of dislocation motion
995 real(pReal), dimension(3,3) :: &
996 my_F, & !< my total deformation gradient
997 neighbor_F, & !< total deformation gradient of my neighbor
998 my_Fe, & !< my elastic deformation gradient
999 neighbor_Fe, & !< elastic deformation gradient of my neighbor
1000 Favg !< average total deformation gradient of me and my neighbor
1001 real(pReal), dimension(3) :: &
1002 normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration
1003 normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration
1004 normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration
1005 normal_me2neighbor_defConf
1007 area, & !< area of the current interface
1008 transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point
1009 lineLength, & !< dislocation line length leaving the current interface
1012 ph = material_phaseat(1,el)
1013 if (timestep <= 0.0_preal)
then
1014 plasticstate(ph)%dotState = 0.0_preal
1018 associate(prm => param(instance), &
1019 dst => microstructure(instance), &
1020 dot => dotstate(instance), &
1021 stt => state(instance))
1027 rho = getrho(instance,of,ip,el)
1030 rho0 = getrho0(instance,of,ip,el)
1031 my_rhosgl0 = rho0(:,sgl)
1033 forall (s = 1:ns, t = 1:4) v(s,t) = plasticstate(ph)%state(iv(s,t,instance),of)
1034 gdot = rhosgl(:,1:4) * v * spread(prm%burgers,2,4)
1036 # 1056 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
1041 tau(s) = math_tensordot(mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,of)
1042 if (abs(tau(s)) < 1.0e-15_preal) tau(s) = 1.0e-15_preal
1045 dlower = prm%minDipoleHeight
1046 dupper(:,1) = prm%mu * prm%burgers/(8.0_preal * pi * (1.0_preal - prm%nu) * abs(tau))
1047 dupper(:,2) = prm%mu * prm%burgers/(4.0_preal * pi * abs(tau))
1049 where(dneq0(sqrt(sum(abs(rho(:,edg)),2)))) &
1050 dupper(:,1) = min(1.0_preal/sqrt(sum(abs(rho(:,edg)),2)),dupper(:,1))
1051 where(dneq0(sqrt(sum(abs(rho(:,scr)),2)))) &
1052 dupper(:,2) = min(1.0_preal/sqrt(sum(abs(rho(:,scr)),2)),dupper(:,2))
1054 dupper = max(dupper,dlower)
1058 rhodotmultiplication = 0.0_preal
1059 isbcc:
if (lattice_structure(ph) == lattice_bcc_id)
then
1060 forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_preal)
1061 rhodotmultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) &
1062 * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s)
1064 rhodotmultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) &
1065 * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s)
1070 rhodotmultiplication(:,1:4) = spread( &
1071 (sum(abs(gdot(:,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) &
1072 * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4)
1075 forall (s = 1:ns, t = 1:4) v0(s,t) = plasticstate(ph)%state0(iv(s,t,instance),of)
1079 rhodotflux = 0.0_preal
1080 if (.not. phase_localplasticity(material_phaseat(1,el)))
then
1083 if (any( abs(gdot) > 0.0_preal &
1084 .and. prm%CFLfactor * abs(v0) * timestep &
1085 > ipvolume(ip,el) / maxval(iparea(:,ip,el))))
then
1086 # 1116 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
1087 plasticstate(ph)%dotState = ieee_value(1.0_preal,ieee_quiet_nan)
1095 m(1:3,:,1) = prm%slip_direction
1096 m(1:3,:,2) = -prm%slip_direction
1097 m(1:3,:,3) = -prm%slip_transverse
1098 m(1:3,:,4) = prm%slip_transverse
1100 my_f = f(1:3,1:3,1,ip,el)
1101 my_fe = matmul(my_f, math_inv33(fp(1:3,1:3,1,ip,el)))
1103 neighbors:
do n = 1,nipneighbors
1105 neighbor_el = ipneighborhood(1,n,ip,el)
1106 neighbor_ip = ipneighborhood(2,n,ip,el)
1107 neighbor_n = ipneighborhood(3,n,ip,el)
1108 np = material_phaseat(1,neighbor_el)
1109 no = material_phasememberat(1,neighbor_ip,neighbor_el)
1111 opposite_neighbor = n + mod(n,2) - mod(n+1,2)
1112 opposite_el = ipneighborhood(1,opposite_neighbor,ip,el)
1113 opposite_ip = ipneighborhood(2,opposite_neighbor,ip,el)
1114 opposite_n = ipneighborhood(3,opposite_neighbor,ip,el)
1116 if (neighbor_n > 0)
then
1117 neighbor_instance = phase_plasticityinstance(material_phaseat(1,neighbor_el))
1118 neighbor_f = f(1:3,1:3,1,neighbor_ip,neighbor_el)
1119 neighbor_fe = matmul(neighbor_f, math_inv33(fp(1:3,1:3,1,neighbor_ip,neighbor_el)))
1120 favg = 0.5_preal * (my_f + neighbor_f)
1125 neighbor_v0 = 0.0_preal
1135 if (neighbor_n > 0)
then
1136 if (phase_plasticity(material_phaseat(1,neighbor_el)) == plasticity_nonlocal_id .and. &
1137 any(compatibility(:,:,:,n,ip,el) > 0.0_preal))
then
1139 forall (s = 1:ns, t = 1:4)
1140 neighbor_v0(s,t) = plasticstate(np)%state0(iv(s,t,neighbor_instance),no)
1141 neighbor_rhosgl0(s,t) = max(plasticstate(np)%state0(irhou(s,t,neighbor_instance),no),0.0_preal)
1144 where (neighbor_rhosgl0 * ipvolume(neighbor_ip,neighbor_el) ** 0.667_preal < prm%significantN &
1145 .or. neighbor_rhosgl0 < prm%significantRho) &
1146 neighbor_rhosgl0 = 0.0_preal
1147 normal_neighbor2me_defconf = math_det33(favg) * matmul(math_inv33(transpose(favg)), &
1148 ipareanormal(1:3,neighbor_n,neighbor_ip,neighbor_el))
1149 normal_neighbor2me = matmul(transpose(neighbor_fe), normal_neighbor2me_defconf) &
1150 / math_det33(neighbor_fe)
1151 area = iparea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
1152 normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me)
1156 topp = t + mod(t,2) - mod(t+1,2)
1157 if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_preal &
1158 .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_preal )
then
1159 linelength = neighbor_rhosgl0(s,t) * neighbor_v0(s,t) &
1160 * math_inner(m(1:3,s,t), normal_neighbor2me) * area
1161 where (compatibility(c,:,s,n,ip,el) > 0.0_preal) &
1162 rhodotflux(:,t) = rhodotflux(1:ns,t) &
1163 + linelength/ipvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_preal
1164 where (compatibility(c,:,s,n,ip,el) < 0.0_preal) &
1165 rhodotflux(:,topp) = rhodotflux(:,topp) &
1166 + linelength/ipvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_preal
1181 if (opposite_n > 0)
then
1182 if (phase_plasticity(material_phaseat(1,opposite_el)) == plasticity_nonlocal_id)
then
1184 normal_me2neighbor_defconf = math_det33(favg) &
1185 * matmul(math_inv33(transpose(favg)),ipareanormal(1:3,n,ip,el))
1186 normal_me2neighbor = matmul(transpose(my_fe), normal_me2neighbor_defconf) &
1188 area = iparea(n,ip,el) * norm2(normal_me2neighbor)
1189 normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor)
1193 if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_preal )
then
1194 if (v0(s,t) * neighbor_v0(s,t) >= 0.0_preal)
then
1195 transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_preal)
1197 transmissivity = 0.0_preal
1199 linelength = my_rhosgl0(s,t) * v0(s,t) &
1200 * math_inner(m(1:3,s,t), normal_me2neighbor) * area
1201 rhodotflux(s,t) = rhodotflux(s,t) - linelength / ipvolume(ip,el)
1202 rhodotflux(s,t+4) = rhodotflux(s,t+4) &
1203 + linelength / ipvolume(ip,el) * (1.0_preal - transmissivity) &
1204 * sign(1.0_preal, v0(s,t))
1220 rhodotsingle2dipoleglide(:,2*c-1) = -2.0_preal * dupper(:,c) / prm%burgers &
1221 * ( rhosgl(:,2*c-1) * abs(gdot(:,2*c)) &
1222 + rhosgl(:,2*c) * abs(gdot(:,2*c-1)) &
1223 + abs(rhosgl(:,2*c+4)) * abs(gdot(:,2*c-1)))
1225 rhodotsingle2dipoleglide(:,2*c) = -2.0_preal * dupper(:,c) / prm%burgers &
1226 * ( rhosgl(:,2*c-1) * abs(gdot(:,2*c)) &
1227 + rhosgl(:,2*c) * abs(gdot(:,2*c-1)) &
1228 + abs(rhosgl(:,2*c+3)) * abs(gdot(:,2*c)))
1230 rhodotsingle2dipoleglide(:,2*c+3) = -2.0_preal * dupper(:,c) / prm%burgers &
1231 * rhosgl(:,2*c+3) * abs(gdot(:,2*c))
1233 rhodotsingle2dipoleglide(:,2*c+4) = -2.0_preal * dupper(:,c) / prm%burgers &
1234 * rhosgl(:,2*c+4) * abs(gdot(:,2*c-1))
1236 rhodotsingle2dipoleglide(:,c+8) = abs(rhodotsingle2dipoleglide(:,2*c+3)) &
1237 + abs(rhodotsingle2dipoleglide(:,2*c+4)) &
1238 - rhodotsingle2dipoleglide(:,2*c-1) &
1239 - rhodotsingle2dipoleglide(:,2*c)
1244 rhodotathermalannihilation = 0.0_preal
1246 rhodotathermalannihilation(:,c+8) = -2.0_preal * dlower(:,c) / prm%burgers &
1247 * ( 2.0_preal * (rhosgl(:,2*c-1) * abs(gdot(:,2*c)) + rhosgl(:,2*c) * abs(gdot(:,2*c-1))) &
1248 + 2.0_preal * (abs(rhosgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhosgl(:,2*c+4)) * abs(gdot(:,2*c-1))) &
1249 + rhodip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c))))
1252 if (lattice_structure(ph) == lattice_fcc_id) &
1253 forall (s = 1:ns, prm%colinearSystem(s) > 0) &
1254 rhodotathermalannihilation(prm%colinearSystem(s),1:2) = - rhodotathermalannihilation(s,10) &
1255 * 0.25_preal * sqrt(stt%rho_forest(s,of)) * (dupper(s,2) + dlower(s,2)) * prm%edgeJogFactor
1259 rhodotthermalannihilation = 0.0_preal
1260 selfdiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kb * temperature))
1261 vclimb = prm%atomicVolume * selfdiffusion * prm%mu &
1262 / ( kb * temperature * pi * (1.0_preal-prm%nu) * (dupper(:,1) + dlower(:,1)))
1263 forall (s = 1:ns, dupper(s,1) > dlower(s,1)) &
1264 rhodotthermalannihilation(s,9) = max(- 4.0_preal * rhodip(s,1) * vclimb(s) / (dupper(s,1) - dlower(s,1)), &
1265 - rhodip(s,1) / timestep - rhodotathermalannihilation(s,9) &
1266 - rhodotsingle2dipoleglide(s,9))
1268 rhodot = rhodotflux &
1269 + rhodotmultiplication &
1270 + rhodotsingle2dipoleglide &
1271 + rhodotathermalannihilation &
1272 + rhodotthermalannihilation
1274 # 1325 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
1277 if ( any(rho(:,mob) + rhodot(:,1:4) * timestep < -prm%atol_rho) &
1278 .or. any(rho(:,dip) + rhodot(:,9:10) * timestep < -prm%atol_rho))
then
1285 plasticstate(ph)%dotState = ieee_value(1.0_preal,ieee_quiet_nan)
1287 dot%rho(:,of) = pack(rhodot,.true.)
1288 dot%gamma(:,of) = sum(gdot,2)
1293 end subroutine plastic_nonlocal_dotstate
1302 module subroutine plastic_nonlocal_updatecompatibility(orientation,instance,i,e)
1304 type(rotation),
dimension(1,discretization_nIP,discretization_nElem),
intent(in) :: &
1306 integer,
intent(in) :: &
1320 real(pReal),
dimension(2,param(instance)%sum_N_sl,param(instance)%sum_N_sl,nIPneighbors) :: &
1323 my_compatibilitySum, &
1326 logical,
dimension(param(instance)%sum_N_sl) :: &
1328 type(rotation) :: mis
1330 ph = material_phaseat(1,e)
1332 associate(prm => param(instance))
1336 my_compatibility = 0.0_preal
1337 forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_preal
1339 neighbors:
do n = 1,nipneighbors
1340 neighbor_e = ipneighborhood(1,n,i,e)
1341 neighbor_i = ipneighborhood(2,n,i,e)
1343 neighbor_phase = material_phaseat(1,neighbor_e)
1345 if (neighbor_e <= 0 .or. neighbor_i <= 0)
then
1348 forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%surfaceTransmissivity)
1349 elseif (neighbor_phase /= ph)
then
1355 if (.not. phase_localplasticity(neighbor_phase) .and. .not. phase_localplasticity(ph)) &
1356 forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_preal
1357 elseif (prm%grainboundaryTransmissivity >= 0.0_preal)
then
1360 if (material_texture(1,i,e) /= material_texture(1,neighbor_i,neighbor_e) .and. &
1361 (.not. phase_localplasticity(neighbor_phase))) &
1362 forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%grainboundaryTransmissivity)
1373 mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e))
1374 myslipsystems:
do s1 = 1,ns
1375 neighborslipsystems:
do s2 = 1,ns
1376 my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
1377 mis%rotate(prm%slip_normal(1:3,s2))) &
1378 * abs(math_inner(prm%slip_direction(1:3,s1), &
1379 mis%rotate(prm%slip_direction(1:3,s2))))
1380 my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
1381 mis%rotate(prm%slip_normal(1:3,s2)))) &
1382 * abs(math_inner(prm%slip_direction(1:3,s1), &
1383 mis%rotate(prm%slip_direction(1:3,s2))))
1384 enddo neighborslipsystems
1386 my_compatibilitysum = 0.0_preal
1387 belowthreshold = .true.
1388 do while (my_compatibilitysum < 1.0_preal .and. any(belowthreshold))
1389 thresholdvalue = maxval(my_compatibility(2,:,s1,n), belowthreshold)
1390 nthresholdvalues = real(count(my_compatibility(2,:,s1,n) >= thresholdvalue),preal)
1391 where (my_compatibility(2,:,s1,n) >= thresholdvalue) belowthreshold = .false.
1392 if (my_compatibilitysum + thresholdvalue * nthresholdvalues > 1.0_preal) &
1393 where (abs(my_compatibility(:,:,s1,n)) >= thresholdvalue) &
1394 my_compatibility(:,:,s1,n) = sign((1.0_preal - my_compatibilitysum)/nthresholdvalues,&
1395 my_compatibility(:,:,s1,n))
1396 my_compatibilitysum = my_compatibilitysum + nthresholdvalues * thresholdvalue
1399 where(belowthreshold) my_compatibility(1,:,s1,n) = 0.0_preal
1400 where(belowthreshold) my_compatibility(2,:,s1,n) = 0.0_preal
1407 compatibility(:,:,:,:,i,e) = my_compatibility
1411 end subroutine plastic_nonlocal_updatecompatibility
1417 module subroutine plastic_nonlocal_results(instance,group)
1419 integer,
intent(in) :: instance
1420 character(len=*),
intent(in) :: group
1424 associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance))
1425 outputsloop:
do o = 1,
size(prm%output)
1426 select case(trim(prm%output(o)))
1427 case(
'rho_sgl_mob_edg_pos')
1428 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_mob_edg_pos,
'rho_sgl_mob_edg_pos', &
1429 'positive mobile edge density',
'1/m²')
1430 case(
'rho_sgl_imm_edg_pos')
1431 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_imm_edg_pos,
'rho_sgl_imm_edg_pos',&
1432 'positive immobile edge density',
'1/m²')
1433 case(
'rho_sgl_mob_edg_neg')
1434 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_mob_edg_neg,
'rho_sgl_mob_edg_neg',&
1435 'negative mobile edge density',
'1/m²')
1436 case(
'rho_sgl_imm_edg_neg')
1437 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_imm_edg_neg,
'rho_sgl_imm_edg_neg',&
1438 'negative immobile edge density',
'1/m²')
1440 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_dip_edg,
'rho_dip_edg',&
1441 'edge dipole density',
'1/m²')
1442 case(
'rho_sgl_mob_scr_pos')
1443 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_mob_scr_pos,
'rho_sgl_mob_scr_pos',&
1444 'positive mobile screw density',
'1/m²')
1445 case(
'rho_sgl_imm_scr_pos')
1446 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_imm_scr_pos,
'rho_sgl_imm_scr_pos',&
1447 'positive immobile screw density',
'1/m²')
1448 case(
'rho_sgl_mob_scr_neg')
1449 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_mob_scr_neg,
'rho_sgl_mob_scr_neg',&
1450 'negative mobile screw density',
'1/m²')
1451 case(
'rho_sgl_imm_scr_neg')
1452 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_sgl_imm_scr_neg,
'rho_sgl_imm_scr_neg',&
1453 'negative immobile screw density',
'1/m²')
1455 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_dip_scr,
'rho_dip_scr',&
1456 'screw dipole density',
'1/m²')
1458 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%rho_forest,
'rho_forest',&
1459 'forest density',
'1/m²')
1461 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%v_edg_pos,
'v_edg_pos',&
1462 'positive edge velocity',
'm/s')
1464 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%v_edg_neg,
'v_edg_neg',&
1465 'negative edge velocity',
'm/s')
1467 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%v_scr_pos,
'v_scr_pos',&
1468 'positive srew velocity',
'm/s')
1470 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%v_scr_neg,
'v_scr_neg',&
1471 'negative screw velocity',
'm/s')
1473 if(prm%sum_N_sl>0)
call results_writedataset(group,stt%gamma,
'gamma',&
1474 'plastic shear',
'1')
1476 if(prm%sum_N_sl>0)
call results_writedataset(group,dst%tau_pass,
'tau_pass',&
1477 'passing stress for slip',
'Pa')
1482 end subroutine plastic_nonlocal_results
1488 subroutine stateinit(ini,phase,NipcMyPhase)
1490 type(tInitialParameters) :: &
1492 integer,
intent(in) :: &
1504 real(pReal),
dimension(2) :: &
1512 real(pReal),
dimension(NipcMyPhase) :: &
1515 instance = phase_plasticityinstance(phase)
1516 associate(stt => state(instance))
1518 if (ini%rhoSglRandom > 0.0_preal)
then
1519 do e = 1,discretization_nelem
1520 do i = 1,discretization_nip
1521 if (material_phaseat(1,e) == phase) volume(material_phasememberat(1,i,e)) = ipvolume(i,e)
1524 totalvolume = sum(volume)
1525 minimumipvolume = minval(volume)
1526 densitybinning = ini%rhoSglRandomBinning / minimumipvolume ** (2.0_preal / 3.0_preal)
1529 meandensity = 0.0_preal
1530 do while(meandensity < ini%rhoSglRandom)
1531 call random_number(rnd)
1532 phasemember = nint(rnd(1)*real(nipcmyphase,preal) + 0.5_preal)
1533 s = nint(rnd(2)*real(sum(ini%N_sl),preal)*4.0_preal + 0.5_preal)
1534 meandensity = meandensity + densitybinning * volume(phasemember) / totalvolume
1535 stt%rhoSglMobile(s,phasemember) = densitybinning
1538 do e = 1, nipcmyphase
1539 do f = 1,
size(ini%N_sl,1)
1540 from = 1 + sum(ini%N_sl(1:f-1))
1541 upto = sum(ini%N_sl(1:f))
1543 noise = [math_samplegaussvar(0.0_preal, ini%rhoSglScatter), &
1544 math_samplegaussvar(0.0_preal, ini%rhoSglScatter)]
1545 stt%rho_sgl_mob_edg_pos(s,e) = ini%rhoSglEdgePos0(f) + noise(1)
1546 stt%rho_sgl_mob_edg_neg(s,e) = ini%rhoSglEdgeNeg0(f) + noise(1)
1547 stt%rho_sgl_mob_scr_pos(s,e) = ini%rhoSglScrewPos0(f) + noise(2)
1548 stt%rho_sgl_mob_scr_neg(s,e) = ini%rhoSglScrewNeg0(f) + noise(2)
1550 stt%rho_dip_edg(from:upto,e) = ini%rhoDipEdge0(f)
1551 stt%rho_dip_scr(from:upto,e) = ini%rhoDipScrew0(f)
1558 end subroutine stateinit
1564 subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
1566 integer,
intent(in) :: &
1567 c, & !< dislocation character (1:edge, 2:screw)
1569 real(pReal),
intent(in) :: &
1571 real(pReal),
dimension(param(instance)%sum_N_sl),
intent(in) :: &
1572 tau, & !< resolved external shear stress (without non Schmid effects)
1573 tauNS, & !< resolved external shear stress (including non Schmid effects)
1575 real(pReal),
dimension(param(instance)%sum_N_sl),
intent(out) :: &
1577 dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
1581 ns, & !< short notation for the total number of active slip systems
1586 tauEff, & !< effective shear stress
1587 tPeierls, & !< waiting time in front of a peierls barriers
1588 tSolidSolution, & !< waiting time in front of a solid solution obstacle
1589 vViscous, & !< viscous glide velocity
1590 dtPeierls_dtau, & !< derivative with respect to resolved shear stress
1591 dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress
1592 meanfreepath_S, & !< mean free travel distance for dislocations between two solid solution obstacles
1593 meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers
1594 jumpWidth_P, & !< depth of activated area
1595 jumpWidth_S, & !< depth of activated area
1596 activationLength_P, & !< length of activated dislocation line
1597 activationLength_S, & !< length of activated dislocation line
1598 activationVolume_P, & !< volume that needs to be activated to overcome barrier
1599 activationVolume_S, & !< volume that needs to be activated to overcome barrier
1600 activationEnergy_P, & !< energy that is needed to overcome barrier
1601 activationEnergy_S, & !< energy that is needed to overcome barrier
1602 criticalStress_P, & !< maximum obstacle strength
1603 criticalStress_S, & !< maximum obstacle strength
1606 associate(prm => param(instance))
1610 dv_dtauns = 0.0_preal
1613 if (abs(tau(s)) > tauthreshold(s))
then
1618 taueff = max(0.0_preal, abs(tauns(s)) - tauthreshold(s))
1619 meanfreepath_p = prm%burgers(s)
1620 jumpwidth_p = prm%burgers(s)
1621 activationlength_p = prm%doublekinkwidth *prm%burgers(s)
1622 activationvolume_p = activationlength_p * jumpwidth_p * prm%burgers(s)
1623 criticalstress_p = prm%peierlsStress(s,c)
1624 activationenergy_p = criticalstress_p * activationvolume_p
1625 taurel_p = min(1.0_preal, taueff / criticalstress_p)
1626 tpeierls = 1.0_preal / prm%fattack &
1627 * exp(activationenergy_p / (kb * temperature) &
1628 * (1.0_preal - taurel_p**prm%p)**prm%q)
1629 if (taueff < criticalstress_p)
then
1630 dtpeierls_dtau = tpeierls * prm%p * prm%q * activationvolume_p / (kb * temperature) &
1631 * (1.0_preal - taurel_p**prm%p)**(prm%q-1.0_preal) * taurel_p**(prm%p-1.0_preal)
1633 dtpeierls_dtau = 0.0_preal
1638 taueff = abs(tau(s)) - tauthreshold(s)
1639 meanfreepath_s = prm%burgers(s) / sqrt(prm%solidSolutionConcentration)
1640 jumpwidth_s = prm%solidSolutionSize * prm%burgers(s)
1641 activationlength_s = prm%burgers(s) / sqrt(prm%solidSolutionConcentration)
1642 activationvolume_s = activationlength_s * jumpwidth_s * prm%burgers(s)
1643 activationenergy_s = prm%solidSolutionEnergy
1644 criticalstress_s = activationenergy_s / activationvolume_s
1645 taurel_s = min(1.0_preal, taueff / criticalstress_s)
1646 tsolidsolution = 1.0_preal / prm%fattack &
1647 * exp(activationenergy_s / (kb * temperature)* (1.0_preal - taurel_s**prm%p)**prm%q)
1648 if (taueff < criticalstress_s)
then
1649 dtsolidsolution_dtau = tsolidsolution * prm%p * prm%q * activationvolume_s / (kb * temperature) &
1650 * (1.0_preal - taurel_s**prm%p)**(prm%q-1.0_preal)* taurel_s**(prm%p-1.0_preal)
1652 dtsolidsolution_dtau = 0.0_preal
1656 taueff = abs(tau(s)) - tauthreshold(s)
1657 mobility = prm%burgers(s) / prm%viscosity
1658 vviscous = mobility * taueff
1663 v(s) = sign(1.0_preal,tau(s)) &
1664 / (tpeierls / meanfreepath_p + tsolidsolution / meanfreepath_s + 1.0_preal / vviscous)
1665 dv_dtau(s) = v(s)**2.0_preal * (dtsolidsolution_dtau / meanfreepath_s + mobility /vviscous**2.0_preal)
1666 dv_dtauns(s) = v(s)**2.0_preal * dtpeierls_dtau / meanfreepath_p
1672 end subroutine kinetics
1679 function getrho(instance,of,ip,el)
1681 integer,
intent(in) :: instance, of,ip,el
1682 real(pReal),
dimension(param(instance)%sum_N_sl,10) :: getRho
1684 associate(prm => param(instance))
1686 getrho = reshape(state(instance)%rho(:,of),[prm%sum_N_sl,10])
1689 getrho(:,mob) = max(getrho(:,mob),0.0_preal)
1690 getrho(:,dip) = max(getrho(:,dip),0.0_preal)
1692 where(abs(getrho) < max(prm%significantN/ipvolume(ip,el)**(2.0_preal/3.0_preal),prm%significantRho)) &
1704 function getrho0(instance,of,ip,el)
1706 integer,
intent(in) :: instance, of,ip,el
1707 real(pReal),
dimension(param(instance)%sum_N_sl,10) :: getRho0
1709 associate(prm => param(instance))
1711 getrho0 = reshape(state0(instance)%rho(:,of),[prm%sum_N_sl,10])
1714 getrho0(:,mob) = max(getrho0(:,mob),0.0_preal)
1715 getrho0(:,dip) = max(getrho0(:,dip),0.0_preal)
1717 where(abs(getrho0) < max(prm%significantN/ipvolume(ip,el)**(2.0_preal/3.0_preal),prm%significantRho)) &
1722 end function getrho0
1724 end submodule plastic_nonlocal