The corresponding files can be obtained from:

Cohesive zone modelling of crack propagationΒΆ

In [153]:
from dolfin import *
from mshr import *
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt


domain = Rectangle(Point(0, 0), Point(60, 30)) - Circle(Point(20, 10), 5, 40) - Circle(Point(40, 20), 5, 40)
mesh = generate_mesh(domain, 50)
plot(mesh)

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 60)
Left().mark(facets, 1)
Right().mark(facets, 2)
dS = Measure("dS", domain=mesh, subdomain_data=facets, metadata={"quadrature_degree": 3})
ds = Measure("ds", domain=mesh, subdomain_data=facets, metadata={"quadrature_degree": 3})
dx = Measure("dx", domain=mesh, metadata={"quadrature_degree": 3})


n = FacetNormal(mesh)
t = as_vector([-n[1], n[0]])
../../_images/demo_cohesive_zone_cohesive_zone_crack_1_0.png
In [154]:
V_u = VectorFunctionSpace(mesh, "DG", 1)
V_cont = VectorFunctionSpace(mesh, "CG", 1)
V_d = FunctionSpace(mesh, "CR", 1)
In [155]:
E = 3.09e3
nu = 0.35
lmbda = Constant(E*nu/(1+nu)/(1-2*nu))
mu = Constant(E/2/(1+nu))

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2*mu*eps(v)
In [174]:
sig_max = 75.
tau_max = 75.
GcI = 0.3
GcII = 0.3
kappa = GcII/GcI
beta = 1.
delta_0n = GcI/sig_max/exp(1)
delta_0t = GcII/tau_max/sqrt(exp(1)/2.)
# compression penalty
K = 1e10

ppos = lambda x: (x+abs(x))/2.

def effective_opening(v, n):
    #A = outer(n, n) + beta**2/kappa*(Identity(2)-outer(n,n))
    #return sqrt(dot(v, dot(A, v)))
    return sqrt(dot(v, v))
def That(v, n):
    vn = dot(v,n)*n
    vt = v-vn
    return beta**2/kappa*vt + vn
def T(v, w_max, n):
    w = effective_opening(v, n)
    #return GcI/delta_0n**2*exp(-w_max/delta_0n)*That(v, n)
    return  GcI/delta_0n**2*exp(-w_max/delta_0n)*v

def damage(w_max):
    d = w_max/delta_0n
    return 1-(1+d+d**2/2.)*exp(-d)
In [202]:
u = Function(V_u, name="Displacement")
u_, du = TestFunction(V_u), TrialFunction(V_u)
d = Function(V_d, name="Damage")
w_max = Function(V_d, name="Maximum opening")

DU = Expression(("t", "0"), t=0, degree=0)
bc = [DirichletBC(V_u, Constant((0., 0.)), Left(), method="geometric"),
      DirichletBC(V_u, DU, Right(), method="geometric")]

gamma = 10*(3*lmbda+2*mu)
h = CellVolume(mesh)
ih = 2/(h("+")+h("-"))
ih = 1/FacetArea(mesh)

bulk_energy = 0.5*inner(sigma(u), eps(u))*dx
interfacial_energy = 0.5*inner(T(jump(u), avg(w_max), n("+")), jump(u))*dS - \
                          inner(T(u, w_max, n("+")), u)*ds(1)  - \
                          inner(T(DU-u, w_max, n("+")), DU-u)*ds(2)  #- \
                          #inner(T(jump(u), avg(w_max), n("+")) - avg(dot(sigma(u),n)), jump(u))*dS + \
                          #inner(T(u, w_max, n("+")) - dot(sigma(u),n), u)*ds(1)  + \
                          #inner(T(DU-u, w_max, n("+")) - dot(sigma(u),n), DU-u)*ds(2)
stabilization = gamma*ih("+")*dot(jump(u), jump(u))*dS #+ ih*gamma*dot(u, u)*ds(1) + ih*gamma*dot(u-DU, u-DU)*ds(2)
total_energy = bulk_energy + interfacial_energy #+ stabilization

D_total_energy = derivative(total_energy, u, u_)
F = replace(D_total_energy, {u: du})
D2_total_energy = derivative(D_total_energy, u, du)
assemble(D_total_energy)
Out[202]:
<dolfin.cpp.la.Vector; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Vector > *' at 0x7f588f0be330> >
In [205]:

w_max = Function(V_d)
w_max.interpolate(Expression("fabs(x[0]-30) <= 1 ? d : 0", d=delta_0n, degree=1))
#p = plot(project(GcI/delta_0n**2*exp(-w_max/delta_0n), FunctionSpace(mesh, "DG", 0)))
#plt.colorbar(p)
solve(lhs(F) == rhs(F), u, bc)
plot(u[0])
--- Instant: compiling ---
Out[205]:
<matplotlib.tri.tricontour.TriContourSet at 0x7f588a997390>
../../_images/demo_cohesive_zone_cohesive_zone_crack_6_2.png
In [191]:
def facet_project(v, V):
    v_, dv = TestFunction(V), TrialFunction(V)
    a_form = inner(avg(v_), avg(dv))*dS #+ inner(v_, dv)*ds(1)
    L_form = inner(avg(v_), v)*dS #+ inner(v_, v_bd)*ds(1)
    u  = Function(V)
    A = assemble(a_form, keep_diagonal=True)
    A.ident_zeros() # Regularize the matrix
    b = assemble(L_form)
    solve(A, u.vector(), b)
    return u

problem = NonlinearVariationalProblem(D_total_energy, u, bc, D2_total_energy)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-4
prm['newton_solver']['relative_tolerance'] = 1E-6
prm['newton_solver']['maximum_iterations'] = 25

w_max = Function(V_d)
Nincr = 20
Force = np.zeros((Nincr+1,))
loading = np.linspace(0, 5., Nincr+1)
for (i, t) in enumerate(loading[1:]):
    DU.t = t
    #solve(lhs(F) == rhs(F), u, bc)
    solver.solve()
    #solve(D2_total_energy == rhs(F), u, bc, J=D2_total_energy)
    w_max.assign(facet_project(sqrt(dot(jump(u), jump(u))), V_d))
    d.assign(facet_project(damage(avg(w_max)), V_d))
    print np.mean(d.vector().get_local())
    print w_max.vector().norm("l2")
    Force[i+1] = assemble(sigma(u)[0,0]*ds(2))
    p = plot(d)
    plt.colorbar(p)
    plt.show()

plt.plot(loading, Force)
plt.show()
                  Calling FFC just-in-time (JIT) compiler, this may take some time.
1.3195907205498413e-06
0.00233082911417
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_1.png
1.022675825638074e-05
0.00465877015097
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_3.png
3.357114987097635e-05
0.00698823501118
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_5.png
7.733587984695249e-05
0.00931754030195
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_7.png
0.00014746064085835057
0.0116601659365
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_9.png
0.00024724374190114806
0.0139764700224
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_11.png
0.00038262850736624725
0.016305385396
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_13.png
0.0005552091763995347
0.0186466329133
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_15.png
0.0007712545273155069
0.0209713916367
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_17.png
0.0010294470566391872
0.0232976044461
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_19.png
0.0013342551467160624
0.0256248753411
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_21.png
0.0016836634166018313
0.0279529400447
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_23.png
0.0020958685597487745
0.0302816154723
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_25.png
0.0025518803596576134
0.0326107707921
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_27.png
0.0030611804119767562
0.0349403100339
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_29.png
0.0036330316307105247
0.0372701612078
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_31.png
0.004251887379701253
0.0396002692569
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_33.png
0.004926634877869133
0.041930591357
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_35.png
0.005657796217323516
0.0442686898544
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_37.png
0.0064457478960519435
0.0465952088922
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_39.png
../../_images/demo_cohesive_zone_cohesive_zone_crack_7_40.png
In [142]:
delta_0n
Out[142]:
0.0014715177646857694
In [9]:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-9-78bbf88282cd> in <module>()
      1 W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=1, quad_scheme='default')
----> 2 W0 = FunctionSpace(mesh, W0e, restriction="facet")

/usr/lib/python2.7/dist-packages/dolfin/functions/functionspace.pyc in __init__(self, *args, **kwargs)
    195                                  + str(args[0]))
    196             elif len(args) == 2:
--> 197                 self._init_from_ufl(*args, **kwargs)
    198             else:
    199                 self._init_convenience(*args, **kwargs)

TypeError: _init_from_ufl() got an unexpected keyword argument 'restriction'
In [ ]:
GcI/delta_0n**2