{"id":457660,"date":"2026-05-06T09:09:33","date_gmt":"2026-05-06T09:09:33","guid":{"rendered":"https:\/\/innovationspace.ansys.com\/forum\/forums\/topic\/the-damage-evolution-law-in-the-holmquist-johnson-cook-hjc-model\/"},"modified":"2026-05-06T09:09:33","modified_gmt":"2026-05-06T09:09:33","slug":"the-damage-evolution-law-in-the-holmquist-johnson-cook-hjc-model","status":"publish","type":"topic","link":"https:\/\/innovationspace.ansys.com\/forum\/forums\/topic\/the-damage-evolution-law-in-the-holmquist-johnson-cook-hjc-model\/","title":{"rendered":"The damage evolution law in the Holmquist-Johnson-Cook (HJC) model"},"content":{"rendered":"<p>&lt;p&gt;Hello everyone, I am a current master&#8217;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!&lt;\/p&gt;&lt;p&gt;&nbsp;&lt;\/p&gt;&lt;p&gt;HJC CODE BEGIN&lt;\/p&gt;&lt;p&gt;&nbsp; &nbsp; &nbsp; subroutine umat41 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,&lt;br&gt;&nbsp; &nbsp; &nbsp;1 temper,failel,crv,cma)&lt;br&gt;c&lt;br&gt;c******************************************************************&lt;br&gt;c| &nbsp;HJC Model Implementation with Plastic Volumetric Strain Damage |&lt;br&gt;c******************************************************************&lt;br&gt;c&lt;br&gt;&nbsp; &nbsp; &nbsp; include &#8216;nlqparm'&lt;br&gt;&nbsp; &nbsp; &nbsp; include &#8216;iounits.inc'&lt;br&gt;&nbsp; &nbsp; &nbsp; common\/bk06\/idmmy,iaddp,ifil,maxsiz,ncycle,ctime(2,30)&lt;br&gt;&nbsp; &nbsp; &nbsp; character*5 etype&lt;br&gt;&nbsp; &nbsp; &nbsp; dimension cm(*),eps(*),sig(*),hsv(*)&lt;br&gt;&nbsp; &nbsp; &nbsp; logical failel&lt;br&gt;&nbsp; &nbsp; &nbsp; real epsp, dt1, capa, tt, temper&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;c &nbsp; &nbsp; Variable definitions&lt;br&gt;&nbsp; &nbsp; &nbsp; real G, K_elastic, A, B, C, N, FC, T, EPSILON_RATE, EFMIN, SFMAX &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; real PC, UC, UL, PL, D1, D2, K1, K2, K3&lt;br&gt;&nbsp; &nbsp; &nbsp; real pre_sij(6), hydro_strain_Inc, eij_Inc(6), tr_sij(6)&lt;br&gt;&nbsp; &nbsp; &nbsp; real tr_sij_J2, Q_trial, mu_old, mu_new, mu_max, mu_plastic&lt;br&gt;&nbsp; &nbsp; &nbsp; real P_new, P_max, P_old_stress, P_star, T_star, T_current&lt;br&gt;&nbsp; &nbsp; &nbsp; real eps_dot_star, strainRate_eq, Yield_star, Yield_HJC&lt;br&gt;&nbsp; &nbsp; &nbsp; real eps_p_f, strainOld, Damage, e_vol_inc, d_mu_p&lt;br&gt;&nbsp; &nbsp; &nbsp; real smallNum, k_value, true_sij(6), P_sum_check&lt;br&gt;&nbsp; &nbsp; &nbsp; real K_unload, F_interp&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; real e_plaStrain_Inc, mu_bar, mu_plaStrain&lt;br&gt;&nbsp; &nbsp; &nbsp; parameter (zero = 0.0d0, half = 0.5d0, one = 1.0d0,&lt;br&gt;&nbsp; &nbsp; &nbsp;1 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; two = 2.0d0, three = 3.0d0)&lt;br&gt;&nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; smallNum = 1.0d-20&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;c &nbsp; &nbsp; Get material parameters&lt;br&gt;&nbsp; &nbsp; &nbsp; G = cm(1)&lt;br&gt;&nbsp; &nbsp; &nbsp; K_elastic = cm(2)&lt;br&gt;&nbsp; &nbsp; &nbsp; A = cm(3)&lt;br&gt;&nbsp; &nbsp; &nbsp; B = cm(4)&lt;br&gt;&nbsp; &nbsp; &nbsp; C = cm(5)&lt;br&gt;&nbsp; &nbsp; &nbsp; N = cm(6)&lt;br&gt;&nbsp; &nbsp; &nbsp; FC = cm(7)&lt;br&gt;&nbsp; &nbsp; &nbsp; T = cm(8)&lt;br&gt;&nbsp; &nbsp; &nbsp; EPSILON_RATE = cm(9)&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; EFMIN = cm(10)&lt;br&gt;&nbsp; &nbsp; &nbsp; SFMAX = cm(11)&lt;br&gt;&nbsp; &nbsp; &nbsp; PC = cm(12)&lt;br&gt;&nbsp; &nbsp; &nbsp; UC = cm(13)&lt;br&gt;&nbsp; &nbsp; &nbsp; PL = cm(14)&lt;br&gt;&nbsp; &nbsp; &nbsp; UL = cm(15)&lt;br&gt;&nbsp; &nbsp; &nbsp; D1 = cm(16)&lt;br&gt;&nbsp; &nbsp; &nbsp; D2 = cm(17)&lt;br&gt;&nbsp; &nbsp; &nbsp; K1 = cm(18)&lt;br&gt;&nbsp; &nbsp; &nbsp; K2 = cm(19)&lt;br&gt;&nbsp; &nbsp; &nbsp; K3 = cm(20)&lt;\/p&gt;&lt;p&gt;&nbsp; &nbsp; &nbsp; if (etype.eq.&#8217;solid&#8217;.or.etype.eq.&#8217;shl_t&#8217;) then&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Step 1: Extract history variables and preprocessing&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;&nbsp; &nbsp; &nbsp; P_old_stress &nbsp; &nbsp;= hsv(1) &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; mu_old &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;= hsv(3) &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; strainOld &nbsp; &nbsp; &nbsp; = hsv(5) &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; Damage &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;= hsv(14)&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; mu_max &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;= hsv(15)&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; P_max &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; = hsv(16)&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; mu_plaStrain &nbsp; &nbsp;= hsv(2)&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;c &nbsp; &nbsp; Calculate equivalent strain rate&lt;br&gt;&nbsp; &nbsp; &nbsp; strainRate_eq = sqrt(two\/three * (eps(1)**2 + eps(2)**2 + &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp;1 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;eps(3)**2 + half*(eps(4)**2 + eps(5)**2 +&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp;2 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;eps(6)**2))) \/ max(dt1, smallNum)&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; eps_dot_star = strainRate_eq \/ max(EPSILON_RATE, smallNum)&lt;br&gt;&nbsp; &nbsp; &nbsp; if (eps_dot_star .lt. one) eps_dot_star = one&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; Elastic predictor (deviatoric stress)&lt;br&gt;&nbsp; &nbsp; &nbsp; pre_sij(1:3) = sig(1:3) + P_old_stress&lt;br&gt;&nbsp; &nbsp; &nbsp; pre_sij(4:6) = sig(4:6)&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; e_vol_inc = (eps(1) + eps(2) + eps(3)) \/ three&lt;br&gt;&nbsp; &nbsp; &nbsp; eij_Inc(1:3) = eps(1:3) &#8211; e_vol_inc&lt;br&gt;&nbsp; &nbsp; &nbsp; eij_Inc(4:6) = eps(4:6)&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; tr_sij(1:3) = pre_sij(1:3) + two * G * eij_Inc(1:3)&lt;br&gt;&nbsp; &nbsp; &nbsp; tr_sij(4:6) = pre_sij(4:6) + G * eij_Inc(4:6)&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; tr_sij_J2 = half * (tr_sij(1)**2 + tr_sij(2)**2 + tr_sij(3)**2)&lt;br&gt;&nbsp; &nbsp; &nbsp;1 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;+ tr_sij(4)**2 + tr_sij(5)**2 + tr_sij(6)**2&lt;br&gt;&nbsp; &nbsp; &nbsp; Q_trial = sqrt(three * tr_sij_J2)&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Step 2: Equation of state (EOS) and pressure parameter update&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;&nbsp; &nbsp; &nbsp; hydro_strain_Inc = &#8211; (eps(1) + eps(2) + eps(3)) &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; mu_new = mu_old + hydro_strain_Inc &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;c &nbsp; &nbsp; [Improvement] Precompute dynamic unloading and tensile bulk modulus K_unload&nbsp;&lt;br&gt;c &nbsp; &nbsp; Reference: Interpolate between K_elastic and K1 based on compaction level&lt;br&gt;&nbsp; &nbsp; &nbsp; if (mu_max .le. UC) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;K_unload = K_elastic&lt;br&gt;&nbsp; &nbsp; &nbsp; else if (mu_max .lt. UL) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;F_interp = (mu_max &#8211; UC) \/ (UL &#8211; UC)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;K_unload = (one &#8211; F_interp) * K_elastic + F_interp * K1&lt;br&gt;&nbsp; &nbsp; &nbsp; else&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;K_unload = K1&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; end if&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; if (mu_new .ge. mu_max .and. mu_new .gt. zero) then&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; ! &#8212; Loading path (skeleton curve) &#8212;&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; if (mu_new .le. UC) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; P_new = K_elastic * mu_new&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; else if (mu_new .le. UL) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; P_new = PC + (mu_new &#8211; UC)*(PL &#8211; PC)\/(UL &#8211; UC)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; else&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mu_bar = (mu_new &#8211; UL) \/ (one + UL)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; P_new = PL + K1 * mu_bar + K2 * mu_bar**2 + K3 * mu_bar**3&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; end if&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; mu_max = mu_new&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; P_max &nbsp;= P_new&lt;br&gt;&nbsp; &nbsp; &nbsp; else&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; ! &#8212; Unloading and tension &#8212;&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; if (mu_new .lt. zero) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ! Tension: use degraded unloading modulus according to literature&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; P_new = K_unload * mu_new&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; else&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ! Compressive unloading: linear return from P_max using dynamic K_unload&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; P_new = P_max &#8211; K_unload * (mu_max &#8211; mu_new)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; if (P_new .lt. zero) P_new = zero&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; end if&lt;br&gt;&nbsp; &nbsp; &nbsp; end if&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; Tension cutoff&lt;br&gt;&nbsp; &nbsp; &nbsp; T_current = T * (one &#8211; Damage)&lt;br&gt;&nbsp; &nbsp; &nbsp; if (P_new .lt. -T_current) P_new = -T_current&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; [Improvement] Calculate plastic volumetric strain increment d_mu_p&nbsp;&lt;br&gt;c &nbsp; &nbsp; Since modulus is dynamic K_unload, the elastic volumetric strain subtraction should also use current modulus&lt;br&gt;&nbsp; &nbsp; &nbsp; d_mu_p = zero&lt;br&gt;&nbsp; &nbsp; &nbsp; if (mu_new .gt. UC .and. hydro_strain_Inc .gt. zero) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;d_mu_p = hydro_strain_Inc &#8211; (P_new &#8211; P_old_stress)\/K_unload&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;if (d_mu_p .lt. zero) d_mu_p = zero&lt;br&gt;&nbsp; &nbsp; &nbsp; end if&lt;br&gt;&nbsp; &nbsp; &nbsp; mu_plaStrain = mu_plaStrain + d_mu_p&lt;br&gt;c &nbsp; &nbsp; Update P* and T* immediately (for denominator and HSV6 calculation)&lt;br&gt;&nbsp; &nbsp; &nbsp; P_star = P_new \/ max(FC, smallNum)&lt;br&gt;&nbsp; &nbsp; &nbsp; T_star = T \/ max(FC, smallNum)&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Step 3: Damage denominator calculation (NaN protection)&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;&nbsp; &nbsp; &nbsp; P_sum_check = P_star + T_star&lt;br&gt;&nbsp; &nbsp; &nbsp; if (P_sum_check .le. 1.0d-8) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; eps_p_f = EFMIN&lt;br&gt;&nbsp; &nbsp; &nbsp; else&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; eps_p_f = D1 * (P_sum_check)**D2&lt;br&gt;&nbsp; &nbsp; &nbsp; end if&lt;br&gt;&nbsp; &nbsp; &nbsp; if (eps_p_f .lt. EFMIN) eps_p_f = EFMIN &nbsp;&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Step 4: Yield surface and radial return (fix multiplication order of SFMAX and strain rate)&lt;br&gt;c &nbsp; &nbsp; Step 4: Failure fuse (advanced protection)&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; If damage exceeds 0.999d0, trigger fuse early to avoid numerical noise&lt;br&gt;&nbsp; &nbsp; &nbsp; if (Damage .gt. 0.999999d0) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Damage = 1.0d0 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Yield_HJC = 0.001d0 * FC &nbsp; ! residual strength&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;e_plaStrain_Inc = zero&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;k_value = Yield_HJC \/ max(Q_trial, smallNum)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;true_sij(1:6) = tr_sij(1:6) * k_value&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; goto 999&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; endif&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; Yield surface calculation&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Yield_star = A * (one &#8211; Damage) + B * (P_star**N)* (one&lt;br&gt;&nbsp; &nbsp; &nbsp;1 &nbsp; &nbsp; &nbsp;+ C * log(eps_dot_star))&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;if (Yield_star .gt. SFMAX) Yield_star = SFMAX&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; &nbsp; &nbsp;Add safety cutoff: prevent Yield_star from becoming negative due to numerical fluctuations&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;if (Yield_star .lt. zero) Yield_star = zero&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; Yield_HJC = Yield_star * FC&lt;\/p&gt;&lt;p&gt;&nbsp; &nbsp; &nbsp; if (Q_trial .gt. Yield_HJC) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; e_plaStrain_Inc = (Q_trial &#8211; Yield_HJC) \/ (three * G)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; strainOld = strainOld + e_plaStrain_Inc&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; k_value = Yield_HJC \/ max(Q_trial, smallNum)&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; true_sij(1:6) = tr_sij(1:6) * k_value&lt;br&gt;&nbsp; &nbsp; &nbsp; else&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; e_plaStrain_Inc = zero&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; true_sij(1:6) = tr_sij(1:6)&lt;br&gt;&nbsp; &nbsp; &nbsp; end if&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Step 5: Damage accumulation (numerical tolerance completion)&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Only accumulate if damage has not reached the set &#8216;quasi-failure threshold'&lt;br&gt;&nbsp; &nbsp; &nbsp; if (Damage .lt. 0.999999d0) then&lt;br&gt;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Damage = Damage + (e_plaStrain_Inc + d_mu_p) \/ eps_p_f&lt;br&gt;&nbsp; &nbsp; &nbsp; end if&lt;\/p&gt;&lt;p&gt;c &nbsp; &nbsp; [Core modification]: Force numerical closure&lt;br&gt;c &nbsp; &nbsp; Once damage exceeds the threshold of 0.999999, directly set to 1.0&lt;br&gt;c &nbsp; &nbsp; This solves the precision error concern and ensures complete red in contour&lt;br&gt;&nbsp; &nbsp; &nbsp; if (Damage .ge. 0.999999d0) Damage = one&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;c &nbsp; &nbsp; Step 6: Update history variables (meeting user requirements: hsv2, hsv6)&lt;br&gt;c &nbsp; &nbsp; =====================================================================&lt;br&gt;999 &nbsp; continue&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; sig(1:3) = true_sij(1:3) &#8211; P_new&lt;br&gt;&nbsp; &nbsp; &nbsp; sig(4:6) = true_sij(4:6)&lt;br&gt;&nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(1) &nbsp;= P_new &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(2) &nbsp;= mu_plaStrain &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(3) &nbsp;= mu_new &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(4) &nbsp;= e_plaStrain_Inc &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(5) &nbsp;= strainOld &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(6) &nbsp;= D1 * (P_sum_check)**D2 &nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(7:12) = eps(1:6)&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(13) = eps_dot_star &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(14) = Damage &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(15) = mu_max &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; hsv(16) = P_max &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&lt;br&gt;&nbsp; &nbsp; &nbsp; epsp &nbsp; &nbsp;= strainOld&lt;\/p&gt;&lt;p&gt;&nbsp; &nbsp; &nbsp; end if&lt;br&gt;&nbsp; &nbsp; &nbsp; return&lt;br&gt;&nbsp; &nbsp; &nbsp; end&lt;\/p&gt;&lt;p&gt;HJC CODE END&lt;\/p&gt;&lt;p&gt;&nbsp;&lt;\/p&gt;&lt;p&gt;LS-DYNA-HJC&nbsp;&lt;\/p&gt;&lt;p&gt;<img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/innovationspace.ansys.com\/forum\/wp-content\/uploads\/sites\/2\/2026\/05\/06-05-2026-1778058347-MAT111_FRONT.png\" alt=\"\" width=\"692\" height=\"319\" \/>&lt;\/p&gt;&lt;p&gt;<img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/innovationspace.ansys.com\/forum\/wp-content\/uploads\/sites\/2\/2026\/05\/06-05-2026-1778058365-MAT111_TOP_DAMAGE.png\" alt=\"\" width=\"834\" height=\"505\" \/>&lt;\/p&gt;&lt;p&gt;<img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/innovationspace.ansys.com\/forum\/wp-content\/uploads\/sites\/2\/2026\/05\/06-05-2026-1778058374-MAT111_BOTTOM_DAMAGE.png\" alt=\"\" width=\"836\" height=\"547\" \/>&lt;\/p&gt;&lt;p&gt;&nbsp;&lt;\/p&gt;&lt;p&gt;UMAT-HJC&lt;\/p&gt;&lt;p&gt;<img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/innovationspace.ansys.com\/forum\/wp-content\/uploads\/sites\/2\/2026\/05\/06-05-2026-1778058432-mceclip0.png\" width=\"833\" height=\"534\" \/>&lt;\/p&gt;&lt;p&gt;<img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/innovationspace.ansys.com\/forum\/wp-content\/uploads\/sites\/2\/2026\/05\/06-05-2026-1778058449-mceclip1.png\" width=\"828\" height=\"554\" \/>&lt;\/p&gt;&lt;p&gt;<img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/innovationspace.ansys.com\/forum\/wp-content\/uploads\/sites\/2\/2026\/05\/06-05-2026-1778058504-mceclip2.png\" width=\"828\" height=\"410\" \/>&lt;\/p&gt;&lt;p&gt;&nbsp;&lt;\/p&gt;<\/p>\n","protected":false},"template":"","class_list":["post-457660","topic","type-topic","status-publish","hentry"],"aioseo_notices":[],"acf":[],"custom_fields":[{"0":{"_bbp_forum_id":["27814"],"_bbp_topic_id":["457660"],"_bbp_author_ip":["2001:da8:d800:146:4d05:a375:c950:5a35"],"_bbp_last_reply_id":["457688"],"_bbp_last_active_id":["457688"],"_bbp_last_active_time":["2026-05-06 12:22:38"],"_bbp_reply_count":["4"],"_bbp_reply_count_hidden":["0"],"_bbp_voice_count":["2"],"_bbp_engagement":["684302","58821"],"_btv_view_count":["86"],"_bbp_topic_status":["unanswered"],"_bbp_subscription":["684302","58821"]},"test":"hscjmail-ustc-edu-cn"}],"_links":{"self":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics\/457660","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics"}],"about":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/types\/topic"}],"version-history":[{"count":0,"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics\/457660\/revisions"}],"wp:attachment":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/media?parent=457660"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}