FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
Elastoplastic.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
7  use hecmw_util
8  use mmaterial
10 
11  implicit none
12 
13  real(kind=kreal), private, parameter :: id(6,6) = reshape( &
14  & (/ 2.d0/3.d0, -1.d0/3.d0, -1.d0/3.d0, 0.d0, 0.d0, 0.d0, &
15  & -1.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0, 0.d0, 0.d0, 0.d0, &
16  & -1.d0/3.d0, -1.d0/3.d0, 2.d0/3.d0, 0.d0, 0.d0, 0.d0, &
17  & 0.d0, 0.d0, 0.d0, 0.5d0, 0.d0, 0.d0, &
18  & 0.d0, 0.d0, 0.d0, 0.d0, 0.5d0, 0.d0, &
19  & 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.5d0/), &
20  & (/6, 6/))
21 
22 contains
23 
25  subroutine calelastoplasticmatrix( matl, sectType, stress, istat, extval, plstrain, D, temperature )
26  type( tmaterial ), intent(in) :: matl
27  integer, intent(in) :: secttype
28  real(kind=kreal), intent(in) :: stress(6)
29  real(kind=kreal), intent(in) :: extval(:)
30  real(kind=kreal), INTENT(IN) :: plstrain
31  integer, intent(in) :: istat
32  real(kind=kreal), intent(out) :: d(:,:)
33  real(kind=kreal), optional :: temperature
34 
35  integer :: i,j,ytype
36  logical :: kinematic
37  real(kind=kreal) :: dum, dj1(6), dj2(6), dj3(6), a(6), de(6,6), g, dlambda
38  real(kind=kreal) :: c1,c2,c3, back(6)
39  real(kind=kreal) :: j1,j2,j3, fai, sita, harden, khard, da(6), devia(6)
40 
41  ytype = getyieldfunction( matl%mtype )
42  if( ytype==3 ) then
43  call uelastoplasticmatrix( matl, stress, istat, extval, d )
44  return
45  endif
46  if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
47  kinematic = iskinematicharden( matl%mtype )
48  khard = 0.d0
49  if( kinematic ) then
50  back(1:6) = extval(2:7)
51  khard = calkinematicharden( matl, extval(1) )
52  endif
53  if( present( temperature ) ) then
54  call calelasticmatrix( matl, secttype, de, temperature )
55  else
56  call calelasticmatrix( matl, secttype, de )
57  endif
58 
59  j1 = (stress(1)+stress(2)+stress(3))
60  devia(1:3) = stress(1:3)-j1/3.d0
61  devia(4:6) = stress(4:6)
62  if( kinematic ) devia = devia-back
63  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
64  dot_product( devia(4:6), devia(4:6) )
65 
66  d(:,:) = de(:,:)
67  if( istat == 0 ) return ! elastic state
68 
69  !derivative of J2
70  dj2(1:3) = devia(1:3)
71  dj2(4:6) = 2.d0*devia(4:6)
72  dj2 = dj2/( 2.d0*dsqrt(j2) )
73  if( present(temperature) ) then
74  harden = calhardencoeff( matl, extval(1), temperature )
75  else
76  harden = calhardencoeff( matl, extval(1) )
77  endif
78 
79  select case (ytype)
80  case (0) ! Mises or. Isotropic
81  a(1:6) = devia(1:6)/dsqrt(2.d0*j2)
82  g = de(4,4)
83  dlambda = extval(1)-plstrain
84  c3 = dsqrt(3.d0*j2)+3.d0*g*dlambda !trial mises stress
85  c1 = 6.d0*dlambda*g*g/c3
86  dum = 3.d0*g+khard+harden
87  c2 = 6.d0*g*g*(dlambda/c3-1.d0/dum)
88 
89  do i=1,6
90  do j=1,6
91  d(i,j) = de(i,j) - c1*id(i,j) + c2*a(i)*a(j)
92  enddo
93  enddo
94 
95  return
96  case (1) ! Mohr-Coulomb
97  fai = matl%variables(m_plconst3)
98  j3 = devia(1)*devia(2)*devia(3) &
99  +2.d0* devia(4)*devia(5)*devia(6) &
100  -devia(6)*devia(2)*devia(6) &
101  -devia(4)*devia(4)*devia(3) &
102  -devia(1)*devia(5)*devia(5)
103  sita = -3.d0*dsqrt(3.d0)*j3/( 2.d0*(j2**1.5d0) )
104  if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) then
105  c1 = 0.d0
106  c2 = dsqrt(3.d0)
107  c3 = 0.d0
108  else
109  if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
110  sita = asin( sita )/3.d0
111  c2 = cos(sita)*( 1.d0*tan(sita)*tan(3.d0*sita) + sin(fai)* &
112  ( tan(3.d0*sita)-tan(sita )/dsqrt(3.d0) ) )
113  c1 = sin(fai)/3.d0
114  c3 = dsqrt(3.d0)*sin(sita)+cos(sita)*sin(fai)/(2.d0*j2*cos(3.d0*sita))
115  endif
116  ! deirivative of j1
117  dj1(1:3) = 1.d0
118  dj1(4:6) = 0.d0
119  ! deirivative of j3
120  dj3(1) = devia(2)*devia(3)-devia(5)*devia(5)+j2/3.d0
121  dj3(2) = devia(1)*devia(3)-devia(6)*devia(6)+j2/3.d0
122  dj3(3) = devia(1)*devia(2)-devia(4)*devia(4)+j2/3.d0
123  dj3(4) = 2.d0*(devia(5)*devia(6)-devia(3)*devia(4))
124  dj3(5) = 2.d0*(devia(4)*devia(6)-devia(1)*devia(5))
125  dj3(6) = 2.d0*(devia(4)*devia(5)-devia(2)*devia(6))
126  a(:) = c1*dj1 + c2*dj2 + c3*dj3
127  case (2) ! Drucker-Prager
128  fai = matl%variables(m_plconst3)
129  ! deirivative of j1
130  dj1(1:3) = 1.d0
131  dj1(4:6) = 0.d0
132  a(:) = fai*dj1(:) + dj2(:)
133  end select
134 
135  da = matmul( de, a )
136  dum = harden + khard+ dot_product( da, a )
137  do i=1,6
138  do j=1,6
139  d(i,j) = de(i,j) - da(i)*da(j)/dum
140  enddo
141  enddo
142 
143  end subroutine
144 
146  real(kind=kreal) function cal_equivalent_stress(matl, stress, extval)
147  type( tmaterial ), intent(in) :: matl
148  real(kind=kreal), intent(in) :: stress(6)
149  real(kind=kreal), intent(in) :: extval(:)
150 
151  integer :: ytype
152  logical :: kinematic
153  real(kind=kreal) :: eqvs, sita, fai, j1,j2,j3, devia(6)
154  real(kind=kreal) :: back(6)
155  kinematic = iskinematicharden( matl%mtype )
156  if( kinematic ) back(1:6) = extval(2:7)
157 
158  ytype = getyieldfunction( matl%mtype )
159  j1 = (stress(1)+stress(2)+stress(3))
160  devia(1:3) = stress(1:3)-j1/3.d0
161  devia(4:6) = stress(4:6)
162  if( kinematic ) devia = devia-back
163  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
164  dot_product( devia(4:6), devia(4:6) )
165 
166  select case (ytype)
167  case (0) ! Mises or. Isotropic
168  eqvs = dsqrt( 3.d0*j2 )
169  case (1) ! Mohr-Coulomb
170  fai = matl%variables(m_plconst1)
171  j3 = devia(1)*devia(2)*devia(3) &
172  +2.d0* devia(4)*devia(5)*devia(6) &
173  -devia(6)*devia(2)*devia(6) &
174  -devia(4)*devia(4)*devia(3) &
175  -devia(1)*devia(5)*devia(5)
176  sita = -3.d0*dsqrt(3.d0)*j3/( 2.d0*(j2**1.5d0) )
177  if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) sita=sign(1.d0, sita)
178  if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
179  sita = asin( sita )/3.d0
180  eqvs = (cos(sita)-sin(sita)*sin(fai)/dsqrt(3.d0))*dsqrt(j2) &
181  +j1*sin(fai)/3.d0
182  case (2) ! Drucker-Prager
183  eqvs = dsqrt(j2)
184  case default
185  eqvs = -1.d0
186  end select
187 
188  cal_equivalent_stress = eqvs
189  end function
190 
192  real(kind=kreal) function cal_mises_strain( strain )
193  real(kind=kreal), intent(in) :: strain(6)
194  cal_mises_strain = 2.d0*dot_product( strain(1:3), strain(1:3) )
195  cal_mises_strain = cal_mises_strain+ dot_product( strain(4:6), strain(4:6) )
196  cal_mises_strain = dsqrt( cal_mises_strain/3.d0 )
197  end function
198 
200  real(kind=kreal) function calhardencoeff( matl, pstrain, temp )
201  type( tmaterial ), intent(in) :: matl
202  real(kind=kreal), intent(in) :: pstrain
203  real(kind=kreal), intent(in), optional :: temp
204 
205  integer :: htype
206  logical :: ierr
207  real(kind=kreal) :: s0, s1,s2, ef, ina(2)
208 
209  calhardencoeff = -1.d0
210  htype = gethardentype( matl%mtype )
211  select case (htype)
212  case (0) ! Linear hardening
213  calhardencoeff = matl%variables(m_plconst2)
214  case (1) ! Multilinear approximation
215  if( present(temp) ) then
216  ina(1) = temp; ina(2)=pstrain
217  call fetch_tablegrad( mc_yield, ina, matl%dict, calhardencoeff, ierr )
218  ! print *, ina, calHardenCoeff; pause
219  else
220  ina(1)=pstrain
221  call fetch_tablegrad( mc_yield, ina(1:1), matl%dict, calhardencoeff, ierr )
222  endif
223  ! print *,1, calHardenCoeff; pause
224  case (2) ! Swift
225  s0= matl%variables(m_plconst1)
226  s1= matl%variables(m_plconst2)
227  s2= matl%variables(m_plconst3)
228  calhardencoeff = s1*s2*( s0+pstrain )**(s2-1)
229  case (3) ! Ramberg-Osgood
230  s0= matl%variables(m_plconst1)
231  s1= matl%variables(m_plconst2)
232  s2= matl%variables(m_plconst3)
233  if( present(temp) ) then
234  ef = calcurryield( matl, pstrain, temp )
235  else
236  ef = calcurryield( matl, pstrain )
237  endif
238  calhardencoeff = s1*(ef/s1)**(1.d0-s2) /(s0*s2)
239  case(4) ! Prager
240  calhardencoeff = 0.d0
241  case(5) ! Prager+linear
242  calhardencoeff = matl%variables(m_plconst2)
243  end select
244  end function
245 
247  real(kind=kreal) function calkinematicharden( matl, pstrain )
248  type( tmaterial ), intent(in) :: matl
249  real(kind=kreal), intent(in) :: pstrain
250 
251  integer :: htype
252  htype = gethardentype( matl%mtype )
253  select case (htype)
254  case(4, 5) ! Prager
255  calkinematicharden = matl%variables(m_plconst3)
256  case default
257  calkinematicharden = 0.d0
258  end select
259  end function
260 
262  real(kind=kreal) function calcurrkinematic( matl, pstrain )
263  type( tmaterial ), intent(in) :: matl
264  real(kind=kreal), intent(in) :: pstrain
265 
266  integer :: htype
267  htype = gethardentype( matl%mtype )
268  select case (htype)
269  case(4, 5) ! Prager
270  calcurrkinematic = matl%variables(m_plconst3)*pstrain
271  case default
272  calcurrkinematic = 0.d0
273  end select
274  end function
275 
277  real(kind=kreal) function calcurryield( matl, pstrain, temp )
278  type( tmaterial ), intent(in) :: matl
279  real(kind=kreal), intent(in) :: pstrain
280  real(kind=kreal), intent(in), optional :: temp
281 
282  integer :: htype
283  real(kind=kreal) :: s0, s1,s2, ina(2), outa(1)
284  logical :: ierr
285  calcurryield = -1.d0
286  htype = gethardentype( matl%mtype )
287 
288  select case (htype)
289  case (0, 5) ! Linear hardening, Linear+Parger hardening
290  calcurryield = matl%variables(m_plconst1)+matl%variables(m_plconst2)*pstrain
291  case (1) ! Multilinear approximation
292  if( present(temp) ) then
293  ina(1) = temp; ina(2)=pstrain
294  call fetch_tabledata(mc_yield, matl%dict, outa, ierr, ina)
295  else
296  ina(1) = pstrain
297  call fetch_tabledata(mc_yield, matl%dict, outa, ierr, ina(1:1))
298  endif
299  if( ierr ) stop "Fail to get yield stress!"
300  calcurryield = outa(1)
301  case (2) ! Swift
302  s0= matl%variables(m_plconst1)
303  s1= matl%variables(m_plconst2)
304  s2= matl%variables(m_plconst3)
305  calcurryield = s1*( s0+pstrain )**s2
306  case (3) ! Ramberg-Osgood
307  s0= matl%variables(m_plconst1)
308  s1= matl%variables(m_plconst2)
309  s2= matl%variables(m_plconst3)
310  if( pstrain<=s0 ) then
311  calcurryield = s1
312  else
313  calcurryield = s1*( pstrain/s0 )**(1.d0/s2)
314  endif
315  case (4) ! Parger hardening
316  calcurryield = matl%variables(m_plconst1)
317  end select
318  end function
319 
321  real(kind=kreal) function calyieldfunc( matl, stress, extval, temp )
322  type( tmaterial ), intent(in) :: matl
323  real(kind=kreal), intent(in) :: stress(6)
324  real(kind=kreal), intent(in) :: extval(:)
325  real(kind=kreal), intent(in), optional :: temp
326 
327  integer :: ytype
328  logical :: kinematic
329  real(kind=kreal) :: eqvs, sita, eta, fai, j1,j2,j3, f, devia(6)
330  real(kind=kreal) :: pstrain, back(6)
331 
332  f = 0.0d0
333 
334  kinematic = iskinematicharden( matl%mtype )
335  if( kinematic ) back(1:6) = extval(2:7)
336 
337  pstrain = extval(1)
338  ytype = getyieldfunction( matl%mtype )
339  j1 = (stress(1)+stress(2)+stress(3))
340  devia(1:3) = stress(1:3)-j1/3.d0
341  devia(4:6) = stress(4:6)
342  if( kinematic ) devia = devia-back
343 
344  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
345  dot_product( devia(4:6), devia(4:6) )
346  if( present(temp) ) then
347  eqvs = calcurryield( matl, pstrain, temp )
348  else
349  eqvs = calcurryield( matl, pstrain )
350  endif
351 
352  select case (ytype)
353  case (0) ! Mises or. Isotropic
354  f = dsqrt( 3.d0*j2 ) - eqvs
355  case (1) ! Mohr-Coulomb
356  fai = matl%variables(m_plconst3)
357  j3 = devia(1)*devia(2)*devia(3) &
358  +2.d0* devia(4)*devia(5)*devia(6) &
359  -devia(2)*devia(6)*devia(6) &
360  -devia(3)*devia(4)*devia(4) &
361  -devia(1)*devia(5)*devia(5)
362  sita = -3.d0*dsqrt(3.d0)*j3/( 2.d0*(j2**1.5d0) )
363  if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) sita=sign(1.d0, sita)
364  if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
365  sita = asin( sita )/3.d0
366  f = (cos(sita)-sin(sita)*sin(fai)/dsqrt(3.d0))*dsqrt(j2) &
367  +j1*sin(fai)/3.d0 - eqvs*cos(fai)
368  case (2) ! Drucker-Prager
369  eta = matl%variables(m_plconst3)
370  f = dsqrt(j2) + eta*j1 - eqvs*matl%variables(m_plconst4)
371  end select
372 
373  calyieldfunc = f
374  end function
375 
377  subroutine backwardeuler( matl, stress, plstrain, istat, fstat, temp )
378  use m_utilities, only : eigen3
379  type( tmaterial ), intent(in) :: matl
380  real(kind=kreal), intent(inout) :: stress(6)
381  real(kind=kreal), intent(in) :: plstrain
382  integer, intent(inout) :: istat
383  real(kind=kreal), intent(inout) :: fstat(:)
384  real(kind=kreal), intent(in), optional :: temp
385 
386  real(kind=kreal), parameter :: tol =1.d-3
387  integer, parameter :: maxiter = 5
388  real(kind=kreal) :: dlambda, f, mat(3,3)
389  integer :: i,ytype, maxp(1), minp(1), mm
390  real(kind=kreal) :: youngs, poisson, pstrain, dum, ina(1), ee(2)
391  real(kind=kreal) :: j1,j2,j3, h, kh, kk, dd, yd, g, k, devia(6)
392  real(kind=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3)
393  real(kind=kreal) :: sita, fai, trialprn(3)
394  real(kind=kreal) :: fstat_bak(7)
395  logical :: kinematic, ierr
396  real(kind=kreal) :: betan, back(6)
397 
398  f = 0.0d0
399 
400  ytype = getyieldfunction( matl%mtype )
401  if( ytype==3 ) then
402  call ubackwardeuler( matl, stress, istat, fstat )
403  return
404  endif
405 
406  pstrain = plstrain
407  if(iskinematicharden( matl%mtype ))fstat_bak(2:7)= fstat(8:13)
408  fstat_bak(1) = plstrain
409  if( present(temp) ) then
410  f = calyieldfunc( matl, stress, fstat_bak, temp )
411  else
412  f = calyieldfunc( matl, stress, fstat_bak )
413  endif
414  if( dabs(f)<tol ) then ! yielded
415  istat = 1
416  return
417  elseif( f<0.d0 ) then ! not yielded or unloading
418  istat =0
419  return
420  endif
421  istat = 1 ! yielded
422  kh = 0.d0; kk=0.d0; betan=0.d0; back(:)=0.d0
423 
424  kinematic = iskinematicharden( matl%mtype )
425  if( kinematic ) then
426  back(1:6) = fstat(8:13)
427  betan = calcurrkinematic( matl, pstrain )
428  endif
429 
430  j1 = (stress(1)+stress(2)+stress(3))/3.d0
431  devia(1:3) = stress(1:3)-j1
432  devia(4:6) = stress(4:6)
433  if( kinematic ) devia = devia-back
434  yd = cal_equivalent_stress(matl, stress, fstat)
435 
436  if( present(temp) ) then
437  ina(1) = temp
438  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
439  else
440  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr )
441  endif
442  if( ierr ) then
443  stop " fail to fetch young's modulus in elastoplastic calculation"
444  else
445  youngs = ee(1)
446  poisson = ee(2)
447  endif
448  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
449  g = youngs/ ( 2.d0*(1.d0+poisson) )
450  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
451  dlambda = 0.d0
452 
453  if( ytype==0 ) then ! Mises or. Isotropic
454  do i=1,maxiter
455  if( present(temp) ) then
456  h= calhardencoeff( matl, pstrain+dlambda, temp )
457  else
458  h= calhardencoeff( matl, pstrain+dlambda )
459  endif
460  if( kinematic ) then
461  kh = calkinematicharden( matl, pstrain+dlambda )
462  endif
463  dd= 3.d0*g+h+kh
464  dlambda = dlambda+f/dd
465  if( dlambda<0.d0 ) then
466  dlambda = 0.d0
467  istat=0; exit
468  endif
469  if( present(temp) ) then
470  dum = calcurryield( matl, pstrain+dlambda, temp )
471  else
472  dum = calcurryield( matl, pstrain+dlambda )
473  endif
474  if( kinematic ) then
475  kk = calcurrkinematic( matl, pstrain+dlambda )
476  endif
477  f = yd-3.d0*g*dlambda-dum -(kk-betan)
478  if( dabs(f)<tol*tol ) exit
479  enddo
480  pstrain = pstrain+dlambda
481  if( kinematic ) then
482  kk = calcurrkinematic( matl, pstrain )
483  fstat(2:7) = back(:)+(kk-betan)*devia(:)/yd
484  endif
485  devia(:) = (1.d0-3.d0*dlambda*g/yd)*devia(:)
486  stress(1:3) = devia(1:3)+j1
487  stress(4:6) = devia(4:6)
488  stress(:)= stress(:)+back(:)
489  elseif(ytype==1) then ! Mohr-Coulomb
490  fai = matl%variables(m_plconst3)
491 
492  ! do j=1,MAXITER
493  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
494  dot_product( devia(4:6), devia(4:6) )
495  j3 = devia(1)*devia(2)*devia(3) &
496  +2.d0* devia(4)*devia(5)*devia(6) &
497  -devia(6)*devia(2)*devia(6) &
498  -devia(4)*devia(4)*devia(3) &
499  -devia(1)*devia(5)*devia(5)
500  sita = -3.d0*dsqrt(3.d0)*j3/( 2.d0*(j2**1.5d0) )
501  if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) sita=sign(1.d0, sita)
502  if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
503  sita = asin( sita )/3.d0
504  do mm=1,6
505  if( dabs(stress(mm))<1.d-10 ) stress(mm)=0.d0
506  enddo
507  call eigen3( stress, prnstre, prnprj )
508  trialprn = prnstre
509  maxp = maxloc( prnstre )
510  minp = minloc( prnstre )
511  mm = 1
512  if( maxp(1)==1 .or. minp(1)==1 ) mm =2
513  if( maxp(1)==2 .or. minp(1)==2 ) mm =3
514  do i=1,maxiter
515  if( present(temp) ) then
516  h= calhardencoeff( matl, pstrain, temp )
517  else
518  h= calhardencoeff( matl, pstrain )
519  endif
520  dd= 4.d0*g*( 1.d0+sin(fai)*sin(sita)/3.d0 )+4.d0*k &
521  *sin(fai)*sin(sita)+4.d0*h*cos(fai)*cos(fai)
522  dlambda = dlambda+f/dd
523  if( 2.d0*dlambda*cos(fai)<0.d0 ) then
524  if( cos(fai)==0.d0 ) stop "Math error in return mapping"
525  dlambda = 0.d0
526  istat=0; exit
527  endif
528  dum = pstrain + 2.d0*dlambda*cos(fai)
529  if( present(temp) ) then
530  yd = calcurryield( matl, dum, temp )
531  else
532  yd = calcurryield( matl, dum )
533  endif
534  f = prnstre(maxp(1))-prnstre(minp(1))+ &
535  (prnstre(maxp(1))+prnstre(minp(1)))*sin(fai)- &
536  (4.d0*g*(1.d0+sin(fai)*sin(sita)/3.d0)+4.d0*k*sin(fai) &
537  *sin(sita))*dlambda-2.d0*yd*cos(fai)
538  if( dabs(f)<tol ) exit
539  enddo
540  pstrain = pstrain + 2.d0*dlambda*cos(fai)
541  prnstre(maxp(1)) = prnstre(maxp(1))-(2.d0*g*(1.d0+sin(fai)/3.d0) &
542  + 2.d0*k*sin(fai) )*dlambda
543  prnstre(minp(1)) = prnstre(minp(1))+(2.d0*g*(1.d0-sin(fai)/3.d0) &
544  - 2.d0*k*sin(fai) )*dlambda
545  prnstre(mm) = prnstre(mm)+(4.d0*g/3.d0-2.d0*k)*sin(fai)*dlambda
546 
547  tstre(:,:) = 0.d0
548  tstre(1,1)= prnstre(1); tstre(2,2)=prnstre(2); tstre(3,3)=prnstre(3)
549  mat= matmul( prnprj, tstre )
550  mat= matmul( mat, transpose(prnprj) )
551  stress(1) = mat(1,1)
552  stress(2) = mat(2,2)
553  stress(3) = mat(3,3)
554  stress(4) = mat(1,2)
555  stress(5) = mat(2,3)
556  stress(6) = mat(3,1)
557  elseif(ytype==2) then ! Drucker-Prager
558  fai = matl%variables(m_plconst3)
559  dum = matl%variables(m_plconst4)
560  do i=1,maxiter
561  if( present(temp) ) then
562  h= calhardencoeff( matl, pstrain, temp )
563  else
564  h= calhardencoeff( matl, pstrain )
565  endif
566  dd= g+k*fai*fai+h*dum*dum
567  dlambda = dlambda+f/dd
568  if( dum*dlambda<0.d0 ) then
569  if( dum==0.d0 ) stop "Math error in return mapping"
570  dlambda = 0.d0
571  istat=0; exit
572  endif
573  if( present(temp) ) then
574  f = calcurryield( matl, pstrain+dum*dlambda, temp )
575  else
576  f = calcurryield( matl, pstrain+dum*dlambda )
577  endif
578  f = yd-g*dlambda+fai*(j1-k*fai*dlambda)- dum*f
579  if( dabs(f)<tol*tol ) exit
580  enddo
581  pstrain = pstrain+dum*dlambda
582  devia(:) = (1.d0-g*dlambda/yd)*devia(:)
583  j1 = j1-k*fai*dlambda
584  stress(1:3) = devia(1:3)+j1
585  stress(4:6) = devia(4:6)
586  end if
587 
588  fstat(1) = pstrain
589  end subroutine backwardeuler
590 
592  subroutine updateepstate( gauss )
593  use mmechgauss
594  type(tgaussstatus), intent(inout) :: gauss ! status of curr gauss point
595  gauss%plstrain= gauss%fstatus(1)
596  if(iskinematicharden(gauss%pMaterial%mtype)) then
597  gauss%fstatus(8:13) =gauss%fstatus(2:7)
598  endif
599  end subroutine
600 
601 end module m_elastoplastic
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
This module provides functions for elastic material.
subroutine calelasticmatrix(matl, sectType, D, temp)
Calculate isotropic elastic matrix.
This module provide functions for elastoplastic calculation.
real(kind=kreal) function calcurryield(matl, pstrain, temp)
This function calcualtes current yield stress.
subroutine updateepstate(gauss)
Clear elatoplastic state.
subroutine backwardeuler(matl, stress, plstrain, istat, fstat, temp)
This subroutine does backward-Euler return calculation.
subroutine calelastoplasticmatrix(matl, sectType, stress, istat, extval, plstrain, D, temperature)
This subroutine calculates elastoplastic constitutive relation.
real(kind=kreal) function cal_equivalent_stress(matl, stress, extval)
This subrouitne calculate equivalent stress.
real(kind=kreal) function calyieldfunc(matl, stress, extval, temp)
This function calcualtes yield state.
real(kind=kreal) function cal_mises_strain(strain)
This subrouitne calculate equivalent stress.
real(kind=kreal) function calkinematicharden(matl, pstrain)
This function calcualtes kinematic hardening coefficient.
real(kind=kreal) function calhardencoeff(matl, pstrain, temp)
This function calcualtes hardening coefficient.
real(kind=kreal) function calcurrkinematic(matl, pstrain)
This function calcualtes state of kinematic hardening.
This module provides aux functions.
Definition: utilities.f90:6
subroutine eigen3(tensor, eigval, princ)
Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using Jacobi iteration adapted from numeri...
Definition: utilities.f90:106
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer function getyieldfunction(mtype)
Get type of yield function.
Definition: material.f90:294
integer(kind=kint), parameter m_plconst4
Definition: material.f90:93
integer(kind=kint), parameter m_plconst1
Definition: material.f90:90
integer function gethardentype(mtype)
Get type of hardening.
Definition: material.f90:306
integer(kind=kint), parameter d3
Definition: material.f90:76
integer(kind=kint), parameter m_plconst2
Definition: material.f90:91
character(len=dict_key_length) mc_yield
Definition: material.f90:121
logical function iskinematicharden(mtype)
If it is a kinematic hardening material?
Definition: material.f90:318
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:119
integer(kind=kint), parameter m_plconst3
Definition: material.f90:92
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
Stucture to management all material relates data.
Definition: material.f90:144
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
subroutine uelastoplasticmatrix(matl, stress, istat, fstat, D)
This subroutine read in used-defined material properties tangent This subroutine calculates elastopla...
Definition: uyield.f90:9
subroutine ubackwardeuler(matl, stress, istat, fstat)
This subroutine does backward-Euler return calculation.
Definition: uyield.f90:20