{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-dimensional phase-field model (Allen-Cahn equation)\n", "This python code was developed by Yamanaka research group of Tokyo University of Agriculture and Technology in August 2019." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem\n", "Let us consider an example: growth of phase B in phase A. \n", "This notebook shows two-dimensional phase-field simulation of the growth of phase B using Allen-Cahn equation. \n", "\n", "For details of the model, please see reference (*S. M. Allen and J. W. Cahn, \"A Microscopic Theory for Antiphase Boundary Motion and Its Application to Antiphase Domain Coarsening,\" Acta Metallurgica, Vol. 27 (1979), pp. 1085–1095.*).\n", "## Formulation\n", "### 1. Order parameter\n", "An non-conserved order parameter $\\phi$ (hereafter, we call it as \"phase-field variable\".) is used. The phase-field variable $\\phi$ is defined as $\\phi = 0$ in phase A and $\\phi = 1$ in phase B. The phase-field variable $\\phi$ smoothly changes from 0 to 1 in an \"diffuse\" interfacial region. \n", "### 2. Total free energy\n", "The total Gibbs free energy of the system is defined by\n", "$$\n", "G = \\int_{V} (g_{chem}(\\phi) + g_{doub}(\\phi) + g_{grad}(\\nabla\\phi)) dV\n", "$$\n", "where, $g_{chem}$, $g_{doub}$, and $g_{grad}$ are the chemical free energy, the double-well potential energy and the gradient energy densities, respectively. In this source code, these energy densities are defined as follows: \n", "$$\n", "g_{chem} = p(\\phi)g_{B} + (1-p(\\phi))g_{A}\n", "$$\n", "$$\n", "g_{doub} = Wq(\\phi)\n", "$$\n", "$$\n", "g_{grad} = \\frac{a^{2}}{2}|\\nabla\\phi|^{2}\n", "$$\n", "where $g_{A}$ and $g_{B}$ are Gibbs free energy densities of pure phase A and pure phase B, respectively. Here, $g_{A}$ and $g_{B}$ are assumed to be constant. \n", "\n", "$p(\\phi)$ and $q(\\phi)$ are an interpolation function and the double-well potential function. These functions are often defined as follows: \n", "$$\n", "p(\\phi) = \\phi^{2}(3-2\\phi)\n", "$$\n", "$$\n", "q(\\phi) = \\phi^{2}(1-\\phi)^{2}\n", "$$\n", "\n", "$W$ and $a$ are the height of double-well potential energy and the gradient energy coefficient, respectively, which are given by\n", "$$\n", "W = \\frac{6\\sigma\\delta}{b}\n", "$$\n", "$$\n", "a = \\sqrt{\\frac{3\\sigma b}{\\delta}}\n", "$$\n", "where $\\sigma$ is the interfacial energy of an interface between phase A and B. $\\delta$ is the thickness (length in one-dimension) of the diffuse interface. $b$ is given as $b = 2\\tanh^{-1}(1-2\\lambda)$ where $\\lambda$ is a constant which is used for determining the thickness of interfacial region. In this source code, we use $\\lambda$ = 0.1 and then $b$ = 2.2. \n", "\n", "### 3. Time evolution equation (Allen-Cahn equation)\n", "The time evolution of the phase-field variable (that describes the migration of the interface) is given by Allen-Cahn equation by assuming that the total free energy of the system, $G$, decreases monotonically with time (i.e. The second law of thermodynamics): \n", "$$\n", "\\frac{\\partial \\phi}{\\partial t} = -M_{\\phi}\\frac{\\delta G}{\\delta \\phi}=-M_{\\phi}\\left(\\frac{\\partial p(\\phi)}{\\partial \\phi}(g_{B}-g_{A})+W\\frac{\\partial q(\\phi)}{\\partial\\phi}-a^{2}\\nabla^{2}\\phi\\right)\n", "$$\n", "where $M_{\\phi}$ is the mobility of phase-field variable and is given by: \n", "$$\n", "M_{\\phi} = \\frac{\\sqrt{2W}}{6a}M\n", "$$\n", "Here, $M$ is a \"physical\" mobility of interface between phase A and B. \n", "The term $g_{B}-g_{A}$ in the Allen-Cahn equation corresponds to the driving force of the growth of phase B. \n", "\n", "Note that the Euler-Lagrange equation given by the following equation is used to calculate the functional derivative of $G$ with respect to the phase-field variable as:\n", "$$\n", "\\frac{\\delta G}{\\delta \\phi}=\\frac{\\partial g}{\\partial \\phi}-\\nabla\\cdot\\frac{\\partial g}{\\partial (\\nabla \\phi)}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Programming " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### import libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib nbagg\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from matplotlib import animation\n", "import math" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### set parameters and physical values " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "nx, ny = 32, 32 # number of computational grids\n", "dx, dy = 0.5e-6, 0.5e-6 # spacing of computational grid [m]\n", "eee = 1.0e+6 # driving force of growth of phase B: g_A - g_B [J/m3]\n", "sigma = 1.0 # ineterfacial energy [J/m2]\n", "delta = 4.*dx # interfacial thickness [m]\n", "amobi = 4.e-14 # interfacial mobilitiy [m4/(Js)]\n", "ram = 0.1 # paraneter which deternines the interfacial area\n", "bbb = 2.*np.log((1.+(1.-2.*ram))/(1.-(1.-2.*ram)))/2. # The constant b = 2.1972" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### calculate phase-field parameters ($a, W$ and $M_{\\phi}$)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "aaa = np.sqrt(3.*delta*sigma/bbb) # gradient energy coefficient \"a\"[(J/m)^(1/2)]\n", "www = 6.*sigma*bbb/delta # potential height W [J/m3]\n", "pmobi = amobi*math.sqrt(2.*www)/(6.*aaa) # mobility of phase-field [m3/(Js)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### define time increment and total number of time steps" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "dt = dx*dx/(5.*pmobi*aaa*aaa)/2 # time increment for a time step [s]\n", "nsteps = 500 # total number of time step" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### declare arrays for phase field variable" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "p = np.zeros((nsteps,nx,ny)) # phase-field variable" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### set initial nuclei of B phase\n", "The initial nuclei of B phase (region of $\\phi = 1$) is located at the origin of the computational domain. \n", "The shape of the nuclei is assumed to be a circle. \n", "\n", "The initial profile of the phase-field variable is calculated using the equilibrium profile of the phase-field variable as: \n", "$$\n", "\\phi = \\frac{1}{2}\\left(1-\\tanh \\frac{\\sqrt{2W}}{2a}(r-r_{0}) \\right)\n", "$$\n", "where $r_{0}$ denotes the initial position of the interface. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "r_nuclei = 5.*dx # radius of the initial B phase\n", "for i in range(0,nx):\n", " for j in range(0,ny):\n", " r = np.sqrt( (i *dx)**2 +(j*dy)**2 ) - r_nuclei\n", " p[0,i,j] = 0.5*(1.-np.tanh(np.sqrt(2.*www)/(2.*aaa)*r))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### define function for solving Allen-Cahn equation\n", "Allen-Cahn equation is discretized by simple finite difference method in this program. \n", "\n", "The 1st-order Euler method is used for time-integration. The 2nd-order central finite difference method is used for spatial derivatives. The discretized Allen-Cahn equation is given as: \n", "$$\n", "\\phi^{t+\\Delta t}_{i,j} = \\phi^{t}_{i,j} + M_{\\phi}\n", "\\left[ 4W\\phi^{t}_{i,j}\\left(1-\\phi^{t}_{i,j}\\right)\\left(\\phi^{t}_{i,j}-\\frac{1}{2}+\\frac{3}{2W}(g_{A}-g_{B})\\right)\n", "+a^{2}\\left( \\frac{\\phi^{t}_{i+1,j}-2\\phi^{t}_{i,j}+\\phi^{t}_{i-1,j}}{(\\Delta x)^{2}} \n", " + \\frac{\\phi^{t}_{i,j+1}-2\\phi^{t}_{i,j}+\\phi^{t}_{i,j-1}}{(\\Delta y)^{2}}\\right)\n", "\\right]\\Delta t\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def do_timestep(p):\n", " for t in range(nsteps-1):\n", " for j in range(ny):\n", " for i in range(nx):\n", " ip = i + 1\n", " im = i - 1\n", " jp = j + 1\n", " jm = j - 1\n", " if ip > nx - 1:\n", " ip = nx -1\n", " if im < 0:\n", " im = 0\n", " if jp > ny - 1:\n", " jp = ny -1\n", " if jm < 0:\n", " jm = 0\n", " p[t+1,i,j] = p[t,i,j] + pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3./(2.*www)*eee)+ aaa*aaa*((p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Time integration of Allen-Cahn equation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('