Benchmark Problem 8: Homogeneous Nucleation¶

See the journal publication entitled “Phase Field Benchmark Problems for Nucleation” for more details about the nucleation benchmark problem. An in-depth analysis focused on this benchmark was also published in the journal publication entitled “On the Effect of Nucleation Undercooling on Phase Transformation Kinetics”. The analysis in this later publication led to several minor changes to the benchmark that are not present in the first journal publication. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

Table of Contents¶

• Overview

• Model Formulation

  • The phase field model in classical nucleation theory

  • The properties of the classical nucleus

  • Avrami plots

• Benchmark 8 Specification

  • Part (a): single seed

  • Part (b): multiple initial seeds

  • Part (c): multiple seeds at random times

• Submission Guidelines

  • Part (a)

  • Part (b)

  • Part (c)

• Results

• References

Overview¶

Phase field modeling of nucleation has a long history and is covered in a number of reviews (Gránásy 2002, Castro 2003, Simmons 2004, Gránásy 2007, Warren 2009, Heo 2014, Gránásy 2019). The problem formulation of crystallization of an ideal pure liquid cooled below its melting point starts with homogeneous nucleation, a process in which the internal fluctuations of the undercooled liquid lead to the formation of crystal-like seeds able to grow to macroscopic sizes. The nucleation can be assisted by the presence of surfaces (container walls, foreign particles, etc.), in which case the process is termed heterogeneous nucleation. We note that homogeneous nucleation is an idealized formulation, and it is unlikely homogeneous nucleation occurs, due to impurities present in experimental apparatus. However, creating benchmark problems for homogeneous nucleation is still needed since the focus in nucleation research often lies on the underlying homogeneous nucleation, as it is the basis for advancing theoretical approaches to the much more complex phenomena governing heterogeneous nucleation. This benchmark (benchmark 8) targets homogeneous nucleation while benchmark 9 targets heterogeneous nucleation.

There are two main modeling approaches to introduce nuclei into a metastable system: the Langevin noise method (Kubo 1966, Gránásy 2002) and the explicit nucleation method (Simmons 2000, Shi 2019) . We focus on the explicit method, where the critical nucleus size and the nucleation energy are determined by the classical nucleation theory. We break down our consideration of homogeneous nucleation into three parts. Part (a) considers the simple case of single-seed homogeneous nucleation. We explore how the particle size influences the evolution of the nucleus when the thermodynamic driving force is close to the critical value where we can observe whether the particle grows, dissolves or remains stationary. In parts (b) and (c) we consider multiple-seed homogeneous nucleation, with nuclei appearing at fixed time $t=0$ or at random times distributed uniformly, respectively, with the latter part illustrating a scenario with constant nucleation rate. The Benchmark Problem probes the differences between transformation kinetics of the two cases and summarizes them using the Avrami plots based on the Johnson-Mehl-Avrami-Kolmogorov (JMAK) theory (Johnson 1939, Avrami 1939, Avrami 1940, Avrami 1941, Kolmogorov 1937).

Model Formulation¶

The phase field model in classical nucleation theory¶

For this benchmark problems we use the simplest possible phase field model with a single non-conserved phase field $\phi$, which describes an isothermal pure substance with one liquid ($\phi=0$) and one solid ($\phi=1$) phase. The free energy of this system is

set $g(\phi)=\phi^2(1-\phi)^2$, a simple double well function with minima at $\phi=0$ and $\phi=1$, $w$ controls the double-well barrier height, $\epsilon^2$ is the gradient energy coefficient, let $p(\phi)=\phi^3(10-15\phi+6\phi^2)$, which ensures that $p(0)=p'(0)=p'(1)=0$ and $p(1)=1$, and $\Delta f$ is the driving force for solidification at the simulation temperature ($\Delta f$ is positive below the melting point). The time evolution of $\phi$ is given by the Allen-Cahn equation (Allen 1979)

where $M$ is the mobility parameter. We will restrict the problem to two dimensions (2D). For a planar interface, one can show that the equilibrium ($\Delta f$=0) solid-liquid interface profile with the interface centered at $x=x_0$ is given by

where $x$ is the perpendicular distance from the interface. We use this expression as the initial condition when we introduce a nucleus. We can also obtain the width $\ell$ and the free energy of the interface, $\gamma$, as

and

Choosing a characteristic length scale $\xi=\ell$ and a characteristic time scale $\tau=1/(Mw)$, we can obtain a nondimensional form of the problem (the nondimensional quantities are denoted by tildes),

and

with $\widetilde{\Delta f}=\Delta f/w$.

The properties of the classical nucleus¶

The explicit nucleation method that we employ introduces nuclei into a metastable system, with their radii and nucleation barrier determined from classical nucleation theory. Classical nucleation theory views crystallite fluctuations appearing in the undercooled liquid as small spherical domains of the bulk crystalline phase bounded by a mathematically sharp solid-liquid interface. For a 2D system, the free energy of a circular solid particle of radius $r$ is

where $\Delta f$ is the nucleation driving force that is used in our phase field model, and $\gamma$ is the free energy of the interface. The free energy of the particle is a balance between the energy cost in forming the solid-liquid interface, and the free energy from the driving force which is released when the crystalline particle forms. Once the rate of change of free energy with respect to particle size becomes negative, the particle can grow. Taking the derivative of $\Delta G(r)$ we get the rate of change of free energy with respect to radius as

By setting $d\Delta G/dr$ to zero, we obtain the critical radius $r^*$ as

The corresponding critical nucleation free energy is

Using the same units as before, the nondimensional forms of these quantities are

and

respectively.

Avrami plots¶

Avrami plots describe how solids transform from one phase to another at constant temperature. In particular, they can describe the kinetics of nucleation. The Avrami plots come from the JMAK theory (Johnson 1939, Avrami 1939, Avrami 1940, Avrami 1941, Kolmogorov 1937), which makes a number of assumptions and simplifications: (i) nucleation occurs randomly and homogeneously over the entire un-transformed portion of the material, (ii) the growth rate is constant and does not depend on the extent of transformation, (iii) the particles have convex shape with the same orientation, and (iv) the size of the system is infinite (both in space and time).

The basis of the Avrami plots is the transformed fraction (solid fraction) $Y(t)$ vs. time. According to the JMAK theory,

where $K$ is a constant depending on the nucleation and growth rates, $n=d+1$ for continuous nucleation and $n=d$ if nucleation happens only at $t=0$, and $d$ is the number of spatial dimensions. If we plot $\log(-\log(1-Y))$ vs. $\log(t)$, then for the JMAK kinetics (Eq. \ref{eqn:Xt}) we get a straight line with slope $n$.

With this straight line, a slope value can be obtained by selecting a fitting range from the values in the range of $-2 < \log(-\log(1-Y)) < 0$ and using linear regression. This range corresponds to the volume fractions $0.02 < Y < 0.9$ which represents a major portion of the transformation. This fitting range enables a consistent fitting range to be applied in all parts which involve Avrami plots (Mancias 2022).

Benchmark 8 Specification¶

From here on, we will use only the nondimensional forms of the phase field (Eqs. \ref{eqn:tFphi}, \ref{eqn:dphidtt}, and \ref{eqn:phir}) and nucleation (Eqs. \ref{eqn:trc} and \ref{eqn:tdeltaGc}) equations, but we will drop the tildes.

We consider 2D simulation domains centered at $x=y=0$. Periodic boundary conditions are applied in order to emulate JMAK theory assumptions of infinite domain and spatiotemporal homogeneity of nucleation and growth.

Part (a): single seed¶

In this problem we examine the morphology change of the nucleus for different initial radii in a domain of $100 \times 100$ units. The driving force is set to $\Delta f = \sqrt{2}/30$, which corresponds to a critical radius of $r^*=5$ (Eq. \ref{eqn:trc}). Next, we place a circular seed at the center of the domain. Incorporating a nucleus into a diffuse interface phase field model leads to a small offset from a classical, sharp-interface model of a nucleus. To account for this, we incorporate a diffuse-interface seed using the profile $\phi(r)$ given by modifying Eq. \ref{eqn:phir} to 2D where $r=\sqrt{x^2+y^2}$, with a radius $r_0=r^*$. This seed is the diffuse-interface approximation of the classical sharp interface nucleus corresponding to the given $\Delta f$, and therefore it should be fairly close to (an unstable) equilibrium. Part (a) is then defined by the Allen-Cahn equation

for the phase field $\phi$, and the three different initial conditions: in the first case with the seed radius $r_0$ corresponding exactly to the critical one (critical nucleus), and in the other two cases, slightly below and above the critical radius (subcritical and supercritical nuclei):

where $r_0=\{r^*, 0.99r^*, 1.01r^*\}$. Figure 1 shows the computational domain with an initial seed of radius $r=r^*$. The time evolution of the system is then followed for times up to $t = 200$ units, and the solid fraction $Y$ and the total free energy $F$ are plotted as functions of time. Finally, the problem also includes a convergence test with respect to mesh (spatial resolution). The closer the initial radius $r$ is to the critical radius $r^*$, the more sensitive the numerical integration will be to round-off errors. Therefore, a convergence test is included for $r_0=1.01r^*$ using successive runs each halving the average spatial mesh size $\Delta x$ between consecutive runs. At least 3 runs are required to estimate the rate of convergence based on the $L^2$ error with the finest mesh being used as a gold standard. The simulations should be run to $t=200$ and the final phase field (at $t=200$) for each successive run should be submitted.

To submit results to the PFHub website, see the submission guidelines below.

Part (b): multiple initial seeds¶

The second part of the homogeneous nucleation benchmark problem focuses on the kinetics of nucleation using Avrami plots and compares them to JMAK theory. To obtain reasonable statistics for these two parts, the simulation volume needs to be larger so it can encompass a larger number of smaller nuclei. We therefore increase the domain size to $500 \times 500$ units of length, and the driving force to $\Delta f = 1/(3 \sqrt{2})$, which corresponds to a critical radius of $r^*=1$. Random initial positions $\mathbf{r}_i$ of 25 supercritical seeds $i, i=1,\ldots,25$ are generated with $r_0=2.2r^*$ drawn from a uniform distribution on the $500\times500$ domain. Periodic boundary conditions are applied. The distribution of the phase field $\phi$ is the sum of the phase fields $\phi_i$ with profiles

from the different seeds $i$. Overlaps between different seeds are managed by setting $\phi=1$ in all regions where the sums of $\phi_i >1$. After placing initial seeds and adjusting the initial phase field $\phi$ so that $\phi\leq1$ everywhere, the simulation is run up to total time $t=200$, at which time the whole domain is transformed to a solid ($\phi=1$). This is repeated at least 10 times, each time with different random seeds, and the total free energy, time evolution of the solid fraction $Y$, discrete particle count $N$ (the number of disjoint regions with $\phi=1$) and Avrami plots are generated from all simulations; from the Avrami plot the exponent $n$ is estimated and compared to the JMAK theory. The number of simulation repeats you can perform in a “reasonable” time will depend on your code and hardware: 10 is required for statistical relevance, but more is better.

Both Part (b) and Part (c) involve a random number generator. Because different executions will in general generate different random numbers, specific individual solutions will in general not be reproduced. However, statistics such as the average slope of Avrami plots can be meaningfully compared.

To submit results to the PFHub website, see the submission guidelines below.

Part (c): multiple seeds at random times¶

For the third part of the homogeneous benchmark problem, the computational domain is first expanded to $1000 \times1000$. Instead of inserting nuclei at fixed time $t = 0$ as in part two of the problem, in this case 100 random nucleation times $t_i$ are generated, $i=1,\ldots,100$, drawn from a uniform distribution in the interval $t_i \in [0, 600)$ with centers $\mathbf{r}_i$ drawn from a uniform distribution on the $1000\times1000$ domain. The driving force is the same as Part (b): $\Delta f=1/(3 \sqrt{2})$, which corresponds to a critical radius $r^*=1$. The initial radii are also kept the same at $r_0 = 2.2r^*$ and the boundary conditions are periodic.

Again, $\phi$ is set to unity in regions of overlaps of nuclei. The simulation is performed up to times $t=600$, and then repeated for a total of ten (10) simulations for statistical significance. The total free energy, time evolution of the solid fraction $Y$, and discrete particle count $N$ are plotted as functions of time, and an Avrami plot is generated for each of the 10 runs.

To submit results to the PFHub website, see the submission guidelines below.

Submission Guidelines¶

Please follow these upload instructions.

In addition to that specified, further data to upload can include a Youtube video, images of the nuclei at different times, or the entire phase field variable at different times. This data can be uploaded directly to the website or stored at a secondary location and only the link to the data provided. This data is not required, but will help others view your work.

Part (a)¶

Part (a) requires running three simulations with a single seed with different values of $r_0$ using the model formulation defined above. The three required simulations are

1. $r_0 = 0.99 r^*$

2. $r_0 = r^*$

3. $r_0 = 1.01 r^*$

Each simulation should be run to $t=200$ and the following data should be collected.

• The solid fraction in the domain (the integral $\int_V \phi \mathrm{d}V / \int_V \mathrm{d}V$) during the transient (not necessarily at every time step but for ~100 data points).

• The free energy integrated over the whole domain, $\mathcal{F}$, during the transient.

• The phase field, $\phi$, at each grid point in the domain at $t=200$.

One PFHub upload is required for all the simulations. The solid fraction and free energy data should be included in the same file for each simulation. The file name needs to be data_{i}.{ext} depending on the simulation number, i, (1 to 3) and the file extension, ext, either csv, json or tsv. The CSV (JSON or TSV are also fine) files should have a format similar to the example below and with the same column headers.

time,fraction,energy
0.000000000000000000e+00,8.214174263100689974e-03,3.579809478015085311e+00
1.000000000000000056e-01,8.205027223394931180e-03,3.577696619091880859e+00
2.000000000000000111e-01,8.197361418448430651e-03,3.576379039721957032e+00
...

The phase field data file should be named phase_field_{i}.{ext} and have the following column headers.

x,y,phi
0.0,1.0,0.5
0.0,1.5,0.49
...

where x and y is a position in the domain and phi is the value at that position.

Part (b)¶

Part (b) requires running ten simulations with multiple seeds randomly positioned at $t=0$, see the definition above. The simulations are run to $t=200$. The following data should be collected.

• The solid fraction in the domain (the integral $\int_V \phi \mathrm{d}V / \int_V \mathrm{d}V$) during the transient (not necessarily at every time step but for ~100 data points).

• The free energy integrated over the whole domain, $\mathcal{F}$, during the transient.

• The discrete particle count, $N$, during the transient.

One PFHub upload is required with ten data files, but with a total of 30 entries in the “Data Files” section. The solid fraction, free energy and particle count can be included in the same data file for each simulation. The entries should be labeled as solid_fraction_{i}, free_energy_{i} or particle_count_{i} depending on the simulation number (1 to 10). The CSV (JSON or TSV are also fine) files should have a format similar to the following.

time,fraction,energy,particle_count
0.000000000000000000e+00,8.214174263100689974e-03,3.579809478015085311e+00,24
1.000000000000000056e-01,8.205027223394931180e-03,3.577696619091880859e+00,24
2.000000000000000111e-01,8.197361418448430651e-03,3.576379039721957032e+00,23
...

Part (c)¶

Part (c) requires running 10 simulations with seeds randomly appearing between $t=0$ and $t=600$, see the definition above. The simulations are run to $t=600$. The following data should be collected.

• The solid fraction in the domain (the integral $\int_V \phi \mathrm{d}V / \int_V \mathrm{d}V$) during the transient (not necessarily at every time step but for ~100 data points).

• The free energy integrated over the whole domain, $\mathcal{F}$, during the transient.

• The discrete particle count, $N$, during the transient.

One PFHub upload is required with ten data files, but with a total of 30 entries in the “Data Files” section. The solid fraction, free energy and particle count can be included in the same data file for each simulation. The entries should be labeled as solid_fraction_{i}, free_energy_{i} or particle_count_{i} depending on the simulation number (1 to 10). The CSV (JSON or TSV are also fine) files should have a format similar to the following.

time,fraction,energy,particle_count
0.000000000000000000e+00,8.214174263100689974e-03,3.579809478015085311e+00,0
1.000000000000000056e-01,8.205027223394931180e-03,3.577696619091880859e+00,1
2.000000000000000111e-01,8.197361418448430651e-03,3.576379039721957032e+00,1
...

Please use the [upload form]({{ site.baseurl }}/simulations/upload_form/) to upload your results.

Results¶

Results from this benchmark problem are displayed on the [simulation result page]({{ site.baseurl }}/simulations/8a.1) for different codes.

References¶

[^] [^]

[^]

[^]

[^]

[^]

[^]

[^]

[^]

[^]

[^]

[^] [^]

[^] [^]

[^] [^]

[^] [^]

[^] [^]

[^]

[^]

Appendix¶

The following code is given as a resource for those implementing the benchmark to help count the particles in Part (b) and Part (c). One possibility as outlined here is to use scipy.ndimage.label. This returns a tuple with a labeled array of unique features as the first item and a count of features as the second item. In the case of the phase field the values will require thresholding to 0 or 1 before using scipy.ndimage.label.

import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt

def make_particles(shape, n_particles, max_radius, seed=99):
    """Make an array with random particles
    
    Args:
      shape: shape of array
      n_particles: number of particles in arrau
      max_radius: maximum possible radius of particles
      seed: random seed
      
    Returns: numpy array mask with particles
    """
    np.random.seed(seed)
        
    def xy_grid(shape):
        """Generate a grid of x, y values in a single array
        
        The returned array has shape (len(shape),) + shape + (1,). The last axis is
        for the operations with particles.
        """
        linspace = lambda x: np.linspace(0, x - 1, x) + 0.5
        return np.array(np.meshgrid(*map(linspace, shape), indexing="ij"))[..., None]
    
    # first axis is dimensions
    # second and third axes are shape of domain
    # last axis is the number of particles. 
    centers = (np.random.random((len(shape), n_particles)) * np.array(shape)[:, None])[:, None, None]
    radii = (np.random.random(n_particles) * max_radius)[None, None]
    return numpy.any(((xy_grid(shape) - centers)**2).sum(0) < radii**2, axis=-1)

arr = make_particles((100, 200), 20, 10)
plt.imshow(arr)
print('Actual number of particles:', scipy.ndimage.label(arr)[1])