HJC模型的损伤演化问题？
浏览：46
大家好，我正在进行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
