DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
constitutive_plastic_kinehardening.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_kinehardening.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_kinehardening.f90"
5 !--------------------------------------------------------------------------------------------------
11 !--------------------------------------------------------------------------------------------------
12 submodule(constitutive) plastic_kinehardening
13 
14  type :: tparameters
15  real(pReal) :: &
16  gdot0 = 1.0_preal, &
17  n = 1.0_preal
18  real(pReal), allocatable, dimension(:) :: &
19  theta0, & !< initial hardening rate of forward stress for each slip
20  theta1, & !< asymptotic hardening rate of forward stress for each slip
21  theta0_b, & !< initial hardening rate of back stress for each slip
22  theta1_b, & !< asymptotic hardening rate of back stress for each slip
23  tau1, &
24  tau1_b
25  real(pReal), allocatable, dimension(:,:) :: &
26  interaction_slipslip
27  real(pReal), allocatable, dimension(:,:,:) :: &
28  P, &
29  nonSchmid_pos, &
30  nonSchmid_neg
31  integer :: &
32  sum_N_sl, & !< total number of active slip system
33  of_debug = 0
34  logical :: &
35  nonSchmidActive = .false.
36  character(len=pStringLen), allocatable, dimension(:) :: &
37  output
38  end type tparameters
39 
40  type :: tkinehardeningstate
41  real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
42  crss, & !< critical resolved stress
43  crss_back, & !< critical resolved back stress
44  sense, & !< sense of acting shear stress (-1 or +1)
45  chi0, & !< backstress at last switch of stress sense
46  gamma0, & !< accumulated shear at last switch of stress sense
47  accshear
48  end type tkinehardeningstate
49 
50 !--------------------------------------------------------------------------------------------------
51 ! containers for parameters and state
52  type(tParameters), allocatable, dimension(:) :: param
53  type(tKinehardeningState), allocatable, dimension(:) :: &
54  dotState, &
55  deltaState, &
56  state
57 
58 contains
59 
60 
61 !--------------------------------------------------------------------------------------------------
64 !--------------------------------------------------------------------------------------------------
65 module subroutine plastic_kinehardening_init
66 
67  integer :: &
68  Ninstance, &
69  p, o, &
70  NipcMyPhase, &
71  sizeState, sizeDeltaState, sizeDotState, &
72  startIndex, endIndex
73  integer, dimension(:), allocatable :: &
74  N_sl
75  real(pReal), dimension(:), allocatable :: &
76  xi_0, & !< initial resistance against plastic flow
77  a
78  character(len=pStringLen) :: &
79  extmsg = ''
80 
81  write(6,'(/,a)') ' <<<+- plastic_'//plasticity_kinehardening_label//' init -+>>>'; flush(6)
82 
83  ninstance = count(phase_plasticity == plasticity_kinehardening_id)
84  if (iand(debug_level(debug_constitutive),debug_levelbasic) /= 0) &
85  write(6,'(a16,1x,i5,/)') '# instances:',ninstance
86 
87  allocate(param(ninstance))
88  allocate(state(ninstance))
89  allocate(dotstate(ninstance))
90  allocate(deltastate(ninstance))
91 
92  do p = 1, size(phase_plasticityinstance)
93  if (phase_plasticity(p) /= plasticity_kinehardening_id) cycle
94  associate(prm => param(phase_plasticityinstance(p)), &
95  dot => dotstate(phase_plasticityinstance(p)), &
96  dlt => deltastate(phase_plasticityinstance(p)), &
97  stt => state(phase_plasticityinstance(p)),&
98  config => config_phase(p))
99 
100  prm%output = config%getStrings('(output)',defaultval=emptystringarray)
101 
102 
103 
104 
105 
106 
107 
108 !--------------------------------------------------------------------------------------------------
109 ! slip related parameters
110  n_sl = config%getInts('nslip',defaultval=emptyintarray)
111  prm%sum_N_sl = sum(abs(n_sl))
112  slipactive: if (prm%sum_N_sl > 0) then
113  prm%P = lattice_schmidmatrix_slip(n_sl,config%getString('lattice_structure'),&
114  config%getFloat('c/a',defaultval=0.0_preal))
115 
116  if(trim(config%getString('lattice_structure')) == 'bcc') then
117  a = config%getFloats('nonschmid_coefficients',defaultval = emptyrealarray)
118  if(size(a) > 0) prm%nonSchmidActive = .true.
119  prm%nonSchmid_pos = lattice_nonschmidmatrix(n_sl,a,+1)
120  prm%nonSchmid_neg = lattice_nonschmidmatrix(n_sl,a,-1)
121  else
122  prm%nonSchmid_pos = prm%P
123  prm%nonSchmid_neg = prm%P
124  endif
125  prm%interaction_SlipSlip = lattice_interaction_slipbyslip(n_sl, &
126  config%getFloats('interaction_slipslip'), &
127  config%getString('lattice_structure'))
128 
129  xi_0 = config%getFloats('crss0', requiredsize=size(n_sl))
130  prm%tau1 = config%getFloats('tau1', requiredsize=size(n_sl))
131  prm%tau1_b = config%getFloats('tau1_b', requiredsize=size(n_sl))
132  prm%theta0 = config%getFloats('theta0', requiredsize=size(n_sl))
133  prm%theta1 = config%getFloats('theta1', requiredsize=size(n_sl))
134  prm%theta0_b = config%getFloats('theta0_b', requiredsize=size(n_sl))
135  prm%theta1_b = config%getFloats('theta1_b', requiredsize=size(n_sl))
136 
137  prm%gdot0 = config%getFloat('gdot0')
138  prm%n = config%getFloat('n_slip')
139 
140  ! expand: family => system
141  xi_0 = math_expand(xi_0, n_sl)
142  prm%tau1 = math_expand(prm%tau1, n_sl)
143  prm%tau1_b = math_expand(prm%tau1_b, n_sl)
144  prm%theta0 = math_expand(prm%theta0, n_sl)
145  prm%theta1 = math_expand(prm%theta1, n_sl)
146  prm%theta0_b = math_expand(prm%theta0_b,n_sl)
147  prm%theta1_b = math_expand(prm%theta1_b,n_sl)
148 
149 !--------------------------------------------------------------------------------------------------
150 ! sanity checks
151  if ( prm%gdot0 <= 0.0_preal) extmsg = trim(extmsg)//' gdot0'
152  if ( prm%n <= 0.0_preal) extmsg = trim(extmsg)//' n_slip'
153  if (any(xi_0 <= 0.0_preal)) extmsg = trim(extmsg)//' crss0'
154  if (any(prm%tau1 <= 0.0_preal)) extmsg = trim(extmsg)//' tau1'
155  if (any(prm%tau1_b <= 0.0_preal)) extmsg = trim(extmsg)//' tau1_b'
156 
157  !ToDo: Any sensible checks for theta?
158  else slipactive
159  xi_0 = emptyrealarray
160  allocate(prm%tau1,prm%tau1_b,prm%theta0,prm%theta1,prm%theta0_b,prm%theta1_b,source=emptyrealarray)
161  allocate(prm%interaction_SlipSlip(0,0))
162  endif slipactive
163 
164 !--------------------------------------------------------------------------------------------------
165 ! allocate state arrays
166  nipcmyphase = count(material_phaseat == p) * discretization_nip
167  sizedotstate = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl
168  sizedeltastate = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl
169  sizestate = sizedotstate + sizedeltastate
170 
171  call material_allocateplasticstate(p,nipcmyphase,sizestate,sizedotstate,sizedeltastate)
172 
173 !--------------------------------------------------------------------------------------------------
174 ! state aliases and initialization
175  startindex = 1
176  endindex = prm%sum_N_sl
177  stt%crss => plasticstate(p)%state (startindex:endindex,:)
178  stt%crss = spread(xi_0, 2, nipcmyphase)
179  dot%crss => plasticstate(p)%dotState(startindex:endindex,:)
180  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_xi',defaultval=1.0_preal)
181  if(any(plasticstate(p)%atol(startindex:endindex) < 0.0_preal)) extmsg = trim(extmsg)//' atol_xi'
182 
183  startindex = endindex + 1
184  endindex = endindex + prm%sum_N_sl
185  stt%crss_back => plasticstate(p)%state (startindex:endindex,:)
186  dot%crss_back => plasticstate(p)%dotState(startindex:endindex,:)
187  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_xi',defaultval=1.0_preal)
188 
189  startindex = endindex + 1
190  endindex = endindex + prm%sum_N_sl
191  stt%accshear => plasticstate(p)%state (startindex:endindex,:)
192  dot%accshear => plasticstate(p)%dotState(startindex:endindex,:)
193  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_gamma',defaultval=1.0e-6_preal)
194  if(any(plasticstate(p)%atol(startindex:endindex) < 0.0_preal)) extmsg = trim(extmsg)//' atol_gamma'
195  ! global alias
196  plasticstate(p)%slipRate => plasticstate(p)%dotState(startindex:endindex,:)
197 
198  o = plasticstate(p)%offsetDeltaState
199  startindex = endindex + 1
200  endindex = endindex + prm%sum_N_sl
201  stt%sense => plasticstate(p)%state (startindex :endindex ,:)
202  dlt%sense => plasticstate(p)%deltaState(startindex-o:endindex-o,:)
203 
204  startindex = endindex + 1
205  endindex = endindex + prm%sum_N_sl
206  stt%chi0 => plasticstate(p)%state (startindex :endindex ,:)
207  dlt%chi0 => plasticstate(p)%deltaState(startindex-o:endindex-o,:)
208 
209  startindex = endindex + 1
210  endindex = endindex + prm%sum_N_sl
211  stt%gamma0 => plasticstate(p)%state (startindex :endindex ,:)
212  dlt%gamma0 => plasticstate(p)%deltaState(startindex-o:endindex-o,:)
213 
214  plasticstate(p)%state0 = plasticstate(p)%state ! ToDo: this could be done centrally
215 
216  end associate
217 
218 !--------------------------------------------------------------------------------------------------
219 ! exit if any parameter is out of range
220  if (extmsg /= '') call io_error(211,ext_msg=trim(extmsg)//'('//plasticity_kinehardening_label//')')
221 
222  enddo
223 
224 
225 end subroutine plastic_kinehardening_init
226 
227 
228 !--------------------------------------------------------------------------------------------------
230 !--------------------------------------------------------------------------------------------------
231 pure module subroutine plastic_kinehardening_lpanditstangent(lp,dlp_dmp,mp,instance,of)
232 
233  real(pReal), dimension(3,3), intent(out) :: &
234  Lp
235  real(pReal), dimension(3,3,3,3), intent(out) :: &
236  dLp_dMp
237 
238  real(pReal), dimension(3,3), intent(in) :: &
239  Mp
240  integer, intent(in) :: &
241  instance, &
242  of
243 
244  integer :: &
245  i,k,l,m,n
246  real(pReal), dimension(param(instance)%sum_N_sl) :: &
247  gdot_pos,gdot_neg, &
248  dgdot_dtau_pos,dgdot_dtau_neg
249 
250  lp = 0.0_preal
251  dlp_dmp = 0.0_preal
252 
253  associate(prm => param(instance))
254 
255  call kinetics(mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
256 
257  do i = 1, prm%sum_N_sl
258  lp = lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
259  forall (k=1:3,l=1:3,m=1:3,n=1:3) &
260  dlp_dmp(k,l,m,n) = dlp_dmp(k,l,m,n) &
261  + dgdot_dtau_pos(i) * prm%P(k,l,i) * prm%nonSchmid_pos(m,n,i) &
262  + dgdot_dtau_neg(i) * prm%P(k,l,i) * prm%nonSchmid_neg(m,n,i)
263  enddo
264 
265  end associate
266 
267 end subroutine plastic_kinehardening_lpanditstangent
268 
269 
270 !--------------------------------------------------------------------------------------------------
272 !--------------------------------------------------------------------------------------------------
273 module subroutine plastic_kinehardening_dotstate(mp,instance,of)
274 
275  real(pReal), dimension(3,3), intent(in) :: &
276  Mp
277  integer, intent(in) :: &
278  instance, &
279  of
280 
281  real(pReal) :: &
282  sumGamma
283  real(pReal), dimension(param(instance)%sum_N_sl) :: &
284  gdot_pos,gdot_neg
285 
286 
287  associate(prm => param(instance), stt => state(instance), dot => dotstate(instance))
288 
289  call kinetics(mp,instance,of,gdot_pos,gdot_neg)
290  dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
291  sumgamma = sum(stt%accshear(:,of))
292 
293 
294  dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
295  * ( prm%theta1 &
296  + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumgamma/prm%tau1) &
297  * exp(-sumgamma*prm%theta0/prm%tau1) &
298  )
299 
300  dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
301  ( prm%theta1_b + &
302  (prm%theta0_b - prm%theta1_b &
303  + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
304  ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
305  )
306 
307  end associate
308 
309 end subroutine plastic_kinehardening_dotstate
310 
311 
312 !--------------------------------------------------------------------------------------------------
314 !--------------------------------------------------------------------------------------------------
315 module subroutine plastic_kinehardening_deltastate(mp,instance,of)
316 
317  real(pReal), dimension(3,3), intent(in) :: &
318  Mp
319  integer, intent(in) :: &
320  instance, &
321  of
322 
323  real(pReal), dimension(param(instance)%sum_N_sl) :: &
324  gdot_pos,gdot_neg, &
325  sense
326 
327  associate(prm => param(instance), stt => state(instance), dlt => deltastate(instance))
328 
329  call kinetics(mp,instance,of,gdot_pos,gdot_neg)
330  sense = merge(state(instance)%sense(:,of), & ! keep existing...
331  sign(1.0_preal,gdot_pos+gdot_neg), & ! ...or have a defined
332  deq0(gdot_pos+gdot_neg,1e-10_preal)) ! current sense of shear direction
333 
334 # 338 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_kinehardening.f90"
335 
336 !--------------------------------------------------------------------------------------------------
337 ! switch in sense of shear?
338  where(dneq(sense,stt%sense(:,of),0.1_preal))
339  dlt%sense (:,of) = sense - stt%sense(:,of) ! switch sense
340  dlt%chi0 (:,of) = abs(stt%crss_back(:,of)) - stt%chi0(:,of) ! remember current backstress magnitude
341  dlt%gamma0(:,of) = stt%accshear(:,of) - stt%gamma0(:,of) ! remember current accumulated shear
342  else where
343  dlt%sense (:,of) = 0.0_preal
344  dlt%chi0 (:,of) = 0.0_preal
345  dlt%gamma0(:,of) = 0.0_preal
346  end where
347 
348  end associate
349 
350 end subroutine plastic_kinehardening_deltastate
351 
352 
353 !--------------------------------------------------------------------------------------------------
355 !--------------------------------------------------------------------------------------------------
356 module subroutine plastic_kinehardening_results(instance,group)
357 
358  integer, intent(in) :: instance
359  character(len=*), intent(in) :: group
360 
361  integer :: o
362 
363  associate(prm => param(instance), stt => state(instance))
364  outputsloop: do o = 1,size(prm%output)
365  select case(trim(prm%output(o)))
366  case('resistance')
367  if(prm%sum_N_sl>0) call results_writedataset(group,stt%crss,'xi_sl', &
368  'resistance against plastic slip','Pa')
369  case('backstress') ! ToDo: should be 'tau_back'
370  if(prm%sum_N_sl>0) call results_writedataset(group,stt%crss_back,'tau_back', &
371  'back stress against plastic slip','Pa')
372  case ('sense')
373  if(prm%sum_N_sl>0) call results_writedataset(group,stt%sense,'sense_of_shear', &
374  'tbd','1')
375  case ('chi0')
376  if(prm%sum_N_sl>0) call results_writedataset(group,stt%chi0,'chi0', &
377  'tbd','Pa')
378  case ('gamma0')
379  if(prm%sum_N_sl>0) call results_writedataset(group,stt%gamma0,'gamma0', &
380  'tbd','1')
381  case ('accumulatedshear')
382  if(prm%sum_N_sl>0) call results_writedataset(group,stt%accshear,'gamma_sl', &
383  'plastic shear','1')
384  end select
385  enddo outputsloop
386  end associate
387 
388 end subroutine plastic_kinehardening_results
389 
390 
391 !--------------------------------------------------------------------------------------------------
393 ! stress.
395 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
396 ! have the optional arguments at the end.
397 !--------------------------------------------------------------------------------------------------
398 pure subroutine kinetics(Mp,instance,of, &
399  gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
400 
401  real(pReal), dimension(3,3), intent(in) :: &
402  Mp
403  integer, intent(in) :: &
404  instance, &
405  of
406 
407  real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
408  gdot_pos, &
409  gdot_neg
410  real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
411  dgdot_dtau_pos, &
412  dgdot_dtau_neg
413 
414  real(pReal), dimension(param(instance)%sum_N_sl) :: &
415  tau_pos, &
416  tau_neg
417  integer :: i
418 
419  associate(prm => param(instance), stt => state(instance))
420 
421  do i = 1, prm%sum_N_sl
422  tau_pos(i) = math_tensordot(mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
423  tau_neg(i) = merge(math_tensordot(mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
424  0.0_preal, prm%nonSchmidActive)
425  enddo
426 
427  where(dneq0(tau_pos))
428  gdot_pos = prm%gdot0 * merge(0.5_preal,1.0_preal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
429  * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
430  else where
431  gdot_pos = 0.0_preal
432  end where
433 
434  where(dneq0(tau_neg))
435  gdot_neg = prm%gdot0 * 0.5_preal & ! only used if non-Schmid active, always 1/2
436  * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
437  else where
438  gdot_neg = 0.0_preal
439  end where
440 
441  if (present(dgdot_dtau_pos)) then
442  where(dneq0(gdot_pos))
443  dgdot_dtau_pos = gdot_pos*prm%n/tau_pos
444  else where
445  dgdot_dtau_pos = 0.0_preal
446  end where
447  endif
448  if (present(dgdot_dtau_neg)) then
449  where(dneq0(gdot_neg))
450  dgdot_dtau_neg = gdot_neg*prm%n/tau_neg
451  else where
452  dgdot_dtau_neg = 0.0_preal
453  end where
454  endif
455  end associate
456 
457 end subroutine kinetics
458 
459 end submodule plastic_kinehardening
config
Reads in the material configuration from file.
Definition: config.f90:13
constitutive
elasticity, plasticity, internal microstructure state
Definition: constitutive.f90:10