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

some comments for discussion

This commit is contained in:
Martin Diehl 2015-04-07 07:42:46 +00:00
parent ab89c84ec8
commit 6e650a42b4
1 changed files with 13 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
processing/misc/yieldSurface.py
View File

@ -42,7 +42,6 @@ def principalStress(p):
sin = np.sin; cos = np.cos sin = np.sin; cos = np.cos
I1,I2,I3 = invariant(p) I1,I2,I3 = invariant(p)
third = 1.0/3.0
I1s3I2= (I1**2 - 3.0*I2)**0.5 I1s3I2= (I1**2 - 3.0*I2)**0.5
numer = 2.0*I1**3 - 9.0*I1*I2 + 27.0*I3 numer = 2.0*I1**3 - 9.0*I1*I2 + 27.0*I3
denom = I1s3I2**(-3.0) denom = I1s3I2**(-3.0)
@ -311,8 +310,12 @@ def Drucker(eqStress, paras, sigmas, mFix, criteria, Jac = False):
if mFix[0]: p = mFix[1] if mFix[0]: p = mFix[1]
else: p = paras[-1] else: p = paras[-1]
I1,I2,I3 = invariant(sigmas) I1,I2,I3 = invariant(sigmas)
#I = invariant(sigmas)
#J = np.zeros([3])
J2 = I1**2/3.0 - I2 J2 = I1**2/3.0 - I2
#J[1] = I[0]**2/3.0 - I[1]
J3 = I1**3/13.5 - I1*I2/3.0 + I3 J3 = I1**3/13.5 - I1*I2/3.0 + I3
#J[2] = I[0]**3/13.5 - I[0]*I[1]/3.0 + I[2] etc.
J2_3p = J2**(3.0*p) J2_3p = J2**(3.0*p)
J3_2p = J3**(2.0*p) J3_2p = J3**(2.0*p)
left = J2_3p - C_D*J3_2p left = J2_3p - C_D*J3_2p
@ -355,7 +358,9 @@ def Hill1979(eqStress,paras, sigmas, mFix, criteria, Jac = False):
coeff = paras[0:6] coeff = paras[0:6]
s1,s2,s3 = principalStresses(sigmas) s1,s2,s3 = principalStresses(sigmas)
# s= principalStresses(sigmas)
diffs = np.array([s2-s3, s3-s1, s1-s2, 2.0*s1-s2-s3, 2.0*s2-s3-s1, 2.0*s3-s1-s2])**2 diffs = np.array([s2-s3, s3-s1, s1-s2, 2.0*s1-s2-s3, 2.0*s2-s3-s1, 2.0*s3-s1-s2])**2
#diffs = np.array([s[1]-s[2], s[2]-s[0], etc ... s1-s2, 2.0*s1-s2-s3, 2.0*s2-s3-s1, 2.0*s3-s1-s2])**2
diffsm = diffs**(m/2.0) diffsm = diffs**(m/2.0)
base = np.dot(coeff,diffsm) base = np.dot(coeff,diffsm)
r = base**(1.0/m)/eqStress #left = base**mi r = base**(1.0/m)/eqStress #left = base**mi
@ -424,7 +429,7 @@ def Barlat1989(eqStress, paras, sigmas, mFix, criteria, Jac=False):
Anisotropic: a, c, h, p, m; m is optional Anisotropic: a, c, h, p, m; m is optional
''' '''
a, c, h, p = paras[0:4] a, c, h, p = paras[0:4]
if mFix[0]: m = mFix[1] if mFix[0]: m = mFix[1] #???
else: m = paras[-1] else: m = paras[-1]
s11 = sigmas[0]; s22 = sigmas[1]; s12 = sigmas[3] s11 = sigmas[0]; s22 = sigmas[1]; s12 = sigmas[3]
@ -796,6 +801,7 @@ def KarafillisBoyce(eqStress, paras, sigmas, mFix, criteria, Jac=False):
jb1 = dphi1db1*(drds*a*phi1**(a-1)*alpha ) jb1 = dphi1db1*(drds*a*phi1**(a-1)*alpha )
jb2 = dphi2db2*(drds*a*phi2**(a-1)*(1.0-alpha)) jb2 = dphi2db2*(drds*a*phi2**(a-1)*(1.0-alpha))
jc1 = np.sum([dphi1dP[i]*dPdc[i] for i in xrange(3)],axis=0)*drds*a*phi1**(a-1.0)*alpha jc1 = np.sum([dphi1dP[i]*dPdc[i] for i in xrange(3)],axis=0)*drds*a*phi1**(a-1.0)*alpha
#jc1 = np.sum([dphi1dP[0:3]*dPdc[0:3])*drds*a*phi1**(a-1.0)*alpha
jc2 = np.sum([dphi2dQ[i]*dQdc[i] for i in xrange(3)],axis=0)*drds*a*phi2**(a-1.0)*(1.0-alpha) jc2 = np.sum([dphi2dQ[i]*dQdc[i] for i in xrange(3)],axis=0)*drds*a*phi2**(a-1.0)*(1.0-alpha)
jalpha = drds * (phi1**a - phi2**a) jalpha = drds * (phi1**a - phi2**a)
@ -1202,8 +1208,10 @@ def doSim(delay,thread):
upper,lower = table.data[line,9],table.data[line-1,9] # values for linear interpolation upper,lower = table.data[line,9],table.data[line-1,9] # values for linear interpolation
stress = np.array(table.data[line-1,0:9] * (upper-threshold)/(upper-lower) + \ stress = np.array(table.data[line-1,0:9] * (upper-threshold)/(upper-lower) + \
table.data[line ,0:9] * (threshold-lower)/(upper-lower)).reshape(3,3) # linear interpolation of stress values table.data[line ,0:9] * (threshold-lower)/(upper-lower)).reshape(3,3) # linear interpolation of stress values
#stress = 0.5*(stress+stress.T) # symmetrise
#for the mapping, a fuction from DAMASK (33to6) simplifies
dstrain= np.array(table.data[line,10:] - table.data[line-1,10:]).reshape(3,3) dstrain= np.array(table.data[line,10:] - table.data[line-1,10:]).reshape(3,3)
#
yieldStress[i,0]= stress[0,0]; yieldStress[i,1]=stress[1,1]; yieldStress[i,2]=stress[2,2] yieldStress[i,0]= stress[0,0]; yieldStress[i,1]=stress[1,1]; yieldStress[i,2]=stress[2,2]
yieldStress[i,3]=(stress[0,1] + stress[1,0])/2.0 # 0 3 5 yieldStress[i,3]=(stress[0,1] + stress[1,0])/2.0 # 0 3 5
yieldStress[i,4]=(stress[1,2] + stress[2,1])/2.0 # * 1 4 yieldStress yieldStress[i,4]=(stress[1,2] + stress[2,1])/2.0 # * 1 4 yieldStress
@ -1239,7 +1247,7 @@ def loadcaseNo():
N_simulations+=1 N_simulations+=1
return N_simulations return N_simulations
def converged(): def converged(): # is there any meaningfull stopping criterion?
global N_simulations global N_simulations
if N_simulations < options.max: if N_simulations < options.max:
return False return False
@ -1262,6 +1270,7 @@ parser.add_option('-g','--geometry', dest='geometry', type='string',
help='name of the geometry file [%default]', metavar='string') help='name of the geometry file [%default]', metavar='string')
parser.add_option('-c','--criterion', dest='criterion', choices=fitCriteria.keys(), parser.add_option('-c','--criterion', dest='criterion', choices=fitCriteria.keys(),
help='criterion for stopping simulations [%default]', metavar='string') help='criterion for stopping simulations [%default]', metavar='string')
# best/worse fitting? Stopping?
parser.add_option('-f','--fitting', dest='fitting', choices=thresholdParameter, parser.add_option('-f','--fitting', dest='fitting', choices=thresholdParameter,
help='yield criterion [%default]', metavar='string') help='yield criterion [%default]', metavar='string')
parser.add_option('-y','--yieldvalue', dest='yieldValue', type='float', nargs=3, parser.add_option('-y','--yieldvalue', dest='yieldValue', type='float', nargs=3,