HJC模型的損傷演化問題？
瀏覽：570
大家好，我正在進行HJC模型的二次開發，我想從最簡單的HJC模型做起，目前單單元已經通過驗證了，但是在質量塊撞擊薄板模型中的損傷演化與應力應變更新與LS-DYNA自帶的HJC模型的不同，我查閱了許多文獻與問了許多人，嘗試了許多方法，改了很多遍，還是改的不一樣，源碼在此，請各位幫我看看有什么問題，我感覺主要問題應該在徑向回歸與損傷累計那塊代碼出現問題。請各位幫我看看，謝謝大家了。
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 變量定義
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 獲取屬性參數
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 第一步：提取歷史變量與預處理
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 計算等效應變率
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 彈性預測 (偏應力)
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 第二步：狀態方程 (EOS) 與 壓力參數更新
c =====================================================================
hydro_strain_Inc = - (eps(1) + eps(2) + eps(3))
mu_new = mu_old + hydro_strain_Inc
c [改進] 預先計算動態卸載與受拉體積模量 K_unload
c 文獻依據：根據材料壓實程度，在 K_elastic 和 K1 之間插值
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
! --- 加載階段 (骨架曲線) ---
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
! --- 卸載與受拉階段 ---
if (mu_new .lt. zero) then
! 受拉狀態：使用插值退化后的模量，符合文獻中對于過渡區受拉的處理方式
P_new = K_unload * mu_new
else
! 受壓卸載：從 P_max 按照動態 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 拉伸截斷
T_current = T * (one - Damage)
if (P_new .lt. -T_current) P_new = -T_current
c [改進] 計算塑性體積應變增量 d_mu_p
c 既然模量改為了動態的 K_unload，彈性部分的體積應變扣除也應使用當前模量
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 及時更新 P* 和 T* (用于計算分母和 HSV6)
P_star = P_new / max(FC, smallNum)
T_star = T / max(FC, smallNum)
c =====================================================================
c 第三步：損傷分母計算 (防 NaN 保護)
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 第四步：屈服面與徑向回歸 (修復 SFMAX 與應變率的乘法順序)
c 第四步：失效熔斷 (進階防御版)
c =====================================================================
c 只要損傷超過了 0.999d0，就提前執行熔斷，不給它產生數值噪音的機會
if (Damage .gt. 0.999999d0) then
Damage = 1.0d0 ! 強行補齊，讓它在云圖上變紅
Yield_HJC = 0.001d0 * FC !殘余強度
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 =====================================================================
c 屈服面計算
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 增加安全截斷：防止由于數值波動導致 Yield_star 算出負數
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 第五步：損傷累加 (數值容差補齊版)
c =====================================================================
c 只有在損傷還未達到我們設定的“準失效閾值”時，才進行累加
if (Damage .lt. 0.999999d0) then
Damage = Damage + (e_plaStrain_Inc + d_mu_p) / eps_p_f
end if
c 【核心修改】：數值強制閉合
c 一旦損傷跨過了 0.999999 這個門檻，直接手動“填平”到 1.0
c 這樣既解決了你擔心的精度誤差問題，又確保了云圖紅得徹底
if (Damage .ge. 0.999999d0) Damage = one
c =====================================================================
c 第六步：更新歷史變量 (滿足用戶要求：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
