sigaboyce.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: ARRUDA BOYCE strain energy potential
C   compute stresses using MATB which is the 
C   left cauchy green tensor [b]=[F][FT]
!||====================================================================
!||    sigaboyce    ../engine/source/materials/mat/mat100/sigaboyce.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 sigaboyce(
     1                NEL , MATB ,C1,C2,C3,
     3                C4  ,C5   ,MU,BETA,D ,
     4                SIG ,BI1  ,JDET ,FLAG_MUL,
     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) :: C1,C2,C3,C4,C5,MU,BETA,D ,
     .                        COEFR, BETAF,COEFM          
      my_real, DIMENSION(NEL, 3,3) , INTENT(IN) ::  MATB 
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)  :: 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     
     .   J2THIRD(NEL),I1(NEL),JTHIRD(NEL),WW (NEL),
     .   dphidi1(nel) ,dphidj(nel) ,eta(nel)
C----------------------------------------------------------------
      eta(1:nel) = one !MULLINS DAMAGE FACTOR
      DO i = 1,nel
         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)))    
        !INVARIANT BAR         
         i1(i) = matb(i,1,1)+matb(i,2,2)+matb(i,3,3)
         IF(jdet(i)>zero) THEN
          jthird(i)  = exp((-third )*log(jdet(i)))
          j2third(i) = jthird(i)**2       
         ELSE
          jthird(i)  = zero
          j2third(i) = zero
         ENDIF
        !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 *( c1 * (bi1(i)- three)+c2 * beta *   (bi1(i)**2- nine)
     .   +                                   c3 * beta**2 *(bi1(i)**3- three**3)    
     .   +                                   c4 * beta**3 *(bi1(i)**4- three**4)
     .   +                                   c5 * beta**4 *(bi1(i)**5- three**5) )
 
       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)* two*mu *( c1 + two * c2 * beta * bi1(i)
     .   +                           three* c3 *(beta *bi1(i))**2    
     .   +                          four* c4 *(beta *bi1(i))**3 
     .   +                          five  * c5 *(beta *bi1(i))**4)/max(em20,jdet(i))
 
        dphidj(i)  = d * (  jdet(i) - one /max(em20,jdet(i)) ) !d=1/d
      ENDDO
      DO i = 1,nel
         sig(i,1,1) =  dphidi1(i)* (matb(i,1,1)-third*bi1(i))
     .       +   dphidj(i)        
         sig(i,2,2) =  dphidi1(i)* (matb(i,2,2)-third*bi1(i))
     .       +   dphidj(i)        
         sig(i,3,3) =  dphidi1(i)* (matb(i,3,3)-third*bi1(i))
     .       +   dphidj(i)        
         sig(i,1,2) =  dphidi1(i)*matb(i,1,2)
         sig(i,2,3) =  dphidi1(i)*matb(i,2,3)
         sig(i,3,1) =  dphidi1(i)*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
 
      RETURN
      END
C   Hyperelastic model: ARRUDA BOYCE strain  energy potential