Poroplastic Data Fitting



In [1]:

    
import numpy as np
from numpy import *
from bokeh import *
from bokeh.plotting import *
output_notebook()
from matmodlab2 import *
from pandas import read_excel
from scipy.optimize import leastsq
diff = lambda x: np.ediff1d(x, to_begin=0.)
trace = lambda x, s='SIG': x[s+'11'] + x[s+'22'] + x[s+'33']
RTJ2 = lambda x: sqrt(((x['SIG11']-x['SIG22'])**2 + 
                       (x['SIG22']-x['SIG33'])**2 + 
                       (x['SIG33']-x['SIG22'])**2)/6.)









    





    
        
        Loading BokehJS ...
    






    














    



Setting up the Matmodlab notebook environment

Summary

In the cells to follow, the following material parameters were found

$$\begin{align} B_0 &= 14617807286.8\\ B_1 &= 40384983097.2\\ B_2 &= 385649437.858\\ P_0 & = −164761936.257 \\ P_1 & = 3.20119273834e−10\\ P_2 & = 7.39166987894e−18\\ P_3 & = 0.0983914345654\\ G_1 & = 9647335534.93 \\ G_2 & = 2.3838775292e−09 \\ G_3 & = −7.40942609805e−07\\ \end{align}$$

Read in the Data

Read in the hydrostatic data and compute derived values.



In [2]:

    
df = read_excel('porodata.xlsx', sheetname='hydrostatic')
df['EV'] = trace(df, 'STRAIN')
df['I1'] = trace(df, 'SIG')
df['dEV'] = diff(df['EV'])
df['dI1'] = diff(df['I1'])

Hydrostatic Response

Elastic Unloading Curve

Plot the pressure vs. volume strain curve and determine the section in which elastic unloading occurs



In [3]:

    
plot = figure(x_axis_label='Volume Strain', y_axis_label='Pressure')
plot.circle(-df['EV'], -df['I1']/3.)
plot.text(-df['EV'], -df['I1']/3.,
    text=range(len(df)),text_color="#333333",
    text_align="left", text_font_size="5pt")
show(plot)

It appears that the unloading occurs at data point 101 and continues until the end of the data. This curve will be used to fit the bulk modulus parameters. Below, scipy is used to optimize the parameters to the curve.



In [4]:

    
kfun = lambda B0, B1, B2, I1: B0 + B1 * exp(-B2 / abs(I1))
def kmm_bulk(x, fac, I1, K):
    B0, B1, B2 = x * fac
    return K - kfun(B0, B1, B2, I1)



In [5]:

    
imax = 101
df1 = df.iloc[imax:].copy()
K = np.array(df1['dI1'] / 3. / df1['dEV'])
b0 = np.array((K[-1], K[0] - K[-1], 1e9))
fac = 1e9
B, icov = leastsq(kmm_bulk, b0/fac, args=(fac, df1['I1'], K))
B0, B1, B2 = B * fac
B0, B1, B2









    Out[5]:





(14617807276.85973, 40384983106.145226, 385649437.68283093)



In [6]:

    
plot = figure(x_axis_label='Bulk Modulus', y_axis_label='Pressure')
plot.circle(-df1['I1']/3., K)
plot.line(-df['I1']/3., kfun(B0, B1, B2, df['I1']), color='red')
show(plot)

Poro response

With the bulk response determined, find the porosity parameters



In [7]:

    
df['EP'] = df['I1'] / 3. / kfun(B0, B1, B2, df['I1']) - df['EV']
p3 = max(df['EP'])
df['PORO'] = p3 - df['EP']
plot = figure(x_axis_label='Plastic Strain', y_axis_label='Pressure')
plot.circle(df['EP'], -df['I1']/3.)
show(plot)



In [8]:

    
plot = figure(x_axis_label='Pressure', y_axis_label='PORO')
df2 = df.iloc[:imax].copy()
plot.circle(-df2['I1']/3., df2['PORO'])
show(plot)



In [9]:

    
def pfun(P0, P1, P2, P3, I1):
    xi = -I1 / 3. + P0
    return P3 * exp(-(P1 + P2 * xi) * xi)
    
def kmm_poro(x, fac, I1, P):
    p0, p1, p2, p3 = asarray(x) * fac
    return P - pfun(p0, p1, p2, p3, I1)



In [10]:

    
p0 = (1, 1, 1, p3)
fac = np.array([1e8, 1e-10, 1e-18, 1])
p, icov = leastsq(kmm_poro, p0, args=(fac, df2['I1'], df2['PORO']))
P0, P1, P2, P3 = p * fac
P0, P1, P2, P3









    Out[10]:





(-28446120.517447289,
 -1.6950842119972055e-09,
 7.3916733442463017e-18,
 0.08958961426627815)



In [11]:

    
plot = figure(x_axis_label='Pressure', y_axis_label='PORO')
plot.circle(-df2['I1']/3., df2['PORO'], legend='Data')
plot.line(-df2['I1']/3., pfun(P0, P1, P2, P3, df2['I1']), color='red', legend='Fit')
show(plot)

Shear Response



In [12]:

    
keys = (2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 22.5, 30.0)
colors = ('red', 'blue', 'orange', 'purple', 
          'green', 'black', 'magenta', 'teal', 'cyan')
df2 = {}
p = figure(x_axis_label='I1', y_axis_label='Sqrt[J2]')
p1 = figure(x_axis_label='Axial Strain', y_axis_label='Axial Stress')
for (i, key) in enumerate(keys):
    key = 'txc p={0:.01f}MPa'.format(key)
    x = read_excel('porodata.xlsx', sheetname=key)
    x['I1'] = trace(x, 'SIG')
    x['RTJ2'] = RTJ2(x)
    df2[key] = x
    p.circle(-df2[key]['I1'], df2[key]['RTJ2'], legend=key[4:], color=colors[i])
    
    # determine where hydrostatic preload ends
    j = nonzero(x['SIG11'] - x['SIG22'])[0]
    E0, S0 = df2[key]['STRAIN11'][j[0]], df2[key]['SIG11'][j[0]]
    p1.circle(-df2[key]['STRAIN11'][j]+E0, -df2[key]['SIG11'][j]+S0,
              legend=key[4:], color=colors[i])

p.legend.orientation = 'horizontal'
show(p1)
show(p)

The axial stress versus axial strain plot shows that the response is linear, meaning that the elastic modulus is constant.



In [13]:

    
key = 'txc p=2.5MPa'
j = nonzero(df2[key]['SIG11'] - df2[key]['SIG22'])[0]
df3 = df2[key].iloc[j].copy()
E0, S0 = df3['STRAIN11'].iloc[0], df3['SIG11'].iloc[0]
EF, SF = df3['STRAIN11'].iloc[-1], df3['SIG11'].iloc[-1]
E = (SF - S0) / (EF - E0)
print('{0:E}'.format(E))









    



2.372314E+10

The shear modulus can now be determined



In [14]:

    
G = lambda I1: 3 * kfun(B0, B1, B2, I1) * E / (9 * kfun(B0, B1, B2, I1) - E)
gfun = lambda g0, g1, g2, rtj2: g0 * (1 - g1 * exp(-g2 * rtj2)) / (1 - g1)
def kmm_shear(x, fac, rtj2, G):
    g0, g1, g2 = asarray(x) * fac
    return G - gfun(g0, g1, g2, rtj2)



In [15]:

    
g = asarray(G(df3['I1']))
g0 = (g[0], .0001, 0)
fac = 1.
g, icov = leastsq(kmm_shear, g0, args=(fac, RTJ2(df3), g))
G0, G1, G2 = g * fac
G0, G1, G2









    Out[15]:





(9647335536.3045807, 2.3838678121801188e-09, -7.4094283435851223e-07)



In [16]:

    
p2 = figure(x_axis_label='Sqrt[J2]', y_axis_label='Shear Modulus')
p2.circle(RTJ2(df3), G(df3['I1']))
p2.line(RTJ2(df3), gfun(G0, G1, G2, RTJ2(df3)), color='red')
show(p2)



In [ ]: