Implementing solution as Dirichlet BC

variational formulation
1

Hello everyone,

I am solving a phase-field model for crack growth problem. For that, I am applying the load in small time-steps, and the problem is quasi-static. Now, I want to define a growing crack region as Dirichlet BC, so that I can implement the region of crack (based on the solution of loading step i) as Dirichlet BC for solving the loading step i+1.
I have the simple code as below. In the below code, I want to know how I can extract the region where |p| > 0.95 (p is obtained in line 43 pnew.assign(p)) and use that as Dirichlet BC for future loading steps?

from dolfin import *
mesh = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5), 20, 20)
W = VectorFunctionSpace(mesh, 'CG', 1)
p, q, dp = Function(W), TestFunction(W), TrialFunction(W)
u, v, du = Function(W), TestFunction(W), TrialFunction(W)
unew, pold, pnew = Function(W), Function(W) , Function(W)
  
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")

def crack_region(x, on_boundary):
    return abs(x[1]) < 10e-03 and x[0] <= 0

bcbot = DirichletBC(W, Constant((0.0, 0)), bot)
bctop = DirichletBC(W, Constant((0.002, 0)), top)
bc_u = [bcbot, bctop]
bc_phi = [DirichletBC(W, Constant((0.0, 1)), crack_region)]

Gc, l, lmbda, mu=  1, 0.01, 100e3, 80e3
def W0(u):
    F = Identity(len(u)) + grad(u)
    C = F.T*F 
    Ic, J = tr(C), det(F)
    E = (mu/2)*(Ic - 2) - mu*ln(J) + (lmbda/2)*(ln(J))**2
    return E

def total_energy(u,d):
    E = ((1- conditional(lt(dot(d,d), 0.001),0.0, sqrt(dot(d,d))))**2)*W0(u) +\
        Gc* ( dot(d,d)/(2*l) + (l/2)*inner(grad(d), grad(d)) )  
    return E
      
E_du = derivative(total_energy(u, pold) * dx , u, v)   
E_phi = derivative(total_energy(unew, p) * dx, p, q) 
J_phi  = derivative(E_phi, p, dp)
J_u = derivative(E_du, u, du)    
p_disp = NonlinearVariationalProblem(E_du, u, bc_u, J_u)
p_phi = NonlinearVariationalProblem(E_phi, p, bc_phi ,J_phi)
solver_disp = NonlinearVariationalSolver(p_disp)
solver_phi = NonlinearVariationalSolver(p_phi)
solver_disp.solve()
unew.assign(u) 
solver_phi.solve()
pnew.assign(p)

Thank you so much for your time!
Mary

2

The usual condition used for irreversibility of crack growth in phase field models is an inequality constraint of p_{n+1} \geq p_{n} on the entire domain (i.e., the phase field cannot decrease anywhere in the domain), not an equality constraint on a p-dependent subset of the domain. You could implement this constraint using a PETScTAOSolver for the minimization of E_phi. Using a PETScTAOSolver with bound constraints defined by Functions is demonstrated (albeit in a different application context) by this example.

3 Likes
3

Thank you so much for your answer. Actually, based on my problem, I would like to apply the irreversibility on a p-dependent subset of the domain. So, I would be so thankful if you could recommend me any method to implement the grown crack region as Dirichlet BC for the next loading step.
Thank you very much!

4

You could probably do this within an augmented Lagrangian iteration, as discussed for irreversibility in phase field methods here. The method of the linked paper enforces the constraint on the full domain, but it could be adapted to limit constraint enforcement to a subset of the domain, by including a UFL conditional in the penalty term.

1 Like