From an old, wise engineering book I've got the following problem: 
There is 2D infinite periodic structure with round holes of radius r with tension applied to it. The relation of maximum stress in the body to the applied one should follow nice 1/x like curve, see the picture. While I'm trying to reproduce this example, I'm getting anything but it.
I use displacement based weak formulation for elasticity with periodic boundary conditions remaping top to bottom and left to right, with exception of corners (0,0), (0,1), (1,0), (1,1).
For problem to be well-posed I've fixed displacement vector in one point (0,0) to be equal (0,0).
Am I right? Is fixing one point enough? I still can rotate the unit square though. Do I need to fix an additional point?
The tension is applied via Neumann boundary conditions with tension applied in the following way:
g = Expression(('1.0','0.0'), degree = 2)
dot(g,v)*ds(2) + dot(-g,v)*ds(1)
This is the place I'm most doubted about applying tension correctly. The rest of the code follows (it produces almost constant max stress equal to 0.4).
from __future__ import print_function
import pickle
from dolfin import *
from dolfin_adjoint import *
n = 64
mesh = UnitSquareMesh(n, n)
p = Constant(5)
eps = Constant(1.0e-2)
class HoleExpression(Expression):
def eval(self, values, x):
values[0] = 1.0
if (x[0]-0.5)**2 + (x[1]-0.5)**2 < rr**2:
values[0] = 0
return values
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((near(x[0],0) or near(x[1], 0)) and not ((near(x[0], 0) and (near(x[1], 0))) or (near(x[0], 0) and (near(x[1], 1))) or (near(x[0], 1) and (near(x[1], 0))) or (near(x[0], 1) and (near(x[1], 1))) ))
def map(self, x, y):
if near(x[0], 1):
y[0] = x[0] - 1.
y[1] = x[1]
elif near(x[1], 1):
y[0] = x[0]
y[1] = x[1] - 1.
else:
y[0] = -1000
y[1] = -1000
cor1 = CompiledSubDomain("near(x[0], 0.0) && near(x[1], 0.0)")
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
right.mark(boundary_parts, 1)
left.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]
norm = FacetNormal(mesh)
tang = as_vector([norm[1], -norm[0]])
g = Expression(('1.0','0.0'), degree = 2)# traction
A = FunctionSpace(mesh, "CG", 1)
PP = VectorFunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())
PN = VectorFunctionSpace(mesh, "CG", 1)
bc1 = DirichletBC(PP, Constant((0.0, 0.0)), cor1, method='pointwise')
bcs = [bc1]
parameters["std_out_all_processes"] = False
def k(a):
return eps + (1 - eps) * a**p
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
def sigma(rho, v):
return 2.0*mu*rho*sym(grad(v)) + lmbda*rho*tr(sym(grad(v)))*Identity(len(v))
def get_u_star(a):
u = TrialFunction(PP)
v = TestFunction(PP)
frm = inner(sigma(k(a), u), sym(grad(v)))*dx
rhs = dot(g,v)*ds(2) + dot(-g,v)*ds(1)
u = Function(PP)
problem = LinearVariationalProblem(frm, rhs, u, bcs)
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] ="bicgstab"
solver.solve()
return u
if __name__ == "__main__":
_sigma_max = []
ligament_eff = []
TTT = TensorFunctionSpace(mesh, "CG", 1)
for _i in range(10):
rr = 0.05*_i
le = 1-2*rr
ligament_eff.append(le)
a = interpolate(HoleExpression(degree = 2), A, name="Control")
u_star = get_u_star(a)
sg = sigma(k(a), u_star)
sigma_TTT = project(sg, TTT)
_sm = max(sigma_TTT.vector().array())
_sigma_max.append(_sm)
print(rr, le, _sm)
pickle.dump([_sigma_max, ligament_eff], open("data.pkl", "w"))
