Multi-Phase-Field Model for Surface/Phase-Boundary Diffusion

TL;DR

This paper extends the multi-phase-field model to include capillarity-driven diffusion along surfaces and phase boundaries, enabling the study of phenomena like thermal grooving and nano-cluster annealing with new insights into transient surface structures.

Contribution

The authors develop a generalized multi-phase-field model incorporating capillarity effects for surface and phase-boundary diffusion, applied to thermal grooving and nano-cluster annealing.

Findings

Observation of decaying oscillations in thermal grooving profiles.

Identification of transient crater-like structures during nano-cluster annealing.

Model successfully captures capillarity-driven diffusion phenomena.

Abstract

The multi-phase-field approach is generalized to treat capillarity-driven diffusion parallel to the surfaces and phase-boundaries, i.e. the boundaries between a condensed phase and its vapor and the boundaries between two or multiple condensed phases. The effect of capillarity is modeled via curvature-dependence of the chemical potential whose gradient gives rise to diffusion. The model is used to study thermal grooving on the surface of a polycrystalline body. Decaying oscillations of the surface profile during thermal grooving, postulated by Hillert long ago but reported only in few studies so far, are observed and discussed. Furthermore, annealing of multi-nano-clusters on a deformable free surface is investigated using the proposed model. Results of these simulations suggest that the characteristic crater-like structure with an elevated perimeter, observed in recent experiments, is a transient non-equilibrium state during the annealing process.

Tables

Table 1. Table 1: Polynomial coefficients of Eq. ( 23 ).

	$a_0$	$-7.803$	$\times {10}^{-1}$
	$a_1$	$10.000$	$\times {10}^{-1}$
	$a_2$	$-2.886$	$\times {10}^{-1}$
	$a_4$	$8.130$	$\times {10}^{-3}$
	$a_6$	$-2.004$	$\times {10}^{-4}$
	$a_8$	$3.625$	$\times {10}^{-6}$
	$a_{10}$	$-4.969$	$\times {10}^{-8}$
	$a_{12}$	$5.339$	$\times {10}^{-10}$
	$a_{14}$	$-4.655$	$\times {10}^{-12}$

Table 2. Table 2: The dependency of the contact angle between the surface of two grains on the grain boundary energy, σ GB subscript 𝜎 GB \sigma_{\text{GB}} . The surface energy is kept constant at σ surface = 1 subscript 𝜎 surface 1 \sigma_{\text{surface}}=1 J/m 2 . Simulation results are listed together with the prediction of the von Neumann’s relation, Eq. ( 24 ).

	${\sigma }_{\text{GB}}$ (J/m2)	${\theta }_{\text{von Neumann}}{(}^{\circ })$	${\theta }_{\text{simulation}}{(}^{\circ })$
	0.25	97.2	96.4
	0.50	104.5	103.4
	0.75	112.0	109.6
	1.00	120.0	118.0
	1.25	128.7	125.3
	1.50	138.6	134.7

Equations

$$\ln \left(\frac{p_V}{p_0}\right)=\frac{\sigma V_m}{k_{B}T}\kappa$$

$$\mu \approx -k_{B}T\ln \left(\frac{k_{B}T}{p_V{\lambda }^3}\right)$$

$$\mu =-k_{B}T\ln \left(\frac{k_{B}T}{p_0{\lambda }^3}\right)+\sigma V_m\kappa$$

$$\mathbf{j}=-\tilde{\mu }c{\nabla }_s\mu =-\tilde{\mu }\sigma c{V}_m{\nabla }_s\kappa$$

$$\dot{c}=\frac{D_s\sigma c{V}_m}{k_{B}T}{\nabla }_s^2\kappa$$

$$v_n=\frac{D_s\sigma c{V}_m^2}{k_{B}T}{\nabla }_s^2\kappa$$

$$v_n=-\sigma M\kappa$$

$$v_n=\sigma \left(\tilde{M}{\nabla }_s^2\kappa -M\kappa \right)$$

$$\cos \theta =\frac{{\sigma }_{\beta \gamma }-{\sigma }_{\beta \alpha }}{{\sigma }_{\gamma \alpha }}$$

$$\frac{{\sigma }_{\beta \gamma }}{\sin {\theta }_{\alpha }}=\frac{{\sigma }_{\alpha \gamma }}{\sin {\theta }_{\beta }}=\frac{{\sigma }_{\alpha \beta }}{\sin {\theta }_{\gamma }}.$$

$$F\left(\left\{\varphi \right\}\right)=\sum _{\alpha =1}^N\sum _{\beta =\alpha +1}^N{\int }_{\Omega }{f}_{\alpha \beta }\left({\varphi }_{\alpha },{\varphi }_{\beta }\right)dV$$

$$f_{\alpha \beta }=\frac{4{\sigma }_{\alpha \beta }}{\eta }\left(-\frac{{\eta }^2}{{\pi }^2}\nabla {\varphi }_{\alpha }\cdot \nabla {\varphi }_{\beta }+\left|{\varphi }_{\alpha }{\varphi }_{\beta }\right|\right)$$

$${\dot{\varphi }}_{\alpha }=-\frac{1}{N}\sum _{\beta =1}^N\nabla \cdot {\tilde{\mathbf{M}}}_{\alpha \beta }\nabla \left(\frac{\delta }{\delta {\varphi }_{\alpha }}-\frac{\delta }{\delta {\varphi }_{\beta }}\right)F=-\frac{1}{N}\sum _{\beta =1}^N\nabla \cdot {\tilde{\mathbf{M}}}_{\alpha \beta }\nabla {\psi }_{\alpha \beta }$$

$${\psi }_{\alpha \beta }=\left(\frac{\delta }{\delta {\varphi }_{\alpha }}-\frac{\delta }{\delta {\varphi }_{\beta }}\right)F=I_{\alpha \beta }-I_{\beta \alpha }+\sum _{\gamma =1,\gamma \ne \alpha ,\gamma \ne \beta }^N\left(I_{\alpha \gamma }-I_{\beta \gamma }\right)$$

$$I_{\alpha \beta }=\frac{\partial f_{\alpha \beta }}{\partial {\varphi }_{\alpha }}-\nabla \frac{\partial f_{\alpha \beta }}{\partial \nabla {\varphi }_{\alpha }}$$

$$=\frac{4{\sigma }_{\alpha \beta }}{\eta }\left(\frac{{\eta }^2}{{\pi }^2}{\nabla }^2{\varphi }_{\beta }+\frac{\partial \left|{\varphi }_{\alpha }{\varphi }_{\beta }\right|}{\partial {\varphi }_{\alpha }}\right)$$

$${\mathbf{n}}_{\alpha \beta }=\frac{{\varphi }_{\alpha }\nabla {\varphi }_{\beta }-{\varphi }_{\beta }\nabla {\varphi }_{\alpha }}{\Vert {\varphi }_{\alpha }\nabla {\varphi }_{\beta }-{\varphi }_{\beta }\nabla {\varphi }_{\alpha }\Vert }$$

$${\tilde{\mathbf{M}}}_{\mathbf{\alpha \beta }}={\tilde{M}}_{\alpha \beta }\left(1-a_{\alpha \beta }{\mathbf{n}}_{\alpha \beta }{\mathbf{n}}_{\alpha \beta }\right)$$

$$a_{\alpha \beta }={\mathrm{limit}}_1\left(\frac{{\eta }^2}{{\pi }^2}\frac{\left|\nabla {\varphi }_{\alpha }\cdot \nabla {\varphi }_{\beta }\right|}{\left|{\varphi }_{\alpha }\right|\left|{\varphi }_{\beta }\right|}\right)$$

$${\mathrm{limit}}_1\left(x\right)=\{\begin{array}{ll}1& x>1\\ x& \text{else}.\end{array}$$

$${\tilde{\mathbf{M}}}_{\alpha \beta }={\tilde{M}}_{\alpha \beta }1$$

$${\varphi }_{\text{eq}}=\{\begin{array}{ll}1& \text{if }r<R-\frac{\eta }{2}\\ \frac{1}{2}-\frac{1}{2}\sin \left(\frac{\pi \left(r-R\right)}{\eta }\right)& \text{if }R-\frac{\eta }{2}\le r\le R+\frac{\eta }{2}\\ 0& \text{else},\end{array}$$

$$y_M\left(x,t\right)=m(Bt{)}^{\frac{1}{4}}\sum _{n=0}{a}_n{\left[\frac{x}{(Bt{)}^{\frac{1}{4}}}\right]}^n+\mathcal{O}\left(m^2\right),$$

$$\theta =\pi -\arctan \left(\frac{\sqrt{4{\sigma }_{\text{Surface}}^2-{\sigma }_{\text{GB}}^2}}{{\sigma }_{\text{GB}}}\right)$$

$$\psi =\frac{8{\sigma }_{\alpha \beta }}{\eta }\left[\frac{{\eta }^2}{{\pi }^2}{\nabla }^2\varphi +2\varphi -1\right]$$

$${\mathbf{F}}_1=\frac{m}{\sqrt{2}}{\left(-1,-1,0\right)}^T$$

$${\mathbf{F}}_2=\frac{m}{\sqrt{2}}{\left(1,-1,0\right)}^T$$

$${\mathbf{F}}_3=\frac{m}{\sqrt{2}}{\left(0,1,-1\right)}^T$$

$${\mathbf{F}}_4=\frac{m}{\sqrt{2}}{\left(0,1,1\right)}^T$$

Full text

I Introduction

Diffusion in materials results from random motion of their constituents (atoms/molecules/voids). Aside from self-diffusion of atoms or molecules within a homogeneous body in thermal equilibrium, diffusion is one of the most frequent transport mechanisms driven by weak deviations from the equilibrium state, such as a gradient in concentration or, more generally, chemical potential. This latter aspect becomes particularly relevant at surfaces/interfaces with spatially varying curvature, which, due to Gibbs-Thompson relation, exhibit a gradient in chemical potential along the surface/phase boundary. It is noteworthy that, while the capillarity driven diffusion at free surfaces is rather well investigated (see e.g. Herring (1999); Mullins (1957); Gugenberger et al. (2008)), much less attention has been payed to diffusion at the boundaries between two condensed phases. This work focuses on this issue using the phase-field method. In order to distinguish this type of diffusive transport from grain-boundary diffusion, in which boundaries are faster tracks for diffusion of solute atoms, the term ‘phase-boundary diffusion’ is introduced here to express the capillarity-driven diffusion in boundaries. A first investigation of a conserved phase-field model has been carried out by Caginalp in Caginalp (1988, 1990). Later, a conserved phase-field method has been used to simulate the evolution of elastically stressed films Rätz et al. (2006); Yeon et al. (2006). These models, however, are limited to surface diffusion between two phases. A comparison of the available phase-field models for surface diffusion is presented in Gugenberger et al. (2008). In multiphase multi-grain systems, diffusion in phase-boundaries requires a treatment of multiple conserved phases at triple and higher-order junctions.

In this work, we present a conserved multi-phase-field method which is able to describe the capillarity-driven surface/phase-boundary diffusion in such complex topologies, where multiple junctions may be present. The current method is based on the works reviewed by Gugenberger et al. Gugenberger et al. (2008) on the two-phase surface diffusion models and the non-conserved multi-phase-field model Steinbach (2009, 2013). For the sake of simplicity, here we consider immiscible phases without bulk diffusion and phase transformation, and solely focus on the effects of surface and phase-boundary diffusion. The present approach can be used to investigate complex systems for which it is not possible to obtain an analytic solution or which are simply too large to be investigated by atomistic methods. One of these possible scenarios is the solid-phase sintering of many particles where we can reach beyond the time and length scales of atomistic simulations. As an example, we apply our model to study surface and phase-boundary diffusion of nano-clusters during annealing on a free surface. This is similar to the experimental work presented in Köhler et al. (2013). Furthermore, the construction of the current model allows the combination of surface/boundary diffusion with dissipative interface diffusion Steinbach et al. (2012); Zhang and Steinbach (2012), elasticity Steinbach and Apel (2006); Mosler et al. (2014) and chemo-mechanical coupling Kamachali et al. (2014); Darvishi Kamachali and Schwarze (2017) effect in the same way that they are applied in non-conserved multi-phase-field methods.

The paper is organized as follows. In Sec. I.1, a short introduction into the thermodynamics of phase-boundary diffusion is given. The equilibrium solution at triple junctions is discussed in Secs. I.2 and I.3. The phase-boundary diffusion model is presented in Sec. II. Section III compiles the obtained simulation results. After a brief discussion of the equilibrium interface profile in Sec. III.1, thermal grooving is addressed in Sec. III.2. The late time profiles obtained from these simulations are then used in Sec. III.3 to check whether the method correctly recovers the von Neumann’s triangle relation. Section III.4 highlights that the equilibrium behavior at the junctions is strongly dependent on the conservation constraints. The evolution of multiple nano-clusters on a deformable surface is discussed in Sec. III.5. A summary compiles the main findings in Sec. IV.

I.1 Thermodynamics of Surface Diffusion

We follow an approach similar to Mullins Mullins (1957) and consider a body of single species surrounded by its vapor phase at constant temperature. The surface is assumed to be in equilibrium with a finite curvature, $\kappa$, and the interface energy, $\sigma$. The curvature is here considered as the sum of the principal curvatures, $\kappa_{1}$ and $\kappa_{2}$, with $\kappa=\kappa_{1}+\kappa_{2}$. The vapor pressure, $p_{V}$, is given by (assuming ideal gas and incompressible solid),

[TABLE]

, where $p_{0}$ is the vapor pressure for the planar surface, $k_{B}$ is the Boltzmann constant, $V_{m}$ is the molecular volume of the solid phase and $T$ is the temperature. The chemical potential, $\mu$, of an ideal gas can be approximated as

[TABLE]

. In this relation, $\lambda$ is the thermal de Broglie wave length. Inserting the vapor pressure, $p_{V}$, from Eq. (1) into Eq. (2) results in an expression for the dependency of the chemical potential on curvature,

[TABLE]

. Considering a body of arbitrary curvature, this leads to a flux of particles, $\boldsymbol{\mathbf{j}}$, from regions of high curvature to regions of low curvature,

[TABLE]

, where $\tilde{\mu}$ is the mobility of the particles and $c$ the concentration of particles per unit area. Here we treat the surface, $s$, as an isosurface of the concentration field with a constant value of $c$. In a more general concept of multiple species, the interface energy may depend on the composition and thus the concentrations may vary alongside the surface as well. Equation (4) describes the flux of particles parallel to the surface. The local increase of concentration of particles per unit area, $\dot{c}$, is obtained by taking the surface divergence of $\boldsymbol{\mathbf{j}}$ so that

[TABLE]

, where $D_{s}=\tilde{\mu}k_{B}T$ is the surface diffusion coefficient. By multiplying $\dot{c}$ with the molecular volume, the surface normal velocity can be obtained,

[TABLE]

. Equation (6) is the basic equation which characterizes the dynamics of surface diffusion. In non-conserved processes, e.g. phase transition, evaporation and condensation, the surface velocity is often described by a linear frictional ansatz

[TABLE]

, where $M$ is a mobility coefficient (see Mullins (1957)). Usually, conserved surface diffusion and non-conserved phase transition processes are studied separately. In principle, however, one can think of situations in which a combination of both effects determines the surface normal velocity,

[TABLE]

, in which a second mobility coefficient has been introduced, $\tilde{M}$. Equation (8) states that the significance of the two contributions is determined by the coefficients $\tilde{M}$ and $M$, as well as the length scale of the process. For smaller scales, i.e. larger curvatures, the effect of surface diffusion becomes increasingly more important. It is to be noted that the dynamics of surface diffusion, Eq. (6), for the solid/vapor interface is derived under the assumption of an ideal gas. This is a good approximation for a capillary driven diffusion at a solid/vapor interface at elevated temperatures.

I.2 Not Deformable Surface

A common example of three phases in contact is a liquid droplet, $\alpha$, on a flat non-deformable solid surface, $\beta$, surrounded by a gas phase, $\gamma$. By minimizing the surface energy, one can obtain Young’s law for the contact angle,

[TABLE]

. One can see that Eq. (9) does not have a solution for all possible values of the interface energies. The droplet can either detach form the surface, $\frac{\sigma_{\beta\gamma}-\sigma_{\beta\alpha}}{\sigma_{\gamma\alpha}}\leq-1$, or can completely wet the surface, $\frac{\sigma_{\beta\gamma}-\sigma_{\beta\alpha}}{\sigma_{\gamma\alpha}}\geq 1$. More frequently, however, one has to deal with the case of partial wetting, described by the intermediate values, $-1<\frac{\sigma_{\beta\gamma}-\sigma_{\beta\alpha}}{\sigma_{\gamma\alpha}}<1$.

I.3 Deformable Surface

The above consideration becomes more complex for three deformable phases in contact. At a triple point in equilibrium, the sum of all forces acting on the point is zero and the force vectors form a triangle (see Fig. 1b). Then the law of sines can be applied to express the relation between the magnitude of the forces and their angles, $\theta_{\alpha}$, $\theta_{\beta}$ and $\theta_{\gamma}$, to each other,

[TABLE]

This construction for the interface energies is referred to as the von Neumann’s triangle relation (see e.g. Neumann (1894); Lester (1961); Shanahan (1987)). A comprehensive derivation of Eq. (10) can be found in section 8.2 of the seminal textbook on molecular theory of capillarity by J.S. Rowlinson and B. Widom Rowlinson and Widom (1982). In the case of a junction between more than three phases, the situation becomes more complex. The general approach, however, is the same and is based on the idea of interface energy minimization.

II The Phase-Boundary Diffusion Model

In the following, we present a model for phase-boundary diffusion which is driven by the minimization of the local interface energy. The present approach combines an existing non-conserved multi-phase-field model with the conserved phase-field models for surface diffusion at two-phase-boundaries. The description of the referred multi-phase-field model is detailed in two articles Steinbach (2009) and Steinbach (2013). Furthermore, a review of surface diffusion models is given in Gugenberger et al. (2008). In comparison to existing phase-field models for surface diffusion, the present model can describe the capillary driven diffusion of multiple phases. Since the phase-boundary diffusion is the only process considered here, the volume of each individual phase is conserved.

II.1 The Interface Free Energy Density Functional

It is convenient to start with a suitable description of the free energy, $F$. The free energy depends on the spacial configuration of $N$ phase-fields $\phi_{\alpha}\left(\boldsymbol{\mathbf{x}},t\right)$ with $\alpha\in\left[1,N\right]$, the spatial position $\boldsymbol{\mathbf{x}}$ and the time $t$. Each phase-field $\phi_{\alpha}$ indicates the location of a thermodynamic phase $\alpha$, in a way that $\phi_{\alpha}=1$ means phase $\alpha$ is present and $\phi_{\alpha}=0$ that it is not. Additionally, a sum constrained for the phase-fields is used, $\sum_{\alpha=1}^{N}\phi_{\alpha}=1$. In the case of phase-boundary diffusion, the phase-field value may be interpreted as an indicator function of a thermodynamic phase with a fixed composition and density.

Following the multi-phase-field approach in Steinbach (2009, 2013), the total free energy can be written as a functional $F\left(\left\{\phi_{\alpha}\right\}\right)$. This free energy functional is defined through the free energy densities, $f_{\alpha\beta}$, of the interface between the phases $\alpha$ and $\beta$. The sum of all interface energy densities is integrated over the volume of interest, $\Omega$,

[TABLE]

. One may add other energy densities into the concept, in order to include other physical effects. The energy densities of the interfaces are proposed in a pair-wise manner between existing phase-fields,

[TABLE]

. Here, $\eta$ defines the width of the transition region between two phase-fields, e.g., $\phi_{\alpha}$ and $\phi_{\beta}$. $\left|\phi_{\alpha}\phi_{\beta}\right|$ is called double obstacle potential and restricts the value of the phase-field to the interval $\left[0,1\right]$. Figure 2 shows a plot of the double obstacle potential which can be simplified to $\left|\phi\left(1-\phi\right)\right|$ with $\phi=\phi_{\alpha}=1-\phi_{\beta}$ in the case of only two phase-fields. As it can be seen from Fig. 2, the use of absolute value is necessary to ensure that $\phi=0$ and $\phi=1$ correspond to two equilibrium solutions.

II.2 The Dynamics of Phase-Boundary Diffusion

The dynamic equations for phase-boundary diffusion can now be derived by using the variational derivative of the free energy functional with respect to the phase-field,

[TABLE]

. Here, $\boldsymbol{\tilde{M}}_{\alpha\beta}=\boldsymbol{\tilde{M}}_{\beta\alpha}$ is a symmetric mobility tensor of the interface between the phases $\alpha$ and $\beta$. It is of fourth-rank with $\tilde{M}_{ji\alpha\beta}=\tilde{M}_{ij\alpha\beta}$, where $i$ and $j$ indicate Cartesian coordinates, $i,j\in\{x,y,z\}$. Equation (13) ensures that each phase-field, $\phi_{\alpha}$, is conserved. Moreover, since $\sum^{N}_{\alpha=1}\dot{\phi}_{\alpha}=0$, it follows that the sum rule, $\sum_{\alpha=1}^{N}\phi_{\alpha}(t)=1$, is fulfilled at all times, $t$, provided that it is valid at $t=0$. In Eq. (13), $\psi_{\alpha\beta}$ is a generalized diffusion potential, defined via

[TABLE]

. The function $I_{\alpha\beta}$ is a shorthand notation for the functional derivative of the free interface density $f_{\alpha\beta}$, {dgroup}

[TABLE]

.

[TABLE]

.

It should be remarked here that $I_{\alpha\beta}\neq I_{\beta\alpha}$. In general, the interface energy, $\sigma_{\alpha\beta}\left(\boldsymbol{\mathbf{n}}_{\alpha\beta}\right)$, is a function of the interface orientation which can be determined by using the interface normal vector,

[TABLE]

. This leads to additional terms in Eq. (16) such as the Herring torque Herring (1999). For the sake of simplicity, however, this study is restricted to the case of a constant interface energy. A comment is necessary here regarding the use of a three dimensional $\nabla$ operator in Eq. (13) instead of the surface gradient $\nabla_{s}$ (see Eq. (5)). This is a consequence of the fact that, in diffuse interface approaches such as the present multi-phase-field model, a strictly two dimensional interface does not exist. Rather, it is replaced by an interface layer with a finite thickness, $\eta$. In the limit of a sharp interface $(\eta\rightarrow 0)$, however, Eq. (13) approaches Eq. (5). This is ensured by an appropriate choice of the mobility tensor, $\boldsymbol{\tilde{M}}_{\alpha\beta}$, in such a way as to effectively restrict the diffusion flux to the directions tangential to the interface. Diffusion along the normal direction is only allowed in order to keep the phase-field profile stable and vanishes in the limit of a sharp interface. An asymptotic analysis of multiple surface diffusion models is performed in Gugenberger et al. (2008). Here we set,

[TABLE]

. In this way, the tensorial mobility, $\boldsymbol{\tilde{M}}_{\alpha\beta}$, is restricted to the tangential plane of the interface by using the interface normal $\boldsymbol{\mathbf{n}}_{\alpha\beta}$. In Eq. (18), $\tilde{M}_{\alpha\beta}$ is the scalar magnitude of $\boldsymbol{\tilde{M}}_{\alpha\beta}$ and $\boldsymbol{1}$ is the unit tensor. Additionally, the function $a_{\alpha\beta}$ interpolates between a purely isotropic mobility ($a_{\alpha\beta}=0$) and a pure tangential one ($a_{\alpha\beta}=1$). Thus, a flux across the interface is generated only if the interface is not in its equilibrium shape,

[TABLE]

. The operator $\operatorname{limit}_{1}$ restricts the interpolation function $a_{\alpha\beta}$ to the interval $\left[0,1\right]$ with

[TABLE]

If the sharp interface limit is not of major interest, it is tempting to use the simpler isotropic mobility tensor,

[TABLE]

. However, it is shown in Gugenberger et al. (2008) that an isotropic mobility tensor is less suitable to model the dynamics of surface diffusion. Therefore, with the exception of a single test shown in Fig. 3, all the simulations reported here are performed using the non-isotropic mobility tensor, Eq. (18).

II.3 Simulation Details

It is noteworthy that the derivatives of the obstacle potential, $\partial\left|\phi_{\alpha}\right|/\partial\phi_{\alpha}$ and $\partial\left|\phi_{\beta}\right|/\partial\phi_{\beta}$, are not continuous at $\phi=0$ and $\phi=1$. Therefore, a so called bent-cable model Chiu (2002) is used in order to piecewise interpolate between the linear derivative of the double obstacle and the non-linear interpolation at the minima.

The model is implemented inside the open source software project OpenPhase (www.OpenPhase.de). A finite difference scheme with a 27-point stencil Spotz and Carey (1995) is used for the discretization of the Laplacian in Eq. (16). This high number of stencil points is used in order to avoid numerical errors. We also tested the method with a 7-point stencil and the results showed no significant difference here. This observation is corroborated by similar studies with a focus on curvature evaluation vak . The divergence and the gradient in Eq. (13) are both discretized with a 3-point first order central finite difference scheme. A first order explicit Euler method is used for the discretization of the time derivative. Because phase-boundary diffusion becomes more significant on small length scales, a spacial discretization of $\Delta x=10^{-9}$m and temporal discretization of $\Delta t=10^{-13}$s is considered. The surface/interface energies used in this study are of the order of $\sigma\sim 1$ J/m2, which is typical for metallic systems. The mobility coefficient of Eq. (13) is set to $\tilde{M}_{\alpha\beta}=10^{-25}\text{m}^{2}/\text{Js}$. This value of the interface mobility does not correspond to a real physical system but is necessary to keep the algorithm stable. Results obtained in this work are thus of generic rather than material specific nature. To better reflect this feature, all lengths and times are reported below in units of $\Delta x$ and $\Delta t$. If not otherwise stated, the interface width is set to $\eta=10\Delta x$. In order to calculate the right hand side of Eq. (14), three additional boundary grid cells for $\phi$ are used around the computation domain. One way to determine these boundary cells is to use periodic boundary conditions. The volume of the phases is in this case conserved up to the machine’s roundoff error.

III Results and Discussion

The phase-boundary diffusion model is first tested with regard to the equilibrium interface profile for a spherically symmetric phase-field in Sec. III.1. A detailed study of thermal grooving is presented in Sec. III.2. Assuming that the late time profile obtained from these simulations is a good approximation for equilibrium shape of a three phase boundary, it is examined in Sec. III.3 whether these simulations are consistent with the von Neumann’s triangle relation. As an example for the influence of periodic boundaries, we simulate a 3D tessellation with octahedra in Sec. III.4. In Sec. III.5, the annealing of nano-clusters on a deformable, initially flat, surface is investigated. The results obtained from these simulations are discussed in the context of the available literature.

III.1 Interface Profile

As a very first test, we consider relaxation of the interface towards its equilibrium profile, similar to the test performed in Gugenberger et al. (2008). In the case of only two phases and a planar interface, the analytical solution of the interface profile is known at equilibrium (see App. A). We show numerically that the phase-field profile obtained from the planar interface is also a good approximation for the profile of a spherically symmetric phase-field. For two phase-fields with $\phi_{\alpha}=1-\phi_{\beta}=\phi$,

[TABLE]

in which, $r$ is the radial coordinate and $R$ the radius of the sphere.

Equation (22) is used as initialization for the phase-field, but with a wider interface width ($\eta_{\text{initial}}=2\eta$). Figure 3 shows the value of $\phi$ along the center line of the sphere for different time steps. In the current set-up, the phase-field at the initial time, $t=0$, is already in a spherical shape, but the interface is wider than the equilibrium solution. Thus, diffusion is expected to occur only in the direction normal to the interface. A survey of Fig. 3 reveals that isotropic and non-isotropic tensorial mobilities result in different relaxation rates towards the equilibrium profile. This is closely related to the fact that, in the case of a non-isotropic mobility tensor, the matrix-elements are chosen such that the transport is mainly allowed along the tangential direction, while it is rather restricted along the normal direction. In the present set-up, where the driving force acts only along the normal direction, this leads to a slower diffusion as compared to an isotropic mobility tensor, where the diffusion process is not restricted along any direction.

III.2 An Example of Thermal Grooving

Mullins has investigated the mechanisms of surface diffusion and evaporation-condensation which lead to the formation of surface grooves in a polycrystal that is heated up to elevated temperatures Mullins (1957). An initially flat surface with a grain boundary perpendicular to the surface is considered. For surface diffusion, a first order approximation of the surface profile evolution,

[TABLE]

has been obtained for small slopes $m$ and with the parameter $B=\tilde{\mu}\sigma cV_{m}^{2}$ (see Mullins (1957)). The values for $a_{n}$ are given in Tab. 1. The profile of Eq. (23) is shape invariant and its amplitude, $Bt$, can not be changed without stretching it along the $x$-axis. This solution suggests that the groove profile increases its amplitude over time without changing its shape.

In order to examine the analytically obtained results of Mullins, we consider a quasi-two dimensional set-up, where three phase-fields are initialized so that they form a flat surface with an interface perpendicular to the surface (Fig. 4a,b). The evolution of the surface profile is shown in Fig. 4c. Two of the three phase-fields represent the solid grains. The third phase-field stands for the surrounding gas phase, $\phi_{\text{gas}}$. The surface profile is taken as the contour of the gas phase with $\phi_{\text{gas}}=0.5$, which is reconstructed with a fourth-order polynomial and a bisection method. This is necessary to obtain the equilibrium angles with a suitable accuracy for the following analysis.

Between the initial and equilibrium configurations, the triple junction is moving downwards, causing the interfaces to bend upwards. Additionally, one can see that the maximum of this bended interface travels outwards. In a private communication to Mullins Mullins (1957), Hillert pointed out that since the flow of matter immediately beyond the point of inflection is toward the origin, and since the curve has a fixed shape with enlarging size, there must be oscillations of the surface profile about the $x$-axis. These oscillations are not predicted in the original paper by Mullins Mullins (1957) who did not exclude this possibility but suspected that it would be difficult to observe these oscillations, due to the decreasing amplitude of the surface profile Mullins (1957). The results obtained in the seminal work of Mullins, who assumes the shape invariance of the surface profile, has been confirmed by numerical solution of the underlying partial differential equations Spotz and Carey (1995); Zhang and Wong (2002). Recently, the above mentioned non-monotonic surface profile has been reported in experiments in the case of a tungsten polycrystal (see Sachenko et al. (2004) and references therein). As seen in Fig. 4, the results of the current multi-phase-field model for surface diffusion provide an independent numerical evidence for the oscillatory propagation of the surface profile during thermal grooving. These results are also in qualitative agreement with two dimensional calculations based on a variational approach Hackl et al. (2013).

A comparison of Eq. (23) with the obtained surface profile (see Fig. 4) shows that Eq. (23) is a good approximation for early times of the evolution and near the triple junction. However, our simulation reveals that, the more the wave front propagates, the flatter its profile becomes in comparison to Eq. (23). A profile similar to the one obtained in this study and the profile obtained by Mullins albeit without oscillations has been experimentally observed in Gladstone et al. (2001).

Since the oscillation pattern increases in size with time, it is more probable that they can be observed in the later stages of thermal grooving. One can however, also not fully exclude the possible role of thermal fluctuations in damping the oscillations. In order to investigate this aspect for the time evolution of the phase-field variable, Eq. (13) shall be extended to a Langevin-type equation including the effect of thermal noise. In addition, on the surface of a polycrystalline body which maintains multi grooving junctions, oscillations interfere with one another. Thus, spacing between the junctions may also play an important role. This feature is not considered in the present simulations as well. A detailed investigation of the effects arising from thermal noise and multiple grooving junctions is left for future work.

III.3 Von Neumann’s triangle relation

In this section, we benchmark our results versus Young’s law at triple junctions. For this purpose, we take the example of thermal grooving and investigate the dependency of the equilibrium groove angle on the grain boundary energy. Here, the equilibrium state is defined as the steady and quasi time independent profile, which establishes during late stages of thermal grooving. In a sharp interface picture, the groove angle is defined as the angle between the left and right tangent lines of the surface profile at the triple junction. In close similarity to this, we use the right and left tangents of the corresponding dual phase contour lines near the triple junction. In order to improve numerical accuracy, a high order polynomial interpolation is used to calculate the tangent lines.

Since Eq. (10) is derived with the assumption of straight/planar interfaces, the accuracy will depend on the boundary conditions in the simulation box. In order to account for this fact, the simulation of thermal grooving has been repeated for different lengths of the computation domain along the $x$-direction.

We find that it is possible to obtain a relation between the size of the computation domain and the equilibrium angle (see Fig. 5a). With this knowledge, one can minimize numerical errors and estimate the effect of the boundary. If one assumes that the surface energy, $\sigma_{\text{surface}}$, is the same for both grains, then the groove angle, $\theta$, can be calculated using Eq. (10),

[TABLE]

. The angles calculated by the simulations with the present method, and those predicted by Eq. (24), are listed in Tab. 2. The same data is plotted in Fig. 5b, highlighting the close agreement between numerically obtained results and Eq. (24).

III.4 Quadruple Junctions

Although the above simulations of thermal grooving are performed in three dimensions, they are essentially equivalent to a two dimensional situation due to the cylindrical symmetry of the set-up (no $z$-dependence). For multi-grain (multi-phase) systems in 3D, however, vertex points (quadruple and higher order junctions) may also exist in which more than three grains meet. In order to investigate this type of scenario, a multi-grain set-up is constructed by dividing the simulation domain into eight cubes such that only triple lines and quadruple points form between the grains (see Fig. 6a). By surface minimization, one can see that the grains arrange themselves into truncated octahedra (see Fig. 6c). The equilibrium angle between the triple lines at a quadruple junction is either $90^{\circ}$, when the lines are connected by a square plane, or $120^{\circ}$, when the lines are connected by a hexagonal plane. We find that the dihedral angle at a triple line between two hexagonal planes of truncated octahedron is roughly $109^{\circ}$ and the angle between a hexagonal and a square plane is $\approx 125^{\circ}$. These are well in line with the values of $109^{\circ}28^{^{\prime}}16^{"}$ and $125^{\circ}15^{^{\prime}}51^{"}$ expected from geometrical consideration. However, from a generalization of the von Neumann’s triangle relation to quadruple junctions, one would expect that the system does not form truncated octahedra and that all dihedral angles have the same value. The difference between the simulation results and the prediction of a generalized von Neumann’s triangle relation for quadruple junctions is probably due to the conservation constraint in the periodic set-up of simulations.

It is quite simple, however, to show that the sum of forces which act alongside the triple lines on a quadruple point is zero (see App. B).

III.5 Droplets on a Deformable Surface

An important application of the phase-boundary diffusion model proposed in this work is the simulation of nano-clusters on free surfaces. In particular, the model allows to study the entrenching of nano-clusters, reported in recent experiments Köhler et al. (2013).

In order to study the influence of the interface energy on the equilibrium contact angles between the nano-clusters and the free surface, we have placed a spherical droplet on a surface and observed its wetting dynamics. The equilibrium configurations for different interface energies are shown in Fig. 7. A complete wetting of the surface, beyond Eq. (10), is visible in Fig. 7c.

The equilibrium configurations shown in Fig. 7 do not resemble the geometry of the nano-clusters observed in Köhler et al. (2013) where elevated perimeters have been observed which surround the ‘crater-like’ entrenchments of the nano-clusters. This perimeter can be explained as intermediate state of the entrenching process, in a sense that the entire system is not in equilibrium. This interpretation is confirmed by our simulations, shown in Fig. 8, where four initially connected droplets on surface are simulated with identical interface energies. One can clearly see that the droplets separate from each other and entrench into the surface. The perimeter is visualized in Fig. 8b. Thus, the present model is not only able to reproduce experimentally observed results but also allows to uncover the transient character of the elevated perimeter in the entrenching process. A full investigation of the multi-nano-clusters on free surfaces is in progress.

IV Summary

A conserved multi-phase-field model is proposed to describe the physical effect of surface/phase-boundary diffusion. Starting from the non-conserved multi-phase-field model Steinbach (2009), a pairwise continuity equation is proposed for temporal evolution of the surface/phase-boundaries by curvature-driven diffusion.

The model is applied to investigate the dynamics of the surface profile during thermal grooving. Simulation results indicate that the shape invariant solution of Mullins Mullins (1957) is a good approximation to the early stages of the thermal grooving process. Moreover, new evidence is provided for the oscillations of the surface profile, first proposed by Hillert. Additional simulations of thermal grooving at variable interface energies are then performed in order to gain further confidence on the reliability of the present approach. This is achieved by a check of the late time profiles obtained from these simulations against the von Neumann’s triangle relation for equilibrium angles at a triple junction.

As a forecast on future applications of the proposed model, the annealing dynamics of nano-clusters on initially flat surfaces is investigated. The reported crater-like structure with an elevated perimeter Köhler et al. (2013) is found to be a transient non-equilibrium state during nano-cluster annealing, thus shedding light onto this complex process from a dynamic perspective. The proposed multi-phase-field method is generic and can be used to study any type of surface and phase-boundary phenomena which contains conserved fields.

Appendix A Static Equilibrium Solution

In the case of two phase-fields with $\phi_{\alpha}=1-\phi_{\beta}=\phi$, the generalized diffusion potential introduced in Eq. (14) simplifies to

[TABLE]

, where we assumed $\phi\left(1-\phi\right)\in\left[0,1\right]$ so that the partial derivate of the obstacle potential (2nd term in Eq. (16)) reduces to $\partial\left|\phi\left(1-\phi\right)\right|/\partial\phi=2\phi-1$.

Using this information, one finds that the equilibrium solution of Eq. (13) for a planar interface, situated at $x=0$ between two phases, is identical to Eq. (22) with the difference that the radial distance $r$ shall be replaced by the Cartesian coordinate along the direction normal to the interface, say $x$.

Appendix B Forces at a quadruple junction of a tessellation with truncated octahedra

If one considers a quadruple in a tessellation with truncated octahedra as in Fig. 6c, one sees that four triple lines are connected to that quadruple junction. The configuration of the phase-fields at each of these triple lines is the same. So it is reasonable to assume that the magnitude, $m$, of the forces alongside the triple line is the same and only their orientation differs. These forces can be written as follows: {dgroup}

[TABLE]

[TABLE]

[TABLE]

[TABLE]

.

One can easily check that the sum of $\boldsymbol{\mathbf{F}}_{1}$, $\boldsymbol{\mathbf{F}}_{2}$, $\boldsymbol{\mathbf{F}}_{3}$ and $\boldsymbol{\mathbf{F}}_{4}$ is zero, and the angles between the forces are either $90^{\circ}$ or $120^{\circ}$.

Acknowledgements.

Financial support by the German Research Foundation DFG within Priority Program SPP1713 under the grants DA 1655/1-1, STE 1116/20-1 and the DFG-project VA 205/17-1 is gratefully acknowledged.