Phase field method in DCB test

General
1

Hi everyone,

I am reproducing a DCB test with phase field method with fracture and I am getting problems with the application of the displacements. Both (upper left and lower left corners) should have equal magnitude but different sign.

After running the code and visualizing the displacements I can see how the left area of the bar deforms downwards. Without expressing logic with the displacements of the DCB test.

I have tried assigning a variable “load_value” in the while loop in case the definition of “t” is causing problems by leaving both loads in the same direction. But I could not correct the problem either.

# Initial load value
load_value = 0.0

# Staggered scheme
while t <= 1.0:
t += deltaT
if t >= 0.7:
deltaT = 0.01
load_value += t * u_r
load_top.t = load_value
load_bot.t = -load_value
iter = 0
err = 1

How can I correct this possible error in the code structure?

Thanks in advance

# Preliminaries and mesh
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import math

# Mesh
vertices=[Point(0.035,0),Point(0,-0.00001),Point(0,-0.0021),Point(0.15,-0.0021),Point(0.15,0.0021),Point(0,0.0021),Point(0,0.00001)]
domain=Polygon(vertices)
mesh=generate_mesh(domain, 250)
plot(mesh)
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 J/m2
Gc =  98.2
l= (0.0042/16)*5
#lambda GPa
lmbda = 15.4838e9
#shear modulus GPa
mu = 4.26e9

# 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
# Right boundary
class Right(SubDomain):
      def inside(self, x, on_boundary):
          return near(x[0], 0.15)

class point(SubDomain):
      def inside(self, x, on_boundary):
          return near(x[0],0) and near(x[1],0.0021, 1e-2)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0) and near(x[1], -0.0021, 1e-2)

#Boundary segments
right = Right()
pt = point()
pt_bot = Bottom()

def Crack(x):
    return abs(x[1]) < 1e-06 and x[0] <= 0.035
load_top = Expression("t", t = 0.0, degree=1)
load_bot = Expression("-t", t = 0.0, degree=1)

bctop= DirichletBC(W.sub(1), load_top, pt, method="pointwise")
bcbot= DirichletBC(W.sub(1), load_bot, pt_bot, method="pointwise")
bcright = DirichletBC(W, Constant((0,0)), right)
bc_u = [bctop, bcbot, bcright]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
pt.mark(boundaries,1)

ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)

# Variational form
unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
E_du = ((1.0-pold)**2)*inner(grad(v),sigma(u))*dx
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
            *inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew, bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)

# Initialization of the iterative procedure and output requests
t = 0

#displacement 10 mm
u_r = 0.01
deltaT  = 0.1
tol = 1e-3
conc_f = File ("./ResultsDir_DCB/phi.pvd")
fname = open('ForcevsDisp_DCB.txt', 'w')

# Staggered scheme
while t<=1.0:
    t += deltaT
    if t >=0.7:
        deltaT = 0.01
    #adding loads top & bot
    load_top.t=t*u_r
    load_bot.t=t*u_r
    iter = 0
    err = 1

    while err > tol:
        iter += 1
        solver_disp.solve()
        solver_phi.solve()
        err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
        err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
        err = max(err_u,err_phi)
        
        uold.assign(unew)
        pold.assign(pnew)
        Hold.assign(project(psi(unew), WW))

        if err < tol:
		
            print ('Iterations:', iter, ', Total time', t)

            if round(t*1e4) % 10 == 0:
                conc_f << pnew

                Traction = dot(sigma(unew),n)
                fy = Traction[1]*ds(1)
                fname.write(str(t*u_r) + "\t")
                fname.write(str(assemble(fy)) + "\n")
	    	    
fname.close()
print ('Simulation completed') 

#Deformação
eps=TensorFunctionSpace(mesh,"CG", degree=1)
deformação=Function(eps)
deformação.assign(project(epsilon(unew),eps))
#File('Bend_Test/deformacion.pvd')<<deformação
File('ResultsDir_DCB/deformacion.pvd')<<deformação

#Tensão
Vsig = TensorFunctionSpace(mesh, "CG", degree=1)
tensões = Function(Vsig)
tensões.assign(project(sigma(unew),Vsig))
File('ResultsDir_DCB/tensiones.pvd')<<tensões

#Deslocamento
File('ResultsDir_DCB/desplazamiento.pvd')<<unew