Phase field modelling for a notched plate with hole

Hello everyone

I am trying to model an experiment from a paper, specifically the 4.6 experiment.

The code I am using is:

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

# Mesh
domain=Rectangle(Point(0,0),Point(65,120)) - Circle(Point(20,20),5) - Circle(Point(20,100),5) - Circle(Point(36.5,51),10) - Rectangle(Point(0,65),Point(10,65.5))
mesh=generate_mesh(domain, 100)

# 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.28
l = 0.1
lmbda = 1.94e3
mu = 2.45e3


# 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],120) && on_boundary")
bot = CompiledSubDomain("near(x[1],0) && on_boundary")

load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0,0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]
bc_phi = []
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)

# 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
u_r = 0.5
deltaT  = 0.5
tol = 1e-3
conc_f = File ("./Resultados_Ambati/phi.pvd")
fname = open('ForcevsDisp_Ambati.txt', 'w')

# Staggered scheme
while t<=100:
    t += deltaT
    load.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")
                x.append(str(t*u_r))
                y.append(str(round(assemble(fy)*1e-3,1)))
                print("Force:", str(round(assemble(fy)*1e-3,1)),'kN')
                print("Displacement:",str(t*u_r),'mm')
               
fname.close()
print ('Simulation completed') 

As I made a line of code that allows you to view the load and displacement at each iteration, I realized that the results were not matching the paper’s. For example, in my model, at 0.25 mm displacement, the applied load is 0.7 kN, which is clearly wrong, because the max load is 0.6 kN. So I would like to know why the results of my model are not correct.

My suspicion is the applied load, because on the machine where the test was performed the way it is applied is different. Is it possible to apply the load using phenics as it is applied in the testing machine?