Import the modules
In [ ]:
import scipy as sp
import numpy as np
import openpnm as op
import matplotlib.pyplot as plt
import openpnm.models.geometry as gm
import openpnm.topotools as tt
%matplotlib inline
np.random.seed(10)
Set the workspace loglevel to not print anything
In [ ]:
ws = op.Workspace()
ws.settings["loglevel"] = 50
In [ ]:
%run shared_funcs.ipynb
We can run multiple times as the network sizes are randomly generated between a given range we can obtain an average
In [ ]:
x_values = []
y_values = []
for ensemble in range(10):
x_ensemble, y_ensemble = simulation(n=8)
x_values.append(x_ensemble)
y_values.append(y_ensemble)
In [ ]:
x_values = np.asarray(x_values).flatten()
y_values = np.asarray(y_values).flatten()
plt.figure()
from matplotlib.font_manager import FontProperties
fontP = FontProperties()
fontP.set_size('small')
wu_average_x_values = [0.004, 0.021, 0.052, 0.081, 0.129, 0.162, 0.186, 0.219, 0.261,
0.286, 0.324, 0.363, 0.42, 0.478, 0.531, 0.586, 0.64, 0.698, 0.747, 0.802]
wu_average_y_values = [0.118, 0.113, 0.105, 0.096, 0.085, 0.078, 0.07, 0.062, 0.054, 0.049, 0.04,
0.033, 0.027, 0.02, 0.012, 0.006, 0.003, 0.002, 0.002, 0.002]
p1, = plt.plot(x_values, y_values, 'ko')
p2, = plt.plot(wu_average_x_values, wu_average_y_values, 'ro')
plt.title('normalized diffusivity versus saturation')
plt.xlabel('saturation')
plt.ylabel(r'$\frac{D_e}{D_b}$')
#plt.ylim([0, .15])
plt.xlim([0, 1])
plt.legend([p1, p2],
[r'$\frac{D_e}{D_b} = f(\epsilon, \phi)g(s, \phi)$' + '\n' + r'$X = 1.8$' +
'\n' + r'$Z_t = 2.0$' + '\n' + r'$Z_i = 4.0$' + '\n' + r'$\beta = 1.0$' + '\n' + r'$n = 14$', "Wu's results"])
plt.show()
And finally extract the g(S) function for relative diffusivity.
In [ ]:
plt.figure()
normalize_factor = max(y_values)
g_values = y_values / normalize_factor
wu_saturation = [0.004, 0.066, 0.0930, .119, 0.14, 0.175, 0.209, 0.24, 0.282, 0.32, 0.371, 0.413,
0.464, 0.517, 0.605, 0.672, 0.761, 0.831, 0.898, 0.948, 0.996]
wu_g_values = [0.986, 0.838, 0.758, 0.701, 0.651, 0.576, 0.516, 0.456, 0.39, 0.335, 0.268, 0.221,
0.171, 0.111, 0.067, 0.04, 0.019, 0.007, 0.003, 0.003, 0.003]
p1, = plt.plot(x_values, g_values, 'ko')
p2, = plt.plot(wu_saturation, wu_g_values, 'ro')
plt.title('g(s) versus saturation')
plt.xlabel('saturation')
plt.ylabel('g(s)')
plt.legend([p1, p2],
["our values", "Wu's values (fitted curve)"], loc='center left', bbox_to_anchor=(1, 0.5), prop = fontP)
plt.show()