1 Introduction

The numerical simulation of phenomena like phase growth has gained popularity over the last decades, and even more so in the last years when it became feasible to make predictions about the failure behavior of various components, e.g., in microelectronics. There are essentially two very different approaches to modeling the domain of reactive diffusion. First, the phase field method [8] has to be mentioned as the most elegant approach. Other than to directly solve the usual conservation/diffusion formalism the evolution equation is obtained from solving a variational minimization problem containing also non-conserved phase describing variables. The functional to be minimized is usually constructed using empirical functions that are at least three times continuously differentiable such as a double-well potential [8]. In the simplest cases, the variable of interest (for example the concentration) is just penalized for obtaining values beyond the usual equilibrium concentrations. Additionally, the interface curvature could be penalized (see for example “Cahn–Hilliard Equation” [4]). This concept can even be expanded to model microstructural grain growth (see “Allen–Cahn Equation” [4]). On the other hand, we have, as indicated before, methods that aim at directly solving conservation/diffusion evolution equations, taking the phase-varying properties into account by introducing them as concentration-dependent functions. Therefore, no empirical functions like the double-well potential have to be introduced and the equations solely include formalisms obtained directly from thermodynamics. These methods shall loosely be called sharp interface models, which is of course due to their tendency to create sharp transitions between phases. This is in strong contrast to the PFM that is characterized by smooth concentration transitions in interface domains, where the width of the interface can be controlled by acting on empirical model parameters. Although sharp interface methods promise higher efficiency in the treatment of complex microstructures, special care must be taken to choose the right numerical method in the implementation. For the Galerkin method, there occur numerical stability problems when trying to solve reactive diffusion problems with the methodology described in this paper. To get around these stability problems, the Finite Difference Method, the Finite Volume Method (FVM) and the Discontinuous Galerkin Method can be employed. The DGM is of special interest for scientists that are familiar with the FEM or want to solve coupled multi-physics problems that require maximum flexibility regarding solution methods. This is because the DGM was designed to fit the FEM-framework.

1.1 Properties of the fluxes in reactive diffusion

Consider a system containing n substitutional components with their respective site fractionsFootnote 1 given by \(y_i\) for \(i = 0,1,..,n\), where \(y_0\) corresponds to vacant (empty) sites in the lattice (vacancies). The flux \(\underline{j}_k\) of a species k is, according to Eq. (1), given by the sum of the gradients of chemical potentials \(\mu _i\) scaled with so-called Onsager coefficients \(L_{ik}\) [8, 10, 11, 16, 36].

$$\begin{aligned} \underline{j}_k&:= \sum _{i=1}^{n} L_{ik} \nabla \mu _i \;\;\;\; k = 1, .. , n \end{aligned}$$
(1)

In general, the chemical potentials and Onsager coefficients can be dependent on the site fractions of all components as well as the pressure and temperature. For the derivations in this paper, it is only necessary to consider their dependency on the composition:

$$\begin{aligned} L_{ik}&:= L_{ik}(y_0, y_1, ... , y_n) , \end{aligned}$$
(2)
$$\begin{aligned} \mu _{i}&:= \mu _{i}(y_0, y_1, ... , y_n) , \end{aligned}$$
(3)

Additionally, the restrictions (4) and (5) must hold at all times.

$$\begin{aligned} \sum _{i=0}^{n} y_i = 1, \end{aligned}$$
(4)
$$\begin{aligned} \sum _{i=0}^{n} \underline{j}_i = \underline{0}, \end{aligned}$$
(5)

The chemical potential of component i can, according to a derivation in [32], be computed from the molar Gibbs Free Energy \(G_m\) as

$$\begin{aligned} \mu _{i} := G_m + \sum _{k=1}^{n-1} (\delta _{ik} - y_k) \frac{\partial G_m}{\partial y_k} \;\;\;\; i = 1, .., n , \end{aligned}$$
(6)

with the assumption that \(y_n\) is, according to Eq. (4), expressed in terms of the other site fractions and does not appear as independent variable:

$$\begin{aligned} y_{n} := 1 - \sum _{i=0}^{n-1} y_{i} \end{aligned}$$
(7)

In (7), \(\delta _{ik}\) denotes the Kronecker Delta. As multiple phases are considered, the molar Gibbs Free Energy has to be phase dependent as well. In this article, the special case of diffusional phase transformations is considered, following the works of Svoboda [34, 35, 37]. This implies that the lattice conversion is instantaneous and that the phase is completely determined by the composition of the lattice. This assumption is reasonable for slow diffusion as observed in solid metals, where the transformations are purely diffusion driven. For a general approach, it has to be acknowledged that the lattice conversion also takes up time which is considered by the PFM [8]. The derivation is continued assuming a system consisting of p phases named \(\phi _{l}\) for \(l = 1, 2,.., p\). It is furthermore assumed that the molar Gibbs Free Energy of every phase is known and given in terms of the composition:

$$\begin{aligned} G_{m}^{\phi _\ell } := G_{m}^{\phi _\ell }(y_{1}^{\phi _\ell },y_{2}^{\phi _\ell }, .., y_{n}^{\phi _\ell }) \end{aligned}$$
(8)

The newly introduced variables \(y_{1}^{\phi _\ell }, y_{2}^{\phi _\ell },.., y_{n}^{\phi _\ell }\) refer to the site fractions of the components explicitly associated with phase \(\phi _\ell \). These variables have to satisfy the following equations:

$$\begin{aligned} y_i = \sum _{\ell =1}^{p} y_{i}^{\phi _\ell } \phi _\ell \;\;\;\;\;\; i = 1, .., n \end{aligned}$$
(9)

Where \(\phi _\ell \) are phase describing functions (as known from the PFM) that obey the constrains:

$$\begin{aligned}&0 \le \phi _\ell \le 1 \;\;\;\;\;\; \ell = 1, .., p \end{aligned}$$
(10)
$$\begin{aligned}&\sum _{\ell =1}^{p} \phi _\ell = 1 \end{aligned}$$
(11)

Wherever \(\phi _{\ell } = 1\), the lattice consists entirely of phase \(\phi _{\ell }\) and so on. Wherever multiple phase describing functions are nonzero, an interface is assumed. Using the introduced functions, the total Gibbs Free Energy of the multi-component, multi-phase system can be formulated:

$$\begin{aligned} \mathcal {G} := \int _{\Omega } \; \frac{1}{\Omega _{m}} \sum _{\ell =1}^{p} \phi _\ell G_{m}^{\phi _\ell } \; \textrm{d}\mathbf {\Omega } \end{aligned}$$
(12)

Note that \(\Omega \) denotes the mathematical domain of the whole system while \(\Omega _{m}\) refers to the molar volume of components. In this article, the molar volume is assumed to be constant for all components and phases, respectively. There are, however, general formulations that review reactive diffusion with phase-dependent molar volumes [36]. The Gibbs Free Energy functional is the starting point for derivations for both the PFM as well as the reactive diffusion methodology developed by Svoboda [35,36,37]. While the flux defined in Eq. (1) follows from maximization of the entropy production rateFootnote 2 the evolution of phase describing functions \(\phi _{\ell }\) strives to minimize the Gibbs Free Energy functional. The main idea of Svoboda’s methodology is to solve for the phase describing functions beforehand, therefore, assuming instantaneous phase transitions and viewing phases entirely determined by composition. Therefore, the phase describing functions are not needed in the presented method. This is achieved by solving the following optimization problem:

$$\begin{aligned} G_{m}(y_1, y_2, .., y_n) := \min _{ y_{i}^{\phi _{\ell }}, \phi _\ell } \mathcal {G} \end{aligned}$$
(13)
Fig. 1
figure 1

Molar Gibbs free energy curves of three phases in a binary system

Full size image
Fig. 2
figure 2

Chemical potentials of component 1 and 2 for the system depicted in Fig. 1. The domains where a certain phase is stable are colored according to the color scheme in Fig. 1

Full size image

Subject to the constraints (9), (10) and (11), where (9) must be incorporated to the functional by using Lagrange multipliers. For the inequality constraints (10), special methods for highly nonlinear systems have to be applied such as homotopy methods. This is the usual approach to compute phase diagrams also referred to as CALPHADFootnote 3 methodology [33]. In Figs. 1 and 2, an example for the computation of chemical potentials from Gibbs Free Energy curves using Eqs. (7) and (13) is depicted. Solving optimization problem (13) can be interpreted geometrically as constructing a convex hull curve around the molar Gibbs Free Energy curves seen in Fig. 1. Therefore, it is also possible to calculate the chemical potentials by the so-called common-tangent construction [16]. It can be thermodynamically favorable to have multiple phases at different fractions at a given chemical composition. Regions where this is the case are the interfaces between two phases or junctions among multiple phases. These regions are characterized by the fact that the overall chemical potential of a component is constant, as can be seen in Fig. 2. Everywhere else, that is where a certain phase is clearly favorable, the chemical potential of the preferable phase is active [25]. For the overall system, we can at least conclude that the chemical potential of a component is a \(C^0\) continuous function of site fractions. Therefore, its derivative is a discontinuous function. Due to Eq. (1), it can be concluded that fluxes are discontinuous for the reviewed model. This is the main reason why numerical treatment with the Galerkin method will fail, and its discontinuous counterpart has to be employed.

2 Methodology

2.1 The discontinuous Galerkin method for reactive diffusion

The most fundamental fact about diffusion of conserved quantities can be stated as follows: The total change in the amount of a conserved species \(y\) in the domain \(\Omega \) is equal to the amount that enters or leaves the domain across its boundary \(\partial \Omega \). In mathematical terms, this can be expressed by the following integral equation:

$$\begin{aligned} \int _\Omega \dot{y} \; \textrm{d}\pmb {x} + \int _{\partial \Omega } \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} = 0, \end{aligned}$$
(14)

To derive the diffusion equation from statement (14), the divergence theorem and the fundamental lemma of the calculus of variations [13] have to be applied in order to achieve a localizationFootnote 4 of the conservation law. In the discontinuous Galerkin method, the same fundamental mathematical considerations as in equation (14) apply. Other than in the derivation of the variational form of the diffusion equation for the GM, however, the divergence theorem must not be applied as it requires the flux to be a sufficiently smooth (differentiable) function [13]. As was shown in the previous chapter, this is a property the flux cannot satisfy in reactive diffusion (as it is discontinuous across phase boundaries). This is the reason the GM does not suffice to solve reactive diffusion problems as discussed in this paper. The DGM uses Eq. (14) directly, but it is not only applied to the whole domain \(\Omega \) but to every finite element \(\Omega \) consists of. In the following, the terms “finite element” will be replaced by the more modern designation “cell”. The domain of an individual cell e is denoted \(\Omega _e\) with its boundary \(\partial \Omega _e\).

$$\begin{aligned} \int _{\Omega _e} \dot{y} \; \textrm{d}\pmb {x} + \int _{\partial \Omega _e} \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {S} = 0, \end{aligned}$$
(15)

The total amount that flows into or out of the cell through its boundary causes a change in the amount \(y\) that remains in the cell. The localization in the DGM is therefore achieved by requiring the species to be conserved with respect to each cell. By assuming the site fraction to be constant over this cell, we also imply that a specific function space for the trial- and test functions is to be used. Let \(D^0\) denote the function space of DG functions of degree zero. \(D^0\) contains all piecewise constant and globally discontinuous functions. That means constant over cells and discontinuous at the cell boundaries, also called “facets”. Every cell is assigned just one trial and one test function which obtains the value one in its assigned cell and the value zero everywhere else. Special attention is needed to consider the boundary terms. It is clear that every cell has its boundary \(\partial \Omega _e\), but the mathematical boundary of the domain \(\partial \Omega \) consists in general of “facets” [20] of multiple cells. Therefore, certain cell boundaries will be part of the mathematical boundary of the problem as well. Facets of cells are therefore always shared between two cells or are part of the mathematical boundary. To increase readability, the boundary of individual cells with shared facets shall be named \(\Gamma \) as depicted in Fig. 3. Requiring Eq. (15) to hold for every cell can now be expressed by the following variational formulation:

$$\begin{aligned} \int _{\Omega } \dot{y} v \; \textrm{d}\pmb {x} + \int _{\partial \Omega } v \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} + \int _{\Gamma } [v \underline{j}] \; \textrm{d}\pmb {S} = 0 \;\;\;\; \forall v \in D^0, \end{aligned}$$
(16)
Fig. 3
figure 3

Schematic domain displaying the different parameterizations of the mathematical boundary and the internal cell boundaries

Full size image

To account for the discontinuities at the internal boundaries \(\Gamma \), the jump operator “\([\cdot ]\)” is introduced [20, 38]:

$$\begin{aligned}{}[v]:= (v^+ - v^-) \underline{n}, \end{aligned}$$
(17)

The jump operator can be interpreted as a numerical tool that allows to account for the discontinuous “jump” when integrating over the internal cell facets \(\Gamma \) as depicted in Fig. 4. Imagine the individual facet \(\Gamma _{-}^{+}\) to be the shared facet of two cells \(e^+\) and \(e^-\). When integrating over \(\Gamma _{-}^{+}\), the function values on both sides of the facet have to be considered:

$$\begin{aligned} \int _{\Gamma _{-}^{+}} ( f^+ v^+ - f^- v^- )\;\underline{n} \cdot \textrm{d}{} {\textbf {S}} = \int _{\Gamma _{-}^{+}} [f v] \; \textrm{d}S, \end{aligned}$$
(18)
Fig. 4
figure 4

Schematic of two neighboring cells \(e^{-}\) and \(e^{+}\) that share the Facet \(\Gamma _{-}^{+}\)

Full size image

For the flux between two cells, we choose, according to Eq. (1), \(\underline{j} = L(y) \nabla \mu (y) \), where \(L(\cdot )\) and \(\mu (\cdot )\) are arbitrary nonlinear functions. On the mathematical boundary, we keep \(\underline{j}\) to represent the inflow or outflow to be applied by the user as von Neumann boundary condition.

$$\begin{aligned} \int _{\Omega } \dot{y} v \; \textrm{d}\pmb {x} + \int _{\partial \Omega } v \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} + \int _{\Gamma } [v L \nabla \mu ] \; \textrm{d}\pmb {S} = 0 \;\;\;\; \forall v \in D^0, \end{aligned}$$
(19)

To complete the derivation, another special operator has to be introduced, namely the average operator “\(<\cdot>\)”, which is defined as follows [38]:

$$\begin{aligned} <v>:= \frac{1}{2} (v^+ + v^-), \end{aligned}$$
(20)

It computes, as the name suggests, the average value of two cells in their shared facet. This operator comes into play when an expression containing the jump operator shall be simplified. For this, the jump identity is introduced [38]:

$$\begin{aligned}{}[uv]:= [u]<v> + <u>[v], \end{aligned}$$
(21)

Application to variational form (19) yields:

$$\begin{aligned}{} & {} \int _{\Omega } \dot{y} v \; \textrm{d}\pmb {x} + \int _{\partial \Omega } v \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} + \int _{\Gamma } [v]<L> \;<\nabla \mu> \cdot \; \underline{n} \; +<v> [L]<\nabla \mu> \cdot \; \underline{n} \; \nonumber \\{} & {} \quad +<v> \; <L> [\nabla \mu ] \; \textrm{d}\pmb {S} = 0 \;\;\;\; \forall v \in D^0, \end{aligned}$$
(22)

Since the Onsager coefficient L is in general a nonlinear but continuous function of the site fraction y, it must also behave continuous in a weak sense. This is achieved by setting \([L]=0\). As was pointed out in the introduction, the equilibrium conditions require the chemical potential to be constant over interfaces. Therefore, we additionally require \([\nabla \mu ]=0\) to hold over all facets and drop the corresponding terms.

$$\begin{aligned} \int _{\Omega } \dot{y} v \; \textrm{d}\pmb {x} + \int _{\partial \Omega } v \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} + \int _{\Gamma } [v]<L> \; <\nabla \mu > \cdot \; \underline{n} \; \textrm{d}\pmb {S} = 0 \;\;\;\; \forall v \in D^0, \end{aligned}$$
(23)

While it is clear that \(\nabla \mu \) must vanish in the cells due to the type of trial function used, it is, however, up to the user to define the behavior of \(\nabla \mu \) at the facets. As there occur jumps in the site fraction between cells, it would be unreasonable to neglect this term. In order to use this method, however, the behavior of this function has to be further investigated. Therefore, we concentrate again on an individual facet between two cells \(e^+\) and \(e^-\) and replace the gradient with a difference quotient:

$$\begin{aligned} \int _{\Gamma _{-}^{+}}< \nabla \mu> \cdot \; \underline{n} \; \textrm{d}\pmb {S} \approx \int _{\Gamma _{-}^{+}} \; < \frac{\mu (y^+) - \mu (y^-)}{l_e} > \; \textrm{d}\pmb {S}, \end{aligned}$$
(24)

Here, \(l_e\) denotes the cell length. As we already required \([\nabla \mu ]=0\) to hold, it is clear, that \(<\nabla \mu > \; =\nabla \mu \) must be true:

$$\begin{aligned} \int _{\Gamma _{-}^{+}} \nabla \mu \cdot \underline{n} \; \textrm{d}\pmb {S} \approx \int _{\Gamma _{-}^{+}} \; \frac{\mu (y^+) - \mu (y^-)}{l_e} \; \textrm{d}\pmb {S}, \end{aligned}$$
(25)

Definition (25) can be simplified to yield the jump operation:

$$\begin{aligned} \int _{\Gamma _{-}^{+}} \nabla \mu \cdot \underline{n} \; \textrm{d}\pmb {S} \approx \int _{\Gamma _{-}^{+}} \; \frac{1}{l_e} [\mu ] \; \textrm{d}\pmb {S}, \end{aligned}$$
(26)

While fluxes are not allowed to have jumps across facets, the mole fraction \(y\) is, at least in this method allowed to have discontinuities. This is a widely used simplification for diffusion problems (see symmetric interior penalty method (SIP/DGM) [6, 7, 30]). The fundamental cause for the flux between two cells is therefore the difference in the site fraction of the cells. Assuming the relationship between the gradient and the jump also explains why this method gives almost identical results as the finite difference method. Substituting this relationship gives the final variational formulation.

$$\begin{aligned} \int _{\Omega } \dot{y} v \; \textrm{d}\pmb {x} + \int _{\partial \Omega } v \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} + \int _{\Gamma } \frac{<L>}{l_e} \; [v][\mu ] \; \textrm{d}\pmb {S} = 0 \;\;\;\; \forall v \in D^0, \end{aligned}$$
(27)

Interestingly, the strong form of the diffusion equation is actually never used in this derivation as this requires the application of the divergence theorem that can necessarily only be applied to smooth functions. Note that regular diffusion problems (simple diffusion problems where the fluxes are evoked by Fick’s law \(\underline{j} = D \nabla y\), where \(D\) is the constant diffusion coefficient) are very well investigated by the DGM community. There even exists a variational formulation for these problems that is independent of the degree of the used trial and test functions [38]. This formulation yields the following variational formulation if functions of degree zero are utilized, such as in the derivation above:

$$\begin{aligned} \int _{\Omega } \dot{y} v \; \textrm{d}\pmb {x} + \int _{\partial \Omega } v \underline{j} \cdot \underline{n} \; \textrm{d}\pmb {s} + \int _{\Gamma } \frac{D}{l_e} \; [v][y] \; \textrm{d}\pmb {S} = 0 \;\;\;\; \forall v \in D^0, \end{aligned}$$
(28)

The similarities to the variational form for reactive diffusion with the fluxes in canonical form come as no surprise. By choosing \(L=D=const.\) and \(\mu (y)=y\), the variational forms become identical, as expected. Although this is a new variational formulation, it can therefore be assumed that the assumptions made are correct and lead to valid results. For completeness, the reader is also referred to [14] where a similar case was studied. Another notable mention is [9] where the so-called harmonic average instead of the average operator of this work was employed to solve the advection–diffusion equation with locally vanishing diffusivity. Note that the discontinuous Galerkin Method for diffusion (28) is, if forward Euler time integration is used, in fact so closely related to the “forward time, central space” diffusion scheme of Finite Difference Methods that the stability criterion of the latter can be applied to find a stable time increment [38]:

$$\begin{aligned} 0 \le \frac{D \Delta t}{(l_e)^2} \le \frac{1}{2} , \end{aligned}$$
(29)

3 Results

3.1 Vacancy-controlled diffusion of a binary, three-phase system in a two-dimensional domain

To generate the following solutions to the discussed variational problems, the legacy version of FEniCSFootnote 5 has been used. FEniCS [2, 20] is developed by many scientists and consists of the modules dolfin [22, 23], ffc [17, 21, 24] and fiat [18, 19]. This finite element software can be used with a python 3 [39] interface. All graphics have been generated with Paraview [1, 3]. To demonstrate the capabilities of the presented method, a microstructural evolution problem is considered. For a detailed derivation of the model, the reader is referred to [36]. A first implementation of the model can be found in [34] where the equations were solved using a finite difference scheme on a one-dimensional domain. The chosen system is inspired by the binary Cu-Sn system, although the simulation will not use real thermodynamical data and material parameters. This system is of technical importance since tin-based alloys are to replace lead-based solders in various microelectronic devices [15, 31]. These solders are used to attach and connect all kinds of microelectronic components to circuit boards and are often exposed to multiple stresses such as temperature fluctuations and mechanical loads. Due to the small scale of solders and the constant attempts of further downsizing, effects like phase growth cannot be neglected to fully understand the damage mechanisms. The two components of the model are given in terms of their site fractions \(y_1\) and \(y_2\), respectively. Additionally, the site fraction of vacancies \(y_{0}\) is considered. Using (4) as constraint, the system is fully defined by using \(y_{0}\) and \(y_{1}\) as variables since \(y_{2} = 1 - y_{1} - y_{0}\) (or \(y_{2} \approx 1 - y_{1}\), since \(y_{0} \ll y_{1}\)). Due to the incorporation of the vacancy mechanism, the model has several special features. Since vacant sites in the lattice are just empty space, they are not conserved and can emerge and vanish at sources and sinks for vacancies. Therefore, the evolution equations of the model are not strictly diffusion equations but a more general form of transport equations:

$$\begin{aligned} \dot{y}_0&= -\Omega _{m} \nabla \cdot \underline{j}_0 + (1-y_0)\alpha \end{aligned}$$
(30)
$$\begin{aligned} \dot{y}_{1}&= -\Omega _{m} \nabla \cdot \underline{j}_1 - \alpha y_1 \end{aligned}$$
(31)

Here, \(\alpha \) denotes the rate of generation and annihilation of vacancies. For a stress free system, it is defined as

$$\begin{aligned} \alpha := -\frac{\mu _{0}}{K \Omega _{m}}, \end{aligned}$$
(32)

where \(\mu _{0}\) is the chemical potential of vacancies and K is the bulk viscosity as introduced in [12]. In general, K is variable over the domain, obtaining different values at grain- and phase-boundaries, but for the sake of simplicity it is assumed to be constant in this example. The chemical potential of vacancies is defined as

$$\begin{aligned} \mu _{0} := RT \ln \left( \frac{y_{0}}{y_{0}^{\text {eq}}}\right) , \end{aligned}$$
(33)

where R and T are gas constant and temperature, respectively, and \(y_{0}^{\text {eq}}\) is the equilibrium site fraction of vacancies. The Onsager coefficients are computed according to

$$\begin{aligned} L_{ik} := A_k \delta _{ik} + \frac{(1-f)}{f} \frac{A_i A_k}{\sum _j A_j}, \end{aligned}$$
(34)

with f being the so-called vacancy correlation factor (\(f=0.7815\) for fcc and \(f=0.7272\) for bcc alloys) [34]. The coefficients \(A_k\) are defined as:

$$\begin{aligned} A_k := \frac{y_{0}}{y_{0}^{\text {eq}}} y_k \frac{D_k}{\Omega _{m} RT} \end{aligned}$$
(35)

\(D_k\) is the Diffusion coefficient of component k in the bulk. In general, the diffusion coefficients can be phase dependent, but again we assume them constant for the sake of demonstration. Having the equations of the model defined, we continue to solve them by the derived DG approach. The variational formulations of Eqs. (30) and (31) read:

$$\begin{aligned}&\int _{\Omega } (y_{0}^{t} - y_0) \; v \; \textrm{d} \pmb {x} + \sum _{i=1}^{2} \sum _{j=1}^{2} \int _{\Gamma } \Delta t \; \Omega _{m} \frac{<L_{ij}>}{l_e} \; [\mu _{j}][v] \; \textrm{d} \pmb {S} + \int _{\Omega } \Delta t \; (1-y_0) \; \alpha \; v \; \textrm{d}\pmb {x} = 0 \;\;\;\; \forall v \in D^0 \end{aligned}$$
(36)
$$\begin{aligned}&\int _{\Omega } (y_{1}^{t} - y_1) \; w \; \textrm{d} \pmb {x} + \sum _{j=1}^{2} \int _{\Gamma } \Delta t \; \Omega _{m} \frac{<L_{1j}>}{l_e} \; [\mu _{j}][w] \; \textrm{d} \pmb {S} - \int _{\Omega } \Delta t \; y_1 \; \alpha \; w \; \textrm{d} \pmb {x} = 0 \;\;\;\; \forall w \in D^0 \end{aligned}$$
(37)
Fig. 5
figure 5

Initial condition of \(y_1\). The mesh consists of 73441 quadrilateral finite elements

Full size image
Fig. 6
figure 6

Initial condition of the chemical potential of vacancies. \(y_0(t=0)=y_{0}^{\text {eq}}\)

Full size image
Fig. 7
figure 7

Site fraction of component one (\(y_1\)) after 0.3 [ms]. Clearly, two new phases \(\phi _2\) and \(\phi _3\) have formed

Full size image
Fig. 8
figure 8

Site fraction of component one (\(y_1\)) after 0.6 [ms]

Full size image
Fig. 9
figure 9

Chemical potential of vacancies after 0.3 [ms]. There is an excess of vacancies in the inclusions due to the Kirkendall effect

Full size image
Fig. 10
figure 10

Chemical potential of vacancies after 0.6 [ms]

Full size image
Fig. 11
figure 11

Site fraction of component one (\(y_1\)) after 0.6 [ms]. The mesh consists of 11881 quadrilateral finite elements

Full size image
Fig. 12
figure 12

Site fraction of component one (\(y_1\)) after 0.6 [ms]. The mesh consists of 4624 quadrilateral finite elements

Full size image

These two nonlinear variational problems were solved using a fully coupled backward Euler scheme on a two-dimensional unit square domain using a total of 73441 quadrilateral finite elements. The calculation was performed using the following parameters: , , , , , [–] and [–]. The Diffusion coefficient of component two, , was chosen twice as large as the diffusion coefficient of component one in order to trigger the Kirkendall effect [34]. The initial condition of \(y_1\) can be viewed in Fig. 5. The initial distribution was created using a NeperFootnote 6 generated Microstructure with 200 grains. Hundred of these grains were randomly chosen and assigned with an initial site fraction of 0.03 [–], while the other grains were set to a value of 0.97 [–]. As driving force of diffusion, the chemical potentials of Fig. 2 are used. The initial distribution of vacancies (\(y_0\)) was chosen to be constant over the domain with its equilibrium site fraction as can be seen in Fig. 6. Note that the chemical potential of vacancies was plotted instead of just the site fraction to clearly point out where shortage and exess of vacancies are located. In Figs. 7 and 8, the evolution of the microstructure is depicted at different timesteps using the site fraction \(y_1\). As can be seen, the phases \(\phi _2\) and \(\phi _3\) from Fig. 1 form at the interface between the initially unmixed components. The Cu–Sn system is also known to form two intermediate phases, often accompanied by so-called Kirkendall voids [34]. These voids or pores form due to vacancy accumulation triggered by the Kirkendall effect. Since Sn-atoms are known to have a higher diffusion coefficient than Cu-atoms, the flux of Sn-atoms out of Sn-rich regions can not be fully balanced by the Cu-atoms leading to an accumulation of vacancies in Sn-crystals (this is due to the flux balance equation) (5)). This causes the Kirkendall effect as can be observed in Figs. 9 and 10, where red and blue regions indicate excess and shortage of vacancies, respectively. Another observation that can be made in Fig. 9 is that regions with high interface curvature are at higher risk of Kirkendall voiding. As time progresses, it can be seen in Fig. 10 that the magnitude of the chemical potential of vacancies is generally sinking, which indicates that the system is approaching the final equilibrium state. To address the theme of mesh dependence, the calculation was also carried out with two coarser meshes. The resulting distribution of \(y_1\) at time \(t=0.0006\) [s] is given in Fig. 11 for a mesh with a total of 11881 quadrilateral finite elements and in Fig. 12 for a mesh consisting of just 4624 finite elements. While the method is stable and convergent for all mesh sizes it is obvious that the \(DG^0\) formulation benefits from fine meshing since the phase boundaries will be resolved more accurately.

4 Conclusion

In the present paper, the DGM was used to derive a new variational formulation (27) for reactive diffusion problems. The main advantage of the presented method over the PFM is the fact that thermodynamical principles can be used accurately and no empirical functions have to be introduced to smooth the phase transitions. The method also shows excellent properties in terms of stability and flexibility. The ability to model effects like the vacancy mechanism of solid body diffusion helps to obtain a better understanding of phenomena like Kirkendall voids and other kinds of vacancy related damage mechanisms. With its simplicity and robustness, the method can handle the added complexity of the vacancy mechanism well, promising to get hold of other ambitious phenomena in the future, such as segregation and trapping models. On the downside, it must be mentioned that for the presented method it is not yet possible to include interface energies to the Gibbs Free Energy, which is possible with the PFM. Furthermore, fine meshing is necessary to resolve detailed microstructures which is of course due to the properties of piecewise constant interpolation functions. Other than for the PFM where the focus lies on resolving the interface area in all its detail the presented method aims at simulating more macroscopic effects of diffusion on the overall bulk, thus idealizing the interface and representing it sharply. The DGMs discussed in this paper are restricted to the special case of trial and test functions of degree zero, e.g., piecewise constant functions. For some applications, it might be worth considering higher degree approximations. This case, however, would require deriving new variational formulations.