Here are a few tutorials meant for beginners:

Title: Propogating waves of self-assembly in organosilane monolayers

Authors: Jack F Douglas, Kirill Efimenko, Daniel A Fischer, Fredrick R Phelan, and Jan Genzer

Source: PNAS, June 19, 2007, Vol. 104, No. 25, pp. 10324-10329


Wavefronts associated with reaction–diffusion and self-assembly processes are ubiquitous in the natural world. For example, propagating fronts arise in crystallization and diverse other thermodynamic ordering processes, in polymerization fronts involved in cell movement and division, as well as in the competitive social interactions and population dynamics of animals at much larger scales. Although it is often claimed that self-sustaining or autocatalytic front propagation is well described by mean-field “reaction–diffusion” or “phase field” ordering models, it has recently become appreciated from simulations and theoretical arguments that fluctuation effects in lower spatial dimensions can lead to appreciable deviations from the classical mean-field theory (MFT) of this type of front propagation. The present work explores these fluctuation effects in a real physical system. In particular, we consider a high-resolution near-edge x-ray absorption fine structure spectroscopy (NEXAFS) study of the spontaneous frontal self-assembly of organosilane (OS) molecules into self-assembled monolayer (SAM) surface-energy gradients on oxidized silicon wafers. We find that these layers organize from the wafer edge as propagating wavefronts having well defined velocities. In accordance with two-dimensional simulations of this type of front propagation that take fluctuation effects into account, we find that the interfacial widths w(t) of these SAM self-assembly fronts exhibit a power-law broadening in time, w(t) {\approx} t^{\beta}, rather than the constant width predicted by MFT. Moreover, the observed exponent values accord rather well with previous simulation and theoretical estimates. These observations have significant implications for diverse types of ordering fronts that occur under confinement conditions in biological or materials-processing contexts.

Title: A discretization of the n-dimensional Laplacian for a dimension-independent stability limit

Author(s): Anand Kumar

Source: Proceedings: Mathematical, Physical and Engineering Sciences, Vol. 457, No. 2015, (Nov. 8, 2001), pp. 2667-2674.

A discretization of the n-dimensional Laplacian employing all the 3^n points of the finite-difference stencil has been presented. Using this discretization, the stability limit of the heat conduction and the wave equations are found to become dimension-independent \frac{1} {2} and 1, respectively.

Title: Spectral implementation of an adaptive moving mesh method for phase-field equations

Authors: W M Feng, P Yu, S Y Hu, Z K Liu, Q Du and L-Q Chen

Source: Journal of Computational Physics, Vo. 220, Issue 1, 20 December 2006, pp. 498-510.


Phase-field simulations have been extensively applied to modeling microstructure evolution during various materials processes. However, large-scale simulations of three-dimensional (3D) microstructures are still computationally expensive. Among recent efforts to develop advanced numerical algorithms, the semi-implicit Fourier spectral method is found to be particularly efficient for systems involving long-range interactions as it is able to utilize the fast Fourier transforms (FFT) on uniform grids. In this paper, we report our recent progress in making grid points spatially adaptive in the physical domain via a moving mesh strategy, while maintaining a uniform grid in the computational domain for the spectral implementation. This approach not only provides more accurate treatment at the interfaces requiring higher resolution, but also retains the numerical efficiency of the semi-implicit Fourier spectral method. Numerical examples using the new adaptive moving mesh semi-implicit Fourier spectral method are presented for both two and three space dimensional microstructure simulations, and they are compared with those obtained by other methods. By maintaining a similar accuracy, the proposed method is shown to be far more efficient than the existing methods for microstructures with small ratios of interfacial widths to the domain size.

Title: Numerical study of the Cahn-Hilliard equation in one, two and three dimensions

Authors: E V L de Mello and Otton Teixeira da Silveira Filho

Source: Physica A: Statistical and theoretical physics, Vol. 347, 1 March 2005, pp. 429-443.


The Cahn–Hilliard (CH) equation is related with a number of interesting physical phenomena like the spinodal decomposition, phase separation and phase ordering dynamics. On the other hand this equation is very stiff and the difficulty to solve it numerically increases with the dimensionality and therefore, there are several published numerical studies in one dimension (1D), dealing with different approaches, and much fewer in two dimensions (2D). In three dimensions (3D) there are very few publications, usually concentrate in some specific result without the details of the used numerical scheme. We present here a stable and fast conservative finite difference scheme to solve the CH with two improvements: a splitting potential into an implicit and explicit in time part and the use of free boundary conditions. We show that gradient stability is achieved in one, two and three dimensions with large time marching steps than normal methods.

In one of my recent posts I noted the peridynamics formulation for the study of phase transformations and wondered how it might be related to phase field models. Over at iMechanica, in the comments section, there are plenty of discussions on the formulation; specifically, here are the sections contrasting the method with phase field methods:

Peridynamics is, very loosely speaking, a continuum analog of pair-potentials in molecular dynamics … It may be a good way to model long-range effects that are important at small scales. We found that the long-range effects provide some interesting mechanics during nucleation that one doesn’t see in phase-field type models (Section 5 and Appendix B of [Journal of the Mechanics and Physics of
Solids, 54 (2006) 1811–1842, and is available at doi:10.1016/j.jmps.2006.04.001] ) …The dispersion relation is similar to atomistic (discrete) models of solids (Section 12 in [1]). In peridynamics, at short wavelengths, the velocity of the waves goes to zero. In atomistic phonon calculations as well, the velocity goes to zero when the wavelength approaches the lattice spacing. This is in sharp contrast to classical continuum (sharp-interface) and phase-field models. In MD models, one doesn’t put dissipation into the atomic potentials, but one has short waves with small velocity generated as a phase boundary (or other defect) moves. These are the source of dissipation when one examines the motion at macroscopic levels and are lumped as heat, roughly. We find a very similar mechanism in peridynamics, where we have dissipation at macroscopic scales without putting any dissipation in at the microscopic scale. This is again unlike phase-field modeling where one needs to explicitly put in viscosity at the microscopic level to obtain dissipation at the macroscopic scale. All of this seems to make peridynamics a interesting tool to examine dynamics processes.

Also, in Stewart Silling’s original paper formulating peridynamics [1], his introduction and discussion sections give a perspective of the theory.

[1] Silling, S.A., 2000. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys.Solids 48, 175–209.

Finally, here are the relevant sections of the paper itself on the question of phase field versus peridynamics:

In classical continuum theory, these phase transforming materials have been modeled using an energy that has multiple minima, each minimum corresponding to a particular phase or variant of martensite. In a dynamic, or even quasistatic, setting, the usual constitutive information, strain energy density as a function of strain, is insufficient to determine a unique solution. For example, even simple Riemann problems with a single phase or twin boundary in the initial conditions allow a one-parameter family of solutions. Therefore, we require further material information to pick the physically correct solution. Abeyaratne and Knowles, 1990 and Abeyaratne and Knowles, 1991b have proposed that this extra information may be specified in the form of a nucleation criterion and a kinetic relation.

The nucleation criterion determines whether a new phase will nucleate from a single phase. The kinetic relation determines the kinetics or the rules that govern the evolution of the phase boundary. It relates the velocity to a thermodynamic driving force, these being conjugate variables in the dissipation (or entropy) inequality. The driving force is related to Eshelby’s idea of the force acting on a defect (Eshelby, 1956 and Eshelby, 1975). The nucleation criterion and the kinetic relation provide uniqueness and well-posedness to initial-boundary value problems. Physically, they can be thought of as a macroscopic remnant of the lattice level atomic motion from one energy well to another that is lost in the continuum theory. However, a systematic derivation from a microscopic theory as well as experimental confirmation remain a topic of active research.

Another approach to overcome the inability of classical continuum mechanics to model the kinetics of phase transformations is to regularize or augment the theory, notably by adding a strain gradient (capillarity) and viscosity to the constitutive relation. This augmented theory leads to a unique solution for the motion of phase boundaries (Abeyaratne and Knowles, 1991a and Truskinovsky, 1993). Further, Abeyaratne and Knowles (1991a) have shown a correspondence between such methods and the kinetic relation. However, nucleation is incompletely explored in this theory, and computational evidence suggests that it is in fact quite difficult. Further, this theory leads to fourth-order equations which are difficult to deal with computationally: they are stiff and one needs smooth elements in the finite element method (see for example, Kloucek and Luskin, 1994 and Dondl and Zimmer, 2004).

There is a closely related phase-field approach (see for example, Artemev et al., 2001 and Wang et al., 1994) in the infinitesimal strain setting. Here, one uses the transformation strain as an internal variable or order parameter, considers the free energy density as a function of this order parameter and uses linear elasticity to penalize the incompatibility in this internal variable field. This leads to a second-order equation which is computationally attractive. However, the equilibrium and the dynamics can be different from that of the regularized theories described earlier (Bhattacharya, 2003). The connection between this theory and kinetic relations remains unexplored (Killough, 1998 has some discussion on this question), nucleation remains difficult and most studies are quasistatic.

Time for some more reading!

Update: A link to peridynamic based fracture code; via Biswajit.

Title: Conservative nonlinear difference scheme for the Cahn-Hilliard equation (Parts I and II)

Authors: S M Choo and S K Chung (Part I); S M Choo, S K Chung and K I Kim(Part II)

Source: Part I — Computers and Mathematics with applications, Vol. 36, Issue 7, October 1998, pp. 31-39; Part II — Computers and Mathematics with applications, Vol . 39, Issues 1-2, January 2000, pp. 229-243.


Part I:

Numerical solutions for the Cahn-Hilliard equation is considered using the Crank-Nicolson type finite difference method. Existence of the solution for the difference scheme has been shown by Brouwer fixed-point theorem. Stability, convergence and error analysis of the scheme are shown. We also show that the scheme preserves the discrete mass, even though the linearized scheme in [1] is conditionally stable and does not preserve the mass.

Part II:

A nonlinear conservative difference scheme is considered for the two-dimensional Cahn-Hilliard equation. Existence of the solution for the finite difference scheme has been shown and the corresponding stability, convergence, and error estimates are discussed. We also show that the scheme preserves the discrete total mass computationally as well as analytically.