Hello there,
I’m working with brittle fracture formulation, specifically I’m basing my studies from this papper: https://www.researchgate.net/publication/338897618_Phase_field_fracture_implementation_in_FEniCS. I changed a bit the code, my problem is basically, my problem is the following:
from __future__ import print_function
from dolfin import *
from mshr import *
from dolfin import File
import matplotlib.pyplot as plt
import numpy as np
from dolfin import plot
domain = Rectangle(Point(-0.5, -0.5), Point(0.5, 0.5)) - Rectangle(Point(-0.5, -0.0001), Point(0, 0.0001))
unitary_circle1 = Circle(Point(0.0, 0.0), 0.020)
unitary_circle2 = Circle(Point(0.5, 0.2), 0.09)
unitary_circle3 = Circle(Point(0.1, 0.4), 0.025)
unitary_circle4 = Circle(Point(0.26, 0.71), 0.05)
domain_with_hole = domain - unitary_circle1 - unitary_circle2 - unitary_circle3 - unitary_circle4
mesh = generate_mesh(domain_with_hole, 50)
mesh = refine(mesh)
plot(mesh)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Mesh Plot')
plt.show()
# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
# Introduce manually the material parameters
Gc = 2.7
l = 0.015
lmbda = 121.1538e3
mu = 80.7692e3
# Constituive functions
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
def psi(u):
return 0.5*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)
# Boundary conditions
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
def Crack(x):
return abs(x[1]) < 1e-03 and x[0] <= 0.0
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
And the output:
But what I tried to do is basically, in this unitary_circles be a “hole” on the plate and in these circles, be a region where the function doesn’t act on the material, and have zero boundary conditions near of these “holes”. Also, I thought that doing it would change the properties of crack propagation, but is really similar with the one without holes. What I did wrong? Thank by your attention.
