Crack Not Propagating in Phase Field Simulation (Hydrogen Embrittlement) - FEniCS

General
,
1

I am trying to simulate crack propagation in a single-edge notched plate subjected to tension in a hydrogen environment using the phase-field approach in FEniCS. However, despite applying incremental displacement, the crack is not propagating in my simulation.
Observed Issues:

  • The phase field (ϕ) remains unchanged throughout the simulation.
  • Load-displacement curve does not show failure.
  • History variable H updates, but crack evolution is not happening.

Possible Causes:

  • Incorrect phase-field equation formulation?
  • Issue with hydrogen-dependent fracture energy?
  • Need finer mesh resolution or fracture length scale l?

Would appreciate any suggestions! Should I modify something or debug in a specific way?

Please check my code, Thanks in advance!

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# Geometry and mesh parameters
length, height = 1.0, 1.0

mesh = RectangleMesh(Point(0, 0), Point(length, height), 100, 100)  # Instead of 50, 50


# Marking initial crack (Horizontal crack at y=0.5, from x=0 to x=0.5)
class InitialCrack(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5) and (0.0 <= x[0] <= 0.5)

crack = InitialCrack()

# Material properties
E, nu = 210e3, 0.3
Gc0, l = 2, 0.015
D = 0.0127
kappa = 1e-8
chi = 0.89  # Hydrogen damage coefficient

# Elasticity parameters
mu = E / (2*(1 + nu))
lmbda = E*nu / ((1 + nu)*(1 - 2*nu))

# Mesh and function spaces
V_u = VectorFunctionSpace(mesh, 'CG', 2)
V_phi = FunctionSpace(mesh, 'CG', 1)
V_c = FunctionSpace(mesh, 'CG', 1)

# Functions
u, phi, C = Function(V_u), Function(V_phi), Function(V_c)
v, dphi, dc = TestFunction(V_u), TestFunction(V_phi), TestFunction(V_c)

# Boundary conditions
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], height)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class CrackBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5) and (0.0 <= x[0] <= 0.5)

top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()
crack_boundary = CrackBoundary()

bc_u_bottom = DirichletBC(V_u, Constant((0.0, 0.0)), bottom_boundary)
bc_phi = DirichletBC(V_phi, Constant(1.0), crack_boundary)
bc_C = DirichletBC(V_c, Constant(1.0), 'on_boundary')

# Initial conditions
C0 = 5.0e-7  # 0.5 ppm in FEniCS units
phi_old = interpolate(Constant(1.0), V_phi)
C_old = interpolate(Constant(C0), V_c)
H = Function(V_phi)  # Initialize H with the same function space as phi

# Define hydrogen-dependent fracture energy
Gc = project(Gc0 * (1 - chi * (C_old / (C_old + 0.03))), V_phi)

# Define variational forms
def eps(u):
    return sym(grad(u))

def sigma(u, phi):
    return ((1 - phi)**2 + kappa) * (lmbda * tr(eps(u)) * Identity(2) + 2 * mu * eps(u))

# Time-stepping parameters
num_steps = 1000
dt = 1000  # Time step in seconds
max_disp = 0.05  # Maximum displacement to track full crack propagation
t = 0
disp_increment = 5e-6  # Correct per the paper

# Output files
vtkfile_phi = File('phi.pvd')
vtkfile_C = File('C.pvd')
vtkfile_u = File('u.pvd')

# Tracking Load vs Displacement
load_displacement = []

dx = Measure("dx", domain=mesh)  # Ensure dx is correctly defined

for step in range(num_steps):
    t += dt
    disp = step * disp_increment  # Incremental displacement
    bc_u_top = DirichletBC(V_u.sub(1), Constant(disp), top_boundary)
    bcs_u = [bc_u_bottom, bc_u_top]
    
    # Step 1: Solve elasticity problem
    elasticity_eq = inner(sigma(u, phi_old), eps(v))*dx
    solve(elasticity_eq == 0, u, bcs_u,
          J=derivative(elasticity_eq, u),
          solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
    
    # Step 2: Compute history variable
    psi = 0.5 * inner(sigma(u, phi_old), eps(u))
    H_new = project(psi, V_phi)
    H.assign(project(conditional(ge(H_new, H), H_new, H), V_phi))  # Ensure crack irreversibility
    #H.assign(max(H, H_new))  # Ensure crack irreversibility
    print(f"Step: {step}, Max H: {H.vector().max()}")

    # Step 3: Update Gc with latest C values
    theta = C_old / (C_old + 0.03)  # Hydrogen concentration factor
    Gc = Gc0 * (1 - chi * theta)
    
    # Step 4: Solve phase-field equation
    phi_eq = ((Gc/l)*phi*dphi + Gc*l*dot(grad(phi), grad(dphi))
              - 2*(1 - phi)*H*dphi)*dx
    solve(phi_eq == 0, phi, bc_phi,
          J=derivative(phi_eq, phi),
          solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
    
    # Step 5: Compute hydrostatic stress AFTER phase-field update
    sigma_H = project(tr(sigma(u, phi)) / 3, V_c)
    
    # Step 6: Solve hydrogen diffusion equation
    C_eq = ((C - C_old)/dt)*dc*dx + dot(D*grad(C), grad(dc))*dx \
        - (D*C/(8.314*300))*2e-6*dot(grad(sigma_H), grad(dc))*dx
    solve(C_eq == 0, C, bc_C,
          J=derivative(C_eq, C),
          solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
    
    # Step 7: Update solutions
    phi_old.assign(phi)
    C_old.assign(C)
    
    # Step 8: Track Load vs Displacement
    load = assemble(inner(sigma(u, phi), eps(u)) * dx)
    load_displacement.append((disp, load))
    
    # Step 9: Save results
    vtkfile_phi << (phi, t)
    vtkfile_C << (C, t)
    vtkfile_u << (u, t)
    
    print(f"Step: {step}, Disp: {disp}, Load: {load}")

# Plot Load vs Displacement
plt.figure()
disp_vals, load_vals = zip(*load_displacement)
plt.plot(disp_vals, load_vals, '-o')
plt.xlabel("Displacement (mm)")
plt.ylabel("Load (N)")
plt.title("Load vs Displacement Curve")
plt.grid(True)
plt.show()

print("Simulation completed successfully.")