ruser36.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/.
!||====================================================================
!||    ruser36      ../engine/source/elements/spring/ruser36.F
!||--- called by ------------------------------------------------------
!||    rforc3       ../engine/source/elements/spring/rforc3.F
!||--- calls      -----------------------------------------------------
!||    get_u_func   ../engine/source/user_interface/ufunc.F
!||    get_u_mid    ../engine/source/user_interface/upidmid.F
!||    get_u_mnu    ../engine/source/user_interface/upidmid.F
!||    get_u_pid    ../engine/source/user_interface/upidmid.F
!||    get_u_pnu    ../engine/source/user_interface/upidmid.F
!||====================================================================
      SUBROUTINE ruser36(
     1             NEL     ,IPROP   ,UVAR   ,NUVAR  ,FR_WAVE,
     2             FX      ,FY      ,FZ     ,XMOM   ,YMOM   ,
     3             ZMOM    ,E       ,OFF    ,STIFM  ,STIFR  ,
     4             VISCM   ,VISCR   ,MASS   ,XINER  ,DT     ,
     5             XL      ,VX      ,RY1    ,RZ1    ,RX     ,
     6             RY2     ,RZ2     )
C-------------------------------------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C----------------------------------------------------------
C   d u m m y   a r g u m e n t s   a n d   f u n c t i o n
C----------------------------------------------------------
      INTEGER NEL,NUVAR,IPROP,
     .        GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
     .        KFUNC,KMAT,KPROP
      MY_REAL
     .   UVAR(NUVAR,*),DT ,
     .   FX(*), FY(*), FZ(*), E(*), VX(*),MASS(*) ,XINER(*),
     .   ry1(*), rz1(*), off(*), xmom(*), ymom(*),
     .   zmom(*), rx(*), ry2(*), rz2(*),xl(*),
     .   stifm(*) ,stifr(*) , viscm(*) ,viscr(*) ,fr_wave(*) ,
     .   get_u_mat, get_u_geo, get_u_func
      EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
     .         get_u_mat,get_u_geo, get_u_func
      parameter(kfunc=29)
      parameter(kmat=31)
      parameter(kprop=33)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I,KCST,N0,
     . imat1,iprop1,
     . imat2,iprop2
      INTEGER J,N,NINDX,NMAX,INDEX(MVSIZ),IFUNC
      MY_REAL 
     .        FAC,RHO,AREA,IXX,IYY,IZZ,IMYZ,YOUNG,G,
     .        area1,ixx1,iyy1,izz1,rho1,young1,g1,
     .        area2,ixx2,iyy2,izz2,rho2,young2,g2
      my_real r_y1,r_y2,r_z1,r_z2,r_y,r_z,r_x,facx,facy,facz,
     .        ay1,ay2,az1,az2,ay,az,by1,by2,bz1,bz2,by,bz,g3,
     .        cx1,cx2,cx,area_y(mvsiz),area_z(mvsiz),
     .        yld(mvsiz),dre,drg,dpla_i,dpla_j(mvsiz),dfe,dfg,
     . err,f,df,yld_i,c1,fn,fm,y1,y2,
     .        pnx,pnx2,pny,pny2,pnz,pnz2,pmx,pmx2,pmy,pmy2,pmz,pmz2,
     . nx2,ny2,nz2,mx2,my2,mz2,tempy,tempz,ktran,krot,
     .        nx(mvsiz),ny(mvsiz),nz(mvsiz),youngm(mvsiz),gm(mvsiz),
     .        mx(mvsiz),my(mvsiz),mz(mvsiz),svm(mvsiz),h(mvsiz),
     .        epsvxx(mvsiz),epsvxy(mvsiz),epsvxz(mvsiz),
     .        epscbx(mvsiz),epscby(mvsiz),epscbz(mvsiz),
     .        signxx(mvsiz),signxy(mvsiz),signxz(mvsiz),
     . momnxx(mvsiz),momnyy(mvsiz),momnzz(mvsiz),
     . gama(mvsiz),gamac(mvsiz),xl2(mvsiz),devol,
     .        pc1,dc1,pr1,ps1,dc2,pr2,ps2,dc,pr,ps,ueq,dtemp,
     .        sig01,sig02, sig0,hpla,h1,h2,m1,m2,m,dsigy,cc1,
     .        fac1,fac2
      DATA nmax/10/
C-----------------------------------------------
      kcst = 0
C
      iprop1 = get_u_pnu(1,iprop,kprop) 
      area1  = get_u_geo(2,iprop1)
      ixx1   = get_u_geo(3,iprop1)
      iyy1   = get_u_geo(4,iprop1)
      izz1   = get_u_geo(5,iprop1)
      r_y1= get_u_geo(6,iprop1)
      r_z1= get_u_geo(7,iprop1)
      iprop2 = get_u_pnu(2,iprop,kprop) 
      area2  = get_u_geo(2,iprop2)
      ixx2   = get_u_geo(3,iprop2)
      iyy2   = get_u_geo(4,iprop2)
      izz2   = get_u_geo(5,iprop2)
      imat1  = get_u_pnu(1,iprop1,kmat) 
      ifunc  = get_u_mnu(1,imat1,kfunc)
      young1 = get_u_mat(7,imat1)
      g1     = get_u_mat(6,imat1)
      rho1   = get_u_mat(0,imat1)
      ay1    = get_u_mat(8,imat1)
      az1    = get_u_mat(9,imat1)
      by1    = get_u_mat(10,imat1)
      bz1    = get_u_mat(11,imat1)
      cx1    = get_u_mat(12,imat1)
      dc1    = get_u_mat(13,imat1)
      pr1    = get_u_mat(14,imat1)
      ps1    = get_u_mat(15,imat1)
      sig01  = get_u_mat(16,imat1)
      h1     = get_u_mat(17,imat1)
      m1     = get_u_mat(18,imat1)
      fac1   = get_u_mat(19,imat1)
      imat2  = get_u_pnu(1,iprop2,kmat)
      young2 = get_u_mat(7,imat2)
      g2     = get_u_mat(6,imat2)
      rho2   = get_u_mat(0,imat2)
      ay2    = get_u_mat(8,imat2)
      az2    = get_u_mat(9,imat2)
      by2    = get_u_mat(10,imat2)
      bz2    = get_u_mat(11,imat2)
      cx2    = get_u_mat(12,imat2)
      dc2    = get_u_mat(13,imat2)
      pr2    = get_u_mat(14,imat2)
      ps2    = get_u_mat(15,imat2)
      sig02  = get_u_mat(16,imat2)
      h2     = get_u_mat(17,imat2)
      m2     = get_u_mat(18,imat2)
      r_y2   = get_u_geo(6,iprop2)
      r_z2   = get_u_geo(7,iprop2)
      rho    = half*(rho1+rho2)
      area   = half*(area1+area2)
      fac    = rho*area
      ixx    = half*(ixx1+ixx2)
      iyy    = half*(iyy1+iyy2)
      izz    = half*(izz1+izz2)
      r_y = half*(r_y1+r_y2)
      r_z = half*(r_z1+r_z2)
      r_x = half*(r_y+r_z)
      cx  = half*(cx1+cx2)
      ay  = half*(ay1+ay2)
      az  = half*(az1+az2)
      by  = half*(by1+by2)
      bz  = half*(bz1+bz2)     
      imyz   = max(iyy,izz)
      young  = half*(young1+young2)
      g      = half*(g1+g2)      
      facx   = r_x/max(ixx,em20)
      facy   = r_y/max(iyy,em20)
      facz   = r_z/max(izz,em20)
      tempy  = g*area/max(twelve*young*iyy,em20)
      tempz  = g*area/max(twelve*young*izz,em20)
      g3     = three*g
      n0     = 11
      dc     = half*(dc1+dc2)
      pr     = half*(pr1+pr2)
      ps     = half*(ps1+ps2)
      sig0   = half*(sig01+sig02)
      hpla   = half*(h1+h2)
      m     = half*(m1+m2)     
C
C----  UVAR(N0,I) -> PLA , N0+1 -> SY0 , N0+2 -> EXX_P , N0+3 -> EYY_P
C----  UVAR(N0+4,I) -> ENDOMMAGEMENT
C
      DO i=1,nel
        youngm(i)  = young*(1-uvar(n0+4,i))
        gm(i)      = g*(1-uvar(n0+4,i))
C== WE FIRST CALCULATE THE VELOCITIES OF GENERALIZED DEFORMATIONS ====
        dtemp = one/max(xl(i),em20)
        epsvxx(i) = vx(i)*dtemp
        epsvxy(i) = -half*(rz1(i) + rz2(i))
        epsvxz(i) = half*(ry1(i) + ry2(i))
        epscbx(i) = r_x*rx(i)*dtemp
        epscby(i) = r_y*(ry2(i) - ry1(i))*dtemp
        epscbz(i) = r_z*(rz2(i) - rz1(i))*dtemp
C== ON CALCULE CONTRAINTES GENERALISEES ====
        xl2(i)    =  xl(i)*xl(i)
        area_y(i) = area/(one +tempy*xl2(i))    
        area_z(i) = area/(one +tempz*xl2(i))
        signxx(i) = fx(i)/max(area,em20)+youngm(i)*epsvxx(i)*dt
        signxy(i) = fy(i)/max(area_y(i),em20)+gm(i)*epsvxy(i)*dt
        signxz(i) = fz(i)/max(area_z(i),em20)+gm(i)*epsvxz(i)*dt
        momnxx(i) = xmom(i)*facx+gm(i)*epscbx(i)*dt
        momnyy(i) = ymom(i)*facy+youngm(i)*epscby(i)*dt
        momnzz(i) = zmom(i)*facz+youngm(i)*epscbz(i)*dt
      ENDDO
C
C--- UVAR(11,I)-> DEF PLASTIQUE EQUIVALENTE -
C--- UVAR(N0,I)-> PLA, N0+1->SY0; N0+2->EXX_P; N0+3->EYY_P;
C--- UVAR(N0+3,I)-> COURBURE PLASTIQUE Y,Z UVAR(N0+2,I)-> COURBURE PLASTIQUE X -
C
      DO i=1,nel
        IF (ifunc /= 0) THEN
          y1= fac1*get_u_func(ifunc,uvar(n0,i),y2)
          y2= fac1*y2   
          yld(i) = max(y1,em20)                
          m=1
        ELSE
          y1 = sig0+hpla*(uvar(n0,i))**m
          y2 = hpla
        ENDIF
        yld(i) = max(y1,em20)
        h(i)   = max(y2,zero)
        uvar(n0+1,i)=min(yld(i),uvar(n0+1,i))
C
        pc1 = abs(six*youngm(i)*(uvar(n0+2,i)/uvar(n0+1,i)))
        c1=one-fourth*exp(-pc1)
        gamac(i)=cx*three_over_4/max(c1,em20)
C
        pc1 = six*youngm(i)*(uvar(n0+3,i)/uvar(n0+1,i))
        cc1 =  (three*pi/sixteen)**2
        c1=one - (one - cc1)*exp(-pc1)
        gama(i)=cc1/max(c1,em20)
C
        nx(i)=signxx(i)
        ny(i)=ay*signxy(i)
        nz(i)=az*signxz(i)
        mx(i)=gamac(i)*momnxx(i)
        my(i)=by*gama(i)*momnyy(i)
        mz(i)=bz*gama(i)*momnzz(i)
        svm(i)=nx(i)*nx(i)+three*(ny(i)*ny(i)+nz(i)*nz(i)+mx(i)*mx(i))+
     1         my(i)*my(i)+mz(i)*mz(i)
        svm(i)=sqrt(svm(i)) 
      ENDDO
        fac2 = g/max(em20,young)
      DO  i=1,nel
C
C       TIME STEP
C
        dtemp    = one/max(xl(i),em20)
        ktran    = max(area,fac2*area_y(i),fac2*area_z(i))*dtemp
        krot     = four *max(iyy*dtemp,izz*dtemp)
        stifm(i) = youngm(i)*ktran
        stifr(i) = max( gm(i)*ixx*dtemp,youngm(i)*krot)
        viscm(i) = zero
        viscr(i) = zero
        mass(i)   = xl(i)*fac
        xiner(i)  = xl(i)*rho*max(ixx,imyz+area*xl2(i)/twelve)
      ENDDO
C
      nindx=0
      DO i=1,nel
        IF (svm(i) > yld(i) .AND. off(i) == one) THEN
          nindx=nindx+1
          index(nindx)=i
        ENDIF
      ENDDO
C
      IF (nindx /= 0) THEN
        DO j=1,nindx
          i=index(j)
          dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
        ENDDO
C
C WARNING, THE VERSION OF THE SPRING MODEL (Q.ZENG/A.COMBESCURE/C.BENOIT)
C NE CONSIDERE PAS D'ERREUR CUMULEE SUR PLUSIEURS ELEMENTS...
C
C --> EXISTE VERSION AVEC TEST MULTI-ELEMENTS (ERR_MULTI/RUSER29.F)
C
C
C
        n=0
        err =2*em4
        DO WHILE (err > em4 .AND. n < nmax)
          DO j=1,nindx
            i=index(j)
            dpla_i=dpla_j(i)
            IF (ifunc /= 0) THEN
              yld_i =yld(i)+h(i)*dpla_i
              ELSEIF (uvar(n0,i) > em20) THEN
              yld_i =yld(i)+h(i)*m*dpla_i*(uvar(n0,i)+dpla_i)**(m-1)
              ELSE
              yld_i =yld(i)
              ENDIF
            dre =youngm(i)*dpla_i/yld_i
            drg =g*dpla_i/yld_i
            pnx  =one/(one+dre)
            pny  =one/(one+three*drg*ay**2)
            pnz  =one/(one+three*drg*az**2)
            pmx  =one/(one+three*drg*gamac(i)**2)
            pmy  =one/(one+dre*(by*gama(i))**2)
            pmz  =one/(one+dre*(bz*gama(i))**2)    
            pnx2=pnx*pnx
            pny2=pny*pny
            pnz2=pnz*pnz
            pmx2=pmx*pmx
            pmy2=pmy*pmy
            pmz2=pmz*pmz
            nx2=nx(i)*nx(i)
            ny2=ny(i)*ny(i)
            nz2=nz(i)*nz(i)
            mx2=mx(i)*mx(i)
            my2=my(i)*my(i)
            mz2=mz(i)*mz(i)
            fn=nx2*pnx2+3.*(ny2*pny2+nz2*pnz2)
            fm=3.*mx2*pmx2+my2*pmy2+mz2*pmz2
            f = fn+fm-yld_i*yld_i
            dfe = nx2*pnx2*pnx+my2*pmy2*pmy*(by*gama(i))**2
     .                    +mz2*pmz2*pmz*(bz*gama(i))**2
            dfg = ny2*pny2*pny*ay**2+nz2*pnz2*pnz*az**2
     .                    +mx2*pmx2*pmx*gamac(i)**2
            IF (ifunc /= 0) THEN
              dsigy = -h(i)*yld_i*two
              ELSEIF (uvar(n0,i)>em20) THEN 
              dsigy = -(h(i)*yld_i*m*(dpla_i+uvar(n0,i))**(m-1))*two
            ELSE
              dsigy = zero
              ENDIF
            df =(-dfe*(youngm(i)-dre*h(i))/yld_i
     .           -nine*dfg*(g-drg*h(i))/yld_i)*two
            df=df+dsigy
            err=abs(f/df)
            dpla_j(i)=max(zero,dpla_i-f/df)
          ENDDO ! DO J=1,NINDX
          n=n+1
        ENDDO ! DO WHILE (ERR > EM4 .AND. N < NMAX)
C
C     CONTRAINTES PLASTIQUEMENT ADMISSIBLES
C
        DO j=1,nindx
          i=index(j)
          uvar(n0,i) = uvar(n0,i) + dpla_j(i)
          dpla_i=dpla_j(i)
C INCREMENT N+1
          IF (ifunc /= 0) THEN
            yld_i =yld(i)+h(i)*dpla_i
          ELSEIF (uvar(n0,i) > em20) THEN
            yld_i =yld(i)+h(i)*m*dpla_i*(uvar(n0,i))**(m-1)
          ELSE
            yld_i =yld(i)
          ENDIF
          c1    = dpla_i/yld_i
          dre =youngm(i)*c1
          drg =g*c1
          pnx  =one/(one+dre)
          pny  =one/(one+three*drg*ay**2)
          pnz  =one/(one+three*drg*az**2)
          pmx  =one/(one+three*drg*gamac(i)**2)
          pmy  =one/(one+dre*(by*gama(i))**2)
          pmz  =one/(one+dre*(bz*gama(i))**2)
C
          signxx(i) = signxx(i)*pnx
          signxy(i) = signxy(i)*pny
          signxz(i) = signxz(i)*pnz
          momnxx(i) = momnxx(i)*pmx
          momnyy(i) = momnyy(i)*pmy
          momnzz(i) = momnzz(i)*pmz
C
          uvar(n0+2,i) = uvar(n0+2,i)+three*c1*gamac(i)**2*abs(momnxx(i))    
          ueq =  c1*sqrt( signxx(i)**2 
     1                 +  gama(i)**4
     1                 *(by**4*momnyy(i)**2 + bz**4*momnzz(i)**2) )
          uvar(n0+3,i) = uvar(n0+3,i) + ueq      
C
C ENDOMAGEMENT
C
          IF (uvar(n0,i) > ps) THEN
             devol = dc*dpla_i/(pr - ps)
          ELSE
               devol = zero
          ENDIF
C
          uvar(n0+4,i) = uvar(n0+4,i)+devol
C
          IF (uvar(n0,i) >= pr .OR. uvar(n0+4,i) >= dc) THEN
            off(i) = zero
            uvar(n0+4,i) = dc
          ENDIF  
        ENDDO ! DO J=1,NINDX
      ENDIF ! IF (NINDX /= 0)
C
C  TRANSLATE FORCES AND MOMENTS :
C
      DO i=1,nel
        fx(i) = signxx(i)*area*off(i)
        fy(i) = signxy(i)*area_y(i)*off(i)
        fz(i) = signxz(i)*area_z(i)*off(i)
        xmom(i) = momnxx(i)*off(i)/facx
        ymom(i) = momnyy(i)*off(i)/facy
        zmom(i) = momnzz(i)*off(i)/facz
      ENDDO
C---
      RETURN
      END