Thermodynamics and Diffusion in Multiphase-Field Modeling (OpenPhase)

This chapter introduces the thermodynamic and diffusion treatment in OpenPhase for multiphase-field simulations.

---

1. Governing Equation for Phase Evolution

In the multiphase-field method, each phase ${\varphi }_{\alpha }$ evolves according to

$$\begin{array}{}\text{(1)}& {\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}$$

Here:

• $M_{\alpha \beta }$: interface mobility,
• ${\sigma }_{\alpha \beta }$: interface energy,
• $I_{\gamma }$: geometrical function of interface orientation,
• $\mathrm{\Delta }{\overline{g}}_{\alpha \beta }$: thermodynamic driving force (see Section 2).

This formulation couples interfacial thermodynamics with kinetics: the interfacial term (surface energy differences) and the thermodynamic driving force (undercooling) together determine how the phase fields evolve.

Where does $\mathrm{\Delta }{\overline{g}}_{\alpha \beta }$ come from?
It stems from compositional undercooling relative to equilibrium, derived in Section 2, and is consistent with the partitioning rules discussed in Sections 3–5.

---

2. Thermodynamic Driving Force

The thermodynamic driving force between two phases $\alpha$ and $\beta$ is given by the undercooling:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{g}_{\alpha \beta }=\mathrm{\Delta }{S}_{\alpha \beta }\mathrm{\Delta }{T}_{\alpha \beta },\mathrm{\Delta }{T}_{\alpha \beta }=m_{\alpha \beta }\mathrm{\Delta }{c}_{\alpha \beta }=m_{\alpha \beta }(c_{\alpha }-c_{\alpha }^{eq}),\end{array}$$

with:

• $\mathrm{\Delta }{S}_{\alpha \beta }$: entropy difference,
• $m_{\alpha \beta }$: slope of the liquidus/solidus line,
• $\mathrm{\Delta }{c}_{\alpha \beta }$: deviation from equilibrium concentration.

Thus, the driving force is controlled by the undercooling relative to equilibrium.
The equal-undercooling condition is the cornerstone of the partitioning relations used in Section 3 and generalized to multiphase junctions in Section 4.

Engineering note. The quantity $m_{\alpha \beta }$ is direction-specific (ordered pair) and its ratio forms the partitioning coefficient $K_{\alpha \beta }=m_{\alpha \beta }/m_{\beta \alpha }$ used explicitly in Sections 4.1–4.2.

---

3. Partitioning at Interfaces

At an $\alpha$–$\beta$ interface, solute concentrations in the phases are determined by:

1. Mass conservation:

$$\begin{array}{}\text{(3)}& {\varphi }_{\alpha }{c}_{\alpha }+{\varphi }_{\beta }{c}_{\beta }=c_{Total},\end{array}$$

2. Partitioning condition (equal undercooling):

$$\begin{array}{}\text{(4)}& m_{\alpha \beta }(c_{\alpha }-{\overline{c}}_{\alpha \beta })=m_{\beta \alpha }(c_{\beta }-{\overline{c}}_{\alpha \beta }),\end{array}$$

with ${\overline{c}}_{\alpha \beta }$ the intersection point of the tangent slopes.

This system of linear equations yields unique solutions for $c_{\alpha }$ and $c_{\beta }$.
These pairwise relations are consistent with the driving-force definition in Section 2.

Tip. In binaries, this reproduces the classical lever-rule geometry on the phase diagram.
In practice, OpenPhase resolves pairwise partitioning exactly on two-phase interfaces and uses the schemes of Section 4 at multi-junctions.

---

4. Multiphase Junctions

When more than two phases meet, enforcing all pairwise partitioning conditions from Section 3 makes the system overdetermined.
OpenPhase provides three schemes to resolve this:

4.1 Two-Step Scheme

1. Compute the total deviation from equilibrium:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{c}_{Total}=c_{Total}-\sum _{\alpha }{\varphi }_{\alpha }{c}_{\alpha }^{eq}.\end{array}$$

2. Redistribute solute according to:

$$\begin{array}{}\text{(6)}& c_{\alpha }=c_{\alpha }^{eq}+\mathrm{\Delta }{c}_{Total}{\left[\sum _{\beta \ne \alpha }{\varphi }_{\beta }\right]}^{-1}\sum _{\beta \ne \alpha }\frac{{\varphi }_{\beta }}{{\varphi }_{\alpha }+K_{\beta \alpha }{\varphi }_{\beta }}.\end{array}$$

where $K_{\alpha \beta }=m_{\alpha \beta }/m_{\beta \alpha }$ (see Section 2).

Rationale (engineering intuition). Two-Step is a rule-of-mixtures redistribution of the global deviation $\mathrm{\Delta }{c}_{Total}$, weighted by neighboring phase fractions and partition coefficients.

4.2 Eiken Scheme

Based on ^[1]:

$$\begin{array}{}\text{(7)}& c_{\alpha }=\frac{c_{Total}-\sum _{\beta }{\varphi }_{\beta }(c_{\beta }^{eq}-c_{\alpha }^{eq}{K}_{\alpha \beta })}{\sum _{\beta }{\varphi }_{\beta }{K}_{\alpha \beta }}.\end{array}$$

Rationale. Eiken embeds the slope consistency more directly via $K_{\alpha \beta }$ and generally yields better thermodynamic fidelity than Two-Step, especially when slopes are highly asymmetric.

4.3 Optimization Scheme

Does not rely on pairwise partitioning, but minimizes deviation from equilibrium:

$$\begin{array}{}\text{(8)}& min\frac{1}{2}\sum _{\alpha }\sum _{\beta \ne \alpha }{m}_{\alpha \beta }(c_{\alpha }-c_{\alpha }^{eq}{)}^2,\end{array}$$

subject to global mass conservation:

$$\begin{array}{}\text{(9)}& c_{Total}=\sum _{\alpha }{\varphi }_{\alpha }{c}_{\alpha }.\end{array}$$

Rationale. This enforces a global, consistent solution across all phases. It is the most robust in complex multi-junctions and multicomponent systems.

---

5. Independent Partitioning of Components

In multicomponent alloys, each component partitions independently:

• For component $i$ at an $\alpha$–$\beta$ interface:

$$\begin{array}{}\text{(10)}& {\varphi }_{\alpha }{c}_{\alpha }^i+{\varphi }_{\beta }{c}_{\beta }^i=c_{Total}^i,\end{array}$$

$$\begin{array}{}\text{(11)}& m_{\alpha \beta }^i(c_{\alpha }^i-{\overline{c}}_{\alpha \beta }^i)=m_{\beta \alpha }^i(c_{\beta }^i-{\overline{c}}_{\alpha \beta }^i).\end{array}$$

• Partition coefficients $K_{\alpha \beta }^i=m_{\alpha \beta }^i/m_{\beta \alpha }^i$ are component-specific.

This modular treatment means that solutes A, B, C, … are redistributed independently, with coupling only via the underlying thermodynamics (e.g., CALPHAD databases).
This independence is consistent with the diffusion treatment per component in Section 7 and the reference-element constraint in Section 6.

---

6. Reference Element in OpenPhase

OpenPhase uses a reference element (typically the solvent).

• The concentrations satisfy:

$$\begin{array}{}\text{(12)}& c_{ref}+\sum _i{c}_i=1,\end{array}$$

where $c_{ref}$ is the concentration of the reference element and $c_i$ are the solute components.

• The reference element is not explicitly diffused; it is defined implicitly from the other components.

• This ensures:

  • Mass conservation,
  • Reduction of numerical effort (one fewer diffusion equation),
  • A natural treatment of solvent-based alloys.

This constraint integrates seamlessly with independent partitioning (Section 5) and diffusion (Section 7).

---

7. Diffusion of Phase Concentrations (Fick’s Law)

Once partitioned, solute evolves by diffusion:

• Fick’s first law in phase $\alpha$:

$$\begin{array}{}\text{(13)}& {\mathbf{J}}_{\alpha }^i=-D_{\alpha }^i\mathrm{\nabla }{c}_{\alpha }^i,\end{array}$$

• Fick’s second law:

$$\begin{array}{}\text{(14)}& \frac{\partial c_{\alpha }^i}{\partial t}=\mathrm{\nabla }\cdot (D_{\alpha }^i\mathrm{\nabla }{c}_{\alpha }^i).\end{array}$$

• Mixture concentration:

$$\begin{array}{}\text{(15)}& c^i(\mathbf{x},t)=\sum _{\alpha }{\varphi }_{\alpha }(\mathbf{x},t)c_{\alpha }^i(\mathbf{x},t).\end{array}$$

This ensures solute conservation during diffusion and phase transformation.
During interface motion, the partitioning schemes (Section 4) determine instantaneous interfacial split, while Fickian diffusion relaxes gradients within phases.

---

8. Antitrapping Flux

Why it is needed

The diffuse interface used in simulations is much thicker than in reality.
When the interface moves, this can cause artificial solute trapping—a numerical artifact that would otherwise corrupt the physically correct partitioning from Sections 3–5.

Correction

An antitrapping flux is added:

$$\begin{array}{}\text{(16)}& {\mathbf{J}}^i=-\sum _{\alpha }{\varphi }_{\alpha }{D}_{\alpha }^i\mathrm{\nabla }{c}_{\alpha }^i+{\mathbf{J}}_{at}^i,\end{array}$$

with

$$\begin{array}{}\text{(17)}& {\mathbf{J}}_{at}^i=-aW\frac{\partial {\varphi }_{\alpha }}{\partial t}\hat{\mathbf{n}}(c_{\alpha }^i-c_{\beta }^i).\end{array}$$

• $a$: numerical constant,
• $W$: interface width,
• $\hat{\mathbf{n}}$: interface normal,
• $(c_{\alpha }^i-c_{\beta }^i)$: solute contrast.

This term counteracts non-physical solute storage in the diffuse interface and vanishes in the sharp-interface limit. It works in tandem with Fickian diffusion (Section 7) to maintain correct solute transport.

---

9. Engineering Takeaway

• OpenPhase enforces thermodynamic consistency via partitioning schemes (Sections 3–4).
• Each component is partitioned independently (Section 5), with one reference element determined implicitly (Section 6).
• Diffusion transports solute within each phase (Section 7), and antitrapping flux corrects artifacts from diffuse interfaces (Section 8).
• Engineers can select the partitioning scheme depending on the balance between accuracy and computational cost (Section 9).

Together, these features allow OpenPhase to realistically simulate solidification, precipitation, grain growth, and phase transformations in engineering alloys.

---

10. Implementation in OpenPhase

To make practical use of the concepts described above, OpenPhase provides a set of classes and interfaces that handle thermodynamic data, diffusion, and partitioning.

If expanding the User Driving Force example, the following additional objects are needed:

    Composition                      Cx(OPSettings,InputFile);
    ThermodynamicPropertiesEQP       TP(OPSettings,InputFile);
    DiffusionProperties              DP(OPSettings,InputFile);
    ThermodynamicInterfaceLPD        TQ(OPSettings,InputFile);
    ThermokineticInterface           TK(OPSettings,InputFile);
    ChemicalDiffusion                DF(OPSettings,InputFile);
    EquilibriumPartitioning          EQP(OPSettings,InputFile);

• Composition is a storage class that stores total and phase compositions.
• ThermodynamicPropertiesEQP and DiffusionProperties store the relevant data for thermodynamics such as slopes, equilibrium concentrations, and diffusion coefficients.
• ThermodynamicInterfaceLPD reads thermodynamic data from a linearized phase diagram stored in an OPID file (see OPID Format).
• ThermokineticInterface reads kinetic data in Arrhenius form from an OPID file.
• ChemicalDiffusion handles diffusion and the antitrapping flux (see Section 8).
• EquilibriumPartitioning handles the calculation of phase concentrations and the thermodynamic driving force (see Section 2 and Section 3).

---

10.1 Timeloop Expansion

The timeloop from the UserDrivingForce example can be expanded as follows:

        IP.Set(Phi, Tx, BC);
        DO.CalculatePhaseFieldIncrements(Phi, IP, dG);

        TP.CheckRange(Phi,Cx,Tx);

        TQ.SetEquilibriumData(Phi,Cx,Tx,TP);
        EQP.CalculatePhaseConcentrations(Phi,Cx,Tx,TP,BC);

        TK.SetDiffusionData(Phi,Cx,Tx,DP);

        TP.CalculateDrivingForce(Phi,IP,Cx,Tx,dG);
        dG.Average(Phi, BC);
        dG.MergePhaseFieldIncrements(Phi, IP);

        Phi.NormalizeIncrements(BC, RTC.dt);

        DF.Solve(Phi,Cx,Tx,DP,TP,BC,RTC.dt);

        Phi.MergeIncrements(BC, RTC.dt);

---

10.2 Step-by-Step Explanation

• TP.CheckRange(Phi,Cx,Tx)
  Checks if thermodynamic data needs to be updated because temperature or composition changed significantly.

  • Thresholds are defined in the OPI file:
    • $MDev for mole fraction changes,
    • $TDev for temperature changes,
    • under the section @ThermodynamicPropertiesEQP.

• TQ.SetEquilibriumData(Phi,Cx,Tx,TP)
  Sets the equilibrium phase concentrations for the current total composition and temperature in each grid point.

• EQP.CalculatePhaseConcentrations(Phi,Cx,Tx,TP,BC)
  Calculates the phase concentrations using one of the partitioning schemes (see Section 4).

• TK.SetDiffusionData(Phi,Cx,Tx,DP)
  Obtains diffusion coefficients from the thermokinetic interface.

• TP.CalculateDrivingForce(Phi,IP,Cx,Tx,dG)
  Calculates the thermodynamic contribution to the driving force using the governing equation introduced in Section 1.

• Phi.NormalizeIncrements(BC, RTC.dt)
  Normalizes the phase-field increments to ensure stability and consistency.

• DF.Solve(Phi,Cx,Tx,DP,TP,BC,RTC.dt)
  Solves the diffusion equation, including the antitrapping flux correction (see Section 8) if activated in the OPI file.

• Phi.MergeIncrements(BC, RTC.dt)
  Final update of the phase fields with the computed increments.

---

See also

• Thermodynamic Phase — per-phase thermodynamic data container.
• Sublattice Model — CALPHAD-style Gibbs-energy formulation.
• Grand Potential Solver — the primary consumer of the thermodynamic state.
• Equilibrium Partition (Binary) — the binary-alloy diffusion solver whose partitioning theory is covered on this page.
• Composition — shared component / concentration storage.

References

[1] J. Eiken et al. (2006). Multiphase-field approach for multicomponent alloys with extrapolation scheme for numerical application. Physical Review E, 73, 066122. DOI: 10.1103/PhysRevE.73.066122 ↩