*UNIT_SYSTEM
SI
*SCRIPT_PYTHON
damage.py
*TIME
5.0e-4
#
# --- SENSORS ---
#
*OUTPUT_SENSOR
"center"
10, 1, 0, 0, 0
"right"
11, 1, 0.05, 0, 0
*OUTPUT_SENSOR_EXTENDED
"damage"
55
10, 11, 5.0e-5
"dmg_cl"
"dmg_jc"
"epsp"
#
# --- MESH ---
#
*COMPONENT_BOX
"sample"
1, 1, 20, 2, 2
-0.1, -0.01, -0.01, 0.1, 0.01, 0.01
#
# --- MATERIAL ---
#
*MAT_METAL
1, 7800.0, 210.0e9, 0.3, 0, 1
1
*FUNCTION
1
5.0e8 + 5.0e8*epsp^0.3
*PROP_THERMAL
1, 0.0, 450.0
#
# --- PART ---
#
*PART
"sample"
1, 1
#
# --- INITIAL VELOCITY ---
#
*INITIAL_VELOCITY
P, 1, fcn(123)
*FUNCTION
123
3000*x
*END
#
# COMPARISON OF DAMAGE MODELS
# dmg_cl - Cockcroft-Latham
# dmg_jc - Johnson-Cook
#
import math
# damage parameters
Wc = 1.0e9
d1 = 0.1
d2 = 0.3
d3 = 0.5
d4 = 0.01
d5 = 0.7
eps_0 = 1.0
T_0 = 0.0
T_m = 1500.0
#
# INITIALIZE DAMAGE
def init():
return [0, 0, 0]
#
# UPDATE DAMAGE
def update(dmg_cl, dmg_jc, epsp_old,
t, dt, epsp, epsc, dmg, T, sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_zx,
eps_xx, eps_yy, eps_zz, eps_xy, eps_yz, eps_zx, vx, vy, vz, x, y, z):
# damage parameters
global Wc, d1, d2, d3, d4, d5, eps_0, T_0, T_m
# incremental plastic strain and strain rate
depsp = epsp - epsp_old
rate = depsp / dt
# principal stresses
[s1, s2, s3] = eigval(sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_zx)
# pressure and von Mises stress
p =-(s1+s2+s3)/3.0
sig_vm = math.sqrt(1.5*( (s1+p)**2 + (s2+p)**2 + (s3+p)**2 ))
# damage increment
if (depsp > 0):
dmg_cl = dmg_cl + max(0,s1,s2,s3)*depsp / Wc
eps_f = (d1 + d2*math.exp(d3*p/sig_vm)) * (1.0 + d4*math.log(rate/eps_0)) * (1.0 + d5*(T - T_0)/(T_m - T_0))
dmg_jc = dmg_jc + depsp / eps_f
return [dmg_cl, dmg_jc, epsp]
#
# COMPUTE PRINCIPAL STRESSES
def eigval(sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_zx):
# constant, linear and quadratic terms in characteristic equation f = x^3 + a2*x^2 + a1*x + a0 = 0
a0 =-sig_xx*sig_yy*sig_zz + sig_xx*sig_yz**2 + sig_xy**2*sig_zz + sig_zx**2*sig_yy - 2.0*sig_xy*sig_yz*sig_zx
a1 =-sig_xy**2 - sig_yz**2 - sig_zx**2 + sig_yy*sig_zz + sig_xx*sig_yy + sig_xx*sig_zz
a2 =-(sig_xx+sig_yy+sig_zz)
Q = a1/3.0 - a2**2/9.0
Q = min(0.0,Q)
R = a1*a2/6.0 - a0/2.0 - a2**3/27.0
T =-a2/3.0
if (Q > -1.0e-10):
cc = 0.0
ss = 0.0
else:
S = max(-1.0,min(1.0,R/math.sqrt(-Q**3)))
theta = math.acos(S)
Q2 = 2.0*math.sqrt(-Q)
cc = Q2*math.cos(theta/3.0)
ss = Q2*math.sin(theta/3.0)
# principal stresses
s1 = cc + T
s2 =-cc/2.0 - math.sqrt(3.0)/2.0*ss + T
s3 =-cc/2.0 + math.sqrt(3.0)/2.0*ss + T
return [s1, s2, s3]