Graph of Cahn-Hilliard equation

Question

Hi! there,

I have set up a new initial condition (triangle) in the demonstrative code of Cahn-Hilliard equation. After executing the code, the triangle become circle gradually as time increased.
I don't quite understand the reason why it change to circle.
Does anyone know the reason in Physics or other perspectives?

Thank you!

class InitialConditions(Expression):
def eval(self,values,x):
    if between(x[0],(0.3,0.8)) and between(x[1],(0.3,0.8*x[0]+0.06)):
         values[0] = 1.0
    else:
         values[0] = 0.0
def value_shape(self):
     return (2,)

######################################################################

Class for interfacing with the Newton solver

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

#

lmbda = 1.0e-03 # surface parameter (IValue= 1.0e-02, AT=1.0e-03)
dt = 1.0e-03 # time step (IValue= 5.0e-06, AT=1.0e-03)
theta = 0.7 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson, AT=0.7

#

Form compiler options

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

Create mesh and define function spaces

#

mesh = UnitSquareMesh(64, 64) # IV=96 AT=64
V = FunctionSpace(mesh, "Lagrange", 1)
ME = V*V

#

Define trial and test functions

du = TrialFunction(ME)
q, v = TestFunctions(ME)

Define functions

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

Split mixed functions

dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)

Create intial conditions and interpolate

u_init = InitialConditions()
u.interpolate(u_init)
u0.interpolate(u_init)

Compute the chemical potential df/dc

The f here is different from the dissertation(ie.f = F)

#

c = variable(c)
f = c2*(1-c)2 # IV=100*c2*(1-c)2, AT=c2*(1-c)2
dfdc = diff(f, c)

#

mu_(n+theta)

mu_mid = (1.0-theta)mu0 + thetamu

Weak statement of the equations (Add:(1.0/lmbda)* )

#

L0 = cqdx - c0qdx + dtdot(grad(mu_mid), grad(q))dx
L1 = muvdx - dfdcvdx - lmbdadot(grad(c), grad(v))dx
L = L0 + L1

#

Compute directional derivative about u in the direction of du (Jacobian)

a = derivative(L, u, du)

Create nonlinear problem and Newton solver

IValue: relative_tolerance = 1e-6

#

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

#

Output file

file = File("MYoutput_CH.pvd", "compressed")

Step in time

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

plot(u.split()[0])
interactive()

Answer 1

Otto, the reason the triangle evolves to a circle is due to the fact that the Cahn-Hilliard equation involves a free energy functional F which is the integral of the sum of the "Density gradient free energy" (E_G) and "Homogeneous free energy" (E_H).
So we have: F = \int_{Omega} E_G + E_H d\Omega.

E_G = 0.5 eps^2 abs( grad(c) )^2 (Density gradient free energy)
E_H = f(c) (Homogeneous free energy)
f(c) is a double well function in [-1,1] or [0,1]
M(c) is a mobility function

Let the chemical potential (aka the total free energy variation) \mu be the
derivative of F w.r.t. c:
\mu = df/dc = ... = f'(c) - eps^2 laplace(c)

Note that this happens to appear in the Cahn-Hilliard eqn (CH)
dc/dt = div [ M(c) grad ( - eps^2 laplace(c) + f'(c) ) ]

CH minimizes both the Homogeneous free energy and Density gradient free energy.
The latter is associated with the energy of the interface separating the phases
and as a natural consequence one ends up with a circular interface in 2D.

Kind regards,
Babak S. Hosseini