neo_hook_t.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/.
C   Hyperelastic model: THERMAL DEPENDANT NEO-HOOK
C   compute stresses using MATB which is the 
C   left cauchy green tensor [b]=[F][FT]
!||====================================================================
!||    neo_hook_t   ../engine/source/materials/mat/mat100/neo_hook_t.F
!||--- called by ------------------------------------------------------
!||    sigeps100    ../engine/source/materials/mat/mat100/sigeps100.F90
!||--- calls      -----------------------------------------------------
!||    mullins_or   ../engine/source/materials/fail/mullins_or/fail_mullins_OR_s.F
!||====================================================================
       SUBROUTINE neo_hook_t(
     1                NEL , MATB,SIG,
     4                BI1,JDET ,FLAG_MUL,MU,D,
     5                NVARF,COEFR, BETAF,COEFM  ,UVARF)
C-----------------------------------------------
C   I M P L I C I T   T Y P E S
C-----------------------------------------------
#include "implicit_f.inc"
C----------------------------------------------------------------
C  I N P U T   A R G U M E N T S
C----------------------------------------------------------------
      INTEGER,  INTENT(IN) ::       NEL,FLAG_MUL,NVARF
      my_real,  INTENT(IN) ::  COEFR, BETAF,COEFM        
 
      my_real, DIMENSION(NEL),  INTENT(IN) ::  mu,d     
      my_real, DIMENSION(NEL, 3,3), INTENT(IN) :: matb(nel,3,3)
C      my_real, DIMENSION(NEL), INTENT(IN) :: MU
C----------------------------------------------------------------
C  O U T P U T   A R G U M E N T S
C----------------------------------------------------------------
      my_real, DIMENSION(NEL), INTENT(OUT) ::  bi1,jdet 
      my_real, DIMENSION(NEL, 3,3), INTENT(OUT)   :: sig
C----------------------------------------------------------------
C  I N P U T  O U T P U T   A R G U M E N T S
C----------------------------------------------------------------
      my_real, DIMENSION(NEL,NVARF), INTENT(INOUT) :: uvarf
C----------------------------------------------------------------
C  L O C A L  V A R I B L E S
C----------------------------------------------------------------
      INTEGER    I
 
      my_real     
     .   lpchain(nel), trace(nel),traceb(nel),j2third(nel),
     .   sb1(nel), sb2(nel),sb3(nel),tbnorm(nel),dgamma(nel),i1(nel),
     .   aa,bb,cc,trb2,trb22,invj(nel), 
     .   jthird(nel),j4third(nel),eta(nel),ww(nel) ,
     .   dphidi1(nel)  , dphidj(nel) 
C----------------------------------------------------------------
      eta(1:nel) = one !MULLINS DAMAGE FACTOR
 
      DO i = 1,nel
         !J = RHO0/RHO = RELATIVE VOLUME = DET F (F = GRADIENT OF DEFORMATION)
         jdet(i)=matb(i,1,1)*matb(i,2,2)*matb(i,3,3) -matb(i,1,1)*matb(i,2,3)*matb(i,3,2) -
     .       matb(i,3,3)*matb(i,1,2)*matb(i,2,1) +matb(i,1,2)*matb(i,2,3)*matb(i,3,1) +
     .       matb(i,2,1)*matb(i,3,2)*matb(i,1,3) -matb(i,2,2)*matb(i,3,1)*matb(i,1,3)
         jdet(i)= sqrt(max(em20, jdet(i)))    
 
         !FIRST INVARIANT 
         i1(i) = matb(i,1,1)+matb(i,2,2)+matb(i,3,3)
 
         IF(jdet(i)>zero) THEN
          j2third(i)  = exp((-two_third )*log(jdet(i)))
         ELSE
          j2third(i)  = zero
         ENDIF
      ENDDO
      DO i = 1,nel
        !INVARIANT BAR
        !first invariant deviator BI1:
         bi1(i) = i1(i) * j2third(i) 
      ENDDO
 
 
      !====================================
      !dammge by mullins effect
      !====================================
      IF(flag_mul == 1)THEN
       DO i = 1,nel
         ww(i) =  mu(i) *(bi1(i)-three) /two   
 
       ENDDO
       CALL mullins_or(
     1                nel ,nvarf,  coefr,betaf ,
     2                coefm, ww  , uvarf,eta )
      ENDIF
      !====================================
 
      DO i = 1,nel
        eta(i) = max(min(eta(i),one),em20)
        dphidi1(i) = eta(i)*mu(i)/two
        dphidj(i)  = d(i)* (jdet(i)-one) 
        invj(i)=one/max(em20,jdet(i))
      ENDDO
 
      DO i=1,nel   
         !CAUCHY STRESS
         aa =  eta(i)*mu(i)/max(em20,jdet(i))
 
         sig(i,1,1) =  aa*(matb(i,1,1)-third*i1(i))         
     .       +   dphidj(i)        
         sig(i,2,2) =  aa*(matb(i,2,2)-third*i1(i))               
     .       +   dphidj(i)        
         sig(i,3,3) =  aa*(matb(i,3,3)-third*i1(i))            
     .       +   dphidj(i)        
         sig(i,1,2) =  aa*matb(i,1,2)
         sig(i,2,3) =  aa*matb(i,2,3)
         sig(i,3,1) =  aa*matb(i,3,1)
         sig(i,2,1)=sig(i,1,2)   
         sig(i,3,2)=sig(i,2,3)   
         sig(i,1,3)=sig(i,3,1)   
      ENDDO
C   
      RETURN
      END