Thin interface limit for phase-field models of solidification with local mobility correction

TL;DR

This paper introduces a novel method for deriving mobility corrections in phase-field models of solidification, utilizing local interface information to improve accuracy across a wider range of diffusivity conditions.

Contribution

It presents a new approach that removes systematic errors in phase-field kinetics without relying on small interface Peclet number approximations.

Findings

Broader applicability for higher interface Peclet numbers.

Elimination of systematic errors in phase-field kinetics.

Analytical form for mobility corrections derived from local interface data.

Abstract

A new approach is developed to derive an analytical form for mobility corrections in phase-field models for pure material solidification. Similar to the thin interface limit approach (Karma and Rappel, 1996) it seeks to remove systematic errors in the kinetics of phase-field models that arise due to the diffusivity of the interface. The new approach harnesses local information within the diffuse interface instead of an expansion within the interface peclet number $p_e=\eta/\delta$, the ratio of diffuse interface thickness to diffusion length. Therefore it has no need for an approximation of $p_e$ being small. This results in a broader range of applicability of the approach for higher $p_e$.

Figures

Equations

$$\dot{\varphi }$$

$$\dot{T}$$

$$v=M^{\text{phys}}({\sigma }^{\ast }\kappa +\Delta S(T_m-T_{\text{i}})),$$

$$p_e=\frac{\eta v}{\alpha }$$

$${\varphi }_{\text{tw}}(x,t)=\{\begin{array}{ll}1& \text{for }x<vt\\ \frac{1}{2}+\frac{1}{2}\cos (\frac{\pi }{\eta }(x-vt))& \text{for }vt\le x<vt+\eta \\ 0& \text{for }vt+\eta <x\end{array}$$

$$v=M^{\text{eff}}\Delta S(T_m-T).$$

$$-v{\partial }_{x}T=\alpha {\partial }_x^{2}T-vk{\partial }_x{\varphi }_{\text{tw}}.$$

$${\partial }_{x}T$$

$$\displaystyle=-\frac{v}{\alpha}T+\frac{v}{\alpha}\Bigl{(}T_{s}-\frac{k}{2}\Big{[}1-\cos(\frac{\pi}{\eta}x)\Bigr{]}\Bigr{)}.$$

$$\displaystyle\partial_{x}C(x)=\frac{v}{\alpha}\Bigl{(}T_{s}-\frac{k}{2}\Bigl{[}1-\cos(\frac{\pi}{\eta}x)\Bigr{]}\Bigr{)}\exp(-\frac{v}{\alpha}x).$$

$$T_R(x)=T_{\text{s}}-v{A}_v(x),$$

$$\displaystyle A_{v}(x):=\frac{k}{2v}\Bigl{(}1-\frac{1}{\pi^{2}\alpha^{2}+\eta^{2}v^{2}}\Bigl{[}\pi^{2}\alpha^{2}\exp(-\frac{v}{\alpha}x)$$

$$\displaystyle+\pi\alpha v\eta\sin(\frac{\pi x}{\eta})+\eta^{2}v^{2}\cos(\frac{\pi x}{\eta})\Bigr{]}\Bigr{)}.$$

$$-v\frac{\partial \varphi }{\partial x}=-M(x)\frac{\partial \varphi }{\partial x}\Delta S(T_m-T_R(x)).$$

$$v=\frac{M(x)}{1-M(x)\Delta S{A}_v(x)}\Delta S(T_m-T_{\text{i}}),$$

$$M(x)=\frac{M^{\text{phys}}}{1+M^{\text{phys}}\Delta S{A}_v(x)}$$

$$x=\frac{\eta }{\pi }\text{arcsin}(2\varphi -1).$$

$$M^{\text{eff}}=\frac{M^{\text{phys}}}{1+M^{\text{phys}}\Delta Sk\eta /(8\alpha )}.$$

$$\displaystyle\dot{\phi}_{i}=M^{\text{eff}}\Bigl{(}\sigma^{*}(\frac{\phi_{i+1}+\phi_{i-1}-2\phi_{i}}{\Delta x^{2}}+\frac{\pi^{2}}{\eta^{2}}(\phi_{i}-\frac{1}{2}))+\frac{\pi}{\eta}\sqrt{\phi_{i}(1-\phi_{i})}\Delta S(T_{m}-T_{i})\Bigr{)}$$

$${\dot{\varphi }}_i=$$

$$\displaystyle+M(\phi_{i})\Bigl{(}\frac{\pi}{\eta}\sqrt{\phi_{i}(1-\phi_{i})}\Delta S(T_{m}-T_{i})\Bigr{)}$$

$${\dot{T}}_i=\alpha \frac{T_{i+1}+T_{i-1}-2T_i}{\Delta x^2}+k{\dot{\varphi }}_i.$$

$$\displaystyle v=-\frac{\partial\phi}{\partial t}\Bigl{(}\frac{\partial\phi}{\partial x}\Bigr{)}^{-1}=-\frac{2\Delta x}{\Delta t}(\phi_{i}^{m}-\phi_{i}^{m-1})/(\phi_{i+1}^{m}-\phi_{i-1}^{m}).$$

$$\Delta t<\frac{\Delta x^2}{2\alpha }$$

$$\Delta t<\frac{\Delta x^2}{2{\sigma }^{\ast }{M}^{\text{eff}}}$$

$${\sigma }^{\ast }=\frac{\Delta S\Delta T_{\text{i}}\eta }{\pi {\delta }_{\text{I}}}$$

$$\Delta t={\delta }_s\frac{\Delta x^2}{\alpha }$$

$$M_0={\delta }_s\frac{\Delta x^2}{\Delta t{\sigma }^{\ast }}.$$

$$\dot{\varphi }$$

$$\displaystyle\text{MSD}=\frac{1}{N}\sum_{i=1}^{N}\Bigl{(}\phi_{\text{a}}(x_{i})-\phi^{m}_{i}\Bigr{)}^{2}$$

Full text

1 Introduction

The phase-field method is established as the method of choice when it comes to the simulation of dendritic microstructure evolution and a broad variety of phase transformations, with the prominent example of solidification processes [1]. Its central idea is to describe an interface between two phases implicitly as a transition region of finite thickness, in which an indicator function $\phi$, the phase-field variable, has a smooth transition between the phases. The dynamics of phase transformation can then be described via a partial differential equation for $\phi$ and the property fields under consideration, here the temperature $T$.

When in 1994 the first dendrites were modeled with the phase-field method by Kobayashi [2], the method gained a significant prominence boost and was consequently investigated more thoroughly. It was realized early on, that in the original simulations of Kobayashi the results depend strongly on the numerical discretization [3], as well as on the spatial resolution, which is quantified by the interface peclet number $p_{e}$. For quantitative numerical simulations, however, one needs a remedy for both types of errors.

Removing the first type of errors is an ongoing venture in research until today [4, 5, 6] and requires special treatment of the numerical discretization.

To accomplish the latter one a correction of the systematic errors of the temperature field compared to the sharp interface solution is necessary. This was first done by the so called sharp interface limit developed by Caginalp [7]. Later this approach was improved to the so called thin interface limit, which was established by Karma and Rappell using matched asymptotics in a first order expansion in the peclet number $p_{e}$ [8, 9] and is the present-day state-of-the-art technique.

The local mobility correction approach presented here can be seen as an extension of the thin interface limits by providing a full analytic solution. We also discuss effects of numerical discretization like interface locking. This paper is organized as follows: First we present the theoretical derivation of the local mobility correction, and then test it in numerical simulations. To differentiate between numerical discretization errors and systematic errors due to finite interface width, first simulation results with constant driving force are examined and then a comparison between the new local mobility correction scheme and the established thin interface limit is done.

2 Theoretical Background

In this work a phase-field model, which describes the solidification of a pure material, is investigated. It uses a double obstacle potential and a linearized driving force. Its governing equations are written as

[TABLE]

where $T$ is the temperature field and $\phi$ is the phase-field variable, with $\phi=0$ indicating presence of the liquid phase and $\phi=1$ presence of the solid phase.

The phase-field mobility $M$ is going to be chosen either as a simple numerical parameter $M^{\text{eff}}$, or as a $\phi$-dependent function $M(\phi)$, depending on whether a thin interface limit or a local mobility correction scheme is employed. $\sigma^{*}$ is the interface stiffness, $\eta$ the numerical interface width, $\Delta S$ is the entropy difference between bulk solid and liquid and $T_{m}$ the melting temperature. $\alpha$ is the thermal diffusivity and $k$ is the so called hypercooling temperature, which is related to the materials latent heat $L$, heat capacity $c_{p}$, and density $\rho$ via $k=L/(\rho c_{p})$.

A phase-field model is regarded as quantitatively correct, if it reproduces the kinetic equation

[TABLE]

of a sharp interface model. Since in this work only the 1D-case is investigated, the interface curvature $\kappa$ will be set to zero. The physical mobility $M^{\text{phys}}$ connects the thermodynamic driving force $\Delta S(T_{m}-T_{\text{i}})$, where $T_{\text{i}}$ is the interface temperature, to the interface velocity $v$.

In the thin interface limit approach by Karma and Rappel [9] formulas for the effective interface mobility $M^{\text{eff}}$ are derived, which express it as a function of the desired physical mobility $M^{\text{phys}}$. These formulas assure, that in the limit of vanishing peclet numbers $p_{e}=\eta/\delta\to 0$ solutions of equation (1) also fullfill equation (3). Because of $\delta=\alpha/v$, the interface peclet number

[TABLE]

is also a measure for the interface velocity.

To start let us consider two aspects, that make matching equation (1) and (3) challenging: First of all, the driving force term in the phase-field model $\Delta S(T_{m}-T)$ is locally varying with $T$, whereas in the sharp interface picture the driving force $\Delta S(T_{m}-T_{\text{i}})$ is a well defined single number. Second, due to its locally varying nature, the driving force in the phase-field model has contributions originating from the bulk undercooling and not only the interfacial undercooling. In the present case this would result in higher velocities than physically reasonable. A sketch of this situation is given in figure 1. From this picture it also gets plausible, that in the limit $p_{e}\to 0$ variations in $T$ throughout the interface become small, and the drivingforce $\Delta S(T_{m}-T)$ approaches a constant value throughout the interface. In the case of one dimensional steady state solidification, where the diffusion fluxes into the solid phase vanish, this driving force value will correspond to the temperature $T_{s}$ which is reached at the solid side of the diffuse interface. Hence we identify $\Delta S(T_{m}-T_{s})$ as the physically reasonable driving force of the phase-field model.

In the sharp and thin interface approaches it is shown, that in the limit $p_{e}\to 0$, where $T$ can be assumed to be constant, a solution of equation (1) in one dimension is the traveling wave solution

[TABLE]

with

[TABLE]

For a constant temperature profile it is therefore shown, that by choosing $M^{\text{eff}}=M^{\text{phys}}$ the desired kinetics are reproduced. If the temperature $T$ respectively the driving force is not constant throughout the interface, neither $\phi_{\text{tw}}(x,t)$ will be a solution of the phase-field equation (1), nor will the kinetic relation (3) be reproduced correctly. Similar to the aforementioned approach, our method also relies on the traveling wave solution, but aims to reproduce it in non approximative manner.

The main idea is now to define $M$ in equation (1) as a mobility function $M(\phi)$ in such a way, that the new phase-field equation is still solved by $\phi_{\text{tw}}$ and the corresponding analogon to equation (3) stays valid even for locally varying $T$. This is possible, because $M(\phi)$ evaluates local information via the phase-field variable and compensates local variations of $T$. Therefore we call this approach local mobility correction.

In order to determine the analytical form of $M(\phi)$, we first transform the governing equations (1) and (2) in a moving frame with velocity $v$ via $z=x-vt$ and assume steady state evolution, so that we may replace $\frac{\text{d}}{\text{d}t}$ by $-v\frac{\partial}{\partial z}$. Within the moving frame and inside the interface the traveling wave solution has the form $\phi_{\text{tw}}(z)=\frac{1}{2}+\frac{1}{2}\cos(\frac{\pi}{\eta}z)$. This enables us to first look for a form of the local mobility $M(z)$, which is dependent on the moving frame coordinate $z$ and not on $\phi$. Since the traveling wave solution is monotonous, $M(z)$ can later be expressed as $M(\phi)$ by inverting $\phi_{\text{tw}}(z)$. In the following the moving frame coordinate $z$ will be written as $x$.

Before we can determine the local mobility $M(x)$, we need to know the analytical form of $T(x)$, which is governed by equation (2). Under the above mentioned conditions equation (2) is

[TABLE]

This ordinary differential equation can be solved by first integrating directly with the initial conditions $T(0)=T_{\text{s}}$ and $T^{\prime}(0)=0$, which results in the equation

[TABLE]

Using the Ansatz $T(x)=C(x)\exp(-\frac{v}{\alpha}x)$ we arrive at

[TABLE]

Directly integrating this leads to the temperature profile

[TABLE]

where $A_{v}(x)$ is given by

[TABLE]

Note that $A_{v}$ if finite in the limit $v\to 0$: $\lim_{v\to 0}A_{v}(x)=\frac{k}{2\alpha}(x+\frac{\eta}{\pi}\sin(\frac{\pi x}{\eta}))$.

Inserting $\phi_{\text{tw}}(x)$ and $T_{R}(x)$ in equation (1) results in

[TABLE]

Straightforward rearranging delivers the equation

[TABLE]

where we identified $T_{\text{s}}$ as the interfacial temperature $T_{\text{i}}$ from the sharp interface picture. The prefactor in front of the driving force $\Delta S(T_{m}-T_{\text{i}})$ is the physical interfacial mobility $M^{\text{phys}}$, which is reproduced by our phase-field model. According to this by choosing

[TABLE]

it is guaranteed, that the modified phase-field model reproduces the physical mobility $M^{\text{phys}}$. Note that $A_{v}(x)$ is always positive, as long as $x>0$ holds and this holds since $x\in[0,\eta]$ in the interfacial region. This is reasonable since the locally corrected mobility should always be smaller than the physical one, so that it can correct for the artificially high driving force within the interface.

In order to get a form of the corrected Mobility $M(\phi)$, which is then also valid in the original and not only the moving coordinate frame, and in order to make use of the local mobility correction in real simulations, the $x$-dependence of $M$ is now expressed by the local value of the phase-field variable $\phi$ itself. The inversion of the traveling wave solution delivers

[TABLE]

By inserting (16) into (15), the $\phi$-dependent form $M(\phi)$ of the local mobility correction is easily obtained.

To derive a non local mobility correction in analogy to [9], one can do the analysis above with $M$ defined as an effective numerical mobility $M^{\text{eff}}$. This results in the equation $v\frac{\partial\phi}{\partial x}=M^{\text{eff}}\frac{\partial\phi}{\partial x}\Delta S(T_{m}-T_{R}(x))$. Since this cannot be valid per point, we integrate by $\int_{0}^{\eta}\text{d}x$. Moreover we use the above mentioned limit $v\to 0$ for $A_{v}$, and than by rearranging and identifying the physical mobility as before, we get the equation

[TABLE]

The effective Mobility $M^{\text{eff}}$ corrects for the overestimated driving force like the local Mobility $M(x)$, but only in an average manner.

3 Numerical implementation

3.1 Numerical Methods

In order to test the applicability of local mobility corrections and examine possible benefits, one dimensional phase-field simulations are carried out, both with the local and the non local version of the mobility correction from equations (15) and (17) respectively.

For the numerical implementation of the non local mobility correction the space discretized version

[TABLE]

of the phase-field evolution equation (1) with $M^{\text{eff}}$ chosen according to equation (17) is used. For the numerical implementation of the local mobility correction

[TABLE]

with $M(\phi_{i})$ chosen according to equation (15) is used.

In both cases the space discretized version of the heat conduction equation (2) is

[TABLE]

One important point is the introduction of a constant auxiliary mobility $M_{0}$ in equation (19). This is done, because using the local mobility $M(\phi_{i})$ in front of the stabilization operator $\Delta\phi-\frac{\pi^{2}}{\eta^{2}}(\phi-0.5)$ has shown to cause severe biases in the resulting velocities. Since the stabilization operator pushes the numerical solution towards the traveling wave shape, it appears more appropriate to let it do so in an spatially uniform manner. We will refer to equation (19) as constant stabilization method.

To solve the time evolution an explicit euler scheme is applied: at time $t_{m}=m\Delta t$ the time derivative is approximated via $\dot{\phi}_{i}^{m}\approx(\phi_{i}^{m+1}-\phi_{i}^{m})/\Delta t$, where $m$ and $m+1$ are indicating function values at timepoints $t_{m}$ and $t_{m+1}$.

The calculation of the velocity $v$, which needs to be known in order to calculate the local mobility $M(\phi_{i}^{m})$, is obtained by accessing information of the phase-field variable one time step earlier:

[TABLE]

3.2 Choice of numerical parameters

Model parameters as well as numerical discretization constants have to meet several conditions in order to grant stable simulations. Because of the usage of en explicit time stepping scheme the stability criterion for the laplacians in the phase-field and the temperature field equation [10]

[TABLE]

have to be fullfilled. In order to avoid, that the interfacial region $0<\phi<1$ gets spread by the varying driving forces beyond scales of $\eta$, the condition $\Delta G<\frac{\pi\sigma^{*}}{\eta}$ has proven to be useful. To meet this condition we choose the interfacial energy, which has no physical meaning in the one dimensional problem, to be

[TABLE]

whith $\delta_{\text{I}}=0.8$. $\Delta T_{\text{i}}=T_{m}-T_{\text{i}}$ is the driving force generating interface undercooling.

The smallest lengthscale, which has to be resolved in phase-field simulations, is the interface width $\eta$. The number $n$ of gridpoints distributed across this length scale determines the spatial resolution of the simulation and is typically chosen to be $n=\frac{\eta}{\Delta x}\geq 5$. In this work it is assumed, that, besides all material and process parameters ($\alpha$, $M^{\text{phys}}$, $\Delta S$, $T_{m}$, $\Delta T$), the interface width $\eta$ is given. With knowing $\eta$ and the desired numerical spatial resolution $n$ the gridspacing $\Delta x$ can be easily calculated. The timestep $\Delta t$ is then calculated according to equation (22) via

[TABLE]

with $\delta_{s}=0.2$ to assure numerical stability of the temperature solver. Note that due to condition (24) this already implies numerical stability of the phase-field as long as $p_{e}\leq\pi\delta_{\text{I}}$, whereas for $p_{e}\geq\pi\delta_{\text{I}}$ numerical stability of the latter becomes the more restrictive condition.

The constant $M_{0}$ is chosen such, that in each timestep the dimensionless stabilization operator $\frac{\eta^{2}}{\Delta x^{2}}(\phi_{i+1}+\phi_{i-1}-2\phi_{i})+\pi^{2}(\phi_{i}-\frac{1}{2})$ contributes the highest possible fraction to the phase-field variable $\phi$ without causing numerical instabilities. This can be achieved analogously to equation (25) by setting

[TABLE]

The used simulation parameters were: $T_{m}=1$, $k=0.1$, $\Delta T_{\text{i}}=0.05$, $\eta=1$, $\sigma^{*}=0.01$, $\Delta S=0.6$ and $\alpha=100.0$. The system was initialized by setting $\phi$ as a step function and setting the temperature field to a constant value of $T_{m}-\Delta T_{\text{i}}$ in the solid region and $T_{m}-\Delta T_{\text{i}}-k$ in the liquid region.

4 Results and discussion

4.1 Simulations with constant driving force

Since mobility correction schemes are meant to compensate for the varying driving force in the diffuse interface and not to correct the errors araising from numerical discretization, both types of errors have to be distinguished. For this purpose simulations with constant driving force were performed, in which only the equation

[TABLE]

with constant $\Delta G:=\Delta S\Delta T_{\text{i}}$ was solved. By varying the mobility $M^{\text{phys}}$, which is equal to $M^{\text{eff}}$ in the constant driving force case, the phase-field model is tweaked towards reproducing different theoretical velocities $v_{\text{th}}=M^{\text{phys}}\Delta G$. In figure 2 the difference between actual and theoretical velocity $v-v_{\text{th}}$ is shown against $v_{\text{th}}$. The velocities were measured by tracking the point with $\phi=0.5$. Those results show the known effect, that the numerical discretization is causing a decreased velocity [3]. It gets more severe for lower spatial resolutions and vanishes asymptotically for high spatial resolutions.

To get an idea about how the simulation kinetics are influenced by the constant stabilization method and by different weightings between stabilization operator and driving force term in normal spatial discretization, we compared the reproduced kinetics of both methods.

In the case of normal spatial discretization the driving force and mobility were weighted by a factor $f$ according to $\Delta G\to f\Delta G$ and $M^{\text{eff}}\to M^{\text{eff}}/f$. The theoretically expected velocity for different $f$-factors is the same. For higher $f$ the stabilization operator has less weight relative to the driving force term. Therefore we expect the simulated traveling wave profile to be more distorted. It can be hypothesized that this will lead to errors also in the kinetics.

To test for this hypothesis, the traveling wave profiles were analysed. For this purpose the analytical form of the traveling wave $\phi_{\text{a}}=\frac{1}{2}\cos{\frac{x-x_{0}}{\eta}}+\frac{1}{2}$ was fitted with respect to $x_{0}$ to the numerical solution at several time points. Then the mean squeared deviation between numerical fitted analytical form $\phi_{\text{a}}(x_{i})$ and the actual numerical values $\phi^{m}_{i}$ was calculated

[TABLE]

with $x_{1},...,x_{N}$ belonging to the diffuse interface. The MSD-value is subject to fluctuations in time, but nevertheless it has well defined time averaged values $\overline{\text{MSD}}$, around which fluctuations take place and which deliver a good measure for the reproduction quality of the traveling wave solution. In order to see the correlation between $\overline{\text{MSD}}$ and the relative error in velocity $(v-v_{th})/v_{th}$, the respective values are plotted together in figure 3. Obviously the error in kinetics depends on the error in the traveling wave reproduction, and with reducing the latter one the first one can reduced without higher numerical effort, which is a good argument for the usage of constant stabilization methods. On top of that it can be seen, that the numerical effects are causing greater errors than the changes in interface peclet numbers. Note that for a given factor $f$, given numerical method and given spatial resolution the respective points in figure 3 differ only by their peclet number. Those points are clustering closely to each other.

It should be mentioned, that the constant stabilization methods exhibit artificial behaviour for very small velocities. They tend to trap the traveling wave solution in an artificial stationary solution, where the stabilization operator compensates in every gridpoint the driving force term. However for intermediate velocities this problem doesn’t occur.

4.2 Comparison between local and non local mobility correction

By comparing results for different spatial resolutions $n=\frac{\eta}{dx}$ in figure 4 one important finding can be made: The local mobility correction scheme works almost equally precise as the non local one for poor spatial resolutions like $n=5$. In order to gain an advantage in the precision of kinetics good spatial resolutions like $n=10$ or higher are required. This can be well understood, since we cannot expect from a theory enhanced by local information to deliver better results in simulations, when the local information is numericaly reduced to a minimum.

Moreover this is together with the results in figure 3 a clear indicator, that numerical errors have a stronger impact than errors arising from the diffusivity of the interface. Nevertheless the results obtained with local mobility correction behave more predictable with staying below the constant driving force results, whereas the non local results are strongly influenced by the applied numerical method and discretization in wether their velocities over- or undershoot.

Additionally the local mobility correction scheme is giving more accuracy compared to the non local one for higher peclet numbers. Even for low spatial resolution this effect can be seen. This is expected, because for higher peclet numbers the diffuse interface errors are getting more severe.

5 Conclusion

In this work we developed an extended approach to thin interface limit analyses, which includes local information in the mobility factor and aims to make quantitative agreement between phase-field and sharp interface models possible with higher accuracy for higher peclet numbers. Since the common thin interface limits already mitigate diffuse interface errors very reliably, most of quantitative errors in phase-field simulations are of numerical nature. Consequently our local approach grants significant advantage compared to the common method only under specific conditions: high spatial resolution and moderately high interface peclet numbers. Moreover the local mobility correction scheme as presented here, is applicable for 1D problems and not trivially extendes to 2D or 3D situations, where interface velocities and diffusion lengths can vary along the interface. In the light of recent research, that aims to remove numerical errors in phase-field simulations [5], the local mobility correction might deliver additional benefit in combination with mitigation measures for numerical errors.

Acknowledgments

This work was funded by the Deutsche Forschungsgemeinschaft (DFG) under grants no. Ste1116/10-2.

Declaration of interest - None.