]'s demos! - Home](../_static/logo.png){kind=logo}
Figure 9#
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 = "Fig9-seeds.dat"
mesh_filebasename = "Fig9-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)
dV = dolfin.Measure("dx", domain=mesh)
Vs0 = dolfin.assemble(dolfin.Constant(1.) * dV)
Phis0 = Vs0/mesh_bbox_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}}
Loading#
epsilon_max = 0.5
pf_lst = [0., 0.25, 0.5, 1.]
Quasi-incompressible solid#
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))
Sigma_ter_macro = 2*sympy.diff(W_skel, C_macro)
eps = sympy.symbols("epsilon")
epsilon = sympy.Matrix(
[[eps, 0 ],\
[ 0 , eps]])
F = sympy.eye(2) + epsilon
J = F.det()
C = F.T * F
Sigma_ter_macro = Sigma_ter_macro.subs(list(zip(C_macro, C))) # 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_incomp["parameters"],
vertices=mesh_bbox_vertices,
vol=mesh_bbox_V0,
bbox=mesh_bbox)
lmbda_tilde, mu_tilde = homogenization_problem.get_lambda_and_mu()
Sigma_ter_macro_incomp = Sigma_ter_macro.subs(lmbda_macro, lmbda_tilde).subs(mu_macro, mu_tilde)
epsilon_macro_vals = numpy.linspace(0., epsilon_max, 100)
Sigma_ter_macro_incomp_ = sympy.lambdify(eps, Sigma_ter_macro_incomp[0,0], modules=["numpy"])
Sigma_ter_macro_incomp_vals = Sigma_ter_macro_incomp_(epsilon_macro_vals)
Microscopic model response#
eps_lst = []
Sigma_ter_lst = []
for k_pf, pf in enumerate(pf_lst):
load_params = {}
load_params["U_bar_00_lst"] = [0., epsilon_max]
load_params["sigma_bar_01_lst"] = [0., 0.]
load_params["sigma_bar_10_lst"] = [0., 0.]
load_params["U_bar_11_lst"] = [0., epsilon_max]
load_params["pf_lst"] = [pf, pf]
step_params = {}
step_params["n_steps"] = 2
step_params["Deltat"] = 1.
step_params["dt_ini"] = 0.1
step_params["dt_min"] = 0.01
step_params["dt_max"] = 0.1
res_basename = "Fig9"
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_xx_lst = qois_vals[:, qois_name_list.index("U_bar_XX")]
eps_yy_lst = qois_vals[:, qois_name_list.index("U_bar_YY")]
sigma_xx_lst = qois_vals[:, qois_name_list.index("sigma_bar_XX")]
sigma_yy_lst = qois_vals[:, qois_name_list.index("sigma_bar_YY")]
p_lst = qois_vals[:, qois_name_list.index("p_f")]
eps_lst.append(eps_xx_lst)
Sigma_ter_lst.append(numpy.empty(len(eps_xx_lst)))
for i in range(len(eps_xx_lst)):
U_bar = numpy.diag([eps_xx_lst[i], eps_yy_lst[i]])
F_bar = numpy.eye(2) + U_bar
J_bar = numpy.linalg.det(F_bar)
C_bar = F_bar.T * F_bar
sigma = numpy.diag([sigma_xx_lst[i], sigma_yy_lst[i]])
Sigma = J_bar * numpy.dot(numpy.linalg.inv(F_bar), numpy.dot(sigma, numpy.linalg.inv(F_bar.T)))
Sigma_ter_lst[-1][i] = (Sigma + p_lst[i] * J_bar * numpy.linalg.inv(C_bar))[0,0]
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_{T}~(kPa)$", fontsize=16)
plt.plot(epsilon_macro_vals, Sigma_ter_macro_incomp_vals, "#99000D")
plt.plot(eps_lst[0], Sigma_ter_lst[0], "#084594")
plt.plot(eps_lst[1], Sigma_ter_lst[1], "#2171B5")
plt.plot(eps_lst[2], Sigma_ter_lst[2], "#4292C6")
plt.plot(eps_lst[3], Sigma_ter_lst[3], "#6BAED6")
plt.legend(["Macroscopic model "r"$\bar{\Sigma}_{T}$"", any"r"$~\bar{p}_f$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=0.00~(kPa)$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=0.25~(kPa)$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=0.5~(kPa)$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=1.0~(kPa)$"])
plt.xlim(0., epsilon_max)
plt.ylim(0., 0.12)
plt.savefig("Fig9-incompressible.pdf", bbox_inches="tight")
plt.show()
Compressible solid#
Microscopic model response#
eps_lst = []
Sigma_ter_lst = []
for k_pf, pf in enumerate(pf_lst):
load_params = {}
load_params["U_bar_00_lst"] = [0., epsilon_max]
load_params["sigma_bar_01_lst"] = [0., 0.]
load_params["sigma_bar_10_lst"] = [0., 0.]
load_params["U_bar_11_lst"] = [0., epsilon_max]
load_params["pf_lst"] = [pf, pf]
step_params = {}
step_params["n_steps"] = 2
step_params["Deltat"] = 1.
step_params["dt_ini"] = 0.1
step_params["dt_min"] = 0.01
step_params["dt_max"] = 0.1
res_basename = "Fig9"
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_xx_lst = qois_vals[:, qois_name_list.index("U_bar_XX")]
eps_yy_lst = qois_vals[:, qois_name_list.index("U_bar_YY")]
sigma_xx_lst = qois_vals[:, qois_name_list.index("sigma_bar_XX")]
sigma_yy_lst = qois_vals[:, qois_name_list.index("sigma_bar_YY")]
p_lst = qois_vals[:, qois_name_list.index("p_f")]
eps_lst.append(eps_xx_lst)
Sigma_ter_lst.append(numpy.empty(len(eps_xx_lst)))
for i in range(len(eps_xx_lst)):
U_bar = numpy.diag([eps_xx_lst[i], eps_yy_lst[i]])
F_bar = numpy.eye(2) + U_bar
J_bar = numpy.linalg.det(F_bar)
C_bar = F_bar.T * F_bar
sigma = numpy.diag([sigma_xx_lst[i], sigma_yy_lst[i]])
Sigma = J_bar * numpy.dot(numpy.linalg.inv(F_bar), numpy.dot(sigma, numpy.linalg.inv(F_bar.T)))
Sigma_ter_lst[-1][i] = (Sigma + p_lst[i] * J_bar * numpy.linalg.inv(C_bar))[0,0]
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_{T}~(kPa)$", fontsize=16)
plt.plot(eps_lst[0], Sigma_ter_lst[0], "#084594")
plt.plot(eps_lst[1], Sigma_ter_lst[1], "#2171B5")
plt.plot(eps_lst[2], Sigma_ter_lst[2], "#4292C6")
plt.plot(eps_lst[3], Sigma_ter_lst[3], "#6BAED6")
plt.legend([r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=0.00~(kPa)$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=0.25~(kPa)$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=0.5~(kPa)$", r"Microscopic model $\tilde{\Sigma}_{T}$, $\tilde{p}_f=1.0~(kPa)$"])
plt.xlim(0., epsilon_max)
plt.ylim(0., 0.12)
plt.savefig("Fig9-compressible.pdf", bbox_inches="tight")
plt.show()