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).
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} \).
To answer the following questions,
a) How can we ensure to have about 3 grid points inside the interface?
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.
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?
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).
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.
#!/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
#....................................................
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