In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import weibull_min, norm
from scipy.stats import skew, kurtosis
import pathlib
import subprocess
from IPython.display import FileLink
#os.chdir('/home/kelidas/workspace_git/python_examples/atena_2')
In [2]:
%config InlineBackend.figure_format = 'svg'
%matplotlib inline
Postprocessing výsledků z numerických výpočtů¶
ve složce
C_AtenaResultsjsou soubory C_1.txt až C_64.txt, které obsahují hodnoty z čtyř monitorů zaznamenaných programem Atena (monitory jsou uvozeny hlavičkou a jsou za sebou)velikostí bylo cca 10, simulací (pilotní studie 32 s opakováním dokud se nenašly optimální parametry, finální 64) - tzn. toto vyhodnocení je nutné zopakovat ještě mnohokrát (např. 20-30x)
bohužel nic lepšího než excel jsem neuměl, takže to vypadalo viz
C_AtenaResults\C_W64.xls+ trochu mně pomáhalo C a gnuplotcílem je vykreslit graf všech křivek
posun vs sílase zvýrazněním maxim, kdeposun = (monitor1 - monitor2) * 1e6 [mm]asíla = monitor3najít maximální sílu všech křivek a spočítat střední hodnotu, smodch, šikmost a špičatost
vykreslit histogramy maximální síly a jí příslušející deformace
fit maximální síly pomocí weibullova a normálního rozdělení
In [3]:
data_dir = 'C_AtenaResults'
data_all = []
for i in range(1,65):
c = 0
data = []
with open(os.path.join(data_dir, 'C_%d.txt' % i), 'r', encoding='cp1250') as f:
monitor = np.array([])
for l in f.readlines():
if l == '-----------------\n':
c += 1
continue
if c == 2:
if l == '\n':
c = 0
data.append(monitor)
monitor = np.array([])
else:
v = np.float64(l.split()[1])
monitor = np.append(monitor, v)
data_all.append(np.array(data).T)
In [4]:
data_all[0][:3, :]
Out[4]:
Vykreslení ld křivek se zvýrazněním maxim¶
In [5]:
fig, ax = plt.subplots(figsize=(6,4), dpi=100)
for data in data_all:
displ = (data[:, 0] - data[:, 1]) * 1e6
force = data[:, 3]
ax.plot(displ, force, 'k-', lw=.5)
fidx = np.argmax(force)
ax.plot(displ[fidx], force[fidx], 'ro', lw=.5, ms=2)
ax.set_xlim(0,20)
ax.set_xlabel('displ [mm]')
ax.set_ylabel('force [kN]');
Analýza maxim a jim odpovídajících posunů¶
In [6]:
f_max = np.array([])
d_fmax = np.array([])
for data in data_all:
f = data[:,3]
fm_idx = np.argmax(f)
f_max = np.append(f_max, f[fm_idx])
displ = (data[:, 0] - data[:, 1]) * 1e6
d_fmax = np.append(d_fmax, displ[fm_idx])
In [7]:
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8,4), dpi=100)
pn = norm.fit(f_max)
pw = weibull_min.fit(f_max, floc=0)
rv_norm = norm(*pn)
rv_weib = weibull_min(*pw)
print('parametry normálního rozdělení (mu, var) =', rv_norm.stats('mv'))
print('parametry dvouparametrického weibullova rozdělení (mu, var) =', rv_weib.stats('mv'))
ax1.hist(f_max, density=True)
x = np.linspace(1.5,3.5,1000)
ax1.plot(x, rv_norm.pdf(x), label='norm')
ax1.plot(x, rv_weib.pdf(x), label='weibull')
ax1.set_xlabel(u'maximální síla [kN]')
ax1.set_ylabel(u'relativní četnosti')
ax1.legend(loc='best')
ax2.hist(d_fmax, density=True)
ax2.set_xlabel('deformace [mm]')
fig.tight_layout()
In [8]:
print('{:40} = {:+.6f}'.format(u'střední hodnota maximální síly', np.mean(f_max)))
print('{:40} = {:+.6f}'.format(u'směrodatná odchylka maximální síly', np.std(f_max)))
print('{:40} = {:+.6f}'.format(u'šikmost maximální síly', skew(f_max)))
print('{:40} = {:+.6f}'.format(u'špičatost maximální síly', kurtosis(f_max)))
Pro každý vzorek vygenerovat A4¶
- obsahuje graf
- hodnotu max síly a odpovídajícího posunu
In [9]:
path = pathlib.Path('report') / 'figs'
path.mkdir(parents=True, exist_ok=True)
In [10]:
report = r"""\documentclass[a4paper,10pt,titlepage]{article}
\usepackage[left=2cm,
textwidth=17cm,
textheight=24cm,
top=3cm]{geometry}
\usepackage[utf8]{inputenc}
\usepackage[czech]{babel}
%\usepackage{times}
\usepackage{lmodern}
\usepackage[unicode]{hyperref}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{verbatim}
%\usepackage{layouts}
\title{Report}
\date{\today}
\author{Václav Sadílek}
\maxdeadcycles=200
\begin{document}
\maketitle
\listoffigures
\section{Vzorky}
"""
In [11]:
fig, ax = plt.subplots(figsize=(6,4), tight_layout=True)
for di, data in enumerate(data_all):
displ = (data[:, 0] - data[:, 1]) * 1e6
force = data[:, 3]
fidx = np.argmax(force)
if di == 0:
l1, = ax.plot(displ, force, 'k-', lw=1)
p1, = ax.plot(displ[fidx], force[fidx], 'ro', lw=1, ms=2)
ax.set_xlim(0,20)
ax.set_ylim(0, 3.6)
ax.set_xlabel('displ [mm]')
ax.set_ylabel('force [kN]')
else:
l1.set_xdata(displ)
l1.set_ydata(force)
p1.set_xdata(displ[fidx])
p1.set_ydata(force[fidx])
fig.savefig('report/figs/fig_{:02d}.pdf'.format(di), bbox_inches='tight', pad_inches=.01)
report += r"""
\subsection{{Vzorek č.{}}}
\begin{{figure}}[htb!]
\centering
\includegraphics[width=10cm]{{figs/fig_{:02d}.pdf}}
\caption{{Vzorek č.~{} - maximální síla {:.3f}\,kN a příslušný posun {:.3f}\,mm}}
\label{{fig:vzorek{:02d}}}
\end{{figure}}
\clearpage
""".format(di+1, di, di+1, force[fidx], displ[fidx], di+1)
In [12]:
report += '\end{document}'
fname = 'report/report.tex'
with open(fname, 'w') as f:
f.write(report)
In [13]:
#os.system('pdflatex -output-directory report report/report.tex')
subprocess.Popen(['pdflatex', '-output-directory', 'report', 'report/report.tex'])
FileLink(fname.replace('tex', 'pdf'))
Out[13]: