Phase field: practical

B. Appolaire [1]
Y. Le Bouar [2]

[1] Institut Jean Lamour, Université de Lorraine – CNRS, France
[2] Laboratoire d'Étude des Microstructures, CNRS – ONERA, France

Feb. 2019


Free energy functional


We use the following free energy functional: $$ \begin{equation} F=\int f(\phi)+\frac{\alpha}{2}|\nabla \phi|^2 \ dV \label{Free_energy_functional} \end{equation} $$

where the homogeneous free energy density is $$ \begin{equation} f(\phi)= 16 \, \Delta f \, \phi^2 \ (1-\phi)^2 \label{dble_well} \end{equation} $$

The equilibrium profile of a 1D interface is $$ \begin{equation} \phi_{eq}(x)= \frac12 \left(1+\tanh\left(\frac{2x}{\delta}\right)\right) \label{tanh} \end{equation} $$

where \( \delta=\sqrt{\alpha/(2\Delta f)} \) is the interface width.

The associated interfacial energy is $$ \begin{equation} \gamma = \frac23 \sqrt{2 \,\alpha \, \Delta f} \label{_auto1} \end{equation} $$

Introducing the grid spacing \( d \), the discretized nondimensional free energy \( \tilde{F}={F}/(d^3 \, \Delta f) \) becomes: $$ \begin{equation} \tilde{F} =\sum_n\, \tilde{f}(\phi_n)+\frac{\tilde{\alpha}}{2}|\tilde{\nabla} \phi_n|^2 \label{free_energy_discrete_nondim} \end{equation} $$

where \( \phi_n \) denotes the value of field \( \phi \) at node \( n \), \( n\in\mathbb{N} \), \( \tilde{f}(\phi)= f(\phi)/\Delta f \), \( \tilde{\alpha}=\alpha/(\Delta f \, d^2) \) and \( \tilde{\nabla} \phi \) is the discrete nondimensional gradient operator, for example in 1D \( \phi((n+1)d)-\phi(nd) \).

The nondimensional interface free energy and interfacial width are related to the nondimensional gradient term in the following way: $$ \begin{align} \tilde{\gamma} = \frac{\gamma}{d\,\Delta f} &=\frac{ 2\sqrt{2 }}3 \, \sqrt{\tilde{\alpha}} \label{_auto2}\\ \tilde{\delta} = \frac{\delta}{d} &= \sqrt{\frac{\tilde{\alpha}}{2}} \label{delta_tilde} \end{align} $$

\( \tilde{\delta} \) should be equal to 3 to 5 for numerical reasons (grid pinning).



Allen-Cahn equation


Assuming a classical relaxation dynamics for the order parameter, the evolution equation reads: $$ \begin{equation} \frac{d\phi}{dt} = - L \, \frac{\delta F}{\delta \phi} \label{_auto3} \end{equation} $$

where the functional derivative is: $$ \begin{equation} \frac{\delta F}{\delta \phi}= f'(\phi)-\alpha \Delta \phi \label{_auto4} \end{equation} $$

We assume that the characteristic decay time \( t_{exp} \) for a long wavelength fluctuation around \( \phi=0 \) is know. Using \( f(\phi)\approx \frac{1}{2} f''(0)\,\phi^2 \), and neglecting the heterogeneous term, we have: $$ \begin{equation} \frac{d\phi}{dt} \approx - L \, f''(0) \, \phi \label{_auto5} \end{equation} $$

and finally \( L= 1/ (f''(0) \, t_{exp}) \).

Using the dimensionless time \( \tilde{t}=t/t_{exp} \), the nondimensional Allen-Cahn equation writes: $$ \begin{equation} \frac{d\phi}{d\tilde{t}} = - \tilde L \,\big(\tilde{f}'(\phi)-\tilde{\alpha} \tilde{\Delta} \phi\big) \label{allencahn_nondim} \end{equation} $$

where we have used the kinetic parameter \( \tilde L = {1}/{\tilde{f}''(0)} \) and the discrete laplacian \( \tilde{\Delta}= d^2\Delta \).

As you may know, life is often much easier in Fourier space. The Allen-Cahn equation in Fourier space writes: $$ \begin{equation} \frac{d\hat\phi}{d\tilde{t}} = - \tilde L \,\left(\,\left\{\tilde{f}'(\phi(\tilde{r}))\right\}_{\tilde{k}} + \tilde{k}^2 \, \tilde{\alpha} \, \hat{\phi} \, \right) \label{_auto6} \end{equation} $$

where the Fourier transform of any field \( \psi \) is denoted either as \( \hat{\psi} \) or \( \left\{\psi\right\}_{\tilde{k}} \).

The following explicit first order Euler spectral scheme is already coded: $$ \begin{equation} \hat\phi(\tilde{t}+\tilde{dt}) = \hat\phi(\tilde{t}) - \tilde L \, \tilde{dt} \,\left(\, \left\{\tilde{f}'(\phi(\tilde{r}))\right\}_{\tilde{k},\tilde{t}} + \tilde{k}^2 \, \tilde{\alpha} \, \hat{\phi}(\tilde{t}) \, \right) \label{eulerspectral} \end{equation} $$

where all quantities on the right hand side are calculated at the previous time step \( \tilde{t} \).



Questions


To answer the following questions,

Beware to indent using spaces rather than tabs (one level = 3 spaces) when modifying the python scripts.

a) How can we ensure to have about 3 grid points inside the interface?

Hint.

Use Eq. \eqref{delta_tilde}. You can check your result by playing with parameter \( \tilde{\alpha} \) considering a plane interface (see next question).

b) Compute an equilibrium flat profile and compute the nondimensional interfacial free energy.

Hint.

Set alpha to the right value, and run a calculation considering a planar precipitate in the middle of the box by decommenting init_interface with the proper parameters (and of course commenting all the other functions changing phi0).

c) Start a simulation with a circle and study the evolution of the average phase field. Explain the result.

d) What is the maximum time step that can be used?

Hint.

Proceed by trial and error: indeed, the formal Von Neumann stability analysis cannot be applied because of the nonlinear term. It may be easier to consider first the case of the lamellar precipitate rather than the circular one. Once a critical value for the time step has been obtained form the lamellar precipitate, check it with the circular one.

e) Code a semi implicit spectral scheme and analyze its stability and accuracy (you may use the simulation of a shrinking circle).

Hint.

The semi-implicit scheme is obtained by evaluating \( \hat{\phi} \) on the right hand side of Eq. \eqref{eulerspectral} at time \( t+dt \).

f) What should be modified to turn dissolution into growth?

g) Start the simulation with the constant value of 0.5 for phi plus a small noise of amplitude 0.01. What is the obtained microstructural evolution?

h) Homework: Test of the spatial discretization scheme. I want to repeat a calculation after increasing the number of grid points in the interface, keeping the real width of the interface, the interfacial energy, the decay time and the real box size constant. What shall I do?

i) Homework: Errors related to the real interfacial width. I want to repeat the calculation after decreasing the real interfacial width, but keeping constant the number of grid points in the interface, the interfacial energy, the decay time and the real box size. What shall I do?

j) Code the Cahn-Hilliard equation in a semi-implicit scheme.



The python scripts


#!/usr/bin/env python

'''      Solve the Allen-Cahn & Cahn-Hillard equations

# (The MIT License)

# Copyright (c) 2013, 2017 
#          Benoit Appolaire <benoit.appolaire@onera.fr>
#        & Yann Le Bouar    <yann.lebouar@onera.fr>

# Permission is hereby granted, free of charge, to any person 
# obtaining a copy of this software and associated documentation 
# files (the "Software"), to deal in the Software without restriction, 
# including without limitation the rights to use, copy, modify, merge, 
# publish, distribute, sublicense, and/or sell copies of the Software, 
# and to permit persons to whom the Software is furnished to do so, 
# subject to the following conditions:

# The above copyright notice and this permission notice shall be 
# included in all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, 
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY 
# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, 
# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
'''

#1
#______________________________________________________________________________
#                         Import modules/libraries
#______________________________________________________________________________

import sys              #system basic routines (may be usefull for debugging)
from scipy import *     #scientific library including fast fourier transform
from pylab import *     #plotting a la matlab with matplotlib
import numpy.random as rng #for random numbers
import operator         #export functions for the intrinsic operators


#2
#______________________________________________________________________________
#                     Define a few functions/subroutines
#______________________________________________________________________________


#3
#..............................................................................
def Egrad(phi_,kx_,ky_,alpha):
   '''
      Gradient energy density
   '''
#..............................................................................

   tmp = 0.5*alpha*(kx_**2+ky_**2)*abs(phi_)**2
   return tmp.real
   
#4
#..............................................................................
def g(phi):
   '''
      Double well potential where the energy barrier is 1
   '''
#..............................................................................

   return 16.*phi**2.*(1.-phi)**2.
   
#5
#..............................................................................
def dg(phi):
   '''
      Derivative of the double well potential $g(\phi) = 16*\phi^2(1-\phi)^2$
      i.e. $\frac{dg}{dphi}$ (Non linear source term in the pde)
   '''
#..............................................................................

   return 2.*phi*(1.-phi)*(1.-2.*phi)*16.

#6
#..............................................................................
def ddg(phi):
   '''
      Second Derivative of the double well potential 
      $g(\phi) = 16\phi^2(1-\phi)^2$ i.e. $\frac{d^2g}{dphi^2}$
   '''
#..............................................................................

   return 32.*(1.-6.*phi*(1.-phi))

#7
#..............................................................................
def discretization(L=1.,N=10):
   '''
      Discretize with N nodes the real and reciprocal spaces 
      along one direction with length 2*L 
   '''
#..............................................................................

   return mgrid[-L:L:2*L/N], fftfreq(N,d=L/N)

#8
#..............................................................................
def init_cosx(x_,y_,scalex=1.):
   '''
      Initialization of a field with cosine along x
   '''
#..............................................................................
   
   return (1.+cos(2.*pi*x_/scalex))/2.

#9
#..............................................................................
def init_circle(x_,y_,centerx=0.,centery=0.,radius=1.,val1=1.,val2=0.):
   '''
      Initialization of a field with a disk centered at (centerx,centery)
      and with a given radius
   '''
#..............................................................................
   
   return where( (x_-centerx)**2.+(y_-centery)**2.-radius**2.<=0., val1, val2 )

#10
#..............................................................................
def init_interface(x_,y_,pos1=0,pos2=1,val1=1.,val2=0.):
   '''
      Initialization of a field with two interfaces normal to x
   '''
#..............................................................................
   
   return where( (( x_ >= pos1) & ( x_ < pos2)), val1, val2 )

#11
#..............................................................................
def iterwhite():
   '''
      Generate a sequence of samples of white noise.
      Generates a never-ending sequence of floating-point values.
      
      # Copyright (c) 2011 Leif Johnson <leif@leifjohnson.net>
      #
      # Permission is hereby granted, free of charge, to any person 
      # obtaining a copy of this software and associated documentation 
      # files (the "Software"), to deal in the Software without restriction, 
      # including without limitation the rights to use, copy, modify, merge, 
      # publish, distribute, sublicense, and/or sell copies of the Software, 
      # and to permit persons to whom the Software is furnished to do so, 
      # subject to the following conditions:
      #
      # The above copyright notice and this permission notice shall be 
      # included in all copies or substantial portions of the Software.
   '''
#..............................................................................

   while True:
      for n in rng.randn(100):
         yield n

#12
#..............................................................................
def _asarray(source, shape):
   '''
      # Copyright (c) 2011 Leif Johnson <leif@leifjohnson.net>
      #
      # Permission is hereby granted, free of charge, to any person 
      # obtaining a copy of this software and associated documentation 
      # files (the "Software"), to deal in the Software without restriction, 
      # including without limitation the rights to use, copy, modify, merge, 
      # publish, distribute, sublicense, and/or sell copies of the Software, 
      # and to permit persons to whom the Software is furnished to do so, 
      # subject to the following conditions:
      #
      # The above copyright notice and this permission notice shall be 
      # included in all copies or substantial portions of the Software.
   '''
#..............................................................................

   noise = source()
   if shape is None:
      return noise.next()
   count = reduce(operator.mul, shape)
   return asarray([noise.next() for _ in range(count)]).reshape(shape)

#13
#..............................................................................
def white(shape=None):
   '''
      Generate white noise.

      shape: If given, returns a numpy array of white noise with this shape. 
      If not given, return just one sample of white noise.

      # Copyright (c) 2011 Leif Johnson <leif@leifjohnson.net>
      #
      # Permission is hereby granted, free of charge, to any person 
      # obtaining a copy of this software and associated documentation 
      # files (the "Software"), to deal in the Software without restriction, 
      # including without limitation the rights to use, copy, modify, merge, 
      # publish, distribute, sublicense, and/or sell copies of the Software, 
      # and to permit persons to whom the Software is furnished to do so, 
      # subject to the following conditions:
      #
      # The above copyright notice and this permission notice shall be 
      # included in all copies or substantial portions of the Software.
   '''
#..............................................................................

   return _asarray(iterwhite, shape)



#14
#..............................................................................
def headline():

#..............................................................................

   string  = '# time'
   string += '       fhom       grad       Ftot  '
   string += '       dphi         <phi>    '
   string += '   stored in'
   print string


#..............................................................................
def print_store(t,iplot,nl,nc,step,isteps,kx_,ky_,phi0,phi,phi_,alpha):
   '''
      Monitoring the evolution with global quantities (average, min, max)
      and snaphots of the phase field
   '''
#..............................................................................
      
   fhom    = g(phi.real)
   grad    = Egrad(phi_,kx_,ky_,alpha)
   E_grad  = real(average(grad))
   E_homo  = nmax*nmax*average(fhom)
   E_tot   = E_grad + E_homo

   fname   = 'phi%02d'%iplot
   savez(fname,phi=phi.real)

   subplot(nl,nc,iplot+1)
   imshow(phi.real,interpolation='nearest',origin='lower',cmap=cm.gray,\
                                                        vmin=0.,vmax=1.)
   #contour(phi.real,levels=[0.5],colors='r')
   title('t = %.2f'%(step*pdt))
   xticks([])
   yticks([])

   isteps.append(step)

   string  = '%.2e'%t
   string += '   %.3e  %.3e  %.3e  '%(E_grad, E_homo, E_tot)
   string += '  %.3e    %.3e   '%((phi-phi0).real.max(),average(phi.real))
   string += '  %s'%(fname+'.npz')
   print string

   iplot += 1

   return iplot



#15
#______________________________________________________________________________

#     This is the main program which will be run when invoking the script

if __name__ == "__main__":
#______________________________________________________________________________
   
   #This is only for automatic integration by doconce of the
   # previous if statement into a jupyter notebook cell to run
   dummy_for_ipynb = True


   #16
   #....................................................
   #(0) A few data

   #Phase field gradient (nondimensional)
   alpha   = 1.
   
   #Kinetic coefficient (nondimensional)
   kin     = 1./ddg(0)
   
   print '# kin = ', kin



   #17
   #....................................................
   #(1) Space discretization
   nmax    = 100
   xmax    = 1.*nmax

   x,kx    = discretization(L=xmax,N=nmax)
   y,ky    = discretization(L=xmax,N=nmax)
   
   x_,y_   = meshgrid(x,y)
   kx_,ky_ = meshgrid(kx,ky)

   kx_     = 2*pi*kx_
   ky_     = 2*pi*ky_


   #18
   #....................................................
   #(2) Time discretization
   tmax    = 1.
   pdt     = 1e-3
   nsteps  = int(tmax/pdt)
   steps   = range(1,nsteps+1)

   #   Times when to plot
   nplots  = 10
   tplots  = linspace(0.,tmax,nplots+1)
   #tplots = [80,82,84,86,88,90,92,94,96,98,100,tmax]
   iplot   = 0


   #19
   #....................................................
   #(3) Initial microstructure
   
   phi0    = zeros(shape(x_))

   phi0   += init_interface(x_,y_,pos1=-x.max()/2.,\
                            pos2=x.max()/2.,val1=1.,val2=0.)
   #phi0   += init_circle(x_,y_,centerx=(0.),centery=(0.),\
   #                      radius=(x.max()/4.), val1=1., val2=0.)
   #phi0   += init_circle(x_,y_,centerx=(x.max()/2.),\
   #                      centery=(x.max()/6.),radius=(x.max()/3.),\
   #                      val1=1., val2=0.)
   #phi0    = 0.5 + 0.01*white(shape(x_))
   #phi0   += 0.01*white(shape(x_))


   #20
   phi     = phi0
   phi_    = fft2(phi)
   grad    = zeros(shape(x_))
   fhom    = zeros(shape(x_))
   phi_av  = [average(phi)]
   phi_min = [phi.real.min()]
   phi_max = [phi.real.max()]


   #21
   #....................................................
   #(4) Prepare the printing and final plots
   
   headline()

   nc      = 4
   nl      = (len(tplots)-1)/nc + 1
   fig1    = figure(1,figsize=(14,3*nl))

   step    = 0
   isteps  = []
   iplot   = print_store(step*pdt,iplot,nl,nc,step,isteps,\
                               kx_,ky_,phi0,phi,phi_,alpha)


   #22
   #....................................................
   #(5) Time integration of the PDE

   for step in steps:
      
      # (5a) Update the previous field
      phi0 = phi
      # (5b) Solve in Fourier space using FFT
      dg_  = fft2(dg(phi))
      phi_ = fft2(phi)
      #---------------------  Allen-Cahn (explicit) ----------------------
      phi_ = phi_ - kin*pdt*(dg_+(kx_**2.+ky_**2.)*alpha*phi_)
      #-------------------------------------------------------------------
      # (5c) Go back to real space
      phi  = ifft2(phi_)
      # (That's all !!!!!)
   
      # (5d) Compute quantities for the time evolution plot
      phi_av.append(average(phi.real))
      phi_min.append(phi.real.min())
      phi_max.append(phi.real.max())

      # (5e) Snapshots of the evolution
      if step == int(tplots[iplot]/pdt):
         iplot = print_store(step*pdt,iplot,nl,nc,step,isteps,\
                                   kx_,ky_,phi0,phi,phi_,alpha) 


   #23
   #....................................................
   #(6) Plot the time evolution of the average and extrema

   subplots_adjust(bottom=0.05,left=0.05,right=0.95,top=0.95,wspace=0.4)


   fig2 = figure(2)
   clf()
   axes([0.14,0.12,0.82,0.83])

   trange = array([0]+steps)*pdt
   #phi_av = array(phi_av)*100.

   plot(trange,phi_av, 'k-', label='avg')
   plot(trange,phi_min,'b-', label='min')
   plot(trange,phi_max,'r-', label='max')
   legend( loc='center right' )
   axis([0.,tmax,-0.05,1.05])

   xlabel('time',   fontsize=20)
   ylabel(r'$\phi$',fontsize=20)


   # plot dots corresponding to the snapshots
   for i in isteps:
      plot(trange[i],phi_av[i], 'ko',markersize=10)
      plot(trange[i],phi_min[i],'bo',markersize=10)
      plot(trange[i],phi_max[i],'ro',markersize=10)

   savetxt('time_table.txt',trange[isteps])
   fig1.savefig('snapshots.png')
   fig2.savefig('conc_extrema.pdf')

   print
   print 'Done'
   print
   
   show()

   #24
   #....................................................



Plotting a snapshot

It can be useful to have a look at a particular snapshot that has been saved in an npz file. This can be done using PlotSnapshot.py which plots side by side the 2D field and the profile along an horizontal line at position given by iy (in nodes).


#!/usr/bin/env python

'''   Visualization of npz files generated by AllenCahnHilliard_1.py

# (The MIT License)
#
# Copyright (c) 2013, 2017 
#          Benoit Appolaire <benoit.appolaire@onera.fr>
#        & Yann Le Bouar    <yann.lebouar@onera.fr>
#
# Permission is hereby granted, free of charge, to any person 
# obtaining a copy of this software and associated documentation 
# files (the "Software"), to deal in the Software without restriction, 
# including without limitation the rights to use, copy, modify, merge, 
# publish, distribute, sublicense, and/or sell copies of the Software, 
# and to permit persons to whom the Software is furnished to do so, 
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be 
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, 
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY 
# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, 
# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

'''

#1
import sys
from pylab import *
from numpy import load


#2
#______________________________________________________________________________
#      This is the main program which will be run when invoking the script

if __name__ == "__main__":
#______________________________________________________________________________


   name = 'phi10'

   res = load(name+'.npz')
   phi = res['phi']

   iy  = 50

   #3
   figure(1,figsize=(11,5))
   ax1 = axes([0.05,0.15,0.4,0.8])
   imshow(phi,interpolation='nearest',origin='lower',\
                         cmap=cm.gray,vmin=0.,vmax=1.)
   axhline(iy,color='r')
   xticks([])
   yticks([])

   #4
   ax2 = axes([0.53,0.15,0.44,0.8])
   plot(phi[iy,:],'r-o',linewidth=1.5)
   axis([0,len(phi)-1,-0.01,1.01])
   xlabel('$x$',fontsize=24)
   ylabel('$\phi$',fontsize=24)


   savefig(name+'.png')
   show()

#5

© 2019, B. Appolaire, Y. Le Bouar. CC BY-NC-SA