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/), &
27 real(kind=
kreal),
parameter :: i2(6) = (/ 1.d0, 1.d0, 1.d0, 0.d0, 0.d0, 0.d0 /)
29 integer,
parameter :: VM_ELASTIC = 0
30 integer,
parameter :: VM_PLASTIC = 1
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
38 integer,
parameter :: DP_ELASTIC = 0
39 integer,
parameter :: DP_PLASTIC_SURF = 1
40 integer,
parameter :: DP_PLASTIC_APEX = 2
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
56 integer :: ytype,hdflag_in
59 if(
present(hdflag) ) hdflag_in = hdflag
64 call calelastoplasticmatrix_vm( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
66 call calelastoplasticmatrix_mc( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
70 call uelastoplasticmatrix( matl%variables, stress, istat, extval, plstrain, d, temperature, hdflag_in )
75 subroutine calelastoplasticmatrix_vm( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
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
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)
92 if( secttype /=
d3 ) stop
"Elastoplastic calculation support only Solid element currently"
95 if( istat == vm_elastic )
return
96 if( hdflag == 2 )
return
98 harden = calhardencoeff( matl, extval(1), temperature )
103 back(1:6) = extval(2:7)
104 khard = calkinematicharden( matl, extval(1) )
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) )
114 a(1:6) = devia(1:6)/sqrt(2.d0*j2)
116 dlambda = extval(1)-plstrain
117 c3 = sqrt(3.d0*j2)+3.d0*g*dlambda
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)
124 d(i,j) = d(i,j) - c1*id(i,j) + c2*a(i)*a(j)
127 end subroutine calelastoplasticmatrix_vm
130 subroutine calelastoplasticmatrix_mc( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
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
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)
151 if( secttype /=
d3 ) stop
"Elastoplastic calculation support only Solid element currently"
154 if( istat == mc_elastic )
return
155 if( hdflag == 2 )
return
157 harden = calhardencoeff( matl, extval(1), temperature )
159 k = d(1,1)-(4.d0/3.d0)*g
164 youngs = 9.d0*k*g/(3.d0*k+g)
165 poisson = (3.d0*k-r2g)/(6.d0*k+r2g)
172 sphsps = sinphi*sinpsi
173 r2cosphi = 2.d0*cosphi
174 r4cos2phi = r2cosphi*r2cosphi
176 call eigen3( stress, prnstre, prnprj )
177 m1 = maxloc( prnstre, 1 )
178 m3 = minloc( prnstre, 1 )
180 m1 = 1; m2 = 2; m3 = 3
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)))
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
214 dum = r4cos2phi*harden
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
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
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
260 end subroutine calelastoplasticmatrix_mc
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
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
279 if( secttype /=
d3 ) stop
"Elastoplastic calculation support only Solid element currently"
282 if( istat == dp_elastic )
return
283 if( hdflag == 2 )
return
285 harden = calhardencoeff( matl, extval(1), temperature )
288 k = d(1,1)-(4.d0/3.d0)*g
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) )
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
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)) &
320 c1 = k*(1.d0 - k/(k + alpha*beta*harden))
323 d(i,j) = c1*i2(i)*i2(j)
330 real(kind=
kreal)
function calhardencoeff( matl, pstrain, temp )
332 real(kind=
kreal),
intent(in) :: pstrain
333 real(kind=
kreal),
intent(in) :: temp
337 real(kind=
kreal) :: s0, s1,s2, ef, ina(2)
339 calhardencoeff = -1.d0
345 ina(1) = temp; ina(2)=pstrain
346 call fetch_tablegrad(
mc_yield, ina, matl%dict, calhardencoeff, ierr )
351 calhardencoeff = s1*s2*( s0+pstrain )**(s2-1)
356 ef = calcurryield( matl, pstrain, temp )
357 calhardencoeff = s1*(ef/s1)**(1.d0-s2) /(s0*s2)
359 calhardencoeff = 0.d0
366 real(kind=
kreal)
function calkinematicharden( matl, pstrain )
368 real(kind=
kreal),
intent(in) :: pstrain
374 calkinematicharden = matl%variables(
m_plconst3)
376 calkinematicharden = 0.d0
381 real(kind=
kreal)
function calcurrkinematic( matl, pstrain )
383 real(kind=
kreal),
intent(in) :: pstrain
389 calcurrkinematic = matl%variables(
m_plconst3)*pstrain
391 calcurrkinematic = 0.d0
396 real(kind=
kreal)
function calcurryield( matl, pstrain, temp )
398 real(kind=
kreal),
intent(in) :: pstrain
399 real(kind=
kreal),
intent(in) :: temp
402 real(kind=
kreal) :: s0, s1,s2, ina(2), outa(1)
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)
419 calcurryield = s1*( s0+pstrain )**s2
424 if( pstrain<=s0 )
then
427 calcurryield = s1*( pstrain/s0 )**(1.d0/s2)
435 subroutine backwardeuler( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag )
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
445 integer :: ytype, hdflag_in
448 if(
present(hdflag) ) hdflag_in = hdflag
453 call backwardeuler_vm( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag_in )
455 call backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag_in )
457 call backwardeuler_dp( matl, stress, plstrain, istat, fstat, temp, hdflag_in )
459 call ubackwardeuler( matl%variables, stress, plstrain, istat, fstat, temp, hdflag_in )
464 subroutine backwardeuler_vm( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag )
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
474 real(kind=
kreal),
parameter :: tol =1.d-8
475 integer,
parameter :: maxiter = 10
476 real(kind=
kreal) :: dlambda, f
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)
485 if( kinematic ) back(1:6) = fstat(8:13)
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) )
494 eqvs = dsqrt( 3.d0*j2 )
495 yd = calcurryield( matl, plstrain, temp )
498 if( abs(f/yd)<tol )
then
501 elseif( f<0.d0 )
then
505 if( hdflag == 2 )
return
508 kh = 0.d0; kk=0.d0; betan=0.d0; back(:)=0.d0
509 if( kinematic ) betan = calcurrkinematic( matl, plstrain )
512 call fetch_tabledata(
mc_isoelastic, matl%dict, ee, ierr, ina)
514 stop
" fail to fetch young's modulus in elastoplastic calculation"
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) )
527 h= calhardencoeff( matl, pstrain, temp )
529 kh = calkinematicharden( matl, pstrain )
532 dlambda = dlambda+f/dd
533 if( dlambda<0.d0 )
then
536 istat=vm_elastic;
exit
538 pstrain = plstrain+dlambda
539 yd = calcurryield( matl, pstrain, temp )
541 kk = calcurrkinematic( matl, pstrain )
543 f = eqvs-3.d0*g*dlambda-yd -(kk-betan)
544 if( abs(f/yd)<tol )
exit
550 kk = calcurrkinematic( matl, pstrain )
551 fstat(2:7) = back(:)+(kk-betan)*devia(:)/eqvs
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(:)
560 h= calhardencoeff( matl, pstrain, temp )
561 yd = calcurryield( matl, plstrain, temp )
562 plpotential = -0.5d0*(eqvs-yd)*(eqvs-yd)/(h+3.d0*g)
564 end subroutine backwardeuler_vm
567 subroutine backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag )
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
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)
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
595 r2cosphi = 2.d0*cosphi
597 call eigen3( stress, prnstre, prnprj )
599 m1 = maxloc( prnstre, 1 )
600 m3 = minloc( prnstre, 1 )
602 m1 = 1; m2 = 2; m3 = 3
607 eqvs = prnstre(m1)-prnstre(m3) + (prnstre(m1)+prnstre(m3))*sinphi
608 cohe = calcurryield( matl, plstrain, temp )
609 f = eqvs - r2cosphi*cohe
611 if( abs(f/cohe)<tol )
then
612 istat = mc_plastic_surf
614 elseif( f<0.d0 )
then
618 if( hdflag == 2 )
return
620 istat = mc_plastic_surf
623 call fetch_tabledata(
mc_isoelastic, matl%dict, ee, ierr, ina)
625 stop
" fail to fetch young's modulus in elastoplastic calculation"
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) )
638 sphsps = sinphi*sinpsi
639 r4cos2phi = r2cosphi*r2cosphi
640 c1 = 4.d0*(g*(1.d0+sphsps/3.d0)+k*sphsps)
642 h= calhardencoeff( matl, pstrain, temp )
644 dlambda = dlambda+f/dd
645 if( r2cosphi*dlambda<0.d0 )
then
646 if( cosphi==0.d0 ) stop
"Math error in return mapping"
649 istat = mc_elastic;
exit
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
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
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
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
680 cohe = calcurryield( matl, plstrain, temp )
681 f = eqvs - r2cosphi*cohe
682 fb = eqvsb - r2cosphi*cohe
685 h= calhardencoeff( matl, pstrain, temp )
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
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
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)
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
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'
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
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 )
732 resid = cotphi*cohe-p
733 if( abs(resid/cohe)<tol )
exit
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) )
755 end subroutine backwardeuler_mc
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
767 real(kind=
kreal),
parameter :: tol =1.d-8
768 integer,
parameter :: maxiter = 10
769 real(kind=
kreal) :: dlambda, f
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
774 real(kind=
kreal) :: alpha, beta, depv, factor, resid
780 j1 = (stress(1)+stress(2)+stress(3))
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) )
788 cohe = calcurryield( matl, plstrain, temp )
789 f = eqvst + eta*pt - xi*cohe
791 if( abs(f/cohe)<tol )
then
792 istat = dp_plastic_surf
794 elseif( f<0.d0 )
then
798 if( hdflag == 2 )
return
800 istat = dp_plastic_surf
803 call fetch_tabledata(
mc_isoelastic, matl%dict, ee, ierr, ina)
805 stop
" fail to fetch young's modulus in elastoplastic calculation"
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) )
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"
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
837 if( eqvs>=0.d0 )
then
838 factor = 1.d0-g*dlambda/eqvst
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'
847 cohe = calcurryield( matl, pstrain, temp )
848 resid = beta*cohe - pt
850 h= calhardencoeff( matl, pstrain, temp )
852 depv = depv - resid/dd
853 pstrain = plstrain+alpha*depv
854 cohe = calcurryield( matl, pstrain, temp )
856 resid = beta*cohe - p
857 if( abs(resid/cohe)<tol )
then
867 devia(:) = factor*devia(:)
868 stress(1:3) = devia(1:3)+p
869 stress(4:6) = devia(4:6)
878 gauss%plstrain= gauss%fstatus(1)
880 gauss%fstatus(8:13) =gauss%fstatus(2:7)
integer(kind=4), parameter kreal
This module provides functions for elastic material.
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix.
This module provide functions for elastoplastic calculation.
subroutine calelastoplasticmatrix_dp(matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag)
This subroutine calculates elastoplastic constitutive relation.
subroutine, public backwardeuler(matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag)
This subroutine does backward-Euler return calculation.
subroutine backwardeuler_dp(matl, stress, plstrain, istat, fstat, temp, hdflag)
This subroutine does backward-Euler return calculation for Drucker-Prager.
subroutine, public calelastoplasticmatrix(matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag)
This subroutine calculates elastoplastic constitutive relation.
subroutine, public updateepstate(gauss)
Clear elatoplastic state.
This module provides aux functions.
subroutine eigen3(tensor, eigval, princ)
Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using Jacobi iteration adapted from numeri...
subroutine deriv_general_iso_tensor_func_3d(dpydpx, dydx, eigv, px, py)
Compute derivative of a general isotropic tensor function of one tensor.
This module summarizes all information of material properties.
integer function getyieldfunction(mtype)
Get type of yield function.
integer(kind=kint), parameter m_plconst5
integer(kind=kint), parameter m_plconst4
integer(kind=kint), parameter m_plconst1
integer function gethardentype(mtype)
Get type of hardening.
integer(kind=kint), parameter d3
integer(kind=kint), parameter m_plconst2
character(len=dict_key_length) mc_yield
logical function iskinematicharden(mtype)
If it is a kinematic hardening material?
character(len=dict_key_length) mc_isoelastic
integer(kind=kint), parameter m_plconst3
This modules defines a structure to record history dependent parameter in static analysis.
This subroutine read in used-defined material properties tangent.
subroutine, public uelastoplasticmatrix(matl, stress, istat, fstat, plstrain, D, temp, hdflag)
This subroutine calculates elastoplastic constitutive relation.
subroutine, public ubackwardeuler(matl, stress, plstrain, istat, fstat, temp, hdflag)
This subroutine does backward-Euler return calculation.
Structure to manage all material related data.
All data should be recorded in every quadrature points.