m32plas.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/.
!||====================================================================
!||    m32plas     ../engine/source/materials/mat/mat032/m32plas.F
!||--- called by ------------------------------------------------------
!||    sigeps32c   ../engine/source/materials/mat/mat032/sigeps32c.F
!||--- calls      -----------------------------------------------------
!||    urotov      ../engine/source/airbag/uroto.F
!||====================================================================
      SUBROUTINE m32plas(JFT    ,JLT    ,PM     ,OFF    ,EPSEQ  ,
     .                   IMAT   ,DIR    ,EZZ    ,IPLA   ,DT1C   ,
     .                   DPLA1  ,EPSPD  ,NEL    ,NGL    ,HARDM  ,
     .                   SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
     .                   DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "impl1_c.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
#include      "com08_c.inc"
#include      "units_c.inc"
#include      "scr17_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER JFT, JLT,IMAT,IPLA,NEL
C     REAL
      my_real
     .   PM(NPROPM,*), OFF(*), DIR(NEL,2), EPSEQ(*),
     .   ezz(*),dt1c(*),dpla1(*),epspd(*),hardm(*),
     .   signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
     .   depsxx(mvsiz),depsyy(mvsiz),depsxy(mvsiz),depsyz(mvsiz),
     .   depszx(mvsiz)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER ICC
      INTEGER I,NMAX,J, NINDX, INDEX(MVSIZ),NGL(MVSIZ),
     .        INDX(MVSIZ),N
C     REAL
      my_real
     .   e,ca, ce, cn, nu,
     .   ymax,cm,epdr,epchk(mvsiz),
     .   a11,a22,a1122,a12,
     .   seq,yld,scale, epsp,s11,s22,s12,g3(mvsiz),nu1,
     .   d11, d22, d12, dpla, dtinv, p,a,b,c,a1s1,a2s2,q,umr,
     .   s_11(mvsiz),s_22(mvsiz),s_12(mvsiz),s1,s2,s3,seqh(mvsiz),
     .   h(mvsiz),dpla_i(mvsiz),dpla_j(mvsiz),etse(mvsiz),nu2(mvsiz),
     .   nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),nnu1(mvsiz),q12(mvsiz),
     .   q21(mvsiz),jq(mvsiz),jq2(mvsiz),a_1,a_2,a_3,
     .   st11(mvsiz),st22(mvsiz),st12(mvsiz),axx(mvsiz),ayy(mvsiz),
     .   a_xy(mvsiz),axy(mvsiz),b_1(mvsiz),b_2(mvsiz),b_3(mvsiz),
     .   yldp(mvsiz),dr,p_1,p_2,p_3,pp1,pp2,pp3,f,df,yld_i,a1,
     .   a2,epsmax,sige(mvsiz,5)
      DATA nmax/3/
C-----------------------------------------------
      e    = pm(20,imat)
      nu   = pm(21,imat)
      ca   = pm(38,imat)
      ce   = pm(39,imat)
      cn   = pm(40,imat)
      ymax = pm(42,imat)
      cm   = pm(43,imat)
      epdr = pm(44,imat)
      a11  = pm(45,imat)
      a22  = pm(46,imat)
      a1122= pm(47,imat)
      a12  = pm(48,imat)
      icc  = nint(pm(49,imat))
      a1   = pm(24,imat)
      a2   = pm(25,imat) 
      epsmax = pm(41,imat)
C
      IF (ipla == 0) THEN
C projection radiale 
C-------------------
C     CONTRAINTE ORTHO
C-------------------
        DO i=jft,jlt
          d11 = dir(i,1)*dir(i,1)
          d22 = dir(i,2)*dir(i,2)
          d12 = dir(i,1)*dir(i,2)
          s11 = d11*signxx(i) + d22*signyy(i) + two*d12*signxy(i)
          s22 = d22*signxx(i) + d11*signyy(i) - two*d12*signxy(i)
          s12 = d12*(signyy(i) - signxx(i)) + ( d11 - d22 )*signxy(i)
C----------------------------
C     VITESSE DE DEFORMATION
C----------------------------
          epsp = max(abs(depsxx(i)),abs(depsyy(i)),half*abs(depsxy(i)))
          dtinv= dt1c(i)/max(dt1c(i)**2,em20)
          epsp = epsp*dtinv
          epsp = max(epsp ,epdr)
          epspd(i) = epsp     
          nnu1(i) = nu / (one - nu)    
          nu5(i) = one-nnu1(i) 
C-------------
C     CRITERE
C-------------
          yld = ca*(ce+epseq(i))**cn*epsp**cm
          yld = min(yld ,ymax)
C------------------------------------------
C     CONTRAINTES PLASTIQUEMENT ADMISSIBLES
C------------------------------------------
          seq = sqrt(a11*s11*s11 + a22*s22*s22
     .            -a1122*s11*s22 + a12*s12*s12)
          scale    = min(one,yld /max(seq ,em20))
          signxx(i)=signxx(i)*scale
          signyy(i)=signyy(i)*scale
          signxy(i)=signxy(i)*scale
          dpla     = off(i)* max(zero,(seq -yld )/e)
          ezz(i)   = - nu5(i)*dpla*half*(signxx(i)+signyy(i))/max(em20,yld)
          epseq(i) = epseq(i) + dpla
          dpla1(i) = dpla
        ENDDO
      ELSEIF (ipla == 2) THEN
C projection avec p=cste + mise a zero de s33 => sig sur le critere 
C-------------------
C     CONTRAINTE ORTHO
C-------------------
        DO i=jft,jlt
          d11 = dir(i,1)*dir(i,1)
          d22 = dir(i,2)*dir(i,2)
          d12 = dir(i,1)*dir(i,2)
          s11 = d11*signxx(i) + d22*signyy(i) + two*d12*signxy(i)
          s22 = d22*signxx(i) + d11*signyy(i) - two*d12*signxy(i)
          s12 = d12*(signyy(i) - signxx(i)) + ( d11 - d22 )*signxy(i)
C----------------------------
C     VITESSE DE DEFORMATION
C----------------------------
          epsp = max(abs(depsxx(i)),abs(depsyy(i)),half*abs(depsxy(i)))
          dtinv= dt1c(i)/max(dt1c(i)**2,em20)
          epsp = epsp*dtinv
          epsp = max(epsp ,epdr)
          epspd(i) = epsp          
C-------------
C     CRITERE
C-------------
          yld =ca*(ce+epseq(i))**cn*epsp**cm
          yld = min(yld ,ymax)
C------------------------------------------
C     CONTRAINTES PLASTIQUEMENT ADMISSIBLES
C------------------------------------------
          nu1 = nu/(one - nu)
          nu5(i) = one-nu1 
          g3(i)  = three_half*e/(one + nu)
          p   = -(s11+s22) * third
          q   = (one -nu1)*p
          s11 = s11+q
          s22 = s22+q
          a1s1 = a11*s11
          a2s2 = a22*s22
          a = a1s1*s11 + a2s2*s22 - a1122*s11*s22 + a12*s12*s12
          b = -q*(a1s1 + a2s2 - half*a1122*(s11+s22))
          c = (a11+a22-a1122)*q*q
          seq = sqrt(a+b+b+c)
          c = c - yld*yld
          scale = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
          umr = one - scale
          q = q*umr
          signxx(i) = signxx(i)*scale - q 
          signyy(i) = signyy(i)*scale - q
          signxy(i) = signxy(i)*scale
          dpla     = off(i)*seq*umr/max(em20,g3(i))
          ezz(i)   = -nu5(i)*dpla*half*(signxx(i)+signyy(i)) / yld
          epseq(i) = epseq(i) + dpla
          dpla1(i) = dpla
        ENDDO
C Projection iterative
      ELSEIF (ipla == 1) THEN
C-- -----------------------
C         CONTRAINTES ORHTO
C-------------------------
        DO i=jft,jlt
          d11 = dir(i,1)*dir(i,1)
          d22 = dir(i,2)*dir(i,2)
          d12 = dir(i,1)*dir(i,2)
          s_11(i) = d11*signxx(i) + d22*signyy(i) + two*d12*signxy(i)
          s_22(i) = d22*signxx(i) + d11*signyy(i) - two*d12*signxy(i)
          s_12(i) = d12*(signyy(i) - signxx(i)) + ( d11 - d22 )*signxy(i)
          nnu1(i) = nu / (one - nu)
C----------------------------
C     VITESSE DE DEFORMATION
C----------------------------
          epsp = max(abs(depsxx(i)),abs(depsyy(i)),half*abs(depsxy(i)))
          dtinv= dt1c(i)/max(dt1c(i)**2,em20)
          epsp = epsp*dtinv
          epsp = max(epsp ,epdr)
          epspd(i) = epsp
C--------------------------------
C     HILL CRITERION
C--------------------------------
          s1=a11*s_11(i)*s_11(i)
          s2=a22*s_22(i)*s_22(i)
          s3=a1122*s_11(i)*s_22(i)
          axy(i)=a12*s_12(i)*s_12(i)
          seqh(i)=sqrt(s1+s2-s3+axy(i))  
          ezz(i) = -(depsxx(i)+depsyy(i))*nnu1(i)
          g3(i)  = three_half*e/(one+nu)
          yldp(i) = ca*(ce+epseq(i))**cn*epsp**cm
          yldp(i) = min(yldp(i) ,ymax)
          IF (yldp(i) >= ymax) THEN
            h(i)=zero
          ELSE
            h(i)=ca*cn*(ce+epseq(i))**(cn-one)*epsp**cm
          ENDIF
        ENDDO
C--------------------------------
C     HARDENING MODULUS
C--------------------------------
        DO i=jft,jlt
          hardm(i) = h(i)
        ENDDO
C-------------------------
C     GATHER PLASTIC FLOW
C-------------------------
        nindx=0
        DO i=jft,jlt
          IF (seqh(i) > yldp(i) .AND. off(i) == one) THEN
            nindx=nindx+1
            index(nindx)=i
          ENDIF
        ENDDO
        IF (nindx > 0) THEN
C---------------------------
C    DEP EN CONTRAINTE PLANE
C---------------------------
#include "vectorize.inc"
          DO  j=1,nindx
            i=index(j)
            dpla_j(i)= (seqh(i)-yldp(i))/(g3(i)+h(i))
            etse(i)= h(i)/(h(i)+e)
            nu2(i) = one-nu*nu
            nu3(i) = nu*half
            nu4(i) = half*(one-nu)
            nu5(i) = one-nnu1(i)
            s1=a11*nu*two-a1122
            s2=a22*nu*two-a1122
            s12=a1122-nu*(a11+a22)
            s3=sqrt(nu2(i)*(a11-a22)*(a11-a22)+s12*s12)
            IF (abs(s1) < em20) THEN 
              q12(i)=zero
            ELSE
              q12(i)=-(a11-a22+s3)/s1
            ENDIF
            IF (abs(s2) < em20) THEN 
              q21(i)=zero 
            ELSE
              q21(i)=(a11-a22+s3)/s2
            ENDIF
            jq(i)=one/(one-q12(i)*q21(i))
            jq2(i)=jq(i)*jq(i)
            a=a11*q12(i)
            b=a22*q21(i)
            a_1=(a11+a1122*q21(i)+b*q21(i))*jq2(i)
            a_2=(a22+a1122*q12(i)+a*q12(i))*jq2(i)
            a_3=(a+b)*jq2(i)*two+a1122*(jq2(i)*two-jq(i))
            st11(i)=s_11(i)+s_22(i)*q12(i)
            st22(i)=q21(i)*s_11(i)+s_22(i)
            axx(i)=a_1*st11(i)*st11(i)
            ayy(i)=a_2*st22(i)*st22(i)
            a_xy(i)=a_3*st11(i)*st22(i)
            a=a1122*nu3(i)
            b=s3*jq(i)
            b_1(i)=a22-a-b
            b_2(i)=a11-a+b
            b_3(i)=a12*nu4(i)
          ENDDO ! DO  J=1,NINDX
C
          DO n=1,nmax
#include "vectorize.inc"
            DO  j=1,nindx
              i=index(j)
              dpla_i(i)=dpla_j(i)
              IF (dpla_i(i) > zero) THEN
                yld_i =min(yldp(i)+h(i)*dpla_i(i),ymax)
                dr  =a1*dpla_i(i)/yld_i
                p_1=one/(one+b_1(i)*dr)
                pp1=p_1*p_1
                p_2=one/(one+b_2(i)*dr)
                pp2=p_2*p_2
                p_3=one/(one+b_3(i)*dr)
                pp3=p_3*p_3
                f  = axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
     .              -yld_i*yld_i
                df = -((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
     .                 (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
     .                  axy(i)*pp3*p_3*b_3(i))*(a1-dr*h(i))/yld_i
     .                  -h(i)*yld_i
                dpla_j(i)=max(zero,dpla_i(i)-f*half/df)
              ELSE
                dpla_j(i)=zero
              ENDIF !IF (DPLA_I(I) > ZERO)
            ENDDO ! DO  J=1,NINDX
          ENDDO ! DO N=1,NMAX
C------------------------------------------
C     CONTRAINTES PLASTIQUEMENT ADMISSIBLES
C------------------------------------------
#include "vectorize.inc"
         DO  j=1,nindx
           i=index(j)
           epseq(i) = epseq(i) + dpla_j(i)
           yld_i =yldp(i)+h(i)*dpla_j(i)
           yld_i =min(yld_i ,ymax)
           dr  =a1*dpla_j(i)/yld_i
           p_1=one/(one+b_1(i)*dr)
           p_2=one/(one+b_2(i)*dr)
           p_3=one/(one+b_3(i)*dr)
            s1=st11(i)*p_1
            s2=st22(i)*p_2
            s_11(i)=jq(i)*(s1-s2*q12(i))
            s_22(i)=jq(i)*(s2-s1*q21(i))
            s_12(i)=s_12(i)*p_3
            s1=a11*s_11(i)+a22*s_22(i)
     .        -a1122*(s_11(i)+s_22(i))*half
            ezz(i) = - nu5(i)*dpla_j(i)*s1/yld_i
            s1 =a1122*half
            p_1=a11*s_11(i)-s1*s_22(i)
            p_2=a22*s_22(i)-s1*s_11(i)
            p_3=a12*s_12(i)
            dpla1(i) = dpla_j(i)
         ENDDO
C
         nindx=0
         DO i=jft,jlt
            IF (epseq(i) > epsmax .AND. off(i) == one) THEN
             nindx=nindx+1
             indx(nindx)=i   
             off(i)=zero
            ENDIF
          ENDDO
C 
          DO i=jft,jlt
            sige(i,1) = s_11(i) 
            sige(i,2) = s_22(i) 
            sige(i,3) = s_12(i) 
            sige(i,4) = zero
            sige(i,5) = zero
          ENDDO
C
          CALL urotov(jft,jlt,sige,dir,nel)
C
          DO i=jft,jlt
           signxx(i)= sige(i,1)
           signyy(i)= sige(i,2)
           signxy(i)= sige(i,3)
          ENDDO
        ENDIF  !  IF (NINDX > 0) THEN
c---
        IF (nindx > 0) THEN
          IF (imconv == 1) THEN
            DO i = 1, nindx
#include   "lockon.inc"
              WRITE(iout, 1000) ngl(indx(i))
              WRITE(istdo,1100) ngl(indx(i)),tt
#include   "lockoff.inc"
            ENDDO
          ENDIF
        ENDIF ! IF (NINDX > 0)
      ENDIF  !  IF (IPLA == 0) THEN
C 
 1000 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'-- RUPTURE OF SHELL ELEMENT :',i10,' AT TIME :',g11.4)
C        
      RETURN
      END