Thanks for your response!
When the Newton solver fails to converge and we reduce the time step (Δt), we should NOT update C_old. Otherwise, we are moving forward with a bad solution, making convergence even harder.
Please correct me if i am wrong.
Here is my code for your kind perusal.
Thanks in advance.
from dolfin import *
from mshr import *
import numpy as np
# Define and refine the mesh
mesh = RectangleMesh(Point(-0.5, -0.5), Point(0.5, 0.5), 200, 200, "crossed")
# Refine mesh around the initial horizontal crack (y=0, x <=0)
class CrackDomain(SubDomain):
def inside(self, x, on_boundary):
return abs(x[1]) < 0.01 and x[0] <= 0.0 # Horizontal crack at y=0
crack_domain = CrackDomain()
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
crack_domain.mark(cell_markers, True)
mesh = refine(mesh, cell_markers)
# Function spaces
V = FunctionSpace(mesh, 'CG', 1) # Phase-field
W = VectorFunctionSpace(mesh, 'CG', 1) # Displacement
Q = FunctionSpace(mesh, 'CG', 1) # Hydrogen concentration
# Define trial and test functions
du, v = TrialFunction(W), TestFunction(W)
dphi, w = TrialFunction(V), TestFunction(V)
dC, q = TrialFunction(Q), TestFunction(Q)
# Functions
u, u_old = Function(W), Function(W)
phi, phi_old = Function(V), Function(V)
C, C_old = Function(Q), Function(Q)
Hold = Function(V)
# Material parameters
E = 210e9 # Young's modulus (Pa)
nu = 0.3 # Poisson's ratio
lmbda = E * nu / ((1 + nu) * (1 - 2*nu)) # Lamé parameter λ
mu = E / (2 * (1 + nu)) # Lamé parameter μ
Gc_0 = 2.7 # Initial fracture energy (N/m)
l = 0.001 # Phase-field length scale (m)
D0 = 1e-9 # Baseline diffusion coefficient (m²/s)
VH = 2e-6 # Partial molar volume (m³/mol)
R = 8.314 # Gas constant (J/(mol·K))
T = 300 # Temperature (K)
chi = 0.89 # Hydrogen embrittlement coefficient
Delta_g0 = 30000 # Gibbs free energy difference (J/mol)
# Boundary conditions
bcbot = DirichletBC(W, Constant((0.0, 0.0)), "near(x[1], -0.5) && on_boundary")
tension = Expression(("0.0", "u_r * t / total_time"), u_r=0.007, t=0.0, total_time=3.0, degree=1)
bctop = DirichletBC(W, tension, "near(x[1], 0.5) && on_boundary")
bc_u = [bcbot, bctop]
bc_C = DirichletBC(Q, Constant(0.1), "on_boundary")
class Crack(SubDomain):
def inside(self, x, on_boundary):
return abs(x[1]) < 1e-3 and x[0] <= 0.0
crack_surface = Crack()
bc_phi = DirichletBC(V, Constant(1.0), crack_surface)
# Initial crack
class InitialCrack(UserExpression):
def eval(self, values, x):
if abs(x[1]) < 1e-3 and x[0] <= 0.0:
values[0] = 1.0
else:
values[0] = 0.0
def value_shape(self):
return ()
phi.interpolate(InitialCrack())
# Constitutive relations
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2 * mu * epsilon(u) + lmbda * tr(epsilon(u)) * Identity(len(u))
def psi(u):
return 0.5 * lmbda * (tr(epsilon(u)))**2 + mu * inner(epsilon(u), epsilon(u))
sigma_H = tr(sigma(u))
theta = C / (C + exp(-Delta_g0 / (R * T)))
Gc = Gc_0 * (1 - chi * theta)
D_eff = (1 - chi * theta) * D0
# Variational forms
# Displacement
E_du = ((1.0 - phi_old**2) * inner(sigma(du), grad(v)) - inner(Constant((0.0, 0.0)), v)) * dx
# Phase-field
def H(u_old, u_new, Hold):
psi_new = psi(u_new)
return conditional(gt(psi_new, Hold), psi_new, Hold)
F_phi = (Gc * l * inner(grad(phi), grad(w)) + (Gc / l) * phi * w - 2 * (1 - phi_old) * H(u_old, u, Hold) * w) * dx
J_phi = derivative(F_phi, phi)
# Hydrogen diffusion
F_C = (( (C - C_old)/Constant(0.01) ) * q * dx
+ D_eff * dot(grad(C), grad(q)) * dx
- (D_eff * VH * C) / (R * T) * dot(grad(sigma_H), grad(q)) * dx )
J_C = derivative(F_C, C)
# Solvers
problem_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), u, bc_u)
solver_disp = LinearVariationalSolver(problem_disp)
problem_phi = NonlinearVariationalProblem(F_phi, phi, bc_phi, J=J_phi)
solver_phi = NonlinearVariationalSolver(problem_phi)
prm_phi = solver_phi.parameters
prm_phi['newton_solver']['linear_solver'] = 'gmres'
prm_phi['newton_solver']['preconditioner'] = 'ilu'
prm_phi['newton_solver']['absolute_tolerance'] = 1e-10
prm_phi['newton_solver']['relative_tolerance'] = 1e-8
problem_C = NonlinearVariationalProblem(F_C, C, bc_C, J=J_C)
solver_C = NonlinearVariationalSolver(problem_C)
prm_C = solver_C.parameters
prm_C['newton_solver']['linear_solver'] = 'gmres'
prm_C['newton_solver']['preconditioner'] = 'ilu'
prm_C['newton_solver']['absolute_tolerance'] = 1e-10
prm_C['newton_solver']['relative_tolerance'] = 1e-8
# Time-stepping parameters
t = 0.0
deltaT = 0.001
num_steps = 3000
# Output
conc_file = XDMFFile("results.xdmf")
conc_file.parameters["flush_output"] = True
# Time loop
max_attempts = 5
successful_step = False
current_deltaT = deltaT
for step in range(num_steps):
t += current_deltaT
tension.t = t
attempts = 0
while attempts < max_attempts:
try:
# Solve displacement using u_old (previous step's solution)
solver_disp.solve()
# Update history variable using u and u_old
Hold.assign(project(H(u_old, u, Hold), V))
# Solve phase-field using phi_old (previous) and u (current)
solver_phi.solve()
# Solve hydrogen diffusion using C_old (previous) and current variables
solver_C.solve()
# Update old variables for next time step
u_old.assign(u)
phi_old.assign(phi)
C_old.assign(C)
successful_step = True
break # Break out of the retry loop if all solves succeed
except RuntimeError as e:
print(f"Solver failed at time {t} with deltaT {current_deltaT}. Retrying with reduced time step.")
current_deltaT *= 0.5 # Reduce time step by half
attempts += 1
if not successful_step:
print(f"Failed to converge after {max_attempts} attempts. Terminating simulation.")
break
# Write output
C.rename("HydrogenConcentration", "label")
phi.rename("PhaseField", "label")
conc_file.write(C, t)
conc_file.write(phi, t)
conc_file.close()
print("Simulation completed.")