An extended stochastic diffusion model for ternary mixtures

Mechanics of Materials 56 (2013) 122–130

Contents lists available at SciVerse ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

An extended stochastic diffusion model for ternary mixtures Denis Anders a,⇑, Kerstin Weinberg a a

Chair of Solid Mechanics, University of Siegen, Paul-Bonatz-Straße 9-11, 57068 Siegen, Germany

a r t i c l e

i n f o

Article history: Received 10 June 2012 Received in revised form 31 August 2012 Available online 22 October 2012 Keywords: Ternary mixtures Microstructures Phase-field simulations Random fields

a b s t r a c t This manuscript presents a novel approach to the simulation of anisotropic diffusion phenomena in multicomponent materials. Since all physical processes are subjected to a certain degree of random inhomogeneity, the influence of random phenomena cannot be neglected in modern physical models. Therefore, an operator-scaling anisotropic random field embedded in a multiphase version of the commonly known Cahn–Hilliard phase-field model is presented here in order to describe phase evolution in an exemplary ternary mixture subjected to random phenomena. The arising set of coupled highly non-linear diffusion equations will be solved numerically by means of a B-spline based finite element method. Since the experimental determination of material parameters characterizing the anisotropy of the specimen is cumbersome and error-prone, our approach gives rise to a completely new perspective for the simulation of anisotropic diffusion systems. Ó 2012 Elsevier Ltd. All rights reserved.

1. Motivation Nowadays high-performance multicomponent materials become more and more important in miscellaneous manufacturing branches. Their range of applications reaches from lightweight design of vehicle and aircraft construction, joint materials in microelectronics, coating industry, high-temperature superalloys in gas turbines and combustion chambers, design of contact lenses, cosmetics fabrication, etc., up to applications such as drug design and development of implants in medical engineering. For this reason, high-performance functional materials constitute a crucial precondition for progress in scientific and technological realms. The ongoing course of an extreme miniaturization of functional devices requires a deeper understanding of the material’s microstructure. The morphological structure on the microscopic level profoundly determines the material behavior as a continuous macroscopic system. With the decreasing size of functional devices, the influence of ⇑ Corresponding author. Tel.: +49 271 740 2225; fax: +49 271 740 4644. E-mail addresses: [email protected] (D. Anders), weinberg@ imr.mb.uni-siegen.de (K. Weinberg). 0167-6636/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.mechmat.2012.10.002

microstructural quantities such as distribution of composition, characteristic size and shape of domains, which differ in structure, orientation and chemical composition, increases on their physical and mechanical properties. In order to accentuate their advantageous/useful features, it is extremely important to gain insight in the essential mechanisms of microstructure formation and evolution. In this context the present work offers an innovative view on the description of microstructures. The ultimate microstructural arrangement of a system depends on its configuration and on exterior driving forces. In addition to this, random physical imperfections within the material and random noises in thermodynamic fields influence in essence the microstructural evolution. Therefore, the impact of random phenomena cannot be neglected in modern physical models. To this end, we take up ideas from Anders and Weinberg (2011a) and introduce advanced stochastic processes—so-called operator-scaling anisotropic random processes—in order to adapt simulation results based on deterministic mathematical models to experimental observations. In Section 2 an extended multiphase stochastic Cahn– Hilliard model will be presented in order to address randomly forced phase evolution in multiphase mixtures.

123

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

Additionally a concise introduction into the basic ideas of operator-scaling anisotropic random fields is provided. In Section 3 the temporal and spatial discretization of the employed evolution equations for a representative ternary system are presented. The demanding continuity requirements for spatial discretization will be accomplished by a B-spline based finite element method according to our earlier works, cf. Anders and Weinberg (2011b), Anders et al. (2012) and Anders and Weinberg (2011c). To illustrate the flexibility and versatility of our approach this section is also devoted to representative two-dimensional simulations of several ternary morphologies. Please note that operator-scaling random fields represent a relatively new-found class of random processes which acquired vogue mostly in mathematical literature. However, due to its enormous potential as an essential ingredient in stochastic mathematical and physical modeling the authors expect that operator-scaling random fields will enter gradually engineering applications. 2. Multiphase stochastic diffusion model In this section we will outline the microstructural evolution equations for multicomponent mixtures affected by random fluctuations. The main physical variable to capture the microstructure of a non-uniform multicomponent mixture is the local composition. The composition of a functional material is commonly described by the mass concentration of its components, c1 ; c2 etc. In a binary mixture it simplifies to a single concentration field cðx; t Þ with c  c2 ¼ 1  c1 . Consequently, in a ternary mixture the local composition is characterized by two concentration fields c1 ðx; t Þ and c2 ðx; tÞ. The third composition variable implicitly follows from the conservation of mass. The temporal evolution of these local concentration fields may roughly be divided into two steps. Starting from an almost homogeneous solution at high temperatures (during the preparation of the material), the system is rapidly quenched and kept at a constant temperature level below a critical temperature T crit . In this context T crit denotes the lowest temperature of the stable/homogeneous one-phase regime. Therefore, the critical temperature characterizes the transition point from the one-phase to the multi-phase regime and vice versa. In the configuration below T crit the solution is instable and it commences to separate into equilibrium phases. The driving quantity for microstructural evolution is here the minimization of the Gibbs free energy of the system. In general the Gibbs free energy G might consist of a configurational free energy Wcon , interfacial/surface energy contribution Wint , elastic strain energy Welast and energy terms due to electrostatic interactions Welec :

Gðc1 ; c2 ; . . . ; cm ; T Þ ¼

Z 



Wcon þ Welast þ Welec þ Wint dX:

X

ð1Þ Hereby X denotes the spatial domain of the considered physical system. We will use Eq. (1) as the Ginzburg–Landau free energy for our phase-field calculations. In the following we write W ¼ Wcon þ Welast þ Welec þ Wint in order

to denote the total free energy density (in units of J/m3 ). Hereby the configurational free energy density Wcon determines for the most part the compositions and volume fractions of the equilibrium phases. The other energy contributions affect the shape and mutual arrangement of the separated domains. In order to accentuate the anisotropic effects stemming from random fluctuations the influence of Welast and Welec is neglected. The configurational energy density Wcon is formally composed according to

Wcon ¼

m m X X g 0i ðT; pÞci þ #ðT Þ ci ln ci þ g exc ðT; p; c1 ; . . . ; cm Þ; i¼1

i¼1

ð2Þ where g 0i ðT; pÞ are the Gibbs free energy densities of the pure ith component and the second part accounts for energy contributions from the entropy of mixing. The entropy contribution is multiplied by the temperature depended material parameter #. The first two energy contributions represent the classical Lewis-Randall ideal solution model. For the consideration of reasonable solutions, which generally show a non-ideal behavior, it is crucial to include a socalled excess energy contribution g exc . Usually the excess energy is represented in an expansion according to Redlich and Kister (1948). The multicomponent version of Redlich– Kister expansion for g exc reads

g exc ¼

n X l 1 X 1 X ðlÞ  ci cj vij ci  cj þ g þ ; 2 16i;j6m 3 16i;j6m ijk l¼0 i–j

ð3Þ

i–j–k

where g ijk ¼ ci cj ck K ijk þ   , cf. (O’Connell and Haile, 2005, ðlÞ p. 216). The parameters vij , etc. are independent of composition; but they generally depend on temperature and pressure. In practice, the ternary subscripted and higherorder terms are usually neglected. Therefore, the common form of Redlich–Kister expansion is given by

g exc ¼

n X l 1 X ðlÞ  ci cj vij ci  cj : 2 16i;j6m l¼0

ð4Þ

i–j

If we take into account only first order chemical interactions between the system’s constituents as it is performed in the scope of Flory–Huggins thermodynamics of mixing (Flory, 1942; Huggins, 1942), we have to truncate the Redlich–Kister expansion after the first term in the inner series subscripted by l. This yields the so-called Porter’s equation, which reads

g exc ¼

1 X ð0Þ ci cj vij : 2 16i;j6m

ð5Þ

i–j

Consequently, the formal representation of the configurational energy density can be stated as

Wcon ¼

m m X X 1 X ð0Þ g 0i ðT; pÞci þ #ðT Þ ci ln ci þ ci cj vij : 2 16i;j6m i¼1 i¼1

ð6Þ

i–j

For T > T crit Eq. (6) is a fully convex function. A separation into phases would be energetically disadvantageous in comparison to the homogeneous mixture. Thus, the mixture remains stable to all fluctuations. For T < T crit Eq. (6)

124

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

forms a piecewise concave potential, i.e., it has several local minima and a concave region (spinodal region) in between. The initially homogeneous mixture is here unstable and will decompose. After separation, the evolved phase islands start to grow. This coarsening process is driven by the minimization of the interfacial energy density in such a manner that the size of the emerged phase regions increases at the expense of their number. The interfacial energy

Wint ¼

m X ji i¼1

2

krci k2

ð7Þ

is related to the curvature of the boundary of the segregated phase islands. Throughout this manuscript the differential operator rðÞ characterizes the spatial gradient and kk denotes the Euclidean norm. The material parameters ji are related to surface energy density c and length ~l of the transition regions between the domains of each phase. A comprehensive overview on techniques to incorporate other types of energy densities into the phase-field framework may be found in Qin and Bhadeshia (2010). In classical phase-field formalism, microstructural evolution equations are characterized by the usual model Btype Ginzburg–Landau equations for i ¼ 1; . . . ; m

  @ci dG þ ni ; ¼ r  Mr dci @t

on X  ½0; T 

ð8Þ

where M is a positive-definite matrix containing the kinetic coefficients of the system. The differential operator r  ðÞ symbolizes the divergence and ½0; T  represents the considered time interval. ni denotes here a suitable random noise which will be specified in the next section. dG Please note that the expression dc means here the variai tional derivative of G with respect to the concentration dG variable ci . In literature on materials science the term dc i is often referred to as the i-th chemical potential of the system. In our framework it holds that

dG @ Wcon @ Wint ¼ r dci @ci @ rci

! ð9Þ

including the natural boundary condition rci  n ¼ 0, where n is the unit outward normal vector to @ X. The expressions

@ ðÞ @ci

and

@ ðÞ @ rci

characterize the partial derivative

with respect to the concentration field ci and the concentration gradient rci , respectively. A detailed derivation of Eq. (9) may be found in A. For a ternary mixture, system (8) reduces to a set of two coupled partial differential equations for c1 and c2 (and implicitly c3 ) according to







@c1 ¼ r  Mr @ c1 Wcon ðc1 ; c2 ; c3 Þ  2j1 Dc1  j2 Dc2 þ n1 ; @t    @c2 ¼ r  Mr @ c2 Wcon ðc1 ; c2 ; c3 Þ  j1 Dc1  2j2 Dc2 þ n2 ; @t ð10Þ whereby the second order differential operator DðÞ symbolizes the Laplacian, defined as the divergence of the spatial gradient. Please note that for the sake of transparency the partial derivatives @@cðiÞ are replaced by the symbol @ ci ðÞ.

Due to the conservation of mass, which can be expressed by the condition

c1 þ c2 þ c3 ¼ 1 () c3 ¼ 1  ðc1 þ c2 Þ;

ð11Þ

we need only to solve equations for c1 and c2 in order to obtain the microstructural configuration of the entire system. Without loss of generality we assume an averaged, isotropic and tensor-valued mobility matrix M ¼ MI with constant Onsager coefficient M. Also the gradient coefficients ji are all set equal to an overall coefficient j. Then the system (10) reads @c1 @t @c2 @t

   ¼ r  Mr @ c1 Wcon ðc1 ; c2 ; c3 Þ  2jDc1  jDc2 þ n1 ;    ¼ r  Mr @ c2 Wcon ðc1 ; c2 ; c3 Þ  jDc1  2jDc2 þ n2 : ð12Þ

This system of partial differential equations is usually subjected to the initial condition c1 ðx; 0Þ ¼ c1;0 ðxÞ; c2 ðx; 0Þ ¼ c2;0 ðxÞ together with the common natural and mass conserving boundary conditions

rc1  n ¼ rc2  n ¼ 0;

ð13Þ

dG dG  n ¼ Mr  n ¼ 0 on @ X  ½0; T : Mr dc1 dc2

ð14Þ

2.1. Operator-scaling random fields In order to address the shortcomings of the deterministic Cahn–Hilliard model ignoring all kinds of random phenomena (Cook, 1970; Langer, 1971) extended the original Cahn–Hilliard equation by random thermal fluctuations. Unfortunately these classical stochastic diffusion models are restricted to isotropic spatially and temporally uncorrelated random fields. However, isotropy and short-range dependence are a serious drawback for many physical applications as shown by the studies of Benson et al. (2006), Zhang et al. (2006) who incorporated operatorscaling anisotropic random fields to anisotropic transport phenomena from hydrology and anomalous diffusion. Motivated by these studies we included operator-scaling random fields as initial noises in our multicomponent diffusion model (8). Here we will give some basic ideas about operator-scaling random fields for spatially two-dimensional settings. For ðx; tÞ 2 X  ½0; tmax ; t max 6 T and parameters (Hurst indices) 0 < H1 ; H2 ; H3 < 1 these fields are defined by

nðx; tÞ ¼ Re

Z



R2 R

  eixy  1 eiðtmax tÞs  1  w1 ðyÞ

 w2 ðsÞWðdy  dsÞ;

ð15Þ

where 2 1  H1 þH 4 2 w1 ðyÞ ¼ jy  h1 j2H1 þ jy  h2 j2H2 ;

1

w2 ðsÞ ¼ jsj

1þ2H3 2

:

ð16Þ ð17Þ

The random measure Wðdy  dsÞ acts here as an isotropic complex Gaussian random measure. Vectors h1 ; h2 denote normalized eigenvectors of a spatial scaling matrix/opera1 tor E with given eigenvalues H1 1 ; H 2 , respectively. The

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

field (15) is anisotropic if the spatial scaling relation has different Hurst indices in the—not necessarily orthogonal—directions of eigenvectors h1 ; h2 . The intensity of the field in the course of time is scaled by H3 . During our simulations in the next section we will use the parameters H1 ; H2 ; H3 and the mutual position of vectors h1 ; h2 in order to initiate anisotropic microstructural evolution. With respect to mass conservation in a representative domain X we need to claim for an appropriate random field

Z

nðx; tÞ dx ¼ 0 8t 2 ½0; tmax ;

ð18Þ

X

which we easily get by the transformation

1 nðx; tÞ ¼ nðx; tÞ  jXj

Z

nðy; tÞ dy

8t 2 ½0; t max :

ð19Þ

X

This expression is well-defined since the random field  n exists for the given parameters and has continuous sample path. The term jXj denotes the corresponding Lebesgue measure of X. In order to illustrate that the stochastic simulation parameters H1 ; H2 and h1 ; h2 are not merely academic parameters we want to present some ideas to estimate these quantities from experimental observations. The essence for the qualitative estimation of the spatial scaling parameters needs to be taken from image analysis of micrographs. The basic idea is to cover the observed microstructure with matching ellipses/ellipsoids in the first step. The frequency distribution of the resulting equatorial radii ax ; ay can be associated to values for H1 and H2 in our random field. The orientation of the observed microstructure is captured by the angle between the ellipse and a fixed frame of reference. This angle may be expressed indirectly by the mutual position of the orientation vectors h1 ; h2 . Fig. 1 sketches the estimation procedure by means of two exemplary microstructural settings. In the isotropic case, the resulting equatorial radii are almost equal and the ellipse looks like a sphere. In this case

125

the scaling parameters are also almost equal and the orientation vectors are perpendicular. The anisotropic case is covered by ellipses with strongly different equatorial radii and a non-vanishing ellipse angle. In our example the ellipse angle measures a  45 . This scenario may be characterized by H1 – H2 and/or orientation vectors enclosing a small angle. Actually, this is a rather heuristic approach to the parameter estimation of our stochastic simulation parameters from experiment but it provides a vivid demonstration of the relationships between our random field and microstructural key properties. However, for a more detailed mathematical treatment concerning operator-scaling random fields we refer to Zhang et al. (2006), Bierme et al. (2007). Operator-scaling anisotropic random fields may look very academic at the first glance as they are extremely sophisticated with regard to their formal structure. In the next section we will unfold the great potential of these random fields during the simulation of ternary microstructural formation. 3. Computational analysis of microstructural evolution This section is devoted to representative numerical studies of ternary systems subjected to anisotropic random fluctuations. To this end we have to introduce a consistent numerical discretization of the ternary evolution Eqs. (12). 3.1. Spatial and temporal discretization In every finite element approach the variational/weak form of the evolution equations is the starting point for the discretization procedure. In this setting the variational formulation of the ternary evolution Eqs. (12) is obtained by multiplying each equation with its corresponding variation (dc1 and dc2 , respectively) and a subsequent integration over X. Applying Green’s first identity (integration by parts) and making use of the homogeneous boundary conditions (13) and (14) provides the weak form of the problem: Find c1 ; c2 2 V # H2 ðXÞ such that

Fig. 1. Microstructural evolution in a randomly forced ternary diffusion scenario.

126

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

Z

Z @c1 dc1 dX ¼  Mr@ c1 Wcon ðc1 ;c2 ;c3 Þrdc1 dX X @t X Z Z  M j ð2Dc1 þ Dc2 ÞDdc1 dX þ n1 dc1 dX; X X Z Z @c2 dc2 dX ¼  Mr@ c2 Wcon ðc1 ;c2 ;c3 Þrdc2 dX X @t X Z Z  M j ðDc1 þ 2Dc2 ÞDdc2 dX þ n2 dc2 dX; 8dc1 ;dc2 2 H2 ðXÞ: ð20Þ X

n X c1k Nk k¼1 n X

dch1 ¼

ch2 ¼

n X c2l Nl ; l¼1

/1i Ni

dch2 ¼

i¼1

n X /2j Nj :

ð21Þ

j¼1

Inserting the identities from (21) into the system (20) and applying Galerkin orthogonality gives the semi-discrete residuals

Rci 1 ¼

nel Z [

Xe

e¼1

@ch1 Ni dX þ M @t

Z

Z Xe

@ 2c1 Wcon rch1 rNi dX

con

@ c1 c2 W rch2 rNi dX Z Z þ 2Mj Dch1 DNi dX þ Mj Dch2 DNi dX Xe Xe Z  n1 Ni dX; þM

Xe

ð22Þ

Xe

Rcj 2 ¼

nel Z [

Xe

e¼1

þM

@ch2 Nj dX þ M @t

Z

Xe

con

@ 2c2 W

Z

ch1;nþ1  ch1;n N i dX Dt e¼1 Xe Z h i þM @ 2c1 Wcon rch1 rN i d X

Xe

þM

ð23Þ

The temporal discretization of the semi-discrete residuals (22) and (23) is straightforward. The considered time interval ½0; T  is divided into nt subintervals In ¼ ½tn ; t nþ1  according to

In :

ð24Þ

n¼0

The first order time derivative is approximated by finite differences

@c cnþ1  cn ¼ @t Dt



Xe

þ 2Mj 

Z

@ c1 c2 Wcon rch2

Z Xe

nþh

rN i d X

Dch1;nþh DNi dX þ Mj

Z Xe

n1;nþ1 Ni dX;

Dch2;nþh DNi dX ð26Þ

Xe

Rcj 2 ¼

nel Z [

ch2;nþ1  ch2;n N j dX Dt e¼1 Xe Z  þM @ c2 c1 Wcon rch1 nþh rNj dX ZXe h i þM @ 2c2 Wcon rch2 rN j d X nþh Xe Z Z þ Mj Dch1;nþh DN j dX þ 2Mj Dch2;nþh DNj dX Xe Xe Z  n2;nþ1 Nj dX; ð27Þ Xe

with ½nþh denoting a temporal evaluation at t nþh , which is linearly interpolated as ½nþh ¼ ð1  hÞ½n þ h½nþ1 . The discrete residual statements (26) and (27) constitute a system of highly nonlinear equations. In order to obtain the required temporal updates ch1;nþ1 ; ch2;nþ1 we employ an iterative solution strategy embedded in the Newton Raphson method. The explicit representation of the consistent tangential matrix K required within the linearization of the residual system (26) and (27) is formally derived in B.

In order to illustrate the essential features of the introduced diffusion model, two exemplary settings within a representative ternary system are considered. Throughout our computational studies we employ a set of dimensionless system parameters given by

rch2 rN j dX

Xe

n[ t 1

nþh

Xe

Z

3.2. Numerical results

@ c2 c1 Wcon rch1 rNj dX

Z Z þ Mj Dch1 DNj dX þ 2Mj Dch2 DNj dX Xe Xe Z  n2 Nj dX:

½0; T  ¼

nel Z [

X

According to the authors’ earlier works (Anders and Weinberg, 2011b; Anders et al., 2012; Anders and Weinberg, 2011c) we introduce a finite dimensional subspace V h H2 ðXÞ which we interpret as a span of C 1 -continuous B-spline basis functions. Let n be the dimension of V h , then we are able to represent the discrete solutions ch1 ; ch2 and the discrete test functions dch1 ; dch2 as linear combinations of basis functions N i 2 V h :

ch1 ¼

Rci 1 ¼

ð25Þ

with (equidistant) time step Dt ¼ t nþ1  t n . Performing a temporal discretization by means of Crank–Nicholson scheme yields the entirely discrete residual formulation

j ¼ 2:5  105 ; #ðT Þ ¼ 0:35; ð0Þ ð0Þ vð0Þ 12 ¼ v13 ¼ v23 ¼ 1:2: M ¼ 1;

ð28Þ

These material parameters represent a thermodynamically unstable mixture which will decompose into three equilibrium phases. As initial concentration profiles we employ

c1;0 ¼ 0:20 þ 0:01f1 ;

ð29Þ

c2;0 ¼ 0:19 þ 0:01f2 ;

ð30Þ

where f1 ; f2 denote normally distributed random noises with mean 0 and variance 1. The local concentration variables c1 ; c2 and c3 characterize the local mass concentration of substance A; B and C, respectively. Before we focus on the numerical simulation of microstructural evolution in ternary mixtures we have to specify the classical morphologies observed in ternary systems. Actually, ternary systems tend to arrange in much more

127

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

(a)

(b)

(c)

(d)

(e)

Fig. 2. Geometrically possible morphologies for the solidification of eutectic ternary mixtures: (a) three lamellar phases, (b) two lamellar phases and one fibrous phase, (c) one lamellar phase and two fibrous phases, (d) two fibrous phases embedded in a continuous matrix phase and (e) three fibrous phases in a hexagonal cobblestone arrangement.

Fig. 3. Microstructural evolution in a deterministic ternary diffusion scenario.

Fig. 4. Microstructural evolution in a randomly forced ternary diffusion scenario.

versatile types of morphologies than binary mixtures. For an overview on different types of microstructures in eutectic binary alloys the reader is referred to Anders and Weinberg (2011a). For three-phase microstructures (Ruggiero and Rutter, 1997) offer a verbal classification of five morphologies expected to form during solidification of ternary

eutectic mixtures. They chose a heuristic approach to ternary morphologies by a combination of the geometrically possible arrangements in the binary case, in which the structure can only be either lamellar or fibrous. Fig. 2 provides an overview of these basic morphologies in ternary systems.

128

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

In the first computational setting we illustrate the microstructural evolution without a random forcing. To this end the random fields n1 and n2 are set zero during the simulation. Starting from an initially homogeneous composition the mixture is quenched into a thermodynamically unstable regime favoring a phase separation into three equilibrium phases, see Fig. 3. Due to the small overall mass fractions of substance A and B we observe a formation of two fibrous phases (A-rich and B-rich phase) in a continuous C-rich matrix phase. Additionally it can be seen that the minor phases tend to arrange in so-called stacks with alternating A-rich and Brich phases. The later stages of microstructural evolution are dominated by a minimization of interfacial energy leading to an overall coarsening of the considered microstructure, see Fig. 3. A qualitatively similar microstructure was observed by Kabassis et al. (1984) during aging experiments of a unidirectionally grown eutectic Bi-In-Sn solder alloy in a plane orthogonal to the solidification direction. The observed kind of morphology may be classified as type (d) from Fig. 2. Motivated by the highly anisotropic Bi-In-Sn microstructure longitudinal to the solidification direction detected in Kabassis et al. (1984) we illustrate the performance of our stochastic diffusion model. In many eutectic mixtures the equilibrium phases arrange in elongated lamellar structures along the growth direction during solidification. By means of an anisotropic random forcing during the early stages of phase separation we are able to capture this kind of microstructural arrangement. For this reason we employ an anisotropic random field n1 with spatial scaling given by

H1 ¼ H2 ¼ 0:05;

h1 ¼ e0:23pi ;

h2 ¼ e0:24pi :

ð31Þ

This field is only operative until a dimensionless system time of 0:0025 is reached. The random field n2 impelling the evolution of substance B must be of a similar orientation as n1 . A completely different orientation of n2 is unphysical because the overall system is expected to respond to locally anisotropic random phenomena in such a way that all system components align with a similar orientation. In order to guarantee that two equilibrium phases are not driven to enrich in the same locus of the specimen we simply set n2 ¼ n1 . Our simulation results are presented in Fig. 4. Please note that as long as the random field is operative the computational effort increases by approximately 20% due to the additional evaluation in comparison with the deterministic model. Our computational studies now show a completely different microstructural arrangement. The evolving structure consists of a major C-rich phase in which elongated fibers of the two A-rich and B-rich phases are embedded. To draw a parallel to the studies of Kabassis et al. (1984) we may identify the light brown major phase with the intermetallic matrix phase BiIn2 , the neon green phase with the In-rich b-phase and the red phase with the Sn-rich c-phase, respectively. Although all material parameters used in our model are isotropic an initial anisotropic

random forcing by means of a suitable random field is sufficient to induce a long term anisotropic microstructural formation. This effect clearly indicates the great potential of our stochastic diffusion model to capture versatile morphologies in multicomponent mixtures. 4. Conclusion Summing up the obtained results, it is apparent that our diffusion model is capable of reproducing isotropic as well as anisotropic diffusion events in multicomponent diffusive systems. Actually, it is possible to employ classical diffusion models with anisotropic material parameters in order to achieve similar results. Such models include anisotropic diffusion tensors and anisotropic surface/gradient energies, see McFadden et al. (1993), Qin and Bhadeshia (2009). However, these models require an exact determination of the anisotropic material parameters for various physical settings. Unfortunately, the experimental determination of such material parameters characterizing the anisotropy of the specimen is not only expensive but also cumbersome and error-prone. In contrast to this substantial deficiency, our approach requires only a rough parameter estimation of the employed noise term taken from few micrographs. Therefore, in the authors’ opinion it is only a matter of time until smart stochastic diffusion models will gain acceptance in engineering applications. Appendix A. The system’s chemical potential as variational derivative of G In Section 2 the i-th chemical potential of a multiphase mixture was introduced as the first variational derivative of the Gibbs free energy G with respect to the concentration field ci . In this appendix we intend to elaborate on some theoretical basics of variational calculus and to provide a comprehensive derivation of the explicit representation in Eq. (9). Variational calculus usually deals with functionals and therefore it operates on infinite-dimensional spaces. In our context, the calculus of variations can be employed to find an unknown function that minimizes or maximizes a considered functional, subjected to certain boundary conditions. Due to this property, many engineering problems are formulated into a variational form. In order to introduce the concept of functional derivative we first need to introduce the idea of Gâteuax derivative which is a generalization of the concept of directional derivative in differential calculus. Definition 1 (Gâteuax derivative). Suppose X and Y are locally convex topological vector spaces (e.g., Banach spaces), U X is open and F : X # Y. F is called Gâteuax differentiable at x0 2 U, if there exists the limit



F ðx0 þ hv Þ  F ðx0 Þ d ¼ F ðx0 þ hv Þ

h!0 h dh h¼0

lim

¼ dF ðx0 ; v Þ 8v 2 X:

ðA:1Þ

The mapping dF ðx0 ; v Þ : X # Y is called the Gâteuax differential of F at x0 2 U in the direction of v 2 X. Although the Gâteuax differential defines a mapping

129

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

dF ðu; Þ : X # Y

ðA:2Þ

at each point u 2 U, which is homogeneous in the sense that

dF ðu; av Þ ¼ adF ðu; v Þ 8a 2 R;

8v 2 X;

ðA:3Þ

it needs not to be linear and continuous. Due to this deficiency the definition of the first variation of a functional mandates deeper requirements. Definition 2 (The first variation of a functional). Let us consider a sufficiently smooth functional F : X # Y. If there exists dF ðx0 ; v Þ with x0 2 U X; v 2 X and dF ðx0 ; v Þ is a continuous linear operator in v, then dF ðx0 ; v Þ is referred to as the first variation of F at x0 in the direction of v and it is written as

dF ðx0 ; v Þ ¼



dF ðx0 Þ; v ; dx0

ðA:4Þ

! dG @ Wcon @ Wint : ¼ r dci @ci @ rci

ðA:10Þ

This result completes the considerations of the current appendix. Appendix B. Derivation of the consistent tangential matrix K The core of the employed Newton Raphson method is the consistent tangential matrix K required within the linearization of the residual system (26) and (27). For the sake of transparency, it is convenient to split up K into four submatrices, which gives

2 6 K¼4

Kc1 c1

Kc1 c2

c2 c 1

c2 c 2

K

3 7 5:

ðB:1Þ

K

where h; i is a natural inner product on X. If Eq. (A.4) exists for all v 2 X 0 X,

The entries of these sub-matrices are calculated according to

dF : X0 # R dx0

Kcik1 c1 :¼

ðA:5Þ

is a linear functional. In this case dF =dx0 is called functional/variational derivative of F at x0 . Now we are able to derive an explicit representation for the system’s i-th chemical potential, which was introduced as the first variational derivative of the system’s energy functional G. According to our approach, the Gibbs free energy G is generally given by

G¼

Z

e¼1

Z 

Dt

þ Mh þ Mh



Nk Ni dX

Xe

Z

Xe

ðA:6Þ

¼

d dh

Z

¼

X

N k rN i d X

rNk rNi dX

ðB:2Þ

@Rci 1 @c2l Z nel [ ¼ Mh

Kcil1 c2 :¼

e¼1

Wðc1 ; . . . ; ci þ hv ; . . . ; cm ; rc1 ;

þ Mh

!  Z @W @W @ Wcon @ Wint vþ rv dX ¼ vþ r v dX @ci @ rci @ci @ rci X

! ! @ Wcon @ Wint ¼ v r v dX @c @ rc X * ! + @ Wcon @ Wint ¼ r ; ;v @ci @ rci 2

h Xe

Z

Xe

X

. . . rðci þ hv Þ; . . . ; rcm ÞdXjh¼0 Z 

nþh

nþh

Xe

in the absence of elastic strain energy and energy terms due to electrostatic interactions. Here it is appropriate to employ X 0 ¼ L2 ðXÞ, the Banach space of square-integrable functions on X, with the corresponding L2 -inner product. Then we formally get

dG ;v dci

h i @ 2c1 Wcon

Z

X



h i @ 3c1 Wcon rch1

 @ c1 c2 c1 Wcon rch2 nþh Nk rNi dX Xe Z þ 2M jh DNk DNi dX; þ Mh

Wcon ðc1 ; . . . ; cm Þ þ Wint ðrc1 ; . . . ; rcm Þ dX

Z

Xe

Wðc1 ; c2 ; . . . ; cm ; rc1 ; . . . ; rcm Þ dX

X

¼

¼

@Rci 1 @c1k Z nel [ 1

ðA:7Þ

þ Mh

Z

Xe

nþh

N l rN i d X

 @ c1 c2 c2 Wcon rch2 nþh Nl rN i dX  @ c1 c2 Wcon nþh rNl rNi dX

Z

þ Mjh

ðA:8Þ

  i @ c2 @ 2c1 Wcon rch1

DN l DNi dX;

ðB:3Þ

Xe

Z

Kcjk2 c1 :¼ ðA:9Þ

@Rcj 2

@c1k Z nel [ ¼ Mh e¼1

L ðXÞ

where in the third line we employed Green’s first identity together with rci  n ¼ 0 on @ X as natural boundary condition. After the second line the argument was omitted for the sake of transparency. Finally, the functional/variational derivative of G with respect to the concentration field ci reads

þ Mh

Z

Xe

Xe

þ Mh

Z

Xe

þ Mhj

 @ c2 c1 c1 Wcon rch1 nþh Nk rN j dX  @ c2 c1 Wcon nþh rNk rNj dX h   i @ c1 @ 2c2 Wcon rch2

Z

Xe

DN k DNj dX;

nþh

N k rNj dX ðB:4Þ

130

D. Anders, K. Weinberg / Mechanics of Materials 56 (2013) 122–130

Kcjl2 c2 :¼ ¼

@Rcj 2 @c2l Z nel [ 1 e¼1

Dt

þ Mh

Z

Nl Nj dX

Xe



Xe

þ Mh

Z

Xe

þ Mh

Z

Xe

þ 2M jh

@ c2 c1 c2 Wcon rch1

h h

@ 3c2 Wcon rch2 @ 2c2 Wcon

Z

i nþh

nþh

i nþh

Nl rNj dX

N l rN j d X

rNl rNj dX

DNl DNj dX:

ðB:5Þ

Xe

Now all numerical instruments to handle the ternary Cahn–Hilliard system are given. References Anders, D., Weinberg, K., 2011a. Application of operator-scaling anisotropic random fields to binary mixtures. Philos. Mag. 91, 3766–3792. Anders, D., Weinberg, K., 2011b. Numerical simulation of diffusion induced phase separation and coarsening in binary alloys. Int. J. Comput. Mater. Sci. 50, 1359–1364. Anders, D., Weinberg, K., 2011c. A variational approach to the decomposition of unstable viscous fluids and its consistent numerical approximation. ZAMM J. Appl. Math. Mech. 91, 609–629. Anders, D., Weinberg, K., Reichardt, R., 2012. Isogeometric analysis of thermal diffusion in binary blends. Int. J. Comput. Mater. Sci. 52, 182– 188.

Benson, D., Meerschaert, M., Bäumer, B., Scheffler, H., 2006. Aquifer operator-scaling and the effect on solute mixing and dispersion. Water Resour. Res. 42, 1–18. Bierme, H., Meerschaert, M., Scheffler, H., 2007. Operator scaling stable random fields. Stoch. Proc. Appl. 117, 312–332. Cook, H., 1970. Brownian motion in spinodal decomposition. Acta Metall. 18, 297–306. Flory, P., 1942. Thermodynamics of high polymer solutions. J. Chem. Phys. 10, 51–61. Huggins, M., 1942. Theory of solutions of high polymers. J. Am. Chem. Soc. 64, 1712–1719. Kabassis, H., Rutter, J., Winegard, W., 1984. Microstructure of one of the ternary eutectic alloys in the bi-in-sn system. Metall. Trans. A 15A, 1515–1517. Langer, J., 1971. Theory of spinodal decomposition in alloys. Ann. Phys. 65, 53–86. McFadden, G.B., Wheeler, A.A., Braun, R.J., Coriell, S.R., Sekerka, R.F., 1993. Phase-field models for anisotropic interfaces. Phys. Rev. E 48, 2016– 2024. O’Connell, J., Haile, J., 2005. Thermodynamics: Fundamentals for Applications. Cambridge University Press, Cambridge. Qin, R., Bhadeshia, H., 2009. Phase-field model study of the effect of interface anisotropy on the crystal morphological evolution of cubic metals. Acta Mater. 57, 2210–2216. Qin, R.S., Bhadeshia, H.K., 2010. Phase field method. Mater. Sci. Technol. 26, 803–811. Redlich, O., Kister, A., 1948. Thermodynamics of nonelectrolyte solutions. Algebraic representation of thermodynamic properties and the classification of solutions. Ind. Eng. Chem. 40, 345–348. Ruggiero, M., Rutter, J., 1997. Origin of microstructure in the 332 K eutectic of the bi-in-sn system. Mater. Sci. Technol. 13, 5–11. Zhang, Y., Benson, D., Meerschaert, M., LaBolle, E., Scheffler, H., 2006. Random walk approximation of fractional-order multiscaling anomalous diffusion. Phys. Rev. E 74, 026706.