Skip to content
Snippets Groups Projects
Commit 242a19fd authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

Can see droplet ripening with 17 droplets that coalesce into one

parent f6be5614
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Fenics_FloryHugg_Binary.ipynb
+ 23
− 21
View file @ 242a19fd
%% Cell type:markdown id: tags:
# Binary mixture based on Flory Huggins free energy
For $\phi_{tot}$ in a binary mixture, assuming a Flory Huggins free energy density, we obtain (Overleaf):
Reference for FEM: https://www.comsol.com/multiphysics/finite-element-method
$\partial_\mathrm{t}\phi_\mathrm{tot}=k_\mathrm{B}T\Gamma_\mathrm{0} \,\nabla\cdot\left[\left[1-2\chi\phi_\mathrm{tot}(1-\phi_\mathrm{tot})\right]\nabla \phi_\mathrm{tot}-\phi_\mathrm{tot}(1-\phi_\mathrm{tot})\kappa\nabla\nabla^2\phi_\mathrm{tot}\right]$
To solve this equation we follow the [Fenics Cahn-Hilliard example](https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1) (the maths is explained [here](https://fenicsproject.org/docs/dolfinx/dev/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html), for our free energy see overleaf).
$\partial_t \phi_{\mathrm{tot}} = \nabla \cdot \Gamma_{\mathrm{tot}}\nabla\mu_\mathrm{tot}$
$\mu_\mathrm{tot}=\nu\frac{\partial f}{\partial \phi_\mathrm{tot}}=\,k_\mathrm{B}T\left[\ln\phi_\mathrm{tot} - \ln\left(1-\phi_\mathrm{tot}\right)+\chi(1-2\phi_\mathrm{tot})-\kappa\nabla^2 \phi_\mathrm{tot} +1-n \right]$
The **boundary conditions** for the flux of the chemical potential and concentration read:
$\nabla\mu_\mathrm{tot}\cdot n=0 $ on $\partial \Omega\;,$
$\nabla\phi_\mathrm{tot}\cdot n = 0$ on $\partial \Omega$
Cast into the **weak form** and integrating by parts we obtain:
$\int_{\Omega} \partial_t c \;q \; dx = \int_{\partial \Omega} \Gamma_\mathrm{tot}\nabla\mu_\mathrm{tot}q \; ds - \int_\Omega \Gamma_\mathrm{tot} \nabla \mu_\mathrm{tot} \nabla q\;dx = - \int_\Omega \Gamma_\mathrm{tot} \nabla \mu_\mathrm{tot} \nabla q\;dx\;, $
$\int_\Omega \mu_\mathrm{tot} q \; dx = k_\mathrm{B} T \int_\Omega\left(\ln\left(\frac{\phi_\mathrm{tot}}{1-\phi_\mathrm{tot}}\right)+\chi(1-2\phi_\mathrm{tot})\right)q + \kappa\nabla\phi_\mathrm{tot}\nabla q \; dx$
where we dropped the boundary terms (not explicitly written in second equation) due to the boundary conditions/weak form. We also disregard 1-n, because constants don't change the flux in the absence of chemical reactions (an offset doesn't change the solution in this case, try out by adding $-6.0*v*dx$ to L1, for this check out fenics tutorial Poisson equation).
According to the review (Weber et al. 2019, eqn. 2.29) the boundary length scale is
$w = \sqrt{2\kappa/b} = \left( \frac{\kappa}{k_\mathrm{B}T\chi}\right) ^{3/5} \sqrt{ \frac{\chi}{\chi-2} }$
%% Cell type:code id: tags:
``` python
# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins
# as discussed on 28/11/2019
# Example can be found at
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/
# documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1
# Runs with fenics 2019.01
# The resulting .pvd file can be opened using default settings in ParaView
import matplotlib.pyplot as plt
import mshr as ms
import numpy as np
from scipy.optimize import curve_fit
import time
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):
# if (x[0] > -0.1) & (x[0] < 0.1) & (x[1] > -0.1) & (x[1] < 0.1) & (x[2] > -0.1) & (x[2] < 0.1):
# # values[0] = 0.63 + 0.02*(0.5 - random.random())
# values[0] = 0.82
# if (x[0]**2+x[1]**2+x[2]**2 < 0.25):
# values[0] = 0.73
if (self.in_rad(x, [0.5, 0.5, 0.5], 0.05)):
# self.in_rad(x, [0.5, 0.5, 0.5], 0.1) or
# self.in_rad(x, [0.8, 0.8, 0.8], 0.1) or
# self.in_rad(x, [0.2, 0.8, 0.8], 0.1) or
# self.in_rad(x, [0.2, 0.8, 0.2], 0.1) or
# self.in_rad(x, [0.2, 0.2, 0.8], 0.1) or
# self.in_rad(x, [0.8, 0.2, 0.2], 0.1) or
# self.in_rad(x, [0.8, 0.8, 0.2], 0.1) or
# self.in_rad(x, [0.8, 0.8, 0.2], 0.1) or
# self.in_rad(x, [0.8, 0.2, 0.8], 0.1)):
values[0] = 0.75
# if (self.in_rad(x, [0.5, 0.5, 0.5], 0.05)):
# values[0] = 0.75
if (self.in_rad(x, [0.5, 0.5, 0.5], 0.05) or
self.in_rad(x, [0.2, 0.2, 0.2], 0.05) or
self.in_rad(x, [0.8, 0.8, 0.8], 0.05) or
self.in_rad(x, [0.2, 0.8, 0.8], 0.05) or
self.in_rad(x, [0.2, 0.8, 0.2], 0.05) or
self.in_rad(x, [0.2, 0.2, 0.8], 0.05) or
self.in_rad(x, [0.8, 0.2, 0.2], 0.05) or
self.in_rad(x, [0.8, 0.8, 0.2], 0.05) or
self.in_rad(x, [0.8, 0.2, 0.8], 0.05) or
self.in_rad(x, [0.65, 0.65, 0.65], 0.05) or
self.in_rad(x, [0.35, 0.35, 0.35], 0.05) or
self.in_rad(x, [0.65, 0.35, 0.65], 0.05) or
self.in_rad(x, [0.65, 0.65, 0.35], 0.05) or
self.in_rad(x, [0.35, 0.65, 0.65], 0.05) or
self.in_rad(x, [0.65, 0.35, 0.35], 0.05) or
self.in_rad(x, [0.35, 0.35, 0.65], 0.05) or
self.in_rad(x, [0.35, 0.65, 0.35], 0.05)):
values[0] = 0.76
else:
values[0] = 0.268
values[1] = 0.0
# values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
def in_rad(self, x, drop, rad):
if np.sqrt((x[0]-drop[0])**2+(x[1]-drop[1])**2+(x[2]-drop[2])**2)<rad:
return True
else:
return False
# 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
kappa = 0.25*10e-05 # surface parameter
dt = 1.0e-02 # time step
dt = 0.5e-02 # time step
X = 2.2 # Chi Flory Huggins parameter
# time stepping family, e.g.:
# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
theta = 0.5
# Form compiler optionsmai
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Create mesh and build function space
x_mesh =100
y_mesh =100
z_mesh =100
# mesh = UnitSquareMesh.create(x_mesh, y_mesh, CellType.Type.quadrilateral)
# mesh = RectangleMesh(Point(0, 0), Point(1, 0.02), x_mesh, y_mesh)
# mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), x_mesh, y_mesh, z_mesh)
mesh = Mesh('box.xml')
mesh = Mesh('Meshes/box_spheres.xml')
# domain = ms.Sphere(Point(0, 0, 0), 1.0)
# mesh = ms.generate_mesh(domain, 100)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# 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)
dfdc = ln(c/(1-c))+(1-2*X*c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*dot(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
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
# solver.parameters["linear_solver"] = "lu"
solver.parameters["linear_solver"] = 'gmres'
solver.parameters["preconditioner"] = 'ilu'
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
```
%% Cell type:code id: tags:
``` python
# Output file
# file = File("output.pvd", "compressed")
file_c = XDMFFile('c_box.xdmf')
file_c = XDMFFile('c_long.xdmf')
# Step in time
t = 0.0
T = 100*dt
ti = time.time()
while (t < T):
# file << (u.split()[0], t)
file_c.write(u.split()[0], t)
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
print(time.time() - ti)
file_c.close()
```
%% Cell type:code id: tags:
``` python
def tanh_fit(x, a, b, c, d):
return np.tanh((x-a)/b)*c+d
# xdata = np.linspace(0, x_mesh, x_mesh+1)
# ydata = u.compute_vertex_values()[0:x_mesh+1]
xdata = np.linspace(0.1, 0.5, 100)
ydata = np.array([u0([x, 0.5, 0.5])[0] for x in xdata])
popt, pcov = curve_fit(tanh_fit, xdata, ydata)
plt.plot(xdata[65:], ydata[65:])
a, b, c, d = popt
plt.plot(xdata[65:], tanh_fit(xdata[65:], a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5
```
%% Cell type:code id: tags:
``` python
mesh.num_cells()
```
%% Cell type:markdown id: tags:
## Thoughts on stability of algorithm
- relatively robust to grid size changes
- initial conditions seem to have strong influence
- time step also plays a role for Newton solver
- best results with small time step and initial conditions close to final concentrations
%% Cell type:code id: tags:
``` python
```
......
Meshes/box_with_spheres.geo 0 → 100644
+ 44
− 0
View file @ 242a19fd
lc = 0.1;
Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc} ;
Point(3) = {1, 1, 0, lc} ;
Point(4) = {0, 1, 0, lc} ;
Point(5) = {.5, .5, .5, lc} ;
Point(6) = {.2, .2, .2, lc} ;
Point(7) = {.8, .8, .8, lc} ;
Point(8) = {.2, .8, .8, lc} ;
Point(10) = {.8, .2, .8, lc} ;
Point(11) = {.8, .8, .2, lc} ;
Point(12) = {.8, .2, .2, lc} ;
Point(13) = {.2, .8, .2, lc} ;
Point(14) = {.2, .2, .8, lc} ;
Point(15) = {.35, .35, .35, lc} ;
Point(16) = {.65, .65, .65, lc} ;
Point(17) = {.35, .65, .65, lc} ;
Point(18) = {.65, .35, .65, lc} ;
Point(19) = {.65, .65, .35, lc} ;
Point(20) = {.65, .35, .35, lc} ;
Point(21) = {.35, .65, .35, lc} ;
Point(22) = {.35, .35, .65, lc} ;
Line(1) = {1,2} ;
Line(2) = {3,2} ;
Line(3) = {3,4} ;
Line(4) = {4,1} ;
Curve Loop(1) = {4,1,-2,3} ;
Plane Surface(1) = {1} ;
Physical Curve(5) = {1, 2, 4} ;
Physical Surface("My surface") = {1} ;
Extrude {0, 0, 1} { Surface{1}; }
Field[1] = Distance;
Field[1].NodesList = {5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22};
Field[1].NNodesByEdge = 100;
Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc / 20;
Field[2].LcMax = lc;
Field[2].DistMin = 0.05;
Field[2].DistMax = 0.1;
Field[3] = Min;
Field[3].FieldsList = {2};
Background Field = 3;
Physical Volume("The volume", 1) = {1};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment