Extension of Jackson-Hunt analysis for curved solid-liquid interfaces

Accepted Manuscript Extension of Jackson-Hunt analysis for curved solid-liquid interfaces E.S. Nani, Britta Nestler PII: DOI: Reference:

S0022-0248(19)30049-1 https://doi.org/10.1016/j.jcrysgro.2019.01.025 CRYS 24935

To appear in:

Journal of Crystal Growth

Received Date: Revised Date: Accepted Date:

22 October 2018 22 January 2019 23 January 2019

Please cite this article as: E.S. Nani, B. Nestler, Extension of Jackson-Hunt analysis for curved solid-liquid interfaces, Journal of Crystal Growth (2019), doi: https://doi.org/10.1016/j.jcrysgro.2019.01.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Extension of Jackson-Hunt analysis for curved solid-liquid interfaces E. S. Nania,∗, Britta Nestlera,b a Institute

of Applied Materials (IAM-CMS), Karlsruhe Institute of Technology (KIT), Geb. 30.48, Strasse am Forum 7, 76131 Karlsruhe, Germany of Digital Materials Science (IDM), Karlsruhe University of Applied Sciences, Moltkestrasse 30, 76133 Karlsruhe, Germany

b Institute

Abstract In a previous paper, 2D Jackson-Hunt (JH) analysis was extended by accounting for the curvatures of the transformation front for the case of binary systems with stoichiometric phase diagrams. In this study, it is generalized to any given phase diagram. That the steady state solute distribution satisfies zeroth order hypothesis is rigorously established in the present analysis as opposed to the previous theory’s mere self-consistency proof. While the errors incurred due to low p´eclet number approximation, once the zeroth order approximation is implemented, were analyzed in the earlier work, the combined errors are estimated in the current article. The solute distribution of the current formulation is shown to be correct upto first order in p´eclet number as is the relative error in steady state velocity dependence on lamellar spacing and undercooling. For the practical cases of growths under an imposed temperature gradient, the formulation can be readily generalized to rapid solidification situations by recoursing to numerical algorithms for the estimation of solid phase fractions. Keywords: A1. Directional solidification; A1. Diffusion; A1. Jackson-Hunt analysis; A1. Curved solidification fronts; A1. Rapid growth situations.

1. Introduction

boundary conditions. However, solving the full-fledged governing equations proves to be quite formidable. Hence, it has beDetermination of structure-property correlations in materials— come customary to introduce a number of approximations into whereby the response of a given material to thermal, mechanthe governing equations and/or the solutions thereof. The low ical and electromagnetic stimuli can be predicted by knowing p´eclet number approximation of neglecting p´eclet number comthe underlying microstructure—is one of the prime focusses of pared to quantities greater than unity [1, 5, 8–13], any of the materials science. For this enables one to engineer materials boundary condition approximations like the eutectic approxiof desired properties by controlling the microstructure during mation introduced by Jackson and Hunt [1] or the tight coufabrication. However, for the latter, proper insights into the efpling approximation suggested by Donaghey and Tiller [2] or fect of process conditions on the microstructural evolution are the zeroth order approximation implemented by Senninger and essential. For example, in case of eutectic alloys that are proVoorhees [11] or a variant of it employed by Zheng et al [6] are duced through the route of directional solidification, process some of the examples. conditions like melt concentration, undercooling, and speed of Another approximation that has invariably been adopted in growth determine whether the morphology would be rod-like or the JH literature until recently is that of simplifying the geomelamellar, oscillatory or non-oscillatory etc. Furthermore, they try of the boundary. The solid-liquid surface which is curved in also regulate the corresponding features’ dimensions and disgeneral is replaced by a planar one for ease of solving the probtribution. Hence, for efficient designing of such alloys, relalem. Geometry being a strong factor in regulating the growth tionships that exist among the former and the latter are highly kinetics, such a simplification limits the applicability of the thedesirable. ‘Jackson-Hunt like analysis’ is a broad term which ory. Various investigators on a number of occasions have idenrefers to any analytical theory that estimates these relationships tified these limitations and have pointed out the need for gener[1–13] following an approach laid out by K.A. Jackson and J.D. alization. Donaghey and Tiller claimed that the theory would be Hunt in their pioneering work [1]. incomplete without the consideration of curvatures of the transThe approach is as follows: the solute distribution in melt formation front [2]. Roy Elliot ascribed the failure of classical during the steady state eutectic growth is obtained in an analytJH analysis in explaining the experimental data of anomalous ical form and is coupled with the generalized Gibbs-Thomson eutectic growth of Al-Si alloys with flake morphology to the relation leading to the desired result. Therefore, the key step in planar interface assumption [14]. Lahiri and Choudhury [12] this theory is the analytical determination of solute distribution. have attributed the deviations between the classical JH predicThis, in principle, can be achieved by solving the partial diftions and their computational results obtained from phase field ferential equation that governs the same subject to appropriate implementations to the planar interface assumption made in the classical theory. It had earlier been shown that the phase field ∗ Corresponding author results are more accurate compared to the JH theory [15]. Email address: [email protected] (E. S. Nani)

Preprint submitted to Journal of Crystal Growth

January 17, 2019

8

Resolving this issue, very recently, the current authors have been able to develop a theory for binary systems in 2-dimensions that is free of this approximation [13]. Among other results of this study, it has been shown that curvature incorporation modifies the steady-state velocity v.s. lamellar spacing curve at a given undercooling in a direction that reduces the typical deviations reported between analytics and phase-field numerics. However, this extension has only been carried out in the restricted setting of stoichiometric phase diagrams, i.e., for binary eutectic alloys in which both the phases constituting the solid are stoichiometric intermetallics. Further, in Ref. [13], the zeroth order approximation of the kind employed by Senninger and Voorhees and the low p´eclet number approximation were adopted. The errors associated with the latter once the former is implemented were estimated, but not the combined errors. Moreover, it was taken for granted that the unapproximated solute distribution satisfies the zeroth order hypothesis, i.e., that only the constant and the plane wave component of the Fourier series of unapproximated solute distribution are of zeroth order in p´eclet number and all the other components are of higher order. Satisfying this hypothesis is a necessary requirement on the part of the unapproximated solute distribution for the theory to be valid, but it was not rigorously established in Ref. [13]. The objective of the current paper is to present a formulation that overcomes all the foregoing shortcomings of the previous theory. The aforementioned order and error estimates play a crucial role in validating the theory. It is owing to these kind of estimates that Donaghey and Tiller could analyze the accuracy of various flat interface JH theories [2]. They claimed that the eutectic approximation employed by Jackson and Hunt is “not valid” and gives rise to “serious errors” when applied to offeutectic freezing or to systems with markedly different partition coefficients. The reason for this is that the error between the actual and the classical JH solute distribution is finite in such cases. They further argued that the tight coupling approximation is a little better in this respect; for the error is of first order in p´eclet number. Finally, the method adopted by Senninger and Voorhees is even better as their solute distribution is identical to that constructed by Donaghey and Tiller which is correct upto first order in p´eclet number and the errors are of second order. It is evident from the foregoing discussion that even in the present case, to substantiate the theory, it is necessary to perform such error analysis. This has been carried out in a rigorous fashion in the current article. However, the approach followed is entirely different from that employed by Donaghey and Tiller. While they used the explicitly constructed solute distribution for error estimation, we use the concepts from regularity theory. The rest of the paper is organized as follows. Section 2 contains the curvature-informed JH like formulation; the solute distribution and the kinetics relation are derived in this section. Section 3 has proof that the unapproximated steady state solute distribution satisfies zeroth order hypothesis; along with this, the accuracy of the approximate solute distribution and kinetics relation are determined in this section. Section 4 critically evaluates the approach and discusses some ready implications. Section 5 concludes the article. The intermediate steps in various

8

Ω

Ω

α

β

fα λ/2

fβ λ

α

Figure 1: Schematic drawing of one repetitive unit of lamellar eutectic solid composed of phases α and β with phase fractions fα and fβ , respectively, and growing into the liquid Ω∞ with a curved transformation front ∂Ω.

derivations are relegated to the supplementary material [16]. 2. Extension of JH theory for curved surfaces A Jackson-Hunt like analysis that takes into consideration the curvatures of solid-liquid interfaces is presented in this section for binary systems in two dimensions. Such an extension has been carried out for the restricted case of stoichiometric phase diagrams in a previous paper [13]; here, it is generalized to any given phase diagram. We begin by considering the isothermal directional growth of a eutectic alloy at the expense of liquid ahead of it. Since we are working with two spatial dimensions, the growth proceeds with the eutectic solid assuming a lamellar morphology. A schematic of one repetitive unit of the lamellar structure with a curved transformation front is shown in Fig. 1. The width of this representative unit, i.e., the lamellar spacing is denoted by λ. The volume fractions of the phases α and β comprising the eutectic solid are denoted by fα and fβ , respectively. The entire system is assumed to be undercooled from the eutectic temperature (T E ) by an amount ∆T . Further, the transformation is assumed to have reached a steady state in which the growth progresses with a constant speed v. The aim of the theory is to establish a relationship among v, λ and ∆T . As mentioned earlier, the Jackson-Hunt approach for accomplishing this is: (i) To determine the steady state solute distribution in the melt in an analytical form. (ii) To use the above obtained solute distribution in conjunction with the generalized Gibbs-Thomson relation. 2.1. Determination of solute distribution The solute distribution in the melt during the steady state can be evaluated by solving the boundary value problem (BVP)

2

Before carrying on with approximating the BVP, we take a short interlude to study a relation satisfied by the Fourier coefficients of the unapproximated solute distribution which will be used later for improving the accuracy of the approximated one. Consider the (undetermined) Fourier series representation of the unapproximated solute distribution, i.e., the solution of the BVP given by Eqs. (1), (3), (4) and (6):

that governs the same. In a co-ordinate frame attached to the interface, this BVP is given by [13] ! ! ∂2 c ∂2 c vλ ∂c + + = 0, (1) D ∂y ∂x2 ∂y2 ! ∂c→ ds vλ −n · ˆi + ∂c→ −n · ˆj = (cl (x) − c s (x)), (2) ∂x ∂y dx D (x,p(x)) c → 0 as y → ∞, ∂c = 0, ∂x (± 12 ,y)

c = a0 e−y + u(x, y) := a0 e−y +

(3)

am cos(km x)e−qm y ,

(7)

m=1

(4)

where  qm = + 2

where D is the diffusivity of the solute atoms in liquid; x, y are the non-dimensionalized coordinates along the axes perpendicular and parallel to growth direction, respectively, with λ as the scaling parameter; p(x) is the non-dimensionalized inter−n is the unit outward normal vector (pointing into face profile; → the solid phases) to the transformation front; c(x, y) is the concentration (mole fraction) of the solute at the point (x, y) minus the far-field composition c∞ ; cl (x) is the concentration of liquid at the point (x, p(x)) on the interface and hence, by definition, cl (x) = c(x, p(x)) + c∞ ; c s (x) is the concentration of the solid at the point (x, p(x)) on the interface; and s(x) is the non-dimensional length of the interface from origin to the point (x, p(x)). Further, the co-ordinate system is positioned in such a way that the configuration is symmetric about the y−axis. Ω∞ and ∂Ω mark the liquid region and solid-liquid interface, respectively, both corresponding to one lamellar width as shown in the Fig 1. In the rest of the article, we will use  to denote the p´eclet number vλ D. c s (x) is related to cl (x) through various solute distribution coefficients [1, 2, 11]. Denoting c s (x) by cν (x) for all those values of x at which solid phase indexed by ν is in contact with the liquid, this relationship is given by cν (x) = Λν cl (x) − Λν cE + cνE ,

∞ X

r km2 +

2 4

and

km = 2πm.

Plugging this into the boundary condition, Eq. (6), yields   ds = c∞ + a0 e− p(x) + u − cν (x) dx   ds → − =⇒ ∇u · n = c∞ + u + Λν cE − cνE dx − Λν (c∞ + a0 e− p(x) + u). (8)

−n a0 e− p(x) + ∇u · →

Integrating the above equation over one scaled lamellar width gives 

Z1/2 −1/2

Z1/2   u dx =  c∞ + u + Λν cE − cνE dx −1/2

−

Z1/2

Λν (c∞ + a0 e− p(x) + u)dx

−1/2

=⇒ Z1/2   c∞ = Λν (c∞ + a0 e− p(x) + u) − Λν cE + cνE dx.

(5)

(9)

−1/2

where cνE is the equilibrium concentration of ν phase at eutectic temperature; cE is the eutectic point concentration; and Λν is the solute distribution coefficient corresponding to phase ν as defined in Ref. [11]. Making the above substitution in Eq. (2) modifies the boundary condition to ! ∂c→ −n · ˆi + ∂c→ −n · ˆj ds = ((c + c)(1 − Λν )) ∞ ∂x ∂y dx   +  Λν cE − cνE .

Note that the above equation is nothing but the overall mass balance condition. That is, the average concentration in the solid at the transformation front is same as the far-field composition. This can be obtained directly from the governing equations by integrating Eq. (1) on Ω∞ and using Eq. (2). The explicit representation of the solution in terms of Fourier series is not necessary to obtain this result. Carrying out the integration on the r.h.s of above equation and rearranging the terms gives h i Z1/2   cνE fν − cE Λν fν Λν u + a0 e− p(x) c∞ = + dx, (10) (1 − Λν fν ) (1 − Λν fν )

(6)

Ideally, the analytical solution of BVP given by Eqs. (1), (3), (4) and (6) should be used in the generalized Gibbs-Thomson relation to obtain the v − λ − ∆T relation. However, there exists no known methods of obtaining the exact analytical solutions for this BVP. Hence, we modify it to get an approximated BVP, and work with the latter’s solution instead. Of course, care has to be taken that the changes we make do not give rise to serious errors. This is ensured in section 3.

−1/2

where Einstein’s convention of summation over repeated indices is employed. That is, Λν fν = Λα fα + Λβ fβ and cνE fν = cαE fα + cβE fβ . This relation among the Fourier components of solute distribution will be leveraged for constructing more accurate approximate solute distribution. We now go back and introduce modifications into the actual BVP. We retain only the terms that are of zeroth order in p´eclet 3

number in the r.h.s of Eq. 6. As u is the Fourier series of the solute distribution apart from the constant and plane wave term, it is of first order in p´eclet number (section 3.1); hence we neglect it. Eq. (6) then becomes −n ∇u · →

h i ds = Λν cE − cνE + (1 − Λν )c∞ − Λν a0 . dx

Table 1: Two components of the approximated problem (AP)

AI ∇2 u =0 −n ds = hΛν cE − cνE i ∇u · → dx h i + H ν cνE fν − cE Λν fν

(11)

If we further apply the low p´eclet number approximation, i.e., replace ‘qm ’s that occur in u on the l.h.s of the above equation with ‘km ’s, then it is equivalent to asking for a series of the form given below that satisfies the above equation. ∞ X

am cos(km x)e−km y .

u → 0 as y → ∞ ∂u =0 ∂x (± 12 ,y) AII

(12)

∇2 u =0 −n ds = H ν (Λν f ) − Λν  ∇u · → ν dx

m=1

Since u in Eq. (11) is now of the form given in Eq. (12), integrating the former over one (scaled) lamellar width results in the vanishment of the l.h.s. Hence, the following relation between c∞ and a0 is obtained. c∞ =

h i  1 cνE fν − cE Λν fν + a0 Λν fν , ν (1 − Λ fν )

u → 0 as y → ∞ ∂u =0 ∂x (± 12 ,y)

(13)

where, once again, Einstein’s convention is made use of. So far, we have employed the zeroth order and the low p´eclet number approximations. This is identical to the approach taken by Senninger and Voorhees [11]. However, from here on, the two formulations are going to diverge. Senninger and Voorhees solved for all the other Fourier coefficients in terms of plane wave coefficient and far-field concentration through the ‘dotting on both sides’ method. Similar step, in principle, could be performed to relate the latter two, but instead, they used the overall mass balance equation for this purpose. We, on the other hand, are going to use the dotting method, i.e., Eq. (13) to first connect them and then evaluate the other coefficients in terms of the independent parameter. Eliminating c∞ in Eq. (11) using Eq. (13) modifies the former to

The approximate solute distribution, cA , can now be written as cA = c∞ + a0 e−y + (u˜a I + a0 u˜a II ).

Note that the relation Eq. (13) between a0 and c∞ has been used in constructing the approximated problem. However, using the same to eliminate c∞ in Eq. (15) would give an approximate solute distribution correct only upto zeroth order in p´eclet number. Instead, if Eq. (13) is used for this purpose, the resulting approximation would be correct upto first order in p´eclet number. This is proved in section 3.2. Hence an approximate solute distribution correct upto first order in p´eclet number that additionally satisfies the overall mass balance condition identically is cA = L +

−n ds = hΛν cE − cνE i + H ν hcνE f − cE Λν f i ∇u · → ν ν dx  ν ν ν + a0 H (Λ fν ) − Λ , (14) ν

ν

(15)

Z1/2 Λν  u˜ I + a  u˜ II + a e− p(x)  a 0 a 0 (1 − Λν fν )

−1/2

dx

+a0 e−y +  u˜a I + a0  u˜a II

(16) h i where L = cνE fν − cE Λν fν /(1 − Λν fν ). Equivalently, expressing the above equation in terms of c∞ instead of a0 , we have, upto the same order,     Z1/2 Z1/2     −y     e  ν II ν I  cA = c∞ +  J 1 − Λ u˜a dx −  Λ u˜a dx K     K

ν

where H stands for (1 − Λ )/(1 − Λ fν ). In summary, to obtain an approximation of solute distribution, we need to solve for a function of the form given in Eq. (12) that satisfies Eq. (14). That is, we ended up at an entirely new BVP which is composed of the Laplace’s equation and Eqs. (3), (4) and (14). We call this BVP, the ‘Approximated Problem’ (AP). This is further split into two BVPs, AI and AII as given in Table 1. Denoting the solution of AP by ua , it can be expressed as ua = ua I + a0 ua II where ua I and ua II are the solutions of AI and AII, respectively. Further, ua I =  u˜a I and ua II =  u˜a II , where u˜a I and u˜a II are the solutions of AI and AII, respectively, when  = 1. Note that when  is replaced by unity, AI and AII do not depend on  anymore, hence, their solutions, u˜a I and u˜a II , are independent of . Thus, we have proved that ua is of first order in p´eclet number.

−1/2

−1/2

+ u˜a I +

J II ˜ K (17)

4

1/2 h i R where J = (1−Λν fν )c∞ − cνE fν − cE Λν fν and K = Λν e− p(x) dx. −1/2

While the first order accuracy of the above solute distribution for any generic front shape is established in a later section, it can be directly verified from Eq. (17) for the case of planar solid-liquid surfaces. In such a case, i.e., when p(x) = 0, Eq. (17) coincides with the explicit expression of solute distribution developed upto first order in p´eclet number by Donaghey and Tiller [2].

ξ1 =

−1/2

β

∆T β = −m hcX −

cEX iβ

+ Γβ hκiβ ,

(18)

β

∆T = −m hcA − c iβ + Γβ hκiβ . E

Substituting Eq. (16) in Eqs. (20) and (21) gives   ∆T =mα b0 + mα a0 ξ0 + ξ2 + ζ0α + ζ2α   + mα ξ1 + ζ1α + Γα hκiα ,   ∆T = − mβ b0 − mβ a0 ξ0 + ξ2 + ζ0β + ζ2β   − mβ ξ1 + ζ1β + Γβ hκiβ ,

ξ0 =

−1/2

−1/2

Λν ua II dx. (1 − Λν fν )

˜ where, ζ˜1α = hu˜a I iα , ζ1β = hu˜a I iβ , Gα = hp(x)iα − hu˜a II iα , Gβ = hp(x)iβ − hu˜a II iβ , and H = 1/ (1 − Λν fν ). Note that the explicitly evaluated solutions of problems AI and AII for the case  = 1, i.e., u˜a I and u˜a II , occur in the kinetics relation Eq. (25). These can be numerically constructed following the approach given in Ref. [13]. However, for the sake of completeness, this approach is outlined in appendix A. We point out that for the case of stoichiometric phase diagrams, i.e., when Λα = 0 and Λβ = 0, the above relation merges into the one obtained for such systems in Ref. [13]. 3. Order and error estimates In the preceding section, the kinetics relation is obtained from an approximated solute distribution. How good the approximation is defines the reliability of the theory. A finite error between actual and approximated solute distributions is referred by Donaghey and Tiller to be a “serious” one. In principle, such finite errors can give rise to relative errors in velocity calculation that grow to arbitrarily large magnitudes as p´eclet number becomes vanishingly small. Fortunately, in the existing JH theories, such diverging errors do not occur. Nevertheless, finite relative errors in velocity evaluation do arise in formulations that assume eutectic approximation like the original JacksonHunt theory and others. It would be highly desirable to have not just the error in solute distribution, but also the relative error in velocity going to zero with p´eclet number. It will be shown in this section that the approximate solute distribution constructed following the approach of the foregoing section fulfills this. For this, we first obtain some order estimates for the unapproximated solution.

(20) (21)

(22)

(23)

where, b0 = L − cE , ζ0α = he− p(x) iα , ζ0β = he− p(x) iβ , ζ1α = hua I iα , ζ1β = hua I iβ , ζ2α = hua II iα , ζ2β = hua II iβ , Z1/2

ξ2 =

Z1/2

  1 mα Γβ hκiβ + mβ Γα hκiα ∆T − β α (m + m ) v= !  !" ! # Γβ hκiβ − Γα hκiα Gβ − Gα λ mβ mα ˜β α ˜ − b + ζ − ζ 0 1 1 D mβ + mα mβ + mα H (25)

(19)

where cX and cY are concentrations of the components X and Y, respectively, in the liquid layer immediately adjacent to the transformation front; cEX and cYE correspond to the concentration values at the eutectic point; mα is the slope of the α−liquidus at the eutectic point evaluated with respect to component Y; mβ is the slope of the β−liquidus at the eutectic point evaluated with respect to component X; Γα and Γβ are the GibbsThomson coefficients of the α−liquid and β−liquid interfaces respectively; ∆T α and ∆T β are the corresponding average undercoolings; the angle brackets h·iα denote an averaging operation on the quantity placed inside it over α−liquid interface in one lamellar width; the angle brackets h·iβ are similarly defined with the averaging carried out over β−liquid interface; κ is the curvature at any point on the transformation front. It has to be emphasized that the components, X and Y, and the phases, α and β, have to be so named on the phase diagram such that mα and mβ are positive. Without loss of generality, we assume that so far we have been evaluating the distribution of the component X. As the concentrations we refer to are actually the mole fractions, we have cY = 1 − cX and cYE = 1 − cEX . Further, as the solidification is assumed to happen at isothermal conditions at an undercooling of ∆T , Eqs. (18) and (19) become ∆T = mα hcA − cE iα + Γα hκiα ,

Λν ua I dx, (1 − Λν fν )

Solving for a0 by eliminating ∆T in Eqs. (22) and (23) results in h i −(mα + mβ )(b0 + ξ1 ) − mα ζ1α + mβ ζ1β + Γβ hκiβ − Γα hκiα a0 = . (mα + mβ )(ξ0 + ξ2 ) + mα (ζ0α + ζ2α ) + mβ (ζ0β + ζ2β ) (24) Finally, by eliminating a0 from Eqs. (22) and (23), and neglecting terms that are of order  2 or higher we obtain the required v − λ − ∆T relation:

2.2. Determination of kinetics (v − λ − ∆T ) relation Once the solute distribution is determined, the kinetics relation can be obtained by inserting the same in the following averaged generalized Gibbs-Thomson relations. ∆T α = −mα hcY − cYE iα + Γα hκiα ,

Z1/2

3.1. Order estimates for the actual solute distribution While employing the zeroth order approximation, it is assumed that the unapproximated solute distribution satisfies the zeroth order hypothesis. That is, that all the Fourier components of the unapproximated concentration field except the plane wave component and the constant term are of positive order in p´eclet

Λν e− p(x) dx, (1 − Λν fν ) 5

number. Such an assumption has also been employed in Ref. [13]; however, it was justified only by way of a self-consistency argument there. Here, a rigorous proof for the same is provided. Note that c of Eq. (7), i.e., the Fourier series representation of unapproximated solute distribution minus the far-field composition satisfies Eq. (6). Further, the far-field concentration and the plane wave coefficient are related to each other through Eq. (10). Thus, u defined in Eq. (7) satisfies the following equation.

UII, respectively. However, we are interested in those solutions of UI and UII that have vanishing plane wave components. This is because the plane wave component of the unapproximated solute distribution is already accounted for in formulating UP. This means, besides satisfying UI and UII respectively, uu I and uu II should further satisfy the following condition Z1/2

u(x, y) dx = 0

∀ y > max {p(x) : x ∈ [−1/2, 1/2]} .

(27)

−1/2

−n ds =(1 − Λν )u +  hΛν cE − cνE i ∇u · → dx

We now show that uu I is of first order in p´eclet number. Consider the partial differential equation of UI. Multiplying it with u and integrating both sides of the resulting equation on Ω∞ gives Z  Z  ∂u 2 ∇ u udV +  u dV =0 (28) ∂y

Z1/2 h i νE E ν ν + H c fν − c Λ fν + H Λν u dx ν

−1/2

  Z1/2     + a0  H ν Λν e− p(x) dx − Λν e− p(x)   

(26)

Ω∞

=⇒ Z

−1/2

Therefore, the Fourier series of unapproximated solute distribution apart from the constant and the plane wave term satisfies a BVP given by Eqs. (1), (3), (4) and (26). We call this the unapproximated problem (UP). As in the case of AP, we split UP into two BVPs, UI and UII, as shown in Table 2.

Ω∞

Ω∞

∂u u dV = ∂y

Z

X ν u dx + 

∂Ω

Z

(1 − Λν )u2 dx

∂Ω

Z

Λν u dx

∂Ω

Z

H ν u dx (29)

∂Ω

(Divergence Theorem) h i h i where X ν = Λν cE − cνE + H ν cνE fν − cE Λν fν . =⇒ Z Z Z Z ∂u (∇u)2 dV −  u dV ≤ X ν u dx + b2 u2 dx ∂y Ω∞ Ω∞ ∂Ω ∂Ω Z + b1 b2 H u2 dx (30)

UI ∂u =0 ∂y −n ds = hΛν cE − cνE i + H ν hcνE f − cE Λν f i ∇u · → ν ν dx 1/2 Z + (1 − Λν )u + H ν Λν u dx

∇2 u + 

∂Ω

n o n o where b1 = max Λα , Λβ , b2 = max 1 − Λα , 1 − Λβ , H = 1/(1−Λν fν ), and the Cauchy-Schwarz inequality is used (twice). =⇒ Z Z Z ∂u (∇u)2 dV +  (2b2 + 2b1 b2 H − 1) u dV ≤ X ν u dx ∂y

−1/2

u → 0 as y → ∞ ∂u =0 ∂x (± 12 ,y) UII

−n ∇u · →

Z

+

Table 2: Two components of the unapproximated problem (UP)

∇2 u + 

(∇u)2 dV − 

Ω∞

Ω∞

Ω∞

∂Ω

(Integration by parts)

∂u =0 ∂y ds =H ν dx

Z1/2

Renaming 2b2 + 2b1 b2 H − 1 as b3 , we have Z Z Z ∂u 2 (∇u) dV + b3 u dV ≤ X ν u dx ∂y

Λν e− p(x) dx − Λν e− p(x)

Ω∞

−1/2

+ (1 − Λν )u + H ν

Z1/2

=⇒ Z

Λν u dx

(∇u) dV − b3 2

−1/2 Ω∞

u → 0 as y → ∞ ∂u =0 ∂x (± 12 ,y)

Ω∞

v u u tZ Ω∞

∂u ∂y

(31)

∂Ω

!2

Z u2 dV

dV Ω∞

≤

v tZ ∂Ω

(X ν )2 dx

Z u2 dx ∂Ω

(∵ b3 > 0 (see supplementary material [16]))

It can be seen that the solution uu of UP can be written as uu = uu I + a0 uu II where uu I and uu II are the solutions of UI and

(and from Cauchy-Schwarz inequality) 6

=⇒ Z

(∇u)2 dV − b3

v u tZ

Ω∞

3.2. Accuracy of approximate solute distribution, cA We now estimate the accuracy of the approximate solute distribution cA given by Eq. (16). Let the actual solute distribution be denoted by cU . It is straightforward to see that cU is similar in form to Eq. (16) and is given by

Z (∇u)2 dV

Ω∞

u2 dV Ω∞

≤

v tZ

(X ν )2 dx

∂Ω

Z u2 dx ∂Ω

cU =L +

(32)

Z1/2

−1/2

+ a0 e−y + uu I + a0 uu II .

It was shown in Ref. [13] that for functions u that satisfy Eq. (27), the following inequality is valid in the domains of the kind we are concerned with (Eq. (A.4) of Ref. [13]). Z Z u2 dV ≤ Y (∇u)2 dV (33) Ω∞

  Λν uu I + a0 uu II + a0 e− p(x) dx ν (1 − Λ fν ) (39)

Hence cU − cA is Z1/2

   Λν I I II II u − u + a u − u dx u a 0 u a (1 − Λν fν ) −1/2     + uu I − ua I + a0 uu II − ua II . (40)

cU − cA =

Ω∞

n o where Y = max 2, 96y2max with ymax = max {p(x) : x ∈ [−1/2, 1/2]}, n o but can be refined to Y = max 2, 40y2max to give more sharper Therefore, to determine the   order of cU − cA , we need the estibounds. mates for u I − u I and u II − u II . In section 2.1, we demonu

v Z tZ √ Z 2 2 ν (X ) =⇒ 1 − b3 Y (∇u) dV ≤  dx u2 dx 

Ω∞

∂Ω

∂Ω

(34) √ Z =⇒ 1 − b3 Y (∇u)2 dV 

Ω∞

v v u Z tZ √ tZ 2 ν (X ) dx 4 (∇u)2 dV u2 dV ≤  2 Ω∞

∂Ω

Ω∞

(35)

=⇒

Z Ω∞

(Integration by parts and Cauchy-Schwarz inequality) v √ √4 uZ  2X Y t 2 (∇u) dV ≤  (∇u)2 dV (36) √  1 − b3 Y

|uu I − ua I | ≤ k 2

Ω∞

(from Eq. (33)) where X =

rR

(X ν )2 dx.

a

and

|uu II − ua II | ≤ k 2 ,

(41)

That is, we have established that the errors incurred in approximating the solute distribution following the algorithm of section 2 are of second order in p´eclet number. However, there is a subtlety here that has to be attended to. While estimating the difference cU − cA , we have used same a0 for both the unapproximate and approximated solute distributions. However, in the Jackson-Hunt approach, a0 is provided by the theory; and from Eq. (24) it is clear that, for a generic phase diagram, it depends on the solutions of the recasted BVPs. As UI and UII are different from AI and AII, their solutions also differ and hence, the plane wave coefficients of unapproximated and approximated solute distributions differ. This means to evaluate the actual difference, the respective a0 s have to be accordingly accommodated. However, this does not affect the order relation Eq. (42) as both the a0 s are of the same generic form as given

√ √4 √ √4  2X Y  2X Y 3 2   ≤ =⇒ kuk ≤ √  √ L (Ω∞ ) 1 − b3 Y 1 − b3 Y (37)

Following the treatment given in section 3 of Ref. [13], from the above equation, we have |u| ≤ k

u

where k is some positive constant. From Eq. (40) and Eq. (41), we thus have |cU − cA | ≤ k 2 (42)

∂Ω

=⇒ k∇ukL2 (Ω∞ )

a

strated that ua I and ua II are of first order in p´eclet number. Further, in section 3.1, it has been shown that uuI and uuII are similarly of first order. Thus, it follows that uu I − ua I and   uu II − ua II are atleast of first order in p´eclet number. More refined estimates are obtained in the subsection for these   current   quantities. Note that uu I − ua I and uu II − ua II are solutions of BVPs given by UI-AI and UII-AII respectively and further satisfy Eq. (27). These latter BVPs which correspond to the difference of the UP and the AP problems are referred to as difference problems DI and DII, respectively, and are given in Table 3. Following an approach similar to that of the preceding subsection and using the fact that uu I is of first order in p´eclet number, it can be shown that the solutions of DI and DII satisfy the following error bounds

(38)

where k is some positive number. Thus we have shown that uu I is of first order in p´eclet number. That uu II is also of same order can be derived following an identical approach. Since uu = uu I + a0 uu II , it is thus proved that the unapproximated solute distribution satisfies the zeroth order hypothesis. 7

cU should be included in the Gibbs-Thomson relations. However, we have used the approximate solute distribution cA . That is, ideally, Eqs. (44) and (45) must have been worked with instead of Eqs. (20) and (21).

Table 3: Two components of the difference problem (DP)

DI ∇2 u = −  −n ∇u · →

∂uu I ∂y

ds =(1 − Λν )uu I + H ν dx

Z1/2

∆T = mα hcU − cE iα + Γα hκiα , β

∆T = −m hcU − c iβ + Γβ hκiβ .

Λν uu I dx

u → 0 as y → ∞ ∂u =0 ∂x (± 12 ,y)

∆T = mα hcA − cE iα + Γα hκiα + mα hcU − cA iα , β

−1/2

Z1/2   − p(x) ν − Λ e − 1 + H Λν uu II dx ν

E

(46) (47)

where vu is the unapproximated steady state velocity; va is the velocity estimated from the formulation of section 2, i.e., the expression given by Eq. (25), mα χα and mβ χβ are the coefficients of a0 in Eqs. (22) and (23), respectively; and ‘G.F’ stands for geometric factor and is given by the square bracket term in the denominator of Eq. (25). From Eq. (42) we therefore have

−1/2

u → 0 as y → ∞ ∂u =0 ∂x (± 12 ,y)

vu − va = O(). vu

in Eq. (24) with identical zeroth order terms. Hence, Eq. (42) holds absolutely. We now turn to address the question why even though Eq. (13) is used for the entire procedure of constructing AP and thus ua I and ua II , it is not used in the construction of cA (that is for eliminating c∞ from Eq. (15)), but Eq. (10) is used instead. This is because, had the former been used, the error, cU − cA , would have been as given in Eq. (43) instead of Eq. (40). That is, the approximate solute distribution would be correct only upto zeroth order in p´eclet number.

Hence, as the zero in velocity occurs at a finite λ, the kinetics relation, Eq. (25), is correct upto first order in p´eclet number. We point out that though the concentration field obtained by using Eq. (13) in Eq. (15) is correct only upto zeroth order in peclet number, the velocity formula resulting from it is first order accurate. 4. Results and discussion



(1 − Λν fν )     + uu I − ua I + a0 uu II − ua II

(45)

Eliminating a0 from terms hcA −cE iα and hcA −cE iβ of Eqs. (46) and (47) gives   mα mβ χβ hcU − cA iα − χα hcU − cA iβ + O( 2 ) vu = va − (48) !
λ mα mβ α α β β m χ + m χ (G.F) D mβ + mα

∂uu II ∇ u=− ∂y Z1/2   −n ds =H ν ∇u · → Λν e− p(x) − 1 dx + (1 − Λν )uu II dx 2

cU − c A =

β

∆T = −m hcA − c iβ + Γβ hκiβ − m hcU − cA iβ .

DII

Z1/2 Λν u I + a u II + a e− p(x) − 1 u 0 u 0

(44)

The associated error can be estimated as follows. The preceding equations can be re-written as

−1/2



E

The following two tasks have been accomplished in the present study. First, a Jackson-Hunt like formulation accounting for the curvatures of the transformation front is developed for an arbitrary binary phase diagram for two spatial dimensions. Secondly, it is rigorously established that the analysis is first order accurate in p´eclet number. That is, the solute distribution, and the steady state velocity dependence on lamellar spacing and undercooling of the current formulation are correct upto first order in p´eclet number. The key trick employed to prove this is to split the BVP that governs the solute distribution into two BVPs in such a way that the consideration of plane wave components in the recasted BVPs is obviated. The advantage offered by such an exchange becomes transparent from a close analysis of the curvature incorporation carried out for stoichiometric systems in a previous study [13]. For such systems, the plane wave component cancels itself out during the stage of determination of solute distribution. This lubricates the error estimation as it

dx

−1/2

(43)

=O s () In the above equation ‘O s ’ is the ‘sharp order’ ( f is said to be of sharp order of g in a limit if f = O(g) and f , O(g)). The fact that the error is of sharp order of  is easily verified by considering the planar interface case. 3.3. Order estimates for relative error in velocity In this subsection, we evaluate the relative error in velocity incurred when the solute distribution is approximated in the manner described in section 2. To obtain the accurate relationship among v, λ and ∆T , the unapproximated solute distribution 8

permits working with functions that satisfy Eq. (27) (Condition (a) of Ref. [13]). For such functions, it is easy to establish the (Poincar´e) inequality Eq. (33) which immediately yields the estimates we wish to obtain. However, this simplification of the plane wave term effectively cancelling out does not naturally arise in case of a generic phase diagram. Hence, we re-express the governing equations and accordingly split the BVP. The zeroth order and low p´eclet number approximations are adopted in the current treatment implying that the analysis is valid only for small undercoolings where the speed of growth is low. An improvement in this regard is not possible in the case of isothermal growth as these approximations are indispensable within the current framework. The primary reason being the inevitable use of numerical methods for the calculation of concentration fields. The computation of steady state solute distribution involves the determination of certain dot products involving exponential functions (appendix A). However, the presence of curvatures precludes the analytical evaluation of these dot products; they have to be computed numerically. If the low p´eclet number approximation is lifted, the exponential terms will contain the unknown p´eclet number thus preventing the numerical evaluation of the integrals. Further, the traditionally employed low p´eclet number approximation simplifies only the exponentials occurring in the higher order components while the plane wave component continues to carry the unknown p´eclet number. The zeroth order approximation helps remedy this situation. It has to be noted that even though given a separate name, the zeroth order approximation is essentially an assumption about the magnitude of the p´eclet number and hence is valid only at small growth rates. Fortunately, for the experimentally relevant cases of growths under imposed temperature gradients, these approximations can be relaxed as will be shown in a while. In its current form, the curvature informed analysis of Section 2 is too deficient to be of any practical avail. This is because the formulation is developed pretending that the shapes of the transformation fronts are available and the far-field compositions are unknown. However, in practice, it is the regulated ambient conditions that are specified and the resultant morphological features’ characteristics like lamellar spacing and phase fractions are expected to be determined by the theory. Even the inverse problem that we currently solve is limited in its outlook as it is not readily transferable to non-invariant eutectic reactions. In such problems, when the phase fractions are provided as opposed to far field compositions of each of the species, there would be insufficient number of equations to determine the unknowns. Thus the present formulation has to be supplemented with some additional analysis of estimation of phase fractions and the calculation of the self-consistent profile shapes. In the case of directional growths with an imposed temperature gradient, the first of these steps can be accomplished by performing a planar analysis similar to the one of Senninger and Voorhees [11] or Lahiri and Choudhury [12]. Once the phase fractions are available, self-consistent interface shapes can be constructed by a straightforward extension of the iterative process devised for stoichiometric systems in Ref. [13] (see appendix B). Similar program can be carried out in the case of isothermal growth,

however, the simplifications in the form of linearizations which were possible in the directional solidification case [11, 12] cannot be exploited. An important caveat about the above procedure is that the far-field composition is re-predicted and updated in each step of the iterative process. It is important that it converges and remains close to the original input value. This is demonstrated through an example in the following. Table 4: Thermo-physical parameters of an example eutectic system.

Parameter TE σαl (α−liquid surface energy) σβl (β−liquid surface energy) σαβ (α − β interfacial energy) D cαE cβE cE θαl (contact angle by α−liquid surface) θβl (contact angle by β−liquid surface) mα mβ Lα (α phase latent heat of fusion) Lβ (β phase latent heat of fusion) Λα Λβ Γα (= σαl T E /Lα ) Γβ (= σβl T E /Lβ )

Value 1305.9 0.41045 0.21392 0.50895 5.5 × 10−10 0.5005 0.667 0.6115 66.02 38.76 2351.518 1685.746 1.957 × 109 1.367 × 109 0.1 0.2 2.739 × 10−7 2.044 × 10−7

Units K Jm−2 Jm−2 Jm−2 m2 s−1 deg deg K K Jm−3 Jm−3 m1 K m1 K

We consider a eutectic system corresponding to the parameter set given in Table 4. The directional growth of this system under a positive temperature gradient at a speed of v = 2.75 × 10−4 m/s is studied employing the current formulation. The results pertaining to growths with minimum front undercoolings are collected in Table 5 for each of the solid phase fractions listed in the table. Thus, for the initial guess which is a planar surface, a far field composition of c∞ = 0.53389 is predicted when the phase fractions are fα = 4/5 and fβ = 1/5. In the later iterations, this value changes as the reconstructed front shapes change before converging eventually. Looking at the results from a different point of view, i.e. supposing that a far field composition of c∞ = 0.53389 is given as the input, it is clear that the planar theory estimates the above volume fractions for the phases in the growing solid. Performing the subsequent calculations with these phase fractions leads to the re-prediction of the far-field composition in each iteration, however, as the converged c∞ is very close to the originally given one, such a method can be applied for practical cases and for non-invariant eutectic reactions. The deviations of the phase fractions from their lever rule counterparts fαE (defined with respect to eutectic temperature) are also tabulated in the fifth column for each iteration. The undercooling of the converged front in each case is compared to its planar counterpart in the sixth column of the table. Results for various other combinations of Λν s are listed in Table 6. The data of Table 4 is the same as considered in Ref. [13] except for the distribution co-efficients which were zeros there due to the stoichiometric nature of the alloys studied. However, this difference is not the underlying reason for 9

Table 5: Reconstructed profiles at various iterations and the associated far-field compositions and undercoolings tabulated for five different phase fractions. In each case, the results corresponding to the characteristic lamellar spacing (λJH ) are reported.

(1/5, 4/5)

Reconstructed interface profiles

p(x)

Phase fractions ( fα , fβ )

Initial guess 1st iteration 2nd iteration 3rd iteration

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

Iteration

c∞

E fα − fα ×10−3

∆T

∆T converged ∆T planar

0th 1st 2nd 4th 6th 8th 10th

0.63317 0.63321 0.63321 0.63321 0.63321 0.63321 0.63321

3.18062 2.92715 2.93505 2.93876 2.93902 2.93904 2.93904

6.42456 5.93616 5.95137 5.95853 5.95903 5.95907 5.95907

0.92754

0th 1st 2nd 3rd 4th 5th 6th

0.61114 0.61118 0.61118 0.61118 0.61118 0.61118 0.61118

2.12712 1.88240 1.88193 1.88218 1.88221 1.88221 1.88221

5.92142 5.29378 5.29258 5.29323 5.29329 5.29329 5.29329

0.89392

0th 1st 2nd 3rd 4th 5th 6th

0.58357 0.58360 0.58360 0.58360 0.58360 0.58360 0.58360

1.08852 0.92266 0.92044 0.92040 0.92040 0.92040 0.92040

5.64531 4.91969 4.90995 4.90979 4.90979 4.90979 4.90979

0.86971

0th 1st 2nd 4th 6th 8th 10th

0.55597 0.55598 0.55598 0.55598 0.55598 0.55598 0.55598

0.18455 0.12810 0.12478 0.12399 0.12395 0.12394 0.12394

5.58815 4.74877 4.69930 4.68759 4.68695 4.68691 4.68691

0.83872

0th 1st 10th 20th 30th 40th 50th

0.53389 0.53388 0.53387 0.53387 0.53387 0.53387 0.53387

-0.51633 -0.45681 -0.42297 -0.42093 -0.42082 -0.42082 -0.42082

5.79709 4.83430 4.28678 4.25389 4.25209 4.25199 4.25198

0.73346

(1/3, 2/3)

p(x)

x Initial guess 1st iteration 2nd iteration 3rd iteration

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

(1/2, 1/2)

p(x)

x Initial guess 1st iteration 2nd iteration 3rd iteration

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

(2/3, 1/3)

p(x)

x Initial guess 1st iteration 2nd iteration 3rd iteration

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

(4/5, 1/5)

p(x)

x Initial guess 1st iteration 2nd iteration 5th iteration 10th iteration 50th iter.

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

x

10

Table 6: Results for various combinations of Λν coefficients and phase fractions

∆T converged

∆T converged ∆T planar

fαE − fα (×10−3 )

planar

0.63370 0.61150 0.58375 0.55600 0.53380

0.0 0.0 0.0 0.0 0.0

6.42921 5.92104 5.64640 5.58756 4.23032

5.96411 5.29639 4.91080 4.68431 5.79511

0.92766 0.89450 0.86972 0.83835 0.72998

0.63237 0.61063 0.58336 0.55605 0.53422

0.63246 0.61072 0.58340 0.55603 0.53410

7.43023 4.66767 2.35013 -0.21460 -0.18095

6.42462 5.92143 5.64539 5.58825 5.79714

5.95832 5.29207 4.90905 4.68907 4.27380

0.92742 0.89371 0.86957 0.83909 0.73723

0.63239 0.61060 0.58323 0.55580 0.53381

0.63250 0.61071 0.58333 0.55586 0.53385

7.24588 4.73581 2.53195 0.83823 -0.29723

6.41258 5.91497 5.64450 5.59181 5.80434

5.94584 5.28738 4.91061 4.69657 4.30076

0.92721 0.89390 0.86998 0.83990 0.74096

(Λα , Λβ )

( fα , fβ )

(0.0, 0.0)

(1/5, 4/5) (1/3, 2/3) (1/2, 1/2) (2/3, 1/3) (4/5, 1/5)

0.63370 0.61150 0.58375 0.55600 0.53380

(0.4, 0.5)

(1/5, 4/5) (1/3, 2/3) (1/2, 1/2) (2/3, 1/3) (4/5, 1/5)

(0.1, 0.5)

(1/5, 4/5) (1/3, 2/3) (1/2, 1/2) (2/3, 1/3) (4/5, 1/5)

input

c∞ converged

the larger deviations between ∆T converged and ∆T planar observed here. Rather, as pointed out in appendix B, a simplification assumed in Ref. [13] is; the results of Table 6 confirm this.

kinds of phase diagrams [3]. In the current paper, however, we showed that not limited to any particular phase diagram, functions of the form c∞ + a0 e−y + uu I + a0 uu II where uu I and uu II are the solutions of BVPs UI and UII respectively, always satisfy the governing equations when c∞ is chosen appropriately with regard to a0 , uu I and uu II . The existence of the solutions uu I and uu II follows immediately from inequality Eq. (33) and LaxMilgram theorem as shown in Ref. [13]. As a0 is a free variable, we have that the solution to the BVP is non-unique until both the far-field composition and phase fractions are specified or an extra hypothesis like that of an isothermal interface is invoked. It has to noted that in stoichiometric phase diagrams the extra condition is necessary for the unique determination of the solution [13].

Finally, as the speed of advancement of the transformation front is apriori specified in the directional growth case, the p´eclet number is no longer an unknown. Hence, numerical integrations can be performed on integrands carrying the p´eclet number; therefore, the low p´eclet number and the zeroth order approximations can be relaxed. However, an algorithm like the one of Ludwig and Leibbrandt [7] has to be implemented in the initial step for obtaining the phase fractions. The above discussion implies that the growth of stoichiometric eutectics under an imposed temperature gradient is the only problem that can be solved accurately with no restriction on the values of p´eclet number, the details are provided in appendix C.1. The adjusted recipe for non-stoichiometric phase diagrams is presented in appendix C.2. A by-product of the study is that it is definitively proved that in any generic non-stoichiometric phase diagram, the BVP that governs the steady state solute distribution is not sufficient by itself to determine the solute distribution uniquely if only one among the far-field composition and phase fractions is specified. Either both have to be provided or an additional condition like that of an isothermal transformation front has to be imposed for this purpose. Though Donaghey and Tiller claimed to have shown this, they only set up an infinite system of linear equations to determine Fourier coefficients but did not demonstrate if it can be solved for multiple combinations of far-field composition and phase fractions [2]. Trivedi et al are the first ones to show the non-uniqueness, albeit, for the special cases of two

5. Summary and conclusion Two fundamental differences in the governing equations preclude a straightforward extension of curvature incroportion carried out for stoichiometric phase diagrams to non-stoichiometric systems. One, a lack of possibility of direct estimation of phase fractions from the far-field composition which otherwise allows numerical computations that are essential for the determination of solute distribution. Two, the absence of a simplification brought in by the effective cancellation of plane wave component in the boundary condition which otherwise facilitates an easy error estimation. The second of these is tackled by decomposing the BVP appropriately and expressing the error in terms of those associated with the daughter BVPs. While for the first problem, a planar analysis is performed in the initial step to predict phase fractions whose validity for the reconstructed fronts of the later steps is verified a posteriori by showing that the repredicted far-field composition differs little from the input value 11

(only in the fourth decimal place for Λν s as high as 0.5). Thus, a first order accurate, curvature informed JH-like formulation is developed for arbitrary binary phase diagrams. For the practical cases of directional growths regulated through an imposed positive temperature gradient, the analysis is readily generalizable to non-invariant reactions and rapid growth conditions. The entire range of compositions in the miscibility gap can be dealt within the scope of the current theory. However, it is well known that when the departure from eutectic composition is sufficiently high, oscillations set in into the lamellar structures at large enough spacings [17]. A theory that accounts for such solid-solid interface curvatures along with the currently considered solid-liquid surface curvatures is the objective of future investigations. Further, as the phase fractions of the solidifying eutectic become increasingly dissimilar, the growth morphology transitions into a ‘rods embedded in a matrix’ kind of pattern. A unified theory that handles both lamellar and rod growth dynamics while taking into account the surface curvatures is yet to be developed.

integrating the resulting equation over the interval (−1/2, 1/2) gives the following equation ∞ X

am Amn = wn

m=1

where wn stands for wn =

1/2 R

w(x)φn (x)dx. Note that corre-

−1/2

The authors gratefully acknowledge the financial support by the German Research Foundation through project no. NE822/221.

sponding to each natural number ‘n’, we have an equation of the above type. Thus we have an infinite system of equations linear in ‘am ’s for ‘m’ a natural number. It is observed that the finite reductions of this infinite system has solutions that give an approximation to the function w(x) when plugged into the truncated series of Eq (A.1). Further, the approximation is seen to get better as the number of terms in the truncation increases [13]. The ‘am ’s thus determined are substituted in the truncated series of Eq. (12) to obtain an approximation of the solution of the BVP. The above procedure when applied to BVPs AI and AII for  = 1 gives the approximations to u˜a I and u˜a II which are required to compute the kinetics relation. It has to be emphasized that the dot products Amn cannot be evaluated analytically in a closed form except for a very few shapes p(x) of the interface profile. As a result, numerical methods have to be employed to deal with the typical front shapes found in experimental or simulated microstructures.

Appendices

B. Construction of self-consistent profiles

Acknowledgments

I

II

A. Computation of u˜a and u˜a

For determining the kinetics relation, the fulfillment of the Gibbs-Thomson relation on an average (for each phase) only is demanded. The freedom available to extend this to a pointwise basis is the key idea that is exploited for constructing the self-consistent profile shapes. For this, in the case of directional solidification and in the presence of phase fractions, the planar surface is chosen as an initial guess and the solute distribution corresponding to it is obtained as per the prescription of section 2.1. Using this, the unknown a0 and the front undercooling ∆T are evaluated from the averaged Gibbs-Thomson relations as in section 2.2. These are then substituted into the following pointwise Gibbs-Thomson relations where the just obtained solute distribution is utilized but the curvature term κ is treated as an unknown to be determined.

In this appendix, the method of constructing solutions to BVPs AI and AII is summarized. As both BVPs are the same except for the explicit function that goes into the r.h.s of the boundary condition corresponding to boundary ∂Ω, we give the method for a generic r.h.s of the form wν (x). The Fourier series representation of the solution with undetermined coefficients is given by Eq. (12). Substituting it into the boundary condition at ∂Ω gives ∞ X


am km sin(km x) p0 (x) − cos(km x) e−km p(x) = m=1 −wν (x)

=⇒

∞ X m=1

am φm (x) :=

∞ X

∆T = mα (c∞ − cE + a0 e− p(x) + ua (x, p(x))) + Γα κ,

am sin(km x)e−km p(x)

∆T = −mβ (c∞ − cE + a0 e− p(x) + ua (x, p(x))) + Γβ κ.

m=1

= −

Z

x

wν (x)dx =: w(x).

Integrating the above equations yields

(A.1)

0

Zx

Determining the coefficients ‘am ’ is equivalent to constructing the solution. A method that was empirically observed to be successful in approximating the coefficients is the following. Taking the dot product of the above equation with φn (x), i.e., multiplying both sides of the above equation with φn (x) and

xtr.jn.

  λ  ∆T − mα c∞ − cE + a0 e− p(x) + ua (x, p(x)) dx Γα −p0 (x) − sin θαl = q , 1 + (p0 (x))2

12

(B.1)

Zx xtr.jn.

where

  λ  ∆T + mβ c∞ − cE + a0 e− p(x) + ua (x, p(x)) dx Γβ −p0 (x) + sin θβl = q , 1 + (p0 (x))2

ν

w (x) = c∞ − c

=⇒

km2 +

2 4

am φm (x) :=

∞ X

q˜ m sin(km x)e−qm p(x) km m=1 Z x = − wν (x)dx =: w(x). am

(C.1)

As Eq. C.1 is of the form of Eq. A.1, the coefficients am can be obtained in a manner identical to that given in appendix A. Note that the plane wave coefficient a0 is not determined by the above calculations; the averaged Gibbs-Thomson relations serve as the two equations to determine the two unknowns a0 and ∆T . C.2. Non-stoichiometric systems Once the phase fractions are determined implementing the algorithm of Ludwig and Leibbrandt [7], the solutal fields corresponding to the planar interface are split appropriately to give uu I and uu II with the help of the solution of the BVP UI-UII. Next, uu I is substituted into the r.h.s of the following equation which is nothing but the boundary condition of problem UI after re-arrangement of the terms. The l.h.s is treated as an unknown to be determined. −n ∇u · →

h i ds − u = − Λν u + H ν cνE fν − cE Λν fν dx Z1/2 h i ν E νE ν +  Λ c − c + H Λν u dx

(C.2)

−1/2

Integrating the above equation we have ∞ X m=1

am

q˜ m sin(km x)e−qm p(x) = −w(x) km

where w(x) is the r.h.s of Eq. C.2 integrated over (0, x), and for p(x), the reconstructed profile is used. The am s are determined by following the algorithm of appendix A. Similar calculations are performed with the UII problem. The obtained solutions are used in re-predicting c∞ , ∆T and the profile shape. In this manner, by substituting the solute distribution of the previous profile in only a part of the boundary condition, the concentration fields associated with the new profile are evaluated and the process is repeated. Alternatively, instead of working with two BVPs as above, the same method can be applied on Eq. 8. However, c∞ and a0 s have to be treated as unknowns alternatingly. While one of them is being treated as an unknown, the value estimated from the previous iteration is used for the other. The convergence of these iterative methods will be tested and reported in a future work.


am (−km ) sin(km x) p0 (x) + qm cos(km x) e−qm p(x)

m=1

am cos(km )e−qm p(x) + (c∞ − cνE )

m=1

=⇒ ∞ X

and

r

0

C.1. Stoichiometric systems The problem of finding the self-consistent profiles and the corresponding solute distributions and undercoolings when the far-field composition and speed of growth are specified can be solved accurately in case of stoichiometric phase diagrams for any p´eclet number. As the concentration in the solid phases is constant in such systems, specifying the far-field composition is equivalent to providing the phase fractions. Hence, the “twostep” route taken for non-stoichiometric systems is not necessary and the associated discrepancies between the input c∞ , and the re-predicted ones, however small, do not arise. The method of obtaining the solute distribution in the low p´eclet number regime is detailed in Ref. [13]. The adjustments for rapid solidification situations are given here. From Eq. 8, we have for stoichiometric systems:

∞ X

∞ X m=1

C. Solute distribution at large p´eclet numbers

=

 q˜ m = − 2

(B.2)

where, θαl , θβl are the acute angles made by the α−liquid and β−liquid surfaces with the horizontal, respectively, and xtr.jn. is the x-coordinate corresponding to the triple junction. The updated profile shape can be obtained from the above equations by computing p(x). The a0 , ∆T , and ua corresponding to the new profile serve as the input for the next iteration and the process is continued until convergence. It has to be noted that the calculation of profile shapes from Eqs. B.1 and B.2 is very sensitive to the simplifying assumptions that can be invoked with regard to the integrands. For instance, if the plane wave e− p(x) is replaced with its phase averages ζ0α and ζ0β in Eqs. B.1 and B.2, respectively, then for ( fα , fβ ) = (4/5, 1/5), the reconstructed profiles would be as given in Table B.1 instead of those reported in Table 5 at the same lamellar spacing. As the predicted undercoolings are strongly dependent on the profile shapes, the deviations between the values for the converged and the planar fronts also differ. In Ref. [13], the plane wave is considered only upto first order in p´eclet number; this is further replaced by its phase averages for interface shape computations. This is the reason behind the more curvedness and the lower undercoolings of the converged fronts in the current article.

∞ X

νE


am km sin(km x) p0 (x) + q˜ m cos(km x) e−qm p(x) = −wν (x)

m=1

13

Table B.1: The fifth row case of Table 5 reproduced by using an approximated plane wave component for the interface shape calculations

p(x)

Reconstructed interface profiles Initial guess 1st iteration 2nd iteration 3rd iteration 4th iteration

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

x

References [1] K. Jackson, J. Hunt, Lamellar and rod eutectic growth, AIME Met. Soc. Trans. 236 (1966) 1129–1142. [2] L. Donaghey, W. Tiller, On the diffusion of solute during the eutectoid and eutectic transformations, part I, Mater. Sci. Eng. 3 (4) (1968) 231–239. [3] R. Trivedi, P. Magnin, W. Kurz, Theory of eutectic growth under rapid solidification conditions, Acta Metall. 35 (4) (1987) 971–980. [4] P. Magnin, R. Trivedi, Eutectic growth: A modification of the Jackson and Hunt theory, Acta Metall. Mater. 39 (4) (1991) 453–467. [5] T. Himemiya, T. Umeda, Three-phase planar eutectic growth models for a ternary eutectic system, Mater. Trans., JIM 40 (7) (1999) 665–674. [6] L. Zheng, D. Larson Jr., H. Zhang, Revised form of Jackson–Hunt theory: application to directional solidification of MnBi/Bi eutectics, J. Cryst. Growth 209 (1) (2000) 110–121. [7] A. Ludwig, S. Leibbrandt, Generalised Jackson–Hunt model for eutectic solidification at low and large peclet numbers and any binary eutectic phase diagram, Mater. Sci. Eng., A 375 (2004) 540–546. [8] A. Choudhury, M. Plapp, B. Nestler, Theoretical and numerical study of lamellar eutectic three-phase growth in ternary alloys, Phys. Rev. E 83 (5) (2011) 051608. [9] K. Ankit, A. Choudhury, C. Qin, S. Schulz, M. McDaniel, B. Nestler, Theoretical and numerical study of lamellar eutectoid growth influenced by volume diffusion, Acta Mater. 61 (11) (2013) 4245–4253. [10] A. Catalina, P. Voorhees, R. Huff, A. Genau, A model for eutectic growth in multicomponent alloys, in: IOP Conference Series: Materials Science and Engineering, Vol. 84, IOP Publishing, 2015, p. 012085. [11] O. Senninger, P. W. Voorhees, Eutectic growth in two-phase multicomponent alloys, Acta Mater. 116 (2016) 308–320. [12] A. Lahiri, A. Choudhury, Revisiting Jackson-Hunt calculations: Unified theoretical analysis for generic multi-phase growth in a multi-component system, Acta Mater. 133 (2017) 316 – 332. [13] E. Nani, B. Nestler, K. Ankit, Analyzing the cooperative growth of intermetallic phases with a curved solidification front, Acta Mater. 159 (2018) 135–149. [14] R. Elliott, Eutectic solidification processing: crystalline and glassy alloys, Elsevier, 2013. [15] R. Folch, M. Plapp, Quantitative phase-field modeling of two-phase growth, Phys. Rev. E 72 (1) (2005) 011602. [16] See supplemental material for the intermediate steps of various derivations that are skipped in the main article. [17] A. Karma, A. Sarkissian, Morphological instabilities of lamellar eutectics, Metall. Mater. Trans. A 27 (3) (1996) 635–656.

14

Iteration

∆T

0th 1st 2nd 4th 6th 8th 10th

5.79709 4.83430 5.07946 5.05156 5.05131 5.05130 5.05130

∆T converged ∆T planar

0.87135

Highlights: • Steady state transport equations have non-unique solutions under practical input data • Up to 27% smaller undercoolings are found compared to classical JH theory predictions • While the far field composition of the converged and the flat profiles is almost same • The changes are not abrupt going from stoichiometric to non-stoichiometric systems