FrontISTR  5.7.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 !-------------------------------------------------------------------------------
6 module m_elastoplastic
7  use hecmw_util
8  use mmaterial
9  use m_elasticlinear
10  use muyield
11 
12  implicit none
13 
14  private
15  public :: calelastoplasticmatrix
16  public :: backwardeuler
17  public :: updateepstate
18 
19  real(kind=kreal), parameter :: id(6,6) = reshape( &
20  & (/ 2.d0/3.d0, -1.d0/3.d0, -1.d0/3.d0, 0.d0, 0.d0, 0.d0, &
21  & -1.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0, 0.d0, 0.d0, 0.d0, &
22  & -1.d0/3.d0, -1.d0/3.d0, 2.d0/3.d0, 0.d0, 0.d0, 0.d0, &
23  & 0.d0, 0.d0, 0.d0, 0.5d0, 0.d0, 0.d0, &
24  & 0.d0, 0.d0, 0.d0, 0.d0, 0.5d0, 0.d0, &
25  & 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.5d0/), &
26  & (/6, 6/))
27  real(kind=kreal), parameter :: i2(6) = (/ 1.d0, 1.d0, 1.d0, 0.d0, 0.d0, 0.d0 /)
28 
29  integer, parameter :: VM_ELASTIC = 0
30  integer, parameter :: VM_PLASTIC = 1
31 
32  integer, parameter :: MC_ELASTIC = 0
33  integer, parameter :: MC_PLASTIC_SURF = 1
34  integer, parameter :: MC_PLASTIC_RIGHT = 2
35  integer, parameter :: MC_PLASTIC_LEFT = 3
36  integer, parameter :: MC_PLASTIC_APEX = 4
37 
38  integer, parameter :: DP_ELASTIC = 0
39  integer, parameter :: DP_PLASTIC_SURF = 1
40  integer, parameter :: DP_PLASTIC_APEX = 2
41 
42 contains
43 
45  subroutine calelastoplasticmatrix( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
46  type( tmaterial ), intent(in) :: matl
47  integer, intent(in) :: secttype
48  real(kind=kreal), intent(in) :: stress(6)
49  real(kind=kreal), intent(in) :: extval(:)
50  real(kind=kreal), intent(in) :: plstrain
51  integer, intent(in) :: istat
52  real(kind=kreal), intent(out) :: d(:,:)
53  real(kind=kreal), intent(in) :: temperature
54  integer(kind=kint), intent(in), optional :: hdflag
55 
56  integer :: ytype,hdflag_in
57 
58  hdflag_in = 0
59  if( present(hdflag) ) hdflag_in = hdflag
60 
61  ytype = getyieldfunction( matl%mtype )
62  select case (ytype)
63  case (0)
64  call calelastoplasticmatrix_vm( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
65  case (1)
66  call calelastoplasticmatrix_mc( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
67  case (2)
68  call calelastoplasticmatrix_dp( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
69  case (3)
70  call uelastoplasticmatrix( matl%variables, stress, istat, extval, plstrain, d, temperature, hdflag_in )
71  end select
72  end subroutine calelastoplasticmatrix
73 
75  subroutine calelastoplasticmatrix_vm( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
76  type( tmaterial ), intent(in) :: matl
77  integer, intent(in) :: secttype
78  real(kind=kreal), intent(in) :: stress(6)
79  real(kind=kreal), intent(in) :: extval(:)
80  real(kind=kreal), intent(in) :: plstrain
81  integer, intent(in) :: istat
82  real(kind=kreal), intent(out) :: d(:,:)
83  real(kind=kreal), intent(in) :: temperature
84  integer(kind=kint), intent(in) :: hdflag
85 
86  integer :: i,j
87  logical :: kinematic
88  real(kind=kreal) :: dum, a(6), g, dlambda
89  real(kind=kreal) :: c1,c2,c3, back(6)
90  real(kind=kreal) :: j1,j2, harden, khard, devia(6)
91 
92  if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
93 
94  call calelasticmatrix( matl, secttype, d, temperature, hdflag=hdflag )
95  if( istat == vm_elastic ) return
96  if( hdflag == 2 ) return
97 
98  harden = calhardencoeff( matl, extval(1), temperature )
99 
100  kinematic = iskinematicharden( matl%mtype )
101  khard = 0.d0
102  if( kinematic ) then
103  back(1:6) = extval(2:7)
104  khard = calkinematicharden( matl, extval(1) )
105  endif
106 
107  j1 = (stress(1)+stress(2)+stress(3))
108  devia(1:3) = stress(1:3)-j1/3.d0
109  devia(4:6) = stress(4:6)
110  if( kinematic ) devia = devia-back
111  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
112  dot_product( devia(4:6), devia(4:6) )
113 
114  a(1:6) = devia(1:6)/sqrt(2.d0*j2)
115  g = d(4,4)
116  dlambda = extval(1)-plstrain
117  c3 = sqrt(3.d0*j2)+3.d0*g*dlambda !trial mises stress
118  c1 = 6.d0*dlambda*g*g/c3
119  dum = 3.d0*g+khard+harden
120  c2 = 6.d0*g*g*(dlambda/c3-1.d0/dum)
121 
122  do i=1,6
123  do j=1,6
124  d(i,j) = d(i,j) - c1*id(i,j) + c2*a(i)*a(j)
125  enddo
126  enddo
127  end subroutine calelastoplasticmatrix_vm
128 
130  subroutine calelastoplasticmatrix_mc( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
131  use m_utilities, only : eigen3,deriv_general_iso_tensor_func_3d
132  type( tmaterial ), intent(in) :: matl
133  integer, intent(in) :: secttype
134  real(kind=kreal), intent(in) :: stress(6)
135  real(kind=kreal), intent(in) :: extval(:)
136  real(kind=kreal), intent(in) :: plstrain
137  integer, intent(in) :: istat
138  real(kind=kreal), intent(out) :: d(:,:)
139  real(kind=kreal), intent(in) :: temperature
140  integer(kind=kint), intent(in) :: hdflag
141 
142  real(kind=kreal) :: g, k, harden, r2g, r2gd3, r4gd3, r2k, youngs, poisson
143  real(kind=kreal) :: phi, psi, cosphi, sinphi, cotphi, sinpsi, sphsps, r2cosphi, r4cos2phi
144  real(kind=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3), prnstra(3)
145  integer(kind=kint) :: m1, m2, m3
146  real(kind=kreal) :: c1,c2,c3, ca1, ca2, ca3, cam, cap, cd1, cd2, cd3, cdiag, coffd
147  real(kind=kreal) :: ck1, ck2, ck3
148  real(kind=kreal) :: dum, da, db, dc, dd, detinv
149  real(kind=kreal) :: dpsdpe(3,3)
150 
151  if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
152 
153  call calelasticmatrix( matl, secttype, d, temperature, hdflag=hdflag )
154  if( istat == mc_elastic ) return
155  if( hdflag == 2 ) return
156 
157  harden = calhardencoeff( matl, extval(1), temperature )
158  g = d(4,4)
159  k = d(1,1)-(4.d0/3.d0)*g
160  r2g = 2.d0*g
161  r2k = 2.d0*k
162  r2gd3 = r2g/3.d0
163  r4gd3 = 2.d0*r2gd3
164  youngs = 9.d0*k*g/(3.d0*k+g)
165  poisson = (3.d0*k-r2g)/(6.d0*k+r2g)
166 
167  phi = matl%variables(m_plconst3)
168  psi = matl%variables(m_plconst4)
169  sinphi = sin(phi)
170  cosphi = cos(phi)
171  sinpsi = sin(psi)
172  sphsps = sinphi*sinpsi
173  r2cosphi = 2.d0*cosphi
174  r4cos2phi = r2cosphi*r2cosphi
175 
176  call eigen3( stress, prnstre, prnprj )
177  m1 = maxloc( prnstre, 1 )
178  m3 = minloc( prnstre, 1 )
179  if( m1 == m3 ) then
180  m1 = 1; m2 = 2; m3 = 3
181  else
182  m2 = 6 - (m1 + m3)
183  endif
184 
185  c1 = 4.d0*(g*(1.d0+sphsps/3.d0)+k*sphsps)
186  if( istat==mc_plastic_surf ) then
187  dd= c1 + r4cos2phi*harden
188  cd1 = (r2g*(1.d0+sinpsi/3.d0) + r2k*sinpsi)/dd
189  cd2 = (r4gd3-r2k)*sinpsi/dd
190  cd3 = (r2g*(1.d0-sinpsi/3.d0) - r2k*sinpsi)/dd
191  cap = 1.d0+sinphi/3.d0
192  cam = 1.d0-sinphi/3.d0
193  ck1 = 1.d0-2.d0*cd1*sinphi
194  ck2 = 1.d0+2.d0*cd2*sinphi
195  ck3 = 1.d0+2.d0*cd3*sinphi
196  dpsdpe(m1,m1) = r2g*( 2.d0/3.d0-cd1*cap)+k*ck1
197  dpsdpe(m1,m2) = (k-r2gd3)*ck1
198  dpsdpe(m1,m3) = r2g*(-1.d0/3.d0+cd1*cam)+k*ck1
199  dpsdpe(m2,m1) = r2g*(-1.d0/3.d0+cd2*cap)+k*ck2
200  dpsdpe(m2,m2) = r4gd3*( 1.d0-cd2*sinphi)+k*ck2
201  dpsdpe(m2,m3) = r2g*(-1.d0/3.d0-cd2*cam)+k*ck2
202  dpsdpe(m3,m1) = r2g*(-1.d0/3.d0+cd3*cap)+k*ck3
203  dpsdpe(m3,m2) = (k-r2gd3)*ck3
204  dpsdpe(m3,m3) = r2g*( 2.d0/3.d0-cd3*cam)+k*ck3
205  else if( istat==mc_plastic_apex ) then
206  cotphi = cosphi/sinphi
207  dpsdpe(:,:) = k*(1.d0-(k/(k+harden*cotphi*cosphi/sinpsi)))
208  else ! EDGE
209  if( istat==mc_plastic_right ) then
210  c2 = r2g*(1.d0+sinphi+sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
211  else if( istat==mc_plastic_left ) then
212  c2 = r2g*(1.d0-sinphi-sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
213  endif
214  dum = r4cos2phi*harden
215  da = c1 + dum
216  db = c2 + dum
217  dc = db
218  dd = da
219  detinv = 1.d0/(da*dd-db*dc)
220  ca1 = r2g*(1.d0+sinphi/3.d0)+r2k*sinpsi
221  ca2 = (r4gd3-r2k)*sinpsi
222  ca3 = r2g*(1.d0-sinpsi/3.d0)-r2k*sinpsi
223  cdiag = k+r4gd3
224  coffd = k-r2gd3
225  if( istat==mc_plastic_right ) then
226  dpsdpe(m1,m1) = cdiag+ca1*(db-dd-da+dc)*(r2g+(r2k+r2gd3)*sinphi)*detinv
227  dpsdpe(m1,m2) = coffd+ca1*(r2g*(da-db)+((db-dd-da+dc)*(r2k+r2gd3)+(dd-dc)*r2g)*sinphi)*detinv
228  dpsdpe(m1,m3) = coffd+ca1*(r2g*(dd-dc)+((db-dd-da+dc)*(r2k+r2gd3)+(da-db)*r2g)*sinphi)*detinv
229  dpsdpe(m2,m1) = coffd+(ca2*(dd-db)+ca3*(da-dc))*(r2g+(r2k+r2gd3)*sinphi)*detinv
230  dpsdpe(m2,m2) = cdiag+(ca2*((r2k*(dd-db)-(db*r2gd3+dd*r4gd3))*sinphi+db*r2g) &
231  & +ca3*((r2k*(da-dc)+(da*r2gd3+dc*r4gd3))*sinphi-da*r2g))*detinv
232  dpsdpe(m2,m3) = coffd+(ca2*((r2k*(dd-db)+(db*r4gd3+dd*r2gd3))*sinphi-dd*r2g) &
233  & +ca3*((r2k*(da-dc)-(da*r4gd3+dc*r2gd3))*sinphi+dc*r2g))*detinv
234  dpsdpe(m3,m1) = coffd+(ca2*(da-dc)+ca3*(dd-db))*(r2g+(r2k+r2gd3)*sinphi)*detinv
235  dpsdpe(m3,m2) = coffd+(ca2*((r2k*(da-dc)+(da*r2gd3+dc*r4gd3))*sinphi-da*r2g) &
236  & +ca3*((r2k*(dd-db)-(db*r2gd3+dd*r4gd3))*sinphi+db*r2g))*detinv
237  dpsdpe(m3,m3) = cdiag+(ca2*((r2k*(da-dc)-(da*r4gd3+dc*r2gd3))*sinphi+dc*r2g) &
238  & +ca3*((r2k*(dd-db)+(db*r4gd3+dd*r2gd3))*sinphi-dd*r2g))*detinv
239  else if( istat==mc_plastic_left ) then
240  dpsdpe(m1,m1) = cdiag+(ca1*((r2k*(db-dd)-(db*r4gd3+dd*r2gd3))*sinphi-dd*r2g) &
241  & +ca2*((r2k*(da-dc)-(da*r4gd3+dc*r2gd3))*sinphi-dc*r2g))*detinv
242  dpsdpe(m1,m2) = coffd+(ca1*((r2k*(db-dd)+(db*r2gd3+dd*r4gd3))*sinphi+db*r2g) &
243  & +ca2*((r2k*(da-dc)+(da*r2gd3+dc*r4gd3))*sinphi+da*r2g))*detinv
244  dpsdpe(m1,m3) = coffd+(ca1*(db-dd)+ca2*(da-dc))*(-r2g+(r2k+r2gd3)*sinphi)*detinv
245  dpsdpe(m2,m1) = coffd+(ca1*((r2k*(dc-da)+(da*r4gd3+dc*r2gd3))*sinphi+dc*r2g) &
246  & +ca2*((r2k*(dd-db)+(db*r4gd3+dd*r2gd3))*sinphi+dd*r2g))*detinv
247  dpsdpe(m2,m2) = cdiag+(ca1*((r2k*(dc-da)-(da*r2gd3+dc*r4gd3))*sinphi-da*r2g) &
248  & +ca2*((r2k*(dd-db)-(db*r2gd3+dd*r4gd3))*sinphi-db*r2g))*detinv
249  dpsdpe(m2,m3) = coffd+(ca1*(dc-da)+ca2*(dd-db))*(-r2g+(r2k+r2gd3)*sinphi)*detinv
250  dpsdpe(m3,m1) = coffd+ca3*((r2k*(-db+dd+da-dc)+(db-da)*r4gd3+(dd-dc)*r2gd3)*sinphi+(dd-dc)*r2g)*detinv
251  dpsdpe(m3,m2) = coffd+ca3*((r2k*(-db+dd+da-dc)+(da-db)*r2gd3+(dd-dc)*r4gd3)*sinphi+(da-db)*r2g)*detinv
252  dpsdpe(m3,m3) = cdiag+ca3*(-db+dd+da-dc)*(-r2g+(r2k+r2gd3)*sinphi)*detinv
253  endif
254  endif
255  ! compute principal elastic strain from principal stress
256  prnstra(1) = (prnstre(1)-poisson*(prnstre(2)+prnstre(3)))/youngs
257  prnstra(2) = (prnstre(2)-poisson*(prnstre(1)+prnstre(3)))/youngs
258  prnstra(3) = (prnstre(3)-poisson*(prnstre(1)+prnstre(2)))/youngs
259  call deriv_general_iso_tensor_func_3d(dpsdpe, d, prnprj, prnstra, prnstre)
260  end subroutine calelastoplasticmatrix_mc
261 
263  subroutine calelastoplasticmatrix_dp( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
264  type( tmaterial ), intent(in) :: matl
265  integer, intent(in) :: secttype
266  real(kind=kreal), intent(in) :: stress(6)
267  real(kind=kreal), intent(in) :: extval(:)
268  real(kind=kreal), intent(in) :: plstrain
269  integer, intent(in) :: istat
270  real(kind=kreal), intent(out) :: d(:,:)
271  real(kind=kreal), intent(in) :: temperature
272  integer(kind=kint), intent(in) :: hdflag
273 
274  integer :: i,j
275  real(kind=kreal) :: dum, a(6), dlambda, g, k
276  real(kind=kreal) :: j1,j2, eta, xi, etabar, harden, devia(6)
277  real(kind=kreal) :: alpha, beta, c1, c2, c3, c4, ca, devia_norm
278 
279  if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
280 
281  call calelasticmatrix( matl, secttype, d, temperature, hdflag=hdflag )
282  if( istat == dp_elastic ) return ! elastic state
283  if( hdflag == 2 ) return
284 
285  harden = calhardencoeff( matl, extval(1), temperature )
286 
287  g = d(4,4)
288  k = d(1,1)-(4.d0/3.d0)*g
289 
290  eta = matl%variables(m_plconst3)
291  xi = matl%variables(m_plconst4)
292  etabar = matl%variables(m_plconst5)
293 
294  if( istat==dp_plastic_surf ) then
295  j1 = (stress(1)+stress(2)+stress(3))
296  devia(1:3) = stress(1:3)-j1/3.d0
297  devia(4:6) = stress(4:6)
298  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
299  dot_product( devia(4:6), devia(4:6) )
300 
301  devia_norm = sqrt(2.d0*j2)
302  a(1:6) = devia(1:6)/devia_norm
303  dlambda = extval(1)-plstrain
304  ca = 1.d0 / (g + k*eta*etabar + xi*xi*harden)
305  dum = sqrt(2.d0)*devia_norm
306  c1 = 4.d0*g*g*dlambda/dum
307  c2 = 2.d0*g*(2.d0*g*dlambda/dum - g*ca)
308  c3 = sqrt(2.d0)*g*ca*k
309  c4 = k*k*eta*etabar*ca
310  do j=1,6
311  do i=1,6
312  d(i,j) = d(i,j) - c1*id(i,j) + c2*a(i)*a(j) &
313  - c3*(eta*a(i)*i2(j) + etabar*i2(i)*a(j)) &
314  - c4*i2(i)*i2(j)
315  enddo
316  enddo
317  else ! istat==DP_PLASTIC_APEX
318  alpha = xi/etabar
319  beta = xi/eta
320  c1 = k*(1.d0 - k/(k + alpha*beta*harden))
321  do j=1,6
322  do i=1,6
323  d(i,j) = c1*i2(i)*i2(j)
324  enddo
325  enddo
326  endif
327  end subroutine calelastoplasticmatrix_dp
328 
330  real(kind=kreal) function calhardencoeff( matl, pstrain, temp )
331  type( tmaterial ), intent(in) :: matl
332  real(kind=kreal), intent(in) :: pstrain
333  real(kind=kreal), intent(in) :: temp
334 
335  integer :: htype
336  logical :: ierr
337  real(kind=kreal) :: s0, s1,s2, ef, ina(2)
338 
339  calhardencoeff = -1.d0
340  htype = gethardentype( matl%mtype )
341  select case (htype)
342  case (0) ! Linear hardening
343  calhardencoeff = matl%variables(m_plconst2)
344  case (1) ! Multilinear approximation
345  ina(1) = temp; ina(2)=pstrain
346  call fetch_tablegrad( mc_yield, ina, matl%dict, calhardencoeff, ierr )
347  case (2) ! Swift
348  s0= matl%variables(m_plconst1)
349  s1= matl%variables(m_plconst2)
350  s2= matl%variables(m_plconst3)
351  calhardencoeff = s1*s2*( s0+pstrain )**(s2-1)
352  case (3) ! Ramberg-Osgood
353  s0= matl%variables(m_plconst1)
354  s1= matl%variables(m_plconst2)
355  s2= matl%variables(m_plconst3)
356  ef = calcurryield( matl, pstrain, temp )
357  calhardencoeff = s1*(ef/s1)**(1.d0-s2) /(s0*s2)
358  case(4) ! Prager
359  calhardencoeff = 0.d0
360  case(5) ! Prager+linear
361  calhardencoeff = matl%variables(m_plconst2)
362  end select
363  end function
364 
366  real(kind=kreal) function calkinematicharden( matl, pstrain )
367  type( tmaterial ), intent(in) :: matl
368  real(kind=kreal), intent(in) :: pstrain
369 
370  integer :: htype
371  htype = gethardentype( matl%mtype )
372  select case (htype)
373  case(4, 5) ! Prager
374  calkinematicharden = matl%variables(m_plconst3)
375  case default
376  calkinematicharden = 0.d0
377  end select
378  end function
379 
381  real(kind=kreal) function calcurrkinematic( matl, pstrain )
382  type( tmaterial ), intent(in) :: matl
383  real(kind=kreal), intent(in) :: pstrain
384 
385  integer :: htype
386  htype = gethardentype( matl%mtype )
387  select case (htype)
388  case(4, 5) ! Prager
389  calcurrkinematic = matl%variables(m_plconst3)*pstrain
390  case default
391  calcurrkinematic = 0.d0
392  end select
393  end function
394 
396  real(kind=kreal) function calcurryield( matl, pstrain, temp )
397  type( tmaterial ), intent(in) :: matl
398  real(kind=kreal), intent(in) :: pstrain
399  real(kind=kreal), intent(in) :: temp
400 
401  integer :: htype
402  real(kind=kreal) :: s0, s1,s2, ina(2), outa(1)
403  logical :: ierr
404  calcurryield = -1.d0
405  htype = gethardentype( matl%mtype )
406 
407  select case (htype)
408  case (0, 5) ! Linear hardening, Linear+Parger hardening
409  calcurryield = matl%variables(m_plconst1)+matl%variables(m_plconst2)*pstrain
410  case (1) ! Multilinear approximation
411  ina(1) = temp; ina(2)=pstrain
412  call fetch_tabledata(mc_yield, matl%dict, outa, ierr, ina)
413  if( ierr ) stop "Fail to get yield stress!"
414  calcurryield = outa(1)
415  case (2) ! Swift
416  s0= matl%variables(m_plconst1)
417  s1= matl%variables(m_plconst2)
418  s2= matl%variables(m_plconst3)
419  calcurryield = s1*( s0+pstrain )**s2
420  case (3) ! Ramberg-Osgood
421  s0= matl%variables(m_plconst1)
422  s1= matl%variables(m_plconst2)
423  s2= matl%variables(m_plconst3)
424  if( pstrain<=s0 ) then
425  calcurryield = s1
426  else
427  calcurryield = s1*( pstrain/s0 )**(1.d0/s2)
428  endif
429  case (4) ! Parger hardening
430  calcurryield = matl%variables(m_plconst1)
431  end select
432  end function
433 
435  subroutine backwardeuler( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag )
436  type( tmaterial ), intent(in) :: matl
437  real(kind=kreal), intent(inout) :: stress(6)
438  real(kind=kreal), intent(in) :: plstrain
439  integer, intent(inout) :: istat
440  real(kind=kreal), intent(inout) :: fstat(:)
441  real(kind=kreal), intent(inout) :: plpotential
442  real(kind=kreal), intent(in) :: temp
443  integer(kind=kint), intent(in), optional :: hdflag
444 
445  integer :: ytype, hdflag_in
446 
447  hdflag_in = 0
448  if( present(hdflag) ) hdflag_in = hdflag
449 
450  ytype = getyieldfunction( matl%mtype )
451  select case (ytype)
452  case (0)
453  call backwardeuler_vm( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag_in )
454  case (1)
455  call backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag_in )
456  case (2)
457  call backwardeuler_dp( matl, stress, plstrain, istat, fstat, temp, hdflag_in )
458  case (3)
459  call ubackwardeuler( matl%variables, stress, plstrain, istat, fstat, temp, hdflag_in )
460  end select
461  end subroutine backwardeuler
462 
464  subroutine backwardeuler_vm( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag )
465  type( tmaterial ), intent(in) :: matl
466  real(kind=kreal), intent(inout) :: stress(6)
467  real(kind=kreal), intent(in) :: plstrain
468  integer, intent(inout) :: istat
469  real(kind=kreal), intent(inout) :: fstat(:)
470  real(kind=kreal), intent(inout) :: plpotential
471  real(kind=kreal), intent(in) :: temp
472  integer(kind=kint), intent(in) :: hdflag
473 
474  real(kind=kreal), parameter :: tol =1.d-8
475  integer, parameter :: maxiter = 10
476  real(kind=kreal) :: dlambda, f
477  integer :: i
478  real(kind=kreal) :: youngs, poisson, pstrain, ina(1), ee(2)
479  real(kind=kreal) :: j1, j2, h, kh, kk, dd, eqvs, yd, g, k, devia(6)
480  logical :: kinematic, ierr
481  real(kind=kreal) :: betan, back(6)
482 
483  plpotential = 0.d0
484  kinematic = iskinematicharden( matl%mtype )
485  if( kinematic ) back(1:6) = fstat(8:13)
486 
487  j1 = (stress(1)+stress(2)+stress(3))
488  devia(1:3) = stress(1:3)-j1/3.d0
489  devia(4:6) = stress(4:6)
490  if( kinematic ) devia = devia-back
491  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
492  dot_product( devia(4:6), devia(4:6) )
493 
494  eqvs = dsqrt( 3.d0*j2 )
495  yd = calcurryield( matl, plstrain, temp )
496  f = eqvs - yd
497 
498  if( abs(f/yd)<tol ) then ! yielded
499  istat = vm_plastic
500  return
501  elseif( f<0.d0 ) then ! not yielded or unloading
502  istat = vm_elastic
503  return
504  endif
505  if( hdflag == 2 ) return
506 
507  istat = vm_plastic ! yielded
508  kh = 0.d0; kk=0.d0; betan=0.d0; back(:)=0.d0
509  if( kinematic ) betan = calcurrkinematic( matl, plstrain )
510 
511  ina(1) = temp
512  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
513  if( ierr ) then
514  stop " fail to fetch young's modulus in elastoplastic calculation"
515  else
516  youngs = ee(1)
517  poisson = ee(2)
518  endif
519  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
520  g = youngs/ ( 2.d0*(1.d0+poisson) )
521  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
522 
523  dlambda = 0.d0
524  pstrain = plstrain
525 
526  do i=1,maxiter
527  h= calhardencoeff( matl, pstrain, temp )
528  if( kinematic ) then
529  kh = calkinematicharden( matl, pstrain )
530  endif
531  dd= 3.d0*g+h+kh
532  dlambda = dlambda+f/dd
533  if( dlambda<0.d0 ) then
534  dlambda = 0.d0
535  pstrain = plstrain
536  istat=vm_elastic; exit
537  endif
538  pstrain = plstrain+dlambda
539  yd = calcurryield( matl, pstrain, temp )
540  if( kinematic ) then
541  kk = calcurrkinematic( matl, pstrain )
542  endif
543  f = eqvs-3.d0*g*dlambda-yd -(kk-betan)
544  if( abs(f/yd)<tol ) exit
545  ! if( i==MAXITER ) then
546  ! stop 'ERROR: BackwardEuler_VM: convergence failure'
547  ! endif
548  enddo
549  if( kinematic ) then
550  kk = calcurrkinematic( matl, pstrain )
551  fstat(2:7) = back(:)+(kk-betan)*devia(:)/eqvs
552  endif
553  devia(:) = (1.d0-3.d0*dlambda*g/eqvs)*devia(:)
554  stress(1:3) = devia(1:3)+j1/3.d0
555  stress(4:6) = devia(4:6)
556  stress(:)= stress(:)+back(:)
557 
558  fstat(1) = pstrain
559 
560  h= calhardencoeff( matl, pstrain, temp ) !a
561  yd = calcurryield( matl, plstrain, temp ) !b
562  plpotential = -0.5d0*(eqvs-yd)*(eqvs-yd)/(h+3.d0*g)
563 
564  end subroutine backwardeuler_vm
565 
567  subroutine backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag )
568  use m_utilities, only : eigen3
569  type( tmaterial ), intent(in) :: matl
570  real(kind=kreal), intent(inout) :: stress(6)
571  real(kind=kreal), intent(in) :: plstrain
572  integer, intent(inout) :: istat
573  real(kind=kreal), intent(inout) :: fstat(:)
574  real(kind=kreal), intent(in) :: temp
575  integer(kind=kint), intent(in) :: hdflag
576 
577  real(kind=kreal), parameter :: tol =1.d-8
578  integer, parameter :: maxiter = 10
579  real(kind=kreal) :: dlambda, f, mat(3,3)
580  integer :: i, m1, m2, m3
581  real(kind=kreal) :: youngs, poisson, pstrain, ina(1), ee(2)
582  real(kind=kreal) :: h, dd, eqvs, cohe, g, k
583  real(kind=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3)
584  real(kind=kreal) :: phi, psi, trialprn(3)
585  logical :: ierr
586  real(kind=kreal) :: c1, c2, cs1, cs2, cs3
587  real(kind=kreal) :: sinphi, cosphi, sinpsi, sphsps, r2cosphi, r4cos2phi, cotphi
588  real(kind=kreal) :: da, db, dc, depv, detinv, dlambdb, dum, eps, eqvsb, fb
589  real(kind=kreal) :: pt, p, resid
590 
591  phi = matl%variables(m_plconst3)
592  psi = matl%variables(m_plconst4)
593  sinphi = sin(phi)
594  cosphi = cos(phi)
595  r2cosphi = 2.d0*cosphi
596 
597  call eigen3( stress, prnstre, prnprj )
598  trialprn = prnstre
599  m1 = maxloc( prnstre, 1 )
600  m3 = minloc( prnstre, 1 )
601  if( m1 == m3 ) then
602  m1 = 1; m2 = 2; m3 = 3
603  else
604  m2 = 6 - (m1 + m3)
605  endif
606 
607  eqvs = prnstre(m1)-prnstre(m3) + (prnstre(m1)+prnstre(m3))*sinphi
608  cohe = calcurryield( matl, plstrain, temp )
609  f = eqvs - r2cosphi*cohe
610 
611  if( abs(f/cohe)<tol ) then ! yielded
612  istat = mc_plastic_surf
613  return
614  elseif( f<0.d0 ) then ! not yielded or unloading
615  istat = mc_elastic
616  return
617  endif
618  if( hdflag == 2 ) return
619 
620  istat = mc_plastic_surf ! yielded
621 
622  ina(1) = temp
623  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
624  if( ierr ) then
625  stop " fail to fetch young's modulus in elastoplastic calculation"
626  else
627  youngs = ee(1)
628  poisson = ee(2)
629  endif
630  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
631  g = youngs/ ( 2.d0*(1.d0+poisson) )
632  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
633 
634  dlambda = 0.d0
635  pstrain = plstrain
636 
637  sinpsi = sin(psi)
638  sphsps = sinphi*sinpsi
639  r4cos2phi = r2cosphi*r2cosphi
640  c1 = 4.d0*(g*(1.d0+sphsps/3.d0)+k*sphsps)
641  do i=1,maxiter
642  h= calhardencoeff( matl, pstrain, temp )
643  dd= c1 + r4cos2phi*h
644  dlambda = dlambda+f/dd
645  if( r2cosphi*dlambda<0.d0 ) then
646  if( cosphi==0.d0 ) stop "Math error in return mapping"
647  dlambda = 0.d0
648  pstrain = plstrain
649  istat = mc_elastic; exit
650  endif
651  pstrain = plstrain + r2cosphi*dlambda
652  cohe = calcurryield( matl, pstrain, temp )
653  f = eqvs - c1*dlambda - r2cosphi*cohe
654  if( abs(f/cohe)<tol ) exit
655  ! if( i==MAXITER ) then
656  ! stop 'ERROR: BackwardEuler_MC: convergence failure'
657  ! endif
658  enddo
659  cs1 =2.d0*g*(1.d0+sinpsi/3.d0) + 2.d0*k*sinpsi
660  cs2 =(4.d0*g/3.d0-2.d0*k)*sinpsi
661  cs3 =2.d0*g*(1.d0-sinpsi/3.d0) - 2.d0*k*sinpsi
662  prnstre(m1) = prnstre(m1)-cs1*dlambda
663  prnstre(m2) = prnstre(m2)+cs2*dlambda
664  prnstre(m3) = prnstre(m3)+cs3*dlambda
665  eps = (abs(prnstre(m1))+abs(prnstre(m2))+abs(prnstre(m3)))*tol
666  if( prnstre(m1) < prnstre(m2)-eps .or. prnstre(m2) < prnstre(m3)-eps ) then
667  ! return mapping to EDGE
668  prnstre = trialprn
669  dlambda = 0.d0
670  dlambdb = 0.d0
671  if( (1.d0-sinpsi)*prnstre(m1) - 2*prnstre(m2) + (1.d0+sinpsi)*prnstre(m3) > 0) then
672  istat = mc_plastic_right
673  eqvsb = prnstre(m1)-prnstre(m2) + (prnstre(m1)+prnstre(m2))*sinphi
674  c2 = 2.d0*g*(1.d0+sinphi+sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
675  else
676  istat = mc_plastic_left
677  eqvsb = prnstre(m2)-prnstre(m3) + (prnstre(m2)+prnstre(m3))*sinphi
678  c2 = 2.d0*g*(1.d0-sinphi-sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
679  endif
680  cohe = calcurryield( matl, plstrain, temp )
681  f = eqvs - r2cosphi*cohe
682  fb = eqvsb - r2cosphi*cohe
683  pstrain = plstrain
684  do i=1,maxiter
685  h= calhardencoeff( matl, pstrain, temp )
686  dum = r4cos2phi*h
687  da = c1 + dum
688  db = c2 + dum
689  dc = db
690  dd = da
691  detinv = 1.d0/(da*dd-db*dc)
692  dlambda = dlambda + detinv*( dd*f - db*fb)
693  dlambdb = dlambdb + detinv*(-dc*f + da*fb)
694  pstrain = plstrain + r2cosphi*(dlambda+dlambdb)
695  cohe = calcurryield( matl, pstrain, temp )
696  f = eqvs - c1*dlambda - c2*dlambdb - r2cosphi*cohe
697  fb = eqvsb - c2*dlambda - c1*dlambdb - r2cosphi*cohe
698  if( (abs(f)+abs(fb))/(abs(eqvs)+abs(eqvsb)) < tol ) exit
699  ! if( i==MAXITER ) then
700  ! stop 'ERROR: BackwardEuler_MC: convergence failure(2)'
701  ! endif
702  enddo
703  if( istat==mc_plastic_right ) then
704  prnstre(m1) = prnstre(m1)-cs1*(dlambda+dlambdb)
705  prnstre(m2) = prnstre(m2)+cs2*dlambda+cs3*dlambdb
706  prnstre(m3) = prnstre(m3)+cs3*dlambda+cs2*dlambdb
707  else
708  prnstre(m1) = prnstre(m1)-cs1*dlambda+cs2*dlambdb
709  prnstre(m2) = prnstre(m2)+cs2*dlambda-cs1*dlambdb
710  prnstre(m3) = prnstre(m3)+cs3*(dlambda+dlambdb)
711  endif
712  eps = (abs(prnstre(m1))+abs(prnstre(m2))+abs(prnstre(m3)))*tol
713  if( prnstre(m1) < prnstre(m2)-eps .or. prnstre(m2) < prnstre(m3)-eps ) then
714  ! return mapping to APEX
715  prnstre = trialprn
716  istat = mc_plastic_apex
717  if( sinphi==0.d0 ) stop 'ERROR: BackwardEuler_MC: phi==0.0'
718  if( sinpsi==0.d0 ) stop 'ERROR: BackwardEuler_MC: psi==0.0'
719  depv = 0.d0
720  cohe = calcurryield( matl, plstrain, temp )
721  cotphi = cosphi/sinphi
722  pt = (stress(1)+stress(2)+stress(3))/3.d0
723  resid = cotphi*cohe - pt
724  pstrain = plstrain
725  do i=1,maxiter
726  h= calhardencoeff( matl, pstrain, temp )
727  dd= cosphi*cotphi*h/sinpsi + k
728  depv = depv - resid/dd
729  pstrain = plstrain + cosphi*depv/sinpsi
730  cohe = calcurryield( matl, pstrain,temp )
731  p = pt-k*depv
732  resid = cotphi*cohe-p
733  if( abs(resid/cohe)<tol ) exit
734  ! if( i==MAXITER ) then
735  ! stop 'ERROR: BackwardEuler_MC: convergence failure(3)'
736  ! endif
737  enddo
738  prnstre(m1) = p
739  prnstre(m2) = p
740  prnstre(m3) = p
741  endif
742  endif
743  tstre(:,:) = 0.d0
744  tstre(1,1)= prnstre(1); tstre(2,2)=prnstre(2); tstre(3,3)=prnstre(3)
745  mat= matmul( prnprj, tstre )
746  mat= matmul( mat, transpose(prnprj) )
747  stress(1) = mat(1,1)
748  stress(2) = mat(2,2)
749  stress(3) = mat(3,3)
750  stress(4) = mat(1,2)
751  stress(5) = mat(2,3)
752  stress(6) = mat(3,1)
753 
754  fstat(1) = pstrain
755  end subroutine backwardeuler_mc
756 
758  subroutine backwardeuler_dp( matl, stress, plstrain, istat, fstat, temp, hdflag )
759  type( tmaterial ), intent(in) :: matl
760  real(kind=kreal), intent(inout) :: stress(6)
761  real(kind=kreal), intent(in) :: plstrain
762  integer, intent(inout) :: istat
763  real(kind=kreal), intent(inout) :: fstat(:)
764  real(kind=kreal), intent(in) :: temp
765  integer(kind=kint), intent(in) :: hdflag
766 
767  real(kind=kreal), parameter :: tol =1.d-8
768  integer, parameter :: maxiter = 10
769  real(kind=kreal) :: dlambda, f
770  integer :: i
771  real(kind=kreal) :: youngs, poisson, pstrain, xi, ina(1), ee(2)
772  real(kind=kreal) :: j1,j2,h, dd, eqvst, eqvs, cohe, g, k, devia(6), eta, etabar, pt, p
773  logical :: ierr
774  real(kind=kreal) :: alpha, beta, depv, factor, resid
775 
776  eta = matl%variables(m_plconst3)
777  xi = matl%variables(m_plconst4)
778  etabar = matl%variables(m_plconst5)
779 
780  j1 = (stress(1)+stress(2)+stress(3))
781  pt = j1/3.d0
782  devia(1:3) = stress(1:3)-pt
783  devia(4:6) = stress(4:6)
784  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
785  dot_product( devia(4:6), devia(4:6) )
786 
787  eqvst = sqrt(j2)
788  cohe = calcurryield( matl, plstrain, temp )
789  f = eqvst + eta*pt - xi*cohe
790 
791  if( abs(f/cohe)<tol ) then ! yielded
792  istat = dp_plastic_surf
793  return
794  elseif( f<0.d0 ) then ! not yielded or unloading
795  istat = dp_elastic
796  return
797  endif
798  if( hdflag == 2 ) return
799 
800  istat = dp_plastic_surf
801 
802  ina(1) = temp
803  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
804  if( ierr ) then
805  stop " fail to fetch young's modulus in elastoplastic calculation"
806  else
807  youngs = ee(1)
808  poisson = ee(2)
809  endif
810  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
811  g = youngs/ ( 2.d0*(1.d0+poisson) )
812  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
813 
814  dlambda = 0.d0
815  pstrain = plstrain
816 
817  do i=1,maxiter
818  h= calhardencoeff( matl, pstrain, temp )
819  dd= g+k*etabar*eta+h*xi*xi
820  dlambda = dlambda+f/dd
821  if( xi*dlambda<0.d0 ) then
822  if( xi==0.d0 ) stop "Math error in return mapping"
823  dlambda = 0.d0
824  pstrain = plstrain
825  istat=0; exit
826  endif
827  pstrain = plstrain+xi*dlambda
828  cohe = calcurryield( matl, pstrain, temp )
829  eqvs = eqvst-g*dlambda
830  p = pt-k*etabar*dlambda
831  f = eqvs + eta*p- xi*cohe
832  if( abs(f/cohe)<tol ) exit
833  ! if( i==MAXITER ) then
834  ! stop 'ERROR: BackwardEuler_DP: convergence failure'
835  ! endif
836  enddo
837  if( eqvs>=0.d0 ) then ! converged
838  factor = 1.d0-g*dlambda/eqvst
839  else ! return mapping to APEX
840  istat = dp_plastic_apex
841  if( eta==0.d0 ) stop 'ERROR: BackwardEuler_DP: eta==0.0'
842  if( etabar==0.d0 ) stop 'ERROR: BackwardEuler_DP: etabar==0.0'
843  alpha = xi/etabar
844  beta = xi/eta
845  depv=0.d0
846  pstrain = plstrain
847  cohe = calcurryield( matl, pstrain, temp )
848  resid = beta*cohe - pt
849  do i=1,maxiter
850  h= calhardencoeff( matl, pstrain, temp )
851  dd= alpha*beta*h + k
852  depv = depv - resid/dd
853  pstrain = plstrain+alpha*depv
854  cohe = calcurryield( matl, pstrain, temp )
855  p = pt-k*depv
856  resid = beta*cohe - p
857  if( abs(resid/cohe)<tol ) then
858  dlambda=depv/etabar
859  factor=0.d0
860  exit
861  endif
862  ! if( i==MAXITER ) then
863  ! stop 'ERROR: BackwardEuler_DP: convergence failure(2)'
864  ! endif
865  enddo
866  endif
867  devia(:) = factor*devia(:)
868  stress(1:3) = devia(1:3)+p
869  stress(4:6) = devia(4:6)
870 
871  fstat(1) = pstrain
872  end subroutine backwardeuler_dp
873 
875  subroutine updateepstate( gauss )
876  use mmechgauss
877  type(tgaussstatus), intent(inout) :: gauss ! status of curr gauss point
878  gauss%plstrain= gauss%fstatus(1)
879  if(iskinematicharden(gauss%pMaterial%mtype)) then
880  gauss%fstatus(8:13) =gauss%fstatus(2:7)
881  endif
882  end subroutine
883 
884 end module m_elastoplastic
m_utilities::deriv_general_iso_tensor_func_3d
subroutine deriv_general_iso_tensor_func_3d(dpydpx, dydx, eigv, px, py)
Compute derivative of a general isotropic tensor function of one tensor.
Definition: utilities.f90:552
m_elastoplastic::backwardeuler
subroutine, public backwardeuler(matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag)
This subroutine does backward-Euler return calculation.
Definition: Elastoplastic.f90:436
mmaterial::mc_yield
character(len=dict_key_length) mc_yield
Definition: material.f90:142
mmaterial::mc_isoelastic
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:140
muyield::ubackwardeuler
subroutine, public ubackwardeuler(matl, stress, plstrain, istat, fstat, temp, hdflag)
This subroutine does backward-Euler return calculation.
Definition: uyield.f90:45
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
mmaterial::iskinematicharden
logical function iskinematicharden(mtype)
If it is a kinematic hardening material?
Definition: material.f90:343
mmaterial::m_plconst2
integer(kind=kint), parameter m_plconst2
Definition: material.f90:99
muyield
This subroutine read in used-defined material properties tangent.
Definition: uyield.f90:7
mmaterial::m_plconst5
integer(kind=kint), parameter m_plconst5
Definition: material.f90:102
mmaterial::m_plconst3
integer(kind=kint), parameter m_plconst3
Definition: material.f90:100
m_fstr::eps
real(kind=kreal) eps
Definition: m_fstr.f90:142
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
m_elastoplastic
This module provide functions for elastoplastic calculation.
Definition: Elastoplastic.f90:6
m_elastoplastic::updateepstate
subroutine, public updateepstate(gauss)
Clear elatoplastic state.
Definition: Elastoplastic.f90:876
m_utilities::eigen3
subroutine eigen3(tensor, eigval, princ)
Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using Jacobi iteration adapted from numeri...
Definition: utilities.f90:106
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:84
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_elasticlinear
This module provides functions for elastic material.
Definition: ElasticLinear.f90:6
muyield::uelastoplasticmatrix
subroutine, public uelastoplasticmatrix(matl, stress, istat, fstat, plstrain, D, temp, hdflag)
This subroutine calculates elastoplastic constitutive relation.
Definition: uyield.f90:32
mmaterial::m_plconst1
integer(kind=kint), parameter m_plconst1
Definition: material.f90:98
mmaterial::getyieldfunction
integer function getyieldfunction(mtype)
Get type of yield function.
Definition: material.f90:319
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
mmaterial::tmaterial
Structure to manage all material related data.
Definition: material.f90:167
mmaterial::m_plconst4
integer(kind=kint), parameter m_plconst4
Definition: material.f90:101
m_elastoplastic::backwardeuler_dp
subroutine backwardeuler_dp(matl, stress, plstrain, istat, fstat, temp, hdflag)
This subroutine does backward-Euler return calculation for Drucker-Prager.
Definition: Elastoplastic.f90:759
m_elastoplastic::calelastoplasticmatrix_dp
subroutine calelastoplasticmatrix_dp(matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag)
This subroutine calculates elastoplastic constitutive relation.
Definition: Elastoplastic.f90:264
mmaterial::gethardentype
integer function gethardentype(mtype)
Get type of hardening.
Definition: material.f90:331
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
m_elastoplastic::calelastoplasticmatrix
subroutine, public calelastoplasticmatrix(matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag)
This subroutine calculates elastoplastic constitutive relation.
Definition: Elastoplastic.f90:46
m_elasticlinear::calelasticmatrix
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix.
Definition: ElasticLinear.f90:16