Figure 8

Figure 8#

Imports#

import dolfin
import matplotlib
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 = "Fig8-seeds.dat"
mesh_filebasename = "Fig8-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_bbox_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#

nus = 0.2
mat_params_E1 = {"model":"CGNH", "parameters":{"E":0.5, "nu":nus}}
mat_params_E2 = {"model":"CGNH", "parameters":{"E":1.0, "nu":nus}}
mat_params_E3 = {"model":"CGNH", "parameters":{"E":2.0, "nu":nus}}

Es = 1.
mat_params_nu1 = {"model":"CGNH", "parameters":{"E":Es, "nu":0.1}}
mat_params_nu2 = {"model":"CGNH", "parameters":{"E":Es, "nu":0.2}}
mat_params_nu3 = {"model":"CGNH", "parameters":{"E":Es, "nu":0.4}}

Loading#

epsilon_max = 0.5

Symbolic expression of the macroscopic model response#

lmbda_macro, mu_macro = sympy.symbols("lambda, mu")

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_macro/4) * (J_macro**2 - 1 - 2*sympy.ln(J_macro))
W_skel += (mu_macro/2) * (I_C_macro - 2 - 2*sympy.ln(J_macro))

pf = sympy.symbols("pf")

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

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

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 = Sigma_macro.subs(pf, 0.)
sigma_macro = F * Sigma_macro * F.T / J

Sensitivity to solid Young modulus#

Computing macroscopic model response#

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_E1["parameters"],
    vertices=mesh_bbox_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde_E1, mu_tilde_E1 = homogenization_problem.get_lambda_and_mu()
sigma_macro_E1 = sigma_macro.subs(lmbda_macro, lmbda_tilde_E1).subs(mu_macro, mu_tilde_E1)

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_E2["parameters"],
    vertices=mesh_bbox_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde_E2, mu_tilde_E2 = homogenization_problem.get_lambda_and_mu()
sigma_macro_E2 = sigma_macro.subs(lmbda_macro, lmbda_tilde_E2).subs(mu_macro, mu_tilde_E2)

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_E3["parameters"],
    vertices=mesh_bbox_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde_E3, mu_tilde_E3 = homogenization_problem.get_lambda_and_mu()
sigma_macro_E3 = sigma_macro.subs(lmbda_macro, lmbda_tilde_E3).subs(mu_macro, mu_tilde_E3)

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

sigma_macro_E1_ = sympy.lambdify(eps, sigma_macro_E1[0,0], modules=["numpy"])
sigma_macro_vals_E1 = sigma_macro_E1_(epsilon_vals)

sigma_macro_E2_ = sympy.lambdify(eps, sigma_macro_E2[0,0], modules=["numpy"])
sigma_macro_vals_E2 = sigma_macro_E2_(epsilon_vals)

sigma_macro_E3_ = sympy.lambdify(eps, sigma_macro_E3[0,0], modules=["numpy"])
sigma_macro_vals_E3 = sigma_macro_E3_(epsilon_vals)

Computing 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 = "Fig8"
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_E1,
    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_glob_E1 = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_glob_E1 = qois_vals[:, qois_name_list.index("sigma_bar_XX")]

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_E2,
    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_glob_E2 = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_glob_E2 = qois_vals[:, qois_name_list.index("sigma_bar_XX")]

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_E3,
    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_glob_E3 = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_glob_E3 = qois_vals[:, qois_name_list.index("sigma_bar_XX")]

Generating plot#

# From https://stackoverflow.com/a/67870930

# define an object that will be used by the legend
class MulticolorPatch(object):
    def __init__(self, colors):
        self.colors = colors

# define a handler for the MulticolorPatch object
class MulticolorPatchHandler(object):
    def legend_artist(self, legend, orig_handle, fontsize, handlebox):
        width, height = handlebox.width, handlebox.height
        patches = []
        for i, c in enumerate(orig_handle.colors):
            patches.append(
                plt.Rectangle([width/len(orig_handle.colors)*i-handlebox.xdescent, -handlebox.ydescent],
                width/len(orig_handle.colors),
                height,
                facecolor=c,
                edgecolor="none"))
        patch = matplotlib.collections.PatchCollection(patches, match_original=True)
        handlebox.add_artist(patch)
        return patch

# ------ choose some colors
colors1 = ["#084594", "#99000D"]
colors2 = ["#4292C6", "#EF3B2C"]
colors3 = ["#9ECAE1", "#FC9272"]

# ------ create a dummy-plot (just to show that it works)
f, ax = plt.subplots()
ax.plot(eps_glob_E1, sigma_glob_E1, c=colors1[0])
ax.plot(eps_glob_E2, sigma_glob_E2, c=colors2[0])
ax.plot(eps_glob_E3, sigma_glob_E3, c=colors3[0])
ax.plot(epsilon_vals, sigma_macro_vals_E1, c=colors1[1])
ax.plot(epsilon_vals, sigma_macro_vals_E2, c=colors2[1])
ax.plot(epsilon_vals, sigma_macro_vals_E3, c=colors3[1])

# ------ get the legend-entries that are already attached to the axis
h, l = ax.get_legend_handles_labels()

# ------ append the multicolor legend patches
h.append(MulticolorPatch(colors1))
l.append(r"$E_s = 0.5~kPa$")

h.append(MulticolorPatch(colors2))
l.append(r"$E_s = 1~kPa$")

h.append(MulticolorPatch(colors3))
l.append(r"$E_s = 2~kPa$")

# ------ create the legend
f.legend(h, l, loc="upper left", handler_map={MulticolorPatch: MulticolorPatchHandler()}, bbox_to_anchor=(0.125,0.875))

plt.xlim(0., epsilon_max)
plt.ylim(0., 0.175)
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~(kPa)$", fontsize=16)
plt.savefig("Fig8-Es.pdf", bbox_inches="tight")
plt.show()

Sensitivity to solid Poisson coefficient#

Computing macroscopic model response#

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_nu1["parameters"],
    vertices=mesh_bbox_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde_nu1, mu_tilde_nu1 = homogenization_problem.get_lambda_and_mu()
sigma_macro_nu1 = sigma_macro.subs(lmbda_macro, lmbda_tilde_nu1).subs(mu_macro, mu_tilde_nu1)

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_nu2["parameters"],
    vertices=mesh_bbox_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde_nu2, mu_tilde_nu2 = homogenization_problem.get_lambda_and_mu()
sigma_macro_nu2 = sigma_macro.subs(lmbda_macro, lmbda_tilde_nu2).subs(mu_macro, mu_tilde_nu2)

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_nu3["parameters"],
    vertices=mesh_bbox_vertices,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde_nu3, mu_tilde_nu3 = homogenization_problem.get_lambda_and_mu()
sigma_macro_nu3 = sigma_macro.subs(lmbda_macro, lmbda_tilde_nu3).subs(mu_macro, mu_tilde_nu3)

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

sigma_macro_nu1_ = sympy.lambdify(eps, sigma_macro_nu1[0,0], modules=["numpy"])
sigma_macro_vals_nu1 = sigma_macro_nu1_(epsilon_vals)

sigma_macro_nu2_ = sympy.lambdify(eps, sigma_macro_nu2[0,0], modules=["numpy"])
sigma_macro_vals_nu2 = sigma_macro_nu2_(epsilon_vals)

sigma_macro_nu3_ = sympy.lambdify(eps, sigma_macro_nu3[0,0], modules=["numpy"])
sigma_macro_vals_nu3 = sigma_macro_nu3_(epsilon_vals)

Computing 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 = "Fig8"
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_nu1,
    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_glob_nu1 = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_glob_nu1 = qois_vals[:, qois_name_list.index("sigma_bar_XX")]

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_nu2,
    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_glob_nu2 = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_glob_nu2 = qois_vals[:, qois_name_list.index("sigma_bar_XX")]

dmech.run_HollowBox_MicroPoroHyperelasticity(
    dim=2,
    mesh=mesh,
    displacement_perturbation_degree=2,
    quadrature_degree=4,
    mat_params=mat_params_nu3,
    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_glob_nu3 = qois_vals[:, qois_name_list.index("U_bar_XX")]
sigma_glob_nu3 = qois_vals[:, qois_name_list.index("sigma_bar_XX")]

Generating plot#

# From https://stackoverflow.com/a/67870930

# ------ create a dummy-plot (just to show that it works)
f, ax = plt.subplots()
ax.plot(eps_glob_nu1, sigma_glob_nu1, c=colors1[0])
ax.plot(eps_glob_nu2, sigma_glob_nu2, c=colors2[0])
ax.plot(eps_glob_nu3, sigma_glob_nu3, c=colors3[0])
ax.plot(epsilon_vals, sigma_macro_vals_nu1, c=colors1[1])
ax.plot(epsilon_vals, sigma_macro_vals_nu2, c=colors2[1])
ax.plot(epsilon_vals, sigma_macro_vals_nu3, c=colors3[1])

# ------ get the legend-entries that are already attached to the axis
h, l = ax.get_legend_handles_labels()

# ------ append the multicolor legend patches
h.append(MulticolorPatch(colors1))
l.append(r"$\nu_s = 0.1$")

h.append(MulticolorPatch(colors2))
l.append(r"$\nu_s = 0.2$")

h.append(MulticolorPatch(colors3))
l.append(r"$\nu_s = 0.4$")

# ------ create the legend
f.legend(h, l, loc="upper left", handler_map={MulticolorPatch: MulticolorPatchHandler()}, bbox_to_anchor=(0.125,0.875))

plt.xlim(0., epsilon_max)
plt.ylim(0., 0.175)
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~(kPa)$", fontsize=16)
plt.savefig("Fig8-Es.pdf", bbox_inches="tight")
plt.show()

Figure 8

Contents

Figure 8#

Imports#

Parameters#

Geometry#

Solid material properties#

Loading#

Symbolic expression of the macroscopic model response#

Sensitivity to solid Young modulus#

Computing macroscopic model response#

Computing microscopic model response#

Generating plot#

Sensitivity to solid Poisson coefficient#

Computing macroscopic model response#

Computing microscopic model response#

Generating plot#