FrontISTR  5.9.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
509  if( kinematic ) then
510  betan = calcurrkinematic( matl, plstrain ) ! keep back = alpha_n (loaded above) so it accumulates at fstat update and is restored into stress
511  else
512  back(:)=0.d0
513  endif
514 
515  ina(1) = temp
516  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
517  if( ierr ) then
518  stop " fail to fetch young's modulus in elastoplastic calculation"
519  else
520  youngs = ee(1)
521  poisson = ee(2)
522  endif
523  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
524  g = youngs/ ( 2.d0*(1.d0+poisson) )
525  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
526 
527  dlambda = 0.d0
528  pstrain = plstrain
529 
530  do i=1,maxiter
531  h= calhardencoeff( matl, pstrain, temp )
532  if( kinematic ) then
533  kh = calkinematicharden( matl, pstrain )
534  endif
535  dd= 3.d0*g+h+kh
536  dlambda = dlambda+f/dd
537  if( dlambda<0.d0 ) then
538  dlambda = 0.d0
539  pstrain = plstrain
540  istat=vm_elastic; exit
541  endif
542  pstrain = plstrain+dlambda
543  yd = calcurryield( matl, pstrain, temp )
544  if( kinematic ) then
545  kk = calcurrkinematic( matl, pstrain )
546  endif
547  f = eqvs-3.d0*g*dlambda-yd -(kk-betan)
548  if( abs(f/yd)<tol ) exit
549  ! if( i==MAXITER ) then
550  ! stop 'ERROR: BackwardEuler_VM: convergence failure'
551  ! endif
552  enddo
553  if( kinematic ) then
554  kk = calcurrkinematic( matl, pstrain )
555  fstat(2:7) = back(:)+(kk-betan)*devia(:)/eqvs
556  endif
557  devia(:) = (1.d0-3.d0*dlambda*g/eqvs)*devia(:)
558  stress(1:3) = devia(1:3)+j1/3.d0
559  stress(4:6) = devia(4:6)
560  stress(:)= stress(:)+back(:)
561 
562  fstat(1) = pstrain
563 
564  h= calhardencoeff( matl, pstrain, temp ) !a
565  yd = calcurryield( matl, plstrain, temp ) !b
566  plpotential = -0.5d0*(eqvs-yd)*(eqvs-yd)/(h+3.d0*g)
567 
568  end subroutine backwardeuler_vm
569 
571  subroutine backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag )
572  use m_utilities, only : eigen3
573  type( tmaterial ), intent(in) :: matl
574  real(kind=kreal), intent(inout) :: stress(6)
575  real(kind=kreal), intent(in) :: plstrain
576  integer, intent(inout) :: istat
577  real(kind=kreal), intent(inout) :: fstat(:)
578  real(kind=kreal), intent(in) :: temp
579  integer(kind=kint), intent(in) :: hdflag
580 
581  real(kind=kreal), parameter :: tol =1.d-8
582  integer, parameter :: maxiter = 10
583  real(kind=kreal) :: dlambda, f, mat(3,3)
584  integer :: i, m1, m2, m3
585  real(kind=kreal) :: youngs, poisson, pstrain, ina(1), ee(2)
586  real(kind=kreal) :: h, dd, eqvs, cohe, g, k
587  real(kind=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3)
588  real(kind=kreal) :: phi, psi, trialprn(3)
589  logical :: ierr
590  real(kind=kreal) :: c1, c2, cs1, cs2, cs3
591  real(kind=kreal) :: sinphi, cosphi, sinpsi, sphsps, r2cosphi, r4cos2phi, cotphi
592  real(kind=kreal) :: da, db, dc, depv, detinv, dlambdb, dum, eps, eqvsb, fb
593  real(kind=kreal) :: pt, p, resid
594 
595  phi = matl%variables(m_plconst3)
596  psi = matl%variables(m_plconst4)
597  sinphi = sin(phi)
598  cosphi = cos(phi)
599  r2cosphi = 2.d0*cosphi
600 
601  call eigen3( stress, prnstre, prnprj )
602  trialprn = prnstre
603  m1 = maxloc( prnstre, 1 )
604  m3 = minloc( prnstre, 1 )
605  if( m1 == m3 ) then
606  m1 = 1; m2 = 2; m3 = 3
607  else
608  m2 = 6 - (m1 + m3)
609  endif
610 
611  eqvs = prnstre(m1)-prnstre(m3) + (prnstre(m1)+prnstre(m3))*sinphi
612  cohe = calcurryield( matl, plstrain, temp )
613  f = eqvs - r2cosphi*cohe
614 
615  if( abs(f/cohe)<tol ) then ! yielded
616  istat = mc_plastic_surf
617  return
618  elseif( f<0.d0 ) then ! not yielded or unloading
619  istat = mc_elastic
620  return
621  endif
622  if( hdflag == 2 ) return
623 
624  istat = mc_plastic_surf ! yielded
625 
626  ina(1) = temp
627  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
628  if( ierr ) then
629  stop " fail to fetch young's modulus in elastoplastic calculation"
630  else
631  youngs = ee(1)
632  poisson = ee(2)
633  endif
634  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
635  g = youngs/ ( 2.d0*(1.d0+poisson) )
636  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
637 
638  dlambda = 0.d0
639  pstrain = plstrain
640 
641  sinpsi = sin(psi)
642  sphsps = sinphi*sinpsi
643  r4cos2phi = r2cosphi*r2cosphi
644  c1 = 4.d0*(g*(1.d0+sphsps/3.d0)+k*sphsps)
645  do i=1,maxiter
646  h= calhardencoeff( matl, pstrain, temp )
647  dd= c1 + r4cos2phi*h
648  dlambda = dlambda+f/dd
649  if( r2cosphi*dlambda<0.d0 ) then
650  if( cosphi==0.d0 ) stop "Math error in return mapping"
651  dlambda = 0.d0
652  pstrain = plstrain
653  istat = mc_elastic; exit
654  endif
655  pstrain = plstrain + r2cosphi*dlambda
656  cohe = calcurryield( matl, pstrain, temp )
657  f = eqvs - c1*dlambda - r2cosphi*cohe
658  if( abs(f/cohe)<tol ) exit
659  ! if( i==MAXITER ) then
660  ! stop 'ERROR: BackwardEuler_MC: convergence failure'
661  ! endif
662  enddo
663  cs1 =2.d0*g*(1.d0+sinpsi/3.d0) + 2.d0*k*sinpsi
664  cs2 =(4.d0*g/3.d0-2.d0*k)*sinpsi
665  cs3 =2.d0*g*(1.d0-sinpsi/3.d0) - 2.d0*k*sinpsi
666  prnstre(m1) = prnstre(m1)-cs1*dlambda
667  prnstre(m2) = prnstre(m2)+cs2*dlambda
668  prnstre(m3) = prnstre(m3)+cs3*dlambda
669  eps = (abs(prnstre(m1))+abs(prnstre(m2))+abs(prnstre(m3)))*tol
670  if( prnstre(m1) < prnstre(m2)-eps .or. prnstre(m2) < prnstre(m3)-eps ) then
671  ! return mapping to EDGE
672  prnstre = trialprn
673  dlambda = 0.d0
674  dlambdb = 0.d0
675  if( (1.d0-sinpsi)*prnstre(m1) - 2*prnstre(m2) + (1.d0+sinpsi)*prnstre(m3) > 0) then
676  istat = mc_plastic_right
677  eqvsb = prnstre(m1)-prnstre(m2) + (prnstre(m1)+prnstre(m2))*sinphi
678  c2 = 2.d0*g*(1.d0+sinphi+sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
679  else
680  istat = mc_plastic_left
681  eqvsb = prnstre(m2)-prnstre(m3) + (prnstre(m2)+prnstre(m3))*sinphi
682  c2 = 2.d0*g*(1.d0-sinphi-sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
683  endif
684  cohe = calcurryield( matl, plstrain, temp )
685  f = eqvs - r2cosphi*cohe
686  fb = eqvsb - r2cosphi*cohe
687  pstrain = plstrain
688  do i=1,maxiter
689  h= calhardencoeff( matl, pstrain, temp )
690  dum = r4cos2phi*h
691  da = c1 + dum
692  db = c2 + dum
693  dc = db
694  dd = da
695  detinv = 1.d0/(da*dd-db*dc)
696  dlambda = dlambda + detinv*( dd*f - db*fb)
697  dlambdb = dlambdb + detinv*(-dc*f + da*fb)
698  pstrain = plstrain + r2cosphi*(dlambda+dlambdb)
699  cohe = calcurryield( matl, pstrain, temp )
700  f = eqvs - c1*dlambda - c2*dlambdb - r2cosphi*cohe
701  fb = eqvsb - c2*dlambda - c1*dlambdb - r2cosphi*cohe
702  if( (abs(f)+abs(fb))/(abs(eqvs)+abs(eqvsb)) < tol ) exit
703  ! if( i==MAXITER ) then
704  ! stop 'ERROR: BackwardEuler_MC: convergence failure(2)'
705  ! endif
706  enddo
707  if( istat==mc_plastic_right ) then
708  prnstre(m1) = prnstre(m1)-cs1*(dlambda+dlambdb)
709  prnstre(m2) = prnstre(m2)+cs2*dlambda+cs3*dlambdb
710  prnstre(m3) = prnstre(m3)+cs3*dlambda+cs2*dlambdb
711  else
712  prnstre(m1) = prnstre(m1)-cs1*dlambda+cs2*dlambdb
713  prnstre(m2) = prnstre(m2)+cs2*dlambda-cs1*dlambdb
714  prnstre(m3) = prnstre(m3)+cs3*(dlambda+dlambdb)
715  endif
716  eps = (abs(prnstre(m1))+abs(prnstre(m2))+abs(prnstre(m3)))*tol
717  if( prnstre(m1) < prnstre(m2)-eps .or. prnstre(m2) < prnstre(m3)-eps ) then
718  ! return mapping to APEX
719  prnstre = trialprn
720  istat = mc_plastic_apex
721  if( sinphi==0.d0 ) stop 'ERROR: BackwardEuler_MC: phi==0.0'
722  if( sinpsi==0.d0 ) stop 'ERROR: BackwardEuler_MC: psi==0.0'
723  depv = 0.d0
724  cohe = calcurryield( matl, plstrain, temp )
725  cotphi = cosphi/sinphi
726  pt = (stress(1)+stress(2)+stress(3))/3.d0
727  resid = cotphi*cohe - pt
728  pstrain = plstrain
729  do i=1,maxiter
730  h= calhardencoeff( matl, pstrain, temp )
731  dd= cosphi*cotphi*h/sinpsi + k
732  depv = depv - resid/dd
733  pstrain = plstrain + cosphi*depv/sinpsi
734  cohe = calcurryield( matl, pstrain,temp )
735  p = pt-k*depv
736  resid = cotphi*cohe-p
737  if( abs(resid/cohe)<tol ) exit
738  ! if( i==MAXITER ) then
739  ! stop 'ERROR: BackwardEuler_MC: convergence failure(3)'
740  ! endif
741  enddo
742  prnstre(m1) = p
743  prnstre(m2) = p
744  prnstre(m3) = p
745  endif
746  endif
747  tstre(:,:) = 0.d0
748  tstre(1,1)= prnstre(1); tstre(2,2)=prnstre(2); tstre(3,3)=prnstre(3)
749  mat= matmul( prnprj, tstre )
750  mat= matmul( mat, transpose(prnprj) )
751  stress(1) = mat(1,1)
752  stress(2) = mat(2,2)
753  stress(3) = mat(3,3)
754  stress(4) = mat(1,2)
755  stress(5) = mat(2,3)
756  stress(6) = mat(3,1)
757 
758  fstat(1) = pstrain
759  end subroutine backwardeuler_mc
760 
762  subroutine backwardeuler_dp( matl, stress, plstrain, istat, fstat, temp, hdflag )
763  type( tmaterial ), intent(in) :: matl
764  real(kind=kreal), intent(inout) :: stress(6)
765  real(kind=kreal), intent(in) :: plstrain
766  integer, intent(inout) :: istat
767  real(kind=kreal), intent(inout) :: fstat(:)
768  real(kind=kreal), intent(in) :: temp
769  integer(kind=kint), intent(in) :: hdflag
770 
771  real(kind=kreal), parameter :: tol =1.d-8
772  integer, parameter :: maxiter = 10
773  real(kind=kreal) :: dlambda, f
774  integer :: i
775  real(kind=kreal) :: youngs, poisson, pstrain, xi, ina(1), ee(2)
776  real(kind=kreal) :: j1,j2,h, dd, eqvst, eqvs, cohe, g, k, devia(6), eta, etabar, pt, p
777  logical :: ierr
778  real(kind=kreal) :: alpha, beta, depv, factor, resid
779 
780  eta = matl%variables(m_plconst3)
781  xi = matl%variables(m_plconst4)
782  etabar = matl%variables(m_plconst5)
783 
784  j1 = (stress(1)+stress(2)+stress(3))
785  pt = j1/3.d0
786  devia(1:3) = stress(1:3)-pt
787  devia(4:6) = stress(4:6)
788  j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) + &
789  dot_product( devia(4:6), devia(4:6) )
790 
791  eqvst = sqrt(j2)
792  cohe = calcurryield( matl, plstrain, temp )
793  f = eqvst + eta*pt - xi*cohe
794 
795  if( abs(f/cohe)<tol ) then ! yielded
796  istat = dp_plastic_surf
797  return
798  elseif( f<0.d0 ) then ! not yielded or unloading
799  istat = dp_elastic
800  return
801  endif
802  if( hdflag == 2 ) return
803 
804  istat = dp_plastic_surf
805 
806  ina(1) = temp
807  call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
808  if( ierr ) then
809  stop " fail to fetch young's modulus in elastoplastic calculation"
810  else
811  youngs = ee(1)
812  poisson = ee(2)
813  endif
814  if( youngs==0.d0 ) stop "YOUNG's ratio==0"
815  g = youngs/ ( 2.d0*(1.d0+poisson) )
816  k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
817 
818  dlambda = 0.d0
819  pstrain = plstrain
820 
821  do i=1,maxiter
822  h= calhardencoeff( matl, pstrain, temp )
823  dd= g+k*etabar*eta+h*xi*xi
824  dlambda = dlambda+f/dd
825  if( xi*dlambda<0.d0 ) then
826  if( xi==0.d0 ) stop "Math error in return mapping"
827  dlambda = 0.d0
828  pstrain = plstrain
829  istat=0; exit
830  endif
831  pstrain = plstrain+xi*dlambda
832  cohe = calcurryield( matl, pstrain, temp )
833  eqvs = eqvst-g*dlambda
834  p = pt-k*etabar*dlambda
835  f = eqvs + eta*p- xi*cohe
836  if( abs(f/cohe)<tol ) exit
837  ! if( i==MAXITER ) then
838  ! stop 'ERROR: BackwardEuler_DP: convergence failure'
839  ! endif
840  enddo
841  if( eqvs>=0.d0 ) then ! converged
842  factor = 1.d0-g*dlambda/eqvst
843  else ! return mapping to APEX
844  istat = dp_plastic_apex
845  if( eta==0.d0 ) stop 'ERROR: BackwardEuler_DP: eta==0.0'
846  if( etabar==0.d0 ) stop 'ERROR: BackwardEuler_DP: etabar==0.0'
847  alpha = xi/etabar
848  beta = xi/eta
849  depv=0.d0
850  pstrain = plstrain
851  cohe = calcurryield( matl, pstrain, temp )
852  resid = beta*cohe - pt
853  do i=1,maxiter
854  h= calhardencoeff( matl, pstrain, temp )
855  dd= alpha*beta*h + k
856  depv = depv - resid/dd
857  pstrain = plstrain+alpha*depv
858  cohe = calcurryield( matl, pstrain, temp )
859  p = pt-k*depv
860  resid = beta*cohe - p
861  if( abs(resid/cohe)<tol ) then
862  dlambda=depv/etabar
863  factor=0.d0
864  exit
865  endif
866  ! if( i==MAXITER ) then
867  ! stop 'ERROR: BackwardEuler_DP: convergence failure(2)'
868  ! endif
869  enddo
870  endif
871  devia(:) = factor*devia(:)
872  stress(1:3) = devia(1:3)+p
873  stress(4:6) = devia(4:6)
874 
875  fstat(1) = pstrain
876  end subroutine backwardeuler_dp
877 
879  subroutine updateepstate( gauss )
880  use mmechgauss
881  type(tgaussstatus), intent(inout) :: gauss ! status of curr gauss point
882  gauss%plstrain= gauss%fstatus(1)
883  if(iskinematicharden(gauss%pMaterial%mtype)) then
884  gauss%fstatus(8:13) =gauss%fstatus(2:7)
885  endif
886  end subroutine
887 
888 end module m_elastoplastic
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
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
m_elasticlinear::calelasticmatrix
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix hdflag: 0=full, 1=deviatoric only (edge-smoothed,...
Definition: ElasticLinear.f90:17
m_elastoplastic
This module provide functions for elastoplastic calculation.
Definition: Elastoplastic.f90:6
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
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
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:763
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_elastoplastic::updateepstate
subroutine, public updateepstate(gauss)
Clear elatoplastic state.
Definition: Elastoplastic.f90:880
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
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:143
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:628
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::getyieldfunction
integer function getyieldfunction(mtype)
Get type of yield function.
Definition: material.f90:321
mmaterial::m_plconst5
integer(kind=kint), parameter m_plconst5
Definition: material.f90:102
mmaterial::m_plconst4
integer(kind=kint), parameter m_plconst4
Definition: material.f90:101
mmaterial::m_plconst1
integer(kind=kint), parameter m_plconst1
Definition: material.f90:98
mmaterial::gethardentype
integer function gethardentype(mtype)
Get type of hardening.
Definition: material.f90:333
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:84
mmaterial::m_plconst2
integer(kind=kint), parameter m_plconst2
Definition: material.f90:99
mmaterial::mc_yield
character(len=dict_key_length) mc_yield
Definition: material.f90:142
mmaterial::iskinematicharden
logical function iskinematicharden(mtype)
If it is a kinematic hardening material?
Definition: material.f90:345
mmaterial::mc_isoelastic
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:140
mmaterial::m_plconst3
integer(kind=kint), parameter m_plconst3
Definition: material.f90:100
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
muyield
This subroutine read in used-defined material properties tangent.
Definition: uyield.f90:7
muyield::uelastoplasticmatrix
subroutine, public uelastoplasticmatrix(matl, stress, istat, fstat, plstrain, D, temp, hdflag)
This subroutine calculates elastoplastic constitutive relation.
Definition: uyield.f90:32
muyield::ubackwardeuler
subroutine, public ubackwardeuler(matl, stress, plstrain, istat, fstat, temp, hdflag)
This subroutine does backward-Euler return calculation.
Definition: uyield.f90:45
mmaterial::tmaterial
Structure to manage all material related data.
Definition: material.f90:167
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:16