# High-fidelity level-set modeling of diffusive solid-state phase
transformations for polycrystalline materials

Nitish Chandrappa
Mines Paris, PSL University
Centre for material forming (CEMEF), UMR CNRS
06904 Sophia Antipolis, France
<EMAIL_ADDRESS>
Marc Bernacki
Mines Paris, PSL University
Centre for material forming (CEMEF), UMR CNRS
06904 Sophia Antipolis, France
<EMAIL_ADDRESS>

###### Abstract

The formation of microstructures in metallic alloys during hot metal forming
involves simultaneous metallurgical complex phenomena. Traditional high-
fidelity numerical frameworks used on the polycrystalline scale tend to focus
on single-phase microstructures or isolate phase transformations from grain
boundary migration mechanisms. The level-set method is highlighted as
effective in proposing a global framework for modeling multiphase
polycrystalline materials and diffusive solid-state phase transformations.
This framework includes novel techniques for efficient large-scale
microstructural representation, strong coupling with ThermoCalc software for
real-time thermodynamic data, application for ternary alloys and beyond by
taking solute drag aspects, and the use of advanced nucleation models.
Numerous applications are then illustrated.

_K_ eywords Full-field method, level-set, microstructural evolution, diffusive
phase transformation, austenite decomposition

##  1 Introduction

The microstructural evolution in metallic alloys in the context of hot metal
forming typically involves a complex interplay of multiple phenomena that
occur simultaneously. Industrial metallurgical products are often subjected to
complex thermomechanical treatments or processing conditions such as high
plastic deformation at elevated temperatures. Under these conditions, solid-
state phase transformations within alloyed materials become inevitable, along
with other concurrent phenomena. Currently, most state-of-the-art numerical
predictions primarily focus on single-phase microstructural evolution or
exclusively model phase transformation, often neglecting the consideration of
other coexisting phenomena. Although delving into individual phenomena
provides valuable insights, focusing solely on singular aspects might
inadvertently lead to oversights and limitations in comprehending the holistic
behavior of these materials. Such a narrow focus could lead to overlooking
crucial interactions that may influence the alloy’s mechanical strength,
thermal stability, or its susceptibility. Neglecting the concomitant nature of
these transformation aspects could limit our ability to predict and optimize
material performance under diverse operating conditions.

Advances in computational resources have paved the way for more intricate
models (such as atomistic [1, 2, 3] and full-field mesoscopic models [4])
based on a complete and explicit representation of the microstructure. These
models enable to closely follow the topological evolution of the
microstructure during a transformation and have the potential to capture
complex evolution aspects, and hence their predictive capabilities are much
wider. Full-field models (FFM) provide a good equilibrium between the high-
level description of the microstructure and the computational demands. Among
the various approaches within the scope of FFM, each numerical method has its
unique scope of application, with some better suited to specific scenarios
than others. For modeling recrystallization and grain growth under high
plastic deformation, the level-set method [5, 6, 7, 8] stands out as a robust
choice and has been effectively utilized. At the mesoscopic scale, Phase-field
methods [9, 10, 11, 12, 13, 14] have been well established as the popular
choice of method to simulate diffusive solid-state phase transformation
(DSSPT) due to its thermodynamic consistency. However, the potential of level-
set methods remains largely unexplored in the realm of DSSPT. The developments
of Bzowski et al. [15] is one of the only few works based on level-set method
for DSSPT. In these works, one level-set function per grain is considered, and
to model solute diffusion, the diffusion equation is resolved in the parent
phase while that in the product phase his considered as negligible.

The authors thus proposed a global level-set framework in [16] aimed at
thoroughly exploring the capabilities of using the level-set method to
simulate DSSPT in a finite element (FE) context. Within this framework, the
grain coloration technique, proposed in [17, 18] was used to limit the number
of level-set functions required to describe a large-scale microstructure. The
diffusion aspects are modeled in both the parent and the product phase.
Illustrations for phase transformation in simple binary alloyed
microstructures were demonstrated. The objective of this work is to generalize
the previously proposed numerical framework for ternary alloys and beyond by
taking solute drag aspects [19, 20] into account. This work also aims to
integrate nucleation models into the framework, to enrich boundary conditions
to broaden the scope of application, and to develop an efficient numerical
coupling with ThermoCalc software [21] to seamlessly extract the necessary
thermodynamic data.

An overview of the numerical formulation proposed in [16] is provided in the
next section, followed by the improvements and the generalizations
incorporated. A benchmark analysis of the numerical model, followed by
illustrations in large scale microstructures, are presented in the third
section. The model’s capabilities in simulating other diffusive solid-state
phenomena are also demonstrated. The concluding remarks and the prospective
works are then discussed in the last part.

##  2 Numerical formulation

The growth regime of a microstructural evolution involving diffusive phase
transformation is influenced by two coupled processes: (i) solute partitioning
between the phases involved in transformation, and (ii) subsequent migration
of the interface network along with other contributions to the interface
kinetics. As outlined in the previous publication [16], level-set (LS)
functions are employed to represent and govern the interface migration
aspects. In the context of a polycrystalline microstructure, a set of global
level set (GLS) functions, $\varphi_{i}$, are used to efficiently represent
and classify large number of grains based on the grain coloration theorem [17,
18]. However, in the context of multiphase microstructure, secondary LS
functions are necessary to distinguish the regions composed of different
phases. The $\varphi_{i}$ GLS functions constitute the primary basis of LS
representation, while the level-set functions used to classify the phases
forms the secondary basis. Fig.1a gives an illustration of a biphasic
microstructure with the color code indicating the family of GLS functions used
to depict different grains. Meanwhile, Fig.1b indicates the regions composed
of $\alpha$ and $\gamma$ phases, and the secondary basis $\varphi_{\alpha}$ LS
function employed to characterize them.

(a) Classification of different grains

(b) Classification of regions of different phases

Figure 1: Biphasic microstructure indicating the underlying the level-set
representation

The presence of different phases introduce material discontinuities.
Consequently, within the classical level-set framework, resolving solute
redistribution may require separate treatment in the domains of different
phases by explicitly taking into account the necessary jump conditions at the
interphase boundaries. Incorporating these jump conditions necessitates
localizing the evolving interphase boundaries at each resolution step, a task
facilitated by the secondary level-set function ($\varphi_{\alpha}$). As
iterated in [16], rather than relying on the sharp interface characteristics
(in theory) provided by the LS representation, we adopt a diffuse interface
hypothesis across the interphase boundaries while addressing solute diffusion
within the domain. This approach aims to smooth out discontinuities and
integrate a unified equation for modeling solute partitioning throughout the
domain. It’s worth noting that the secondary basis level-set function
$\varphi_{\alpha}$ is initially constructed from the input phase
characteristic function ($\chi_{\alpha}$) [16]. Then, as the evolution
proceeds, $\varphi_{\alpha}$ is updated accordingly. The transition from
$\varphi_{\alpha}$ to a diffuse interface representation ($\phi$) is achieved
using a hyperbolic tangent relation between them, which is extensively
explained in our previous publication [16]:

$\phi=\frac{1}{2}tanh\left(\frac{3\varphi_{\alpha}}{\eta}\right)+\frac{1}{2}.$
(1)

The $\eta$ in the above relation is a diffuse interface thickness parameter.

After resolving the solute diffusion equation, the resulting solute
distribution influences the kinetics of interface migration, alongside other
inherent effects depending on the various considered phenomena. The resulting
motion of the multi-phase grain interface network are then governed by the
transport of primary basis GLS functions ($\varphi_{i}$), based on their
classical level-set description.

###  2.1 Solute redistribution

Upon establishing a diffuse interface description, $\forall i\in\\{C,X\\}$, a
total concentration field variable ($x_{i}$) can be represented as a
continuous field between the concentrations in the individual phases
($x_{i}^{\alpha},x_{i}^{\gamma}$) for each solute $i$:

$x_{i}=\phi x_{i}^{\alpha}+(1-\phi)x_{i}^{\gamma}.$ (2)

As detailed in [16], through the imposition of solute flux continuity for each
phase and assuming a constant concentration ratio between the phases at the
interface throughout, the following unified diffusion equation is derived,
adhering to Fick’s laws of diffusion:

$\frac{\partial x_{i}}{\partial
t}=\bm{\nabla}\cdot\left\\{D^{*}(\phi)\left[\bm{\nabla}x_{i}-\frac{x_{i}(k_{p}^{i}-1)}{1+\phi(k_{p}^{i}-1)}\bm{\nabla}\phi\right]\right\\}.$
(3)

With further simplifications, the above Eq.(3) could be reformulated into the
following flux form:

$\frac{\partial x_{i}}{\partial
t}=\bm{\nabla}\cdot\left[D^{*}(\phi)\bm{\nabla}x_{i}-x_{i}\bm{\mathcal{A}}(\phi)\right].$
(4)

A Convective-Diffusive-Reactive (CDR) formulation is then obtained by
expanding the above equation:

$\frac{\partial x_{i}}{\partial
t}+\left(\bm{\mathcal{A}}-\bm{\nabla}D^{*}\right)\cdot\bm{\nabla}x_{i}-D^{*}\Delta
x_{i}+\mathcal{R}x_{i}=0.$ (5)

In the above CDR formulation:

  * •

$D^{*}$ is called the mixed diffusivity, expressed as:

$D^{*}(\phi)=\frac{D_{\gamma}^{i}+\phi(k_{p}^{i}D^{i}_{\alpha}-D^{i}_{\gamma})}{1+\phi(k_{p}^{i}-1)},$
(6)

where $D^{i}_{\alpha}$ and $D^{i}_{\gamma}$ represent the diffusivity of the
$i$th solute element in $\alpha$ and $\gamma$ phases, respectively.
$k_{p}^{i}$ is the constant concentration ratio respected by the
redistribution of the solute atoms between the parent and the product phase at
the interphase. This ratio is assumed to be equal to the equilibrium
partitioning ratio ($k_{p_{eq}}^{i}$), corresponding to the conditions at the
instant $t$:

$k_{p}^{i}\left(t\right)=\frac{x_{i}^{\alpha}}{x_{i}^{\gamma}}\approx
k_{p_{eq}}^{i}\left(t\right)=\frac{\left.x_{i}^{\alpha}\right|_{eq}}{\left.x_{i}^{\gamma}\right|_{eq}},$
(7)

where $\left.x^{\alpha}_{i}\right|_{eq}$ and
$\left.x^{\gamma}_{i}\right|_{eq}$ are the equilibrium concentrations of
$\alpha$ and $\gamma$ phases respectively at temperature $T$ corresponding to
an instant $t$.

  * •

$\bm{\mathcal{A}}$ is a part of the advective coefficient and is expressed as:

$\bm{\mathcal{A}}(\phi)=\frac{D^{*}(\phi)(k_{p}^{i}-1)}{1+\phi(k_{p}^{i}-1)}\bm{\nabla}\phi.$
(8)

  * •

$\mathcal{R}$ represents the reactive term and is defined as:

$\mathcal{R}=\bm{\nabla}\cdot\bm{\mathcal{A}}.$ (9)

For a test function, $\psi\in H^{1}(\Omega)$, the generic P1 finite element
weak formulation of Eq.(5) could be simplified to the following form:

$\begin{split}\int_{\Omega}\frac{\partial x_{i}}{\partial t}\psi\
d\Omega+\int_{\Omega}\bm{\mathcal{A}}\cdot\bm{\nabla}x_{i}\psi\
d\Omega+\int_{\Omega}D^{*}\bm{\nabla}\psi\cdot\bm{\nabla}x_{i}\ d\Omega\\\
+\int_{\Omega}\mathcal{R}x_{i}\psi\ d\Omega-\int_{\partial\Omega}\psi
D^{*}\bm{\nabla}x_{i}\cdot\bm{n}_{b}\ dS=0,\end{split}$ (10)

where $\bm{n}_{b}$ is the outward unit normal to the domain boundary and the
last boundary integral term is subject to the imposed boundary conditions.

In our previous works [16], we employed homogeneous Neumann conditions on
$x_{i}$ for solute mass conservation, i.e.,
$\left.\bm{\nabla}x_{i}\cdot\bm{n}_{b}\right|{\partial\Omega}=0$. At first
glance, these assumptions may seem reasonable, focusing on pure solute
diffusion within the computational domain without any influx or outflux of
solute atoms. However, these conditions are incomplete and are only valid in a
specific scenario. The mathematical formulation in Eq.(4) suggests that
applying homogeneous Neumann-type boundary conditions on $x_{i}$ conserves
solely the diffusive flux ($-D^{*}\bm{\nabla}x_{i}$), while the advective flux
($x_{i}\bm{\mathcal{A}}$) remains non-conserved.

If we apply the Reynolds’s transport theorem on $x_{i}$ for a fixed domain
boundary, for mass conservation:

$\frac{d}{dt}\int_{\Omega}x_{i}d\Omega=\int_{\Omega}\frac{\partial
x_{i}}{\partial t}d\Omega=0.$ (11)

Plugging Eq.(4), and applying the divergence theorem, we obtain:

$\frac{d}{dt}\int_{\Omega}x_{i}d\Omega=\int_{\Omega}\frac{\partial
x_{i}}{\partial
t}d\Omega=\int_{\partial\Omega}\left[D^{*}(\phi)\bm{\nabla}x_{i}-x_{i}\bm{\mathcal{A}}(\phi)\right]\cdot\bm{n}_{b}dS=0.$
(12)

Thus, the above equation holds true for:

$\left.\left[-D^{*}(\phi)\bm{\nabla}x_{i}+x_{i}\bm{\mathcal{A}}(\phi)\right]\cdot\bm{n}_{b}\right|_{\partial\Omega}=0.$
(13)

The aforementioned Robin-type condition constitutes the correct and complete
boundary condition for Eq.5.

Equation (8) implies that $\bm{\mathcal{A}}\rightarrow 0$ far from a phase
interface where $\bm{\nabla}\phi\rightarrow 0$. In this scenario, Eq.(13)
simplifies to a pure Neumann-type boundary condition, i.e.,
$\left.-D^{*}\bm{\nabla}x_{i}\cdot\bm{n}_{b}\right|_{\partial\Omega}=0$.
Therefore, a pure Neumann-type boundary condition is a particular case,
applicable only when the local boundary region does not interact with any
phase interface. In cases where phase interfaces approach certain sections of
the boundary region, the local advective flux should also be conserved through
a Robin type boundary condition as in Eq.(13).

Thus, for the resolution of the solute redistribution, there’s a continuous
local switch between homogeneous Neumann-type and Robin-type boundary
conditions. This switch is based on the interaction of phase interfaces with
the domain boundaries during the phase transformation.

Applying the generic Robin-type boundary condition to the weak form in
Eq.(10), we eventually obtain:

$\int_{\Omega}\frac{\partial x_{i}}{\partial t}\psi\
d\Omega+\int_{\Omega}\bm{\mathcal{A}}\cdot\bm{\nabla}x_{i}\psi\
d\Omega+\int_{\Omega}D^{*}\bm{\nabla}\psi\cdot\bm{\nabla}x_{i}\
d\Omega+\int_{\Omega}\mathcal{R}x_{i}\psi\
d\Omega-\int_{\partial\Omega}x_{i}\psi\bm{\mathcal{A}}\cdot\bm{n}_{b}\ dS=0.$
(14)

###  2.2 Interface migration

For a polycrystalline microstructure with $N_{g}$ grains whose interface
network is represented by a set of $N_{LS}$ GLS functions ($\varphi_{i}$), the
following convective equations need to be resolved for modeling the interface
migration:

$\begin{cases}\frac{\partial\varphi_{i}}{\partial
t}+\bm{v}\cdot\bm{\nabla}\varphi_{i}=0\\\
\varphi_{i}(\bm{x},t=0)=\varphi_{i}^{0}(\bm{x})\end{cases}\forall
i\in\\{1,2,...,N_{LS}\\}.$ (15)

The velocity field ($\bm{v}$) in this application is based as a function of
various driving pressures affecting the microstructural evolution.

In [16], the following Convective-Diffusive level-set transport equations were
formulated for grain/phase interface network migration for each of the $ith$
GLS function ($\varphi_{i}$) during a $\gamma\rightarrow\alpha$ phase
transformation:

$\begin{gathered}\frac{\partial\varphi_{i}}{\partial t}+\left[\bm{v_{\Delta
G_{\gamma\alpha}}}+\bm{v_{\llbracket\mathcal{E}\rrbracket}}\right]_{i}\cdot\bm{\nabla}\varphi_{i}-\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\sigma_{l}\right]\Delta\varphi_{i}=0\qquad\forall
i\in\\{1,2,...,N_{LS}\\},\end{gathered}$ (16)

where $\mathscr{S}=\\{\gamma\gamma,\gamma\alpha,\alpha\alpha\\}$ with
$\chi_{\gamma\gamma},\chi_{\gamma\alpha},\text{ and }\chi_{\alpha\alpha}$
being the interface characteristic functions, $M$ is the interface mobility,
and $\sigma$ is the interfacial energy. The convective velocity components are
defined as:

$\begin{gathered}\bm{v_{\Delta
G_{\gamma\alpha}}}=\chi_{\gamma\alpha}M_{\gamma\alpha}\Delta
G_{\gamma\rightarrow\alpha}\bm{n}\\\
\bm{v_{\llbracket\mathcal{E}\rrbracket}}=\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\llbracket\mathcal{E}\rrbracket_{l}\right]\bm{n}.\end{gathered},$
(17)

with $\bm{n}$ being the outward unit normal to a migrating interface.

$\Delta G_{\gamma\rightarrow\alpha}$ corresponds to the driving pressure for
phase transformation of $\gamma$ phase into $\alpha$ phase, and
$\llbracket\mathcal{E}\rrbracket$ constitutes a driving pressure linked to the
jump in stored energy ($\mathcal{E}$) due to plastic deformation.

The formulation in Eq.(16) corresponds to a scenario where the net driving
pressure ($P$) is assumed to be composed of the phase transformation
contribution, the stored energy component and the capillarity effects
($-\kappa\sigma$):

$P=\Delta
G_{\gamma\rightarrow\alpha}+\llbracket\mathcal{E}\rrbracket-\kappa\sigma.$
(18)

In a biphasic polycrystalline context, accounting for different kinds of
interfaces, this net driving pressure is rewritten in a generalized form as:

$\begin{gathered}P=\chi_{\gamma\alpha}\Delta
G_{\gamma\rightarrow\alpha}+\sum_{l\in\mathscr{S}}\chi_{l}\left(\llbracket\mathcal{E}\rrbracket_{l}-\kappa\sigma_{l}\right).\end{gathered}$
(19)

As a reminder, the interface migration velocity ($\bm{v}$) concerning a
microstructural evolution at the mesoscopic scale is assumed to be a product
of the interface mobility and the net local driving pressures describing the
various involved phenomena at the interface [22]: $\bm{v}=MP\bm{n}$. The
interface mobility is also expressed in a generalized form accounting for
different types of interfaces such that,
$M=\sum_{l\in\mathscr{S}}\chi_{l}M_{l}$.

In the context of Ternary alloys and beyond, comprising of slow diffusing
substitutional elements, an additional driving pressure contribution due to
the solute drag effects linked to these slow diffusing elements could be taken
into account. Accounting for this contribution, Eq.(19) could be modified as:

$\begin{gathered}P=\chi_{\gamma\alpha}\Delta
G_{\gamma\rightarrow\alpha}+\sum_{l\in\mathscr{S}}\chi_{l}\left(\llbracket\mathcal{E}\rrbracket_{l}-\kappa\sigma_{l}+\left[P_{SD}\right]_{l}\right),\end{gathered}$
(20)

where $P_{SD}$ is the solute drag pressure associated to the drag effects.

Upon prescribing the above modified interface kinetics into the transport
equations (Eqs.(15)), and with further simplifications using the metric
properties of a LS function, we obtain the following form $\forall
i\in\\{1,2,...,N_{LS}\\}$:

$\begin{gathered}\frac{\partial\varphi_{i}}{\partial t}+\left[\bm{v_{\Delta
G_{\gamma\alpha}}}+\bm{v_{\llbracket\mathcal{E}\rrbracket}}\right]_{i}\cdot\bm{\nabla}\varphi_{i}-\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\left[P_{SD}\right]_{l}-\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\sigma_{l}\right]\Delta\varphi_{i}=0,\end{gathered}$
(21)

where,

$\bm{n}_{i}=-\frac{\bm{\nabla}\varphi_{i}}{\|\bm{\nabla}\varphi_{i}\|}=-\bm{\nabla}\varphi_{i}\implies\kappa_{i}=\bm{\nabla}\cdot\bm{n}_{i}=-\Delta\varphi_{i}.$
(22)

The driving pressure due to solute drag contributes passively by opposing the
actively available interface kinetics (linked to the principally active
driving pressures), as the solutes tend to bind to the migrating interface due
to segregation. Thus, the magnitude of the solute drag pressure ($\Delta
G_{SD}=\lvert P_{SD}\rvert$) is generally assumed to be a function of the
interface migration velocity itself, $\Delta G_{SD}=f(\bm{v})$. This has the
effect of rendering the overall interface kinetics description non-linear. The
dependence on the velocity is such that, in the vicinity of zero velocity, or
at sufficiently high interfacial velocities, the drag pressure should
theoretically vanish. The latter is due to the fact that at high velocities,
there is not enough scope or time for the solute atoms to interact with the
rapidly migrating interface, thus reducing the potential to impart any
retardation effects of significance. For modeling the solute drag driving
pressure ($P_{SD}$), Cahn’s solute drag model [23] provides a simplified
description for the magnitude of the drag pressure ($\Delta G_{SD}$):

$\Delta
G_{SD}\left(\bm{v}\right)=\frac{\alpha_{C}x_{X}^{0}\left\lVert\bm{v}\right\rVert}{1+\beta_{C}^{2}\left\lVert\bm{v}\right\rVert^{2}},$
(23)

where $\alpha_{C}$ and $\beta_{C}$ are parameters defined as functions of
temperature, interfacial solute diffusivity, interface width and the binding
energy. Cahn proposed the following definitions for these parameters:

$\alpha_{C}=4N_{V}k_{B}T\int_{-\infty}^{+\infty}\frac{\sinh^{2}\left(\frac{E(z)}{2k_{B}T}\right)}{D^{X}_{\Gamma}(z)}\
dz,$ (24)

$\frac{\alpha_{C}}{\beta_{C}^{2}}=\frac{N_{V}}{k_{B}T}\int_{-\infty}^{+\infty}\left(\frac{\partial
E}{\partial z}\right)^{2}D^{X}_{\Gamma}(z)\ dz.$ (25)

Even though Cahn’s model lacks physical sense across interphase boundaries and
is more apt for grain boundaries, the simplicity in its description is
attractive for numerical implementation, especially in a FE framework,
considering the non-linearity of a solute drag pressure on the velocity field.
Cahn’s simplified drag pressure is capable of capturing the general trends
expected for solute drag effects whether it is for grain interfaces or phase
interfaces as illustrated in [24]. In first approximation, we shall thus
consider Cahn’s simplified description in this work to model solute drag
pressure.

The solute drag driving pressure component can be written as:

$P_{SD}\left(\bm{v},\bm{n}\right)=-\Delta
G_{SD}\frac{(\bm{v}\cdot\bm{n})}{\lVert\bm{v}\rVert}=-\frac{\alpha_{C}x_{X}^{0}(\bm{v}\cdot\bm{n})}{1+\beta_{C}^{2}\left\lVert\bm{v}\right\rVert^{2}}.$
(26)

With the application of Cahn’s solute drag model even for phase interfaces,
for the parameters,
$\left.\left(\alpha_{C},\beta_{C}\right)\right|_{l=\gamma\alpha}$, instead of
using the analytical definitions provided by Cahn in Eq.(24) and Eq.(25), the
idea is to consider them as a set of temperature dependent model parameters to
be fitted in accordance with relevant experimental results. This avoids the
complexity of precisely quantifying some of the physical parameters required
by their analytical definitions, especially across phase interfaces. So, by
selectively controlling the two parameters ($\alpha_{C},\beta_{C}$), we could
attempt to converge to expected magnitudes for the drag pressure across phase
interfaces.

Since the solute drag pressure renders the velocity field non-linear. if we
consider a fully implicit time discretization scheme for the resolution of
Eq.(21), we will obtain a non-linear formulation. Considering Euler implicit
time discretization, Eq.(21) could be discretized as:

$\frac{\partial\varphi_{i}}{\partial
t}\approx\frac{\varphi_{i}^{k+1}-\varphi_{i}^{k}}{\Delta
t}=\mathcal{F}(\varphi_{i}^{k+1}),\qquad\forall i\in\\{1,2,...,N_{LS}\\},$
(27)

where

$\mathcal{F}(\varphi_{i}^{k+1})=\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\sigma_{l}\right]\Delta\varphi_{i}^{k+1}+\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\left[P_{SD}(\bm{v}^{k+1})\right]_{l}-\left[\bm{v_{\Delta
G_{\gamma\alpha}}}+\bm{v_{\llbracket\mathcal{E}\rrbracket}}\right]_{i}\cdot\bm{\nabla}\varphi_{i}^{k+1},$
(28)

is a time dependent operator, with $k$ being the index of time stepping.

If one considers $P_{SD}^{k+1}$ in the second term in Eq.(28), one obtains
from Eq.(26),

$P_{SD}(\bm{v}^{k+1})=-\frac{\alpha_{C}x_{X}^{0}(\bm{v}\cdot\bm{n})^{k+1}}{1+\beta_{C}^{2}\left\lVert\bm{v}^{k+1}\right\rVert^{2}}.$
(29)

Since the interface migration velocity is normal to the interface, at any
position $\bm{x}$, for the signed distance $i$th GLS function, such that
$\varphi_{i}^{k+1}(\bm{x})\geq 0$, it is possible to express:

$(\bm{v}\cdot\bm{n})^{k+1}=\frac{\varphi_{i}^{k+1}-\varphi_{i}^{k}}{\Delta
t}\approx\frac{\partial\varphi_{i}}{\partial t},$ (30)

and,

$\left\lVert\bm{v}^{k+1}\right\rVert=\frac{\left|\varphi_{i}^{k+1}-\varphi_{i}^{k}\right|}{\Delta
t}.$ (31)

So, $P_{SD}^{k+1}$ would then become:

$P_{SD}^{k+1}=-\frac{\alpha_{C}x_{X}^{0}\left(\frac{\varphi_{i}^{k+1}-\varphi_{i}^{k}}{\Delta
t}\right)}{1+\beta_{C}^{2}\left(\frac{\varphi_{i}^{k+1}-\varphi_{i}^{k}}{\Delta
t}\right)^{2}}.$ (32)

Clearly the term in the denominator due to
$\left\lVert\bm{v}^{k+1}\right\rVert$ manifests the operator $\mathcal{F}$ in
Eq.(28) and hence the time discretized formulation in Eq.(27) non-linear. To
avoid this non-linearity and simplify the resolution procedure, as seen in the
works of [25], the following explicit scheme based assumption is made:

$\left\lVert\bm{v}^{k+1}(\bm{x})\right\rVert^{2}\approx\left\lVert\bm{v}^{k}(\bm{x})\right\rVert^{2}=\left(\frac{\varphi_{i}^{k}(\bm{x})-\varphi_{i}^{k-1}(\bm{x})}{\Delta
t}\right)^{2}=\left\lVert\bm{v}^{old}(\bm{x})\right\rVert^{2}\quad\text{ for
}i\mid\varphi_{i}^{k}(\bm{x})\geq 0,$ (33)

where $\bm{v}^{k}$ or $\bm{v}^{old}$ is the net migration velocity field
prescribed to compute the GLS solution $\varphi_{i}^{k}$ at the current state.
Even though this formulation deviates slightly from the fully implicit time
scheme, for the small ranges of time steps expected to be considered in the
GLS resolution, the errors accumulated are expected to be of the small order
while considerably simplifying the resolution procedure, thanks to the
linearization. Taking this into account, Eq.(32) could be rewritten as:

$P_{SD}^{k+1}\approx-\left(\frac{\alpha_{C}x_{X}^{0}}{1+\beta_{C}^{2}\left\lVert\bm{v}^{old}\right\rVert^{2}}\right)\left(\frac{\varphi_{i}^{k+1}-\varphi_{i}^{k}}{\Delta
t}\right).$ (34)

This would then linearize the operator $\mathcal{F}$ in Eq.(28) and hence the
time discretized formulation in Eq.(27). Although not exactly, but such a
strategy of splitting the operator $\mathcal{F}$ into implicit and explicit
parts weakly mimics the ideologies of IMplicit-EXplicit (IMEX) time
integration methods [26, 27]. Substituting Eq.(34) into the time discretized
formulation and going back to the continuum description in Eq.(21), we obtain:

$\begin{gathered}\frac{\partial\varphi_{i}}{\partial t}+\left[\bm{v_{\Delta
G_{\gamma\alpha}}}+\bm{v_{\llbracket\mathcal{E}\rrbracket}}\right]_{i}\cdot\bm{\nabla}\varphi_{i}+\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\left(\frac{\alpha_{C}x_{X}^{0}}{1+\beta_{C}^{2}\left\lVert\bm{v}^{old}\right\rVert^{2}}\right)_{l}\frac{\partial\varphi_{i}}{\partial
t}=\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\sigma_{l}\right]\Delta\varphi_{i}.\end{gathered}$
(35)

By imposing,

$\mathscr{M}_{SD}=1+\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\left(\frac{\alpha_{C}x_{X}^{0}}{1+\beta_{C}^{2}\left\lVert\bm{v}^{old}\right\rVert^{2}}\right)_{l},$
(36)

we obtain:

$\begin{gathered}\mathscr{M}_{SD}\frac{\partial\varphi_{i}}{\partial
t}+\left[\bm{v_{\Delta
G_{\gamma\alpha}}}+\bm{v_{\llbracket\mathcal{E}\rrbracket}}\right]_{i}\cdot\bm{\nabla}\varphi_{i}=\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\sigma_{l}\right]\Delta\varphi_{i}.\end{gathered}$
(37)

Eventually, we can express the above equation into the following linearized
convective-diffusive (CD) formulation $\forall i\in\\{1,2,...,N_{LS}\\}$:

$\begin{gathered}\frac{\partial\varphi_{i}}{\partial
t}+\frac{1}{\mathscr{M}_{SD}}\left[\bm{v_{\Delta
G_{\gamma\alpha}}}+\bm{v_{\llbracket\mathcal{E}\rrbracket}}\right]_{i}\cdot\bm{\nabla}\varphi_{i}-\frac{1}{\mathscr{M}_{SD}}\left[\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\sigma_{l}\right]\Delta\varphi_{i}=0.\end{gathered}$
(38)

In the following, we shall refer $\mathscr{M}_{SD}$ as the solute drag
pressure coefficient. When the solute drag effects are absent or negligible,
$\mathscr{M}_{SD}\rightarrow 1$, while when they are non-negligible,
$\mathscr{M}_{SD}>1$. In contrast to the previous formulation (Eq.(16))
without the solute drag aspects, the new formulation only seems to involve a
factor, $\mathscr{M}_{SD}$, which has the global effect of heterogeneously
lowering the interface mobility (or increasing the drag resistance) for
migration as a function of the local interface kinetics.

In [16], it has been highlighted that a better description for $\bm{v_{\Delta
G_{\gamma\alpha}}}$ is essential in the context of polycrystals comprising of
multi-junctions, and thus the following description was proposed:

$\bm{v_{\Delta
G_{\gamma\alpha}}}(\bm{x},t)=\sum_{i=1}^{N_{LS}}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{N_{LS}}\chi_{G_{i}}M_{ij}\exp{\left(-\beta_{e}|\varphi_{j}|\right)}\chi_{\gamma\alpha}\Delta
G_{\gamma\rightarrow\alpha}\mathscr{F}_{s}(-\bm{n}_{j}),$ (39)

However, the use of all the $N_{LS}$ level-set functions is not necessary
unlike in the case of stored energy component, since $\Delta
G_{\gamma\rightarrow\alpha}$ is valid only on the phase interfaces and not on
the entire grain/phase interface network. So, an alternative and simple
approach could be proposed that takes into account only the secondary basis
level-set function representing all the $\alpha$ phase grains,
$\varphi_{\alpha}$:

$\bm{v_{\Delta
G_{\gamma\alpha}}}(\bm{x},t)=M_{\gamma\alpha}\chi_{\gamma\alpha}\Delta
G_{\gamma\rightarrow\alpha}(\bm{x},t)(-\bm{\nabla}\varphi_{\alpha}).$ (40)

#### Computation of the solute drag pressure coefficient, $\mathscr{M}_{SD}$:

To compute the solute drag pressure coefficient in Eq.(36):

$\mathscr{M}_{SD}=1+\sum_{l\in\mathscr{S}}\chi_{l}M_{l}\left(\frac{\alpha_{C}x_{X}^{0}}{1+\beta_{C}^{2}\left\lVert\bm{v}^{old}\right\rVert^{2}}\right)_{l}\quad\forall
l\in\\{\gamma\gamma,\gamma\alpha,\alpha\alpha\\},$

we need the respective laws/ expressions governing
$\left\\{\alpha_{C}(T),\beta_{C}(T)\right\\}_{l}$ and $M_{l}(\bm{x},t,T)$
provided, and $\chi_{l}$ and $\lVert\bm{v}^{old}\rVert$ need to be computed.
The computation of $\chi_{l}$ functions is already discussed in [16]. In order
to compute $\lVert\bm{v}^{old}\rVert$, let us consider its discrete
definition:

$\left\lVert\bm{v}^{old}(\bm{x})\right\rVert=\frac{\left\lvert\varphi_{i}^{k}(\bm{x})-\varphi_{i}^{k-1}(\bm{x})\right\rvert}{\Delta
t}\quad\text{ for }i\mid\varphi_{i}^{k}(\bm{x})\geq 0,$ (41)

where $\varphi_{i}^{k}$ and $\varphi_{i}^{k-1}$ are the GLS solutions
available at the current and the previous states of time, respectively. In
practice, the above definition is reformulated as:

$\left\lVert\bm{v}^{old}(\bm{x})\right\rVert=\max\limits_{i\in\\{1,...,N_{LS}\\}}\left(\mathcal{I}_{i}^{k}(\bm{x})\frac{\left\lvert\varphi_{i}^{k}(\bm{x})-\varphi_{i}^{k-1}(\bm{x})\right\rvert}{\Delta
t}\right),$ (42)

where $\mathcal{I}_{i}^{k}(\bm{x})$ is used as an indicator function
associated to the $i$th GLS function such that $\varphi_{i}^{k}(\bm{x})\geq
0$, i.e.,

$\mathcal{I}_{i}^{k}(\bm{x})=\begin{cases}1&\text{if}\quad\varphi_{i}^{k}(\bm{x})\geq
0\\\ 0&\text{otherwise}\end{cases}.$ (43)

At the initial state of time ($k=0$), $\lVert v^{old}\rVert=0$, as
$\varphi_{i}^{0}=\varphi_{i}^{-1}$ The solute drag pressure coefficient needs
to be computed before the next resolution step of Eq.(38).

###  2.3 ThermoCalc Coupling

As described in [16], the driving pressure related to the diffusive phase
transformation ($\Delta G_{\gamma\rightarrow\alpha}$) is prescribed through a
linearization of the phase diagram using a thermodynamic database such as
ThermoCalc [21]. The driving pressure is subsequently expressed as a function
of the current temperature, and the phase interface concentrations as follows:

$\begin{gathered}\Delta
G_{\gamma\rightarrow\alpha}(T,x^{\alpha}_{C},x^{\gamma}_{C})=\Delta
S^{\gamma\alpha}\left[(T^{R}-T)+0.5m^{\alpha/(\alpha+\gamma)}_{A-C}\left(x^{\alpha}_{C}-\left.x^{\alpha}_{C}\right|_{R}\right)\right.\\\
+\left.0.5m^{\gamma/(\alpha+\gamma)}_{A-C}\left(x^{\gamma}_{C}-\left.x^{\gamma}_{C}\right|_{R}\right)\right],\end{gathered}$
(44)

where $C$ is any interstitial element (Carbon for instance). $T^{R}$ is the
reference temperature at which the phase diagram is locally linearized, with
$\left(\left.x^{\alpha}_{C}\right|_{R},\left.x^{\gamma}_{C}\right|_{R}\right)$
being the corresponding equilibrium concentrations.
$m^{\alpha/(\alpha+\gamma)}_{A-C}$, and $m^{\gamma/(\alpha+\gamma)}_{A-C}$ are
the local slopes of the linearized phase boundaries, and $\Delta
S^{\gamma\alpha}$ is the entropy difference between the two phase at the
temperature $T$.

From Eqs.(2) and (7), the above description is rewritten in the form of the
total concentration variable:

$\begin{gathered}\Delta G_{\gamma\rightarrow\alpha}(T,x_{C})=\Delta
S^{\gamma\alpha}\left[T^{R}-T+0.5m^{\alpha/(\alpha+\gamma)}_{A-C}\left(\frac{k_{p}^{C}x_{C}}{1+\phi(k_{p}^{C}-1)}-\left.x^{\alpha}_{C}\right|_{R}\right)\right.\\\
+\left.0.5m^{\gamma/(\alpha+\gamma)}_{A-C}\left(\frac{x_{C}}{1+\phi(k_{p}^{C}-1)}-\left.x^{\gamma}_{C}\right|_{R}\right)\right]\end{gathered}$
(45)

In the context of ternary alloys ($A-C-X$), if we assume ortho-equilibrium
(complete local-equilibrium) between the phase for both the solutes $C$ and
$X$, the driving pressure could be expressed using the same principles of
phase diagram linearization [13, 24] as follows:

$\begin{gathered}\Delta G_{\gamma\rightarrow\alpha}(T,x_{C},x_{X})=\Delta
S^{\gamma\alpha}\left(T^{R}-T\right)+0.5\Delta
S^{\gamma\alpha}\sum_{q\in\\{C,X\\}}\left[m^{\alpha/(\alpha+\gamma)}_{A-q}\left(\frac{k_{p}^{q}x_{q}}{1+\phi(k_{p}^{q}-1)}-\left.x^{\alpha}_{q}\right|_{R}\right)\right.\\\
+\left.m^{\gamma/(\alpha+\gamma)}_{A-q}\left(\frac{x_{q}}{1+\phi(k_{p}^{q}-1)}-\left.x^{\gamma}_{q}\right|_{R}\right)\right].\end{gathered}$
(46)

Ortho-equilibrium hypothesis is feasible if both $C$ and $X$ are interstitial
solutes or if the rates of diffusion of $C$ and $X$ are fairly consistent.

However, if the ternary alloy is composed of a slow-diffusing substitutional
element ($X$) along with a faster diffusing interstitial element ($C$), ortho-
equilibrium may not be feasible as it could be difficult to simultaneously
satisfy the mass balances for both $C$, and $X$. In such a scenario, a
constrained phase equilibria is assumed. Within this research framework, the
para-equilibrium (PE) hypothesis [9] is selected as the constrained phase
equilibrium model. According to PE conditions, the slow-diffusing element $X$
is presumed to be completely immobile, allowing only the fast-diffusing
element $C$ to redistribute between the involved phases. Thus, local
equilibrium is maintained solely for the $C$ element. This hypothesis is
generally applicable when the diffusion of $X$ element is not possible within
the experimental timeframe or the selected heat treatment. The driving
pressure is then derived by extracting the hypothetical phase diagrams
obtained under the PE conditions and subsequently linearizing them as
previously:

$\begin{gathered}\Delta G_{\gamma\rightarrow\alpha}^{PE}(T,x_{C})=\Delta
S^{\gamma\alpha}_{PE}\left[T^{R}-T+0.5m^{\alpha/(\alpha+\gamma)}_{[A-C]_{PE}}\left(\frac{k_{p}^{C}x_{C}}{1+\phi(k_{p}^{C}-1)}-\left.x^{\alpha}_{C}\right|_{R}^{PE}\right)\right.\\\
+\left.0.5m^{\gamma/(\alpha+\gamma)}_{[A-C]_{PE}}\left(\frac{x_{C}}{1+\phi(k_{p}^{C}-1)}-\left.x^{\gamma}_{C}\right|_{R}^{PE}\right)\right],\end{gathered}$
(47)

where $m^{\alpha/(\alpha+\gamma)}_{[A-C]_{PE}}$ and
$m^{\gamma/(\alpha+\gamma)}_{[A-C]_{PE}}$ are the slopes of the linearized
imaginary solvus lines constructed under PE assumptions, while
$\left.x^{\alpha}_{C}\right|_{R}^{PE}$ and
$\left.x^{\gamma}_{C}\right|_{R}^{PE}$ are the corresponding constrained
equilibrium concentrations at $T^{R}$.

In previous works [16], the interaction with ThermoCalc was manual, limiting
the number of reference temperatures for data extraction. Furthermore, the
ThermoCalc console or its graphical user interface does not facilitate the
seamless calculation of certain thermodynamic data under para-equilibrium
conditions. These difficulties were overcome by establishing a coupling
between the numerical model and the TQ-interface (SDK) of ThermoCalc to enable
automatic computation of relevant thermodynamic data over a large set of
reference points under both ortho-equilibrium as well as para-equilibrium
conditions. Fig.2 gives a general overview of the ThermoCalc coupling
framework using the TQ-interface. For a given alloy system, the desired
thermodynamic conditions for data extraction are passed as inputs from the
numerical model to the ThermoCalc coupling interface. The coupling model then
provides a matrix of desired thermodynamic data for a given range of
temperatures over a large set of points. The ThermoCalc data extraction is
typically performed as a pre-processing step in the numerical framework. The
discrete and rich thermodynamic dataset generated is then utilized in the
computational loop with the assistance of an interpolation model, enabling the
interpolation of necessary thermodynamic data for any temperature at each time
step. This approach implies local piecewise linearization of the phase diagram
at the current temperature, $T^{R}=T$. This is crucial for accurately
capturing strong topological evolution of the solvus surfaces (or lines) in a
phase diagram, particularly for non-isothermal or continuous cooling
transformations involving a long thermal path
($\left|T_{initial}-T_{final}\right|$). This helps to enrich the driving
pressure description that is based on strong assumptions, and potentially
minimize errors. Interested readers are encouraged to refer [24] for further
specifications on the coupling.

Figure 2: General overview of the ThermoCalc coupling with the numerical model

###  2.4 Nucleation model

Precisely modeling nucleation events at the mesoscopic scale is complex and
often more suitable for much lower scales of modeling, such as atomistic
simulations. The literature offers a wide array of hypotheses and parameters
aimed at characterizing nucleation events in phase transformation. These
parameters encompass various aspects such as the nucleation start temperature
($T_{N_{s}}$), nucleation temperature range ($\delta T_{N}$), shield distance
($\delta l_{N}$), shield time ($\delta t_{N}$), parameters related to nucleus
shape, and composition-based factors influenced by local parent phase
features, among others. Additionally, adhering to the Continuous Nucleation
Theory (CNT) introduces supplementary parameters such as incubation time and
frequency factor. However, due to experimental constraints in observing
nucleation events, quantifying most of these parameters independently presents
significant challenges. Therefore, in this study, we will employ a simplified
approach, aiming for potential refinements in subsequent research. The nucleus
shape will be assumed to be spherical, and only a select few parameters from
the aforementioned list will be taken into account.

The most favored sites for nucleating the new phase are the grain corners
(triple junctions in 2D, quadruple junctions in 3D), followed by grain edges
and surfaces. At lower cooling rates, grain corners are preferred. However, at
higher cooling rates, the nucleation density is expected to increase, raising
the likelihood of saturating the grain corners for nucleation. This situation
consequently opens up opportunities for nucleation on grain edges and surfaces
[28, 29]. Fig.3 gives a 2D representation of nucleation of $\alpha$ phase on
grain corners and grain boundaries in the parent $\gamma$ phase.

Figure 3: Illustration of spherical nuclei of $\alpha$ phase (blue) nucleating
in a 2D parent $\gamma$ phase

The nucleation start temperature ($T_{N_{s}}$) is one of the parameters used
to qualify the onset of nucleation of the new phase. As the temperature
decreases below the transformation temperature ($T_{A3}$), the potential
driving pressure for initiating the stable phase’s nucleation typically
increases. Considering capillary effects, the likelihood of a critical nucleus
forming near the $T_{A3}$ temperature is relatively low. Conversely, as the
temperature decreases further and more driving pressure becomes available, the
likelihood of nucleation increases. Therefore, it’s logical to assume that
$T_{N_{s}}$ is generally set slightly below the transformation temperature.
This difference, known as the nucleation undercooling ($\Delta
T_{N}=T_{A3}-T_{N_{s}}$), is often determined practically by observing the
temperature at which a certain initial percentage of the product phase
fraction is experimentally detected, as mentioned in literature such as [29].

Nucleation can occur in a site-saturated or continuous manner. In site-
saturated nucleation, all nuclei form simultaneously below the nucleation
start temperature, while continuous nucleation involves nuclei forming over a
certain period. In the context of continuous cooling transformation, the
latter mode of nucleation translates to nucleation occurring over a particular
temperature span, referred as nucleation temperature range ($\delta T_{N}$).
In a study by Offerman et al. [28] investigating nucleation, the existence of
a distinct temperature range is suggested, implying a preference for
continuous mode of nucleation over site-saturated scenarios. They also noted a
correlation between the cooling rate and this temperature range. Specifically,
they observed that $\delta T_{N}$ becomes more pronounced with higher cooling
rates while being relatively smaller at lower cooling rates. $\delta T_{N}$ is
generally assumed to be a modeling parameter [29] adjusted to align with
relevant experimental results.

From a numerical modeling perspective, a site characteristic function
($\chi_{sn}$) is adopted to indicate the sites (and their computational nodes)
available for nucleation:

$\chi_{sn}(\bm{x},t)=\begin{cases}1&\text{ if }\bm{x}\in\text{Choice of
nucleation site}\\\ 0&\text{otherwise}\end{cases},\quad\forall t\mid
T_{N_{s}}-\delta T_{N}<T(t)\leq T_{N_{s}}.$ (48)

If the chosen nucleation site is on the grain boundaries, then
$\chi_{sn}(\bm{x},t)=\chi_{\gamma\gamma}$ qualifies all the eligible nodes on
the $\gamma$-phase grain boundaries as illustrated in green in Fig.4.

Figure 4: Characterization of computational nodes for grain boundary
nucleation in a 2D context

Alternatively, when opting for grain corners as the choice of nucleation site,
the nodes are qualified as follows:

$\epsilon_{MJ}$, $\varphi_{i}$

$\text{count }\leftarrow 0$

$i\leftarrow 0$

for $\bm{x}$ in $\Omega$ do

while $i<N_{LS}$ do

if $\lvert\varphi_{i}(\bm{x})\rvert\leq\epsilon_{MJ}$ then

$\text{count }\leftarrow\text{count }+1$

end if

end while

if $\text{count }\geq d_{i}+1$ then$\triangleright$ $d_{i}$ is the
dimensionality

$\chi_{sn}(\bm{x})\leftarrow 1$

else

$\chi_{sn}(\bm{x})\leftarrow 0$

end if

end for

Here, $\epsilon_{MJ}$ characterizes a small distance threshold around the
multiple junction. In 3D, if a node has at least $4$ GLS functions that meet
the condition, it qualifies as a multi-junction node available for nucleation.
In 2D, the equivalent criterion is fulfilled with $3$ GLS functions. A similar
criterion can be applied to qualify grain edge nodes in 3D, where the
condition is met by $3$ GLS functions. Fig.5 gives a 2D representation of the
small localized region in green around the triple junction satisfying the
above condition and hence the nodes within this zone are provisionally
qualified as eligible sites.

The site characteristic function is however utilized in conjunction with a
shield distance parameter, $\delta l_{N}$. As discussed in [29], shield
distance is an inter-nucleation distance parameter, shielding pre-existing
nuclei from the formation of additional ones in their proximity. From a
physical perspective, this parameter reflects that upon nucleus formation, the
surrounding local region becomes less conducive to further nucleation due to
subsequent changes in local characteristics, such as solute concentration.
This parameter primarily serves as a modeling parameter due to the challenges
in consistently quantifying it, stemming from the lack of substantial physical
underpinning. In the current model, at any instant $t$, $\delta l_{N}$ is
taken as, $\delta l_{N}(t)=k_{sh}r^{*}_{num}(t)$, where $k_{sh}$ is a shield
distance factor, and $r^{*}_{num}$ is the numerically imposed critical radius
of the nucleus.

In many scenarios, the critical radii of nuclei are significantly smaller
compared to the size of the considered domain. Consequently, accurately
capturing these small nuclei poses substantial computational demands, leading
to them being rarely precisely prescribed at their theoretical critical radii.
To address this, a minimum radius is set based on the underlying mesh
resolution, while a maximum size is imposed to prevent uncharacteristically
large nuclei. Thus, the numerical radius of a nucleus ($r^{*}_{num}$) in this
model is determined by the following expression:

$r^{*}_{num}=\min\left(\max\left(r^{*}_{th},k_{1}h_{min}\right),k_{2}h_{min}\right),$
(49)

where $r^{*}_{th}$ is the theoretically estimated critical radius, $h_{min}$
is the minimum size of the underlying mesh resolution, $k_{1}$ is a constant
governing the minimum number of mesh elements to be included along the radius,
and $k_{2}(>k_{1})$ represents another constant limiting the maximum size
permissible for a nucleus. The theoretical critical radius technically depends
on the site of nucleation, however in this context, it is roughly estimated as
follows:

$r^{*}_{th}=\frac{(d_{i}-1)\sigma_{\gamma\alpha}}{\Delta g^{\alpha}_{V}},$
(50)

where $d_{i}$ denotes the dimensionality, and $\Delta g^{\alpha}_{V}(T)$ is
the driving pressure for nucleation at the current temperature. This driving
pressure can be extracted during the pre-processing stage, thanks to the
ThermoCalc coupling, and subsequently interpolated at any temperature as
previously discussed.

Therefore, the shield distance parameter, combined with the site
characteristic function, further filters out nodes that are ineligible for
nucleation (i.e., $\varphi_{\alpha}(\bm{x})\geq-\delta l_{N}$). The shielded
zone and the nodes within it are highlighted in red in Fig.5. So, eventually
one of the nodes within the green zone is shielded and hence marked in red as
ineligible for nucleation.

Figure 5: Characterization of nucleation sites for grain corner nucleation:
eligible nodes (green), shielded ineligible nodes (red)

The functioning of the nucleation model adopted within this framework is
summarized in Fig.6.

Figure 6: General workflow of the nucleation model

Currently, the model adopts a constant nucleation rate for the continuous
nucleation scenario. The implementation of Continuous Nucleation Theory (CNT)
has not been incorporated into this work. The nucleation rate ($\dot{n_{N}}$)
is estimated using the cooling rate $\dot{T}$, the imposed nucleation density
($\rho_{N}$) and the nucleation temperature range as follows:

$\dot{n_{N}}=\frac{\rho_{N}V_{0}\dot{T}}{\delta T_{N}},$ (51)

where $V_{0}$ represents the total volume of the domain (or area $A_{0}$ in
2D). The nucleation density could be approximated using the expected final
fraction of the product phase (e.g., $f^{\alpha}_{f}$) and its anticipated
average grain size ($r_{\alpha}^{f}$) as outlined in [29]:

$\rho_{N}=\frac{3f^{\alpha}_{f}}{4\pi(r_{\alpha}^{f})^{3}}.$ (52)

In 2D, the above relation is modified to
$\rho_{N}=\frac{f^{\alpha}_{f}}{\pi(r_{\alpha}^{f})^{2}}$.

Upon the creation of a nucleus, its introduction triggers a local change in
the phase status, subsequently altering the local diffusivity field. This
alteration allows the solute redistribution equation to naturally adjust the
local composition. Thus, the local composition isn’t manually altered upon
nucleus introduction to avoid tampering with the solute mass conservation. At
the moment of creation, the newly formed nucleus is dynamically integrated
into one of the existing GLS functions, $\varphi_{i}$, or a new function is
created if needed. As the local phase features are updated, the associated
dependent fields, including the secondary basis function $\varphi_{\alpha}$,
and the diffuse interface function ($\phi$), etc., are correspondingly
modified.

###  2.5 Overview of the numerical framework

Fig.7 presents a general outline of the numerical framework highlighting the
key stages and tasks involved. The numerical simulations are carried out in a
FE context using unstructured triangular meshes with a P1 interpolation, and
employing an implicit backward Euler time scheme for the time discretization.
Each system associated with Eq.(10) and the weak formulation of Eqs.(38) is
assembled using typical P1 finite elements with a Streamline Upwind Petrov-
Galerkin (SUPG) stabilization for the convective terms [30]. A particular
adaptive meshing/ remeshing strategy is adopted along the grain and interphase
boundaries, outlined in [16]. The boundary conditions applied to the GLS
functions ($\varphi_{i}$) are classical null Neumann boundary conditions,
while the total solute concentration field is imposed with Robin-type
conditions. It should be highlighted that, post the resolution of the LS
transport equations, a particular numerical treatment [31] is performed on the
$\varphi_{i}$ functions to avoid any kinematic incompatibilities at the
multiple junctions:

$\hat{\varphi_{i}}=\frac{1}{2}\left(\varphi_{i}-\max_{j\neq
i}\varphi_{j}\right)\quad\forall i\in\\{1,...,N_{LS}\\}.$ (53)

It is then followed by a reinitialization procedure of the GLS functions using
a recent efficient strategy [32] to restore the metric properties of a signed
distance function.

Figure 7: General overview of the numerical framework

##  3 Results and discussion

###  3.1 Benchmarking

To validate this numerical framework, the semi-analytical mixed-mode (SAM)
sharp interface model proposed by Bos et al. [33] for diffusive phase
transformation is used. Furthermore, the results are compared with those
obtained from a similar analysis conducted by Mecozzi et al. [34] using their
phase-field numerical-model (PF-NM) framework.

In the context of $\gamma\rightarrow\alpha$ phase transformation, the semi-
analytical formulation assumes an analytical function of position in a semi-
infinite domain to model the solute concentration profile in front of the
$\alpha/\gamma$ interface. The assumed concentration profile is then
constrained by flux and mass balance relations. The diffusion in the product
$\alpha$ phase is assumed to be instantaneous, and hence the concentration
profile within this phase is set to be homogeneous at its equilibrium
concentration. Fig.8 depicts the representation of the concentration profile
in a 1D domain within the sharp interface based semi-analytical model. A
single $\alpha$ phase nucleus is assumed to form in the infinite $\gamma$
phase matrix.

Figure 8: Solute concentration profile assumed in the semi-analytical model by
Bos et al. [33]

The driving pressure for phase transformation is assumed to be proportional to
the deviation of $\gamma$ phase interface concentration
($x^{\gamma}_{\Gamma}$) from its equilibrium concentration
($x^{\gamma}_{eq}$):

$\Delta
G_{\gamma\rightarrow\alpha}=\Upsilon\left(x^{\gamma}_{eq}-x^{\gamma}_{\Gamma}\right),$
(54)

where $\Upsilon$ is a temperature dependent proportionality factor derived
from thermodynamic databases. Interested readers are referred to [33, 24] for
further details on the formulation and the resolution procedure.

#### Simulation conditions:

An isothermal austenite decomposition ($\gamma\rightarrow\alpha$) in a steel
alloy (Fe - $0.1\text{\,}\mathrm{wt\%}$C - $0.5\text{\,}\mathrm{wt\%}$Mn) at
$1000\text{\,}\mathrm{K}$ is considered for benchmarking. This scenario
mirrors the case studied by Mecozzi et al. [34] in their comparison of a
phase-field numerical model with the aforementioned semi-analytical method.
The phase transformation is assumed to occur under the para-equilibrium (PE)
hypothesis. Notably, this study does not consider the solute drag effects
potentially induced by the Mn element, and the capillarity effects at the
phase interface are neglected for this analysis. For this particular case
study, following [34], we impose an additional hypothesis of equal
undercoolings in both the phases when deriving the driving pressure via phase
diagram linearization. This added assumption enables the use of a coherent
thermodynamic description between the sharp interface model and the numerical
model. So, the driving pressure description in Eq.(44) for the numerical model
could be simplified to the following form:

$\Delta G_{\gamma\rightarrow\alpha}=\Delta
S^{\gamma\alpha}m^{\gamma/(\alpha+\gamma)}_{Fe-C}\left(x^{\gamma}_{C}-\left.x^{\gamma}_{C}\right|_{eq}\right),$
(55)

where $x^{\gamma}_{C}$ is the interface concentration in the $\gamma$ phase.
By comparing the above equation with the driving pressure formulation of the
semi-analytical sharp interface model as described in Eq.(54), we can deduce
that, $\Upsilon=-m^{\gamma/(\alpha+\gamma)}_{Fe-C}\Delta S^{\gamma\alpha}$.
The above derivation under the equal undercooling assumption is thoroughly
outlined in [34, 24].

The required mobility, carbon diffusivity, and thermodynamic data have been
adopted from [34] to maintain consistency with the comparison. Therefore, the
temperature dependency of phase interface mobility and carbon diffusivity in
the austenite phase has been considered based on an Arrhenius-type law as
follows:

$\begin{gathered}M_{\gamma\alpha}=3.5\times
10^{17}\exp\left(\frac{-Q_{m}}{RT}\right)\quad
in\quad$\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{4}\text{$\cdot$}{\mathrm{J}}^{-1}\text{$\cdot$}{\mathrm{s}}^{-1}$\\\
D_{\gamma}^{C}=1.5\times
10^{7}\exp\left(\frac{-Q^{C}_{\gamma}}{RT}\right)\quad
in\quad$\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{2}\text{$\cdot$}{\mathrm{s}}^{-1}$\end{gathered},$ (56)

where $Q_{m}=$140\text{\,}\mathrm{kJ}\text{$\cdot$}{\mathrm{mol}}^{-1}$$, and
$Q^{C}_{\gamma}=$142.1\text{\,}\mathrm{kJ}\text{$\cdot$}{\mathrm{mol}}^{-1}$$.

The thermodynamic data extracted according to [34], with a reference
temperature of $T^{R}=$1073\text{\,}\mathrm{K}$$, have been tabulated in
Tab.1. These data have been derived by linearizing the pseudo-binary Fe-C
phase diagrams under the para-equilibrium hypothesis.

$T^{R}$ ($\text{\,}\mathrm{K}$) |  $\left.x_{C}^{\alpha}\right|_{R}$ ($\text{\,}\mathrm{wt\%}$) |  $\left.x_{C}^{\gamma}\right|_{R}$ ($\text{\,}\mathrm{wt\%}$) |  $\Delta S^{\gamma\alpha}$ ($\text{\,}\mathrm{J}\text{$\cdot$}{\mathrm{K}}^{-1}\text{$\cdot$}{\mathrm{\SIUnitSymbolMicro m}}^{-3}$) |  $m^{\alpha/(\alpha+\gamma)}_{Fe-C}$ ($\text{\,}\mathrm{K}\text{$\cdot$}{\mathrm{wt\%}}^{-1}$) |  $m^{\gamma/(\alpha+\gamma)}_{Fe-C}$ ($\text{\,}\mathrm{K}\text{$\cdot$}{\mathrm{wt\%}}^{-1}$)
---|---|---|---|---|---
$1073$ |  $0.009$ |  $0.279$ |  $3.46\times 10^{-13}$ |  $-10250$ |  $-186.2$
Table 1: ThermoCalc data extracted at $T^{R}=$1073\text{\,}\mathrm{K}$$ from
[34] for Fe - $0.1\text{\,}\mathrm{wt\%}$C - $0.5\text{\,}\mathrm{wt\%}$Mn

In this comparison, two scenarios are examined: (i) A pseudo-1D domain
measuring $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in length (with a
relatively small and insignificant breadth of
$0.15\text{\,}\mathrm{\SIUnitSymbolMicro m}$), and (ii) A square 2D domain of
size $100\times 100$ $\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$.
Initially, the domain is entirely austenitic. Upon instantaneous cooling below
the $T_{A3}$ temperature to $1000\text{\,}\mathrm{K}$, a ferrite nucleus is
allowed to grow at the center of the domain for a duration of
$200\text{\,}\mathrm{s}$. In the LS numerical-model (LS-NM), the initial
ferrite nucleus radius ($r_{\alpha}^{i}$) is set to
$0.25\text{\,}\mathrm{\SIUnitSymbolMicro m}$ due to numerical restrictions
associated with the FE mesh size (whereas it starts from
$0\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in the SAM). For the LS-NM, the
initial ferrite nucleus composition is established at
$0.001832\text{\,}\mathrm{wt\%}$C, which is the critical nucleus composition
at $T=$1000\text{\,}\mathrm{K}$$ under the para-equilibrium (PE) constraints,
extracted from ThermoCalc. However, in the SAM, the nucleus composition is
directly set to the equilibrium composition for the ferrite phase
($0.016013\text{\,}\mathrm{wt\%}$C) under PE, assuming instantaneous diffusion
in this phase. Throughout the numerical simulation, a uniform time step of
$\Delta t=$0.01\text{\,}\mathrm{s}$$ is utilized. An adaptive
meshing/remeshing strategy is employed with a diffuse interface thickness of
$\eta=$0.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ for the 1D case and
$\eta=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ for the 2D case. The mesh
resolution within the diffuse interface ($h_{min}$) is specified as
$5\text{\,}\mathrm{nm}$ for the 1D case and $70\text{\,}\mathrm{nm}$ for the
2D case.

#### Results:

Figs.9 and Figs.10 illustrate the phase distribution at the initial state and
at the end of $t=$200\text{\,}\mathrm{s}$$, obtained by the LS-NM for the 1D
and the 2D cases, respectively.

(a) $t=$0\text{\,}\mathrm{s}$$

(b) $t=$200\text{\,}\mathrm{s}$$

Figure 9: Phase distribution obtained by the level-set based numerical model -
1D case

(a) $t=$0\text{\,}\mathrm{s}$$

(b) $t=$200\text{\,}\mathrm{s}$$

Figure 10: Phase distribution obtained by the level-set based numerical model
- 2D case

Fig.11 depicts a comparison between the LS-NM and the SAM concerning the
evolution of carbon concentration profiles ahead of the $\alpha/\gamma$ phase
interface within the austenite side at different time intervals for the 1D
scenario. The solid lines denote the numerical solution, whereas the open
circular markers denote the semi-analytical solution. Similarly, Fig.12
presents this comparison for the 2D case.

Figure 11: Evolution of the carbon concentration profiles in the austenite
phase - 1D case Figure 12: Evolution of the carbon concentration profiles in
the austenite phase - 2D case analyzed along the A-A line described in Fig.10

The curves in Fig.11 and Fig.12 demonstrate good agreement at different
instants between the LS-NM and the SAM in predicting the solute concentration
profiles in the austenite phase for both the 1D and the 2D cases (analyzed
along the A-A line described in Fig.10). In the 2D case, a slight deviation
between the two models at $t=$200\text{\,}\mathrm{s}$$ is noticeable. This
discrepancy likely stems from the effects of boundary conditions in the
numerical model as the phase interface approaches the boundaries of the 2D
domain. In contrast, in the 1D case, the phase interface predicted by the
numerical model at $t=$200\text{\,}\mathrm{s}$$ is further away from the
boundaries, thus unaffected by boundary conditions. It’s important to note
that while the semi-analytical model operates within a semi-infinite domain,
the numerical model operates within a bounded domain with imposed boundary
conditions. Analyzing these concentration profiles reveals the mixed-mode
nature of the transformation kinetics. It’s clear that the interfacial
concentration doesn’t immediately attain the local equilibrium concentration,
leading to finite interface migration. At the same time, diffusion in the bulk
of the austenite phase doesn’t occur instantaneously, resulting in a gradient
between the interfacial concentration and the bulk concentration.

In Figs.13, the evolution of the radius of the ferrite phase nucleus between
the SAM, LS-NM, and the PF-NM is compared. On the other hand, the evolution of
the average interfacial carbon concentration in the austenite phase across
these different models is compared in Figs.14. The PF-NM results used in these
plots were sourced from [34], where an identical analysis was performed. These
curves clearly show that the predictions from the LS-NM closely match those of
the sharp interface SAM, in contrast to the PF-NM, which predicts slower
transformation kinetics. Although the PF-NM predicts relatively lower
interfacial carbon concentrations, which theoretically should translate to
more driving pressure during the transformation, this does not lead to faster
interface kinetics. In the phase-field method, the migration of interface
kinetics is formulated using an effective or numerical interface mobility,
which depends on the diffuse interface thickness parameter and the physical
interface mobility. As highlighted in [34], during the phase-field simulation,
the effective interface thickness is not constant and is lower than the
initially imposed value. In contrast to the phase-field method, the level-set
numerical method relies solely on the physical interface mobility to dictate
interface migration kinetics. Since the migration of level-sets is
theoretically synonymous to that of the sharp interface method, only the
velocity field over the iso-zero contour of the level-set
($\varphi_{\alpha}=0$) is of prime importance. This differs from PF-NM, where
kinetics are considered across the entire diffuse interface. This sharp
interface characteristic of LS-NM in interface migration could perhaps reflect
the closer resemblance observed in its predictions compared to the fully
diffuse approach of PF-NM.

(a) 1D case

(b) 2D case

Figure 13: Comparison of the evolution of phase interface predicted by
different models

(a) 1D case

(b) 2D case

Figure 14: Comparison of the evolution of interfacial carbon concentration in
the austenite side as predicted by different models

However, Mecozzi et al. [34] demonstrate in the context of PF-NM that
enhancing mesh resolution and ensuring an adequate number of discretization
points within the diffuse interface indeed leads to better convergence in
their model’s predictions.

###  3.2 2D biphasic polycrystal cases

#### Binary alloy: A step towards a realistic large scale case

We now examine austenite decomposition in a binary steel alloy (Fe -
$0.02\text{\,}\mathrm{wt\%}$C) within a large-scale microstructure, emulating
the complexity often encountered in industrial settings. At atmospheric
pressure, the austenitization temperature for this alloy is approximately
$T_{A3}=$1175\text{\,}\mathrm{K}$$. We consider an initial polycrystalline
microstructure in a square domain of size $1\text{\,}\mathrm{mm}$, situated
slightly above the $T_{A3}$ temperature with
$T^{i}=$1176\text{\,}\mathrm{K}$$. At this stage, the microstructure is
composed entirely of austenite grains. Fig.15 illustrates the morphology of
the initial microstructure, featuring $1592$ austenite grains. The initial
grain morphology has been obtained by respecting the grain size distribution
(evaluated for each grain as the radius of a circle of equivalent area)
depicted in Fig.16 with an average grain radius of about
$14\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

Figure 15: Illustration of the initial austenitic microstructure with $1592$
grains (the grain interfaces are highlighted in black) Figure 16: Grain size
distribution of the initial large-scale microstructure

The microstructure is cooled at a constant rate of
$\dot{T}=$-10\text{\,}\mathrm{K}\text{$\cdot$}{\mathrm{s}}^{-1}$$ to
$T^{f}=$976\text{\,}\mathrm{K}$$. Following the cooling phase, the
microstructure is maintained at this temperature for an additional $80$
seconds to ensure complete transformation. It is then rapidly reheated to
$1100\text{\,}\mathrm{K}$ at a rate of
$20\text{\,}\mathrm{K}\text{$\cdot$}{\mathrm{s}}^{-1}$ to enhance grain
boundary mobility and thus promote grain growth effects. After the reheating,
the microstructure is held at $1100\text{\,}\mathrm{K}$ for the remaining
duration of the simulation until $t=$1000\text{\,}\mathrm{s}$$. The thermal
path corresponding to this scenario is illustrated in Fig.17.

Figure 17: Thermal path imposed during the transformation: the broken axis
representation indicates that the same temperature is maintained during that
period

The interface mobility and carbon diffusivity data are adopted from [35, 13,
36] with Arrhenius type law for temperature dependence:

$\begin{gathered}M=6\times
10^{17}\exp{\left(\frac{-140000}{RT}\right)},\quad\text{in}\quad$\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{4}\text{$\cdot$}{\mathrm{J}}^{-1}\text{$\cdot$}{\mathrm{s}}^{-1}$\\\
D_{\alpha}^{C}=2.2\times
10^{8}\exp\left(\frac{-122500}{RT}\right),\quad\text{in}\quad$\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{2}\text{$\cdot$}{\mathrm{s}}^{-1}$\\\ D_{\gamma}^{C}=1.5\times
10^{7}\exp\left(\frac{-142100}{RT}\right),\quad\text{in}\quad$\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{2}\text{$\cdot$}{\mathrm{s}}^{-1}$\end{gathered}.$ (57)

The value for the interface energy is taken following [37, 12], i.e.
$\sigma_{\gamma\alpha}=$1.0\text{\times}{10}^{-6}\text{\,}\mathrm{J}\text{$\cdot$}{\mathrm{mm}}^{-2}$$.
In this study, mobility and interfacial energy are assumed to be homogeneous,
isotropic, and identical across phase and grain interfaces of both phases. It
is however crucial to emphasize that the generalized kinetic description
employed in this numerical model enables the seamless integration of any
heterogeneity or anisotropy aspects if required.

Taking into account the available driving pressures for ferrite nucleation at
different temperatures below the austenitization temperature, and factoring in
the capillarity effects, the nucleation start temperature is estimated to be
$T_{N_{s}}=$1166\text{\,}\mathrm{K}$$. This temperature is chosen to increase
the probability of ferrite nuclei entering the growth regime of the
transformation. This delay in nucleation could also loosely represent the
notion of an incubation time. In this scenario, nucleation is configured to
occur solely at the triple junctions (grain corners), continuously, until
these sites are saturated, and spanning over a temperature range of $\delta
T_{N}=$40\text{\,}\mathrm{K}$$. Under these conditions, a total of $2207$
ferrite nuclei are generated.

The resolution time step is fixed at $\Delta t=$0.02\text{\,}\mathrm{s}$$. An
adaptive mesh with periodic remeshing is utilized. The local mesh size within
the diffuse phase interface is set as
$h_{min}=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, and the diffuse phase
interface thickness is set as $\eta=$8\text{\,}\mathrm{\SIUnitSymbolMicro m}$$
($\approx 20h_{min}$).

In Fig.18, the time evolution of the ferrite ($\alpha$) fraction is
illustrated. It is evident that the ferrite fraction converges to unity,
consistent with the ThermoCalc prediction indicating complete transformation
to a ferritic microstructure below $990\text{\,}\mathrm{K}$.

Figure 18: Time evolution of the ferrite fraction during the transformation

Figs.19 present the snapshots of the transforming microstructure at different
intervals during the cooling stage. The associated carbon redistribution maps
are illustrated in Figs.20. The continuous mode of nucleation results in a
dispersion in the eventual ferrite grain size distribution. Early-formed
nuclei gain an initial advantage in growth over subsequently appearing ones.
Moreover, the first few nuclei experience minimal hard impingement until new
neighboring nuclei gradually emerge. The presence of local nuclei clusters
also contribute to the dispersion in ferrite grain size due to the enhanced
hard impingement among the ferrite grains. Nucleation limited to grain corners
and the relatively high nucleation density due to rapid cooling tend to yield
globally finer equiaxed ferrite grains as evident from the Figs.19. The carbon
enrichment in the austenite phase is noticeable as the transformation
proceeds. The tiny dark spots observed in the Fig.20d is characterized by the
formation of minuscule trapped austenite islands or pockets. Due to the high
cooling rate, there is a rapid evolution resulting in a significant influx of
carbon atoms into those austenite grains, which are thoroughly surrounded by a
large number of ferrite grains. This substantial carbon enrichment
consequently reduces the driving pressure for phase transformation in these
areas, as local equilibrium is rapidly achieved. This results in trapped
austenite, characterized by visible solute gradients. Given sufficient time,
these solute gradients may diffuse into the bulk of the grains, enabling the
additional growth of ferrite grains that extend into these confined austenite
islands.

(a) $t=$1\text{\,}\mathrm{s}$,T=$1166\text{\,}\mathrm{K}$$

(b) $t=$3\text{\,}\mathrm{s}$,T=$1146\text{\,}\mathrm{K}$$

(c) $t=$5\text{\,}\mathrm{s}$,T=$1126\text{\,}\mathrm{K}$$

(d) $t=$20\text{\,}\mathrm{s}$,T=$976\text{\,}\mathrm{K}$$

Figure 19: Snapshots of austenite decomposition into ferrite in a large-scale
microstructure, at different instants till the end of cooling

(a) $t=$1\text{\,}\mathrm{s}$,T=$1166\text{\,}\mathrm{K}$$

(b) $t=$3\text{\,}\mathrm{s}$,T=$1146\text{\,}\mathrm{K}$$

(c) $t=$5\text{\,}\mathrm{s}$,T=$1126\text{\,}\mathrm{K}$$

(d) $t=$20\text{\,}\mathrm{s}$,T=$976\text{\,}\mathrm{K}$$

Figure 20: Snapshots of carbon evolution between the phases at different
instants in a large-scale microstructure till the conclusion of cooling Figure
21: Configuration obtained post complete transformation to ferrite at the end
of $t=$100\text{\,}\mathrm{s}$$: Phase distribution (left), Carbon
concentration field (right)

Figs.21 represent the state of the ferritic microstructure obtained after the
transformation stage (at $t=$100\text{\,}\mathrm{s}$$). The triple junctions
formed by the ferrite grains in Fig.21 are clearly far from the expected
equilibrium configuration ($120\text{\,}\mathrm{\SIUnitSymbolDegree}$ for
isotropic interfacial properties). Given the right activation, the system
could further reduce the excess free energy, primarily existing as surface
energy, through the grain growth (GG) phenomena. However, at temperatures as
low as $976\text{\,}\mathrm{K}$, the interface mobility remains low, and the
effects of GG are negligible within the timescales of interest. Therefore, to
initiate the necessary activation, the microstructure is reheated to a higher
temperature ($1100\text{\,}\mathrm{K}$) and maintained there to enhance the GG
effects. It should be remarked that the GG stage of the simulation is resolved
at a higher time step ($0.5\text{\,}\mathrm{s}$ in this case).

(a) $t=$200\text{\,}\mathrm{s}$,T=$1100\text{\,}\mathrm{K}$$

(b) $t=$400\text{\,}\mathrm{s}$,T=$1100\text{\,}\mathrm{K}$$

(c) $t=$700\text{\,}\mathrm{s}$,T=$1100\text{\,}\mathrm{K}$$

(d) $t=$1000\text{\,}\mathrm{s}$,T=$1100\text{\,}\mathrm{K}$$

Figure 22: Snapshots of grain growth effects observed in ferrite grains at
different instants after reheating

The influence of GG phenomenon on the resulting ferrite grain structure is
portrayed in Figs.22. The growth of larger grains at the expense of smaller
ones, aimed at minimizing the interfacial area, is clearly evident. Notably,
the triple junctions approach the equilibrium angles
($120\text{\,}\mathrm{\SIUnitSymbolDegree}$) characteristic of isotropic GG.
The differences in the ferrite grain size distributions obtained before the
onset of GG regime, and after GG effects at $t=$1000\text{\,}\mathrm{s}$$ are
showcased in Figs.23. As the GG effects become prominent, the grain size
distribution becomes more dispersed, with the average grain size gradually
increasing.

(a) At the beginning of GG stage ($t=$100\text{\,}\mathrm{s}$$)

(b) At $t=$1000\text{\,}\mathrm{s}$$

Figure 23: Ferrite grain size distribution during the grain growth (GG) stage

The initial and the final microstructures resulting from this thermal
treatment are depicted in Figs.24 using grain coloration. This scenario
produces a finer microstructure after both the transformation and the
subsequent grain growth. The current model does not account for the reverse
phase transformation during the reheating stage. However, since the driving
pressure description remains consistent, adapting the existing framework to
control the direction of the phase transformation should be relatively
straightforward. This could be a potential avenue for future research,
involving additional complexities such as nucleation for the reverse
transformation.

(a) Initial microstructure

(b) Final microstructure

Figure 24: Comparison of the microstructure morphologies before and after the
complete thermal treatment: represented by grain coloration

#### Ternary alloy: Influence of solute drags effects

We shall now examine diffusive phase transformation in a ternary alloy
(Fe-$0.1\text{\,}\mathrm{wt\%}$C-$0.5\text{\,}\mathrm{wt\%}$Mn) in the
presence of a substitutional Mn element. The influence of Mn on the
transformation kinetics are incorporated through solute drag aspects under
para-equilibrium hypothesis. For this case, a simple microstructure
illustrated in Fig.25 is considered with $317$ austenite grains.

Figure 25: Grain coloration adopted for the initial austenitic grain
morphology

For this alloy, the calculated austenitization temperature under the PE
hypothesis is approximately $1118\text{\,}\mathrm{K}$. The initial temperature
is set to this $T_{A3}^{PE}$ temperature. The microstructure is cooled to a
final temperature of $T^{f}=$910\text{\,}\mathrm{K}$$ at a cooling rate of
$-10\text{\,}\mathrm{K}\text{$\cdot$}{\mathrm{s}}^{-1}$. At end of cooling,
the microstructure is maintained at the final temperature for an additional
$30\text{\,}\mathrm{s}$ to provide ample time for phase transformation and the
associated solute redistribution. The interfacial energy, mobility,
diffusivity data, and most numerical parameters are fixed the same as in the
previous case. The time step is chosen as $\Delta t=$0.02\text{\,}\mathrm{s}$$
for this case. The ferrite nucleation has been limited to grain corners. For
simplicity, the number of ferrite nuclei stochastically introduced have been
restricted to $300$ for this case. The nucleation start temperature is set to
$T_{N_{s}}=$1112\text{\,}\mathrm{K}$$, and the range is assumed as $\delta
T_{N}=$30\text{\,}\mathrm{K}$$.

To model solute drag aspects, the simplified description provided by Cahn’s
model [23] is adopted as detailed in the second section. In practical
scenarios, the objective is to utilize the two parameters ($\alpha,\beta$)
introduced by Cahn’s solute drag pressure as fitting model parameters to agree
with the experimental results, rather than employing the analytical
expressions provided by Cahn in Eqs.(24) and (25). However, for demonstration
purposes in this context, we utilize the analytical expressions based on
literature [38] to establish some parameters in Eqs.(24) and (25), aiming for
as realistic a representation as possible. Figs.26 provide a visualization of
the Cahn’s parameters computed in this simulation, depicting their variation
with temperature. The graph illustrating the variation of $\beta_{C}^{2}$ in
Fig.26b indicates a rise in solute drag pressure at elevated temperatures,
gradually weakening as the temperature decreases.

(a) $\alpha_{C}(T)$

(b) $\beta_{C}^{2}(T)$

Figure 26: Evolution of Cahn’s solute drag parameters with temperature for the
chosen case

Fig.27 showcases the evolution of the ferrite phase fraction over time for two
scenarios: one without any solute drag effects, and the other incorporating
solute drag effects. It could be inferred that, enhanced drag effects at
higher temperatures limit the phase fraction initially. However, with
sufficient time, both the scenarios eventually converge to a similar ferrite
fraction value as depicted in Fig.27, since at lower temperatures the solute
drag effects are significantly reduced.

Figure 27: Evolution of ferrite fraction with time for the difference cases

Figs.28, and 29 compare the phase distribution and the solute distribution
during the transformation between the scenarios with and without the solute
drag effects. Both the scenarios globally result in equiaxed grains. However,
in the instance incorporating solute drag effects, the growth of ferrite
nuclei in the beginning at higher temperatures is sluggish due to strong drag
effects. Generally, this could limit the potential dispersion in grain size
distribution caused by continuous nucleation, resulting in a relatively
uniform distribution of grains. However, as evident in Figs.31, in this case,
the difference are minimal, possibly due to a smaller temperature range of
nucleation for the considered cooling rate. The comparison between the final
microstructure morphologies illustrated in Figs.30 demonstrates similar grain
characteristics, including a comparable mean ferrite grain size. The
noticeable distinctions primarily lie in the transformation kinetics between
the two scenarios for the assumed solute drag parameters in this illustration.

(a) Without solute drag effects

(b) With solute drag effects

Figure 28: Snapshots of phase evolution with and without the consideration of
solute drag effects

(a) Without solute drag effects

(b) With solute drag effects

Figure 29: Snapshots of carbon distribution with and without the consideration
of solute drag effects

(a) No solute drag effects

(b) With solute drag effects

Figure 30: Final phase distribution obtained with and without solute drag
effects

(a) Without solute drag effects

(b) With solute drag effects

Figure 31: Final grain size distributions obtained for the two cases

###  3.3 Simulation of other solid-state diffusive phenomena: Particle
coarsening

After a first-order diffusive phase transformation where second-phase
particles (SPPs) separate within a matrix phase, local equilibrium might be
achieved between the matrix and the SPPs. However, this state may not
represent the system’s minimum free energy configuration, particularly when
the matrix phase consists of a dispersed or multi-modal distribution of SPPs.
The excess free energy contributed by the interfacial energy of the of SPPs
renders the system thermodynamically unfavorable. Consequently, the system
attempts to liberate this excess free energy by dissolving smaller particles
in favor of larger ones through diffusional mass transport. The curvature of
SPPs alter the local equilibrium concentrations with the matrix phase due to
the Gibbs-Thomson effect. The multi-modal distribution of SPPs thus creates a
gradient in concentration which drives the diffusional mass transport between
a region of higher concentration and lower concentration. This phenomenon is
commonly referred to as particle coarsening or Ostwald ripening (OR) [39, 40].
In the absence of elastic misfits between the two phases, the kinetics of
particle coarsening are governed by the capillarity induced driving pressure
and the driving pressure due to the phase transformation. It should thus be
possible to simulate this phenomenon using the existing numerical framework.

To simulate this phenomenon, we’ve chosen to use a hypothetical scenario:
representing the SPPs as the ferrite phase ($\alpha$) and the matrix as the
austenite phase ($\gamma$). It’s crucial to highlight that despite this
hypothetical setup not mirroring the conventional appearance of Ostwald
ripening, the inherent kinetics and dynamics remain identical. The primary
objective is to showcase the numerical model’s capability in simulating
various diffusive solid-state phenomena, including particle coarsening,
utilizing the same kinetic framework. Therefore, our simulation involves an
austenite phase matrix within a $1\text{\,}\mathrm{mm}$-sized domain, wherein
the second phase consists of ferrite, forming a bimodal distribution of
particles, as illustrated in Fig.32. The initial configuration is represented
in Fig.33 through the grain coloration employed for this case.

Figure 32: Initial size distribution of second-phase particles (SPPs) Figure
33: Initial grain coloration of the matrix-SPPs setup

The initial condition is slightly offset ($T^{i}=$1162\text{\,}\mathrm{K}$$)
from the isothermal simulation temperature of $T=$1160\text{\,}\mathrm{K}$$
for the Fe-$0.02\text{\,}\mathrm{wt\%}$C alloy. This proximity ensures that
the system is nearly at local equilibrium for phase transformation and can
swiftly transition to the particle coarsening regime. The resolution time step
is fixed at $\Delta t=$0.5\text{\,}\mathrm{s}$$. The mobility and carbon
diffusivity data are consistent with previous cases. However, considering the
long time scale of this phenomenon, for illustrative purposes, the interfacial
energy is set to a higher value
($\sigma_{\gamma\alpha}=$2\text{\times}{10}^{-6}\text{\,}\mathrm{J}\text{$\cdot$}{\mathrm{mm}}^{-2}$$)
to accentuate the coarsening process. The other numerical parameters remain
the same as the previous cases. The simulation is run for a duration of
$10000\text{\,}\mathrm{s}$.

Figure 34: Snapshots of coarsening and dissolution of second phase particles
(SPPs)

Figs.34 present a series of snapshots illustrating the particle evolution at
different instants. The system is anticipated to reach local equilibrium after
the initial phase transformation at approximately $50\text{\,}\mathrm{s}$.
Thus, the particle coarsening regime is assumed to commence from
$t=$50\text{\,}\mathrm{s}$$. The expected coarsening of larger particles is
visibly apparent, as the smaller particles dissolve over time. It’s important
to note that physical coalescence aspects are not accounted for in this
simulation. Whenever two particles approach close proximity, the local
diffusional flux and the interaction with capillarity pressure seem to flatten
their interfaces.

Fig.35 demonstrates the particle size distribution obtained at the end of the
simulation. It seems evident that, in comparison to Fig.32, there has been an
increase in the frequency of larger particles over time as a result of
coarsening.

Figure 35: Final size distribution of second-phase particles (SPPs)

In the Fig.36, the kinetics of the particle size evolution are compared with
the modern Ostwald ripening theories [41]. These theories are based on a
generic expression:

$\bar{R}_{\alpha}^{3}-\bar{R}_{\alpha}^{3}(0)=kt,$ (58)

where $\bar{R}_{\alpha}$ is the mean particle size, $\bar{R}_{\alpha}^{3}(0)$
is the initial mean particle size, and $k$ is a parameter characterized based
on the concerned theory. It appears that the simulated mean particle size can
be fit into the aforementioned expression, resembling the quantitative
description provided by the OR theories. Therefore, the hypothetical scenario
we’ve adopted demonstrates kinetic characteristics akin to those of a typical
Ostwald ripening phenomenon. This serves as evidence of the numerical model’s
seamless potential to simulate a range of diffusive solid-state phenomena.

Figure 36: Comparison with modern Ostwald ripening theories

##  4 Conclusions and perspectives

A generalization of the global level-set numerical framework used to primarily
simulate diffusive solid-state phase transformation in metallic
polycrystalline materials is presented. The generalized framework incorporates
solute drag aspects into the transformation kinetics in the presence of any
substitutional solutes. A coupling with a thermodynamic database was
established to provide seamless extraction of necessary thermodynamic data for
the transformation. The level-set numerical model was benchmarked against a
state-of-the-art sharp interface semi-analytical model, alongside comparisons
with predictions from a corresponding phase-field numerical model. The
potential of the proposed numerical framework to replicate the phase
transformation behavior in complex large-scale microstructures with thousands
of grains, relevant to industrial settings was showcased. An illustrative case
considering solute drag aspects in a ternary alloy was presented, outlining
the model’s ability to reproduce sluggish transformation kinetics. The
adjustable solute drag parameters offer additional flexibility to control the
transformation kinetics, to be in better agreement with the experimental
transformation curves. The numerical model’s versatility to simulate other
diffusive solid-state phenomena such as particle coarsening, without any
modifications to the existing formulation was demonstrated. The generalized
nature of the kinetic framework alongside the adaptive meshing capabilities
offer diverse scope to seamlessly integrate other complex evolution aspects
into the model, including simulations in 3D.

The positive outcomes of this work are expected to guide future research,
potentially expanding the range of modeling approaches available alongside the
established phase-field method in the realm of diffusive solid-state phase
transformation. A comprehensive experimental analysis to validate the proposed
numerical model is planned. Such an analysis would enable improvements to the
current nucleation model, potentially refining it based on experimental
observations. In the proposed numerical framework, solute diffusion is assumed
homogeneous within the bulk of a grain and across the grains of the same
phase. Consequently, exploring the impact of short-circuit diffusion effects
during phase transformation emerges as a potential avenue for future work.
Integrating the effects of elastic/plastic accommodation due to volume misfits
at the interphase boundaries into the current kinetic framework could further
enhance our understanding of interface kinetics, providing a more
comprehensive description of the physics involved. The global objective of
this work was to establish a kinetic framework capable of encompassing various
evolution aspects simultaneously, notably in a high plastic deformation
context. The existing numerical tools and concepts devised for
recrystallization could be integrated into the current framework, enabling the
simulation of dynamic recrystallization in a multiphase context and
potentially considering, at the mesoscopic scale, the modeling of realistic
industrial thermomechanical paths for multiphase materials. Another
perspective of these developments concerns the potential
improvement/optimization of state-of-the-art equations used in macroscale code
to describe type II residual-stresses due to phase transformations mainly
during quenching and surface hardening treatments.

## Acknowledgements

The authors thank ArcelorMittal, Aperam, Ascometal, Aubert & Duval, CEA,
Constellium, Framatome, and Safran companies and the ANR for their financial
support through the DIGIMU consortium and RealIMotion ANR industrial Chair
(Grant No. ANR-22-CHIN-0003).

## References

  * [1] Y. Jin, A. Khachaturyan, Atomic density function theory and modeling of microstructure evolution at the atomic scale, Journal of applied physics 100 (1) (2006) 013519.
  * [2] Y. Mishin, M. Asta, J. Li, Atomistic modeling of interfaces and their impact on microstructure and properties, Acta Materialia 58 (4) (2010) 1117–1151.
  * [3] M. Biglari, Simulation of the austenite-ferrite transformation; effect of applied stress (2013).
  * [4] K. Janssens, D. Raabe, E. Kozeschnik, M. Miodownik, B. Nestler, Computational materials engineering: an introduction to microstructure evolution, Academic Press, 2010.
  * [5] M. Bernacki, R. Logé, T. Coupez, Level set framework for the finite-element modelling of recrystallization and grain growth in polycrystalline materials, Scripta Materialia 64 (6) (2011) 525–528.
  * [6] J. Fausty, N. Bozzolo, D. Pino Muñoz, M. Bernacki, A novel level-set finite element formulation for grain growth with heterogeneous grain boundary energies, Materials & Design 160 (2018) 578–590.  doi:https://doi.org/10.1016/j.matdes.2018.09.050.
URL https://www.sciencedirect.com/science/article/pii/S0264127518307524

  * [7] K. Alvarado, S. Florez, B. Flipon, N. Bozzolo, M. Bernacki, A level set approach to simulate grain growth with an evolving population of second phase particles, Modelling and Simulation in Materials Science and Engineering 29 (3) (2021) 035009\.  doi:10.1088/1361-651X/abe0a7.
URL https://hal.archives-ouvertes.fr/hal-03463482

  * [8] L. Maire, Full field and mean field modeling of dynamic and post-dynamic recrystallization in 3D – Application to 304L steel, Theses, Université Paris sciences et lettres (Nov. 2018).
URL https://pastel.archives-ouvertes.fr/tel-02418810

  * [9] M. Hillert, J. Ågren, On the definitions of paraequilibrium and orthoequilibrium, Scripta Materialia 50 (5) (2004) 697–699.
  * [10] D.-H. Yeon, P.-R. Cha, J.-K. Yoon, A phase field study for ferrite–austenite transitions under paraequilibrium, Scripta Materialia 45 (2001) 661–668.  doi:10.1016/S1359-6462(01)01077-6.
  * [11] G. Pariser, P. Schaffnit, I. Steinbach, W. Bleck, Simulation of the $\gamma-\alpha$ transformation using the phase-field method, Steel Research 72 (2001) 354–360.  doi:10.1002/srin.200100130.
  * [12] C. Huang, D. Browne, S. McFadden, A phase-field simulation of austenite to ferrite transformation kinetics in low carbon steels, Acta Materialia 54 (1) (2006) 11–21.  doi:10.1016/j.actamat.2005.08.033.
  * [13] M. Mecozzi, Phase field modelling of the austenite to ferrite transformation in steels, Ph.D. thesis, Delf University of Technology (2006).
  * [14] M. Militzer, M. G. Mecozzi, J. Sietsma, S. Zwaag, Three-dimensional phase field modelling of the austenite-to-ferrite transformation, Acta Materialia 54 (2006) 3961–3972.  doi:10.1016/j.actamat.2006.04.029.
  * [15] K. Bzowski, L. Rauch, M. Pietrzyk, Application of statistical representation of the microstructure to modeling of phase transformations in dp steels by solution of the diffusion equation, Procedia Manufacturing 15 (2018) 1847–1855.
  * [16] N. Chandrappa, M. Bernacki, A level-set formulation to simulate diffusive solid/solid phase transformation in polycrystalline metallic materials-application to austenite decomposition in steels, Computational Materials Science 216 (2023) 111840.
  * [17] B. Scholtes, M. Shakoor, A. Settefrati, P.-O. Bouchard, N. Bozzolo, M. Bernacki, New finite element developments for the full field modeling of microstructural evolutions using the level-set method, Computational Materials Science 109 (2015) 388–398.  doi:10.1016/j.commatsci.2015.07.042.
  * [18] B. Scholtes, Development of an efficient level set framework for the full field modeling recrystallization in 3d, Ph.D. thesis, Université Paris sciences et lettres (2016).
  * [19] C. Hutchinson, H. Zurob, Y. Brechet, The growth of ferrite in fe-cx alloys: The role of thermodynamics, diffusion, and interfacial conditions, Metallurgical and Materials Transactions A 37 (2006) 1711–1720.
  * [20] G. Purdy, J. Ågren, A. Borgenstam, Y. Bréchet, M. Enomoto, T. Furuhara, E. Gamsjager, M. Gouné, M. Hillert, C. Hutchinson, et al., Alemi: A ten-year history of discussions of alloying-element interactions with migrating interfaces, Metallurgical and Materials Transactions A 42 (2011) 3703–3718.
  * [21] Thermo-calc software, http://www.thermocalc.se/.
  * [22] J. Christian, The theory of transformations in metals and alloys, Newnes, 2002.
  * [23] J. Cahn, The impurity-drag effect in grain boundary motion, Acta metallurgica 10 (9) (1962) 789–798.
  * [24] N. Chandrappa, A global level-set based full-field numerical framework to simulate diffusive solid-state phase transformation in polycrystalline metallic materials, Ph.D. thesis, Université Paris sciences et lettres (2024).
  * [25] J. Furstoss, C. Petit, A. Tommasi, C. Ganino, D. Muñoz, M. Bernacki, On the role of solute drag in reconciling laboratory and natural constraints on olivine grain growth kinetics, Geophysical Journal International 224 (2) (2021) 1360–1370.
  * [26] M. Crouzeix, Une méthode multipas implicite-explicite pour l’approximation des équations d’évolution paraboliques, Numerische Mathematik 35 (3) (1980) 257–276.
  * [27] U. Ascher, S. Ruuth, R. Spiteri, Implicit-explicit runge-kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics 25 (2-3) (1997) 151–167.
  * [28] S. Offerman, N. Van Dijk, J. Sietsma, S. Grigull, E. Lauridsen, L. Margulies, H. Poulsen, M. Rekveldt, S. Van der Zwaag, Grain nucleation and growth during phase transformations, Science 298 (5595) (2002) 1003–1005.
  * [29] M. Mecozzi, M. Militzer, J. Sietsma, S. Van der Zwaag, The role of nucleation behavior in phase-field simulations of the austenite to ferrite transformation, Metallurgical and Materials Transactions A 39 (2008) 1237–1247.
  * [30] A. Brooks, T. Hughes, Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations, Computer Methods in Applied Mechanics and Engineering 32 (1) (1982) 199 – 259.  doi:https://doi.org/10.1016/0045-7825(82)90071-8.
  * [31] B. Merriman, J. Bence, S. Osher, Motion of multiple junctions: a level set approach, Journal of Computational Physics 112 (1994) 334–363.
  * [32] M. Shakoor, B. Scholtes, P.-O. Bouchard, M. Bernacki, An efficient and parallel level set reinitialization method–application to micromechanics and microstructural evolutions, Applied Mathematical Modelling 39 (23-24) (2015) 7291–7302.
  * [33] C. Bos, J. Sietsma, A mixed-mode model for partitioning phase transformations, Scripta materialia 57 (12) (2007) 1085–1088.
  * [34] M. Mecozzi, J. Eiken, M. Apel, J. Sietsma, Quantitative comparison of the phase-transformation kinetics in a sharp-interface and a phase-field model, Computational Materials Science 50 (6) (2011) 1846–1853.
  * [35] G. Krielaart, S. Van der Zwaag, Kinetics of $\gamma$→ $\alpha$ phase transformation in fe-mn alloys containing low manganese, Materials science and technology 14 (1) (1998) 10–18.
  * [36] D. Lide, CRC handbook of chemistry and physics, Vol. 85, CRC press, 2004.
  * [37] I. Loginova, J. Odqvist, G. Amberg, J. Ågren, The phase-field approach and solute drag modeling of the transition to massive $\gamma$→ $\alpha$ transformation in binary fe-c alloys, Acta Materialia 51 (5) (2003) 1327–1339.
  * [38] F. Fazeli, M. Militzer, Application of solute drag theory to model ferrite formation in multiphase steels, Metallurgical and Materials Transactions A 36 (2005) 1395–1405.
  * [39] W. Ostwald, Blocking of ostwald ripening allowing long-term stabilization, Phys. Chem 37 (1901) 385.
  * [40] P. Voorhees, The theory of ostwald ripening, Journal of Statistical Physics 38 (1-2) (1985) 231–252, copyright: Copyright 2007 Elsevier B.V., All rights reserved.  doi:10.1007/BF01017860.
  * [41] A. Baldan, Review progress in ostwald ripening theories and their applications to nickel-base superalloys part i: Ostwald ripening theories, Journal of materials science 37 (2002) 2171–2202.