1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
| # calculate the mean, alternating and equivalent stresses
# and the service life and the damage
def calcul(k,f):
maxStress=0
for idx, value in enumerate(Stress) :
S11 = value.data[0]
S22 = value.data[1]
S12 = value.data[3]
e = Stress1[idx]
S_11 = e.data[0]
S_22 = e.data[1]
S_12 = e.data[3]
# calculate the alternating stress tensor and its eigenvalues
D11 = (S11-S_11)/2
D22 = (S22-S_22)/2
D12 = (S12-S_12)/2
Dyn1 = ((D11+D22)-sqrt(D11**2+D22**2+4*D12**2-2*D11*D22))/2
Dyn2 = ((D11+D22)+sqrt(D11**2+D22**2+4*D12**2-2*D11*D22))/2
# calculate the mean stress tensor and its eigenvalues
St11 = (S11+S_11)/2
St22 = (S22+S_22)/2
St12 = (S12+S_12)/2
Stat1 = ((St11+St22)+sqrt(St11**2+St22**2+4*St12**2-2*St11*St22))/2
Stat2 = ((St11+St22)-sqrt(St11**2+St22**2+4*St12**2-2*St11*St22))/2
# calculate the alternating stress and the mean stress
if abs(Dyn1)>=abs(Dyn2):
sig_dyn = abs(Dyn1)
else:
sig_dyn = abs(Dyn2)
if abs(Stat1)>=abs(Stat2):
signe = Stat1/abs(Stat1)
sig_stat = signe*sqrt(Stat1**2-Stat1*Stat2+Stat2**2)
else:
signe = Stat2/abs(Stat2)
sig_stat = signe*sqrt(Stat2**2-Stat1*Stat2+Stat1**2)
# calculate the maximum and minimum stress and the ratio
Max = sig_dyn+sig_stat
Min = sig_stat-sig_dyn
R = Min/Max
# calculate the equivalent stress
alpha_G = abs(sig_D_mean/((sig_dyn/sig_stat)+(sig_D_mean/Rm)))*\
((1+(sig_dyn/sig_stat)**2)/(sig_dyn**2+sig_stat**2))**0.5
alpha_S = abs(sig_D_mean/((sig_dyn/sig_stat)+(sig_D_mean/Re)))*\
((1+(sig_dyn/sig_stat)**2)/(sig_dyn**2+sig_stat**2))**0.5
if choix == 'Goodman':
sig_eq = sig_Goodman(Rm, sig_dyn, sig_stat, alpha_G)
elif choix == 'Gerber':
sig_eq = sig_Gerber(Re, sig_dyn, sig_stat)
elif choix == 'Soderberg':
sig_eq = sig_Soderberg(Re, sig_dyn, sig_stat, alpha_S)
else:
print ('Votre choix n appartient pas a la liste precedante')
if sig_eq > maxStress:
maxStress = sig_eq
# determination of the service life
if sig_eq > y[0]:
N_cycle = 0
elif sig_eq <= y[-1]:
N_cycle = 1.0E+10
else:
i=0
len_y=len(y)
while (i < len_y - 1):
if (sig_eq < y[i]) and (sig_eq > y[i+1]):
N_cycle = x[i+1]-((y[i+1]-sig_eq)*(x[i+1]-x[i])/(y[i+1]-y[i]))
break
i=i+1
# calculate the damage
Damage = T_cycle[k-1]/N_cycle
# projection of the alternating stresses on the DTA and FTA curves
if (R<-1) or (R>1):
propa_DTA1=-1.0E-01
propa_DTA2=-1.0E-01
propa_FTA=-1.0E-01
else:
j =0
len_DTA1 = len(DTA1)
while (j < len_DTA1 - 1):
if (R > Ratio[j]) and (R < Ratio[j+1]):
sig_DTA1 = DTA1[j+1]-((Ratio[j+1]-R)*(DTA1[j+1]-DTA1[j])/(Ratio[j+1]-Ratio[j]))
propa_DTA1 = sig_dyn-sig_DTA1
sig_DTA2 = DTA2[j+1]-((Ratio[j+1]-R)*(DTA2[j+1]-DTA2[j])/(Ratio[j+1]-Ratio[j]))
propa_DTA2 = sig_dyn-sig_DTA2
sig_FTA = FTA[j+1]-((Ratio[j+1]-R)*(FTA[j+1]-FTA[j])/(Ratio[j+1]-Ratio[j]))
propa_FTA = sig_dyn-sig_FTA
break
j=j+1
f.write('Instance:\t %s\tElement:\t %6d\tSig_stat:\t %E\tSig_dyn:\t %E\
\tMax:\t %E\tMin:\t %E\tR:\t %E\tSig_eq:\t %E\tN_cycle:\t %E\
\tMaxStress:\t %E\tDamage:\t %E\tpropa_DTA1:\t %E\tpropa_DTA2:\t %E\
\tpropa_FTA:\t %E\n'%(value.instance.name,\
value.elementLabel,sig_stat,\
sig_dyn,Max,Min,R,sig_eq,N_cycle,
maxStress,Damage,propa_DTA1,propa_DTA2,propa_FTA))
return maxStress
#-----------------------------------------------------------GAG-----------------------------------------------------------------
# determine the stress value at the frame number 2 and 11 of the second step
# in the "SKIN" element set (GAG)
stressField = lcf.steps['Load'].frames[2].fieldOutputs['S']
stressField1 = lcf.steps['Load'].frames[11].fieldOutputs['S']
topCenter = lcf.rootAssembly.instances['V20110406-1-1'].elementSets['SKIN']
Stress = stressField.getSubset(region=topCenter, position=CENTROID).values
Stress1 = stressField1.getSubset(region=topCenter, position=CENTROID).values
maxStresses = []
maxStresses.append(calcul(nb_loadcase+1, f10))
# etc. for each call to calcul()
|
Partager