Figure 7

Imports

import dolfin
import numpy
import scipy.optimize

from scipy.interpolate import interp1d

import dolfin_mech                    as dmech
import matplotlib.pyplot              as plt
import micro_poro_structure_generator as gen

Importing experimental data

smith_PV_inflation_gamma_0   = numpy.load('smith_PV_inflation_gamma_0.npy')
p_smith_PV_inflation_gamma_0 = smith_PV_inflation_gamma_0[:, 0]
v_smith_PV_inflation_gamma_0 = smith_PV_inflation_gamma_0[:, 1]

smith_PV_deflation_gamma_0   = numpy.load('smith_PV_deflation_gamma_0.npy')
p_smith_PV_deflation_gamma_0 = smith_PV_deflation_gamma_0[:, 0]
v_smith_PV_deflation_gamma_0 = smith_PV_deflation_gamma_0[:, 1]

Defining geometry

seeds_filename    = "Fig7-seeds.dat"
mesh_filebasename = "Fig7-mesh"
qois_filename     = "Fig7-qois.dat"
res_basename      = "Fig7"

domain_y  = 0.1 * 0.8
domain_x  = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.012 * 0.8

gen.generate_seeds_semi_regular(
    DoI = 0.,
    row = 1,
    domain_y = domain_y,
    seeds_filename = seeds_filename)
gen.generate_mesh_2D_rectangle_w_voronoi_inclusions(
    mesh_filename = mesh_filebasename,
    seeds_filename = seeds_filename,
    h = thickness,
    lcar = thickness/5,
    domain_x = domain_x,
    domain_y = domain_y,
    shift_y = 0.,
    remove_seeds = True)

mesh = dolfin.Mesh()
dolfin.XDMFFile(mesh_filebasename+".xdmf").read(mesh)
dV = dolfin.Measure("dx",domain=mesh)

coord = mesh.coordinates()
xmax = max(coord[:,0]); xmin = min(coord[:,0])
ymax = max(coord[:,1]); ymin = min(coord[:,1])

V   = (xmax - xmin)*(ymax - ymin)
VS0 = dolfin.assemble(dolfin.Constant(1) * dV)
Vf0 = V - VS0

Loading

load_params = {}
load_params["pf"] = 4
load_params["sigma_bar_00"] = 0.0
load_params["sigma_bar_11"] = 0.0
load_params["sigma_bar_01"] = 0.0
load_params["sigma_bar_10"] = 0.0

Model response

def pressure_computer(params):
    mat_params = {"model":"exponentialneoHookean", "parameters":{"beta1":params[0], "beta2":params[1], "beta3":params[2], "beta4":100*params[0], "alpha":params[3]}}    
    dmech.run_HollowBox_MicroPoroHyperelasticity(
        dim=2,
        mesh=mesh,
        mat_params=mat_params,
        load_params=load_params,
        step_params={"Deltat":1., "dt_ini":0.01, "dt_min":0.005, "dt_max":0.01},
        res_basename=res_basename,
        write_qois_limited_precision=False,
        verbose=1
    )

    qois_vals = numpy.loadtxt(qois_filename)
    qois_name_list = open(qois_filename).readline().split()
    pf_lst = qois_vals[:, qois_name_list.index("p_f") - 1]*10.20
    vf_lst = qois_vals[:, qois_name_list.index("vf") - 1]

    for i in range(1, len(vf_lst)):
        slope = (vf_lst[i] - vf_lst[i - 1])/(pf_lst[i] - pf_lst[i - 1])
        if slope < 0.1 * Vf0:
            break

    vf_asym = vf_lst[i]
    vf_lst = [vf_/vf_asym *100 for vf_ in vf_lst]

    return pf_lst, vf_lst

Plot

params = [0.08855929243285596, 0.011039510924095856, 0.3, 3.409513409247871]
pf_lst_1, vf_lst_1 = pressure_computer(params)
params = [0.08855929243285596, 0.011039510924095856, 0.6281487879627474, 3.409513409247871]
pf_lst_2, vf_lst_2 = pressure_computer(params)
params = [0.08855929243285596, 0.011039510924095856, 0.8, 3.409513409247871]
pf_lst_3, vf_lst_3 = pressure_computer(params)
params = [0.08855929243285596, 0.011039510924095856, 1, 3.409513409247871]
pf_lst_4, vf_lst_4 = pressure_computer(params)

plt.plot()
plt.rc('xtick' , labelsize=14)
plt.rc('ytick' , labelsize=14)
plt.rc('legend', fontsize=12)
plt.xlabel(r'$p_f~(cm H_2O)$', fontsize=16)
plt.ylabel(r'$Volume~(\% TLC)$', fontsize=16)

colors = ["#DEEBF7", "#9ECAE1", "#4292C6", "#084594"]  # light to dark blue
labels = [
    r'Model, $\beta_3 = 300~Pa$',
    r'Model, $\beta_3 = 626~Pa$',
    r'Model, $\beta_3 = 800~Pa$',
    r'Model, $\beta_3 = 1000~Pa$'
]
plt.plot(pf_lst_1[0:80], vf_lst_1[0:80], color=colors[0], label=labels[0])
plt.plot(pf_lst_2[0:80], vf_lst_2[0:80], color=colors[1], label=labels[1])
plt.plot(pf_lst_3[0:80], vf_lst_3[0:80], color=colors[2], label=labels[2])
plt.plot(pf_lst_4[0:80], vf_lst_4[0:80], color=colors[3], label=labels[3])

plt.xlim(0, 25)
plt.ylim(0, 112)
plt.legend(loc = 'lower right', fontsize=12, shadow=True)
plt.savefig('PV_beta3_variation.pdf',bbox_inches='tight')
plt.savefig("PV_beta3_variation.tif", format="tiff", dpi=300, bbox_inches="tight")
plt.show()