Hyperelastic Model Fitting



In [1]:

    
%load_ext autoreload
%autoreload 2
from numpy import *
import numpy as np
from bokeh.plotting import *
from pandas import read_excel
from matmodlab2.fitting.hyperopt import *
output_notebook()









    



Setting up the Matmodlab notebook environment






    





    
        
        Loading BokehJS ...

Experimental Data



In [2]:

    
# uniaxial data
udf = read_excel('Treloar_hyperelastic_data.xlsx', sheetname='Uniaxial')
ud = udf.as_matrix(columns=('Engineering Strain', 'Engineering Stress (MPa)'))

# Biaxial data
bdf = read_excel('Treloar_hyperelastic_data.xlsx', sheetname='Biaxial')
bd = bdf.as_matrix(columns=('Engineering Strain', 'Engineering Stress (MPa)'))

# Pure shear data
sdf = read_excel('Treloar_hyperelastic_data.xlsx', sheetname='Pure Shear')
sd = sdf.as_matrix(columns=('Engineering Strain', 'Engineering Stress (MPa)'))

Uniaxial Data

Find the optimal fit to the uniaxial stress data with a hyperelastic polynomial model. The symbolic constant UNIAXIAL_DATA instructs hyperopt to interpret the input data as coming from a uniaxial stress experiment.



In [3]:

    
uf = hyperopt(UNIAXIAL_DATA, ud[:,0], ud[:,1])
print(uf.summary())









    



            Data type: Uniaxial
Number of data points: 23
     Polynomial order: 4
        I2 dependence: True
           Parameters: C10=9783822746.674, C01=-3260926703.250, C20=877819255.795, C11=-189535803.642, C02=10417237.759, C30=-32317870.391, C21=128676376.268, C12=31082108.870, C03=-127257924.315, C40=8765132.778, C31=7155466.130, C22=-37641058.416, C13=38503721.758, C04=-43648109.681
                Error: 0.02076748143710473

At this point, the optimal parameters have been determined and are accessible with the popt attribute:



In [4]:

    
uf.popt









    Out[4]:





array([  9.78382275e+09,  -3.26092670e+09,   8.77819256e+08,
        -1.89535804e+08,   1.04172378e+07,  -3.23178704e+07,
         1.28676376e+08,   3.10821089e+07,  -1.27257924e+08,
         8.76513278e+06,   7.15546613e+06,  -3.76410584e+07,
         3.85037218e+07,  -4.36481097e+07])

The optimal parameters are also available as a dictionary via the todict method:



In [5]:

    
uf.todict()









    Out[5]:





{'C01': -3260926703.2502761,
 'C02': 10417237.759391822,
 'C03': -127257924.31490304,
 'C04': -43648109.680867009,
 'C10': 9783822746.674469,
 'C11': -189535803.64177945,
 'C12': 31082108.870263141,
 'C13': 38503721.757908188,
 'C20': 877819255.795259,
 'C21': 128676376.26756807,
 'C22': -37641058.415643252,
 'C30': -32317870.390860472,
 'C31': 7155466.1302872244,
 'C40': 8765132.7779212315}

The error in the fit:



In [6]:

    
uf.error









    Out[6]:





0.020767481437104732

Plots are generated with the bp_plot method



In [7]:

    
show(uf.bp_plot())

Biaxial Data

Biaxial data is fit in a similar manner:



In [8]:

    
bf = hyperopt(BIAXIAL_DATA, bd[:,0], bd[:,1])
print(bf.summary())
show(bf.bp_plot())









    



            Data type: Biaxial
Number of data points: 16
     Polynomial order: 4
        I2 dependence: True
           Parameters: C10=48126616.023, C01=-14323537.793, C20=597202131.473, C11=-420483882.811, C02=108970996.079, C30=416916020.562, C21=-85123697.713, C12=-35925284.204, C03=-3410748.604, C40=105722602.577, C31=9671542.604, C22=2942836.173, C13=-267685.707, C04=16785.785
                Error: 0.00990862004186524

Shear Data

Lastly, the shear data is fit



In [9]:

    
sf = hyperopt(SHEAR_DATA, sd[:,0], sd[:,1])
print(sf.summary())
show(sf.bp_plot())









    



            Data type: Shear
Number of data points: 13
     Polynomial order: 3
        I2 dependence: True
           Parameters: C10=-480403082.071, C01=161106124.736, C20=-278863992.929, C11=38177864.867, C02=56699970.680, C30=-29401178.157, C21=-4432.913, C12=2952118.412, C03=-752519.876
                Error: 0.018242925241750776

Comparison of Fits

Examine the error in the shear fit using parameters from the uniaxial fit



In [10]:

    
y1 = sf.eval(overlay=uf)
y2 = sf.eval()
err = sqrt(mean((y1-y2)**2)) / average(abs(y2))
print(err)
show(sf.bp_plot(overlay=[bf, uf]))
show(uf.bp_plot(overlay=[sf]))
show(bf.bp_plot(overlay=[uf, sf]))









    



0.0859995309373

hyperopt2

hyperopt2 attempts to find the model that fits all given data the best.



In [11]:

    
f = hyperopt2(SHEAR_DATA, sd[:,0], sd[:,1], 
              UNIAXIAL_DATA, ud[:,0], ud[:,1],
              BIAXIAL_DATA, bd[:,0], bd[:,1])



In [12]:

    
print(f.summary())









    



            Data type: Uniaxial
Number of data points: 23
     Polynomial order: 4
        I2 dependence: False
           Parameters: C10=430761.746, C20=46252.207, C30=-22262.415, C40=3276.462
                Error: 0.044158206214665915



In [13]:

    
f.error2









    Out[13]:





0.017267590200161907



In [14]:

    
p = f.bp_plot(strain=linspace(0,6.5), points=False)
p.circle(sd[:,0], sd[:,1], color='black', legend='Shear data')
p.circle(bd[:,0], bd[:,1], color='red', legend='Biaxial data')
p.circle(ud[:,0], ud[:,1], color='green', legend='Uniaxial data')
show(p)



In [ ]: