notebooks/Phase field theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c33fb5ec",
   "metadata": {},
   "source": [
    "# Derivation of Governing Equations in Phase Field Method: Allen-Cahn and Cahn-Hilliard Equations\n",
    "\n",
    "Zhengtao Gan, August 9, 2023\n",
    "\n",
    "## Recall variational derivative\n",
    "\n",
    "Many problems in physics and mechanics involve finding extrema, either local minima or maxima. This could relate to identifying the maximum or minimum values of a scalar function, where the techniques of differential calculus, such as taking derivatives, are employed. Alternatively, it may involve finding a function that minimizes or maximizes a certain path or integral. In this latter case, the calculus of variations and the use of variational derivatives are essential tools. These concepts form the basis for understanding a wide range of physical phenomena and are central to fields such as classical mechanics and phase field modeling\n",
    "\n",
    "**Reference**: [Variational Calculus Notes by P. Gutierrez](http://www.nhn.ou.edu/~gut/notes/cm/lect_05.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e008a06",
   "metadata": {},
   "source": [
    "### The stationary value of integral determined by variational derivative\n",
    "\n",
    "#### Variational Calculus\n",
    "\n",
    "Many problems involve determining the extrema of integrals. The general problem can be stated as:\n",
    "\n",
    "Given a definite integral\n",
    "\n",
    "$$ I = \\int_{a}^{b} F(y, y', x) \\, dx $$\n",
    "\n",
    "with boundary conditions\n",
    "\n",
    "$$ f(a) = \\alpha \\quad \\text{and} \\quad f(b) = \\beta $$\n",
    "\n",
    "we aim to find the function $ y = f(x) $ that provides a stationary value to the integral.\n",
    "\n",
    "Assuming that $ y = f(x) $ indeed provides a stationary value, we can evaluate the integral for a slightly modified value of the function:\n",
    "\n",
    "$$ \\bar{f}(x) = f(x) + \\epsilon \\phi(x) $$\n",
    "\n",
    "where $ \\phi(x) $ is continuous, differentiable, and $ \\epsilon $ is an arbitrarily small parameter. Comparing the two functions at a point $ x $, we get:\n",
    "\n",
    "$$ \\delta y = \\bar{f}(x) - f(x) = \\epsilon \\phi(x) $$\n",
    "\n",
    "Note the difference between $ \\delta y $ and $ dy $. The differential $ dy $  corresponds an infinitesimal change in $ f(x) $ caused by an infinitesimal change $dx$, while the variation $ \\delta y $ is an infinitesimal change in $y$ producing a new function.\n",
    "\n",
    "Given the boundary conditions remain fixed:\n",
    "\n",
    "$$ \\delta f(x) |_{x=a} = 0 \\quad \\text{and} \\quad \\delta f(x) |_{x=b} = 0 $$\n",
    "\n",
    "Varying the function $ F(y, y', x) $ with a virtual displacement in $ y $, we obtain:\n",
    "\n",
    "$$ \\delta F(y, y', x) = F(y + \\epsilon \\phi, y' + \\epsilon \\phi', x) - F(y, y', x) = \\epsilon \\left( \\frac{\\partial F}{\\partial y} \\phi + \\frac{\\partial F}{\\partial y'} \\phi' \\right) $$\n",
    "\n",
    "where the higher order terms in the Taylor series are neglected since they approach zero. Subsequently, the variation on the integral becomes:\n",
    "\n",
    "$$ \\delta \\int_{a}^{b} F \\, dx =  \\int_{a}^{b} \\delta F \\, dx = \\epsilon \\int_{a}^{b} \\left( \\frac{\\partial F}{\\partial y} \\phi + \\frac{\\partial F}{\\partial y'} \\phi' \\right) \\, dx $$\n",
    "\n",
    "We calculate the rate of change of the integral by dividing by $ \\epsilon $:\n",
    "\n",
    "$$\n",
    "\\frac{\\delta I}{\\epsilon} = \\int_{a}^{b} \\left( \\frac{\\partial F}{\\partial y} \\phi + \\frac{\\partial F}{\\partial y'} \\phi' \\right) dx \\quad \n",
    "$$\n",
    "\n",
    "\n",
    "To further reduce this integral, we carry out an integration by parts on the second term:\n",
    "\n",
    "$$\n",
    "\\int_{a}^{b} \\frac{\\partial F}{\\partial y'} \\phi' dx = \\left[ \\frac{\\partial F}{\\partial y'} \\phi \\right]_a^b - \\int_{a}^{b} \\frac{d}{dx} \\left( \\frac{\\partial F}{\\partial y'} \\right) \\phi dx \\quad \n",
    "$$\n",
    "\n",
    "The first term drops out since $ \\phi = 0 $ at the endpoint (no variation in the boundary conditions). The rate of change of the integral becomes:\n",
    "\n",
    "$$\n",
    "\\frac{\\delta I}{\\epsilon} = \\int_{a}^{b} \\left( \\frac{\\partial F}{\\partial y} - \\frac{d}{dx} \\left( \\frac{\\partial F}{\\partial y'} \\right) \\right) \\phi dx = 0 \\quad \n",
    "$$\n",
    "\n",
    "Since the function $ \\phi $ is arbitrary, it follows that the rate of change of the integral is zero only if the term in brackets is zero.\n",
    "\n",
    "We deduce that for the integral to be stationary:\n",
    "\n",
    "$$ \n",
    "\\frac{\\partial F}{\\partial y} - \\frac{d}{dx} \\left( \\frac{\\partial F}{\\partial y'} \\right) = 0 \n",
    "$$\n",
    "\n",
    "In summary, the variational derivative of a functional can be expressed as\n",
    "\n",
    "**One-dimensional**:\n",
    "\n",
    "$$\n",
    "F = \\int f(x, \\eta(x), \\eta'(x)) \\, dx\n",
    "$$\n",
    "\n",
    "The variational derivative is given by:\n",
    "\n",
    "$$\n",
    "\\frac{\\delta F}{\\delta \\eta} = \\frac{\\partial f}{\\partial \\eta} - \\frac{d}{dx} \\left( \\frac{\\partial f}{\\partial \\eta'} \\right)\n",
    "$$\n",
    "\n",
    "**Three-dimensional**:\n",
    "\n",
    "$$\n",
    "F = \\int_V f(\\mathbf{r}, \\eta(\\mathbf{r}), \\nabla \\eta(\\mathbf{r})) \\, dV \\quad \n",
    "$$\n",
    "\n",
    "The variational derivative in three dimensions is:\n",
    "\n",
    "$$\n",
    "\\frac{\\delta F}{\\delta \\eta} = \\frac{\\partial f}{\\partial \\eta} - \\frac{d}{dx} \\left( \\frac{\\partial f}{\\partial \\eta_x} \\right) - \\frac{d}{dy} \\left( \\frac{\\partial f}{\\partial \\eta_y} \\right) - \\frac{d}{dz} \\left( \\frac{\\partial f}{\\partial \\eta_z} \\right)\n",
    "$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "45360204",
   "metadata": {},
   "source": [
    "# Governing equation of phase field method\n",
    "\n",
    "Using the concept of variational derivative, we can derive the governing equations used in the phase field method: Allen-Cahn equation and Cahn-Hillard equation. The phase field model is typically governed by the time-dependent **Ginzburg-Landau equation**:\n",
    "The rate of change of the non-conserved $i$-th order parameter $\\eta_i$ is assumed to be linearly proportional to its driving force, or the variational derivative of the total free energy $F$ with respect to the phase-field variable (i.e., the order parameter):\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\eta_i}{\\partial t} = -M_i \\frac{\\delta F}{\\delta \\eta_i}\n",
    "$$\n",
    "\n",
    "Where:\n",
    "- $\\eta_i$: the $i$-th order parameter\n",
    "- $M_i$: the $i$-th grain boundary mobility\n",
    "- $F$: total free energy\n",
    "\n",
    "This equation is derived based on **linear kinetics assumption** and is only valid for systems that are not too far from equilibrium (*Chen, L.Q. and Zhao, Y., 2022. From classical thermodynamics to phase-field method. Progress in Materials Science, 124, p.100868*). \n",
    "\n",
    "The expression for total free energy $F$ is given by:\n",
    "\n",
    "$$\n",
    "F = \\int_V \\left[ f(\\boldsymbol{\\eta}) + \\frac{\\kappa}{2} \\sum_{i=1}^N \\|\\nabla \\eta_i\\|^2 \\right] dV\n",
    "$$\n",
    "\n",
    "With the first term is the local free energy density and the second term is the gradient energy density and:\n",
    "- $V$: volumetric domain\n",
    "- $\\kappa$: gradient energy coefficient\n",
    "- $N$: the number of order parameters\n",
    "\n",
    "By deriving the variational derivative as discussed above, we obtain the **Allen-Cahn equation**:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\eta_i}{\\partial t} = -M_i \\left( \\frac{\\partial f}{\\partial \\eta_i} - \\kappa \\sum_{i=1}^{N} \\nabla^2 \\eta_i \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bef2cf44",
   "metadata": {},
   "source": [
    "#### How is the formulation of free energy obtained?  ####\n",
    "\n",
    "The free energy density is assumed to be a function of an order parameter and its spatial derivatives. While we primarily focus on a single order parameter, the derivation can be extended to multiple order parameters with ease. Broadly, the total free energy, $F$, is given by:\n",
    "\n",
    "$$ F = \\int_V  f(\\eta, \\nabla \\eta, \\nabla^{2}\\eta, ...) dV $$\n",
    "\n",
    "Where:\n",
    "\n",
    "- $ \\nabla \\eta $ represents the gradient of $ \\eta $: $ \\frac{\\partial \\eta}{\\partial x_i} $.\n",
    "- $ \\nabla^{2} \\eta $ is its second derivative, a tensor: $ \\frac{\\partial^2 \\eta}{\\partial x_i \\partial x_j} $.\n",
    "\n",
    "Expanding the free energy density, $ f $, using Taylor series, we get:\n",
    "\n",
    "$$\n",
    "f(\\eta, \\nabla \\eta, \\nabla^{2}\\eta, ...) = f(\\eta) + \n",
    "     \\frac{\\partial f}{\\partial \\nabla \\eta} \\cdot \\nabla \\eta +\n",
    "     \\frac{\\partial f}{\\partial \\nabla^2 \\eta} \\colon \\nabla^2 \\eta +\n",
    "     \\frac{1}{2!}\\frac{\\partial^2 f}{\\partial \\nabla \\eta \\partial \\nabla \\eta} \\colon (\\nabla \\eta)^2 + ...\n",
    "$$\n",
    "\n",
    "Neglecting higher order terms and defining:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\frac{\\partial f}{\\partial \\nabla \\eta} &\\equiv \\alpha_i, \\\\\n",
    "    \\frac{\\partial f}{\\partial \\nabla^2 \\eta} &\\equiv \\beta_{ij}, \\\\\n",
    "    \\frac{\\partial^2 f}{\\partial \\nabla \\eta \\partial \\nabla \\eta} &\\equiv \\gamma_{ij},\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We can rewrite our expansion as:\n",
    "\n",
    "$$\n",
    "f(\\eta, \\nabla \\eta, \\nabla^{2}\\eta) = f(\\eta) + \n",
    "     \\alpha_i \\cdot \\frac{\\partial \\eta}{\\partial x_i} +\n",
    "     \\beta_{ij} \\colon \\frac{\\partial^2 \\eta}{\\partial x_i \\partial x_j} +\n",
    "     \\frac{1}{2!}\\gamma_{ij} \\colon (\\frac{\\partial \\eta}{\\partial x_i})(\\frac{\\partial \\eta}{\\partial x_j})\n",
    "$$\n",
    "\n",
    "Assuming the system has inversion symmetry, $ \\alpha_i = 0 $.  The thrid and forth terms can be combined by reforming the third term using intergration by parts:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\int_V \\beta_{ij} \\colon \\frac{\\partial^2 \\eta}{\\partial x_i \\partial x_j} dV =\n",
    "\\int_V \\beta_{ij} \\colon \\frac{\\partial}{\\partial x_i}(\\frac{\\partial \\eta}{\\partial x_j}) dV =\n",
    "\\left[\\beta_{ij}\\frac{\\partial \\eta}{\\partial x_j}\\right]_s - \\int_V \\frac{\\partial \\beta_{ij}}{\\partial x_i}\\frac{\\partial \\eta}{\\partial x_j} dV = \n",
    "\\left[\\beta_{ij}\\frac{\\partial \\eta}{\\partial x_j}\\right]_s - \\int_V \\frac{\\partial \\beta_{ij}}{\\partial \\eta} \\frac{\\partial \\eta}{x_i} \\frac{\\partial \\eta}{\\partial x_j} dV\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Neglecting the boundary term, the total free energy is:\n",
    "\n",
    "$$ F = \\int_V  f(\\eta) + \\left[ \\frac{\\gamma_{ij}}{2} - \\frac{\\partial \\beta_{ij}}{\\partial \\eta} \\right] \\colon (\\frac{\\partial \\eta}{\\partial x_i})(\\frac{\\partial \\eta}{\\partial x_j})  dV $$\n",
    "\n",
    "Defining:\n",
    "\n",
    "$$\n",
    "\\left[ \\frac{\\gamma_{ij}}{2} - \\frac{\\partial \\beta_{ij}}{\\partial \\eta} \\right] \\equiv \\frac{\\kappa_{ij}}{2},\n",
    "$$\n",
    "\n",
    "We have:\n",
    "\n",
    "$$ F = \\int_V  f(\\eta) + \\frac{\\kappa_{ij}}{2} \\colon (\\frac{\\partial \\eta}{\\partial x_i})(\\frac{\\partial \\eta}{\\partial x_j}) dV $$\n",
    "\n",
    "For isotropic systems dominated by a single parameter, $ \\kappa $, this simplifies to:\n",
    "\n",
    "$$ F = \\int_V  f(\\eta) + \\frac{\\kappa}{2} \\| \\nabla \\eta \\|^2 dV $$\n",
    "\n",
    "For systems with multiple order parameters, we generalize as:\n",
    "\n",
    "$$ F = \\int_V \\left[ f(\\boldsymbol{\\eta}) + \\frac{\\kappa}{2} \\sum_{i=1}^N \\|\\nabla \\eta_i\\|^2 \\right] dV $$\n",
    "\n",
    "**Reference**: [Phase Field Modelling Video](https://www.youtube.com/watch?v=8SUhwZYsc1o)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a1db5bb",
   "metadata": {},
   "source": [
    "#### Cahn-Hillard equation\n",
    "\n",
    "If we assume the free energy density is the function of species concentrations instead of order parameters, we can derive Cahn-Hillard equation, which is a conserative form.\n",
    "The total free energy in term of a species concentration $c$ is expressed as\n",
    "$$ F = \\int_V  f(c, \\nabla c, \\nabla^{2}c, ...) dV $$\n",
    "\n",
    "Following the same procedure discussed above, we can get a simplied free energy expression as \n",
    "\n",
    "$$ F = \\int_V  f(c) + \\frac{\\kappa}{2} \\| \\nabla c \\|^2 dV $$\n",
    "\n",
    "The chemial potential $\\mu$ is the variational derivative of $F$ with respect to $c$\n",
    "\n",
    "$$ \\mu = \\frac{\\delta F}{\\delta c}$$\n",
    "\n",
    "The rate of change of the concentration is determined by the gradient of the chemical potential\n",
    "\n",
    "$$ \\frac{\\partial c}{\\partial t} = -\\nabla J =\\nabla \\cdot M \\nabla \\mu$$\n",
    "\n",
    "where $J$ is the flux of the chemical potential and $M$ is the mobility coefficient. Substuting the expresion of the totoal free energy, we can obtain the conserative **Cahn-Hillard equation**\n",
    "\n",
    "$$ \\frac{\\partial c}{\\partial t} = M \\nabla^2 (\\frac{\\partial f}{\\partial c} - \\kappa \\nabla^2c)$$\n",
    "\n",
    "Reorganzing the derivatives, we can obtain the governing equation for diffusion phenomena, considering the effects of both concentration gradient and curvature:\n",
    "\n",
    "$$ \\frac{\\partial c}{\\partial t} = M \\frac{\\partial^2 f}{\\partial c^2} \\nabla^2c - \\kappa M \\nabla^4c$$\n",
    "\n",
    "It's enlightening to compare the Cahn-Hilliard equation with the traditional diffusion equation deduced from Fick's law: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\frac{\\partial c}{\\partial t} = M \\frac{\\partial^2 f}{\\partial c^2} \\nabla^2c - \\kappa M \\nabla^4c \n",
    "    \\\\\n",
    "    \\frac{\\partial c}{\\partial t} = D \\nabla^2c\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "From the comparison, it's evident that the coefficient $ M \\frac{\\partial^2 f}{\\partial c^2} $ parallels the diffusion coefficient $ D $. Moreover, the Cahn-Hilliard equation introduces an additional term, $ - \\kappa M \\nabla^4c $, which accounts for modulating the effects linked to concentration curvature.\n",
    "\n",
    "**Reference**: [Phase Field Modelling Video](https://www.youtube.com/watch?v=X5k_WtQIeDo&list=PLaikkcmmwwxL4Ujl4wjPoPp_AHjVZgigc&index=66)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd332c49",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}