m2lawp.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    m2lawp       ../engine/source/materials/mat/mat002/m2lawp.F
!||--- called by ------------------------------------------------------
!||    main_beam3   ../engine/source/elements/beam/main_beam3.F
!||====================================================================
      SUBROUTINE m2lawp(
     .    PM      ,FOR     ,MOM     ,EINT    ,GEO     ,
     .    OFF     ,PLA     ,EXX     ,EXY     ,EXZ     ,
     .    KXX     ,KYY     ,KZZ     ,AL      ,FA1     ,
     .    FA2     ,FA3     ,MA1     ,MA2     ,MA3     ,
     .    NEL     ,MAT     ,PID     ,NGL     ,IPM     ,
     .    NUMMAT  ,NUVAR   ,UVAR    ,SIGY)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
#include      "param_c.inc"
#include      "scr17_c.inc"
#include      "com08_c.inc"
#include      "impl1_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER ,INTENT(IN) :: NEL,NUMMAT,NUVAR
      INTEGER ,INTENT(IN) :: MAT(NEL),PID(NEL),NGL(NEL)
      INTEGER ,DIMENSION(NPROPMI,NUMMAT) ,INTENT(IN) :: IPM
      my_real
     .   PM(NPROPM,*), FOR(NEL,3), MOM(NEL,3), EINT(NEL,2), GEO(NPROPG,*),
     .   OFF(*), PLA(*),
     .   al(nel),
     .   exx(nel),exy(nel),exz(nel),kxx(nel),kyy(nel),
     .   kzz(nel),fa1(nel),fa2(nel),fa3(nel),
     .   ma1(nel),ma2(nel),ma3(nel),a1(nel)
      my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) :: uvar
      my_real ,DIMENSION(NEL)       ,INTENT(INOUT) :: sigy
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER INDX(NEL),ICC(NEL),IRTY(NEL)
      INTEGER  I,NPIF, J, NINDX
      my_real
     .   yldtmp(nel),ym(nel), 
     .   b1(nel), b2(nel), b3(nel), degmb(nel),
     .   degfx(nel),esp(nel),rho(nel),g(nel),
     .   dmpm(nel), dmpf(nel), shf(nel),
     .   f1(nel), m1(nel), m2(nel), m3(nel),
     .   degsh(nel), yeq(nel), yld(nel), dwpla(nel),
     .   epmx(nel), dwelm(nel), dwelf(nel),z3(nel),z4(nel),
     .   ca(nel), cb(nel), cn(nel), ymax(nel),
     .   rr(nel), cc(nel), epdr(nel), epsp(nel),
     .   sh(nel), yma2(nel), sh10(nel), sh20(nel), sh0(nel),
     .   sh1(nel), sh2(nel), fact, epif,plap(nel),dpla(nel),
     .   tmp1(nel), tmp2(nel), tmp3(nel), etmp,plap1,vp,asrate,
     .   factphi,yma2phi, sh10p, sh20p, sh1p, sh2p,shp, sh0p,
     .   ys,ysp,gs,gsp,trm1,trm1p,trm2,trm2p,trm3,trm3p,yeq0,yeq1
C-----------------------------------------------
      epif = zero
      npif = 0
C
      IF (impl_s == 0 .OR. idyna > 0) THEN
        DO i=1,nel                    
          dmpm(i)=geo(16,pid(i))*al(i)  
          dmpf(i)=geo(17,pid(i))*al(i)  
        ENDDO
      ELSE
        DO i=1,nel                    
          dmpm(i)=zero  
          dmpf(i)=zero  
        ENDDO 
      ENDIF                                 
C
 
      DO i=1,nel
        vp     = ipm(255,mat(i)) 
        rho(i) = pm( 1,mat(i))
        asrate = min(one, pm(9,mat(i))*dt1)
        g(i)   =pm(22,mat(i))
        ym(i)  =pm(20,mat(i))
C--------------------------
        ca(i)  =pm(38,mat(i))
        cb(i)  =pm(39,mat(i))
        cn(i)  =pm(40,mat(i))
        epmx(i)=pm(41,mat(i))
        ymax(i)=pm(42,mat(i))
        cc(i)  =pm(43,mat(i))
        IF(vp == 1)THEN
          epdr(i) =   max(em20,pm(44,mat(i)))
        ELSE               
          epdr(i) = max(em20,pm(44,mat(i))*dt1)  
        ENDIF
        icc(i) =nint(pm(49,mat(i)))
C-----------------------------
        a1(i)  =geo(1,pid(i))
        b1(i)  =geo(2,pid(i))
        b2(i)  =geo(18,pid(i))
        b3(i)  =geo(4,pid(i))
C        SHF(I) =GEO(13,PID(I))
        shf(i) =geo(37,pid(i))
        epif    = max(epif,cc(i))
        irty(i) = nint(pm(50,mat(i)))
        npif    = npif+irty(i)
        z3(i)   =pm(51,mat(i))
        z4(i)   =pm(52,mat(i))
      ENDDO
C-----------------------------
C     DAMPING terms removed to pforce3
C-----------------------------
      DO i=1,nel
        esp(i) = (eint(i,1)+eint(i,2))/al(i)/a1(i)
      ENDDO
C
      DO i=1,nel
        degmb(i) = for(i,1)*exx(i)
        degsh(i) = for(i,2)*exy(i)+for(i,3)*exz(i)
        degfx(i) = mom(i,1)*kxx(i)+mom(i,2)*kyy(i)+mom(i,3)*kzz(i)
      ENDDO
C
C     CISAILLEMENT TRANSVERSAL CALCULE AVEC K1=12EI/L**2 K2=5/6GA
C
      DO i=1,nel
        sh(i)=five_over_6*g(i)*a1(i)
        yma2(i)=twelve*ym(i)/al(i)**2
        sh10(i)=yma2(i)*b1(i)
        sh20(i)=yma2(i)*b2(i)
        sh0(i)=(one - shf(i))*sh(i)
        sh1(i)=sh0(i)*sh10(i)/(sh(i)+sh10(i)) + shf(i)*sh10(i)
        sh2(i)=sh0(i)*sh20(i)/(sh(i)+sh20(i)) + shf(i)*sh20(i)
C
        f1(i)   =for(i,1)+ exx(i)*a1(i)*ym(i)
        for(i,2)=for(i,2)+ exy(i)*sh2(i)
        for(i,3)=for(i,3)+ exz(i)*sh1(i)
        m1(i)   =mom(i,1)+ kxx(i)*g(i) *b3(i)
        m2(i)   =mom(i,2)+ kyy(i)*ym(i)*b1(i)
        m3(i)   =mom(i,3)+ kzz(i)*ym(i)*b2(i)
      ENDDO
C-------------
C     CRITERE
C-------------
      DO i=1,nel
        yeq(i)= f1(i)*f1(i) + three * a1(i) *
     +                    ( m1(i)*m1(i) / max(b3(i),em20)
     +                    + m2(i)*m2(i) / max(b1(i),em20)
     +                    + m3(i)*m3(i) / max(b2(i),em20) )
        yeq(i)= sqrt(yeq(i))/a1(i)
      ENDDO
C-------------
C     STRAIN RATE (JOHNSON-COOK, ZERILLI-ARMSTRONG)
C-------------
      IF (epif /= zero) THEN
        DO i=1,nel
          IF(vp == 1)THEN
            plap(i) = uvar(i,1)
            plap(i) = max(plap(i),epdr(i))
            epsp(i) = log(plap(i)/epdr(i))
          ELSE
            epsp(i)=abs(degmb(i)+degfx(i))/(yeq(i)+ em20)/a1(i)
            tmp2(i)=abs(degmb(i)+degfx(i))
            tmp3(i)=epsp(i)
            epsp(i)= max(epsp(i),epdr(i))
            epsp(i)= log(epsp(i)/epdr(i))
          ENDIF
        ENDDO
        IF (npif == zero) THEN
          DO i=1,nel
            epsp(i)=(one + cc(i) * epsp(i))
            IF (icc(i) == 1) ymax(i) = ymax(i) * epsp(i)
          ENDDO
        ELSEIF (npif == nel) THEN
          DO i=1,nel
            epsp(i)=(one + cc(i) * epsp(i))
            tmp1(i)=epsp(i)
            epsp(i)= cc(i)*exp((-z3(i)+z4(i) * epsp(i))*esp(i))
            IF(icc(i)==1)ymax(i)= ymax(i) + epsp(i)
            ca(i) = ca(i) + epsp(i)
            epsp(i)=one
          ENDDO
        ELSE
          DO i=1,nel
            IF (irty(i) == 0)THEN
              epsp(i)=(one + cc(i) * epsp(i))
              tmp1(i)=ymax(i)
              IF(icc(i)==1) ymax(i) = ymax(i) * epsp(i)
            ELSE
              epsp(i)=(one  + cc(i) * epsp(i))
              tmp1(i)=epsp(i)
              epsp(i)= cc(i)*exp((-z3(i)+z4(i) * epsp(i))*esp(i))
              IF(icc(i)==1)ymax(i)= ymax(i) + epsp(i)
              ca(i) = ca(i) + epsp(i)
              epsp(i)=one
            ENDIF
          ENDDO
        ENDIF
      ELSE
        DO i=1,nel
          epsp(i)=one
        ENDDO
      ENDIF
C-----------------------------------
C     YIELD
C-----------------------------------
      DO i=1,nel
        yld(i)= ca(i) + cb(i) * exp(cn(i) * log(pla(i)+em30))
        yldtmp(i)=yld(i)
        yld(i)= min(yld(i)*epsp(i),ymax(i))
        rr(i) = min(one,yld(i)/(yeq(i)+ em20))
        sigy(i) = yld(i)
      ENDDO
C
      DO i=1,nel
        f1(i)    = f1(i)*rr(i)
        dwelm(i) =(f1(i)+for(i,1))*(f1(i)-for(i,1))/ym(i)/a1(i)
        degmb(i) = degmb(i) + f1(i)*exx(i)
      ENDDO
C
      DO i=1,nel
        m1(i) = m1(i)*rr(i)
        m2(i) = m2(i)*rr(i)
        m3(i) = m3(i)*rr(i)
        dwelf(i) =(m1(i)+mom(i,1))*(m1(i)-mom(i,1))/ g(i)/b3(i)+
     .            (m2(i)+mom(i,2))*(m2(i)-mom(i,2))/ym(i)/b1(i)+
     .            (m3(i)+mom(i,3))*(m3(i)-mom(i,3))/ym(i)/b2(i)
        degfx(i) = degfx(i)+ m1(i)*kxx(i)+m2(i)*kyy(i)+m3(i)*kzz(i)
      ENDDO
C
      DO i=1,nel
        dwpla(i)= degmb(i)+degfx(i)-dwelm(i)-dwelf(i)
      ENDDO
C-----------------------
C     EPS PLASTIQUE
C-----------------------
      DO i=1,nel
        tmp1(i)=epsp(i)*dwpla(i)/yld(i)
c       PLA(I)=PLA(I)+OFF(I)*MAX(ZERO,0.5*EPSP(I)*DWPLA(I)/YLD(I)/A1(I))
        dpla(i) = max(zero,half*tmp1(i)/a1(i))
        pla(i)  = pla(i)+off(i) * dpla(i)
      ENDDO
      DO i=1,nel
        IF (vp == 1) THEN 
          plap1     = dpla(i) / max(em20,dt1)
          uvar(i,1) = asrate * plap1 + (one - asrate) * plap(i)
        ENDIF
      ENDDO
 
C--------------------------------
C     TEST DE RUPTURE DUCTILE
C--------------------------------
      DO i=1,nel
          IF (off(i) < em01)  off(i) = zero
          IF (off(i) < one )  off(i) = off(i)*four_over_5
      ENDDO
C
      nindx  = 0
      DO i=1,nel
        IF (off(i) < one)      cycle
        IF (pla(i) < epmx(i))  cycle
        off(i)=four_over_5
        idel7nok = 1
        nindx=nindx+1
        indx(nindx)=i
      ENDDO
C
      IF (nindx > 0 .AND. imconv == 1) THEN
        DO j=1,nindx
#include "lockon.inc"
          WRITE(iout, 1000) ngl(indx(j))
          WRITE(istdo,1100) ngl(indx(j)),tt
#include "lockoff.inc"
        ENDDO
      ENDIF
C
      DO i=1,nel
        fa1(i) = f1(i)
        fa2(i) = for(i,2) 
        fa3(i) = for(i,3)  
        ma1(i) = m1(i) 
        ma2(i) = m2(i) 
        ma3(i) = m3(i) 
      ENDDO               
C
 1000 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT :',i10,' AT TIME :',g11.4)
C-----------------------------------------------
      RETURN
      END