The spinodal decomposition phenomenon is a spontaneous separation of an initially homogenous mixture into two distinct regions of different properties (spin-up/spin-down, component A/component B). It is a “barrierless” phase separation process, such that under the right thermodynamic conditions, any fluctuation, no matter how small, will tend to grow. This is in contrast to nucleation, where a fluctuation must exceed some critical magnitude before it will survive and grow. Spinodal decomposition can be described by the “Cahn-Hilliard” equation (also known as “conserved Ginsberg-Landau” or “model B” of Hohenberg & Halperin)
where is a conserved order parameter, possibly representing alloy composition or spin. The double-well free energy function penalizes states with intermediate values of between 0 and 1. The gradient energy term , on the other hand, penalizes sharp changes of . These two competing effects result in the segregation of into domains of 0 and 1, separated by abrupt, but smooth, transitions. The parameters and determine the relative weighting of the two effects and is a rate constant.
We can simulate this process in FiPy with a simple script:
>>> from fipy import *
>>> if __name__ == "__main__": ... nx = ny = 1000 ... else: ... nx = ny = 20 >>> mesh = Grid2D(nx=nx, ny=ny, dx=0.25, dy=0.25) >>> phi = CellVariable(name=r"$\phi$", mesh=mesh)
We start the problem with random fluctuations about
>>> phi.setValue(GaussianNoiseVariable(mesh=mesh, ... mean=0.5, ... variance=0.01))
FiPy doesn’t plot or output anything unless you tell it to:
>>> if __name__ == "__main__": ... viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)
FiPy would automatically interpolate D * a**2 * (1 - 6 * phi * (1 - phi)) onto the Faces, where the diffusive flux is calculated, but we obtain somewhat more accurate results by performing a linear interpolation from phi at Cell centers to PHI at Face centers. Some problems benefit from non-linear interpolations, such as harmonic or geometric means, and FiPy makes it easy to obtain these, too.
>>> PHI = phi.getArithmeticFaceValue() >>> D = a = epsilon = 1. >>> eq = (TransientTerm() ... == DiffusionTerm(coeff=D * a**2 * (1 - 6 * PHI * (1 - PHI))) ... - DiffusionTerm(coeff=(D, epsilon**2)))
Because the evolution of a spinodal microstructure slows with time, we use exponentially increasing time steps to keep the simulation “interesting”. The FiPy user always has direct control over the evolution of their problem.
>>> dexp = -5 >>> elapsed = 0. >>> if __name__ == "__main__": ... duration = 1000. ... else: ... duration = 1000.
>>> while elapsed < duration: ... dt = min(100, exp(dexp)) ... elapsed += dt ... dexp += 0.01 ... eq.solve(phi, dt=dt, solver=DefaultSolver(precon=None)) ... if __name__ == "__main__": ... viewer.plot() ... elif (max(phi.getGlobalValue()) > 0.7) and (min(phi.getGlobalValue()) < 0.3) and elapsed > 10.: ... break
>>> print (max(phi.getGlobalValue()) > 0.7) and (min(phi.getGlobalValue()) < 0.3) True