Figure 10

Figure 10#

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 = "Fig10-seeds.dat"
mesh_filebasename = "Fig10-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_params = {"E":Es, "nu":nus}
mat_params_c = {"model":"CGNH"    , "parameters":mat_params_params}
mat_params_d = {"model":"CGNH_bar", "parameters":mat_params_params}

Loading#

beta_min = 1.
beta_max = 2.
n_steps = 10
beta_lst = numpy.linspace(beta_min, beta_max, n_steps)

Symbolic expression of the macroscopic model response#

Symbolic variables#

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

K_macro = lmbda_macro + 2*mu_macro/3
G_macro = mu_macro

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")

Coupled model expression#

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

Sigma_macro_c  = 2*sympy.diff(W_c, C_macro)
Sigma_macro_c -= pf * J_macro * C_macro.inv()

p_hydro_macro_c = - sympy.trace(Sigma_macro_c.T * C_macro)/3/J_macro

Decoupled model expression#

W_d  = (K_macro/4) * (J_macro**2 - 1 - 2*sympy.ln(J_macro))
W_d += (G_macro/2) * (J_macro**(-2/3) * I_C_macro - 3)

Sigma_macro_d  = 2*sympy.diff(W_d, C_macro)
Sigma_macro_d -= pf * J_macro * C_macro.inv()

p_hydro_macro_d = - sympy.trace(Sigma_macro_d.T * C_macro)/3/J_macro

Computing macroscopic model response#

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

p_hydro_macro_c = p_hydro_macro_c.subs(list(zip(C_macro, C))).subs(pf, 0.) # list & zip should not be needed, cf. https://github.com/sympy/sympy/issues/10589

p_hydro_macro_d = p_hydro_macro_d.subs(list(zip(C_macro, C))).subs(pf, 0.) # list & zip should not be needed, cf. https://github.com/sympy/sympy/issues/10589

homogenization_problem = dmech.HomogenizationProblem(
    dim=2,
    mesh=mesh,
    mat_params=mat_params_params,
    vol=mesh_bbox_V0,
    bbox=mesh_bbox)
lmbda_tilde, mu_tilde = homogenization_problem.get_lambda_and_mu()

p_hydro_macro_c = p_hydro_macro_c.subs(lmbda_macro, lmbda_tilde).subs(mu_macro, mu_tilde)
p_hydro_macro_d = p_hydro_macro_d.subs(lmbda_macro, lmbda_tilde).subs(mu_macro, mu_tilde)

beta_macro_vals = numpy.linspace(beta_min, beta_max, 100)

p_hydro_macro_c_ = sympy.lambdify(beta, p_hydro_macro_c, modules=["numpy"])
p_hydro_macro_c_vals = p_hydro_macro_c_(beta_macro_vals)

p_hydro_macro_d_ = sympy.lambdify(beta, p_hydro_macro_d, modules=["numpy"])
p_hydro_macro_d_vals = p_hydro_macro_d_(beta_macro_vals)

Computing microscopic model response#

load_params = {}
load_params["pf_lst"] = [0.]*len(beta_lst)
load_params["U_bar_00_lst"] = [float(beta - 1) for beta in beta_lst]
load_params["sigma_bar_01_lst"] = [0.]*len(beta_lst)
load_params["sigma_bar_10_lst"] = [0.]*len(beta_lst)
load_params["U_bar_11_lst"] = [float(1./beta-1.) for beta in beta_lst]

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

Coupled model#

res_basename = "Fig10"
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_c,
    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 = qois_vals[:, qois_name_list.index("U_bar_XX")]
beta_micro_c_vals = eps_lst+1
p_hydro_micro_c_vals = qois_vals[:, qois_name_list.index("p_hydro")]

Decoupled model#

res_basename = "Fig10"
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_d,
    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 = qois_vals[:, qois_name_list.index("U_bar_XX")]
beta_micro_d_vals = eps_lst+1
p_hydro_micro_d_vals = qois_vals[:, qois_name_list.index("p_hydro")]

Generating plot#

plt.figure()

plt.rc("xtick", labelsize=14)
plt.rc("ytick", labelsize=14)
plt.rc("legend", fontsize=12)

plt.xlabel(r"$\beta$"          , fontsize=16)
plt.ylabel(r"$p_{hydro}~(kPa)$", fontsize=16)

plt.plot(beta_micro_c_vals, p_hydro_micro_c_vals, "#084594"                    )
plt.plot(beta_micro_d_vals, p_hydro_micro_d_vals, "#084594", linestyle="dashed")

plt.plot(beta_macro_vals, p_hydro_macro_c_vals, "#99000D"                    )
plt.plot(beta_macro_vals, p_hydro_macro_d_vals, "#99000D", linestyle="dashed")

plt.legend(
    ["Microscopic model " r"$\tilde{p}_{hydro}$"" - coupled",\
     "Microscopic model " r"$\tilde{p}_{hydro}$"" - decoupled",\
     "Macroscopic model " r"$\bar{p}_{hydro}$"" - coupled",\
     "Macroscopic model " r"$\bar{p}_{hydro}$"" - decoupled"])

plt.xlim(beta_min, beta_max)
plt.ylim(-0.1, +0.01)

plt.savefig("Fig10-p_hydro.pdf", bbox_inches="tight")
plt.show()
Contents