Hello,
I have extended/modified the fenics19 Cahn-Hilliard demo code to consider
the contact angle boundary condition
\nabla c \cdot n = 1/(eps sqrt(2)) cos( \theta ) (1 - c^2) (contact angle)
eps := parameter related to the interface thickness
\theta := contact angle
As test case we use a unit square that is “filled” in the upper and lower
parts with c = +1, and c = -1, respectively.
This initial condition defines a fluid-fluid interface that is orthogonal to
the domain boundary.
Note that the right hand side of the contact angle b.c. vanishes
for a contact angle of 90 degrees.
In fact the fluid-fluid interface remains orthogonal to the domain boundary
if we do nothing or impose a contact angle b.c. with \theta = 90.
However, specifying some other angle, say 45 degrees, lets the code break down.
It would be very interesting to find out why this happens and how one can fix it.
Any hints are very appreciated.
Kind regards,
Babak S. Hosseini
The code is attached to this message.
# Mixed form (phase field, chemical potential) Cahn-Hilliard solver.
#
# The code is based on the fenics19 Cahn-Hilliard demo and has slightly
# been modified to additionally consider a nonlinear Neumann boundary
# condition for imposing the contact angle boundary condition:
#
# \nabla c \cdot n = 1/(eps sqrt(2)) cos( \theta ) (1 - c^2) (contact angle)
# eps := parameter related to the interface thickness
# \theta := contact angle
#
# As far as modifications of the original code are concerned, the
# double well function has been replaced by
# f = 0.25 * (c-1)**2 * (c+1)**2
# which is suitable for c \in [-1,1]
#
# As test case we use a unit square that is "filled" in the upper and lower
# parts with c = +1, and c = -1, respectively.
# This initial condition defines a fluid-fluid interface that is orthogonal to
# the domain boundary.
# Note that the right hand side of the contact angle b.c. vanishes
# for a contact angle of 90 degrees.
# In fact the fluid-fluid interface remains orthogonal to the domain boundary
# if we do nothing or impose a contact angle b.c. with \theta = 90.
# However, specifying some other angle, say 45 degrees, lets the code break down.
# It would be very interesting to find out why this happens and how one can fix it.
#
#
# Equations referred to are from the following paper:
# Isogeometric Analysis of the Navier-Stokes-Cahn-Hilliard
# equations with application to incompressible two-phase flows.
# Babak S. Hosseini et al.
#
# Run the code as follows:
# $: python3 demo_cahn-hilliard.py
# =================================================================
import random
from dolfin import *
# Class representing the intial conditions
class InitialConditions(UserExpression):
def __init__(self, **kwargs):
random.seed(2 + MPI.rank(MPI.comm_world))
super().__init__(**kwargs)
def eval(self, values, x):
# test case: divided domain (unit square with the upper half filled with +1 and the lower half filled with -1):
minval = -1
maxval = 1
if x[1] <= 0.5:
values[0] = minval
else:
values[0] = maxval
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
iface_eps = 0.04 # inteface thickness (lmbda = eps**2)
st = 1 # surface tension
# Paramters alpha and beta related to the surface tension and diffuse interface width.
# See eqn (24)
aux = 3.0 / (2 * sqrt( 2.0 ));
alpha = aux * st * iface_eps;
beta = aux * st / iface_eps;
dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Create mesh and build function space
mesh = UnitSquareMesh.create(200, 200, CellType.Type.quadrilateral)
#bsh begin
# FIXME: We need P2 Lagrange elements
# The code fails to work with elements of degree 2..very strange
#bsh end
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# Trial and test functions of the space ``ME`` are now defined::
# 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(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
# Compute the chemical potential df/dc
c = variable(c)
# Double well function psi(c) penalizing the mixing of phases.
# Suitable for c \in [-1,1]
# See eqn (12)
f = 0.25 * (c-1)**2 * (c+1)**2
dfdc = diff(f, c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
# Scaling parameter for the mobility function
D = 1
# Mobility function according to eqn (17)
m = D * ((c0)**2 - 1)**2
sqrt_two = 1.41421356237310
pi = 3.14159265358979
pi_over_two = 1.57079632679490
pi_over_four = 0.78539816339744
# Contact angle between the fluid-fluid interface and the domain boundary
angleDeg = 90.0
# Weak statement of the equations:
# ---------------------------------
# Bilinear form for the evolution of the phase field (c). Test function: q
L0 = c * q * dx - c0 * q * dx + dt * m * dot( grad(mu_mid), grad(q) ) * dx
# Bilinear form for the evolution of the chemical potential (mu). Test function: v
L1 = mu * v * dx - beta * dfdc * v * dx - alpha * dot( grad(c), grad(v) ) * dx
# Form for the imposition of the nonlinear Robin b.c. with test function q:
L2 = ( 1.0 / (sqrt_two * iface_eps) * cos( angleDeg * pi / 180.0 ) * (1.0 - c*c) ) * q * ds
#L = L0 + L1 + L2
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
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")
file << (u.split()[0], 0.0)
# Step in time
maxSteps = 10
t = 0.0
T = maxSteps * dt
step = 1
while (t < T):
print( "" )
print( "------- time iter: {} / {} -----".format(step,maxSteps) )
t = t + dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
step = step + 1