creep.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module mcreep
  use hecmw_util
  use mmaterial
  use m_elasticlinear
 
  implicit none
 
contains
 
  subroutine iso_creep(matl, sectType, stress, strain, extval,plstrain,            &
      dtime,ttime,stiffness, temp, hdflag)
    type( tmaterial ), intent(in)    :: matl
    integer, intent(in)              :: sectType
    real(kind=kreal), intent(in)     :: stress(6) 
    real(kind=kreal), intent(in)     :: strain(6) 
    real(kind=kreal), intent(in)     :: extval(:) 
    real(kind=kreal), intent(in)     :: plstrain  
    real(kind=kreal), intent(in)     :: ttime     
    real(kind=kreal), intent(in)     :: dtime     
    real(kind=kreal), intent(out)    :: stiffness(6,6) 
    real(kind=kreal), intent(in)     :: temp
    integer(kind=kint), intent(in), optional :: hdflag
 
    integer :: i, j, hdflag_in
    logical :: ierr
    real(kind=kreal) :: ina(1), outa(3)
    real(kind=kreal) :: xxn, aa
 
    real(kind=kreal) :: c3,e,un,g,stri(6),p,dstri,c4,c5,f,df, eqvs
 
    hdflag_in = 0
    if( present(hdflag) ) hdflag_in = hdflag
 
    aa = 0.0d0
    xxn = 0.0d0
    !
    !     elastic
    !
    ina(1) = temp
    call calelasticmatrix( matl, secttype, stiffness, temp, hdflag=hdflag_in )
 
    if( dtime==0.d0 .or. all(stress==0.d0) ) return
    if( hdflag_in == 2 ) return
    !
    !     elastic constants
    !
    ina(1) = temp
    call fetch_tabledata( mc_isoelastic, matl%dict, outa(1:2), ierr, ina )
 
    if( ierr ) then
      stop "error in isotropic elasticity definition"
    else
      e=outa(1)
      un=outa(2)
    endif
 
    !      Norton
    if( matl%mtype==norton ) then         ! those with no yield surface
      ina(1) = temp
      call fetch_tabledata( mc_norton, matl%dict, outa, ierr, ina )
      xxn=outa(2)
      aa=outa(1)*((ttime+dtime)**(outa(3)+1.d0)-ttime**(outa(3)+1.d0))/(outa(3)+1.d0)
    endif
 
    g=0.5d0*e/(1.d0+un)
    !
    !     creep
    !
    stri(:)=stress(:)
    p=-(stri(1)+stri(2)+stri(3))/3.d0
    do i=1,3
      stri(i)=stri(i)+p
    enddo
    !
    dstri=dsqrt(1.5d0*(stri(1)*stri(1)+stri(2)*stri(2)+stri(3)*stri(3)+    &
      2.d0*(stri(4)*stri(4)+stri(5)*stri(5)+stri(6)*stri(6))) )
    !
    !     unit trial vector
    !
    stri(:)=stri(:)/dstri
 
    eqvs = dstri
    if( eqvs<1.d-10 ) eqvs=1.d-10
    f=aa*eqvs**xxn
    df=xxn*f/eqvs
    !
    !     stiffness matrix
    !
    c3=6.d0*g*g
    c4=c3*plstrain/(dstri+3.d0*g*plstrain)
    c3=c4-c3*df/(3.d0*g*df+1.d0)
    c5=c4/3.d0
 
    do i=1,6
      do j=1,6
        stiffness(i,j) = stiffness(i,j) +c3*stri(i)*stri(j)
      enddo
    enddo
    do i=1,3
      stiffness(i,i) = stiffness(i,i) - c4
      do j=1,3
        stiffness(i,j) = stiffness(i,j) +c5
      enddo
    enddo
    do i=4,6
      stiffness(i,i) = stiffness(i,i) - c4/2.d0
    enddo
 
  end subroutine
 
  subroutine update_iso_creep(matl, sectType, strain, stress, extval, plstrain,         &
      dtime, ttime, temp, hdflag)
    type( tmaterial ), intent(in)    :: matl
    integer, intent(in)              :: secttype
    real(kind=kreal), intent(in)     :: strain(6) 
    real(kind=kreal), intent(inout)  :: stress(6) 
    real(kind=kreal), intent(inout)  :: extval(:) 
    real(kind=kreal), intent(out)    :: plstrain  
    real(kind=kreal), intent(in)     :: ttime     
    real(kind=kreal), intent(in)     :: dtime     
    real(kind=kreal), intent(in)     :: temp
    integer(kind=kint), intent(in), optional :: hdflag
 
    integer :: i, hdflag_in
    logical :: ierr
    real(kind=kreal) :: ina(1), outa(3)
    real(kind=kreal) :: xxn, aa
 
    real(kind=kreal) :: e,un,g,dg,ddg,stri(6),p,dstri,f,df, eqvs
 
    hdflag_in = 0
    if( present(hdflag) ) hdflag_in = hdflag
 
    aa = 0.0d0
    xxn = 0.0d0
 
    if( dtime==0.d0 ) return
    if( hdflag_in == 2 ) return
    !
    !     elastic constants
    !
    ina(1) = temp
    call fetch_tabledata( mc_isoelastic, matl%dict, outa(1:2), ierr, ina )
    if( ierr ) then
      stop "error in isotropic elasticity definition"
    else
      e=outa(1)
      un=outa(2)
    endif
 
    !      Norton
    if( matl%mtype==norton ) then         ! those with no yield surface
      ina(1) = temp
      call fetch_tabledata( mc_norton, matl%dict, outa, ierr, ina )
      if( ierr ) then
        stop "error in isotropic elasticity definition"
      else
        xxn=outa(2)
        aa=outa(1)*((ttime+dtime)**(outa(3)+1.d0)-ttime**(outa(3)+1.d0))/(outa(3)+1.d0)
      endif
    endif
 
    g=0.5d0*e/(1.d0+un)
 
    !
    !     creep
    !
    stri(:)=stress(:)
    p=-(stri(1)+stri(2)+stri(3))/3.d0
    do i=1,3
      stri(i)=stri(i)+p
    enddo
    !
    dstri=dsqrt(1.5d0*(stri(1)*stri(1)+stri(2)*stri(2)+stri(3)*stri(3)+    &
      2.d0*(stri(4)*stri(4)+stri(5)*stri(5)+stri(6)*stri(6))) )
    !
    !     determination of the consistency parameter
    !
    dg=0.d0
    do
      if( matl%mtype==norton ) then
        eqvs = dstri-3.d0*g*dg
        f=aa*eqvs**xxn
        df=xxn*f/eqvs
        ddg = (f-dg)/(3.d0*g*df+1.d0)
        dg = dg+ddg
        if((ddg<dg*1.d-6).or.(ddg<1.d-12)) exit
      endif
    enddo
 
    stri(:) = stri(:)-3.d0*g*dg*stri(:)/dstri
    stress(1:3) = stri(1:3)-p
    stress(4:6) = stri(4:6)
 
    !
    !     state variables
    !
    plstrain= dg
    extval(1)=eqvs
 
  end subroutine
 
  subroutine updateviscostate( gauss )
    use mmechgauss
    type(tgaussstatus), intent(inout) :: gauss  ! status of curr gauss point
 
    gauss%fstatus(2) = gauss%fstatus(2)+gauss%plstrain
  end subroutine
 
end module