SYSCOMMON_SUBROUTINES.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/>.
 
!---------------------Central Difference Scheme---------------------------
 
subroutine central_difference(NDOF, M, M_inv, C, dT, U, U1, U0, P, flag_Minv)
    
    implicit none
 
    integer, intent(in) :: NDOF
    integer, intent(inout) :: flag_Minv
    integer :: idof, j
 
    real*8, dimension(NDOF, NDOF), intent(in) :: m, c
    real*8, dimension(NDOF), intent(in) :: u1, u0, p
    real*8, intent(in) :: dt
 
    real*8, dimension(NDOF, NDOF), intent(inout) :: m_inv
    real*8, dimension(NDOF), intent(inout) :: u
    real*8, dimension(NDOF, NDOF) :: m1, m2, m3, m4
    real*8, dimension(NDOF) :: t
    real*8 :: x, y
    
    ! Check if for some cases we need to change damping matrix (C); then we can pre define these variables (m1, m2, m3, m4); may be we dont need M,M1_inv,C matrices later
    u=0; 
 
    m1=m/dt/dt + c/2./dt    !!! coeff for u(n+1) -> 1/M_inv
    m2=-2.*m/dt/dt          !!! coeff for u(n)
    m3=m/dt/dt-c/2./dt      !!! coeff for u(n-1)
 
    if (flag_minv.eq.0) then
        if (ndof.gt.1) then
            call matinv(m1, m4, ndof)
        elseif (ndof.eq.1) then
            m4 = 1/m1
        endif
        m_inv = m4
        flag_minv = 1
    else
        m4 = m_inv
    endif
 
    do idof= 1,ndof
        x=0;    y=0;
        do j=1, ndof;
            x = x + m2(idof,j)*u1(j);  
            y = y + m3(idof,j)*u0(j);
        enddo
        t(idof) = p(idof) -x -y
    enddo
 
    !For SDOF
    !T=P-m2*U1-m3*U0
    !U=T/m1
 
    do idof=1,ndof
        do j=1,ndof
            u(idof) = u(idof) + m_inv(idof,j)*t(j)
        enddo
    enddo
 
    return
end subroutine central_difference
 
!------------------Inverse of a Matrix------------------------------------
subroutine matinv(A, B, N) 
 
    implicit none
 
    integer(4), intent(in) :: N
    integer(4) :: IS(N),JS(N)
    integer(4) :: I, J, K
 
    real(8), intent (in) :: A(N,N)
    real(8) :: B(N,N)
    real(8) :: D, T
 
    b=a
    do k=1,n
        d=0.0d0
 
        do i=k,n
            do j=k,n
                if(abs(b(i,j))>d) then
                    d=abs(b(i,j)); is(k)=i; js(k)=j
                end if
            enddo
        enddo
 
        do j=1,n
            t=b(k,j); b(k,j)=b(is(k),j); b(is(k),j)=t
        enddo
 
        do i=1,n
            t=b(i,k);   b(i,k)=b(i,js(k));  b(i,js(k))=t
        enddo
 
        b(k,k)=1.d0/b(k,k);
        do j=1,n; if(j.NE.k) b(k,j)=b(k,j)*b(k,k);  enddo
        
        do i=1,n
            if(i.NE.k) then
                do j=1,n; if(j.NE.k) b(i,j)=b(i,j)-b(i,k)*b(k,j); enddo
            endif
        enddo
 
        do i=1,n; if(i.NE.k) b(i,k)=-b(i,k)*b(k,k); enddo
    enddo
 
    do k=n,1,-1
        do j=1,n
            t=b(k,j); b(k,j)=b(js(k),j); b(js(k),j)=t
        enddo
        do i=1,n
            t=b(i,k); b(i,k)=b(i,is(k)); b(i,is(k))=t
        enddo
    enddo
 
    return
end subroutine matinv
 
 
!------------------Const-Law: Non-linear Shear Spring----------------------------
subroutine ksteel02(props,s,e,de,Et,statev,spd, yield, IDeath) !M, ndof
!
    implicit none
    real*8 e0,sy0,eta,mu,gama,esoft,alpha,beta,a_k,omega
    real*8 emax,emin,ert,srt,erc,src,ehc,eh1,dt,dc,eu
    real*8 de,s,e,s0,et,e_unload,sign,sy,evs,eve,epeak,smax,max
    real*8 sres,eres,x,e_slip,s_slip,e_close,s_close,srel,et1
    real*8 smin,spd,strain_end
    
    !real*8 mu
    integer kon, yield, IDeath  !ndof
    !real*8 M(ndof,ndof)!
    real*8 props(10), statev(11)
 
    e0  = props(1) 
    sy0 = props(2) 
    eta = props(3) 
    mu  = props(4) 
    gama= props(5) 
    esoft=props(6) 
    alpha= props(7) 
    beta = props(8)
    a_k= props(9) 
    omega= props(10) 
    emax  = statev(1) !maximum strain
    emin  = statev(2) !minimum strain
    ert   = statev(3) !strain at load reversal toward tension
    srt   = statev(4)
    erc   = statev(5) !strain at load reversal toward compression
    src   = statev(6)
    kon   = nint(statev(7)) !
    ehc   = statev(8) !effective cummulative hysteresis energy
    eh1   = statev(9) !hysteresis energy in a half cycle
    dt    = statev(10) !damage index for tension
    dc    = statev(11) !damage index for compression
      
    eu    = mu * sy0/e0 !characteristic ultimate strain
 
 
    if(omega<=0.) omega=0.5  
    if(a_k<0.) a_k=0. 
    if(eta<=0.) eta=1.d-6;      
    if(esoft>=0.) esoft=-1.d-6  
      
    if (kon.eq.0) then
        emax =  sy0/e0
        emin = -beta*sy0/e0
        if (de.ge.0.0) then
            kon = 1
        else
            kon = 2
        end if
    else if ((kon.eq.1).and.(de.lt.0.0)) then !Load reversal
            kon = 2
            if (s.gt.0.0) then
                erc = e
                src = s
            end if
            ehc = ehc + eh1 * (erc / eu ) ** 2.0
            eh1 = 0.0 !a new half cycle is to begin
            if (e.gt.emax) emax = e
    else if ((kon.eq.2).and.(de.gt.0.0)) then !Load reversal
            kon = 1
            if (s.lt.0.0) then
                ert = e
                srt = s
            endif
            ehc = ehc + eh1 * (ert / eu ) ** 2.0
            eh1 = 0.0 !a new half cycle is to begin
            if (e.lt.emin) emin = e
    end if
!
    s0=s
    s = s + e0 * de
    et = e0
    
    if(a_k>0.) then
        if(s0>0.) e_unload=e0*(abs(emax/(sy0/e0)))**(-a_k)
        if(s0<0.) e_unload=e0*(abs(emin/(sy0/e0)))**(-a_k) 
        if(e_unload<0.1*e0) e_unload=0.1*e0
    else
        e_unload=e0
    end if
 
    if(s0*de<0.) then 
        s = s0 + e_unload * de
        et = e_unload
        if(s*s0<0.) then 
            de=de-s0/e_unload
            s0=1.d-6*sy0*sign(1.d0,s) 
            et=e0
        end if
    end if
 
        
    if ( de .ge. 0.0 .and. s0>=0.) then
        sy = (1.0 - dt) * sy0
        !loading envelope
        ! Hardening
        if(e+de>sy/e0) then
            evs = max( sy + ( e + de - sy/e0) * eta * e0, 0.)
            eve = eta * e0
           if (s .ge. evs) then
              s = evs
              et = eve
              yield=1
           end if
        end if
        ! Softening
        epeak=sy/e0+(alpha-1.)*sy/e0/eta
        if(e+0.5*de>epeak) then
            evs=max(sy*alpha+esoft*e0*(e+de-epeak),0.0*sy)
            if(sy*alpha+esoft*e0*(e+de-epeak)<=0.) ideath=1 ! Complete damage
            eve=esoft*e0
            if (s .ge. evs) then
               s = evs
               et = eve
               yield=1
            end if
        end if
 
        !reloading envelope
        smax = max(sy, sy + (emax - sy/e0) * eta * e0)  
        if(emax>epeak) then
            smax=max(sy*alpha+esoft*e0*(emax-epeak),0.0*sy)
        end if
        sres = 0.02 * smax                         
        eres = ert - (srt - sres) / e_unload                 
 
        x=emax-smax/e0  
        e_slip=gama*emax+(1.-gama)*x
        s_slip=smax*gama 
        e_close=e_slip*omega 
        s_close=(e_close-eres)/(e_slip-eres) * (s_slip-sres) + sres
 
        if (eres .le. emax - smax / e0) then   
            if(e+0.5*de<e_close)  then  
                srel = (e+de-eres)/(e_slip-eres) * (s_slip-sres) + sres
                et1=(s_slip-sres)/(e_slip-eres)
            else                 
                srel=(e+de-e_close)/(emax-e_close)*(smax-s_close)+s_close
                et1 = (smax - s_close) / (emax - e_close)
            end if
            if (s .gt. srel) then
               s = max( srel, 0.)
               et = et1
            end if
        end if
 
    elseif ( de .lt. 0.0 .and. s0<0. ) then
        sy = (1.0 - dc) * sy0 *beta
        !loading envelope
        ! Hardening
        if(e+de<-sy/e0) then
            evs =  min(-sy + ( e + de + sy/e0) * eta * e0,-0.0*sy)
            eve = eta * e0
            if (s .le. evs) then
                s = evs
                et = eve
                yield=1
            end if
        end if
        ! Softening
        epeak=-sy/e0-(alpha-1.)*sy/e0/eta
        if(e+0.5*de<epeak) then
            evs=min(-sy*alpha+esoft*e0*(e+de-epeak),-0.*sy)
            if(-sy*alpha+esoft*e0*(e+de-epeak)>=0.) ideath=1 
            eve=esoft*e0
            if (s .le. evs) then
                s = evs
                et = eve
                yield=1
          end if
        end if
 
        !reloading envelope
        smin = min(-sy, -sy + (emin + sy/e0) * eta * e0)
        if(emin<epeak)  then 
            smin=min(-sy*alpha+esoft*e0*(emin-epeak),0.)
        end if  
        sres = 0.02 * smin 
        eres = erc - (src - sres) /  e_unload
 
        x=emin-smin/e0 
        e_slip=gama*emin+(1.-gama)*x 
        s_slip=smin*gama 
        e_close=e_slip*omega 
        s_close=(e_close-eres)/(e_slip-eres) * (s_slip-sres) + sres
 
        if (eres .ge. emin - smin / e0) then    
            if(e+0.5*de>e_close) then 
                srel = (e+de-eres)/(e_slip-eres) * (s_slip-sres) + sres
                et1=(s_slip-sres)/(e_slip-eres)
            else         
                srel=(e+de-e_close)/(emin-e_close)*(smin-s_close)+s_close
                et1 = (smin - s_close) / (emin - e_close)
            end if
            if (s .lt. srel) then
                s = min(srel, 0.)
                et = et1
            end if
        end if
    end if
 
 
    if (et.ne.e0 .and. et.ne. e_unload) then 
        spd = spd + s * de
        eh1 = eh1 + s * de
        if ( s .ge. 0.0 ) then
            dc = min(ehc /(3.0 * beta* sy0 * eu), 0.7)
        else
            dt = min(ehc /(3.0 * sy0 * eu), 0.7)
        end if
    end if
 
    x=max(sy0, beta*sy0)
    epeak=x/e0+(alpha-1.)*x/e0/eta
    strain_end=epeak+abs(x/(esoft*e0))
    x=max(abs(emax),abs(emin),abs(e+de))
 
    if (ideath==1) then
        s=0; 
        et=1e-6*e0
    end if
    !
    statev(1)   = emax
    statev(2)   = emin
    statev(3)   = ert
    statev(4)   = srt
    statev(5)   = erc
    statev(6)   = src
    statev(7)   = kon
    statev(8)   = ehc
    statev(9)   = eh1
    statev(10)  = dt
    statev(11)  = dc
    return
end subroutine ksteel02