Figure 6

Figure 6#

Imports#

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

import dolfin_mech                    as dmech
import micro_poro_structure_generator as gen

Parameters#

Geometry#

seeds_filename = "Fig6-seeds.dat"
mesh_filebasename = "Fig6-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)

mesh_coord = mesh.coordinates()
mesh_xmax = max(mesh_coord[:,0]); mesh_xmin = min(mesh_coord[:,0])
mesh_ymax = max(mesh_coord[:,1]); mesh_ymin = min(mesh_coord[:,1])
mesh_vertices = numpy.array(
    [[mesh_xmin, mesh_ymin],
     [mesh_xmax, mesh_ymin],
     [mesh_xmax, mesh_ymax],
     [mesh_xmin, mesh_ymax]])
mesh_bbox = [mesh_xmin, mesh_xmax, mesh_ymin, mesh_ymax]
mesh_bbox_V0 = (mesh_xmax-mesh_xmin) * (mesh_ymax-mesh_ymin)

Solid material properties#

Es = 1.
nus = 0.499
mat_params = {"model":"CGNH", "parameters":{"E":Es, "nu":nus}}

Symbolic expressions#

Kinematical and pressure variables#

eps = sympy.symbols("eps")
epsilon = sympy.Matrix(
    [[eps,  0 ],\
     [ 0 , eps]])
F = sympy.eye(2) + epsilon
J = F.det()
C = F.T * F

pf = sympy.symbols("pf")

Linearized model response expression#

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params["parameters"],
    vertices=mesh_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde, mu_tilde = homogenization_problem.get_lambda_and_mu()
kappa_tilde = homogenization_problem.get_kappa()

sigma_lin = lmbda_tilde * sympy.trace(epsilon) * sympy.eye(2) + 2*mu_tilde * epsilon - pf * sympy.eye(2)

Macroscopic model response expression#

C_macro = sympy.MatrixSymbol("C", 2, 2).as_explicit()
I_C_macro = sympy.trace(C_macro)
III_C_macro = sympy.det(C_macro)
J_macro = sympy.sqrt(III_C_macro)

W_skel  = (lmbda_tilde/4) * (J_macro**2 - 1 - 2*sympy.ln(J_macro))
W_skel += (mu_tilde/2) * (I_C_macro - 2 - 2*sympy.ln(J_macro))

Sigma_macro  = 2*sympy.diff(W_skel, C_macro)
Sigma_macro -= pf * J_macro * C_macro.inv()

Sigma_macro = Sigma_macro.subs(list(zip(C_macro, C))) # list & zip should not be needed, cf. https://github.com/sympy/sympy/issues/10589
sigma_macro = F * Sigma_macro * F.T / J

Biaxial strain loading#

epsilon_max = 0.5

Linearized and macroscopic model responses#

epsilon_vals = numpy.linspace(0., epsilon_max, 100)

sigma_lin_ = sympy.lambdify(eps, sigma_lin.subs(pf, 0.)[0,0], modules=["numpy"])
sigma_lin_vals = sigma_lin_(epsilon_vals)
sigma_lin_vals /= Es

sigma_macro_ = sympy.lambdify(eps, sigma_macro.subs(pf, 0.)[0,0], modules=["numpy"])
sigma_macro_vals = sigma_macro_(epsilon_vals)
sigma_macro_vals /= Es

Microscopic model response#

load_params = {}
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
load_params["pf"] = 0.0

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

res_basename = "Fig6"
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,
    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:]
epsilon_micro_vals = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_micro_vals = qois_vals[:, qois_name_list.index("sigma_bar_XX")]
sigma_micro_vals /= Es

Plot#

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"$\sigma/E_s~()$", fontsize=16)
plt.plot(epsilon_vals      , sigma_lin_vals  , "#005824")
plt.plot(epsilon_micro_vals, sigma_micro_vals, "#084594")
plt.plot(epsilon_vals      , sigma_macro_vals, "#99000D")
plt.legend(
    ["Linear model "+ r"$\tilde{\sigma}_{xx},~\tilde{\sigma}_{yy}$",\
     "Microscopic model "+ r"$\tilde{\sigma}_{xx},~\tilde{\sigma}_{yy}$",\
     "Macroscopic model "+ r"$\bar{\sigma}_{xx},~\bar{\sigma}_{yy}$"])
plt.xlim(0., epsilon_max)
plt.ylim(0., 0.15)
plt.savefig("Fig6-eps.pdf", bbox_inches="tight")
plt.show()

Fluid pressure loading#

pf_max = 0.5

Linearized and macroscopic model responses#

pf_vals = numpy.linspace(0., pf_max, 100)

sigma_lin_ = sympy.lambdify(pf, sigma_lin.subs(eps, 0.)[0,0], modules=["numpy"])
sigma_lin_vals = sigma_lin_(pf_vals)
sigma_lin_vals /= Es

sigma_macro_ = sympy.lambdify(pf, sigma_macro.subs(eps, 0.)[0,0], modules=["numpy"])
sigma_macro_vals = sigma_macro_(pf_vals)
sigma_macro_vals /= Es

Microscopic model response#

load_params = {}
load_params["U_bar_00"] = 0.0
load_params["sigma_bar_01"] = 0.0
load_params["sigma_bar_10"] = 0.0
load_params["U_bar_11"] = 0.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

res_basename = "Fig6"
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,
    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_micro_vals = qois_vals[:, qois_name_list.index("p_f")]
sigma_micro_vals = qois_vals[:, qois_name_list.index("sigma_bar_XX")]
sigma_micro_vals /= Es

Plot#

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"$\sigma/{E_s}~()$", fontsize=16)
plt.plot(pf_vals      , sigma_lin_vals  , "#005824")
plt.plot(pf_micro_vals, sigma_micro_vals, "#084594")
plt.plot(pf_vals      , sigma_macro_vals, "#99000D")
plt.legend(
    ["Linear model "+ r"$\tilde{\sigma}_{xx},~\tilde{\sigma}_{yy}$",\
     "Microscopic model "+ r"$\tilde{\sigma}_{xx},~\tilde{\sigma}_{yy}$",\
     "Macroscopic model "+ r"$\bar{\sigma}_{xx},~\bar{\sigma}_{yy}$"])
plt.xlim(0., pf_max)
plt.ylim(-0.5, 0. )
plt.savefig("Fig6-pf.pdf", bbox_inches="tight")
plt.show()

Figure 6

Contents

Figure 6#

Imports#

Parameters#

Geometry#

Solid material properties#

Symbolic expressions#

Kinematical and pressure variables#

Linearized model response expression#

Macroscopic model response expression#

Biaxial strain loading#

Linearized and macroscopic model responses#

Microscopic model response#

Plot#

Fluid pressure loading#

Linearized and macroscopic model responses#

Microscopic model response#

Plot#