fstr_contact_damping.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 m_fstr_contact_damping
  use hecmw
  use elementinfo
  use mcontactdef
  use m_fstr_contact_elem_common
  implicit none
 
  private
  public :: is_damping_enabled
  public :: contact_damping_weight
  public :: getdampingstiffness
  public :: getdampingnodalforce
 
contains
 
  pure logical function is_damping_enabled(contact)
    type(tcontact), intent(in) :: contact
    is_damping_enabled = (contact%damp_alpha > 0.0d0 .and. contact%damp_gact > 0.0d0)
  end function is_damping_enabled
 
  pure real(kind=kreal) function contact_damping_weight(g, gact)
    real(kind=kreal), intent(in) :: g     
    real(kind=kreal), intent(in) :: gact  
    real(kind=kreal) :: t
 
    if (g <= 0.0d0) then
      contact_damping_weight = 1.0d0
    else if (g >= gact) then
      contact_damping_weight = 0.0d0
    else
      t = g / gact
      contact_damping_weight = 1.0d0 - 3.0d0*t*t + 2.0d0*t*t*t
    end if
  end function contact_damping_weight
 
  subroutine getdampingstiffness(cstate, surf, alpha, gact, stiff, smoothing_type)
    type(tcontactstate), intent(in)    :: cstate
    type(tsurfelement), intent(in)     :: surf
    real(kind=kreal), intent(in)       :: alpha         
    real(kind=kreal), intent(in)       :: gact          
    real(kind=kreal), intent(out)      :: stiff(:,:)    
    integer(kind=kint), optional, intent(in) :: smoothing_type
 
    integer(kind=kint) :: nnode, ndof, i, j
    real(kind=kreal) :: bn(l_max_surface_node*3+3)
    real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
    real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))
    real(kind=kreal) :: w, coeff
 
    nnode = size(surf%nodes)
    ndof = nnode*3 + 3
 
    stiff = 0.0d0
 
    ! Compute Bn directly via computeTm_Tt (no need for metric/Ht/Gt)
    call computetm_tt(cstate, surf, 0.0d0, tm, tt, smoothing_type, bn=bn)
 
    ! Compute weight from frozen distance
    w = contact_damping_weight(cstate%distance, gact)
    coeff = alpha * w
 
    if (coeff <= 0.0d0) return
 
    ! Build stiffness: K_damp = coeff * Bn (x) Bn
    do j = 1, ndof
      do i = 1, ndof
        stiff(i,j) = coeff * bn(i) * bn(j)
      enddo
    enddo
 
  end subroutine getdampingstiffness
 
  subroutine getdampingnodalforce(cstate, surf, ndDu, alpha, gact, ctNForce, ctTForce, smoothing_type)
    type(tcontactstate), intent(in)    :: cstate
    type(tsurfelement), intent(in)     :: surf
    real(kind=kreal), intent(in)       :: nddu(:)       
    real(kind=kreal), intent(in)       :: alpha         
    real(kind=kreal), intent(in)       :: gact          
    real(kind=kreal), intent(out)      :: ctnforce(:)   
    real(kind=kreal), intent(out)      :: cttforce(:)   
    integer(kind=kint), optional, intent(in) :: smoothing_type
 
    integer(kind=kint) :: nnode, ndof
    real(kind=kreal) :: bn(l_max_surface_node*3+3)
    real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
    real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))
    real(kind=kreal) :: w, coeff, delta_gn
 
    nnode = size(surf%nodes)
    ndof = nnode*3 + 3
 
    ctnforce = 0.0d0
    cttforce = 0.0d0
 
    ! Compute Bn directly via computeTm_Tt (no need for metric/Ht/Gt/elemcrd)
    call computetm_tt(cstate, surf, 0.0d0, tm, tt, smoothing_type, bn=bn)
 
    ! Compute weight from frozen distance
    w = contact_damping_weight(cstate%distance, gact)
    coeff = alpha * w
 
    if (coeff <= 0.0d0) return
 
    ! Compute normal gap increment: delta_gn = Bn^T * ddu
    delta_gn = dot_product(bn(1:ndof), nddu(1:ndof))
 
    ! Damping force: f_damp = coeff * delta_gn * Bn
    ctnforce(1:ndof) = coeff * delta_gn * bn(1:ndof)
 
    ! Lagrange row (not used, set to 0 for interface compatibility)
    ctnforce((nnode+1)*3+1) = 0.0d0
    cttforce((nnode+1)*3+1) = 0.0d0
 
  end subroutine getdampingnodalforce
 
end module m_fstr_contact_damping