A Phase Field Model for Binary Alloy Solidification with Boundary Interface Intersection

A phase field model for binary alloy solidification with boundary interface intersection is developed. In the phase field model, the heat and solute conservation equations are appropriately modified to account for the presence of heat and solute rejection inside the diffuse interface, and a relaxation boundary condition for the phase field variable is introduced to balance the interface energy and boundary surface energy in the multiphase contact region. The thin interface asymptotic analysis is applied on the phase field model to yield the free interface problem with dynamic contact point condition. Introduction The study of the evolution of the interface structure in alloy solidification near physical boundary is of great technological interest. The solidification near the boundary interface multiphase intersection is governed by heat and solute flux as well as the boundary interface surface energy balance. To simulate the micro-structure evolution during solidification, the phase field approach has been known as a widely accepted technique[2, 8, 10, 11, 12, 16, 17, 18, 20]. This approach has the well known advantage that it avoids explicitly tracking a sharp boundary by introducing an order parameter, i.e., a phase field variable, which takes constant values in the bulk phases and varies smoothly but steeply in a diffuse interface region with small thickness. In order to model the coupled heat and solute diffusion in binary alloy solidification, the heat and solute conservation equations are appropriately modified to account for the presence of heat and solute rejection inside the diffuse interface in the phase field model[17]. However, the above mentioned work didn’t consider the case when the interface intersect with the external physical boundary of the domain, where additional mechanism to balance the interface energy and boundary surface energy in the multiphase contact region should be imposed[19]. To recover the additional mechanism in the multiphase contact region, we consider the solidification of a dilute binary alloy in a channel Ω = {x = (x, y) : x ∈ R, 0 ≤ y ≤ L} with a moving interface Γ(t) separating solid and liquid states, Ω+(t) and Ω−(t), where Ω+(t) denotes the liquid region, Ω−(t) denotes the solid region and Γ(t) is the solid/liquid interface. The evolution of the interface is extensively studied when the interface is a closed curve that never intersect the external boundary. We will concern the case when Γ(t) intersect ∂Ω at {y = 0} and {y = L}, and study the dynamics of the boundary-interface intersection. The motion of the interface is controlled by coupled heat and solute diffusion in the whole domain. We assume the alloy have equal thermal diffusivity α in solid and liquid, the solutal diffusivity is D in liquid and zero in solid. The governing free-interface problem is ∂tU = D∇U (liquid), (1) ∂tθ = α∇θ (liquid and solid), (2) (1 + (1− k)Ui)Vn = −D∂nU |l (mass conservation), (3) Vn = α(∂nθ|s − ∂nθ|l) (heat conservation), (4) International Frontier Science Letters Submitted: 2016-03-15 ISSN: 2349-4484, Vol. 8, pp 9-18 Revised: 2016-04-15 doi:10.18052/www.scipress.com/IFSL.8.9 Accepted: 2016-04-15 © 2016 SciPress Ltd., Switzerland Online: 2016-06-29 SciPress applies the CC-BY 4.0 license to works we publish: https://creativecommons.org/licenses/by/4.0/ θi +Mc∞Ui = −d0κ− βVn (Gibbs-Thomson), (5) where U is the dimensionless concentration (mole fraction of solute) and θ is the temperature relative to the equilibrium freezing temperature scaled by latent heat of fusion and the heat capacity, k is the equilibrium partition coefficient, Vn is the normal velocity of the interface, n is the direction normal to the interface pointing from solid to liquid, M is the scaled magnitude of the liquidus slope, c∞ is the initial concentration of the alloy, d0 is the thermal capillary length, κ is the local curvature of the interface, β is the kinetic coefficient. Note that the subscript i on θ and U is present when these quantities are being evaluated at the interface. In this setting, we need to solve the coupled heat and solute diffuse equations, and the shape and position of the free boundary Γ(t). The evolution of the interface Γ(t) is governed by the above kinetic Gibbs-Thomson relation (5), this equation can be understood similarly as curvature flow equation, which is a second order parabolic type partial differential equation. When the free boundary Γ(t) intersect the external boundary, to complete the formulation, we need to impose conditions on the boundary-interface contact points to incorporate the boundary-interface interaction in the evolution of the interface, thus additional dynamic conditions on the contact points are required[19]. To get some insight on the dynamic contact point condition, we study this problem by phase field model, in which the sharp interface is replaced by a diffuse interface, thus an auxiliary order parameter φ, the phase field, is introduced to indicate the phase. The quantityφ is a continuous variable that takes constants in the bulk phases, say +1 in the solid and−1 in the liquid, and varies smoothly but steeply from−1 to +1 over a thin layer, the diffuse interface. Then a phase field model will be constructed to govern the time evolution of φ, U and θ, which incorporate the interfacial physics of the problem and the diffuse of heat and solute in the system. The standard phase field model, a version of the Model C of Halperin, Hohenberg and Ma [9] was originally derived on the basis of a gradient functional for the free energy of the system, with dynamics that ensure the free energy decreasing with time. The time rate of change of the order parameter is coupled to the energy equation to incorporate the release of latent heat as solidification occurs. This model was further elaborated in Caginalp [3], Fix [8] and Langer [14] etc.. However, the phase field model has the disadvantage that it is hard to use quantitatively because it is often computationally too stringent to choose the interface thickness small enough to resolve the desired sharp interface limit of the phase field model. To be useful for real simulation, phase-field models need to bridge the scale gap between the thickness of the physical solid-liquid interfaces and the typical scale of the microstructures. This is achieved by increasing the orders of magnitude of the interface width such that the width of the diffuse interface is about one order of magnitude smaller than the radius of curvature of the interface but much larger than the real microscopic width of a solid-liquid interface. This procedure magnifies the physical effect that is due to the diffuseness of the interface, therefore, to guarantee precise simulations, all these effects have to be controlled or, if possible, eliminated. The privileged tool to achieve this is the so-called thin-interface limit presented by Karma and Rappel [11, 12], Almgren [1] for pure substances, and by Karma [10] for isothermal solidification of dilute binary alloys: the equations of the phase field model are analysed under the assumption that the interface thickness is much smaller than any other physical length-scale present in the problem, but much larger than the capillary length. The procedure of matched asymptotic expansions then yields the effective interface conditions valid at the macroscale, which contain all effects of the finite interface thickness up to the order to which the expansions are carried out. Note that the isothermal solidification of dilute binary alloys considered in Karma [10], the socalled antitrapping current was introduced to eliminate nonequilibrium effects at the interface, and it was shown that quantitative simulations of alloy solidification are possible with this model [5]. However, the concentrations in dilute binary alloys is small, the thermal diffusivity may be much larger than the solutal diffusivity, the solutal and thermal effect are still comparable and the solute and heat conservation equations need to be solved simultaneously. For this purposes, a phase-field model 10 Volume 8


Introduction
The study of the evolution of the interface structure in alloy solidification near physical boundary is of great technological interest. The solidification near the boundary interface multiphase intersection is governed by heat and solute flux as well as the boundary interface surface energy balance. To simulate the micro-structure evolution during solidification, the phase field approach has been known as a widely accepted technique [2,8,10,11,12,16,17,18,20]. This approach has the well known advantage that it avoids explicitly tracking a sharp boundary by introducing an order parameter, i.e., a phase field variable, which takes constant values in the bulk phases and varies smoothly but steeply in a diffuse interface region with small thickness. In order to model the coupled heat and solute diffusion in binary alloy solidification, the heat and solute conservation equations are appropriately modified to account for the presence of heat and solute rejection inside the diffuse interface in the phase field model [17]. However, the above mentioned work didn't consider the case when the interface intersect with the external physical boundary of the domain, where additional mechanism to balance the interface energy and boundary surface energy in the multiphase contact region should be imposed [19].
To recover the additional mechanism in the multiphase contact region, we consider the solidification of a dilute binary alloy in a channel Ω = {x = (x, y) : x ∈ R, 0 ≤ y ≤ L} with a moving interface Γ(t) separating solid and liquid states, Ω + (t) and Ω − (t), where Ω + (t) denotes the liquid region, Ω − (t) denotes the solid region and Γ(t) is the solid/liquid interface. The evolution of the interface is extensively studied when the interface is a closed curve that never intersect the external boundary. We will concern the case when Γ(t) intersect ∂Ω at {y = 0} and {y = L}, and study the dynamics of the boundary-interface intersection.
The motion of the interface is controlled by coupled heat and solute diffusion in the whole domain. We assume the alloy have equal thermal diffusivity α in solid and liquid, the solutal diffusivity is D in liquid and zero in solid. The governing free-interface problem is where U is the dimensionless concentration (mole fraction of solute) and θ is the temperature relative to the equilibrium freezing temperature scaled by latent heat of fusion and the heat capacity, k is the equilibrium partition coefficient, V n is the normal velocity of the interface, n is the direction normal to the interface pointing from solid to liquid, M is the scaled magnitude of the liquidus slope, c ∞ is the initial concentration of the alloy, d 0 is the thermal capillary length, κ is the local curvature of the interface, β is the kinetic coefficient. Note that the subscript i on θ and U is present when these quantities are being evaluated at the interface.
In this setting, we need to solve the coupled heat and solute diffuse equations, and the shape and position of the free boundary Γ(t). The evolution of the interface Γ(t) is governed by the above kinetic Gibbs-Thomson relation (5), this equation can be understood similarly as curvature flow equation, which is a second order parabolic type partial differential equation. When the free boundary Γ(t) intersect the external boundary, to complete the formulation, we need to impose conditions on the boundary-interface contact points to incorporate the boundary-interface interaction in the evolution of the interface, thus additional dynamic conditions on the contact points are required [19].
To get some insight on the dynamic contact point condition, we study this problem by phase field model, in which the sharp interface is replaced by a diffuse interface, thus an auxiliary order parameter φ, the phase field, is introduced to indicate the phase. The quantity φ is a continuous variable that takes constants in the bulk phases, say +1 in the solid and −1 in the liquid, and varies smoothly but steeply from −1 to +1 over a thin layer, the diffuse interface. Then a phase field model will be constructed to govern the time evolution of φ, U and θ, which incorporate the interfacial physics of the problem and the diffuse of heat and solute in the system.
The standard phase field model, a version of the Model C of Halperin, Hohenberg and Ma [9] was originally derived on the basis of a gradient functional for the free energy of the system, with dynamics that ensure the free energy decreasing with time. The time rate of change of the order parameter is coupled to the energy equation to incorporate the release of latent heat as solidification occurs. This model was further elaborated in Caginalp [3], Fix [8] and Langer [14] etc.. However, the phase field model has the disadvantage that it is hard to use quantitatively because it is often computationally too stringent to choose the interface thickness small enough to resolve the desired sharp interface limit of the phase field model. To be useful for real simulation, phase-field models need to bridge the scale gap between the thickness of the physical solid-liquid interfaces and the typical scale of the microstructures. This is achieved by increasing the orders of magnitude of the interface width such that the width of the diffuse interface is about one order of magnitude smaller than the radius of curvature of the interface but much larger than the real microscopic width of a solid-liquid interface. This procedure magnifies the physical effect that is due to the diffuseness of the interface, therefore, to guarantee precise simulations, all these effects have to be controlled or, if possible, eliminated. The privileged tool to achieve this is the so-called thin-interface limit presented by Karma and Rappel [11,12], Almgren [1] for pure substances, and by Karma [10] for isothermal solidification of dilute binary alloys: the equations of the phase field model are analysed under the assumption that the interface thickness is much smaller than any other physical length-scale present in the problem, but much larger than the capillary length. The procedure of matched asymptotic expansions then yields the effective interface conditions valid at the macroscale, which contain all effects of the finite interface thickness up to the order to which the expansions are carried out.
Note that the isothermal solidification of dilute binary alloys considered in Karma [10], the socalled antitrapping current was introduced to eliminate nonequilibrium effects at the interface, and it was shown that quantitative simulations of alloy solidification are possible with this model [5]. However, the concentrations in dilute binary alloys is small, the thermal diffusivity may be much larger than the solutal diffusivity, the solutal and thermal effect are still comparable and the solute and heat conservation equations need to be solved simultaneously. For this purposes, a phase-field model 10 Volume 8 using a computationally tractable thin-interface analysis that makes it possible to obtain results that are independent of the interface width and that provides the freedom to choose the phase-field parameters such that growth without kinetic effects or solute trapping were set forth in [17].
In this paper, we extend the phase field model in [17] to incorporate the interaction of external boundary and the interface. To recover this surface-interface energy balancing mechanism in the boundary interface intersection, we will construct a phase field model based on a modified free energy functional proposed by Cahn [4] to include the boundary-interface contact energy. Thus a modified gradient flow with relaxational boundary condition were derived. We further illustrate by asymptotic analysis that the relaxational boundary condition leads to a dynamic contact angle relation at the boundary-interface intersection. Compared with the Neumann conditions proposed in [13] or the stationary contact angle conditions proposed in [6], this model with relaxational boundary condition gives the dynamic nature of the contact angle relation in the study of boundary-interface interaction in binary alloy solidification.
In Section 2, we will describe the new phase-field model. Asymptotic analysis of this model and the leading order interface dynamics are presented in Section 3.

Phase field model
We briefly indicate the construction of the phase field model that reduces to the above set of freeinterface equations. In our case of solidification with boundary-interface intersection, the associated free energy functional is assumed to be of the Ginzburg-Landau form with an additional surface energy term in order to account for the boundary-interface interaction [4]: where the small parameter ε is gradient energy coefficient which also characterize the thickness of the interface, f AB (φ, U, θ) is the bulk free energy density of a binary mixture of A and B atoms/molecules and U denotes the dimensionless solute concentration defined as the scaled mole fraction of B, γ(φ) is the surface free energy density function on the external boundary. The square-gradient term is the contribution of the interface. The integrations are understood as in some suitable finite domain under consideration. By variational principle, the dynamical equation for φ can be given by with the relaxational boundary condition where τ, η are the mobilities of the interface and external boundary, respectively, δ is the interfacial energy per unit area, l is the outward normal direction to ∂Ω.
To specify the free energy density, it is worthwhile to keep φ fixed in the bulk phase thus a form of the free energy density including the pure and solute contribution is chosen as here p(φ) = (φ 2 − 1) 2 /4 is a double well potential, the conventional dimensionless parameter λ is the ratio of the interfacial width to capillary length and g(φ) is a smooth, odd function of φ that monotonically increase from g(−1) = −1 to g(1) = 1 in order to incorporate the release of latent heat as solidification occurs.
The equations of heat and solute diffusion are modified by incorporating time rate change of the order parameter φ, thus Equation (10) is equivalent to the mass continuity relation. The gradient solute current consists of two parts. The first part is the standard flux driven by the gradient of chemical potential that reduces to Fick's law of diffusion in the liquid, the secondj at is an "antitrapping" current that is only nonvanishing in the diffuse interface region [10,17], and hence does not affect the bulk thermodynamic properties of the model. This current produces a solute flux from solid to liquid along a direction normal to the interface that counterbalances the solute flux induced by the gradient of chemical potential across a moving interface. The "antitrapping" current can be chosen as The equations (7)- (12) lead to the following final set of phase field model ∂θ ∂t = α∇ 2 θ + 1 2 ∂φ ∂t (15) with the relaxational boundary condition (8). This problem may be solved by properly imposed initial and boundary conditions for U and θ.
Next, we carry out the asymptotic analysis for the phase field model (13)- (15) and derive the leading order behavior, with ε ≪ 1 being assumed, whereby an outer solution valid away from the interface is matched to an inner solution that is valid in the interfacial region [7].
For ε > 0, we define the interface Γ(t; ε) to be the level set {φ = 0}. Our primary attention is on the motion of the interface and its intersection with external boundary. Denote a point on the interface by x 0 (s, t; ε), where s is the arc-length measured from one intersection point. For each point x in a neighborhood of Γ(t; ε), we defined a signed distance function d(x, t; ε) with the same sign as φ(x, t; ε) and the value equal to the distance from x to its projection point x 0 on Γ(t; ε). With those notations, each point (x, t) in a neighborhood of Γ(t; ε) can be prescribed by a new local orthogonal curvilinear coordinate system (d, s, t). The inner expansion for such a point is then described via a stretched variable ρ = d/ε. With a slightly abuse of notations, for each ε, the functions φ, U and θ may be viewed as functions either of (x, t) or of (ρ, s, t). Elementary computations gives us the formula for converting derivatives w.r. t. (x, t) to derivatives w.r. t. (ρ, s, t) as and we further note that ∇d = (d x , d y ), |∇d| = 1. △d = κ is the curvature of the level surface of d, normal translations of Γ. d t = −V is the negative normal velocity of the interface. This change of coordinates is valid near Γ, in particular it is valid within a distance O(ε).

The outer expansions
We firstly consider the outer expansion away from the interface with thus the leading order heat and solute diffuse equations in the outer regions are which recovers Equation (1)-(2) in the bulk region. To solve the free interface problem, we need further to find the transmission conditions on Γ. These will be derived from the inner problem near the interface.

The inner expansions
In the inner region near a point x 0 on the interface, we use the stretched variable ρ = d/ε and expand φ(x, t; ε) = Φ(ρ, s, t; ε) = Φ 0 (ρ, s, t) + εΦ 1 (ρ, s, t) + O(ε 2 ), the resulting equations are solved order by order in ε, with far field boundary conditions that are obtained by matching to the outer solutions. This requires the matching conditions for U . The matching conditions for θ are similar. For φ, we have By our definition that φ = 0 on Γ, we also have the condition that Φ i (0, s, t) = 0 for each i.

Classical asymptotics
We reformulate (13)- (15) in curvilinear coordinates (d, s, t) and rescale d to ρ = d/ε, by using the differential relations (16), to obtain The conventional analysis assumed λ = εΛ to be a small parameter, which requires the interface width assumed to be small compared to the capillary length. This requirement however can be loosen by "thin interface model" [11,12] or "isothermal asymptotics" [1], in which λ is assumed to be of order is still retained, see [15]. We now substitute the inner expansion (20) into (21) and use the matching conditions to get the structure and motion of the interface.

Leading order phase field equation.
By the thin interface assumption ( [10,15]) corresponds to a stationary interface, the leading order phase field equation is this yields for solution at leading order the stationary phase field profile which is monotone decreasing and approaching its limiting value exponentially as ρ → ±∞.

The leading order heat and solute diffuse equations.
Similarly, the leading order heat diffuse equation is Θ 0 ρρ = 0, integrating and using the matching condition to get this result and the thin interface assumption (22) yield both Θ 0 and U 0 are constants independent of ρ.

First order phase field equation. The first order phase field equation is
The solvability condition yields This equation involves first order U and θ thus we need to go further analysis. 4. First order U 1 . Notice that U 0 is independent of ρ, from the second equation in (21), we have .
Take the limit as ρ → +∞ we have which is the mass conservation equation on the interface in (3). One expression for U 1 is 5. First order Θ 1 . Note that Θ 0 is independent of ρ, from the third equation in (21), we have take the integration w.r.t. ρ to get where A is a constant. Pass to the limit as ρ → ±∞ and subtract to have the heat conservation equation on the interface in (4): the expression for Θ 1 is Substitute U 1 in (27) and Θ 1 in (29) into the solvability condition (25): note that Φ 0 ρ is an even function of ρ has been used. By taking the thermal capillary length d 0 and kinetic coefficient β as we have the Gibbs-Thomson condition (5): 6. Boundary-interface intersection. Now let us derive the boundary-interface intersection condition at Γ ∩ {y = 0}. Take the inner expansion of the relaxational boundary condition (8), we obtain from the leading order that where ψ is the dynamic contact angle between the sharp interface Γ and the x-axis. Multiply the above equation by Φ 0 ρ and integrate in ρ over (−∞, +∞), notice the interfacial energy δ = here ψ is the dynamic contact angle between the interface Γ and the x-axis. The Young's equation is recovered by introducing an "equilibrium contact angle" ψ 0 which is defined when the contact point velocity V = 0. The condition at Γ ∩ {y = L} can be similarly derived.

The leading order behavior
We now summarize the results obtained from the above asymptotic analysis, these are Eq. (19) in the bulk, Eq. (26), (28) and (31) on the free interface, together with the contact point condition (32) at Γ ∩ ∂Ω: Note that here the dynamic contact angle ψ is determined as long as the shape and position of the free interface problem are solved. The last condition can be seen as a constraint to the free interface problem when the interface intersect the physical boundary, which implies the dynamic nature of the contact point. If the dimensionless concentration function u 0 takes value 0, which is the case of solidification of pure substance, the above system is similar to the system considered in [19].

Conclusion
In this paper a phase field model for binary alloy solidification with boundary interface intersection is developed. To balance the interface energy and boundary surface energy in the multiphase contact region, a surface free energy term is incorporated in the free energy functional [4], and a relaxation boundary condition for the phase field variable is introduced. By applying thin interface asymptotic analysis [10], the leading order free interface problem is derived with constraint at the contact point to accommodate the dynamic nature of the contact angle. This problem may be completely solved by further specified initial conditions and boundary conditions on physical boundary.