Viscoelastic.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 mviscoelastic
  use hecmw_util
  use mmaterial
  implicit none
 
  private :: hvisc, trs, trsinc
 
contains
 
  real(kind=kreal) function hvisc(x,expx)
    real(kind=kreal), intent(in) :: x       
    real(kind=kreal), intent(in) :: expx    
 
    if(x < 1.d-04) then
 
      hvisc = 1.d0 - 0.5d0*x*(1.d0 - x/3.d0*(1.d0     &
        - 0.25d0*x*(1.d0 - 0.2d0*x)))
 
    else
 
      hvisc = (1.d0 - expx)/x
 
    endif
 
  end function
 
  real(kind=kreal) function trsinc(tn1,tn,mtype, mvar)
    real(kind=kreal), intent(in) :: tn1     
    real(kind=kreal), intent(in) :: tn      
    integer, intent(in)          :: mtype
    real(kind=kreal), intent(in) :: mvar(:) 
 
    real(kind=kreal)  ::  hsn1, hsn, asn1, asn
 
    if(mtype==viscoelastic+2) then   ! Arrhenius
 
      hsn = -mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
      hsn1 = -mvar(2)*( 1.d0/(tn1-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
 
    else   ! WLF
 
      if( mvar(3)+tn1-mvar(1)<=0.d0 ) then
        trsinc= -1.d0; return
      endif
      if( mvar(3)+tn-mvar(1)<=0.d0 ) then
        trsinc= -1.d0; return
      endif
      hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
      hsn1 = mvar(2)*( tn1-mvar(1) )/ (mvar(3)+tn1-mvar(1) )*dlog(10.d0)
 
    endif
 
    if( dabs(hsn1-hsn)<1.d-5 ) then
      trsinc = dexp((hsn1+hsn)*0.5d0)
    else
      asn1 = dexp(hsn1)
      asn = dexp(hsn)
      trsinc = (asn1-asn)/(hsn1-hsn)
    endif
 
  end function
 
  real(kind=kreal) function trs(tn,mtype, mvar)
    real(kind=kreal), intent(in) :: tn      
    integer, intent(in)          :: mtype
    real(kind=kreal), intent(in) :: mvar(:) 
 
    real(kind=kreal)  ::  hsn
 
    if(mtype==viscoelastic+2) then   ! Arrhenius
      hsn = mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
    else   ! WLF
      hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
    endif
 
    trs = dexp(hsn)
 
  end function
 
  !-------------------------------------------------------------------------------
  !-------------------------------------------------------------------------------
  subroutine calviscoelasticmatrix(matl, sectType, dt, D, temp)
    use m_elasticlinear
    type( tmaterial ), intent(in) :: matl
    integer, intent(in)           :: secttype
    real(kind=kreal), intent(in)  :: dt        
    real(kind=kreal), intent(out) :: d(:,:)    
    real(kind=kreal), intent(in)  :: temp
 
    integer   i,j, n
    real(kind=kreal) :: g,gg,k,kg, gfac,exp_n,mu_0,mu_n,dq_n,dtau
    real(kind=kreal) :: ina(1), outa(2), ee, pp, ddt
 
    type(ttable), pointer  :: dicval
    logical :: ierr
 
    d(:,:)=0.d0
 
    ddt = dt
 
    ina(1) = temp
    call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
    if( ierr ) then
      ee = matl%variables(m_youngs)
      pp = matl%variables(m_poisson)
    else
      ee = outa(1)
      pp = outa(2)
    endif
    if( matl%mtype>viscoelastic ) ddt=trs(temp,matl%mtype, matl%variables)*dt
 
    if( ddt<=0.d0 ) then
      call calelasticmatrix( matl, d3, d, temp  )
      return
    endif
 
    !     Set elastic parameters for G (mu) and lambda
 
    g = ee/(2.d0*(1.d0 + pp))
    k = ee/(3.d0*(1.d0 - 2.d0*pp))
 
    !     Set properties for integrating the q_i terms
 
    mu_0 = 0.0d0
    gfac = 0.0d0
 
    call fetch_table( mc_viscoelastic, matl%dict, dicval, ierr )
    if( ierr ) stop "Viscoelastic properties not defined"
 
    do n = 1,dicval%tbrow
      mu_n  = dicval%tbval(1,n)
      dtau  = ddt/dicval%tbval(2,n)
      exp_n = dexp(-dtau)
 
      dq_n = mu_n * hvisc(dtau,exp_n)
      gfac = gfac + dq_n
      mu_0 = mu_0 + mu_n
 
    end do
 
    mu_0 = 1.d0 - mu_0
    gfac = gfac + mu_0
 
    !    Set tangent parameters
    gg = g*gfac
    kg = k - 0.6666666666666666d0*gg
    do j =1,3
      do i = 1,3
        d(i,j) = kg
      end do
      d(j,j) = d(j,j) + 2.d0*gg
    end do
 
    do i = 4,6
      d(i,i) = gg
    end do
 
 
  end subroutine
 
  !-------------------------------------------------------------------------------
  subroutine updateviscoelastic(matl, sectType, eps, sig, vsig, dt, temp, tempn)
    use m_elasticlinear
    type( tmaterial ), intent(in)   :: matl
    integer, intent(in)             :: sectType
    real(kind=kreal), intent(in)    :: eps(6)    
    real(kind=kreal), intent(out)   :: sig(6)    
    real(kind=kreal), intent(inout) :: vsig(:)   
    real(kind=kreal), intent(in)    :: dt        
    real(kind=kreal), intent(in)    :: temp
    real(kind=kreal), intent(in)    :: tempn
 
    integer   i,n
    real(kind=kreal) :: g,k,kth, exp_n,mu_0,mu_n,dq_n,dtau, theta
    real(kind=kreal) :: ina(1), outa(2), ee, pp, ddt
 
    real(kind=kreal) :: devstrain(6), en(6), d(6,6)
    type(ttable), pointer  :: dicval
    logical :: ierr
 
    ddt = dt
 
    ina(1) = temp
    call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
    if( ierr ) then
      ee = matl%variables(m_youngs)
      pp = matl%variables(m_poisson)
    else
      ee = outa(1)
      pp = outa(2)
    endif
    if( matl%mtype>viscoelastic ) ddt=trsinc(temp, tempn, matl%mtype, matl%variables)*dt
 
    if( ddt<=0.d0 ) then
      call calelasticmatrix( matl, secttype, d, temp )
      sig = matmul( d(:,:), eps(:) )
      vsig = 0.d0
      return
    endif
 
    !     Set elastic parameters for G (mu) and lambda
 
    g = ee/(2.d0*(1.d0 + pp))
    k = ee/(3.d0*(1.d0 - 2.d0*pp))
 
    !     Compute volumetric strain and deviatoric components
 
    theta = (eps(1) + eps(2) + eps(3))/3.d0
    do i = 1,3
      devstrain(i) = eps(i) - theta
    end do
    do i = 4,6
      devstrain(i) = eps(i)*0.5d0
    end do
 
    !     Set properties for integrating the q_i terms
 
    sig(:) = 0.0d0
    mu_0 = 0.0d0
 
    call fetch_table( mc_viscoelastic, matl%dict, dicval, ierr )
    if( ierr ) stop "Viscoelastic properties not defined"
 
    en(:) = vsig(12*dicval%tbrow+1:12*dicval%tbrow+6)
 
    do n = 1,dicval%tbrow
      mu_n  = dicval%tbval(1,n)
      dtau  = ddt/dicval%tbval(2,n)
      exp_n = dexp(-dtau)
 
      dq_n = mu_n * hvisc(dtau,exp_n)
      mu_0 = mu_0 + mu_n
 
      do i = 1,6
        vsig((n-1)*12+i+6) = exp_n*vsig((n-1)*12+i) + dq_n*(devstrain(i) - en(i))
        sig(i)  = sig(i) + vsig((n-1)*12+i+6)
      end do
    end do
 
    !     Finish stress update
 
    mu_0 = 1.d0 - mu_0
    do i = 1,6
      sig(i) = 2.d0*g*(mu_0*devstrain(i) + sig(i))
    end do
 
    !     Add elastic bulk term
 
    kth = k*theta*3.0d0
    do i = 1,3
      sig(i) = sig(i) + kth
    end do
 
  end subroutine
 
  subroutine updateviscoelasticstate( gauss )
    use mmechgauss
    type(tgaussstatus), intent(inout) :: gauss  ! status of curr gauss point
 
    integer :: i, n, nrow
    real(kind=kreal) :: thetan, vstrain(6)
    nrow = fetch_tablerow( mc_viscoelastic, gauss%pMaterial%dict )
    do n = 1,nrow
      do i = 1,6
        gauss%fstatus((n-1)*12+i) = gauss%fstatus((n-1)*12+i+6)
      end do
    end do
    vstrain = gauss%strain
    thetan = (vstrain(1) + vstrain(2) + vstrain(3))/3.d0
    do i = 1,3
      gauss%fstatus(nrow*12+i) = vstrain(i) - thetan
    end do
    do i = 4,6
      gauss%fstatus(nrow*12+i) = vstrain(i)*0.5d0
    end do
  end subroutine
 
end module