Hey @conpierce8 !
i tried to combine what you said with this code Periodic homogenization but i still don’t get a perfectly symmetrical effective tensor.
[[ 2.73e+05 1.04e+05 -1.48e-05]
[ 1.04e+05 2.72e+05 -2.63e-05]
[ 7.93e-06 1.41e-05 8.42e+04]]
This is my geometry named 4cells.geo :
Point(1) = {0, 0, 0, 1.0};
Point(2) = {2, 0, 0, 1.0};
Point(3) = {0, 2, 0, 1.0};
Point(4) = {2, 2, 0, 1.0};
Point(5) = {2, 1, 0, 1.0};
Point(6) = {0, 1, 0, 1.0};
Point(7) = {1, 0, 0, 1.0};
Point(8) = {1, 1, 0, 1.0};
Point(9) = {1, 2, 0, 1.0};
//+
Line(1) = {1, 7};
Line(2) = {7, 2};
Line(3) = {6, 8};
Line(4) = {8, 5};
Line(5) = {3, 9};
Line(6) = {9, 4};
Line(7) = {1, 6};
Line(8) = {6, 3};
Line(9) = {7, 8};
Line(10) = {8, 9};
Line(11) = {2, 5};
Line(12) = {5, 4};
//+
Curve Loop(1) = {7, 3, -9, -1};
Plane Surface(1) = {1};
Curve Loop(2) = {2, 11, -4, -9};
Plane Surface(2) = {2};
Curve Loop(3) = {3, 10, -5, -8};
Plane Surface(3) = {3};
Curve Loop(4) = {4, 12, -6, -10};
Plane Surface(4) = {4};
//+
Physical Curve(1) = {7, 8};
Physical Curve(2) = {1, 2};
Physical Curve(3) = {11, 12};
Physical Curve(4) = {6, 5};
//+
Physical Surface(0) = {1};
Physical Surface(1) = {2};
Physical Surface(2) = {3};
Physical Surface(3) = {4};
And this is the fenics code
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Define geometry
fname = "4cells"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
L = 2.
vol = 4.
tol1 =DOLFIN_EPS
class PeriodicBoundary(SubDomain):
def inside(self, y, on_boundary):
#left and front boundary are target domain
return bool((near(y[0],0) and on_boundary) \
or (near(y[1],0,tol1) and on_boundary)) \
and (not(near(y[0],L) or near(y[1],L)))
def map(self,x,y):
#map right and back boundary to target domain
if near(x[0], L) and near(x[1], L):
# Added this condition to deal with the
# case x = y = L
y[0] = x[0] - L
y[1] = x[1] - L
elif near(x[0], L):
y[0] = x[0] - L
y[1] = x[1]
elif near(x[1], L):
y[0] = x[0]
y[1] = x[1] - L
else:
y[0] = -0.5
y[1] = -0.5
#define materials
Em, num = 220e3 , 0.25
Er, nur = 210e3 , 0.3
material_parameters = [(Em, num), (Er, nur),(Em, num), (Er, nur)]
nphases = len(material_parameters)
def eps(v):
return sym(grad(v))
def sigma(v, i, Eps):
E, nu = material_parameters[i]
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
return lmbda*tr(eps(v) + Eps)*Identity(2) + 2*mu*(eps(v)+Eps)
Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary())
V = FunctionSpace(mesh, Ve)
v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
dx = Measure('dx')(subdomain_data=subdomains)
Eps = Constant(((0, 0), (0, 0)))
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
def macro_strain(i):
"""returns the macroscopic strain for the 3 elementary load cases"""
Eps_Voigt = np.zeros((3,))
Eps_Voigt[i] = 1
return np.array([[Eps_Voigt[0], Eps_Voigt[2]/2.],
[Eps_Voigt[2]/2., Eps_Voigt[1]]])
def stress2Voigt(s):
return as_vector([s[0,0], s[1,1], s[0,1]])
Chom = np.zeros((3, 3))
for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
print("Solving {} case...".format(case))
Eps.assign(Constant(macro_strain(j)))
solve(a == L, w, [], solver_parameters={"linear_solver": "cg"})
(v, lamb) = split(w)
Sigma = np.zeros((3,))
for k in range(3):
Sigma[k] = assemble(sum([stress2Voigt(sigma(v, i, Eps))[k]*dx(i) for i in range(nphases)]))/vol
Chom[j, :] = Sigma
print(np.array_str(Chom, precision=2))
lmbda_hom = Chom[0, 1]
mu_hom = Chom[2, 2]
print(Chom[0, 0], lmbda_hom + 2*mu_hom)
E_hom = mu_hom*(3*lmbda_hom + 2*mu_hom)/(lmbda_hom + mu_hom)
nu_hom = lmbda_hom/(lmbda_hom + mu_hom)/2
print("Apparent Young modulus:", E_hom)
print("Apparent Poisson ratio:", nu_hom)
i also noticed that with a model of an 8x8 cells, when i defined the 64 materials for each cell with a Gaussian distribution like this:
E, nu = 210e3 , 0.3
random_E = np.random.normal(E,2,64)
random_nu = np.random.normal(nu,0.03,64)
the code only works if the standard deviation is small (0.03 in this case). and it’s also not symmetrical:
[[-6.46e-03 -3.64e-01 -1.49e+00]
[-1.93e-01 -3.94e-03 1.57e+00]
[ 7.17e-01 -1.94e+00 2.45e-04]]
otherwise it gives me this error : # Solve nonlinear variational problem