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/), &
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
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)
46 if( secttype /=
d3 ) stop
"Elastoplastic calculation support only Solid element currently"
50 back(1:6) = extval(2:7)
53 if(
present( temperature ) )
then
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) )
67 if( istat == 0 )
return
71 dj2(4:6) = 2.d0*devia(4:6)
72 dj2 = dj2/( 2.d0*dsqrt(j2) )
73 if(
present(temperature) )
then
81 a(1:6) = devia(1:6)/dsqrt(2.d0*j2)
83 dlambda = extval(1)-plstrain
84 c3 = dsqrt(3.d0*j2)+3.d0*g*dlambda
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)
91 d(i,j) = de(i,j) - c1*id(i,j) + c2*a(i)*a(j)
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
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) ) )
114 c3 = dsqrt(3.d0)*sin(sita)+cos(sita)*sin(fai)/(2.d0*j2*cos(3.d0*sita))
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
132 a(:) = fai*dj1(:) + dj2(:)
136 dum = harden + khard+ dot_product( da, a )
139 d(i,j) = de(i,j) - da(i)*da(j)/dum
148 real(kind=
kreal),
intent(in) :: stress(6)
149 real(kind=
kreal),
intent(in) :: extval(:)
153 real(kind=
kreal) :: eqvs, sita, fai, j1,j2,j3, devia(6)
154 real(kind=
kreal) :: back(6)
156 if( kinematic ) back(1:6) = extval(2:7)
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) )
168 eqvs = dsqrt( 3.d0*j2 )
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) &
193 real(kind=
kreal),
intent(in) :: strain(6)
202 real(kind=
kreal),
intent(in) :: pstrain
203 real(kind=
kreal),
intent(in),
optional :: temp
207 real(kind=
kreal) :: s0, s1,s2, ef, ina(2)
215 if(
present(temp) )
then
216 ina(1) = temp; ina(2)=pstrain
233 if(
present(temp) )
then
249 real(kind=
kreal),
intent(in) :: pstrain
264 real(kind=
kreal),
intent(in) :: pstrain
279 real(kind=
kreal),
intent(in) :: pstrain
280 real(kind=
kreal),
intent(in),
optional :: temp
283 real(kind=
kreal) :: s0, s1,s2, ina(2), outa(1)
292 if(
present(temp) )
then
293 ina(1) = temp; ina(2)=pstrain
294 call fetch_tabledata(
mc_yield, matl%dict, outa, ierr, ina)
297 call fetch_tabledata(
mc_yield, matl%dict, outa, ierr, ina(1:1))
299 if( ierr ) stop
"Fail to get yield stress!"
310 if( pstrain<=s0 )
then
323 real(kind=
kreal),
intent(in) :: stress(6)
324 real(kind=
kreal),
intent(in) :: extval(:)
325 real(kind=
kreal),
intent(in),
optional :: temp
329 real(kind=
kreal) :: eqvs, sita, eta, fai, j1,j2,j3, f, devia(6)
330 real(kind=
kreal) :: pstrain, back(6)
335 if( kinematic ) back(1:6) = extval(2:7)
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
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
354 f = dsqrt( 3.d0*j2 ) - eqvs
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)
370 f = dsqrt(j2) + eta*j1 - eqvs*matl%variables(
m_plconst4)
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
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)
408 fstat_bak(1) = plstrain
409 if(
present(temp) )
then
414 if( dabs(f)<tol )
then
417 elseif( f<0.d0 )
then
422 kh = 0.d0; kk=0.d0; betan=0.d0; back(:)=0.d0
426 back(1:6) = fstat(8:13)
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
436 if(
present(temp) )
then
438 call fetch_tabledata(
mc_isoelastic, matl%dict, ee, ierr, ina)
443 stop
" fail to fetch young's modulus in elastoplastic calculation"
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) )
455 if(
present(temp) )
then
464 dlambda = dlambda+f/dd
465 if( dlambda<0.d0 )
then
469 if(
present(temp) )
then
477 f = yd-3.d0*g*dlambda-dum -(kk-betan)
478 if( dabs(f)<tol*tol )
exit
480 pstrain = pstrain+dlambda
483 fstat(2:7) = back(:)+(kk-betan)*devia(:)/yd
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
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
505 if( dabs(stress(mm))<1.d-10 ) stress(mm)=0.d0
507 call eigen3( stress, prnstre, prnprj )
509 maxp = maxloc( prnstre )
510 minp = minloc( prnstre )
512 if( maxp(1)==1 .or. minp(1)==1 ) mm =2
513 if( maxp(1)==2 .or. minp(1)==2 ) mm =3
515 if(
present(temp) )
then
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"
528 dum = pstrain + 2.d0*dlambda*cos(fai)
529 if(
present(temp) )
then
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
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
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) )
557 elseif(ytype==2)
then
561 if(
present(temp) )
then
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"
573 if(
present(temp) )
then
578 f = yd-g*dlambda+fai*(j1-k*fai*dlambda)- dum*f
579 if( dabs(f)<tol*tol )
exit
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)
595 gauss%plstrain= gauss%fstatus(1)
597 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)
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.
subroutine eigen3(tensor, eigval, princ)
Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using Jacobi iteration adapted from numeri...
This module summarizes all infomation of material properties.
integer function getyieldfunction(mtype)
Get type of yield function.
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.
Stucture to management all material relates data.
All data should be recorded in every quadrature points.
subroutine uelastoplasticmatrix(matl, stress, istat, fstat, d)
This subroutine read in used-defined material properties tangent This subroutine calculates elastopla...
subroutine ubackwardeuler(matl, stress, istat, fstat)
This subroutine does backward-Euler return calculation.