]'s demos! - Home](../_static/logo.png){kind=logo}
Figure 5#
Imports#
import dolfin
import matplotlib.pyplot as plt
import math
import numpy
import dolfin_mech as dmech
import micro_poro_structure_generator as gen
Parameters#
Es = 1
nus = 0.2
mus = Es/2/(1 + nus)
lmbdas = Es*nus/(1 + nus)/(1 - 2*nus)
ks = (3*lmbdas + 2*mus)/3
mat_params = {"E":Es, "nu":nus}
Computing homogenized parameters for various geometries#
Voronoi inclusions#
lmbda_lst_voronoi = []
mu_lst_voronoi = []
phi_lst_voronoi = []
kappa_lst_voronoi = []
mesh_filebasename = "Fig5-voronoi"
seeds_filename = "Fig5-voronoi-seeds.dat"
domain_y = 1.
domain_x = domain_y * numpy.sqrt(3)/1.5
thickness_lst = [0.01, 0.04, 0.06, 0.09, 0.12, 0.15, 0.17]
for k_thickness, thickness in enumerate(thickness_lst):
print(k_thickness)
Phif0 = 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=False)
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)
homogenization_problem = dmech.HomogenizationProblem(
dim=2,
mesh=mesh,
mat_params=mat_params,
vol=mesh_bbox_V0,
bbox=mesh_bbox)
lmbda_tilde, mu_tilde = homogenization_problem.get_lambda_and_mu()
kappa_tilde = homogenization_problem.get_kappa()
phi_lst_voronoi.append(Phif0)
lmbda_lst_voronoi.append(lmbda_tilde)
mu_lst_voronoi.append(mu_tilde)
kappa_lst_voronoi.append(kappa_tilde)
Hexagonal inclusions#
lmbda_lst_hexagonal = []
mu_lst_hexagonal = []
phi_lst_hexagonal = []
kappa_lst_hexagonal = []
mesh_filebasename = "Fig5-hexagonal"
seeds_filename = "Fig5-hexagonal-seeds.dat"
domain_y = 1.
domain_x = domain_y * numpy.sqrt(3)/1.5/2
gen.generate_seeds_semi_regular(
DoI=0.,
row=1,
domain_y=domain_y,
seeds_filename=seeds_filename)
for k_thickness in range(11):
thickness = (k_thickness+0.5)/20
Phif0 = 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=False)
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)
homogenization_problem = dmech.HomogenizationProblem(
dim=2,
mesh=mesh,
mat_params=mat_params,
vol=mesh_bbox_V0,
bbox=mesh_bbox)
lmbda_tilde, mu_tilde = homogenization_problem.get_lambda_and_mu()
kappa_tilde = homogenization_problem.get_kappa()
phi_lst_hexagonal.append(Phif0)
lmbda_lst_hexagonal.append(lmbda_tilde)
mu_lst_hexagonal.append(mu_tilde)
kappa_lst_hexagonal.append(kappa_tilde)
Circular inclusions#
lmbda_lst_circular = []
mu_lst_circular = []
phi_lst_circular = []
kappa_lst_circular = []
mesh_filebasename = "Fig5-circular"
width = 1.
V0 = math.sqrt(3) * width**2
for Phif0 in numpy.linspace(0.01, 0.85, 42):
r = math.sqrt(Phif0*V0/math.pi/2)
gen.generate_mesh_2D_rectangle_w_circular_inclusions(
mesh_filename=mesh_filebasename,
width=width,
r=r,
lcar=width/40,
shift_x=0.,
shift_y=0.)
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)
homogenization_problem = dmech.HomogenizationProblem(
dim=2,
mesh=mesh,
mat_params=mat_params,
vol=mesh_bbox_V0,
bbox=mesh_bbox)
lmbda_tilde, mu_tilde = homogenization_problem.get_lambda_and_mu()
kappa_ = homogenization_problem.get_kappa()
phi_lst_circular.append(Phif0)
lmbda_lst_circular.append(lmbda_tilde)
mu_lst_circular.append(mu_tilde)
kappa_lst_circular.append(kappa_tilde)
Post-processing#
Adding limit porosities#
k_lst_voronoi = []
k_lst_hexagonal = []
k_lst_circular = []
phi_lst_voronoi.insert(0, 1.)
mu_lst_voronoi.insert(0, 0.)
lmbda_lst_voronoi.insert(0, 0.)
kappa_lst_voronoi.insert(0, 0.)
phi_lst_voronoi.append(0.)
mu_lst_voronoi.append(mus)
lmbda_lst_voronoi.append(lmbdas)
kappa_lst_voronoi.append(ks)
for i in range(len(phi_lst_voronoi)):
k_lst_voronoi.append((3*lmbda_lst_voronoi[i] + 2*mu_lst_voronoi[i])/3)
phi_lst_hexagonal.insert(0, 1.)
mu_lst_hexagonal.insert(0, 0.)
lmbda_lst_hexagonal.insert(0, 0.)
kappa_lst_hexagonal.insert(0, 0.)
phi_lst_hexagonal.append(0.)
mu_lst_hexagonal.append(mus)
lmbda_lst_hexagonal.append(lmbdas)
kappa_lst_hexagonal.append(ks)
for i in range (len(phi_lst_hexagonal)):
k_lst_hexagonal.append((3*lmbda_lst_hexagonal[i] + 2*mu_lst_hexagonal[i])/3)
phi_lst_circular.insert(0, 0.)
mu_lst_circular.insert(0, mus)
lmbda_lst_circular.insert(0, lmbdas)
kappa_lst_circular.insert(0, ks)
phi_lst_circular.append(1.)
mu_lst_circular.append(0.)
lmbda_lst_circular.append(0.)
kappa_lst_circular.append(0.)
for i in range (len(phi_lst_circular)):
k_lst_circular.append((3*lmbda_lst_circular[i] + 2*mu_lst_circular[i])/3)
Computing homogenized parameters for the dilute and differential schemes [Dormieux, 2006]#
def Eq1(mu_hom, phi):
return ((1+4*mus/(3*ks))*(mu_hom/mus)**3)/(2 - (1 - 4*mus/(3*ks))*(mu_hom/mus)**(3/5)) - (1 - phi)**6
def Eq2(mu_hom, k_hom):
return mu_hom/mus - ((-4/3*mu_hom/k_hom + 1)**(5/3))/((-4/3*mus/ks + 1)**(5/3))
def differential_scheme(phi_0):
h = 10e-6
mu_hom = mus
k_hom = ks
res = 1
while (res > 10e-7):
res = Eq1(mu_hom, phi_0)
jac = (Eq1(mu_hom + h, phi_0) - Eq1(mu_hom, phi_0))/h
delta = - res/jac
mu_hom = mu_hom + delta
res = 1
while (res > 10e-7):
res = Eq2(mu_hom, k_hom)
jac = (Eq2(mu_hom, k_hom + h) - Eq2(mu_hom, k_hom))/h
delta = - res/jac
k_hom = k_hom + delta
return mu_hom, k_hom
def k_diff_func(mu_hom):
return 4/3*mu_hom/(1 - (mu_hom/mus)**(3/5)*(1 - 4/3*mus/ks))
mu_diff = []
k_diff = []
mu_dil = []
k_dil = []
for phi in phi_lst_hexagonal:
mu_tilde, k_tilde = differential_scheme(phi)
mu_diff.append(mu_tilde)
k_diff.append(k_diff_func(mu_tilde))
mu_dil.append(mus*(1 - 5*phi*(3*ks + 4*mus)/(9*ks + 8*mus)))
k_dil.append(ks*(1 - phi*(1 + 3*ks/4/mus)))
Normalizing homogenized parameters#
mu_lst_norm_voronoi = []
k_lst_norm_voronoi = []
kappa_lst_norm_voronoi = []
mu_lst_hexagonal_norm = []
k_lst_hexagonal_norm = []
kappa_lst_hexagonal_norm = []
mu_lst_norm_circular = []
k_lst_norm_circular = []
kappa_lst_norm_circular = []
mu_diff_norm = []
k_diff_norm = []
mu_dil_norm = []
k_dil_norm = []
for i in range (len(mu_lst_hexagonal)):
mu_lst_hexagonal_norm.append(mu_lst_hexagonal[i]/mus)
mu_dil_norm.append(mu_dil[i]/mus)
mu_diff_norm.append(mu_diff[i]/mus)
k_lst_hexagonal_norm.append(k_lst_hexagonal[i]/ks)
k_diff_norm.append(k_diff[i]/ks)
k_dil_norm.append(k_dil[i]/ks)
kappa_lst_hexagonal_norm.append(kappa_lst_hexagonal[i]/ks)
for i in range (len(mu_lst_voronoi)):
mu_lst_norm_voronoi.append(mu_lst_voronoi[i]/mus)
k_lst_norm_voronoi.append(k_lst_voronoi[i]/ks)
kappa_lst_norm_voronoi.append(kappa_lst_voronoi[i]/ks)
for i in range (len(mu_lst_circular)):
mu_lst_norm_circular.append(mu_lst_circular[i]/mus)
k_lst_norm_circular.append(k_lst_circular[i]/ks)
kappa_lst_norm_circular.append(kappa_lst_circular[i]/ks)
Generating plots#
phi_dil = numpy.linspace(min(phi_lst_hexagonal),(9*ks + 8*mus)/(3*ks + 4*mus)/5,100)
mu_dil_norm = 1 - 5*phi_dil*(3*ks + 4*mus)/(9*ks + 8*mus)
plt.figure()
plt.rc("xtick", labelsize=14)
plt.rc("ytick", labelsize=14)
plt.rc("legend", fontsize=12)
plt.xlabel(r"$\Phi_{f0}$" , fontsize=16)
plt.ylabel(r"$\tilde{\mu}/\mu_s$", fontsize=16)
plt.plot(phi_lst_circular , mu_lst_norm_circular , "#084594" )
plt.plot(phi_lst_hexagonal, mu_lst_hexagonal_norm, "#084594", linestyle="dashdot")
plt.plot(phi_lst_voronoi , mu_lst_norm_voronoi , "#084594", linestyle="dashed" )
plt.plot(phi_lst_hexagonal, mu_diff_norm , "#005824" )
plt.plot(phi_dil , mu_dil_norm , "#005824", linestyle="dashed" )
plt.xlim(0., 1.)
plt.ylim(0., 1.)
plt.legend(
["Microscopic model - circular inclusions",\
"Microscopic model - hexagonal inclusions",\
"Microscopic model - polygonal inclusions",\
"Differential scheme [Dormieux et al. 2006]",\
"Dilute scheme [Dormieux et al. 2006]"])
plt.text(phi_lst_hexagonal[9]-0.020, mu_lst_hexagonal_norm[9]-0.08, "(a)", fontsize=12)
plt.text(phi_lst_hexagonal[6]-0.015, mu_lst_hexagonal_norm[6]-0.07, "(b)", fontsize=12)
plt.text(phi_lst_hexagonal[4]-0.025, mu_lst_hexagonal_norm[4]-0.04, "(c)", fontsize=12)
plt.text(phi_lst_hexagonal[2]-0.015, mu_lst_hexagonal_norm[2]+0.02, "(d)", fontsize=12)
plt.savefig("Fig5-mu.pdf", bbox_inches="tight")
plt.show()
phi_dil = numpy.linspace(min(phi_lst_hexagonal),4*mus/(4*mus + 3*ks),100)
k_dil_norm = (1 - phi_dil*(1 + 3*ks/4/mus))
plt.figure()
plt.rc("xtick", labelsize=14)
plt.rc("ytick", labelsize=14)
plt.rc("legend", fontsize=12)
plt.xlabel(r"$\Phi_{f0}$" , fontsize = 16)
plt.ylabel(r"$\tilde{k}/k_s$", fontsize = 16)
plt.plot(phi_lst_circular , k_lst_norm_circular , "#084594" )
plt.plot(phi_lst_hexagonal, k_lst_hexagonal_norm, "#084594", linestyle="dashdot")
plt.plot(phi_lst_voronoi , k_lst_norm_voronoi , "#084594", linestyle="dashed" )
plt.plot(phi_lst_hexagonal, k_diff_norm , "#005824" )
plt.plot(phi_dil , k_dil_norm , "#005824", linestyle="dashed" )
plt.xlim(0., 1.)
plt.ylim(0., 1.)
plt.legend(
["Microscopic model - circular inclusions",\
"Microscopic model - hexagonal inclusions",\
"Microscopic model - polygonal inclusions",\
"Differential scheme [Dormieux et al. 2006]",\
"Dilute scheme [Dormieux et al. 2006]"])
plt.text(phi_lst_hexagonal[9]+0.01 , k_lst_hexagonal_norm[9] , "(a)", fontsize=12)
plt.text(phi_lst_hexagonal[6] , k_lst_hexagonal_norm[6]+0.025, "(b)", fontsize=12)
plt.text(phi_lst_hexagonal[4]-0.015, k_lst_hexagonal_norm[4]+0.025, "(c)", fontsize=12)
plt.text(phi_lst_hexagonal[2]-0.015, k_lst_hexagonal_norm[2]+0.03 , "(d)", fontsize=12)
plt.savefig("Fig5-k.pdf", bbox_inches="tight")
plt.show()