material.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 mmaterial
  use hecmw_util
  use m_table
  use table_dicts
  implicit none
 
  ! Following algorithm type
  integer(kind=kint), parameter :: infinitesimal = 0
  integer(kind=kint), parameter :: totallag  = 1
  integer(kind=kint), parameter :: updatelag = 2
 
  ! Following material types. All material type number consists with integer of six digits.
  !   First digit: Indicates physical type
  !                1: mechanical deformation analysis
  !                2: heat conduct analysis
  !                ......
  !   Second digit:
  !     Mechanical analysis
  !     1: Elastic
  !     2: Elastoplastic
  !     3: Hyperelastic
  !     4: Viscoelastic
  !     5: Viscoplastic
  !     6: Incomp newtonian
  !     7: Connector(Spring, Dashpot, Joint etc.)
  !     Heat conductiovity
  !     ......
  !   Third digit:
  !     For elastic or elastoplastic deformation, elastic
  !       0: isotropic                      ie. 110000
  !       1: with transever anisotropity        111000
  !     For hyperelastic deformation
  !       0: Neo-Hooke                          130000
  !       1: Mooney-Rivlin                      131000
  !       2: Arruda-Boyce                       132000
  !     For spring or dashpot, joint etc.
  !       0: spring dof                     ie. 170000
  !       1: spring axial                       171000
  !       2: dashpot dof                        172000
  !       3: dashpot axial                      173000
  !   Fourth digit
  !     For spring_d or dashpot_d, joint_d etc.
  !       k: number of dof param (1st digit)
  !   Fifth digit:
  !     For elastoplastic deformation, hardening law
  !       0: Linear hardening           i.e.  120000
  !       1: Multilinear hardening            120010
  !       2: Swift
  !       3: Ramberg-Osgood
  !       4: linear kinematic
  !       5: combined (linear kinematic + linear isotropic)
  !     For spring_d or dashpot_d, joint_d etc.
  !       k: number of dof param (2nd digit)
  !   Six digit:
  !     For visco-elastoplastic deformation, visco law
  !       0: Norton                     i.e.  150000
  !       1: Striab                           150001
  integer(kind=kint), parameter :: usermaterial          = 100000
 
  integer(kind=kint), parameter :: elastic               = 110000
  integer(kind=kint), parameter :: mn_orthoelastic       = 111000
  integer(kind=kint), parameter :: userelastic           = 112000
 
  integer(kind=kint), parameter :: eplastic              = 120000
 
  integer(kind=kint), parameter :: neohooke              = 130000
  integer(kind=kint), parameter :: mooneyrivlin          = 131000
  integer(kind=kint), parameter :: arrudaboyce           = 132000
  integer(kind=kint), parameter :: userhyperelastic      = 133000
  integer(kind=kint), parameter :: mooneyrivlin_aniso    = 134000
 
  integer(kind=kint), parameter :: viscoelastic          = 140000
  integer(kind=kint), parameter :: norton                = 150000
 
  integer(kind=kint), parameter :: incomp_newtonian      = 160000
  integer(kind=kint), parameter :: connector             = 170000
 
  ! Following section type
  integer(kind=kint), parameter :: d3            = -1
  integer(kind=kint), parameter :: planestress   = 1
  integer(kind=kint), parameter :: planestrain   = 0
  integer(kind=kint), parameter :: axissymetric  = 2
  integer(kind=kint), parameter :: shell         = 3
 
  ! Material constants are saved in an array of size 100 and their physical meaning
  ! correspond to their position in the array
  integer(kind=kint), parameter :: m_youngs  = 1
  integer(kind=kint), parameter :: m_poisson = 2
  integer(kind=kint), parameter :: m_density = 3
 
  ! following plastic constitutive parameter
  integer(kind=kint), parameter :: m_plconst1 = 5
  integer(kind=kint), parameter :: m_plconst2 = 6
  integer(kind=kint), parameter :: m_plconst3 = 7
  integer(kind=kint), parameter :: m_plconst4 = 8
  integer(kind=kint), parameter :: m_plconst5 = 9
  integer(kind=kint), parameter :: m_kinehard = 10
 
  integer(kind=kint), parameter :: m_exapnsion = 20
 
  integer(kind=kint), parameter :: m_alpha_over_mu = 21
 
  integer(kind=kint), parameter :: m_beam_radius = 22
  integer(kind=kint), parameter :: m_beam_angle1 = 23
  integer(kind=kint), parameter :: m_beam_angle2 = 24
  integer(kind=kint), parameter :: m_beam_angle3 = 25
  integer(kind=kint), parameter :: m_beam_angle4 = 26
  integer(kind=kint), parameter :: m_beam_angle5 = 27
  integer(kind=kint), parameter :: m_beam_angle6 = 28
 
  integer(kind=kint), parameter :: m_viscocity = 29
 
  ! additional plastic constitutive parameter
  integer(kind=kint), parameter :: m_plconst6 = 30
  integer(kind=kint), parameter :: m_plconst7 = 31
  integer(kind=kint), parameter :: m_plconst8 = 32
  integer(kind=kint), parameter :: m_plconst9 = 33
  integer(kind=kint), parameter :: m_plconst10 = 34
 
  integer(kind=kint), parameter :: m_damping_rm = 35
  integer(kind=kint), parameter :: m_damping_rk = 36
 
  integer(kind=kint), parameter :: m_spring_dof    = 0
  integer(kind=kint), parameter :: m_spring_axial  = 1
  integer(kind=kint), parameter :: m_dashpot_dof   = 2
  integer(kind=kint), parameter :: m_dashpot_axial = 3
 
  integer(kind=kint), parameter :: m_spring_d_ndoffset  = 0
  integer(kind=kint), parameter :: m_spring_a_ndoffset  = 72
  integer(kind=kint), parameter :: m_dashpot_d_ndoffset = 73
  integer(kind=kint), parameter :: m_dashpot_a_ndoffset = 145
 
  ! Dictionary constants
  character(len=DICT_KEY_LENGTH) :: mc_isoelastic= 'ISOELASTIC'      ! youngs modulus, poisson's ratio
  character(len=DICT_KEY_LENGTH) :: mc_orthoelastic= 'ORTHOELASTIC'  ! ortho elastic modulus
  character(len=DICT_KEY_LENGTH) :: mc_yield = 'YIELD'               ! plastic strain, yield stress
  character(len=DICT_KEY_LENGTH) :: mc_themoexp = 'THEMOEXP'         ! thermo expansion coefficient
  character(len=DICT_KEY_LENGTH) :: mc_orthoexp = 'ORTHOEXP'         ! thermo expansion coefficient
  character(len=DICT_KEY_LENGTH) :: mc_viscoelastic = 'VISCOELASTIC' ! Prony coeff only curr.
  character(len=DICT_KEY_LENGTH) :: mc_norton = 'NORTON'             ! NOrton's creep law
  character(len=DICT_KEY_LENGTH) :: mc_incomp_newtonian = 'INCOMP_FLUID' ! viscocity
  character(len=DICT_KEY_LENGTH) :: mc_spring= 'SPRING'              ! spring
  character(len=DICT_KEY_LENGTH) :: mc_dashpot= 'DASHPOT'            ! dashpot
 
  type tshellmat
    integer(kind=kint)         :: ortho
    real(kind=kreal)           :: ee
    real(kind=kreal)           :: pp
    real(kind=kreal)           :: ee2
    real(kind=kreal)           :: g12
    real(kind=kreal)           :: g23
    real(kind=kreal)           :: g31
    real(kind=kreal)           :: angle
    real(kind=kreal)           :: rho
    real(kind=kreal)           :: alpha
    real(kind=kreal)           :: alpha_over_mu
    real(kind=kreal)           :: weight
  end type tshellmat
 
  type tmaterial
    integer(kind=kint)         :: nlgeom_flag
    integer(kind=kint)         :: mtype
    integer(kind=kint)         :: nfstatus
    character(len=30)          :: name
    real(kind=kreal)           :: variables(200)    
    integer(kind=kint)         :: variables_i(200)
    type(tshellmat), pointer   :: shell_var(:) => null() 
    integer(kind=kint)         :: totallyr
    integer(kind=kint)         :: cdsys_id
    integer(kind=kint)         :: n_table
    real(kind=kreal), pointer  :: table(:)=>null()  
    type(dict_struct), pointer :: dict
    logical :: is_elem_rayleigh_damping
  end type tmaterial
 
  type(tmaterial), allocatable :: materials(:)
 
contains
 
  subroutine initmaterial( material )
    type( tmaterial ), intent(inout) :: material
    material%mtype = -1                  ! not defined yet
    material%nfstatus = 0                ! Default: no status
    material%nlgeom_flag = infinitesimal ! Default: INFINITESIMAL ANALYSIS
    material%variables = 0.d0           ! not defined yet
    material%variables_i =  0            ! not defined yet
    material%totallyr = 0               ! not defined yet
    material%is_elem_Rayleigh_damping = .false. ! not defined yet
 
    call dict_create( material%dict, 'INIT', dict_null )
  end subroutine
 
  subroutine finalizematerial( material )
    type( tmaterial ), intent(inout) :: material
    if( associated(material%table) ) deallocate( material%table )
    if( associated(material%dict) ) call dict_destroy( material%dict )
  end subroutine finalizematerial
 
  subroutine initializematls( nm )
    integer, intent(in) :: nm
    integer :: i
    if( allocated(materials) ) deallocate( materials )
    allocate( materials( nm ) )
    do i=1,nm
      call initmaterial( materials(i) )
    enddo
  end subroutine
 
  subroutine finalizematls()
    integer :: i
    if( allocated( materials ) ) then
      do i=1,size(materials)
        call finalizematerial( materials(i) )
      enddo
      deallocate( materials )
    endif
  end subroutine
 
  integer function fetchdigit( npos, cnum )
    integer, intent(in) :: npos
    integer, intent(in) :: cnum
    integer :: i, idum,cdum,dd
    fetchdigit = -1
    cdum = cnum
    if( npos<=0 .or. npos>6) return
    if( cnum<100000 .or. cnum>999999 ) return
    dd = 100000
    do i=1,npos-1
      idum = cdum/dd
      cdum = cdum-idum*dd
      dd = dd/10
    enddo
    fetchdigit = cdum/10**(6-npos)
  end function
 
  subroutine setdigit( npos, ival, mtype )
    integer, intent(in)              :: npos
    integer, intent(in)              :: ival
    integer, intent(inout)           :: mtype
    integer :: i, idum,cdum, cdum1, dd
    cdum = mtype
    if( npos<=0 .or. npos>6 ) return
    if( ival<0 .or. ival>9 ) return
    dd =100000
    cdum1 = 0
    do i=1,npos-1
      idum = cdum/dd
      cdum1 = cdum1+ idum*dd
      cdum = cdum-idum*dd
      dd=dd/10
    enddo
    cdum1 = cdum1 + ival*dd
    idum = cdum/dd
    cdum = cdum-idum*dd
    dd=dd/10
    do i=npos+1,6
      idum = cdum/dd
      cdum1 = cdum1+ idum*dd
      cdum = cdum-idum*dd
      dd=dd/10
    enddo
    mtype = cdum1
  end subroutine
 
  integer function getelastictype( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    getelastictype = -1
    itype = fetchdigit( 1, mtype )
    if( itype/=1 ) return  ! not defomration problem
    itype = fetchdigit( 2, mtype )
    if( itype/=1 .and. itype/=2 ) return  ! not defomration problem
    getelastictype = fetchdigit( 3, mtype )
  end function
 
  integer function getyieldfunction( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    getyieldfunction = -1
    itype = fetchdigit( 1, mtype )
    if( itype/=1 ) return  ! not defomration problem
    itype = fetchdigit( 2, mtype )
    if( itype/=2 ) return  ! not elstoplastic problem
    getyieldfunction = fetchdigit( 4, mtype )
  end function
 
  integer function gethardentype( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    gethardentype = -1
    itype = fetchdigit( 1, mtype )
    if( itype/=1 ) return  ! not defomration problem
    itype = fetchdigit( 2, mtype )
    if( itype/=2 ) return  ! not elstoplastic problem
    gethardentype = fetchdigit( 5, mtype )
  end function
 
  logical function iskinematicharden( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    iskinematicharden = .false.
    itype = fetchdigit( 5, mtype )
    if( itype==4 .or. itype==5 ) iskinematicharden = .true.
  end function
 
  logical function iselastic( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    iselastic = .false.
    itype = fetchdigit( 2, mtype )
    if( itype==1 ) iselastic = .true.
  end function
 
  logical function iselastoplastic( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    iselastoplastic = .false.
    itype = fetchdigit( 2, mtype )
    if( itype==2 ) iselastoplastic = .true.
  end function
 
  logical function ishyperelastic( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    ishyperelastic = .false.
    itype = fetchdigit( 2, mtype )
    if( itype==3 ) ishyperelastic = .true.
  end function
 
  logical function isviscoelastic( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    isviscoelastic = .false.
    itype = fetchdigit( 2, mtype )
    if( itype==4 ) isviscoelastic = .true.
  end function
 
  subroutine ep2e( mtype )
    integer, intent(inout) :: mtype
    if( .not. iselastoplastic( mtype ) ) return
    call setdigit( 2, 1, mtype )
  end subroutine
 
  integer function getconnectortype( mtype )
    integer, intent(in) :: mtype
    integer :: itype
    getconnectortype = -1
    itype = fetchdigit( 1, mtype )
    if( itype/=1 ) return  ! not defomration problem
    itype = fetchdigit( 2, mtype )
    if( itype/=7 ) return  ! not connector
    getconnectortype = fetchdigit( 3, mtype )
  end function
 
  integer function getnumofspring_dparam( material )
    type( tmaterial ), intent(in) :: material
    getnumofspring_dparam = material%variables_i(m_spring_d_ndoffset+1)
  end function
 
  integer function getnumofspring_aparam( material )
    type( tmaterial ), intent(in) :: material
    getnumofspring_aparam = material%variables_i(m_spring_a_ndoffset+1)
  end function
 
  integer function getnumofdashpot_dparam( material )
    type( tmaterial ), intent(in) :: material
    getnumofdashpot_dparam = material%variables_i(m_dashpot_d_ndoffset+1)
  end function
 
  integer function getnumofdashpot_aparam( material )
    type( tmaterial ), intent(in) :: material
    getnumofdashpot_aparam = material%variables_i(m_dashpot_a_ndoffset+1)
  end function
 
end module