1
# -*- coding: utf-8 -*-
3
Created on Thu Sep 3 23:17:27 2020
7
plot stress vs. strain with Bridgeman correction
12
import matplotlib.pyplot as plt
13
from matplotlib.widgets import Button
15
from scipy import interpolate
16
from scipy.optimize import curve_fit
17
import bottleneck as bn
21
cfg = configparser.ConfigParser()
22
cfg.read("constrain.ini")
23
notched_specimen = cfg.getboolean("DEFAULT", "notched_specimen")
24
factor = cfg.getint("DEFAULT", "nsmooth") # moving average filter -- window width for strain data
25
D0 = cfg.getfloat("DEFAULT", "D0") # initial diameter in cylindric gauge section
26
case = cfg.get("DEFAULT", "case")
27
dt_strain = cfg.getfloat("DEFAULT", "dt_images") # strain acquisition time period
28
offset_strain = cfg.getint("DEFAULT", "offset_images") # number of images which were skipped after trigger point
29
nueps = cfg.getfloat("DEFAULT", "nueps") # plastic Poisson ratio
31
def find_limits(t1, t2):
32
""" return common time range between to time axes
35
low = max(t1[0], t2[0])
36
high = min(t1[-1], t2[-1])
40
def power_law_plasticity(strain):
41
""" return power law plastic strain according to LS-Dyna material model Nr. 18"
47
eps_ela = (E/K)**(1/(n-1))
48
print("elastic strain at yield is: ", eps_ela)
49
sigma = K * (eps_ela + strain)**n
50
return eps_ela + strain, sigma * 1000
52
def load_force_LSDyna():
53
""" read force from LS-Dyna axis-symmetric analysis
55
force_filename = "force.csv"
56
data = np.genfromtxt(force_filename, delimiter=",", skip_header=2)
57
time_force = data[:,0]
58
force = data[:,1] * 2 * np.pi * 1000
59
dt_force = time_force[1] - time_force[0]
60
return time_force, force, dt_force
62
def load_force_Zwick():
63
""" read force from Zwick UTM
65
force_filename = "force.TRA"
66
data = np.genfromtxt(force_filename, delimiter=";", skip_header=16)
67
time_force = data[:,0]
69
dt_force = time_force[1] - time_force[0]
70
return time_force, force, dt_force
72
def load_force_SHTB():
73
""" read force from SHTB data analysis
75
data = np.genfromtxt("time_force_nonshifted.txt", skip_header=0)
76
time_force = data[:,0]
77
force = data[:,3] * 1000 # convert kN to N
78
dt_force = time_force[1] - time_force[0]
80
#plt.plot(time_force, force)
82
return time_force, force, dt_force
85
""" load strain and triaxiality data from conStrain analysis """
87
strain_data = np.genfromtxt("strain.txt")
88
indices_strain = strain_data[:,0]
89
dias = strain_data[:,2]
90
notch_radii = strain_data[:,3]
93
eps_eq = (1./nueps) * np.log(dias[0] / dias)
94
eps_eq = bn.move_mean(eps_eq,window=factor,min_count=1)
95
dias = bn.move_mean(dias,window=factor,min_count=1)
97
if notched_specimen == False:
98
notch_radii = np.where(eps_eq < 0.25, 100*dias, notch_radii)
100
notch_radii = bn.move_mean(notch_radii,window=factor,min_count=1)
101
b = triaxiality.bridgman_function(dias, notch_radii)
102
triaxs = triaxiality.triax_function(dias, notch_radii)
105
return indices_strain, eps_eq, dias, b, triaxs
107
def shift_conStrain(indices_strain, nFrames, dt_strain):
108
""" shift conStrain data by a number of frames to make it coincident in time with the force trigger time """
109
time_strain = (indices_strain + nFrames) * dt_strain
112
time_strain += 120.0/5090 # account for strain gauge being 120 mm positioned after center of specimen
115
def interpolate_to_common_time(time_strain):
116
# construct new time axis
117
dt = min(dt_force, dt_strain)
118
low, high = find_limits(time_force, time_strain)
119
print("common time limits: %g -- %g" % (low, high))
120
n = int((high-low)/dt)
121
time = np.linspace(low, high, num=n, endpoint=True)
122
print("new time axis:", time[0], time[-1])
124
# interpolate strain-based quantities onto new time axis
125
f = interpolate.interp1d(time_strain, eps_eq0)
127
f = interpolate.interp1d(time_strain, b0)
129
f = interpolate.interp1d(time_strain, triaxs0)
131
f = interpolate.interp1d(time_strain, dias0)
134
# interpolate force quantities onto new time axis
135
f = interpolate.interp1d(time_force, force0)
138
return eps_eq_i, diameter_i, force_i, b_i, triax_i, time
142
def compute_stress_measures(force, diameter, b):
143
# nominal (engineering) stress
144
area0 = 0.25 * np.pi * D0**2
145
stress_nominal = force/area0
147
# stress measures based on current area
148
cd_mm = D0 * diameter / diameter[0]# current diameter in mm
149
area = 0.25 * np.pi * cd_mm**2
150
stress_cauchy = force / area
151
stress_eq = stress_cauchy * b
153
return stress_nominal, stress_cauchy, stress_eq
158
time_force, force0, dt_force = load_force_Zwick()
160
time_force, force0, dt_force = load_force_SHTB()
161
elif case == "LS-Dyna":
162
time_force, force0, dt_force = load_force_LSDyna()
164
print("unknown force case, exiting")
167
print("force time axis covers %g to %g" % (time_force[0], time_force[-1]))
168
indices_strain, eps_eq0, dias0, b0, triaxs0 = load_conStrain()
169
time_strain = shift_conStrain(indices_strain, offset_strain, dt_strain)
170
print("strain time axis covers %g to %g" % (time_strain[0], time_strain[-1]))
171
eps_eq, diameter, force, b, triax, time = interpolate_to_common_time(time_strain)
173
stress_nominal, stress_cauchy, stress_eq = compute_stress_measures(force, diameter, b)
175
# first plot window: nominal stress and strain over time
176
fig, ax = plt.subplots(3, figsize=(8,7))
177
ax[0].plot(time, stress_nominal, color="blue", label="stress")
179
ax[0].set_ylabel("nominal stress [MPa]", color="blue")
180
ax[0].set_xlabel("time [ms]")
181
ax[0].set_ylim((0,None))
183
ax2.set_ylabel("true strain [-]", color="red")
184
ax2.plot(time, eps_eq, "r-", label="strain")
185
ax2.set_ylim((0,None))
190
# third plot: triaxiality over time
191
ax[1].plot(time, triax, label="triaxiality", color="blue")
193
ax[1].set_ylabel("triaxiality [-]", color="blue")
194
ax[1].set_xlabel("time")
195
ax[1].set_ylim((0, 1.2))
197
ax32.set_ylabel("bridgman factor [-]", color="red")
198
ax32.plot(time, b, "r-", label=b)
199
ax32.set_ylim((0,1.1))
202
# 3rd plot: stress / strain
203
l_stress_nominal, = ax[2].plot(eps_eq, stress_nominal, lw=2, label="nominal stress")
204
l_stress_cauchy, = ax[2].plot(eps_eq, stress_cauchy, lw=2, label="Cauchy stress")
205
l_stress_eq, = ax[2].plot(eps_eq, stress_eq, "g--", lw=2, label="true eq. stress, Bridgman corrected")
207
if case == "LS-Dyna":
208
x, y = power_law_plasticity(eps_eq)
209
ax[2].plot(x, y, "b--", lw=2, label="material model stress")
212
ax[2].set_xlabel("equivalent strain [-]")
213
ax[2].set_ylabel("stress [MPa]")
215
ax[2].set_xlim(0, None)
216
ax[2].set_ylim(0, None)
225
self.update_stresses()
227
def update_stresses(self):
228
offset = offset_strain + self.ind
229
time_strain = shift_conStrain(indices_strain, offset, dt_strain)
230
self.eps_eq, self.diameter, self.force, self.b, self.triax, time = interpolate_to_common_time(time_strain)
231
self.stress_nominal, self.stress_cauchy, self.stress_eq = compute_stress_measures(self.force, self.diameter, self.b)
233
l_stress_nominal.set_xdata(self.eps_eq)
234
l_stress_nominal.set_ydata(self.stress_nominal)
235
l_stress_cauchy.set_xdata(self.eps_eq)
236
l_stress_cauchy.set_ydata(self.stress_cauchy)
237
l_stress_eq.set_xdata(self.eps_eq)
238
l_stress_eq.set_ydata(self.stress_eq)
239
ax[2].set_title("frame offset: %d, @failure: triax=%3.2f, stress=%4d, strain=%3.2f" % (offset, self.triax[-1], self.stress_eq[-1], self.eps_eq[-1]))
242
def next(self, event):
244
self.update_stresses()
246
def prev(self, event):
248
self.update_stresses()
251
axprev = plt.axes([0.2, 0.01, 0.2, 0.04])
252
axnext = plt.axes([0.6, 0.01, 0.2, 0.04])
253
bnext = Button(axnext, '+ 1 frame')
254
bnext.on_clicked(callback.next)
255
bprev = Button(axprev, '- 1 frame')
256
bprev.on_clicked(callback.prev)
258
plt.subplots_adjust(bottom=0.15, top=1)
262
data = np.column_stack((callback.eps_eq, callback.stress_nominal, callback.stress_cauchy, callback.stress_eq, callback.triax, callback.b))
263
np.savetxt("stress_strain.dat", data, header="strain, nominal_stress, Cauchy_stress, Cauchy_Bridgman_stress, triaxiality, bridgman_factor")
268
#-----------------------
269
#fit Hockett-Sherby plasticity model to Cauchy-Bridgman stress
275
bounds = [(500, 25000),
279
bounds = np.asarray(bounds)
282
p0 = (sig0, sig1, m, n)
283
indices = np.asarray(np.where(eps_eq > 0.015)).ravel()
284
print("indices:", indices)
290
print("start index: ", idx)
291
print("eps_eq:", eps_eq)
293
eps_pl = eps_eq[idx:]
295
popt, pcov = curve_fit(hocket_sherby, eps_pl, stress_eq[idx:], p0=p0, bounds=bounds)
296
fit = hocket_sherby(eps_pl, *popt)
297
format = "%10s %6s %6s %6s %6s"
298
format_num = "%10s %06.2f %06.2f %06.2f %06.2f"
299
print(format % ("parameters", "sig0", "sig1", "m", "n"))
300
print(format_num % ("start", *p0))
301
print(format_num % ("final", *popt))
303
plt.plot(eps_pl, stress_eq[idx:])
304
plt.plot(eps_pl, fit)
308
# # fit Swift-Voce plasticity model to Cauchy-Bridgman stress
317
# bounds = [(1000, 2000),
324
# bounds = np.asarray(bounds)
327
# p0 = (A, eps0, n, k0, Q, b, alpha)
328
# idx = np.where(eps_eq > 0.0025)[0][0]
329
# eps_pl = eps_eq[idx:]
330
# print("start index: ", idx)
331
# popt, pcov = curve_fit(swift_voce, eps_pl, stress_eq[idx:], p0=p0, bounds=bounds)
332
# fit = swift_voce(eps_pl, *popt)
333
# format = "%10s %6s %6s %6s %6s %6s %6s %6s"
334
# format_num = "%10s %06.2f %06.2f %06.2f %06.2f %06.2f %06.2f %06.2f"
335
# print(format % ("parameters", "A", "eps0", "n", "k0", "Q", "b", "alpha"))
336
# print(format_num % ("start", *p0))
337
#print(format_num % ("final", *popt))
added
removed
align_stress_strain.py
