Strong-form meshfree collocation method for non-equilibrium solidification of multi-component alloy

Abstract

This work presents a strong form meshfree collocation method for a multi-phase field model with finite dissipation effects due to rapid solidification. We use the collocation method to simulate and study solidification of a low concentration (0.2 at% Sn) Al–Sn binary alloy system under periodic boundary conditions to address non-equilibrium solidification. Numerical implementation takes place through spatial discretization of the governing equations with the collocation method followed by application of the Crank–Nicolson method to integrate through time. Analysis begins with a benchmark, a simple two-grain case with symmetry in domain size, grain positioning, and boundary conditions to study the behavior of the field equations and key terms embedded within. This occurs by studying field and embedded term values along the axis of symmetry. Solidification analysis is then extended for 10 and 20 grains where upon full solidification, the regions with the highest overall concentrations exist within grain boundary region consisting for four or more adjacent grains. An analysis of alloy solidification over a substrate demonstrates epitaxial nucleation and growth.

Appendices

Appendix 1: Implementation of the Gibbs free energy

The bulk (Gibbs) free energy $f_{\gamma }$ for phase ${\gamma }$ evaluated in this work are based on the experimental thermo-chemical data and methods provided within the COST-507 [34]. Data consists of a collection of coefficients and interaction parameters used to calculate free energies of pure element and multi-component systems. The coefficients are temperature T dependent and are used with the power series to evaluate the molar free energy g for a pure species i in phase ${\gamma }$ with respect to its stable state enthalpy; i.e.

$$\begin{aligned}&g_i^\gamma (T) -h_i^\gamma (298.1K) = a+bT+cT\ln {T}\nonumber \\&\quad +dT^2+eT^3+fT^{-1}+gT^7+hT^{-9}. \end{aligned}$$

(33)

The free energy of multi-component system consists of contributions from each constituent species within the mixture and of mixing between species. Contributions from mixing are further classified as either ideal or non-ideal where the latter is modeled with a Redlich-Kister polynomial that contains temperature dependent interaction parameters $L^\nu$. For a substitutional three-component alloy system, the free energy is given by

$$\begin{aligned} \begin{aligned}G^\gamma (x_i,T) & = \sum _{i=1}^{3}x_ig_i^\gamma +RT\sum _{i=3}^{3}\sum _{j>i}^{3}x_ix_j\sum _{\nu =0}^{3}L_{ij}^\nu (x_i-x_j)^\nu \\&\quad +\sum _{i=1}^{3}\sum _{j>i}^{3}\sum _{k>j}^{3}x_ix_jx_k(x_iL^0_{ijk}+x_jL^1_{ijk}+x_kL^2_{ijk}) \end{aligned} \end{aligned}$$

(34)

where the right hand side terms respectively correspond to the unmixed, ideal mixing, binary non-ideal, and ternary non-ideal mixing terms. Invoking the chain rule to differentiate Eq. (34) with respect to composition yields

$$\begin{aligned} \frac{\partial G^\gamma (x_i,T)}{\partial x_m} = \frac{\partial G^\gamma (x_i,T)}{\partial x_i}\frac{\partial x_i}{\partial x_m} = G^\gamma _{,i} x_{i,m} \end{aligned}$$

(35)

with

$$\begin{aligned} \begin{aligned}\frac{\partial G^\gamma (x_i,T)}{\partial x_m} & = \sum _{i=1}^{3}x_{i,m}g^\gamma _i +RT\sum _{i=1}^{3}(x_{i,m}(\ln x_i+1))\\&\quad + \sum _{i=1}^{3}\sum _{j>i}^{3}(x_{i,m}x_j+x_{j,m}x_i)+\sum _{\nu =0}^{n}L_{ij}^\nu (x_i-x_j)^\nu \\&\quad + \sum _{i=1}^{3}\sum _{j>i}^{3}x_ix_j\sum _{\nu =0}^{n}\nu L_{ij}^\nu (x_i-x_j)^{\nu -1}(x_{i,m}-x_{j,m}) \\&\quad + \sum _{i=1}^{3}\sum _{j>i}^{3} \sum _{k>j}^{3}(x_{i,m}x_jx_k+x_ix_{j,m}x_k+x_ix_jx_{k,m})(x_iL^0_{ijk}\\&\quad +x_jL^1_{ijk}+x_kL^2_{ijk})\\&\quad + \sum _{i=1}^{3}\sum _{j>i}^{3}\sum _{k>j}^{3}x_ix_jx_k(x_{i.m}L^0_{ijk}+x_{j,m}L^1_{ijk}+x_{k,m}L^2_{ijk}). \end{aligned} \end{aligned}$$

(36)

Note that for the constraint, i.e. $\sum x_i =1$ we use the conditional relation to implement Eq. (36) in our computational analysis program:

$$\begin{aligned} x_{i,m} = \left\{ \begin{array}{cc} 1 &{} i=m \\ -1 &{} i \ne m \\ \end{array} \right. \end{aligned}$$

(37)

Figure 10 represents an output of the Al–Sn–Zn system along the Al–Zn, Zn–Sn, and Al–Sn surfaces for the liquid, FCC, BCT and HCP phases at a temperature of 600 K. Thermo-Calc plots are also presented for comparison. The derivatives used to construct the tangent lines for the Al–Zn (lower left plot) at $x_{\text {Al}} = 0.2$ are based on Eq. (36).

Fig. 10

Comparisons of the computed Gibbs energy (a–c) against Thermo-Calc software (d–f) at the same temperature, i.e. $T = 600$ K; the calculated Gibbs energy at along each edge of the Al–Sn–Zn ternary system, i.e. (Al, Sn, Zn = 0), (Al, Sn = 0, Zn) and (Al = 0, Sn, Zn) are shown

Appendix 2: Verification of the discrete differential operators

One way to measure the validity of computed discretized differential operators is through its fundamental properties, i.e. consistency [37]. Let $x_i$ where $x_1 = x$ and $x_2 = y$, and to ease notation let $\Phi _{I}^{x_i} = \{\Phi _{I}^{(1,0)}, \Phi _{I}^{(1,0)}\}$ and $\Phi _{I}^{x_i x_j} = \{\Phi _{I}^{(2,0)}, \Phi _{I}^{(1,1)}, \Phi _{I}^{(0,2)}\}$ represent the first and second order differential operator sets respectively. For any point $p\in \bar{\Omega }$ such that $0<p\le N$ with associated coordinate $x^p_i$ following properties have been numerically verified:

• 1. $\sum _{I=1}^{N} \Phi _{I}^{(0,0)} = 1$,

• 2. $\sum _{I=1}^{N} \Phi (\mathbf{x} ^p)_{I}^{(0,0)}[x_i]_I = x_i$,

• 3. $\sum _{I=1}^{N} \Phi (\mathbf{x} ^p)_{I}^{x_j}[x_i]_I = \delta _{ij}$,

• 4. $\sum _{I=1}^{N} \Phi (\mathbf{x} ^p)_{I}^{x_i x_j}[x_i x_j]_I = 1$    if $i \ne j$,

• 5. $\sum _{I=1}^{N} \Phi (\mathbf{x} ^p)_{I}^{x_i x_j}[x_i x_j]_I = 2$    if $i = j$.

Note that there also exist similar properties for ${\mathbb {R}}^3$ domains and higher ordered derivatives (i.e $m>2$) as well.

Another way to measure the effectiveness of differential operators is through its interpolation error. This may be achieved by benchmarking solution data for a given partial differential equation against data from another known or trusted solution analytical, numerical, or manufactured solution. In this work, we take the latter approach by manufacturing and differentiating two primary fields c(x, y) and $\phi (x,y)$ as shown in Eqs. (38) and (39) respectively to produce a term $\nabla \cdot (\phi (x,y)\nabla c(x,y))$. The manufactured fields are variants of those taken from a helpful report on the method of manufactured solutions [38].

$$\begin{aligned} c(x,y)= &\, {} c_0\left[ 1+\sin ^2\left( \frac{x}{R}\right) \sin ^2\left( \frac{2y}{R}\right) \right] \end{aligned}$$

(38)

$$\begin{aligned}&\quad \phi (x,y) = \phi _0 \left[ \frac{\sqrt{x^2+2y^2}}{R}\right] \end{aligned}$$

(39)

Fig. 11

Comparison of the solution surface for $\nabla \cdot \left( \phi (x,y)\nabla c(x,y) \right)$; the analytical solution surface is derived based on the manufactured fields $\phi (x,y)$ and c(x, y)

Figure 11 presents a comparison between analytical y and approximate solution $y^h$ fields for the term $\nabla (\phi (x,y) D\nabla c(x,y)$ which test all two-field first order derivative products, second order derivatives and combinations thereof; for the computation, an unit square domain is discretized with uniformly distributed 4900 collocation points. The computed discrete L₂ norm error which is given by

$$\begin{aligned} \text {error}_{L_2} = \sqrt{\left( \frac{y-y^h}{y}\right) ^2} \end{aligned}$$

(40)

was measured at a value less than 0.3%. Those who interested in the computational resources to construct differential operators are referred to Sect. 4 of this study.