Figure 7

Figure 7#

Imports#

import dolfin
import matplotlib.pyplot as plt
import numpy
import sympy
import sys

import dolfin_mech                    as dmech
import micro_poro_structure_generator as gen

Parameters#

Geometry#

seeds_filename = "Fig7-seeds.dat"
mesh_filebasename = "Fig7-mesh"

domain_y = 1.
domain_x = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.092

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)

coord = mesh.coordinates()
xmax = max(coord[:,0]); xmin = min(coord[:,0])
ymax = max(coord[:,1]); ymin = min(coord[:,1])
V0 = (xmax-xmin) * (ymax-ymin)
dV = dolfin.Measure("dx", domain=mesh)
Vs0 = dolfin.assemble(dolfin.Constant(1.) * dV)
Phis0 = Vs0/V0

Solid material properties#

Es = 1.
nus_incomp = 0.499
mat_params_incomp = {"model":"CGNH", "parameters":{"E":Es, "nu":nus_incomp}}
nus_comp = 0.2
mat_params_comp = {"model":"CGNH", "parameters":{"E":Es, "nu":nus_comp}}

Biaxial strain loading#

epsilon_max = 0.5

load_params = {}
load_params["pf"] = 0.0
load_params["U_bar_00"] = epsilon_max
load_params["sigma_bar_01"] = 0.0
load_params["sigma_bar_10"] = 0.0
load_params["U_bar_11"] = epsilon_max

step_params = {}
step_params["Deltat"] = 1.
step_params["dt_ini"] = 0.1
step_params["dt_min"] = 0.01
step_params["dt_max"] = 0.1

Quasi-incompressible solid#

res_basename = "Fig7"
qois_filename = res_basename+"-qois.dat"

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_incomp,
    load_params=load_params,
    step_params=step_params,
    res_basename=res_basename,
    verbose=1,
    write_qois_limited_precision=False)

qois_vals = numpy.loadtxt(qois_filename)
qois_name_list = open(qois_filename).readline().split()[1:]
eps_lst_incomp = qois_vals[:, qois_name_list.index("U_bar_XX")]
vs_lst_incomp = qois_vals[:, qois_name_list.index("vs")]
Phi_s_lst_incomp = [vs_/V0 for vs_ in vs_lst_incomp]

Compressible solid#

res_basename = "Fig7"
qois_filename = res_basename+"-qois.dat"

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_comp,
    load_params=load_params,
    step_params=step_params,
    res_basename=res_basename,
    verbose=1,
    write_qois_limited_precision=False)

qois_vals = numpy.loadtxt(qois_filename)
qois_name_list = open(qois_filename).readline().split()[1:]
eps_lst_comp = qois_vals[:, qois_name_list.index("U_bar_XX")]
vs_lst_comp = qois_vals[:, qois_name_list.index("vs")]
Phi_s_lst_comp = [vs_/V0 for vs_ in vs_lst_comp]

Plot#

Phi_s_macro_lst = [Phis0]*len(eps_lst_incomp) # Macroscopic model response

plt.figure()
plt.rc("xtick", labelsize=14)
plt.rc("ytick", labelsize=14)
plt.rc("legend", fontsize=12)
plt.xlabel(r"$E_x,~E_y~()$", fontsize=16)
plt.ylabel(r"$\Phi_s$", fontsize=16)
plt.plot(eps_lst_comp  , Phi_s_lst_comp  , "#084594")
plt.plot(eps_lst_incomp, Phi_s_lst_incomp, "#084594", linestyle="dashed")
plt.plot(eps_lst_incomp, Phi_s_macro_lst , "#99000D")
plt.legend(["Microscopic model - compressible solid", "Microscopic model - quasi-incompressible solid", "Macroscopic model"])
plt.xlim(0., epsilon_max)
plt.ylim(0.18, 0.42)
plt.savefig("Fig7-Phi-s_E.pdf", bbox_inches="tight")
plt.show()

Fluid pressure loading#

pf_max = 0.5

load_params = {}
load_params["U_bar_00"] = 0.
load_params["sigma_bar_01"] = 0.0
load_params["sigma_bar_10"] = 0.0
load_params["U_bar_11"] = 0.
load_params["pf"] = pf_max

step_params = {}
step_params["Deltat"] = 1.
step_params["dt_ini"] = 0.1
step_params["dt_min"] = 0.01
step_params["dt_max"] = 0.1

Quasi-incompressible solid#

res_basename = "Fig7"
qois_filename = res_basename+"-qois.dat"

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_incomp,
    load_params=load_params,
    step_params=step_params,
    res_basename=res_basename,
    verbose=1,
    write_qois_limited_precision=False)

qois_vals = numpy.loadtxt(qois_filename)
qois_name_list = open(qois_filename).readline().split()[1:]
pf_lst_incomp = qois_vals[:, qois_name_list.index("p_f")]
vs_lst_incomp = qois_vals[:, qois_name_list.index("vs")]
Phi_s_lst_incomp = [vs_/V0 for vs_ in vs_lst_incomp]

Compressible solid#

res_basename = "Fig7"
qois_filename = res_basename+"-qois.dat"

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_comp,
    load_params=load_params,
    step_params=step_params,
    res_basename=res_basename,
    verbose=1,
    write_qois_limited_precision=False)

qois_vals = numpy.loadtxt(qois_filename)
qois_name_list = open(qois_filename).readline().split()[1:]
pf_lst_comp = qois_vals[:, qois_name_list.index("p_f")]
vs_lst_comp = qois_vals[:, qois_name_list.index("vs")]
Phi_s_lst_comp = [vs_/V0 for vs_ in vs_lst_comp]

Plot#

Phi_s_macro_lst = [Phis0]*len(pf_lst_incomp) # Macroscopic model response

plt.figure()
plt.rc("xtick", labelsize=14)
plt.rc("ytick", labelsize=14)
plt.rc("legend", fontsize=12)
plt.xlabel(r"$p_f~(kPa)$", fontsize=16)
plt.ylabel(r"$\Phi_s$", fontsize=16)
plt.plot(pf_lst_comp  , Phi_s_lst_comp  , "#084594")
plt.plot(pf_lst_incomp, Phi_s_lst_incomp, "#084594", linestyle="dashed")
plt.plot(pf_lst_incomp, Phi_s_macro_lst , "#99000D")
plt.legend(["Microscopic model - compressible solid", "Microscopic model - quasi-incompressible solid", "Macroscopic model"])
plt.xlim(0., pf_max)
plt.ylim(0.18, 0.42)
plt.savefig("Fig7-Phi-s_p-f.pdf", bbox_inches="tight")
plt.show()
Contents