fail_orthbiquad_s.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/.
!||====================================================================
!||    fail_orthbiquad_s   ../engine/source/materials/fail/orthbiquad/fail_orthbiquad_s.F
!||--- called by ------------------------------------------------------
!||    mmain               ../engine/source/materials/mat_share/mmain.f90
!||    mmain8              ../engine/source/materials/mat_share/mmain8.F
!||    mulaw               ../engine/source/materials/mat_share/mulaw.F90
!||    mulaw8              ../engine/source/materials/mat_share/mulaw8.F90
!||    usermat_solid       ../engine/source/materials/mat_share/usermat_solid.F
!||--- calls      -----------------------------------------------------
!||    finter              ../engine/source/tools/curve/finter.F
!||====================================================================
      SUBROUTINE fail_orthbiquad_s (
     1     NEL      ,NUPARAM  ,NUVAR    ,MFUNC    ,KFUNC    ,
     2     NPF      ,TF       ,TIME     ,TIMESTEP ,UPARAM   ,
     3     NGL      ,DPLA     ,EPSP     ,UVAR     ,OFF      ,
     4     SIGNXX   ,SIGNYY   ,SIGNZZ   ,SIGNXY   ,SIGNYZ   ,SIGNZX   ,
     5     DFMAX    ,TDELE    ,ALDT     )
C!-----------------------------------------------
C!   I m p l i c i t   T y p e s
C!-----------------------------------------------
#include      "implicit_f.inc"
C!---------+---------+---+---+--------------------------------------------
C! VAR     | SIZE    |TYP| RW| DEFINITION
C!---------+---------+---+---+--------------------------------------------
C! NEL     |  1      | I | R | SIZE OF THE ELEMENT GROUP NEL 
C! NUPARAM |  1      | I | R | SIZE OF THE USER PARAMETER ARRAY
C! NUVAR   |  1      | I | R | NUMBER OF FAILURE ELEMENT VARIABLES
C!---------+---------+---+---+--------------------------------------------
C! MFUNC   |  1      | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
C! KFUNC   | NFUNC   | I | R | FUNCTION INDEX not used
C! NPF     |  *      | I | R | FUNCTION ARRAY   
C! TF      |  *      | F | R | FUNCTION ARRAY 
C!---------+---------+---+---+--------------------------------------------
C! TIME    |  1      | F | R | CURRENT TIME
C! TIMESTEP|  1      | F | R | CURRENT TIME STEP
C! UPARAM  | NUPARAM | F | R | USER FAILURE PARAMETER ARRAY
C!---------+---------+---+---+--------------------------------------------
C! SIGNXX  | NEL     | F | W | NEW ELASTO PLASTIC STRESS XX
C! SIGNYY  | NEL     | F | W | NEW ELASTO PLASTIC STRESS YY
C! ...     |         |   |   |
C! ...     |         |   |   |
C!---------+---------+---+---+--------------------------------------------
C! UVAR    |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
C! OFF     | NEL     | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
C!---------+---------+---+---+--------------------------------------------
#include "mvsiz_p.inc"
#include "scr17_c.inc"
#include "units_c.inc"
#include "comlock.inc"
C!-----------------------------------------------
      INTEGER NEL, NUPARAM, NUVAR,NGL(NEL)
      my_real TIME,TIMESTEP,UPARAM(*),
     .   SIGNXX(NEL),SIGNYY(NEL),SIGNZZ(NEL),
     .   SIGNXY(NEL),SIGNYZ(NEL),SIGNZX(NEL),UVAR(NEL,NUVAR),
     .   DPLA(NEL),EPSP(NEL),OFF(NEL),DFMAX(NEL),TDELE(NEL),
     .   aldt(nel)     
C!-----------------------------------------------
C!   VARIABLES FOR FUNCTION INTERPOLATION 
C!-----------------------------------------------
      INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
      my_real FINTER ,TF(*)
      EXTERNAL FINTER
C!        Y = FINTER(IFUNC(J),X,NPF,TF,DF)
C!        Y       : y = f(x)
C!        X       : x
C!        DF      : f'(x) = dy/dx
C!        IFUNC(J): FUNCTION INDEX
C!              J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
C!        NPF,TF  : FUNCTION PARAMETER
C!-----------------------------------------------
C!   L o c a l   V a r i a b l e s
C!-----------------------------------------------
      INTEGER I,K,J,INDX(MVSIZ),NINDX,SEL,NANGLE,IREG,IRATE
      my_real 
     .        DF,P,triaxs,SVM,SXX,SYY,SZZ,EPS_FAIL,
     .        P1X,P1Y,S1X,S1Y,S2Y,A1,B1,C1,REF_EL_LEN,LAMBDA,FAC,
     .        X_1(3),X_2(3),COS2(10,10),EPSD0,CJC,RATE_SCALE,FRATE,
     .        mohr_radius,cos2theta(nel)
      my_real sig1(nel),sig2(nel),mohr_center
      my_real, DIMENSION(:), ALLOCATABLE :: q_x11,q_x12,q_x13,q_x21,q_x22,q_x23,q_inst
C
      DATA  cos2/
     1 1.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,   
     2 0.   ,1.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     3 -1.  ,0.   ,2.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     4 0.   ,-3.  ,0.   ,4.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     5 1.   ,0.   ,-8.  ,0.   ,8.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     6 0.   ,5.   ,0.   ,-20. ,0.   ,16.  ,0.   ,0.   ,0.   ,0.   ,     
     7 -1.  ,0.   ,18.  ,0.   ,-48. ,0.   ,32.  ,0.   ,0.   ,0.   ,     
     8 0.   ,-7.  ,0.   ,56.  ,0.   ,-112.,0.   ,64.  ,0.   ,0.   ,     
     9 1.   ,0.   ,-32. ,0.   ,160. ,0.   ,-256.,0.   ,128. ,0.   ,     
     a 0.   ,9.   ,0.   ,-120.,0.   ,432. ,0.   ,-576 ,0.   ,256. /
C!--------------------------------------------------------------
      !=======================================================================
      ! - INITIALISATION OF COMPUTATION ON TIME STEP
      !=======================================================================
      ! Recovering model parameters
      sel        = int(uparam(2))
      ref_el_len = uparam(3)
      epsd0      = uparam(4)
      cjc        = uparam(5)
      rate_scale = uparam(6)
      nangle     = int(uparam(7))
c
      ! -> Allocation of factors
      ALLOCATE(q_x11(nangle),q_x12(nangle),q_x13(nangle),
     .         q_x21(nangle),q_x22(nangle),q_x23(nangle))
      ! Recovering factors
      IF (sel == 3) THEN
        ALLOCATE(q_inst(nangle))        
        DO j = 1,nangle
          q_x11(j)  = uparam(8  + 7*(j-1))
          q_x12(j)  = uparam(9  + 7*(j-1))
          q_x13(j)  = uparam(10 + 7*(j-1)) 
          q_x21(j)  = uparam(11 + 7*(j-1))
          q_x22(j)  = uparam(12 + 7*(j-1))
          q_x23(j)  = uparam(13 + 7*(j-1))
          q_inst(j) = uparam(14 + 7*(j-1))
        ENDDO 
      ELSE
        DO j = 1,nangle
          q_x11(j)  = uparam(8  + 6*(j-1))
          q_x12(j)  = uparam(9  + 6*(j-1))
          q_x13(j)  = uparam(10 + 6*(j-1)) 
          q_x21(j)  = uparam(11 + 6*(j-1))
          q_x22(j)  = uparam(12 + 6*(j-1))
          q_x23(j)  = uparam(13 + 6*(j-1)) 
        ENDDO
      ENDIF
      IF (sel == 3) sel = 2
c
      ! Recovering functions
      irate = 0
      IF (rate_scale /= zero) THEN
        irate = 1
      ENDIF
      ireg  = 0
      IF (ref_el_len /= zero) THEN
        ireg = irate + 1
      ENDIF       
c      
      ! At initial time, compute the element size regularization factor
      IF (uvar(1,3) == zero .AND. ireg > 0) THEN
        DO i=1,nel   
          lambda    = aldt(i)/ref_el_len 
          fac       = finter(kfunc(ireg),lambda,npf,tf,df) 
          uvar(i,3) = fac
        ENDDO
      ENDIF
c
      ! Checking element failure and recovering user variable
      DO i=1,nel
       IF (off(i) < em01) off(i) = zero
       IF (off(i) < one .AND. off(i) > zero) off(i) = off(i)*four_over_5
      END DO
C
      ! Initialization of variable
      nindx = 0  
c      
      !====================================================================
      ! - LOOP OVER THE ELEMENT TO COMPUTE THE DAMAGE VARIABLE
      !====================================================================       
      DO i=1,nel
c
        ! Failure strain initialization
        eps_fail = zero
c
        ! If the element is not broken
        IF (off(i) ==  one .AND. dpla(i) /= zero) THEN
c
          ! Computation of loading angle
          mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
          mohr_center = (signxx(i)+signyy(i))/two
          IF (mohr_radius > em20) THEN
            cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
          ELSE
            cos2theta(i) = one
          ENDIF
          sig1(i) = mohr_center + mohr_radius
          sig2(i) = mohr_center - mohr_radius
          IF (sig1(i)<zero.OR.((sig2(i)<zero).AND.(sig2(i)<-sig1(i)))) THEN 
            cos2theta(i) = -cos2theta(i)
          ENDIF
c
          ! Computation of the BIQUAD parameters
          x_1(1:3) = zero
          x_2(1:3) = zero
          DO j = 1,nangle
            DO k = 1,j              
              x_1(1) = x_1(1) + q_x11(j)*cos2(k,j)*(cos2theta(i))**(k-1)
              x_1(2) = x_1(2) + q_x12(j)*cos2(k,j)*(cos2theta(i))**(k-1)
              x_1(3) = x_1(3) + q_x13(j)*cos2(k,j)*(cos2theta(i))**(k-1)
              x_2(1) = x_2(1) + q_x21(j)*cos2(k,j)*(cos2theta(i))**(k-1)
              x_2(2) = x_2(2) + q_x22(j)*cos2(k,j)*(cos2theta(i))**(k-1)
              x_2(3) = x_2(3) + q_x23(j)*cos2(k,j)*(cos2theta(i))**(k-1)
            ENDDO          
          ENDDO
c        
          ! Computation of hydrostatic stress, Von Mises stress, and stress triaxiality
          p   = third*(signxx(i) + signyy(i) + signzz(i))
          sxx = signxx(i) - p
          syy = signyy(i) - p
          szz = signzz(i) - p
          svm = half*(sxx**2 + syy**2 + szz**2)
     .          + signxy(i)**2 + signzx(i)**2 + signyz(i)**2
          svm = sqrt(three*svm)
          triaxs = p/max(em20,svm)
          IF (triaxs < -two_third) triaxs = -two_third
          IF (triaxs >  two_third) triaxs =  two_third
c
          ! Computing the plastic strain at failure according to stress triaxiality
          IF (triaxs <= third) THEN
            eps_fail = x_1(1) + x_1(2)*triaxs + x_1(3)*triaxs**2
            IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
          ELSE
            SELECT CASE (sel)
              CASE(1)
                eps_fail = x_2(1) + x_2(2)*triaxs + x_2(3)*triaxs**2
                IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
              CASE(2)
                IF (triaxs <= one/sqr3) THEN                     ! triax < 0.57735
                  p1x      = third
                  p1y      = x_1(1) + x_1(2)*p1x + x_1(3)*p1x**2
                  s1x      = one/sqr3
                  s1y      = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
                  a1       = (p1y - s1y)/(p1x - s1x)**2
                  b1       = -two*a1*s1x
                  c1       = a1*s1x**2 + s1y 
                  eps_fail = c1 + b1*triaxs + a1*triaxs**2
                  IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
                ELSE                                             ! triax > 0.57735
                  p1x      = two*third
                  p1y      = x_2(1) + x_2(2)*p1x + x_2(3)*p1x**2
                  s1x      = one/sqr3
                  s1y      = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
                  a1       = (p1y - s1y)/(p1x - s1x)**2
                  b1       = -two*a1*s1x
                  c1       = a1*s1x**2 + s1y 
                  eps_fail = c1 + b1*triaxs + a1*triaxs**2
                  IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
                ENDIF
            END SELECT
          ENDIF
c
          ! Strain-rate effects
          IF (cjc /= zero .OR. irate /= 0) THEN 
            IF (cjc /= zero) THEN 
              frate = one
              IF (epsp(i) > epsd0) frate = frate + cjc*log(epsp(i)/epsd0)
            ELSEIF (irate /= 0) THEN 
              frate = finter(kfunc(irate),epsp(i)/rate_scale,npf,tf,df) 
            ENDIF
            eps_fail = eps_fail*frate
          ENDIF          
c
          ! Computation of damage variables
          dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
          dfmax(i) = min(one,dfmax(i))
c
          ! Checking element failure using global damage
          IF (dfmax(i) >= one .AND. off(i) == one) THEN
            off(i)      = four_over_5
            nindx       = nindx + 1
            indx(nindx) = i
            tdele(i)    = time    
          ENDIF
        ENDIF 
      ENDDO
c------------------------
c------------------------
      ! Deallocation of table
      IF (ALLOCATED(q_x11))  DEALLOCATE(q_x11)
      IF (ALLOCATED(q_x12))  DEALLOCATE(q_x12)
      IF (ALLOCATED(q_x13))  DEALLOCATE(q_x13)
      IF (ALLOCATED(q_x21))  DEALLOCATE(q_x21)
      IF (ALLOCATED(q_x22))  DEALLOCATE(q_x22)
      IF (ALLOCATED(q_x23))  DEALLOCATE(q_x23)
      IF (ALLOCATED(q_inst)) DEALLOCATE(q_inst)
c------------------------
c------------------------
      IF (nindx > 0) THEN
        DO j=1,nindx
          i = indx(j)     
#include "lockon.inc"
          WRITE(iout, 1000) ngl(i),time
          WRITE(istdo,1100) ngl(i),time
#include "lockoff.inc"
        END DO
      END IF         
c------------------------
 1000 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (ORTHBIQUAD) el#',i10,
     .          ' AT TIME :',1pe12.4)     
 1100 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (ORTHBIQUAD) el#',i10,
     .          ' AT TIME :',1pe12.4)
      END