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