#

# 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]