SDOF_SHEAR_MODEL.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	sdof_shear_model (sid, ndof, gr_acc, direction)
	Computation for SDOF shear model.

Function/Subroutine Documentation

sdof_shear_model()

subroutine sdof_shear_model	(	integer*4	sid,
		integer*4	ndof,
		real*8	gr_acc,
		integer*4	direction
	)

Computation for SDOF shear model.

Author

Aline Herlin

Date

November, 2020

Version

1.0

Parameters

[in]	sID	oscillator ID
[in]	gr_acc	base acceleration
[in]	direction	direction of motion

Definition at line 27 of file SDOF_SHEAR_MODEL.f90.

 
  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

References central_difference(), ksteel02(), linear_elastic(), perfectly_plastic(), speed_sci::sys, and trilinear().

Here is the call graph for this function: