SDOF_SHEAR_MODEL.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
 
subroutine sdof_shear_model (sID, ndof, gr_acc, direction)
 
  use speed_sci
 
  implicit none
 
  integer*4 :: sID, ndof, direction, Usign
  integer*4 :: idof, i, j
  real*8 :: gr_acc
  real*8, dimension(ndof) :: e, de, s, syspeff0, syspeff
 
  syspeff0 = 0
 
  do idof=1,ndof
    syspeff(idof) = gr_acc*sys(sid)%Ms(idof,idof)     !!!mass x (-u_g'') = inertial force from the ground acceleration
  enddo
 
  if (ndof.gt.1) then
    do idof=1,ndof-1
      syspeff0(idof) = syspeff(idof) - sys(sid)%IntForce(idof,direction) + sys(sid)%IntForce(idof+1,direction)     !!! SDOFItF = interaction force Ku(n)
    enddo
  endif
  syspeff0(ndof) = syspeff(ndof) - sys(sid)%IntForce(ndof,direction)
 
  ! returns sys(sID)%tempU = displacement at time n+1 = Relative displacement of structure w.r.t. ground
  if (sys(sid)%StructType.eq.1) then
    call central_difference(ndof, sys(sid)%Ms(1,1), sys(sid)%Ms_inv(1,1), sys(sid)%Cs, sys(sid)%dt, &
                            sys(sid)%tempU(1,direction), sys(sid)%tempU1(1,direction), &
                            sys(sid)%tempU0(1,direction), syspeff0, sys(sid)%flag_Minv )
  else
    call central_difference(ndof, sys(sid)%Ms(1:ndof,1:ndof), sys(sid)%Ms_inv(1:ndof,1:ndof), sys(sid)%SysC(1:ndof,1:ndof), &
                            sys(sid)%dt, sys(sid)%tempU(1:ndof,direction), sys(sid)%tempU1(1:ndof,direction), &
                            sys(sid)%tempU0(1:ndof,direction), syspeff0, sys(sid)%flag_Minv )
    ! if (direction.eq.1) then
    !   write(*,*) 'InputForces: Acc =', gr_acc, 'sysPeff = ', (sysPeff0(idof), idof = 1,sys(1)%NDOF)
    !   write(*,*) 'Input Un-1 = ', (sys(1)%tempU0(idof,1), idof=1,sys(1)%NDOF), 'Un = ', (sys(1)%tempU1(idof,1), idof=1,sys(1)%NDOF)
    !   write(*,*) 'Output Un+1 = ', (sys(1)%tempU(idof,1), idof=1,sys(1)%NDOF)
    ! endif
  endif
 
  ! Modify this for SSI - with MDOF (below code works for SSI-SDOF)
  if(sys(sid)%const_law.eq.3) then
    if (sys(sid)%tempU(1, direction).ge.0.0) usign = 1
    if (sys(sid)%tempU(1, direction).lt.0.0) usign = -1
    if ((abs(sys(sid)%tempU(1, direction)).ge.sys(sid)%EU).or.(sys(sid)%damage(direction).eq.1)) &
        sys(sid)%tempU(1,direction) = usign*sys(sid)%EU
  endif
 
  ! Time Variation of Interstory Drift (IDR)
  sys(sid)%variIDR(1,direction) = sys(sid)%tempU(1,direction) - sys(sid)%IDR(1,direction)
  sys(sid)%IDR(1,direction)=sys(sid)%tempU(1,direction); 
  if (sys(sid)%StructType.eq.2) then
    do j=2,ndof
        sys(sid)%variIDR(j,direction) = (sys(sid)%tempU(j,direction) - sys(sid)%tempU(j-1,direction)) - sys(sid)%IDR(j,direction)
        sys(sid)%IDR(j,direction) = sys(sid)%tempU(j,direction) - sys(sid)%tempU(j-1,direction)
    enddo
    !Max IDR 
    do j=1,ndof
        sys(sid)%MaxIDR(j,direction) = max(sys(sid)%MaxIDR(j,direction), abs(sys(sid)%IDR(j,direction)/sys(sid)%Floor_h))
    enddo
  endif
 
  do idof= 1,ndof
      !!! Relative acceleration
      sys(sid)%tempRA1(idof,direction) = ( sys(sid)%tempU(idof,direction) + sys(sid)%tempU0(idof,direction) - &
                                          2.d0*sys(sid)%tempU1(idof, direction))/sys(sid)%dt2
      !!! Absolute acceleration AH
      sys(sid)%tempA1(idof,direction) = sys(sid)%tempRA1(idof,direction) - gr_acc   
 
      ! Updating displacement Values
      sys(sid)%tempU0(idof,direction)=sys(sid)%tempU1(idof,direction)
      sys(sid)%tempU1(idof,direction)=sys(sid)%tempU(idof,direction)
      sys(sid)%tempU(idof,direction)=0
  enddo
 
  do idof=1,ndof
      e(idof) = sys(sid)%IDR(idof,direction) - sys(sid)%variIDR(idof,direction)       ! IDR at (n)th time step
      de(idof)= sys(sid)%variIDR(idof,direction)    ! change in IDR at (n+1)th timestep
      s(idof) = sys(sid)%IntForce(idof,direction)   ! Shear Force at (n)th timestep
 
      if (sys(sid)%StructType.eq.1) then
          if (sys(sid)%const_law.eq.1) then
            call linear_elastic(sys(sid)%Ks, s, de)      !!! elastic constitutive law
          elseif (sys(sid)%const_law.eq.2) then
            call perfectly_plastic(sys(sid)%Ks, s, de, sys(sid)%FY)      !!! elastoplastic constitutive law
          elseif (sys(sid)%const_law.eq.3) then
            call trilinear(sys(sid)%Ks, sys(sid)%Hs, sys(sid)%Ss, s, de, e, sys(sid)%FY, sys(sid)%FH, sys(sid)%FU, sys(sid)%EY, &
                sys(sid)%EH, sys(sid)%EU, sys(sid)%branch(direction), sys(sid)%damage(direction))      !!! trilinear constitutive law
          endif
      elseif (sys(sid)%StructType.eq.2) then
 
        if (sys(sid)%const_law.eq.1) then
          call linear_elastic(sys(sid)%props(idof,1,1), s(idof), de(idof))
        elseif (sys(sid)%const_law.eq.2) then
          if ((abs(sys(sid)%IDR(idof,direction))/sys(sid)%Floor_h).GT.sys(sid)%dval(idof,4)) then
              sys(sid)%MDOFIDeath(idof,direction) = 1
          endif
          call ksteel02(sys(sid)%props(idof,1:10,1), s(idof), e(idof), de(idof), sys(sid)%MDOFEt(idof,direction), &
                        sys(sid)%MDOFstatev(idof,1:11,direction), sys(sid)%MDOFspd(idof,direction), &
                        sys(sid)%MDOFyield(idof,direction), sys(sid)%MDOFIDeath(idof,direction) ) !sys(sID)%Ms(1:ndof,1:ndof), ndof
        endif
      endif
 
      sys(sid)%IntForce(idof,direction) = s(idof)
 
  enddo
 
  return
end subroutine sdof_shear_model