Phase Field - The Newton's method does not converge issue

mathematics
1

I wrote a program to solve a problem similar to the Cahn-Hilliard equation,
but I encountered an issue where the results do not converge.

Here is the formula:

image
image665×576 19.7 KB

φ represents phase field.
U represents concentration field.

Here is my code:

import math
import random
from dolfin import *

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.2
        values[1] = 0.0
    def value_shape(self):
        return (2,)
    
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

# Model parameters

dt     = 2.0e-04  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

epsilon_4 = 0.05
fold = 4
alpha = 0   
beta = 0 
M_PI = math.pi
Cinf = 0.25
k = 0.1
Cl = 0.74
a1 = 5.0 * math.sqrt(2.0) / 8.0
a2 = 47.0 / 75.0
epsilon = 100
Vp = 5.0e-6
d0 = 1.3e-8
D = 1.0e-9
mCinf = 2
m = mCinf / Cinf
G = 100.0e2
lT = m*(1 - k)*Cl / G
theta_inf = 1.0

Dtelda = a1 * a2 * epsilon
Vptelda = a1 * a2 * Vp * d0 / D * pow(epsilon, 2.0)
lTtelda = lT / d0 / epsilon
lamda = a1 * epsilon
tau = a1 * a2 * pow(epsilon, 3.0) * pow(d0, 2.0)

T0 = -mCinf / k
Tinf = T0 + G * Vp * theta_inf * lTtelda / Vptelda

# Step in time
t = 0.0
T = 50*dt


parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

du     = TrialFunction(ME)
q1,q2  = TestFunctions(ME)
v1,v2  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dmu,  duu2 = split(du)
dphi, duu1 = split(du)

phi, uu1 = split(u)
mu, uu2  = split(u)

phi0, uu10 = split(u0)
mu0, uu20 = split(u0)


# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)


phi = variable(phi)

dphidx = grad(phi)[0]
dphidy = grad(phi)[1]
thetaa = atan(dphidy/dphidx)

AA      = 1 + epsilon_4 * cos(fold * thetaa)
dAAdphi = diff(AA, phi)
f       = abs(dot(1,phi))**2 * AA * dAAdphi
dfdx    = grad(f)[0]
dfdy    = grad(f)[1]
temp    = (grad(phi)[1] - Vp * t) / lT

uu1_mid = (1.0 - theta) * uu10 + theta * uu1
uu2_mid = (1.0 - theta) * uu20 + theta * uu2

L0 = (1-(1-k)*temp)*(AA**2)*((phi - phi0)/dt)*q1*dx - phi*q1*dx + phi**3*q1*dx + lamda*(1-phi**2)*(mu + temp)*q1*dx - dot(grad(uu1_mid), grad(q1))*dx
L1 = uu1*q2*dx - AA**2*dot(grad(phi), grad(q2))*dx + dot(grad(phi), grad(phi))*AA*dAAdphi*q2*dx
L2 = 0.5*((1+k)-(1-k)*phi)*((mu-mu0)/dt)*v1*dx - (1+(1-k)*mu)*0.5*((phi - phi0)/dt)*v1*dx - dot(grad(uu2_mid), grad(v1))*dx
L3 = uu2*v2*dx - D*(1-phi)*0.5*dot(grad(mu), grad(v2))*dx - ((pow(2,0.5))/4)*(1+(1-k)**mu)*((phi - phi0)/dt)*dot(grad(phi), grad(v2))/sqrt(dot(grad(phi), grad(phi)))*dx
L = L0 + L1 + L2 + L3

a = derivative(L, u, du)

problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("output.pvd", "compressed")


while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    file << (u.split()[0], t)

thank you all