Figure 9¶
Imports¶
In [ ]:
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¶
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_gamma_7 = numpy.load('smith_PV_inflation_gamma_7.npy')
p_smith_PV_inflation_gamma_7 = smith_PV_inflation_gamma_7[:, 0]
v_smith_PV_inflation_gamma_7 = smith_PV_inflation_gamma_7[:, 1]
smith_PV_inflation_gamma_16 = numpy.load('smith_PV_inflation_gamma_16.npy')
p_smith_PV_inflation_gamma_16 = smith_PV_inflation_gamma_16[:, 0]
v_smith_PV_inflation_gamma_16 = smith_PV_inflation_gamma_16[:, 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_gamma_7 = numpy.load('smith_PV_deflation_gamma_7.npy')
p_smith_PV_deflation_gamma_7 = smith_PV_deflation_gamma_7[:, 0]
v_smith_PV_deflation_gamma_7 = smith_PV_deflation_gamma_7[:, 1]
smith_PV_deflation_gamma_16 = numpy.load('smith_PV_deflation_gamma_16.npy')
p_smith_PV_deflation_gamma_16 = smith_PV_deflation_gamma_16[:, 0]
v_smith_PV_deflation_gamma_16 = smith_PV_deflation_gamma_16[:, 1]
Loadings¶
In [ ]:
qois_filename = "Fig9-qois.dat"
res_basename = "Fig9"
seeds_filename = "Fig9.dat"
mesh_filebasename = "Fig9-mesh"
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_gamma_0 = {}
load_params_gamma_0["pf"] = 4
load_params_gamma_0["sigma_bar_00"] = 0.0
load_params_gamma_0["sigma_bar_11"] = 0.0
load_params_gamma_0["sigma_bar_01"] = 0.0
load_params_gamma_0["sigma_bar_10"] = 0.0
gamma = 7 * 1e-3
p_ini_7 = 0.3
p_fin_7 = 3.0
load_params_gamma_7 = {}
load_params_gamma_7["pf_lst"] = [p_ini_7, p_fin_7]
load_params_gamma_7["sigma_bar_00_lst"] = [0.0, 0.0]
load_params_gamma_7["sigma_bar_11_lst"] = [0.0, 0.0]
load_params_gamma_7["sigma_bar_01_lst"] = [0.0, 0.0]
load_params_gamma_7["sigma_bar_10_lst"] = [0.0, 0.0]
load_params_gamma_7["gamma_lst"] = [gamma, gamma]
step_params_gamma_7 = {}
step_params_gamma_7["n_steps"] = 2
step_params_gamma_7["Deltat"] = 1.
step_params_gamma_7["dt_ini"] = 0.01
step_params_gamma_7["dt_min"] = 0.0005
step_params_gamma_7["dt_max"] = 0.005
gamma = 16 * 1e-3
p_ini_16 = 1
p_fin_16 = 3.0
load_params_gamma_16 = {}
load_params_gamma_16["pf_lst"] = [p_ini_16, p_fin_16]
load_params_gamma_16["sigma_bar_00_lst"] = [0.0, 0.0]
load_params_gamma_16["sigma_bar_11_lst"] = [0.0, 0.0]
load_params_gamma_16["sigma_bar_01_lst"] = [0.0, 0.0]
load_params_gamma_16["sigma_bar_10_lst"] = [0.0, 0.0]
load_params_gamma_16["gamma_lst"] = [gamma, gamma]
step_params_gamma_16 = {}
step_params_gamma_16["n_steps"] = 2
step_params_gamma_16["Deltat"] = 1.
step_params_gamma_16["dt_ini"] = 0.01
step_params_gamma_16["dt_min"] = 0.0005
step_params_gamma_16["dt_max"] = 0.005
In [ ]:
def diameter_identifier(exp_data, diameter_ratio_initial, bnds, vf_asym):
exp_data = exp_data[(exp_data[:, 0] >= 4.5) & (exp_data[:, 0] <= 15)]
p_exp = exp_data[:, 0]
v_exp = exp_data[:, 1]
def J_cost(diameter_ratio_dimless):
diameter_ratio = diameter_ratio_initial[0]*(1 + diameter_ratio_dimless[0])
print("diameter_ratio: " +str(diameter_ratio))
print("diameter_ratio: " +str(diameter_ratio))
print("diameter_ratio: " +str(diameter_ratio))
print("diameter_ratio: " +str(diameter_ratio))
print("diameter_ratio: " +str(diameter_ratio))
print("diameter_ratio: " +str(diameter_ratio))
domain_y = 0.1 * diameter_ratio
domain_x = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.012 * diameter_ratio
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)
gamma = 7 * 1e-3
p_ini = 0.4
p_fin = 1.5
load_params_gamma_7 = {}
load_params_gamma_7["pf_lst"] = [p_ini, p_fin]
load_params_gamma_7["sigma_bar_00_lst"] = [0.0, 0.0]
load_params_gamma_7["sigma_bar_11_lst"] = [0.0, 0.0]
load_params_gamma_7["sigma_bar_01_lst"] = [0.0, 0.0]
load_params_gamma_7["sigma_bar_10_lst"] = [0.0, 0.0]
load_params_gamma_7["gamma_lst"] = [gamma, gamma]
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_7,
step_params=step_params_gamma_7,
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]
index = numpy.argmin(numpy.abs(pf_lst - p_ini*10.20), axis=None)
vf_lst_gamma_7 = vf_lst[index:]
pf_lst_gamma_7 = pf_lst[index:]
model_interpolator = interp1d(pf_lst_gamma_7, vf_lst_gamma_7, kind='cubic')
JC = 0
for i in range(len(exp_data)):
JC += (v_exp[i] - model_interpolator(p_exp[i]))**2
print("JC: " +str(JC))
return JC
diameter_ratio_dimless = 0
diameter_ratio = (scipy.optimize.minimize(J_cost, x0=diameter_ratio_dimless, method="L-BFGS-B", bounds=bnds, tol=1e-2)).x
diameter_ratio = diameter_ratio_initial[0]*(1 + diameter_ratio[0])
return diameter_ratio
domain_y = 0.1
domain_x = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.012
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
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_0,
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]
params_initial = [0.8]
bnds = [(-0.2, 6)]
diameter_ratio = diameter_identifier(numpy.vstack((smith_PV_deflation_gamma_7, smith_PV_inflation_gamma_7)), params_initial, bnds, vf_asym)
In [ ]:
diameter_ratio = 0.8034860511809191
$D_{alv} = 54~\mu m$¶
Defiding geometry¶
In [ ]:
domain_y = 0.1 * 0.8034860511809191
domain_x = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.012 * 0.8034860511809191
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
Model response¶
In [ ]:
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_0,
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]
vf_lst_D_54_gamma_0 = vf_lst
pf_lst_D_54_gamma_0 = pf_lst
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_7,
step_params=step_params_gamma_7,
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]
index = numpy.argmin(numpy.abs(pf_lst - p_ini_7*10.20), axis=None)
vf_lst_D_54_gamma_7 = vf_lst[index:]
pf_lst_D_54_gamma_7 = pf_lst[index:]
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_16,
step_params=step_params_gamma_16,
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]
index = numpy.argmin(numpy.abs(pf_lst - p_ini_16*10.20), axis=None)
vf_lst_D_54_gamma_16 = vf_lst[index:]
pf_lst_D_54_gamma_16 = pf_lst[index:]
$D_{alv} = 80~\mu m$¶
Defining geometry¶
In [ ]:
domain_y = 0.1 * 1.2
domain_x = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.012 * 1.2
# shift_y = -0.02
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
Model response¶
In [ ]:
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_0,
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]
vf_lst_D_80_gamma_0 = vf_lst
pf_lst_D_80_gamma_0 = pf_lst
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_7,
step_params=step_params_gamma_7,
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]
index = numpy.argmin(numpy.abs(pf_lst - p_ini_7*10.20), axis=None)
vf_lst_D_80_gamma_7 = vf_lst[index:]
pf_lst_D_80_gamma_7 = pf_lst[index:]
gamma = 16 * 1e-3
p_ini_16 = 0.5
p_fin_16 = 3.0
load_params_gamma_16 = {}
load_params_gamma_16["pf_lst"] = [p_ini_16, p_fin_16]
load_params_gamma_16["sigma_bar_00_lst"] = [0.0, 0.0]
load_params_gamma_16["sigma_bar_11_lst"] = [0.0, 0.0]
load_params_gamma_16["sigma_bar_01_lst"] = [0.0, 0.0]
load_params_gamma_16["sigma_bar_10_lst"] = [0.0, 0.0]
load_params_gamma_16["gamma_lst"] = [gamma, gamma]
step_params_gamma_16 = {}
step_params_gamma_16["n_steps"] = 2
step_params_gamma_16["Deltat"] = 1.
step_params_gamma_16["dt_ini"] = 0.01
step_params_gamma_16["dt_min"] = 0.0005
step_params_gamma_16["dt_max"] = 0.005
dmech.run_HollowBox_MicroPoroHyperelasticity(
dim=2,
mesh=mesh,
mat_params=mat_params,
load_params=load_params_gamma_16,
step_params=step_params_gamma_16,
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]
index = numpy.argmin(numpy.abs(pf_lst - p_ini_16*10.20), axis=None)
vf_lst_D_80_gamma_16 = vf_lst[index:]
pf_lst_D_80_gamma_16 = pf_lst[index:]
Plot¶
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_inflation_gamma_0, v_smith_PV_inflation_gamma_0, '#D94801', label='[Smith, 1986], inflation, $\gamma = 0~dyn/cm$')
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_7, v_smith_PV_inflation_gamma_7, '#FD8D3C', label='[Smith, 1986], inflation, $\gamma = 7~dyn/cm$')
plt.plot(p_smith_PV_deflation_gamma_7, v_smith_PV_deflation_gamma_7, '#FD8D3C', label='[Smith, 1986], deflation, $\gamma = 7~dyn/cm$', linestyle='dotted')
plt.plot(p_smith_PV_inflation_gamma_16, v_smith_PV_inflation_gamma_16, '#FDAE6B', label='[Smith, 1986], inflation, $\gamma = 16~dyn/cm$')
plt.plot(p_smith_PV_deflation_gamma_16, v_smith_PV_deflation_gamma_16, '#FDAE6B', label='[Smith, 1986], deflation, $\gamma = 16~dyn/cm$', linestyle='dotted')
plt.plot(pf_lst_D_54_gamma_0, vf_lst_D_54_gamma_0, '#084594', label='Model, $D_{alv} = 54~\mu m$, $\gamma = 0~dyn/cm$', linewidth=2.5)
plt.plot(pf_lst_D_54_gamma_7, vf_lst_D_54_gamma_7,'#4292C6' , label='Model, $D_{alv} = 54~\mu m$, $\gamma = 7~dyn/cm$', linewidth=2.5)
plt.plot(pf_lst_D_54_gamma_16, vf_lst_D_54_gamma_16,'#9ECAE1' , label='Model, $D_{alv} = 54~\mu m$, $\gamma = 16~dyn/cm$', linewidth=2.5)
plt.plot(pf_lst_D_80_gamma_0, vf_lst_D_80_gamma_0, '#005A32', label='Model, $D_{alv} = 80~\mu m$, $\gamma = 0~dyn/cm$', linewidth=2.5)
plt.plot(pf_lst_D_80_gamma_7, vf_lst_D_80_gamma_7,'#41AB5D' , label='Model, $D_{alv} = 80~\mu m$, $\gamma = 7~dyn/cm$', linewidth=2.5)
plt.plot(pf_lst_D_80_gamma_16, vf_lst_D_80_gamma_16,'#A1D99B' , label='Model, $D_{alv} = 80~\mu m$, $\gamma = 16~dyn/cm$', linewidth=2.5)
plt.xlim(0, 25)
plt.legend(loc = 'lower right', fontsize=12, bbox_to_anchor=(1.82, 0.), shadow=True)
plt.savefig('RVE_size_surface_tension.pdf',bbox_inches='tight')
plt.savefig("RVE_size_surface_tension.tif", format="tiff", dpi=300, bbox_inches="tight")