fail_wilkins_c.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_wilkins_c   ../engine/source/materials/fail/wilkins/fail_wilkins_c.F
!||--- called by ------------------------------------------------------
!||    mulawc           ../engine/source/materials/mat_share/mulawc.F90
!||    usermat_shell    ../engine/source/materials/mat_share/usermat_shell.F
!||====================================================================
      SUBROUTINE fail_wilkins_c(
     1     NEL      ,NUPARAM  ,NUVAR    ,UPARAM   ,UVAR     ,
     2     TIME     ,NGL      ,IPG      ,ILAY     ,IPT      ,
     3     SIGNXX   ,SIGNYY   ,SIGNXY   ,SIGNYZ   ,SIGNZX   ,
     4     DPLA     ,FOFF     ,DFMAX    ,TDEL     )
C-----------------------------------------------
C     2D Wilkins Failure model
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include "units_c.inc"
#include "comlock.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 UPARAM  | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
C NUVAR   |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES
C UVAR    |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
C---------+---------+---+---+--------------------------------------------
C TIME    |  1      | F | R | CURRENT TIME
C TIMESTEP|  1      | F | R | CURRENT TIME STEP
C---------+---------+---+---+--------------------------------------------
C SIGNXX  | NEL     | F | R | NEW ELASTO PLASTIC STRESS XX
C SIGNYY  | NEL     | F | R | NEW ELASTO PLASTIC STRESS YY
C ...     |         |   |   |
C DPLA    | NEL     | F |R/W| PLASTIC STRAIN
C---------+---------+---+---+--------------------------------------------
C OFF     | NEL     | F | R | DELETED ELEMENT FLAG (=1. ON, =0. OFF)
C FOFF    | NEL     | I |R/W| DELETED INTEGRATION POINT FLAG (=1 ON, =0 OFF)
C DFMAX   | NEL     | F |R/W| MAX DAMAGE FACTOR 
C TDEL    | NEL     | F | W | FAILURE TIME
C---------+---------+---+---+--------------------------------------------
C NGL                         ELEMENT ID
C IPG                         CURRENT GAUSS POINT (in plane)
C ILAY                        CURRENT LAYER
C IPT                         CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
C---------+---------+---+---+--------------------------------------------
C   I N P U T   A r g u m e n t s
C-----------------------------------------------
      INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,ILAY,IPT
      INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
      my_real ,INTENT(IN) :: TIME
      my_real ,DIMENSION(NEL) ,INTENT(IN) :: DPLA,
     .   signxx,signyy,signxy,signyz,signzx
      my_real,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
C-----------------------------------------------
C   I N P U T   O U T P U T   A r g u m e n t s 
C-----------------------------------------------
      INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: FOFF
      my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: DFMAX
      my_real ,DIMENSION(NEL) ,INTENT(OUT)   :: TDEL
      my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: UVAR
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,J,NINDX
      INTEGER ,DIMENSION(NEL) :: INDX
      my_real :: AL,BE,PC,DC,P,W1,W2,SS,S1,S2,S3,SS0,SS1,SS2,SIG1,SIG2,
     .           A,A1,A2
C=======================================================================
      al = uparam(1)
      be = uparam(2)
      pc = uparam(3)
      dc = uparam(4)
      nindx  = 0  
C-----------------------------------------------
      DO i=1,nel
        IF (foff(i) == 1) THEN                               
          ss1 = half*(signxx(i) + signyy(i))  
          ss2 = half*(signxx(i) - signyy(i))  
          ss0 = sqrt(ss2**2 + signxy(i)**2)     
          sig1 = ss1  + ss0                       
          sig2 = ss1  - ss0                       
          ss = sig1                               
          IF (sig2 > sig1) THEN                     
            sig1 = sig2                              
            sig2  = ss                             
          ENDIF
          p   = third*(signxx(i) + signyy(i))
          ss1 = ss1 - p
          s1  = ss1  + ss0
          s2  = ss1  - ss0
          s3  = -p
          IF (sig2 > zero) THEN
            a1 = s2 / s1
            a2 = s2 / s3
          ELSEIF (sig1 < zero) THEN
            a1 = s1 / s3
            a2 = s1 / s2
          ELSEIF (sig1 /= zero) THEN
            a1 = s3 / s1
            a2 = s3 / s2
          ELSE
            a1 = one
            a2 = one
          ENDIF
          a = max(a1,a2)
          w1 = max(em20,(two - a))**be
          w2 = one - p/pc
          w2 = (max(em20,one/w2))**al
          uvar(i,1) = uvar(i,1) + w1*w2*dpla(i)             
          IF (uvar(i,1) >= dc) THEN
            nindx = nindx + 1
            indx(nindx) = i
            foff(i) = 0
            tdel(i) = time
          ENDIF          
        ENDIF 
      ENDDO      
c
c---  Maximum Damage storing for output : 0 < DFMAX < 1
      DO i=1,nel
        dfmax(i) = min(one,max(dfmax(i),uvar(i,1)/dc))
      ENDDO      
c------------------------
      IF (nindx > 0) THEN
        DO j=1,nindx
          i = indx(j)
#include  "lockon.inc"
          WRITE(iout, 2000) ngl(i),ipg,ilay,ipt
          WRITE(istdo,2100) ngl(i),ipg,ilay,ipt,time
#include  "lockoff.inc"
        END DO
      END IF              
c------------------------
 2000 FORMAT(1x,'FAILURE (WILKINS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
     .       i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3)
 2100 FORMAT(1x,'FAILURE (WILKINS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
     .       i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,1x,'AT TIME :',1pe12.4)
c------------------------
      RETURN
      END