Figure 10 & 11

Imports

import dolfin
import numpy
import sympy

from numpy import linspace
from sympy import lambdify
from shapely.geometry import LineString

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_inflation_air_filled = numpy.load('smith_PV_inflation_air_filled.npy')
p_smith_PV_inflation_air_filled = smith_PV_inflation_air_filled[:, 0]
v_smith_PV_inflation_air_filled = smith_PV_inflation_air_filled[:, 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]

smith_PV_deflation_air_filled = numpy.load('smith_PV_deflation_air_filled.npy')
p_smith_PV_deflation_air_filled = smith_PV_deflation_air_filled[:, 0]
v_smith_PV_deflation_air_filled = smith_PV_deflation_air_filled[:, 1]

Defining geometry

seeds_filename    = "Fig10.dat"
mesh_filebasename = "Fig10-mesh"
qois_filename     = "Fig10-qois.dat"
res_basename      = "Fig10"

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 & Parameters

params = [0.08855929243285596, 0.011039510924095856, 0.6281487879627474, 3.409513378002055]
mat_params = {"model":"exponentialneoHookean", "parameters":{"beta1":params[0], "beta2":params[1], "beta3":params[2], "beta4":100*params[0], "alpha":params[3]}}  
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

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    mat_params=mat_params,
    load_params=load_params,
    step_params={"Deltat":1., "dt_ini":0.1, "dt_min":0.005, "dt_max":0.1},
    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]
S_lst = qois_vals[:, qois_name_list.index("S_area") - 1]
S_hat_gamma_0_lst = [S/S_lst[0] for S in S_lst]

vf_gamma_0_lst = vf_lst
pf_gamma_0_lst = pf_lst


def loading_model(gamma, p_ini):
    gamma = gamma * 1e-3
    load_params = {}
    load_params["pf_lst"] = [p_ini, 3]
    load_params["sigma_bar_00_lst"] = [0.0, 0.0]
    load_params["sigma_bar_11_lst"] = [0.0, 0.0]
    load_params["sigma_bar_01_lst"] = [0.0, 0.0]
    load_params["sigma_bar_10_lst"] = [0.0, 0.0]
    load_params["gamma_lst"] = [gamma, gamma]

    step_params = {}
    step_params["n_steps"] = 2
    step_params["Deltat"] = 1.
    step_params["dt_ini"] = 0.01
    step_params["dt_min"] = 0.0005
    step_params["dt_max"] = 0.005

    dmech.run_HollowBox_MicroPoroHyperelasticity(
        dim=2,
        mesh=mesh,
        mat_params=mat_params,
        load_params=load_params,
        step_params=step_params,
        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]
    vf_lst = [vf_/vf_asym *100 for vf_ in vf_lst]
    S_lst = qois_vals[:, qois_name_list.index("S_area") - 1]
    S_hat_lst = [S/S_lst[0] for S in S_lst]

    index = numpy.argmin(numpy.abs(pf_lst - p_ini*10.20), axis=None)

    return pf_lst[index:], vf_lst[index:], S_hat_lst[index:]

model_response_lst = []
model_response_lst.append([0 , loading_model(0, 0.)])
model_response_lst.append([3 , loading_model(3, 0.1)])
model_response_lst.append([5 , loading_model(5, 0.1)])
model_response_lst.append([7 , loading_model(7, 0.3)])
model_response_lst.append([10, loading_model(10, 0.3)])
model_response_lst.append([28, loading_model(28, 1.6)])
model_response_lst.append([30, loading_model(30, 1.8)])

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)

plt.plot(p_smith_PV_deflation_gamma_0, v_smith_PV_deflation_gamma_0, '#D94801', linestyle='dotted', label='[Smith, 1986], deflation, $\gamma = 0~dyn/cm$')
plt.plot(p_smith_PV_inflation_gamma_0, v_smith_PV_inflation_gamma_0, '#D94801', label='[Smith, 1986], inflation, $\gamma = 0~dyn/cm$')

plt.plot(model_response_lst[0][1][0], model_response_lst[0][1][1], '#084594', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[0][0])+ '$~dyn/cm$')
plt.plot(model_response_lst[1][1][0], model_response_lst[1][1][1], '#2171B5', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[1][0])+ '$~dyn/cm$')
plt.plot(model_response_lst[2][1][0], model_response_lst[2][1][1], '#4292C6', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[2][0])+ '$~dyn/cm$')
plt.plot(model_response_lst[3][1][0], model_response_lst[3][1][1], '#6BAED6', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[3][0])+ '$~dyn/cm$')
plt.plot(model_response_lst[4][1][0], model_response_lst[4][1][1], '#9ECAE1', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[4][0])+ '$~dyn/cm$')
plt.plot(model_response_lst[5][1][0], model_response_lst[5][1][1], '#C6DBEF', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[5][0])+ '$~dyn/cm$')
plt.plot(model_response_lst[6][1][0], model_response_lst[6][1][1], '#DEEBF7', label='Model, $D_{alv} = 54~\mu m$, $\gamma =$' +str(model_response_lst[6][0])+ '$~dyn/cm$')

plt.plot(p_smith_PV_deflation_air_filled, v_smith_PV_deflation_air_filled, '#D94801', label='[Smith, 1986], air-filled, deflation', linestyle='dashdot')
plt.plot(p_smith_PV_inflation_air_filled, v_smith_PV_inflation_air_filled, '#D94801', label='[Smith, 1986], air-filled, inflation', linestyle='dashed')

plt.xlim(0, 25)
plt.ylim(0, 115)
# plt.legend(loc = 'lower right', fontsize=12, bbox_to_anchor=(1.8, 0.), shadow=True)
plt.legend(loc = 'lower right', fontsize=12, bbox_to_anchor=(0.84, -1.04), shadow=True)
plt.savefig('smith_intersect.pdf',bbox_inches='tight')
plt.savefig("smith_intersect_new.tif", format="tiff", dpi=300, bbox_inches="tight")

Intersection points

def intersect(lst, reference):
    p_lst = lst[0]
    v_lst = lst[1]
    S_lst = lst[2]
    pf_reference = reference[:, 0]
    vf_reference = reference[:, 1]
    first_line = LineString(numpy.column_stack((p_lst, v_lst)))
    second_line = LineString(numpy.column_stack((pf_reference, vf_reference)))
    intersection = first_line.intersection(second_line)
    gamma_S = intersection.x
    v = intersection.y

    first_line = LineString(numpy.column_stack((p_lst, S_lst)))
    second_line = LineString(numpy.column_stack((len(p_lst)*[gamma_S], S_lst)))
    intersection = first_line.intersection(second_line)
    return intersection.y, v

for i in range (1, len(model_response_lst)):
    model_response_lst[i].append(intersect(model_response_lst[i][1], smith_PV_deflation_air_filled))

gamma_lst_def   = [sublist[0]/30 for sublist in model_response_lst[1:]]
S_intersect_def = [sublist[2][0] for sublist in model_response_lst[1:]]
gamma_lst_def   = [0] + gamma_lst_def
S_intersect_def = [1] + S_intersect_def

for i in range (2, len(model_response_lst)):
    model_response_lst[i].append(intersect(model_response_lst[i][1], smith_PV_inflation_air_filled))

gamma_lst_inf   = [sublist[0]/30 for sublist in model_response_lst[2:]]
S_intersect_inf = [sublist[3][0] for sublist in model_response_lst[2:]]
gamma_lst_inf   = [0] + gamma_lst_inf
S_intersect_inf = [1] + S_intersect_inf

S_hat = sympy.symbols("S_hat")

gamma_S_iden_def =  1.030747711797792/(1 + (S_hat/2.1708848554526874)**(-14.828598856766776))
gamma_S_iden_inf = 1.0130287663205635/(1 + (S_hat/1.6742366271475184)**(-10.288589574038403))

lam_def = lambdify(S_hat, gamma_S_iden_def, modules=['numpy'])
lam_inf = lambdify(S_hat, gamma_S_iden_inf, modules=['numpy'])

S_hat_vals = linspace(1, 3, 100)
gamma_S_iden_vals_deflation = lam_def(S_hat_vals)
gamma_S_iden_vals_inflation = lam_inf(S_hat_vals)

plt.rc('xtick' , labelsize=14)
plt.rc('ytick' , labelsize=14)
plt.rc('legend', fontsize=12)
plt.xlabel(r'$S^*$'     , fontsize=16)
plt.ylabel(r'$\gamma^*$', fontsize=16)

plt.plot(S_intersect_inf, gamma_lst_inf, 'o', color='#084594')
plt.plot(S_intersect_def, gamma_lst_def, 'o', color='#238B45')
plt.plot(S_hat_vals, gamma_S_iden_vals_inflation, '#084594', linestyle = 'dashed', label='Inflation')
plt.plot(S_hat_vals, gamma_S_iden_vals_deflation, '#238B45', linestyle = 'dashed', label='Deflation')
plt.legend(loc = 'lower right', fontsize=12, shadow=True)
plt.xlim(1, 2.8)
plt.savefig('inf_def_fit.pdf',bbox_inches='tight')
plt.savefig("inf_def_fit.tif", format="tiff", dpi=300, bbox_inches="tight")