Figure 12¶
Imports¶
In [ ]:
import dolfin
import numpy
import dolfin_mech as dmech
import matplotlib.pyplot as plt
import micro_poro_structure_generator as gen
Importing experimental data¶
In [ ]:
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¶
In [ ]:
seeds_filename = "Fig12.dat"
mesh_filebasename = "Fig12-mesh"
qois_filename = "Fig12-qois.dat"
res_basename = "Fig12"
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¶
In [ ]:
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 respone¶
$\gamma = 0$¶
In [ ]:
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
$\gamma = \gamma(S)$¶
In [ ]:
load_params = {}
load_params["pf_lst"] = [0.1, 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"] = [0.03, 0.03]
load_params["tension_params"] = {"surface_dependancy":1, "d1":1.0130287663205635, "d2":1.6742366271475184, "d3":-10.288589574038403}
step_params = {}
step_params["n_steps"] = 2
step_params["Deltat"] = 1.
step_params["dt_ini"] = 0.01
step_params["dt_min"] = 0.001
step_params["dt_max"] = 0.005
phi = 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]
vf_inf_lst = vf_lst
pf_inf_lst = pf_lst
load_params["tension_params"] = {"surface_dependancy":1, "d1":1.030747711797792, "d2":2.1708848554526874, "d3":-14.828598856766776}
phi = 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]
vf_def_lst = vf_lst
pf_def_lst = pf_lst
In [ ]:
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', label='[Smith, 1986], deflation, $\gamma = 0~dyn/cm$', linestyle='dotted')
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(pf_gamma_0_lst, vf_gamma_0_lst, '#084594', label='Model, $\gamma = 0~dyn/cm$')
plt.plot(pf_def_lst, vf_def_lst , '#084594', label='Model, deflation, $\gamma = \gamma(S)$', linestyle='dashdot')
plt.plot(pf_inf_lst, vf_inf_lst , '#084594', label='Model, inflation, $\gamma = \gamma(S)$', linestyle='dashed')
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.savefig('surface_dependent_surface_tension.pdf',bbox_inches='tight')
plt.savefig("surface_dependent_surface_tension.tif", format="tiff", dpi=300, bbox_inches="tight")