Figure 11

Figure 11#

Imports#

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

import dolfin_mech                    as dmech
import micro_poro_structure_generator as gen

Parameters#

Geometry#

seeds_filename = "Fig11-seeds.dat"
mesh_filebasename = "Fig11-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 = [mesh_xmin, mesh_xmax, mesh_ymin, mesh_ymax]
mesh_bbox_V0 = (mesh_xmax-mesh_xmin) * (mesh_ymax-mesh_ymin)

Solid material parameters#

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

Loading#

beta = 1.25
pf_lst = [0.0, 0.2]

Symbolic expression of the macroscopic model response#

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

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

pf = sympy.symbols("p_f")

W_skel  = (lmbda_macro/4) * (J_macro**2 - 1 - 2*sympy.ln(J_macro))
W_skel += (mu_macro/2) * (I_C_macro - 3 - 2*sympy.ln(J_macro))

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

p_hydro_macro = - sympy.trace(Sigma_macro.T * C_macro)/3/J_macro

Sigma_D_macro = Sigma_macro + p_hydro_macro * J_macro * C_macro.inv()
Sigma_VM_macro = sympy.sqrt(3/2 * sympy.trace(Sigma_D_macro.T * Sigma_D_macro))

Computing macroscopic model response#

F = sympy.Matrix(
    [[beta, 0     , 0],\
     [0   , 1/beta, 0],\
     [0   , 0     , 1]])
J = F.det()
C = F.T * F

p_hydro_macro = p_hydro_macro.subs(list(zip(C_macro, C)))
Sigma_VM_macro = Sigma_VM_macro.subs(list(zip(C_macro, C)))

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

p_hydro_macro = p_hydro_macro.subs(lmbda_macro, lmbda_tilde).subs(mu_macro, mu_tilde)
Sigma_VM_macro = Sigma_VM_macro.subs(lmbda_macro, lmbda_tilde).subs(mu_macro, mu_tilde)

p_hydro_macro_lst = []
Sigma_VM_macro_lst = []
for pf_ in pf_lst:
    p_hydro_macro_lst.append(p_hydro_macro.subs(pf, pf_))
    Sigma_VM_macro_lst.append(Sigma_VM_macro.subs(pf, pf_))
p_hydro_macro_lst

Computing microscopic model response#

p_hydro_micro_lst = []
Sigma_VM_micro_lst = []

for pf_ in pf_lst:
    load_params = {}
    load_params["U_bar_00"] = beta-1
    load_params["sigma_bar_01"] = 0.
    load_params["sigma_bar_10"] = 0.
    load_params["U_bar_11"] = 1/beta-1
    load_params["pf"] = pf_

    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 = "Fig11"
    qois_filename = res_basename+"-qois.dat"

    problem = 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,
        add_p_hydro_and_Sigma_VM_FoI=True,
        verbose=1)

    p_hydro_micro_lst.append(problem.get_foi("p_hydro").func.vector().get_local())
    Sigma_VM_micro_lst.append(problem.get_foi("Sigma_VM").func.vector().get_local())

Generating plots#

Hydrostatic pressure distribution#

pf_min = -0.3
pf_max = +0.3
for (pf, p_hydro_micro, p_hydro_macro) in zip(pf_lst, p_hydro_micro_lst, p_hydro_macro_lst):
    plt.figure()
    plt.rc("xtick", labelsize=14)
    plt.rc("ytick", labelsize=14)
    plt.rc("legend", fontsize=12)
    plt.xlabel("Hydrostatic pressure (kPa)", fontsize=16)
    plt.ylabel("Frequency distribution"    , fontsize=16)
    seaborn.histplot(
        p_hydro_micro,
        bins=51,
        binrange=[pf_min, pf_max],
        color="#084594",
        stat="density",
        kde=True,
        kde_kws=dict(cut=3, bw_adjust=0.9),
        alpha=0.4,
        edgecolor=(1, 1, 1, 0.4))
    plt.axvline(
        x=p_hydro_macro,
        color="#99000D",
        label="axvline - full height")
    plt.legend(
        ["Microscopic model",\
        "Macroscopic model"])
    plt.xlim(pf_min, pf_max)
    plt.ylim(0., 17.)
    plt.savefig("Fig11-p_hydro_density_pf="+str(pf)+".pdf", bbox_inches="tight")
    plt.show()

Von Mises stress distribution#

Sigma_VM_min = 0.0
Sigma_VM_max = 0.5
for (pf, Sigma_VM_micro, Sigma_VM_macro) in zip(pf_lst, Sigma_VM_micro_lst, Sigma_VM_macro_lst):
    plt.figure()
    plt.rc("xtick", labelsize=14)
    plt.rc("ytick", labelsize=14)
    plt.rc("legend", fontsize=12)
    plt.xlabel("Von Mises stress (kPa)", fontsize=16)
    plt.ylabel("Frequency distribution", fontsize=16)
    seaborn.histplot(
        Sigma_VM_micro,
        bins=51,
        binrange=[Sigma_VM_min, Sigma_VM_max],
        color="#084594",
        stat="density",
        kde=True,
        kde_kws=dict(cut=0, bw_adjust=0.6),
        alpha=0.4,
        edgecolor=(1, 1, 1, 0.4))
    plt.axvline(
        x=Sigma_VM_macro,
        color="#99000D",
        label="axvline - full height")
    plt.legend(
        ["Microscopic model",\
        "Macroscopic model"])
    plt.xlim(Sigma_VM_min, Sigma_VM_max)
    plt.ylim(0., 20.)
    plt.savefig("Fig11-Sigma_VM_density_pf="+str(pf)+".pdf", bbox_inches="tight")
    plt.show()
Contents