From f3ba88ea0a658412e0126fd6276acbfec682ad3b Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Thu, 28 Nov 2019 17:11:41 +0100
Subject: [PATCH] Amended Cahn Hilliard code in a minimal way to account for
Flory Huggins binary mixture with the corresponding mobility scaling. Seems
to work, but not sure about the volume fraction to concentration conversion.
---
Fenics_Cahn_Hilliard.ipynb | 31 ++++++++++++++++---------------
1 file changed, 16 insertions(+), 15 deletions(-)
diff --git a/Fenics_Cahn_Hilliard.ipynb b/Fenics_Cahn_Hilliard.ipynb
index b38153e..cb709ff 100644
--- a/Fenics_Cahn_Hilliard.ipynb
+++ b/Fenics_Cahn_Hilliard.ipynb
@@ -18,7 +18,13 @@
" random.seed(2 + MPI.rank(MPI.comm_world))\n",
" super().__init__(**kwargs)\n",
" def eval(self, values, x):\n",
- " values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
+ " if x[0] > 0.5:\n",
+ "# values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
+ " values[0] = 0.6\n",
+ " else:\n",
+ " values[0] = 0.35 \n",
+ " values[1] = 0.0\n",
+ "# values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
" values[1] = 0.0\n",
" def value_shape(self):\n",
" return (2,)\n",
@@ -33,14 +39,14 @@
" def J(self, A, x):\n",
" assemble(self.a, tensor=A)\n",
"# Model parameters\n",
- "lmbda = 1.0e-02 # surface parameter\n",
- "dt = 5.0e-06 # time step\n",
+ "lmbda = 2.0e-02 # surface parameter\n",
+ "dt = 5.0e-03 # time step\n",
"theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson\n",
"# Form compiler options\n",
"parameters[\"form_compiler\"][\"optimize\"] = True\n",
"parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
"# Create mesh and build function space\n",
- "mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)\n",
+ "mesh = UnitSquareMesh.create(96, 3, CellType.Type.quadrilateral)\n",
"P1 = FiniteElement(\"Lagrange\", mesh.ufl_cell(), 1)\n",
"ME = FunctionSpace(mesh, P1*P1)\n",
"# Define trial and test functions\n",
@@ -60,12 +66,14 @@
"u0.interpolate(u_init)\n",
"# Compute the chemical potential df/dc\n",
"c = variable(c)\n",
- "f = 100*c**2*(1-c)**2\n",
- "dfdc = diff(f, c)\n",
+ "# f = 100*c**2*(1-c)**2\n",
+ "# dfdc = diff(f, c)\n",
+ "X = 2.2\n",
+ "dfdc = ln(c/(1-c))+(1-2*X*c)\n",
"# mu_(n+theta)\n",
"mu_mid = (1.0-theta)*mu0 + theta*mu\n",
"# Weak statement of the equations\n",
- "L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx\n",
+ "L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*dot(grad(mu_mid), grad(q))*dx\n",
"L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx\n",
"L = L0 + L1\n",
"# Compute directional derivative about u in the direction of du (Jacobian)\n",
@@ -81,20 +89,13 @@
"\n",
"# Step in time\n",
"t = 0.0\n",
- "T = 50*dt\n",
+ "T = 5000*dt\n",
"while (t < T):\n",
" t += dt\n",
" u0.vector()[:] = u.vector()\n",
" solver.solve(problem, u.vector())\n",
" file << (u.split()[0], t)"
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
--
GitLab