Learning Forum

[hscj]

Hello everyone, I am a current master's student. I am currently attempting to develop the secondary application of the HJC constitutive model. The single-unit test has been successfully completed so far. However, in the model where the quality block impacts the thin plate, the damage and stress-strain updates of the custom HJC constitutive model are different from the evolution of the built-in HJC constitutive model of LS-DYNA. I have consulted many people and AI, and also reviewed numerous literature. I have found that the main problem might lie in the radial regression algorithm and the code for damage accumulation. Could everyone please help me identify which part of the code has the issue? Thank you very much!

 

HJC CODE BEGIN

      subroutine umat41 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,
     1 temper,failel,crv,cma)
c
c******************************************************************
c|  HJC Model Implementation with Plastic Volumetric Strain Damage |
c******************************************************************
c
      include 'nlqparm'
      include 'iounits.inc'
      common/bk06/idmmy,iaddp,ifil,maxsiz,ncycle,ctime(2,30)
      character*5 etype
      dimension cm(*),eps(*),sig(*),hsv(*)
      logical failel
      real epsp, dt1, capa, tt, temper
      
c     Variable definitions
      real G, K_elastic, A, B, C, N, FC, T, EPSILON_RATE, EFMIN, SFMAX   
      real PC, UC, UL, PL, D1, D2, K1, K2, K3
      real pre_sij(6), hydro_strain_Inc, eij_Inc(6), tr_sij(6)
      real tr_sij_J2, Q_trial, mu_old, mu_new, mu_max, mu_plastic
      real P_new, P_max, P_old_stress, P_star, T_star, T_current
      real eps_dot_star, strainRate_eq, Yield_star, Yield_HJC
      real eps_p_f, strainOld, Damage, e_vol_inc, d_mu_p
      real smallNum, k_value, true_sij(6), P_sum_check
      real K_unload, F_interp 
      real e_plaStrain_Inc, mu_bar, mu_plaStrain
      parameter (zero = 0.0d0, half = 0.5d0, one = 1.0d0,
     1           two = 2.0d0, three = 3.0d0)
     
      smallNum = 1.0d-20
      
c     Get material parameters
      G = cm(1)
      K_elastic = cm(2)
      A = cm(3)
      B = cm(4)
      C = cm(5)
      N = cm(6)
      FC = cm(7)
      T = cm(8)
      EPSILON_RATE = cm(9) 
      EFMIN = cm(10)
      SFMAX = cm(11)
      PC = cm(12)
      UC = cm(13)
      PL = cm(14)
      UL = cm(15)
      D1 = cm(16)
      D2 = cm(17)
      K1 = cm(18)
      K2 = cm(19)
      K3 = cm(20)

      if (etype.eq.'solid'.or.etype.eq.'shl_t') then

c     =====================================================================
c     Step 1: Extract history variables and preprocessing
c     =====================================================================
      P_old_stress    = hsv(1)  
      mu_old          = hsv(3)  
      strainOld       = hsv(5)  
      Damage          = hsv(14) 
      mu_max          = hsv(15) 
      P_max           = hsv(16) 
      mu_plaStrain    = hsv(2)
      
c     Calculate equivalent strain rate
      strainRate_eq = sqrt(two/three * (eps(1)**2 + eps(2)**2 +    
     1                eps(3)**2 + half*(eps(4)**2 + eps(5)**2 + 
     2                eps(6)**2))) / max(dt1, smallNum)
      
      eps_dot_star = strainRate_eq / max(EPSILON_RATE, smallNum)
      if (eps_dot_star .lt. one) eps_dot_star = one
      

c     Elastic predictor (deviatoric stress)
      pre_sij(1:3) = sig(1:3) + P_old_stress
      pre_sij(4:6) = sig(4:6)
      
      e_vol_inc = (eps(1) + eps(2) + eps(3)) / three
      eij_Inc(1:3) = eps(1:3) - e_vol_inc
      eij_Inc(4:6) = eps(4:6)
      
      tr_sij(1:3) = pre_sij(1:3) + two * G * eij_Inc(1:3)
      tr_sij(4:6) = pre_sij(4:6) + G * eij_Inc(4:6)
      
      tr_sij_J2 = half * (tr_sij(1)**2 + tr_sij(2)**2 + tr_sij(3)**2)
     1          + tr_sij(4)**2 + tr_sij(5)**2 + tr_sij(6)**2
      Q_trial = sqrt(three * tr_sij_J2)

c     =====================================================================
c     Step 2: Equation of state (EOS) and pressure parameter update
c     =====================================================================
      hydro_strain_Inc = - (eps(1) + eps(2) + eps(3))  
      mu_new = mu_old + hydro_strain_Inc                   
      
c     [Improvement] Precompute dynamic unloading and tensile bulk modulus K_unload 
c     Reference: Interpolate between K_elastic and K1 based on compaction level
      if (mu_max .le. UC) then
         K_unload = K_elastic
      else if (mu_max .lt. UL) then
         F_interp = (mu_max - UC) / (UL - UC)
         K_unload = (one - F_interp) * K_elastic + F_interp * K1
      else
         K_unload = K1 
      end if
      
      if (mu_new .ge. mu_max .and. mu_new .gt. zero) then 
        ! --- Loading path (skeleton curve) ---
        if (mu_new .le. UC) then
          P_new = K_elastic * mu_new
        else if (mu_new .le. UL) then
          P_new = PC + (mu_new - UC)*(PL - PC)/(UL - UC)
        else
          mu_bar = (mu_new - UL) / (one + UL)
          P_new = PL + K1 * mu_bar + K2 * mu_bar**2 + K3 * mu_bar**3
        end if
        mu_max = mu_new
        P_max  = P_new
      else 
        ! --- Unloading and tension ---
        if (mu_new .lt. zero) then
          ! Tension: use degraded unloading modulus according to literature
          P_new = K_unload * mu_new
        else 
          ! Compressive unloading: linear return from P_max using dynamic K_unload
          P_new = P_max - K_unload * (mu_max - mu_new)
          if (P_new .lt. zero) P_new = zero
        end if
      end if

c     Tension cutoff
      T_current = T * (one - Damage)
      if (P_new .lt. -T_current) P_new = -T_current

c     [Improvement] Calculate plastic volumetric strain increment d_mu_p 
c     Since modulus is dynamic K_unload, the elastic volumetric strain subtraction should also use current modulus
      d_mu_p = zero
      if (mu_new .gt. UC .and. hydro_strain_Inc .gt. zero) then
         d_mu_p = hydro_strain_Inc - (P_new - P_old_stress)/K_unload
         if (d_mu_p .lt. zero) d_mu_p = zero
      end if
      mu_plaStrain = mu_plaStrain + d_mu_p
c     Update P* and T* immediately (for denominator and HSV6 calculation)
      P_star = P_new / max(FC, smallNum)
      T_star = T / max(FC, smallNum)

c     =====================================================================
c     Step 3: Damage denominator calculation (NaN protection)
c     =====================================================================
      P_sum_check = P_star + T_star
      if (P_sum_check .le. 1.0d-8) then
        eps_p_f = EFMIN
      else
        eps_p_f = D1 * (P_sum_check)**D2
      end if
      if (eps_p_f .lt. EFMIN) eps_p_f = EFMIN  

c     =====================================================================
c     Step 4: Yield surface and radial return (fix multiplication order of SFMAX and strain rate)
c     Step 4: Failure fuse (advanced protection)
c     =====================================================================
c     If damage exceeds 0.999d0, trigger fuse early to avoid numerical noise
      if (Damage .gt. 0.999999d0) then
         Damage = 1.0d0             
         Yield_HJC = 0.001d0 * FC   ! residual strength
         e_plaStrain_Inc = zero
         k_value = Yield_HJC / max(Q_trial, smallNum)
         true_sij(1:6) = tr_sij(1:6) * k_value
          goto 999 
      endif

c     Yield surface calculation
         Yield_star = A * (one - Damage) + B * (P_star**N)* (one
     1      + C * log(eps_dot_star))
         if (Yield_star .gt. SFMAX) Yield_star = SFMAX

c        Add safety cutoff: prevent Yield_star from becoming negative due to numerical fluctuations
         if (Yield_star .lt. zero) Yield_star = zero
         
      
      Yield_HJC = Yield_star * FC

      if (Q_trial .gt. Yield_HJC) then
        e_plaStrain_Inc = (Q_trial - Yield_HJC) / (three * G)
        strainOld = strainOld + e_plaStrain_Inc
        k_value = Yield_HJC / max(Q_trial, smallNum)
        true_sij(1:6) = tr_sij(1:6) * k_value
      else
        e_plaStrain_Inc = zero
        true_sij(1:6) = tr_sij(1:6)
      end if

c     =====================================================================
c     Step 5: Damage accumulation (numerical tolerance completion)
c     =====================================================================
c     Only accumulate if damage has not reached the set 'quasi-failure threshold'
      if (Damage .lt. 0.999999d0) then
            Damage = Damage + (e_plaStrain_Inc + d_mu_p) / eps_p_f
      end if

c     [Core modification]: Force numerical closure
c     Once damage exceeds the threshold of 0.999999, directly set to 1.0
c     This solves the precision error concern and ensures complete red in contour
      if (Damage .ge. 0.999999d0) Damage = one
c     =====================================================================
c     Step 6: Update history variables (meeting user requirements: hsv2, hsv6)
c     =====================================================================
999   continue 
      sig(1:3) = true_sij(1:3) - P_new
      sig(4:6) = true_sij(4:6)
      
      hsv(1)  = P_new                        
      hsv(2)  = mu_plaStrain                
      hsv(3)  = mu_new                       
      hsv(4)  = e_plaStrain_Inc              
      hsv(5)  = strainOld                    
      hsv(6)  = D1 * (P_sum_check)**D2       
      hsv(7:12) = eps(1:6)
      hsv(13) = eps_dot_star          
      hsv(14) = Damage                       
      hsv(15) = mu_max                       
      hsv(16) = P_max                        
      epsp    = strainOld

      end if
      return
      end

HJC CODE END

 

LS-DYNA-HJC 

 

UMAT-HJC

 

[hscj]

My email address is hscj@mail.ustc.edu.cn. If you are interested in the second-level development, please feel free to contact me. I am also planning to explore new technologies such as near-field dynamics and neural network rendering!

[ErKo]

Hi
Troubleshooting and reviewing user defined material is outside what Ansys employees can help with. Please wait for other forum members that have experience with umats and user defined material to chime in and perhaps provide feedback. Thank you

Erko

[hscj]

Hi

Thank you for your reply. I am currently continuing to investigate the causes of the abnormal evolution of the damage. Thank you very much.

CJ

[ErKo]

 

Hi 

As we said troubleshooting and reviewing user defined material is outside what Ansys employees can help with.

Be patient and wait for other forum members (non Ansys empl.) that have experience with umats and user defined material to chime in perhaps and provide feedback. Thank you

Erko