223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

1. Implement Karafillis-Boyce yield criterion (the residum and Jacobian); 

2. Karafillis-Boyce works
This commit is contained in:
Haiming Zhang 2015-02-24 16:46:09 +00:00
parent 2baac6fc10
commit 8b1846e78e
1 changed files with 209 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

271
processing/misc/yieldSurface.py
View File

@ -135,14 +135,6 @@ class Hosford(object):
def jac(self, (sigma0, a), ydata, sigmas): def jac(self, (sigma0, a), ydata, sigmas):
return HosfordBasis(sigma0, 1.0,1.0,1.0, a, sigmas, Jac=True, nParas=2) return HosfordBasis(sigma0, 1.0,1.0,1.0, a, sigmas, Jac=True, nParas=2)
#more to do
# KarafillisAndBoyce
# ---------------------------------------------------------------------------------------------
# isotropic yield surfaces
# ---------------------------------------------------------------------------------------------
class Hill1948(object): class Hill1948(object):
''' '''
residuum of Hill 1948 quadratic yield criterion (eq. 2.48) residuum of Hill 1948 quadratic yield criterion (eq. 2.48)
@ -160,11 +152,6 @@ class Hill1948(object):
for f,g,h,l,m,n in zip(jF, jG, jH, jL, jM, jN): jaco.append([f,g,h,l,m,n]) for f,g,h,l,m,n in zip(jF, jG, jH, jL, jM, jN): jaco.append([f,g,h,l,m,n])
return np.array(jaco) return np.array(jaco)
#more to do
# Hill 1979
# Hill 1990,1993 need special stresses to fit
class generalHosford(object): class generalHosford(object):
''' '''
residuum of Hershey yield criterion (eq. 2.104, sigmas = sigma0) residuum of Hershey yield criterion (eq. 2.104, sigmas = sigma0)
@ -212,6 +199,22 @@ class Yld200418p(object):
d12,d21,d23,d32,d31,d13,d44,d55,d66, m), ydata, sigmas): d12,d21,d23,d32,d31,d13,d44,d55,d66, m), ydata, sigmas):
return Yld200418pBasis(sigma0, c12,c21,c23,c32,c31,c13,c44,c55,c66, return Yld200418pBasis(sigma0, c12,c21,c23,c32,c31,c13,c44,c55,c66,
d12,d21,d23,d32,d31,d13,d44,d55,d66, m, sigmas, Jac=True) d12,d21,d23,d32,d31,d13,d44,d55,d66, m, sigmas, Jac=True)
class KarafillisBoyce(object):
'''
residuum of Karafillis-Boyce yield criterion
'''
def __init__(self, uniaxialStress):
self.stress0 = uniaxialStress
def fun(self, (c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,
b1, b2, a, alpha), ydata, sigmas):
return KarafillisBoyceBasis(self.stress0, c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,
b1, b2, a, alpha, sigmas)
def jac(self, (c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,
b1, b2, a, alpha), ydata, sigmas):
return KarafillisBoyceBasis(self.stress0, c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,
b1, b2, a, alpha, sigmas, Jac=True)
class BBC2003(object): class BBC2003(object):
''' '''
residuum of anisotropic Barlat 1991 yield criterion (eq. 2.37) residuum of anisotropic Barlat 1991 yield criterion (eq. 2.37)
@ -258,7 +261,7 @@ def Cazacu_Barlat3DBasis(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,c,
s1_3, s2_3, s3_3 = s11*s1_2, s22*s2_2, s33*s3_2 s1_3, s2_3, s3_3 = s11*s1_2, s22*s2_2, s33*s3_2
s12_2, s23_2, s31_2 = s12**2, s23**2, s31**2 s12_2, s23_2, s31_2 = s12**2, s23**2, s31**2
d12_2, d23_2, d31_2 = (s11-s22)**2, (s22-s33)**2, (s33-s11)**2 d12_2, d23_2, d31_2 = (s11-s22)**2, (s22-s33)**2, (s33-s11)**2
J20 = ( a1*d12_2 + a2*d23_2 + a3*d31_2 )/6.0 + a4*s12_2 + a5*s23_2 + a6*s31_2 J20 = ( a1*d12_2 + a2*d23_2 + a3*d31_2 )/6.0 + a4*s12_2 + a5*s23_2 + a6*s31_2
J30 = ( (b1 +b2 )*s1_3 + (b3 +b4 )*s2_3 + ( b1+b4-b2 + b1+b4-b3 )*s3_3 )/27.0- \ J30 = ( (b1 +b2 )*s1_3 + (b3 +b4 )*s2_3 + ( b1+b4-b2 + b1+b4-b3 )*s3_3 )/27.0- \
( (b1*s22+b2*s33)*s1_2 + (b3*s33+b4*s11)*s2_2 + ((b1+b4-b2)*s11 + (b1+b4-b3)*s22)*s3_2 )/9.0 + \ ( (b1*s22+b2*s33)*s1_2 + (b3*s33+b4*s11)*s2_2 + ((b1+b4-b2)*s11 + (b1+b4-b3)*s22)*s3_2 )/9.0 + \
@ -269,7 +272,7 @@ def Cazacu_Barlat3DBasis(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,c,
)/3.0 )/3.0
f0 = (J20**3 - c*J30**2)/18.0 f0 = (J20**3 - c*J30**2)/18.0
r = f0**(1.0/6.0)*(3.0/sigma0) r = f0**(1.0/6.0)*(3.0/sigma0)
if not Jac: if not Jac:
return (r - 1.0).ravel() return (r - 1.0).ravel()
else: else:
@ -311,7 +314,7 @@ def Cazacu_Barlat2DBasis(a1,a2,a3,a4,b1,b2,b3,b4,b5,b10,c,
f0 = (J20**3 - c*J30**2)/18.0 f0 = (J20**3 - c*J30**2)/18.0
r = f0**(1.0/6.0)*(3.0/sigma0) r = f0**(1.0/6.0)*(3.0/sigma0)
if not Jac: if not Jac:
return (r - 1.0).ravel() return (r - 1.0).ravel()
else: else:
@ -330,7 +333,6 @@ def Cazacu_Barlat2DBasis(a1,a2,a3,a4,b1,b2,b3,b4,b5,b10,c,
jaco.append(jacv) jaco.append(jacv)
return np.array(jaco) return np.array(jaco)
def DruckerBasis(sigma0, C_D, p, sigmas, Jac=False, nParas=2): def DruckerBasis(sigma0, C_D, p, sigmas, Jac=False, nParas=2):
s11 = sigmas[0]; s22 = sigmas[1]; s33 = sigmas[2] s11 = sigmas[0]; s22 = sigmas[1]; s33 = sigmas[2]
s12 = sigmas[3]; s23 = sigmas[4]; s31 = sigmas[5] s12 = sigmas[3]; s23 = sigmas[4]; s31 = sigmas[5]
@ -371,7 +373,7 @@ def HosfordBasis(sigma0, F,G,H, a, sigmas, Jac=False, nParas=1):
base = F*diff23**a + G*diff31**a + H*diff12**a; expo = 1.0/a base = F*diff23**a + G*diff31**a + H*diff12**a; expo = 1.0/a
left = base**expo left = base**expo
right = 2.0**expo*sigma0 right = 2.0**expo*sigma0
if not Jac: if not Jac:
if nParas == 1: return (left - right).ravel() if nParas == 1: return (left - right).ravel()
else: return (left/right - 1.0).ravel() else: return (left/right - 1.0).ravel()
@ -384,7 +386,7 @@ def HosfordBasis(sigma0, F,G,H, a, sigmas, Jac=False, nParas=1):
drda = sigma0*(2.0**expo)*ln(2.0)*deda drda = sigma0*(2.0**expo)*ln(2.0)*deda
dldb = expo*left/base dldb = expo*left/base
jaco = [] jaco = []
if nParas == 1: # von Mises if nParas == 1: # von Mises
return ones*(-2.0**0.5) return ones*(-2.0**0.5)
elif nParas == 2: # isotropic Hosford elif nParas == 2: # isotropic Hosford
@ -402,7 +404,7 @@ def HosfordBasis(sigma0, F,G,H, a, sigmas, Jac=False, nParas=1):
for jacv in zip(js, jF,jG,jH,ja): for jacv in zip(js, jF,jG,jH,ja):
jaco.append(jacv) jaco.append(jacv)
return np.array(jaco) return np.array(jaco)
def Barlat1991Basis(sigma0, a, b, c, f, g, h, m, sigmas, Jac=False, nParas=2): def Barlat1991Basis(sigma0, a, b, c, f, g, h, m, sigmas, Jac=False, nParas=2):
''' '''
residuum of Barlat 1997 yield criterion residuum of Barlat 1997 yield criterion
@ -507,7 +509,7 @@ def BBC2003Basis(sigma0, a,b,c, d,e,f,g, k, sigmas, Jac=False):
dldl1 = a *k*(l1s**k1)*(2.0*l1) dldl1 = a *k*(l1s**k1)*(2.0*l1)
dldl2 = a *k*(l2s**k1)*(2.0*l2) dldl2 = a *k*(l2s**k1)*(2.0*l2)
dldl3 = (1-a)*k*(l3s**k1)*(2.0*l3) dldl3 = (1-a)*k*(l3s**k1)*(2.0*l3)
dldGama = (dldl1 + dldl2)*b dldGama = (dldl1 + dldl2)*b
dldPsi = (dldl1 - dldl2 + 2.0*dldl3)*c dldPsi = (dldl1 - dldl2 + 2.0*dldl3)*c
@ -551,7 +553,7 @@ def principalStress(p):
return np.array([S1,S2,S3]), np.array([I1,I2,I3]) return np.array([S1,S2,S3]), np.array([I1,I2,I3])
def principalStrs_Der(p, Invariant, s1, s2, s3, s4, s5, s6): def principalStrs_Der(p, Invariant, s1, s2, s3, s4, s5, s6, Karafillis=False):
sin = np.sin; cos = np.cos sin = np.sin; cos = np.cos
I1 = Invariant[0,:]; I2 = Invariant[1,:]; I3 = Invariant[2,:] I1 = Invariant[0,:]; I2 = Invariant[1,:]; I3 = Invariant[2,:]
@ -600,49 +602,86 @@ def principalStrs_Der(p, Invariant, s1, s2, s3, s4, s5, s6):
dI3dp1 = p[2]*p[0] - p[5]**2; dI3dp5 = -2.0*p[5]*p[1] + 2.0*p[3]*p[4] dI3dp1 = p[2]*p[0] - p[5]**2; dI3dp5 = -2.0*p[5]*p[1] + 2.0*p[3]*p[4]
dI3dp2 = p[0]*p[1] - p[3]**2; dI3dp3 = -2.0*p[3]*p[2] + 2.0*p[4]*p[5] dI3dp2 = p[0]*p[1] - p[3]**2; dI3dp3 = -2.0*p[3]*p[2] + 2.0*p[4]*p[5]
dI1dc12 = dI1dp0*(-s2); dI2dc12 = dI2dp0*(-s2); dI3dc12 = dI3dp0*(-s2) # c12 if Karafillis:
dI1dc21 = dI1dp1*(-s1); dI2dc21 = dI2dp1*(-s1); dI3dc21 = dI3dp1*(-s1) # c21 dS1dp0 = dS1dI1*dI1dp0 + dS1dI2*dI2dp0 + dS1dI3*dI3dp0
dI1dc23 = dI1dp1*(-s3); dI2dc23 = dI2dp1*(-s3); dI3dc23 = dI3dp1*(-s3) # c23 dS1dp1 = dS1dI1*dI1dp1 + dS1dI2*dI2dp1 + dS1dI3*dI3dp1
dI1dc32 = dI1dp2*(-s2); dI2dc32 = dI2dp2*(-s2); dI3dc32 = dI3dp2*(-s2) # c32 dS1dp2 = dS1dI1*dI1dp2 + dS1dI2*dI2dp2 + dS1dI3*dI3dp2
dI1dc31 = dI1dp2*(-s1); dI2dc31 = dI2dp2*(-s1); dI3dc31 = dI3dp2*(-s1) # c31 dS1dp3 = dS1dI1*dI1dp3 + dS1dI2*dI2dp3 + dS1dI3*dI3dp3
dI1dc13 = dI1dp0*(-s3); dI2dc13 = dI2dp0*(-s3); dI3dc13 = dI3dp0*(-s3) # c13 dS1dp4 = dS1dI1*dI1dp4 + dS1dI2*dI2dp4 + dS1dI3*dI3dp4
dI1dc44 = dI1dp3* s4 ; dI2dc44 = dI2dp3* s4 ; dI3dc44 = dI3dp3* s4 # c44 dS1dp5 = dS1dI1*dI1dp5 + dS1dI2*dI2dp5 + dS1dI3*dI3dp5
dI1dc55 = dI1dp4* s5 ; dI2dc55 = dI2dp4* s5 ; dI3dc55 = dI3dp4* s5 # c55 dS2dp0 = dS2dI1*dI1dp0 + dS2dI2*dI2dp0 + dS2dI3*dI3dp0
dI1dc66 = dI1dp5* s6 ; dI2dc66 = dI2dp5* s6 ; dI3dc66 = dI3dp5* s6 # c66 dS2dp1 = dS2dI1*dI1dp1 + dS2dI2*dI2dp1 + dS2dI3*dI3dp1
dS2dp2 = dS2dI1*dI1dp2 + dS2dI2*dI2dp2 + dS2dI3*dI3dp2
dS2dp3 = dS2dI1*dI1dp3 + dS2dI2*dI2dp3 + dS2dI3*dI3dp3
dS2dp4 = dS2dI1*dI1dp4 + dS2dI2*dI2dp4 + dS2dI3*dI3dp4
dS2dp5 = dS2dI1*dI1dp5 + dS2dI2*dI2dp5 + dS2dI3*dI3dp5
dS3dp0 = dS3dI1*dI1dp0 + dS3dI2*dI2dp0 + dS3dI3*dI3dp0
dS3dp1 = dS3dI1*dI1dp1 + dS3dI2*dI2dp1 + dS3dI3*dI3dp1
dS3dp2 = dS3dI1*dI1dp2 + dS3dI2*dI2dp2 + dS3dI3*dI3dp2
dS3dp3 = dS3dI1*dI1dp3 + dS3dI2*dI2dp3 + dS3dI3*dI3dp3
dS3dp4 = dS3dI1*dI1dp4 + dS3dI2*dI2dp4 + dS3dI3*dI3dp4
dS3dp5 = dS3dI1*dI1dp5 + dS3dI2*dI2dp5 + dS3dI3*dI3dp5
dS1dc12 = dS1dI1 * dI1dc12 + dS1dI2 * dI2dc12 + dS1dI3 * dI3dc12 dS1dc1 = ( dS1dp0*0.0 + dS1dp1*(s2-s3) + dS1dp2*(s3-s2) )/3.0
dS1dc21 = dS1dI1 * dI1dc21 + dS1dI2 * dI2dc21 + dS1dI3 * dI3dc21 dS1dc2 = ( dS1dp0*(s1-s3) + dS1dp1*0.0 + dS1dp2*(s3-s1) )/3.0
dS1dc23 = dS1dI1 * dI1dc23 + dS1dI2 * dI2dc23 + dS1dI3 * dI3dc23 dS1dc3 = ( dS1dp0*(s1-s2) + dS1dp1*(s2-s1) + dS1dp2*0.0 )/3.0
dS1dc32 = dS1dI1 * dI1dc32 + dS1dI2 * dI2dc32 + dS1dI3 * dI3dc32 dS2dc1 = ( dS2dp0*0.0 + dS2dp1*(s2-s3) + dS2dp2*(s3-s2) )/3.0
dS1dc31 = dS1dI1 * dI1dc31 + dS1dI2 * dI2dc31 + dS1dI3 * dI3dc31 dS2dc2 = ( dS2dp0*(s1-s3) + dS2dp1*0.0 + dS2dp2*(s3-s1) )/3.0
dS1dc13 = dS1dI1 * dI1dc13 + dS1dI2 * dI2dc13 + dS1dI3 * dI3dc13 dS2dc3 = ( dS2dp0*(s1-s2) + dS2dp1*(s2-s1) + dS2dp2*0.0 )/3.0
dS1dc44 = dS1dI1 * dI1dc44 + dS1dI2 * dI2dc44 + dS1dI3 * dI3dc44 dS3dc1 = ( dS3dp0*0.0 + dS3dp1*(s2-s3) + dS3dp2*(s3-s2) )/3.0
dS1dc55 = dS1dI1 * dI1dc55 + dS1dI2 * dI2dc55 + dS1dI3 * dI3dc55 dS3dc2 = ( dS3dp0*(s1-s3) + dS3dp1*0.0 + dS3dp2*(s3-s1) )/3.0
dS1dc66 = dS1dI1 * dI1dc66 + dS1dI2 * dI2dc66 + dS1dI3 * dI3dc66 dS3dc3 = ( dS3dp0*(s1-s2) + dS3dp1*(s2-s1) + dS3dp2*0.0 )/3.0
dS1dc4 = dS1dp3*s4; dS1dc5 = dS1dp4*s5; dS1dc6 = dS1dp5*s6
dS2dc4 = dS2dp3*s4; dS2dc5 = dS2dp4*s5; dS2dc6 = dS2dp5*s6
dS3dc4 = dS3dp3*s4; dS3dc5 = dS3dp4*s5; dS3dc6 = dS3dp5*s6
dS2dc12 = dS2dI1 * dI1dc12 + dS2dI2 * dI2dc12 + dS2dI3 * dI3dc12 return dS1dc1,dS1dc2,dS1dc3,dS1dc4,dS1dc5,dS1dc6,\
dS2dc21 = dS2dI1 * dI1dc21 + dS2dI2 * dI2dc21 + dS2dI3 * dI3dc21 dS2dc1,dS2dc2,dS2dc3,dS2dc4,dS2dc5,dS2dc6,\
dS2dc23 = dS2dI1 * dI1dc23 + dS2dI2 * dI2dc23 + dS2dI3 * dI3dc23 dS3dc1,dS3dc2,dS3dc3,dS3dc4,dS3dc5,dS3dc6
dS2dc32 = dS2dI1 * dI1dc32 + dS2dI2 * dI2dc32 + dS2dI3 * dI3dc32 else:
dS2dc31 = dS2dI1 * dI1dc31 + dS2dI2 * dI2dc31 + dS2dI3 * dI3dc31 dI1dc12 = dI1dp0*(-s2); dI2dc12 = dI2dp0*(-s2); dI3dc12 = dI3dp0*(-s2) # c12
dS2dc13 = dS2dI1 * dI1dc13 + dS2dI2 * dI2dc13 + dS2dI3 * dI3dc13 dI1dc21 = dI1dp1*(-s1); dI2dc21 = dI2dp1*(-s1); dI3dc21 = dI3dp1*(-s1) # c21
dS2dc44 = dS2dI1 * dI1dc44 + dS2dI2 * dI2dc44 + dS2dI3 * dI3dc44 dI1dc23 = dI1dp1*(-s3); dI2dc23 = dI2dp1*(-s3); dI3dc23 = dI3dp1*(-s3) # c23
dS2dc55 = dS2dI1 * dI1dc55 + dS2dI2 * dI2dc55 + dS2dI3 * dI3dc55 dI1dc32 = dI1dp2*(-s2); dI2dc32 = dI2dp2*(-s2); dI3dc32 = dI3dp2*(-s2) # c32
dS2dc66 = dS2dI1 * dI1dc66 + dS2dI2 * dI2dc66 + dS2dI3 * dI3dc66 dI1dc31 = dI1dp2*(-s1); dI2dc31 = dI2dp2*(-s1); dI3dc31 = dI3dp2*(-s1) # c31
dI1dc13 = dI1dp0*(-s3); dI2dc13 = dI2dp0*(-s3); dI3dc13 = dI3dp0*(-s3) # c13
dI1dc44 = dI1dp3* s4 ; dI2dc44 = dI2dp3* s4 ; dI3dc44 = dI3dp3* s4 # c44
dI1dc55 = dI1dp4* s5 ; dI2dc55 = dI2dp4* s5 ; dI3dc55 = dI3dp4* s5 # c55
dI1dc66 = dI1dp5* s6 ; dI2dc66 = dI2dp5* s6 ; dI3dc66 = dI3dp5* s6 # c66
dS3dc12 = dS3dI1 * dI1dc12 + dS3dI2 * dI2dc12 + dS3dI3 * dI3dc12 dS1dc12 = dS1dI1 * dI1dc12 + dS1dI2 * dI2dc12 + dS1dI3 * dI3dc12
dS3dc21 = dS3dI1 * dI1dc21 + dS3dI2 * dI2dc21 + dS3dI3 * dI3dc21 dS1dc21 = dS1dI1 * dI1dc21 + dS1dI2 * dI2dc21 + dS1dI3 * dI3dc21
dS3dc23 = dS3dI1 * dI1dc23 + dS3dI2 * dI2dc23 + dS3dI3 * dI3dc23 dS1dc23 = dS1dI1 * dI1dc23 + dS1dI2 * dI2dc23 + dS1dI3 * dI3dc23
dS3dc32 = dS3dI1 * dI1dc32 + dS3dI2 * dI2dc32 + dS3dI3 * dI3dc32 dS1dc32 = dS1dI1 * dI1dc32 + dS1dI2 * dI2dc32 + dS1dI3 * dI3dc32
dS3dc31 = dS3dI1 * dI1dc31 + dS3dI2 * dI2dc31 + dS3dI3 * dI3dc31 dS1dc31 = dS1dI1 * dI1dc31 + dS1dI2 * dI2dc31 + dS1dI3 * dI3dc31
dS3dc13 = dS3dI1 * dI1dc13 + dS3dI2 * dI2dc13 + dS3dI3 * dI3dc13 dS1dc13 = dS1dI1 * dI1dc13 + dS1dI2 * dI2dc13 + dS1dI3 * dI3dc13
dS3dc44 = dS3dI1 * dI1dc44 + dS3dI2 * dI2dc44 + dS3dI3 * dI3dc44 dS1dc44 = dS1dI1 * dI1dc44 + dS1dI2 * dI2dc44 + dS1dI3 * dI3dc44
dS3dc55 = dS3dI1 * dI1dc55 + dS3dI2 * dI2dc55 + dS3dI3 * dI3dc55 dS1dc55 = dS1dI1 * dI1dc55 + dS1dI2 * dI2dc55 + dS1dI3 * dI3dc55
dS3dc66 = dS3dI1 * dI1dc66 + dS3dI2 * dI2dc66 + dS3dI3 * dI3dc66 dS1dc66 = dS1dI1 * dI1dc66 + dS1dI2 * dI2dc66 + dS1dI3 * dI3dc66
return dS1dc12, dS1dc21, dS1dc23, dS1dc32, dS1dc31, dS1dc13, dS1dc44, dS1dc55, dS1dc66, \ dS2dc12 = dS2dI1 * dI1dc12 + dS2dI2 * dI2dc12 + dS2dI3 * dI3dc12
dS2dc12, dS2dc21, dS2dc23, dS2dc32, dS2dc31, dS2dc13, dS2dc44, dS2dc55, dS2dc66, \ dS2dc21 = dS2dI1 * dI1dc21 + dS2dI2 * dI2dc21 + dS2dI3 * dI3dc21
dS3dc12, dS3dc21, dS3dc23, dS3dc32, dS3dc31, dS3dc13, dS3dc44, dS3dc55, dS3dc66 dS2dc23 = dS2dI1 * dI1dc23 + dS2dI2 * dI2dc23 + dS2dI3 * dI3dc23
dS2dc32 = dS2dI1 * dI1dc32 + dS2dI2 * dI2dc32 + dS2dI3 * dI3dc32
dS2dc31 = dS2dI1 * dI1dc31 + dS2dI2 * dI2dc31 + dS2dI3 * dI3dc31
dS2dc13 = dS2dI1 * dI1dc13 + dS2dI2 * dI2dc13 + dS2dI3 * dI3dc13
dS2dc44 = dS2dI1 * dI1dc44 + dS2dI2 * dI2dc44 + dS2dI3 * dI3dc44
dS2dc55 = dS2dI1 * dI1dc55 + dS2dI2 * dI2dc55 + dS2dI3 * dI3dc55
dS2dc66 = dS2dI1 * dI1dc66 + dS2dI2 * dI2dc66 + dS2dI3 * dI3dc66
dS3dc12 = dS3dI1 * dI1dc12 + dS3dI2 * dI2dc12 + dS3dI3 * dI3dc12
dS3dc21 = dS3dI1 * dI1dc21 + dS3dI2 * dI2dc21 + dS3dI3 * dI3dc21
dS3dc23 = dS3dI1 * dI1dc23 + dS3dI2 * dI2dc23 + dS3dI3 * dI3dc23
dS3dc32 = dS3dI1 * dI1dc32 + dS3dI2 * dI2dc32 + dS3dI3 * dI3dc32
dS3dc31 = dS3dI1 * dI1dc31 + dS3dI2 * dI2dc31 + dS3dI3 * dI3dc31
dS3dc13 = dS3dI1 * dI1dc13 + dS3dI2 * dI2dc13 + dS3dI3 * dI3dc13
dS3dc44 = dS3dI1 * dI1dc44 + dS3dI2 * dI2dc44 + dS3dI3 * dI3dc44
dS3dc55 = dS3dI1 * dI1dc55 + dS3dI2 * dI2dc55 + dS3dI3 * dI3dc55
dS3dc66 = dS3dI1 * dI1dc66 + dS3dI2 * dI2dc66 + dS3dI3 * dI3dc66
return dS1dc12, dS1dc21, dS1dc23, dS1dc32, dS1dc31, dS1dc13, dS1dc44, dS1dc55, dS1dc66, \
dS2dc12, dS2dc21, dS2dc23, dS2dc32, dS2dc31, dS2dc13, dS2dc44, dS2dc55, dS2dc66, \
dS3dc12, dS3dc21, dS3dc23, dS3dc32, dS3dc31, dS3dc13, dS3dc44, dS3dc55, dS3dc66
def Yld200418pBasis(sigma0, c12,c21,c23,c32,c31,c13,c44,c55,c66, def Yld200418pBasis(sigma0, c12,c21,c23,c32,c31,c13,c44,c55,c66,
d12,d21,d23,d32,d31,d13,d44,d55,d66, m, sigmas, Jac=False): d12,d21,d23,d32,d31,d13,d44,d55,d66, m, sigmas, Jac=False):
@ -739,6 +778,105 @@ def Yld200418pBasis(sigma0, c12,c21,c23,c32,c31,c13,c44,c55,c66,
jaco.append(jacv) jaco.append(jacv)
return np.array(jaco) return np.array(jaco)
def KarafillisBoyceBasis(sigma0, c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,
b1, b2, a, alpha , sigmas, Jac=False):
s1 = sigmas[0]; s2 = sigmas[1]; s3 = sigmas[2]
s4 = sigmas[3]; s5 = sigmas[4]; s6 = sigmas[5]
p = np.empty_like(sigmas); q = np.empty_like(sigmas)
p[0] = ( (c12+c13)*s1 - c13*s2 - c12*s3 )/3.0
p[1] = ( (c13+c11)*s2 - c13*s1 - c11*s3 )/3.0
p[2] = ( (c11+c12)*s3 - c12*s1 - c11*s2 )/3.0
p[3] = c14*s4
p[4] = c15*s5
p[5] = c16*s6
q[0] = ( (c22+c23)*s1 - c23*s2 - c22*s3 )/3.0
q[1] = ( (c23+c21)*s2 - c23*s1 - c21*s3 )/3.0
q[2] = ( (c21+c22)*s3 - c22*s1 - c21*s2 )/3.0
q[3] = c24*s4
q[4] = c25*s5
q[5] = c26*s6
plambdas, pInvariant = principalStress(p) # no sort
qlambdas, qInvariant = principalStress(q) # no sort
P1 = plambdas[0,:]; P2 = plambdas[1,:]; P3 = plambdas[2,:]
Q1 = qlambdas[0,:]; Q2 = qlambdas[1,:]; Q3 = qlambdas[2,:]
b1h = b1/2.0; b1h1 = b1h-1.0; b2h = b2/2.0; b2h1 = b2h-1.0
b1i = 1.0/b1; b2i = 1.0/b2
ai = 1.0/a
P2P3s = (P2-P3)**2; Q1s = Q1**2
P3P1s = (P3-P1)**2; Q2s = Q2**2
P1P2s = (P1-P2)**2; Q3s = Q3**2
phi10 = P2P3s**b1h + P3P1s**b1h + P1P2s**b1h
phi20 = Q1s**b2h+Q2s**b2h+Q3s**b2h; rb2 = 3.0**b2/(2.0**b2+2.0)
phi1 = (0.5*phi10)**b1i
phi2 = (rb2*phi20)**b2i
Stress = alpha*phi1**a + (1.0-alpha)*phi2**a; EqStress = Stress**ai
r = EqStress/sigma0 - 1.0
if not Jac:
return r.ravel()
else:
ln = lambda x : np.log(x + 1.0e-32)
drds = (r+1.0)*ai/Stress
drdphi1 = drds* alpha *a*phi1**(a-1.0)
drdphi2 = drds*(1.0-alpha)*a*phi2**(a-1.0)
dsda = alpha*phi1**a*ln(phi1) + (1.0-alpha)*phi2**a*ln(phi2)
dphi1dphi10 = phi1/phi10/b1; dphi2dphi20 = phi2/phi20/b2
dphi1dP1 = dphi1dphi10*b1*( P3P1s**b1h1*(P1-P3) + P1P2s**b1h1*(P1-P2))
dphi1dP2 = dphi1dphi10*b1*( P2P3s**b1h1*(P2-P3) + P1P2s**b1h1*(P2-P1))
dphi1dP3 = dphi1dphi10*b1*( P3P1s**b1h1*(P3-P1) + P2P3s**b1h1*(P3-P2))
dphi2dQ1 = dphi2dphi20*b2*Q1s*b2h1*Q1
dphi2dQ2 = dphi2dphi20*b2*Q2s*b2h1*Q2
dphi2dQ3 = dphi2dphi20*b2*Q3s*b2h1*Q3
dP1dc1,dP1dc2,dP1dc3,dP1dc4,dP1dc5,dP1dc6, \
dP2dc1,dP2dc2,dP2dc3,dP2dc4,dP2dc5,dP2dc6, \
dP3dc1,dP3dc2,dP3dc3,dP3dc4,dP3dc5,dP3dc6= \
principalStrs_Der(p, pInvariant, s1,s2,s3,s4,s5,s6, Karafillis=True)
dQ1dc1,dQ1dc2,dQ1dc3,dQ1dc4,dQ1dc5,dQ1dc6, \
dQ2dc1,dQ2dc2,dQ2dc3,dQ2dc4,dQ2dc5,dQ2dc6, \
dQ3dc1,dQ3dc2,dQ3dc3,dQ3dc4,dQ3dc5,dQ3dc6= \
principalStrs_Der(q, qInvariant, s1,s2,s3,s4,s5,s6, Karafillis=True)
dphi10db1 = ( (P2P3s**b1h)*ln(P2P3s)+(P3P1s**b1h)*ln(P3P1s)+(P1P2s**b1h)*ln(P1P2s) )*0.5
dphi20db2 = ( (P2P3s**b1h)*ln(P2P3s)+(P3P1s**b1h)*ln(P3P1s)+(P1P2s**b1h)*ln(P1P2s) )*0.5
drb2db2 = rb2*ln(3.0) - rb2*ln(2.0)/(1.0+2.0**(1.0-b2))
dphi1db1 = phi1*ln(phi10)*(-b1i*b1i) + b1i*phi1/(0.5*phi10)* 0.5*dphi10db1
dphi2db2 = phi2*ln(phi20)*(-b2i*b2i) + b2i*phi2/(rb2*phi20)*(rb2*dphi20db2 + drb2db2*phi20)
ja = drds*dsda - (r+1.0)*ln(Stress)/a/a #drda
jb1 = drds * ( alpha *a*phi1**(a-1)) * dphi1db1
jb2 = drds * ((1.0-alpha)*a*phi2**(a-1)) * dphi2db2
jalpha = drds * (phi1**a - phi2**a)
jc11 = drdphi1*( dphi1dP1*dP1dc1 + dphi1dP2*dP2dc1 + dphi1dP3*dP3dc1 )
jc12 = drdphi1*( dphi1dP1*dP1dc2 + dphi1dP2*dP2dc2 + dphi1dP3*dP3dc2 )
jc13 = drdphi1*( dphi1dP1*dP1dc3 + dphi1dP2*dP2dc3 + dphi1dP3*dP3dc3 )
jc14 = drdphi1*( dphi1dP1*dP1dc4 + dphi1dP2*dP2dc4 + dphi1dP3*dP3dc4 )
jc15 = drdphi1*( dphi1dP1*dP1dc5 + dphi1dP2*dP2dc5 + dphi1dP3*dP3dc5 )
jc16 = drdphi1*( dphi1dP1*dP1dc6 + dphi1dP2*dP2dc6 + dphi1dP3*dP3dc6 )
jc21 = drdphi2*( dphi2dQ1*dQ1dc1 + dphi2dQ2*dQ2dc1 + dphi2dQ3*dQ3dc1 )
jc22 = drdphi2*( dphi2dQ1*dQ1dc2 + dphi2dQ2*dQ2dc2 + dphi2dQ3*dQ3dc2 )
jc23 = drdphi2*( dphi2dQ1*dQ1dc3 + dphi2dQ2*dQ2dc3 + dphi2dQ3*dQ3dc3 )
jc24 = drdphi2*( dphi2dQ1*dQ1dc4 + dphi2dQ2*dQ2dc4 + dphi2dQ3*dQ3dc4 )
jc25 = drdphi2*( dphi2dQ1*dQ1dc5 + dphi2dQ2*dQ2dc5 + dphi2dQ3*dQ3dc5 )
jc26 = drdphi2*( dphi2dQ1*dQ1dc6 + dphi2dQ2*dQ2dc6 + dphi2dQ3*dQ3dc6 )
jaco = []
for jacv in zip(jc11,jc12,jc13,jc14,jc15,jc16,jc21,jc22,jc23,jc24,jc25,jc26,
jb1,jb2,ja,jalpha):
jaco.append(jacv)
return np.array(jaco)
fittingCriteria = { fittingCriteria = {
'tresca' :{'func' : Tresca, 'tresca' :{'func' : Tresca,
@ -835,6 +973,14 @@ fittingCriteria = {
\n Y, c12,c21,c23,c32,c31,c13,c44,c55,c66,d12,d21,d23,d32,d31,d13,d44,d55,d66,m\n', \n Y, c12,c21,c23,c32,c31,c13,c44,c55,c66,d12,d21,d23,d32,d31,d13,d44,d55,d66,m\n',
'error': 'The standard deviation errors are: ' 'error': 'The standard deviation errors are: '
}, },
'karafillis' :{'func' : KarafillisBoyce,
'num' : 16,'err':np.inf,
'name' : 'Yld200418p',
'paras': 'c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,b1,b2,a,alpha',
'text' : '\nCoefficients of Karafillis-Boyce yield criterion: \
\n c11,c12,c13,c14,c15,c16,c21,c22,c23,c24,c25,c26,b1,b2,a,alpha\n',
'error': 'The standard deviation errors are: '
},
'worst' :{'err':np.inf}, 'worst' :{'err':np.inf},
'best' :{'err':np.inf} 'best' :{'err':np.inf}
} }
@ -917,7 +1063,8 @@ class Criterion(object):
popt, pcov, infodict, errmsg, ierr = \ popt, pcov, infodict, errmsg, ierr = \
leastsqBound (criteria.fun, initialguess, args=(ydata,stress), leastsqBound (criteria.fun, initialguess, args=(ydata,stress),
bounds=bounds, Dfun=criteria.jac, full_output=True) bounds=bounds, Dfun=criteria.jac, full_output=True)
if ierr not in [1, 2, 3, 4]: raise RuntimeError("Optimal parameters not found: " + errmsg) if ierr not in [1, 2, 3, 4]:
raise RuntimeError("Optimal parameters not found: " + errmsg)
if (len(ydata) > len(initialguess)) and pcov is not None: if (len(ydata) > len(initialguess)) and pcov is not None:
s_sq = (criteria.fun(popt, *(ydata,stress))**2).sum()/(len(ydata)-len(initialguess)) s_sq = (criteria.fun(popt, *(ydata,stress))**2).sum()/(len(ydata)-len(initialguess))
pcov = pcov * s_sq pcov = pcov * s_sq