atena

Lang cs
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
#os.chdir('/home/kelidas/workspace_git/python_examples/atena_2')
%config InlineBackend.figure_format = 'svg'
%matplotlib inline

Postprocessing výsledků z numerických výpočtů

  • ve složce C_AtenaResults jsou 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 gnuplot

  • cílem je vykreslit graf všech křivek posun vs síla se zvýrazněním maxim, kde posun = (monitor1 - monitor2) * 1e6 [mm] a síla = monitor3

  • nají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 [2]:
data_dir = ''
data_all = []
for i in range(1,2):
    c = 0
    data = []
    with open(os.path.join(data_dir, 'C_%d.txt' % i), encoding='cp1250') as f:
        monitor = np.array([])
        for line in f.readlines():
            l = line
            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)  
data_all
Out[2]:
[array([[5.600e-08, 2.400e-08, 1.343e-04, 1.119e-02],
        [1.120e-07, 4.800e-08, 2.686e-04, 2.238e-02],
        [1.680e-07, 7.201e-08, 4.029e-04, 3.357e-02],
        [2.240e-07, 9.601e-08, 5.372e-04, 4.476e-02],
        [2.800e-07, 1.200e-07, 6.715e-04, 5.595e-02],
        [3.360e-07, 1.440e-07, 8.058e-04, 6.715e-02],
        [3.920e-07, 1.680e-07, 9.400e-04, 7.834e-02],
        [4.480e-07, 1.920e-07, 1.074e-03, 8.953e-02],
        [5.040e-07, 2.160e-07, 1.209e-03, 1.007e-01],
        [5.600e-07, 2.400e-07, 1.343e-03, 1.119e-01],
        [6.160e-07, 2.640e-07, 1.477e-03, 1.231e-01],
        [6.720e-07, 2.880e-07, 1.612e-03, 1.343e-01],
        [7.280e-07, 3.120e-07, 1.746e-03, 1.455e-01],
        [7.840e-07, 3.360e-07, 1.880e-03, 1.567e-01],
        [8.400e-07, 3.600e-07, 2.014e-03, 1.679e-01],
        [8.960e-07, 3.840e-07, 2.149e-03, 1.791e-01],
        [9.520e-07, 4.080e-07, 2.283e-03, 1.902e-01],
        [1.008e-06, 4.320e-07, 2.417e-03, 2.014e-01],
        [1.064e-06, 4.560e-07, 2.552e-03, 2.126e-01],
        [1.120e-06, 4.800e-07, 2.686e-03, 2.238e-01],
        [1.176e-06, 5.040e-07, 2.820e-03, 2.350e-01],
        [1.232e-06, 5.280e-07, 2.954e-03, 2.462e-01],
        [1.288e-06, 5.520e-07, 3.089e-03, 2.574e-01],
        [1.344e-06, 5.760e-07, 3.223e-03, 2.686e-01],
        [1.400e-06, 6.000e-07, 3.357e-03, 2.798e-01],
        [1.456e-06, 6.240e-07, 3.492e-03, 2.910e-01],
        [1.512e-06, 6.480e-07, 3.626e-03, 3.022e-01],
        [1.568e-06, 6.720e-07, 3.760e-03, 3.133e-01],
        [1.624e-06, 6.961e-07, 3.894e-03, 3.245e-01],
        [1.680e-06, 7.201e-07, 4.029e-03, 3.357e-01],
        [1.736e-06, 7.441e-07, 4.163e-03, 3.469e-01],
        [1.792e-06, 7.681e-07, 4.297e-03, 3.581e-01],
        [1.848e-06, 7.921e-07, 4.432e-03, 3.693e-01],
        [1.904e-06, 8.161e-07, 4.566e-03, 3.805e-01],
        [1.960e-06, 8.401e-07, 4.700e-03, 3.917e-01],
        [2.016e-06, 8.641e-07, 4.835e-03, 4.029e-01],
        [2.072e-06, 8.881e-07, 4.969e-03, 4.141e-01],
        [2.128e-06, 9.121e-07, 5.103e-03, 4.253e-01],
        [2.184e-06, 9.361e-07, 5.237e-03, 4.364e-01],
        [2.240e-06, 9.601e-07, 5.372e-03, 4.476e-01],
        [2.296e-06, 9.841e-07, 5.506e-03, 4.588e-01],
        [2.352e-06, 1.008e-06, 5.640e-03, 4.700e-01],
        [2.408e-06, 1.032e-06, 5.775e-03, 4.812e-01],
        [2.464e-06, 1.056e-06, 5.909e-03, 4.924e-01],
        [2.520e-06, 1.080e-06, 6.043e-03, 5.036e-01],
        [2.576e-06, 1.104e-06, 6.177e-03, 5.148e-01],
        [2.632e-06, 1.128e-06, 6.312e-03, 5.260e-01],
        [2.688e-06, 1.152e-06, 6.446e-03, 5.372e-01],
        [2.744e-06, 1.176e-06, 6.580e-03, 5.484e-01],
        [2.800e-06, 1.200e-06, 6.715e-03, 5.596e-01],
        [2.856e-06, 1.224e-06, 6.849e-03, 5.707e-01],
        [2.912e-06, 1.248e-06, 6.983e-03, 5.819e-01],
        [2.968e-06, 1.272e-06, 7.118e-03, 5.931e-01],
        [3.024e-06, 1.296e-06, 7.252e-03, 6.043e-01],
        [3.080e-06, 1.320e-06, 7.386e-03, 6.155e-01],
        [3.136e-06, 1.344e-06, 7.520e-03, 6.267e-01],
        [3.192e-06, 1.368e-06, 7.655e-03, 6.379e-01],
        [3.248e-06, 1.392e-06, 7.789e-03, 6.491e-01],
        [3.304e-06, 1.416e-06, 7.923e-03, 6.603e-01],
        [3.360e-06, 1.440e-06, 8.058e-03, 6.715e-01],
        [3.416e-06, 1.464e-06, 8.192e-03, 6.827e-01],
        [3.472e-06, 1.488e-06, 8.326e-03, 6.938e-01],
        [3.528e-06, 1.512e-06, 8.460e-03, 7.050e-01],
        [3.584e-06, 1.536e-06, 8.595e-03, 7.162e-01],
        [3.640e-06, 1.560e-06, 8.729e-03, 7.274e-01],
        [3.696e-06, 1.584e-06, 8.863e-03, 7.386e-01],
        [3.752e-06, 1.608e-06, 8.998e-03, 7.498e-01],
        [3.808e-06, 1.632e-06, 9.132e-03, 7.610e-01],
        [3.864e-06, 1.656e-06, 9.266e-03, 7.722e-01],
        [3.920e-06, 1.680e-06, 9.400e-03, 7.834e-01],
        [3.976e-06, 1.704e-06, 9.535e-03, 7.946e-01],
        [4.032e-06, 1.728e-06, 9.669e-03, 8.058e-01],
        [4.088e-06, 1.752e-06, 9.803e-03, 8.169e-01],
        [4.144e-06, 1.776e-06, 9.938e-03, 8.281e-01],
        [4.200e-06, 1.800e-06, 1.007e-02, 8.393e-01],
        [4.256e-06, 1.824e-06, 1.021e-02, 8.505e-01],
        [4.312e-06, 1.848e-06, 1.034e-02, 8.617e-01],
        [4.368e-06, 1.872e-06, 1.047e-02, 8.729e-01],
        [4.424e-06, 1.896e-06, 1.061e-02, 8.841e-01],
        [4.480e-06, 1.920e-06, 1.074e-02, 8.953e-01],
        [4.536e-06, 1.944e-06, 1.088e-02, 9.065e-01],
        [4.592e-06, 1.968e-06, 1.101e-02, 9.177e-01],
        [4.648e-06, 1.992e-06, 1.115e-02, 9.289e-01],
        [4.704e-06, 2.016e-06, 1.128e-02, 9.400e-01],
        [4.760e-06, 2.040e-06, 1.141e-02, 9.512e-01],
        [4.816e-06, 2.064e-06, 1.155e-02, 9.624e-01],
        [4.872e-06, 2.088e-06, 1.168e-02, 9.736e-01],
        [4.928e-06, 2.112e-06, 1.182e-02, 9.848e-01],
        [4.984e-06, 2.136e-06, 1.195e-02, 9.960e-01],
        [5.040e-06, 2.160e-06, 1.209e-02, 1.007e+00],
        [5.096e-06, 2.184e-06, 1.222e-02, 1.018e+00],
        [5.152e-06, 2.208e-06, 1.235e-02, 1.030e+00],
        [5.208e-06, 2.232e-06, 1.249e-02, 1.041e+00],
        [5.264e-06, 2.256e-06, 1.262e-02, 1.052e+00],
        [5.320e-06, 2.280e-06, 1.276e-02, 1.063e+00],
        [5.376e-06, 2.304e-06, 1.289e-02, 1.074e+00],
        [5.432e-06, 2.328e-06, 1.303e-02, 1.086e+00],
        [5.488e-06, 2.352e-06, 1.316e-02, 1.097e+00],
        [5.544e-06, 2.376e-06, 1.329e-02, 1.108e+00],
        [5.600e-06, 2.400e-06, 1.343e-02, 1.119e+00],
        [5.656e-06, 2.424e-06, 1.356e-02, 1.130e+00],
        [5.712e-06, 2.448e-06, 1.370e-02, 1.141e+00],
        [5.768e-06, 2.472e-06, 1.383e-02, 1.153e+00],
        [5.824e-06, 2.496e-06, 1.397e-02, 1.164e+00],
        [5.880e-06, 2.520e-06, 1.410e-02, 1.175e+00],
        [5.936e-06, 2.544e-06, 1.424e-02, 1.186e+00],
        [5.992e-06, 2.568e-06, 1.437e-02, 1.197e+00],
        [6.048e-06, 2.592e-06, 1.450e-02, 1.209e+00],
        [6.104e-06, 2.616e-06, 1.464e-02, 1.220e+00],
        [6.160e-06, 2.640e-06, 1.477e-02, 1.231e+00],
        [6.216e-06, 2.664e-06, 1.491e-02, 1.242e+00],
        [6.272e-06, 2.688e-06, 1.504e-02, 1.253e+00],
        [6.328e-06, 2.712e-06, 1.518e-02, 1.265e+00],
        [6.384e-06, 2.736e-06, 1.531e-02, 1.276e+00],
        [6.440e-06, 2.760e-06, 1.544e-02, 1.287e+00],
        [6.496e-06, 2.784e-06, 1.558e-02, 1.298e+00],
        [6.552e-06, 2.808e-06, 1.571e-02, 1.309e+00],
        [6.608e-06, 2.832e-06, 1.585e-02, 1.321e+00],
        [6.664e-06, 2.856e-06, 1.598e-02, 1.332e+00],
        [6.720e-06, 2.880e-06, 1.612e-02, 1.343e+00],
        [6.776e-06, 2.904e-06, 1.625e-02, 1.354e+00],
        [6.832e-06, 2.928e-06, 1.638e-02, 1.365e+00],
        [6.888e-06, 2.952e-06, 1.652e-02, 1.376e+00],
        [6.944e-06, 2.976e-06, 1.665e-02, 1.388e+00],
        [7.000e-06, 3.000e-06, 1.679e-02, 1.399e+00],
        [7.056e-06, 3.024e-06, 1.692e-02, 1.410e+00],
        [7.112e-06, 3.048e-06, 1.706e-02, 1.421e+00],
        [7.168e-06, 3.072e-06, 1.719e-02, 1.432e+00],
        [7.224e-06, 3.096e-06, 1.732e-02, 1.444e+00],
        [7.280e-06, 3.120e-06, 1.746e-02, 1.455e+00],
        [7.336e-06, 3.144e-06, 1.759e-02, 1.466e+00],
        [7.392e-06, 3.168e-06, 1.773e-02, 1.477e+00],
        [7.448e-06, 3.192e-06, 1.786e-02, 1.488e+00],
        [7.504e-06, 3.216e-06, 1.800e-02, 1.500e+00],
        [7.560e-06, 3.240e-06, 1.813e-02, 1.511e+00],
        [7.616e-06, 3.264e-06, 1.826e-02, 1.522e+00],
        [7.672e-06, 3.288e-06, 1.840e-02, 1.533e+00],
        [7.728e-06, 3.312e-06, 1.853e-02, 1.544e+00],
        [7.784e-06, 3.336e-06, 1.867e-02, 1.556e+00],
        [7.840e-06, 3.360e-06, 1.880e-02, 1.567e+00],
        [7.896e-06, 3.384e-06, 1.894e-02, 1.578e+00],
        [7.952e-06, 3.408e-06, 1.907e-02, 1.589e+00],
        [8.008e-06, 3.432e-06, 1.920e-02, 1.600e+00],
        [8.064e-06, 3.456e-06, 1.934e-02, 1.612e+00],
        [8.120e-06, 3.480e-06, 1.947e-02, 1.623e+00],
        [8.176e-06, 3.504e-06, 1.961e-02, 1.634e+00],
        [8.232e-06, 3.528e-06, 1.974e-02, 1.645e+00],
        [8.288e-06, 3.552e-06, 1.987e-02, 1.656e+00],
        [8.345e-06, 3.575e-06, 2.000e-02, 1.667e+00],
        [8.402e-06, 3.598e-06, 2.013e-02, 1.677e+00],
        [8.459e-06, 3.621e-06, 2.026e-02, 1.688e+00],
        [8.517e-06, 3.644e-06, 2.038e-02, 1.698e+00],
        [8.576e-06, 3.666e-06, 2.049e-02, 1.708e+00],
        [8.636e-06, 3.687e-06, 2.060e-02, 1.717e+00],
        [8.696e-06, 3.708e-06, 2.071e-02, 1.726e+00],
        [8.758e-06, 3.728e-06, 2.081e-02, 1.734e+00],
        [8.821e-06, 3.747e-06, 2.090e-02, 1.741e+00],
        [8.884e-06, 3.766e-06, 2.098e-02, 1.748e+00],
        [8.949e-06, 3.783e-06, 2.106e-02, 1.755e+00],
        [9.016e-06, 3.799e-06, 2.113e-02, 1.761e+00],
        [9.083e-06, 3.814e-06, 2.119e-02, 1.766e+00],
        [9.152e-06, 3.827e-06, 2.124e-02, 1.770e+00],
        [9.223e-06, 3.839e-06, 2.129e-02, 1.774e+00],
        [9.295e-06, 3.850e-06, 2.133e-02, 1.777e+00],
        [9.368e-06, 3.860e-06, 2.136e-02, 1.780e+00],
        [9.442e-06, 3.868e-06, 2.138e-02, 1.782e+00],
        [9.518e-06, 3.875e-06, 2.140e-02, 1.783e+00],
        [9.595e-06, 3.881e-06, 2.141e-02, 1.784e+00],
        [9.673e-06, 3.886e-06, 2.141e-02, 1.784e+00],
        [9.751e-06, 3.889e-06, 2.141e-02, 1.784e+00],
        [9.831e-06, 3.892e-06, 2.141e-02, 1.784e+00],
        [9.911e-06, 3.894e-06, 2.140e-02, 1.783e+00],
        [9.992e-06, 3.896e-06, 2.139e-02, 1.783e+00],
        [1.007e-05, 3.896e-06, 2.138e-02, 1.782e+00],
        [1.015e-05, 3.897e-06, 2.137e-02, 1.781e+00],
        [1.024e-05, 3.897e-06, 2.135e-02, 1.779e+00],
        [1.032e-05, 3.896e-06, 2.134e-02, 1.778e+00],
        [1.040e-05, 3.896e-06, 2.132e-02, 1.777e+00],
        [1.048e-05, 3.895e-06, 2.131e-02, 1.776e+00],
        [1.056e-05, 3.893e-06, 2.129e-02, 1.774e+00],
        [1.065e-05, 3.892e-06, 2.127e-02, 1.773e+00],
        [1.073e-05, 3.891e-06, 2.126e-02, 1.771e+00],
        [1.081e-05, 3.889e-06, 2.124e-02, 1.770e+00],
        [1.089e-05, 3.887e-06, 2.122e-02, 1.769e+00],
        [1.098e-05, 3.885e-06, 2.120e-02, 1.767e+00],
        [1.106e-05, 3.882e-06, 2.119e-02, 1.765e+00],
        [1.114e-05, 3.880e-06, 2.117e-02, 1.764e+00],
        [1.122e-05, 3.877e-06, 2.114e-02, 1.762e+00],
        [1.131e-05, 3.874e-06, 2.112e-02, 1.760e+00],
        [1.139e-05, 3.870e-06, 2.109e-02, 1.758e+00],
        [1.148e-05, 3.866e-06, 2.107e-02, 1.755e+00],
        [1.156e-05, 3.862e-06, 2.103e-02, 1.753e+00],
        [1.165e-05, 3.857e-06, 2.100e-02, 1.750e+00],
        [1.173e-05, 3.852e-06, 2.097e-02, 1.747e+00],
        [1.182e-05, 3.847e-06, 2.093e-02, 1.744e+00],
        [1.190e-05, 3.842e-06, 2.090e-02, 1.741e+00],
        [1.199e-05, 3.836e-06, 2.086e-02, 1.738e+00],
        [1.208e-05, 3.830e-06, 2.082e-02, 1.735e+00],
        [1.216e-05, 3.824e-06, 2.078e-02, 1.732e+00],
        [1.225e-05, 3.818e-06, 2.074e-02, 1.729e+00],
        [1.234e-05, 3.812e-06, 2.070e-02, 1.725e+00],
        [1.242e-05, 3.807e-06, 2.067e-02, 1.722e+00],
        [1.251e-05, 3.801e-06, 2.063e-02, 1.719e+00],
        [1.260e-05, 3.795e-06, 2.059e-02, 1.716e+00],
        [1.268e-05, 3.789e-06, 2.056e-02, 1.713e+00],
        [1.277e-05, 3.784e-06, 2.052e-02, 1.710e+00],
        [1.286e-05, 3.778e-06, 2.048e-02, 1.707e+00],
        [1.294e-05, 3.772e-06, 2.045e-02, 1.704e+00],
        [1.303e-05, 3.767e-06, 2.041e-02, 1.701e+00],
        [1.311e-05, 3.761e-06, 2.038e-02, 1.698e+00],
        [1.320e-05, 3.756e-06, 2.034e-02, 1.695e+00],
        [1.329e-05, 3.750e-06, 2.031e-02, 1.692e+00],
        [1.337e-05, 3.745e-06, 2.028e-02, 1.690e+00],
        [1.346e-05, 3.740e-06, 2.024e-02, 1.687e+00],
        [1.354e-05, 3.735e-06, 2.021e-02, 1.684e+00],
        [1.363e-05, 3.730e-06, 2.018e-02, 1.682e+00],
        [1.371e-05, 3.725e-06, 2.015e-02, 1.679e+00],
        [1.380e-05, 3.720e-06, 2.011e-02, 1.676e+00],
        [1.389e-05, 3.714e-06, 2.008e-02, 1.674e+00],
        [1.397e-05, 3.709e-06, 2.005e-02, 1.671e+00],
        [1.406e-05, 3.704e-06, 2.002e-02, 1.668e+00],
        [1.414e-05, 3.699e-06, 1.999e-02, 1.666e+00],
        [1.423e-05, 3.694e-06, 1.995e-02, 1.663e+00],
        [1.431e-05, 3.689e-06, 1.992e-02, 1.660e+00]])]

Vykreslení ld křivek se zvýrazněním maxim

In [3]:
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]');
In [4]:
data = data_all[0]
displ = (data[:, 0] - data[:, 1]) * 1e6
force = data[:, 3]

fidx = np.argmax(force)
displ /= displ[fidx]
force /= force[fidx]
In [5]:
plt.plot(displ, force)
Out[5]:
[]
In [6]:
import pandas as pd
df = pd.read_csv('http://ws-cheetah.stm.fce.vutbr.cz/vyuka/CD004/priklady/tlakova_zkouska/data.txt', sep=';', header=None)
fcube = df.values.flatten()
In [7]:
force.shape
Out[7]:
(224,)
In [8]:
for i, fi in enumerate(fcube):
    f = force * fi
    d = displ * (np.random.randn() * .1 + 2)
    plt.plot(d,f)
    np.savetxt('krivka_%02d.txt' % i, np.vstack((d, f)).T, header='nerealne krivky pevnosti v tlaku\nstrain;stress', delimiter=';')

Analýza maxim a jim odpovídajících posunů

In [9]:
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 [24]:
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, normed=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, normed=True)

ax2.set_xlabel('deformace [mm]')

fig.tight_layout()

parametry normálního rozdělení (mu, var) = (array(2.729125), array(0.144352515625))
parametry dvouparametrického weibullova rozdělení (mu, var) = (array(2.728093261007631), array(0.15355103015394583))
In [26]:
print u'{:40} = {:+.6f}'.format(u'střední hodnota maximální síly', np.mean(f_max))
print u'{:40} = {:+.6f}'.format(u'směrodatná odchylka maximální síly', np.std(f_max))
print u'{:40} = {:+.6f}'.format(u'šikmost maximální síly', skew(f_max))
print u'{:40} = {:+.6f}'.format(u'špičatost maximální síly', kurtosis(f_max))

střední hodnota maximální síly           = +2.729125
směrodatná odchylka maximální síly       = +0.379938
šikmost maximální síly                   = -0.266260
špičatost maximální síly                 = -0.408317
In [ ]:
In [ ]: