fail_brokmann.F File Reference

#include "implicit_f.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	fail_brokmann (nel, nuparam, nuvar, time, timestep, uparam, ngl, signxx, signyy, signxy, uvar, off, ipt, nindxf, indxf, tdel)

Function/Subroutine Documentation

fail_brokmann()

subroutine fail_brokmann	(	integer	nel,
		integer	nuparam,
		integer	nuvar,
			time,
			timestep,
		intent(in)	uparam,
		integer, dimension(nel)	ngl,
		dimension(nel)	signxx,
		dimension(nel)	signyy,
		dimension(nel)	signxy,
		intent(inout)	uvar,
		dimension(nel)	off,
		integer	ipt,
		integer	nindxf,
		integer, dimension(nel)	indxf,
		dimension(nel)	tdel )

Definition at line 32 of file fail_brokmann.F.

!-----------------------------------------------
!   m o d u l e s
!-----------------------------------------------
      use newman_raju_mod
C-----------------------------------------------
c    Test of Ch. Brokmann criterion in glass, taking account for initial
c    statistically distributed micro cracks
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---------+---------+---+---+--------------------------------------------
      INTEGER NEL,NUPARAM,NUVAR,IPT,NINDXF
      my_real time,timestep
c
      INTEGER, DIMENSION(NEL) :: NGL,INDXF
      my_real, DIMENSION(NEL) :: signxx,signyy,signxy,off,tdel
      my_real, DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
      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,IDEB,NINDX
      INTEGER ,DIMENSION(NEL) :: INDX
c
      my_real :: alpha,alphai,exp_n,exp_m,k_ic,k_th,sig_ini,dsig_ini,xdamp,v0,
     .   da,dc,v_a,v_c,kcm,ktm,vm,sig_cos,sig_cosw,dsig_n,ya,yc,k1_a,k1_c,
     .   cr_ang,fac_m,fac_l,fac_t,fac_lenm,fac_mpa,fac_pi0,fac_pi2
c--------------------------------------------------------------------
c---  state variables for Ch.Brokmann extension
c     UVAR(15) = FAIL_B   :  Failure flag, set = 1 to allow Alter test
c     UVAR(16) = CR_LEN   :  Crack length
c     UVAR(17) = CR_DEPTH :  Crack width
c     UVAR(18) = CR_ANG   :  Random crack angle
c     UVAR(19) = THK0     :  Initial thickness of the shell (in micrometers)
c     UVAR(20) = ALDT0    :  Initial width of the shell (in micrometers)
c     UVAR(21) = SIG_COS  :  Crack opening stress (saved for filtering)
C=======================================================================
      exp_n   = uparam(1)
      k_ic    = uparam(6)
      k_th    = uparam(7)
      v0      = uparam(8)
      alpha   = uparam(10)
      ideb    = nint(uparam(17))
      sig_ini = uparam(30)
      fac_m   = uparam(33) 
      fac_l   = uparam(34)
      fac_t   = uparam(35)
c--------------------------------------------------
c     parameter initialization and unit_conversions
c--------------------------------------------------
      fac_lenm = ep06 * fac_l   ! conversion to micrometers
      fac_mpa  = em6  * fac_m / (fac_l * fac_t**2) ! stress conversion to MPa
      vm       = v0 * fac_l / fac_t  ! conversion to (m/s)
      kcm      = k_ic * sqrt(fac_l)  ! conversion to [stress*sqrt(m)]
      ktm      = k_th * sqrt(fac_l)  ! conversion to [stress*sqrt(m)]
c-----------------------------------
      alphai = one - alpha
      nindx  = 0
c
      dsig_ini = 2.0    ! [MPa/s]
      xdamp    = 25.0   ! hard coded damping exponent
      exp_m    = one / (one + xdamp) 
      fac_pi0  = zero
      fac_pi2  = half
c--------------------------------------------------------------------
      DO i = 1,nel
        IF (off(i) == one .and. uvar(i,15) == zero .and. uvar(i,16) > zero) THEN
          v_a  = zero
          v_c  = zero          
c
c         Crack openenig stress (stress rotated in crack direction)
          cr_ang  = uvar(i,18) * two
          sig_cos = half*(signxx(i)+signyy(i)) 
     .            + cos(cr_ang) * half*(signxx(i)-signyy(i))
     .            + sin(cr_ang) * half*signxy(i)
          sig_cos = sig_cos * alpha + uvar(i,21)*alphai   ! filter using exp average
c
c         take into account high stress rates
          dsig_n = (sig_cos - uvar(i,21)) / max(timestep,em20)
          dsig_n = dsig_n * fac_mpa / fac_t  ! => (MPa/s)
 
          uvar(i,21) = sig_cos ! Save Crack Opening Stress (COS)
          IF (dsig_n > dsig_ini) THEN
            sig_cos = sig_cos * (dsig_ini/abs(dsig_n))**exp_m
          END IF          
c------------------------
c         Residual stress
c------------------------
          sig_cosw = sig_cos - sig_ini   ! Subtracting initial surface stress
          sig_cosw = max(zero , sig_cosw)
c------------------------
c         Test criteria
c------------------------
          CALL newman_raju(uvar(i,17),uvar(i,16),uvar(i,19),uvar(i,20),fac_pi2,ya)
          CALL newman_raju(uvar(i,17),uvar(i,16),uvar(i,19),uvar(i,20),fac_pi0,yc)     
          k1_a = ya * sig_cosw * sqrt(pi*uvar(i,16)*em6)  
          k1_c = yc * sig_cosw * sqrt(pi*uvar(i,17)*em6)
c
c         Check rupture criterion #1
          IF (k1_a >= kcm) THEN
            tdel(i) = time            ! tag to start damage                                                          
            nindx   = nindx + 1                                                           
            nindxf  = nindxf+ 1                                                    
            indx(nindx)   = i                                                           
            indxf(nindxf) = i         ! tag to save crack direction                                                          
            uvar(i,15)    = one        ! set failure flag = 1
          END IF
c         
          IF (k1_a >= ktm .and. k1_a < kcm) v_a = v0*(k1_a/kcm)**exp_n
          IF (k1_c >= ktm .and. k1_a < kcm) v_c = v0*(k1_c/kcm)**exp_n
c
c         Crack growth
          da = v_a * timestep * fac_lenm
          dc = v_c * timestep * fac_lenm
          uvar(i,16) = uvar(i,16) + da
          uvar(i,17) = uvar(i,17) + dc
c
c         Check rupture criterion #2
          k1_a = ya * sig_cos * sqrt(pi*uvar(i,16)*em6)
          IF (k1_a >= kcm .and. uvar(i,15) == zero) THEN
            tdel(i) = time            ! tag to start damage                                                          
            nindx   = nindx + 1                                                           
            nindxf  = nindxf+ 1                                                    
            indx(nindx)   = i                                                           
            indxf(nindxf) = i         ! tag to save crack direction                                                          
            uvar(i,15)    = one        ! set failure flag = 1
          END IF
    
        END IF
      END DO
c-----------------------------------------------
      IF (nindx > 0) THEN
        DO j=1,nindx
          i = indx(j)
#include "lockon.inc"
          WRITE(iout, 1000) ngl(i),ipt,time                                
          WRITE(istdo,1100) ngl(i),ipt,time                               
#include "lockoff.inc"
        ENDDO        
      ENDIF          
c-----------------------------------------------
 1000 FORMAT(1x, ' **** FAILURE due to Brokmann criterion : ELEMENT ',i10,
     .       ' INT POINT',i2,' AT TIME ',1pe12.4)
 1100 FORMAT(1x, ' **** FAILURE due to Brokmann criterion : ELEMENT ',i10,
     .       ' INT POINT',i2,' AT TIME ',1pe12.4)
c-----------
      RETURN