Phase Field

See also: PhaseField class API reference — for the method-by-method interface of the PhaseField C++ class.

The phase-field model is used to simulate processes that can be described by the evolution of domain-level order parameters — solidification, solid-state transformations, grain growth, and related problems. It represents interfaces between sub-domains ${\mathrm{\Omega }}_i\subset \mathrm{\Omega }$ within a larger domain $\mathrm{\Omega }\subset {\mathbb{R}}^d$. These sub-domains, referred to as phases, should not be confused with sub-domains in the context of domain decomposition methods.

Phase fields ${\varphi }_{\alpha }:\mathrm{\Omega }\times \mathbb{R}\to [0,1]$ are auxiliary functions used to distinguish between different phases and to describe their evolution over time. A value of 1 in a phase field indicates the presence (or bulk) of the corresponding phase, while a value of 0 signifies its absence. Transitions between phases occur across diffuse interfaces, represented by smooth variations in the phase field functions ${\varphi }_{\alpha }$.

The phase fields satisfy the constraint:

$$\begin{array}{}\text{(1)}& \sum _{\alpha =1}^{\mathrm{\infty }}{\varphi }_{\alpha }=1\end{array}$$

at every point $x\in \mathrm{\Omega }$ within the domain.

The following sections present a detailed description of the multiphase-field model as implemented in OpenPhase.

The phase-field model in OpenPhase describes the evolution of phase fields ${\varphi }_{\alpha }$ using the following equation:

$$\begin{array}{}\text{(2)}& {\dot{\varphi }}_{\alpha }=\frac{\partial }{\partial t}{\varphi }_{\alpha }(x,t)=-\sum _{\beta =1}^N\frac{{\pi }^2{M}_{\alpha \beta }}{8\eta N}(\frac{\delta F}{\delta {\varphi }_{\alpha }}-\frac{\delta F}{\delta {\varphi }_{\beta }})\end{array}$$

where:

• $M_{\alpha \beta }$ is the interface mobility between phase fields ${\varphi }_{\alpha }$ and ${\varphi }_{\beta }$,
• $F$ is the total free energy, given by:

$$\begin{array}{}\text{(3)}& F={\int }_{\mathrm{\Omega }}(f_{\text{GB}}+\dots )dV\end{array}$$

The number of phase fields contributing at a given position $x$ is defined as:

$$\begin{array}{}\text{(4)}& N(x)=\sum _{\alpha =1}^{\mathrm{\infty }}{\chi }_{T(x)}({\varphi }_{\alpha })\end{array}$$

with the active set of phase fields at $x$:

$$\begin{array}{}\text{(5)}& T(x)=\{\varphi \in C^2(\mathrm{\Omega })|\varphi (x)\ne 0\vee \mathrm{\nabla }\varphi (x)\ne 0\vee {\mathrm{\nabla }}^2\varphi (x)\ne 0\}\end{array}$$

The interfacial free energy density is defined as:

$$\begin{array}{}\text{(6)}& f_{\text{GB}}=\sum _{\alpha \ne \beta }\frac{4{\sigma }_{\alpha \beta }}{\eta }(-\frac{{\eta }^2}{{\pi }^2}\mathrm{\nabla }{\varphi }_{\alpha }\cdot \mathrm{\nabla }{\varphi }_{\beta }+{\varphi }_{\alpha }{\varphi }_{\beta })\end{array}$$

Here:

• $\eta$ is the interface width,
• ${\sigma }_{\alpha \beta }$ is the interfacial energy between phases $\alpha$ and $\beta$.

Defining the interface term:

$$\begin{array}{}\text{(7)}& I_{\alpha }={\mathrm{\nabla }}^2{\varphi }_{\alpha }+\frac{{\pi }^2}{{\eta }^2}{\varphi }_{\alpha }\end{array}$$

the grain boundary contribution to the phase evolution becomes:

$$\begin{array}{}\text{(8)}& {\dot{\varphi }}_{\alpha }^{\text{GB}}=\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N\frac{M_{\alpha \beta }}{N}(\sum _{\begin{array}{c}\gamma =1\\ \gamma \ne \beta \end{array}}^N({\sigma }_{\beta \gamma }-{\sigma }_{\alpha \gamma })I_{\gamma })\end{array}$$

An alternative form:

$$\begin{array}{}\text{(9)}& {\dot{\varphi }}_{\alpha }^{\text{GB}}=\frac{1}{N}\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N{M}_{\alpha \beta }[{\sigma }_{\alpha \beta }(I_{\alpha }-I_{\beta })+\sum _{\begin{array}{c}\gamma =1\\ \gamma \ne \alpha ,\gamma \ne \beta \end{array}}^N({\sigma }_{\beta \gamma }-{\sigma }_{\alpha \gamma })I_{\gamma }]\end{array}$$

Or, in terms of interface fields ${\psi }_{\alpha \beta }^{\text{GB}}$:

$$\begin{array}{}\text{(10)}& {\dot{\varphi }}_{\alpha }^{\text{GB}}=\frac{1}{N}\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N{\dot{\psi }}_{\alpha \beta }^{\text{GB}}\end{array}$$

Here, ${\psi }_{\alpha \beta }^{\text{GB}}$ represents the contribution of the grain boundary interface between phase fields $\alpha$ and $\beta$.

Computational Efficiency

A naive implementation of the phase-field equation would have computational complexity $\mathcal{O}(n^3)$, where $n$ is the number of phase fields. This is impractical for simulations involving many phases.

To address this, active parameter tracking was introduced to reduce computation by discarding insignificant phase field contributions. In models using a double-well potential, the phase-field profile is governed by a $\tanh$ function, so values approach but never exactly reach 0 or 1. Truncating near-zero values in such cases may impact accuracy.

However, for models using a double-obstacle potential (with $\sin$-shaped profiles), we can apply the following:

Remark 1: If both ${\varphi }_k(x)=0$ and ${\mathrm{\nabla }}^2{\varphi }_k(x)=0$, then phase $k$ does not contribute to ${\dot{\psi }}_{ij}^{\text{GB}}$.

This allows safe omission of these terms without affecting results.

Remark 2: It is still possible for ${\dot{\varphi }}_i^n(x)>0$ even when both ${\varphi }_i^n(x)=0$ and ${\mathrm{\nabla }}^2{\varphi }_i^n(x)=0$. This signifies the nucleation of a new phase, which is handled separately in the nucleation section.

Incorporating this into the evolution equation, we obtain:

$$\begin{array}{}\text{(11)}& {\dot{\varphi }}_{\alpha }={\chi }_{\alpha }\frac{1}{N}\left[\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N{\chi }_{\beta }{M}_{\alpha \beta }({\sigma }_{\alpha \beta }(I_{\alpha }-I_{\beta })+\sum _{\begin{array}{c}\gamma =1\\ \gamma \ne \alpha ,\gamma \ne \beta \end{array}}^N({\sigma }_{\beta \gamma }-{\sigma }_{\alpha \gamma })I_{\gamma })\right]\end{array}$$

Here, the indicator function ${\chi }_{\alpha }$ determines whether the phase-field ${\varphi }_{\alpha }$ is active:

$$\begin{array}{}\text{(12)}& {\chi }_{\alpha }=\{\begin{array}{ll}1& \text{if }{\varphi }_{\alpha }\ne 0\vee {\mathrm{\nabla }}^2{\varphi }_{\alpha }\ne 0\\ 0& \text{otherwise}\end{array}\end{array}$$

This optimization allows the solver to skip unnecessary computations while preserving accuracy.

Normal Grain Growth Example

A basic example for the multiphase-field method is the simulation of grain growth, it only considers the interface energies and not other contributions to the driving force. Minimizing the interface energy results in large grains to grow and small grains to shrink. Here is a basic example code for a grain growth simulation

#include "Settings.h"
#include "RunTimeControl.h"
#include "DoubleObstacle.h"
#include "PhaseField.h"
#include "Initializations.h"
#include "BoundaryConditions.h"
#include "InterfaceProperties.h"
#include "Tools/TimeInfo.h"
#include "Tools/MicrostructureAnalysis.h"

using namespace std;
using namespace openphase;

int main(int argc, char *argv[])
{
	std::string InputFile;
    if(argc > 1)
    {
        InputFile = argv[1];
    }
    else
    {
        std::cerr << "No Inputfile provided, trying to use default ProjectInput.opi" << std::endl;
        InputFile = "ProjectInput.opi";
    }
    Settings                    OPSettings;
    OPSettings.ReadInput(InputFile);

    RunTimeControl              RTC(OPSettings, InputFile);
    PhaseField                  Phi(OPSettings, InputFile);
    DoubleObstacle              DO(OPSettings, InputFile);
    InterfaceProperties         IP(OPSettings, InputFile);
    BoundaryConditions          BC(OPSettings, InputFile);

    //generating initial grain structure using Voronoi algorithm
    int number_of_grains = 200;
    size_t GrainsPhase = 0;
    Initializations::VoronoiTessellation(Phi, BC, number_of_grains, GrainsPhase);

    std::cout << "Entering the Time Loop!!!" << std::endl;
    for(RTC.tStep = RTC.tStart; RTC.tStep <= RTC.nSteps; RTC.IncrementTimeStep())
    {
        IP.Set(Phi, BC);
        DO.CalculatePhaseFieldIncrements(Phi, IP);
        Phi.NormalizeIncrements(BC, RTC.dt);
        Phi.MergeIncrements(BC, RTC.dt);

        /// Output to VTK file
        if (RTC.WriteVTK())
        {
            Phi.WriteVTK(OPSettings, RTC.tStep);
            MicrostructureAnalysis::WriteGrainsStatistics(Phi,RTC.tStep);
        }
        /// Output raw data
        if (RTC.WriteRawData())
        {
            Phi.Write(OPSettings, RTC.tStep);
        }
        /// Output to screen
        if (RTC.WriteToScreen())
        {
            double I_En = DO.AverageEnergyDensity(Phi, IP);
            std::string message = ConsoleOutput::GetStandard("Interface energy density", I_En);
            ConsoleOutput::WriteTimeStep(RTC, message);
        }
    }
    return EXIT_SUCCESS;
}

This example shows how to setup a simple simulation using the OpenPhase library. The path to input file (*.opi) for this simulation can be given as an argument when calling the executable of this example. The input file should contain the following information:

@RunTimeControl

$SimTtl   Simulation Title                          : Normal grain growth
$nSteps   Number of Time Steps                      : 2000
$FTime    Output to disk every (tSteps)             : 100
$STime    Output to screen every (tSteps)           : 100
$LUnits   Units of length                           : m
$TUnits   Units of time                             : s
$MUnits   Units of mass                             : kg
$EUnits   Energy units                              : J
$dt       Initial Time Step                         : 1e-5
$nOMP     Number of OpenMP Threads                  : 1
$Restrt   Restart switch (Yes/No)                   : No
$tStart   Restart at time step                      : 0
$tRstrt   Restart output every (tSteps)             : 10000

@GridParameters

$Nx       System Size in X Direction                : 101
$Ny       System Size in Y Direction                : 101
$Nz       System Size in Z Direction                : 101
$dx       Grid Spacing                              : 1e-6
$IWidth   Interface Width (in grid points)          : 5.0

@Settings

$Phase_0  Name of Phase 0                           : Phase1

@InterfaceProperties

$EnergyModel_0_0    Interface energy model          : ISO
$Sigma_0_0          Interface energy                : 0.24

$MobilityModel_0_0  Interface energy model          : ISO
$Mu_0_0             Interface mobility              : 1.0e-7

@BoundaryConditions

$BC0X   X axis beginning boundary condition         : Periodic
$BCNX   X axis far end boundary condition           : Periodic

$BC0Y   Y axis beginning boundary condition         : Periodic
$BCNY   Y axis far end boundary condition           : Periodic

$BC0Z   Z axis beginning boundary condition         : Periodic
$BCNZ   Z axis far end boundary condition           : Periodic

In this RunTimeControl controls anything related to time and time steps, such as the time step dt, the total number of timesteps nSteps as well as the frequency of output to files FTime and to screen STime. The functions RTC.WriteVTK() and RTC.WriteRawData() will return true if and only if

$$\begin{array}{}\text{(13)}& 0\equiv RTC.tStep(\mathrm{mod}FTime),\end{array}$$

same for RTC.WriteToScreen() and STime. Gridparameters is a subclass of Settings and the input is automatically read when the input of Settings is read, therefore Gridparameters does not appear in the c++-code explicitly, but is used by all other classes that contain storages. This is the reason why Settings needs call ReadInput before all other classes are constructed as information such as domain size needs to passed to each class so that the storages of the appropriate size can be allocated in the memory. Settings itself only contains the name of the one phase in this simulation, in general it contains the names of all phases and elements. InterfaceProperties is used to define the interface energy and mobility and their anisotropy, for the simulation of normal grain growth the interface energy and mobility are both isotropic. For this simulation we also use periodic boundary conditions in all directions. Instead of

RunTimeControl              RTC(OPSettings, InputFile);

it also possible to create an object with the default constructor and then call the functions Initialize(Settings& settings) and ReadInput(std::string InputFile) like this

RunTimeControl              RTC;
RTC.Initialize(OPSettings);
RTC.ReadInput(InputFile);

this can be done for all classes in the OpenPhase library that have an Initialize function with Settings as an argument. The initializations class offers a variety of static functions that produce an initialize arrangement of phases for the simulation, in this we use a Voronoi Tesselation. In the timeloop

IP.Set(Phi, BC);
DO.CalculatePhaseFieldIncrements(Phi, IP);
Phi.NormalizeIncrements(BC, RTC.dt);
Phi.MergeIncrements(BC, RTC.dt);

the function IP.Set calculates the interface energy and mobility in each interface point. As this is a simulation with isotropic interface energy and mobility the calculation is trivial. In further simulation anisotropy and temperature dependence of the interface properties will be discussed. DO.CalculatePhaseFieldIncrements calculates the increments $dot{\varphi }_{\alpha }$, depending on the interface properties. Phi.NormalizeIncrements is an auxiliary function that in itself is not present in the phasefield equation, it limits the increment to ${\dot{\varphi }}_{\alpha }$ to make sure that increment ${\varphi }_{\alpha }+dt{\dot{\varphi }}_{\alpha }^n\in [0,1]$. Finally Phi.MergeIncrements applies the increment to phasefields and finalizes the timestep. The other function are just used to output the results of the simulation.

See also: examples

In OpenPhase-main/examples/:

• NormalGG — the canonical normal-grain-growth simulation; the worked example on this page is built from it.
• FacetedGG — grain growth with faceted anisotropy.
• Pearlite, PrecipitationNiTi, Superalloys — solid-state phase transformations.
• SolidificationFeC etc. — solidification (combine PhaseField with Composition and EquilibriumPartitionBinary).
• HeatDiffusion, LatentHeat — phase-field with a coupled thermal field.

Driving Forces

The free energy in the phase-field model can be extended by including additional energy contributions, such as an additional driving force term $f_{\text{AD}}$. The total free energy then becomes:

$$\begin{array}{}\text{(14)}& F={\int }_{\mathrm{\Omega }}(f_{\text{GB}}+f_{\text{AD}})dV\end{array}$$

Using a variational derivation, the phase-field evolution equation is extended to:

$$\begin{array}{}\text{(15)}& {\dot{\varphi }}_{\alpha }=\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N\frac{M_{\alpha \beta }}{N}[\sum _{\begin{array}{c}\gamma =1\\ \gamma \ne \beta \end{array}}^N({\sigma }_{\beta \gamma }-{\sigma }_{\alpha \gamma })I_{\gamma }+\frac{{\pi }^2}{8\eta }\mathrm{\Delta }{g}_{\alpha \beta }]\end{array}$$

with the driving force:

$$\begin{array}{}\text{(16)}& \mathrm{\Delta }{g}_{\alpha \beta }=\frac{\partial f_{\text{AD}}}{\partial {\varphi }_{\beta }}-\frac{\partial f_{\text{AD}}}{\partial {\varphi }_{\alpha }}\end{array}$$

A key strength of the phase-field model in OpenPhase is its rigorous treatment of interface properties. A fundamental assumption of the model is that a steady phase-field contour along the interface normal is maintained throughout the simulation. This ensures a constant interface width, which is essential for correctly modeling interface properties and kinetics.

Depending on the origin of the driving force $\mathrm{\Delta }g$ and the context of the simulation, the driving force may vary across the interface. This variation can distort the phase-field profile and negatively impact interface kinetics. To preserve a correct profile, the driving force must remain constant along the normal direction of the interface and vary only along the interface.

To enforce this, $\mathrm{\Delta }{g}_{\alpha \beta }$ is replaced with $h^{\prime }\mathrm{\Delta }{\overline{g}}_{\alpha \beta }$, where $\mathrm{\Delta }{\overline{g}}_{\alpha \beta }$ is the average driving force constant along the normal direction. The prefactor $h^{\prime }$ is:

$$\begin{array}{}\text{(17)}& h^{\prime }=\frac{8}{\pi }\sqrt{{\varphi }_{\alpha }{\varphi }_{\beta }}\end{array}$$

This form guarantees that the driving force contribution is independent of spatial discretization, similar to the curvature term. Integrating the traveling wave solution of the phase-field equation yields the final expression:

$$\begin{array}{}\text{(18)}& {\dot{\varphi }}_{\alpha }=\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N\frac{M_{\alpha \beta }}{N}[\sum _{\begin{array}{c}\gamma =1\\ \gamma \ne \beta \end{array}}^N({\sigma }_{\beta \gamma }-{\sigma }_{\alpha \gamma })I_{\gamma }+\frac{\pi }{\eta }\sqrt{{\varphi }_{\alpha }{\varphi }_{\beta }}\mathrm{\Delta }{\overline{g}}_{\alpha \beta }]\end{array}$$

This can be decomposed into:

$$\begin{array}{}\text{(19)}& {\dot{\varphi }}_{\alpha }=\frac{1}{N}\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N{\dot{\psi }}_{\alpha \beta }^{\text{GB}}+\frac{1}{N}\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N{M}_{\alpha \beta }\frac{\pi }{\eta }\sqrt{{\varphi }_{\alpha }{\varphi }_{\beta }}\mathrm{\Delta }{\overline{g}}_{\alpha \beta }\end{array}$$

or:

$$\begin{array}{}\text{(20)}& {\dot{\varphi }}_{\alpha }=\frac{1}{N}\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N({\dot{\psi }}_{\alpha \beta }^{\text{GB}}+{\dot{\psi }}_{\alpha \beta }^{\text{AD}})\end{array}$$

and finally:

$$\begin{array}{}\text{(21)}& {\dot{\varphi }}_{\alpha }={\dot{\varphi }}_{\alpha }^{\text{GB}}+{\dot{\varphi }}_{\alpha }^{\text{AD}}\end{array}$$