Extended Larché--Cahn framework for reactive Cahn--Hilliard multicomponent systems
EEXTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARDMULTICOMPONENT SYSTEMS
SANTIAGO P. CLAVIJO , LUIS ESPATH & VICTOR M. CALO Abstract.
At high temperature and pressure, solid diffusion and chemical reactions between rock mineralslead to phase transformations. Chemical transport during uphill diffusion causes phase separation, that is,spinodal decomposition. Thus, to describe the coarsening kinetics of the exsolution microstructure, we derivea thermodynamically consistent continuum theory for the multicomponent Cahn–Hilliard equations whileaccounting for multiple chemical reactions and neglecting deformations. Our approach considers multiplebalances of microforces augmented by multiple constituent content balance equations within an extendedLarch´e–Cahn framework. As for the Larch´e–Cahn framework, we incorporate into the theory the Larch´e–Cahn derivatives with respect to the phase fields and their gradients. We also explain the implications of theresulting constrained gradients of the phase fields in the form of the gradient energy coefficients. Moreover, wederive a configurational balance that includes all the associated configurational fields in agreement with theLarch´e–Cahn framework. We study phase separation in a three-component system whose microstructuralevolution depends upon the reaction-diffusion interactions and to analyze the underlying configurationalfields. This simulation portrays the interleaving between the reaction and diffusion processes and how theconfigurational tractions drive the motion of interfaces.
AMS subject classifications: · · · · · · Contents
1. Introduction 22. Theoretical framework 22.1. Constituent content balances 22.2. Concentrations, phase fields, and microforce balances 32.3. Thermocompatible constitutive relations 32.4. Thermodynamical constraints 32.5. Chemical reaction 53. Configurational fields 54. Dimensionless multicomponent Cahn–Hilliard equations 65. Numerical simulation: merging of circular inclusions 76. Final remarks 147. Acknowledgments 14Appendix A. Thermodynamically consistent continuum theory for the multicomponent Cahn–Hilliard equations 14A.1. Larch´e–Cahn derivatives 14A.2. Thermodynamics 15A.3. Theory of reacting materials 16A.4. Configurational stress and force 17References 18 Ali I. Al-Naimi Petroleum Engineering Research Center, King Abdullah University of Science & Technology(KAUST), Thuwal, 23955-6900, Saudi Arabia. Department of Mathematics, RWTH Aachen University, Pontdriesch 14-16, 52062 Aachen, Germany. School of Earth and Planetary Sciences, Curtin University, Kent Street, Bentley, Perth, WA 6102, Aus-tralia.Curtin Institute for Computation, Curtin University, Kent Street, Bentley, Perth, WA 6102, Australia.Mineral Resources, Commonwealth Scientific and Industrial Research Organisation (CSIRO), 10 Kensington,Perth, WA 6152, Australia.
E-mail address : [email protected] , [email protected], [email protected] . Date : January 27, 2021. a r X i v : . [ phy s i c s . g e o - ph ] J a n EXTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS Introduction
Deep in the Earth, both high temperature and pressure allow for solid diffusion and chemical reactionsbetween rock minerals, which in turn, lead to phase transformations and induced deformation. Significantly,the transport of chemical constituents during uphill diffusion generates phase separation processes as a resultof the geothermal gradient in the crust. The phases that compose these solid solutions of minerals diffuseat different rates, and when considering changes in temperature as a result of uphill diffusion, for instanceduring cooling, phase separation processes such as spinodal decomposition occur. For example, ternaryfeldspars formed by orthoclase, anorthite, and albite show spinodal decomposition during cooling. Thus,such a process controls the coarsening kinetics of the exsolution microstructure [1–3]. Rocks are complexsystems composed of several minerals, grain-boundaries, fractures, and pore space where the chemical andmechanical properties may vary in each direction. We describe each mineral as a component of a solidsolution; this interpretation of the mixture allows us to explain the coupled reactive spinodal decompositionduring exsolution. As a case study, we model a phase merging process driven by interfacial responses coupledwith chemical reactions as a first attempt to understand the dynamics of reactive exsolution by spinodaldecomposition [4]. We derive a thermodynamically-consistent reactive n species Cahn-Hilliard model thatcaptures the dynamics of such interactions while following in detail the configurational forces that drive thiscoupled kinematical process.The multicomponent Cahn–Hilliard model is a useful tool for studying the kinetics of multiphase systemsundergoing phase separation. Most importantly, this model tracks the microstructure evolution of theresulting phases to enhance our understanding of the resulting material properties. To describe the underlyingphysics of this problem, we consider n phase fields representing the concentration of conserved species anduse a set of coupled Cahn–Hilliard equations. This representation leads to a system of n degenerate nonlinearfourth-order parabolic partial differential equations. The degeneracy is due to a nonlinear mobility tensorthat can vanish depending on the phase field values. We assume that there exist n microforce balances, assimilarly proposed by Fried & Gurtin [5, 6], and n mass balances accounting for all the relevant chemicalreactions, as similarly proposed by Clavijo et al. [7]. We then build an extended Larch´e–Cahn framework toaccount for the interdependence between the conserved species. Given the set ϕ = { ϕ , . . . , ϕ n } of species,where n ∈ N , we consider n − ϕ = ϕ \ { ϕ σ } while the σ -th conserved species isused as a reference and determined by ˜ ϕ . Thus, to compute partial derivatives with respect to ϕ α andgrad ϕ α , the dependence among species must be taken into account. We then redefine the partial derivativeof functions depending upon ϕ and grad ϕ where ϕ σ is constrained to derive the multicomponent Cahn–Hilliard equations. Moreover, in defining these partial derivatives, we arrive at a constrained inner producton a constrained space to appropriately define the gradient energy coefficients Γ αβ .The outline of this article is as follows. In section §
2, we introduce the balances of microforces andaugment them with mass balances. In section §
3, we present the configurational forces and their balances,and describe how they drive the interface evolution. In section §
4, we make the equations dimensionless.Section § Theoretical framework
We give a brief overview of the theoretical framework that describes the isothermal evolution of n ≥ B of a three-dimensional point space.2.1. Constituent content balances.
We assume that a mass density (cid:37) α , a diffusive flux α , and a reactivemass supply rate s α characterize the instantaneous state of each constituent α = 1 , . . . , n . Also, we requirethat (cid:37) α , α , and s α evolve subject to a pointwise constituent content balance in the form˙ (cid:37) α = − div α + s α , ∀ ≤ α ≤ n, (1)where a superposed dot denotes partial differentiation with respect to time and div denotes the divergenceon B . Stipulating that the mass supply rates and the diffusive fluxes satisfy constraints of the form n (cid:88) α =1 s α = 0 and n (cid:88) α =1 α = 0 , (2) XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 3 we sum the constituent content balance (1) over α from 1 to n to find that the total mass density (cid:37) = n (cid:88) α =1 (cid:37) α (3)must satisfy ˙ (cid:37) = 0 and, thus, is constant.2.2. Concentrations, phase fields, and microforce balances.
Introducing a concentration ϕ α = (cid:37) α (cid:37) (4)for each species α = 1 , . . . , n , from expressions (3) and (4) together with the requirement that the total massdensity (cid:37) is a fixed constant, we have that the following constraint n (cid:88) α =1 ϕ α = 1 (5)must hold in conjunction with (2). Moreover, from constituent content balance (1) for constituent α =1 , . . . , n , we have that (cid:37) ˙ ϕ α = − div α + s α , ∀ ≤ α ≤ n. (6)Next, let ξ α be the α -th microstress, and π α ( γ α ) field is the α -th internal (external) microforce. Thus,we express the microforce balances of Fried & Gurtin [8, § IV] in its pointwise form asdiv ξ α + π α + γ α = 0 . (7)In the partwise form of expression (7), the surface microtraction is ξ α := ξ α · n .2.3. Thermocompatible constitutive relations.
We introduce constitutive relations for the diffusiveflux α , the reactive mass supply rate s α , the internal microforce density π α , and the microstress ξ α foreach constituent α = 1 , . . . , n , which allow us to close the system of evolution equations for the phase fields ϕ α , α = 1 , . . . , n . These relations must be compatible with the constraints (2) and (5) and with the firstand second laws of thermodynamics, which, since we consider only isothermal processes, combine to yieldan inequality of the form (cid:37) ˙ ψ − n (cid:88) α =1 { ( (cid:37)µ α − π α ) ˙ ϕ α − α · grad µ α + µ α s α + ξ α · grad ˙ ϕ α } ≤ , (8)where ψ is the specific free-energy, and µ α is the chemical potential of constituent α (for details, see Appen-dix A.2). We define the chemical potential using the Coleman–Noll procedure in the next section.2.4. Thermodynamical constraints.
Throughout the derivation of the constitutive relations for the mul-ticomponent Cahn–Hilliard system, we use the Larch´e–Cahn derivative (49) from Appendix A.1. Using theColeman–Noll procedure [9], we find the sufficient conditions to ensure the inequality (8) for arbitrary fields.Thus, a set of paired constitutive equations emerges for each kinematic process. We assume the followingconstitutive dependency of the free energy ψ within the context of isothermal processes ψ := ˆ ψ ( ϕ , grad ϕ ) , (9)which specializes the free-energy (8) as follows n (cid:88) α =1 (cid:40) (cid:37) µ α − π α − (cid:37) ∂ ( σ ) ˆ ψ∂ϕ α (cid:41) ˙ ϕ α + n (cid:88) α =1 (cid:40) ξ α − (cid:37) ∂ ( σ ) ˆ ψ∂ (grad ϕ α ) (cid:41) · grad ˙ ϕ α − n (cid:88) α =1 { α · grad µ α + µ α s α } ≤ . (10) EXTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS
The free-energy imbalance (10) must hold for any arbitrary ˙ ϕ α , grad ˙ ϕ α , and grad µ α fields at a given timeand place. Thus, the following relations must hold π ασ = (cid:37) (cid:32) µ ασ − ∂ ( σ ) ˆ ψ∂ϕ α (cid:33) , (11a) ξ ασ = (cid:37) ∂ ( σ ) ˆ ψ∂ (grad ϕ α ) , (11b) ασ = − n (cid:88) β =1 M αβ grad µ βσ , (11c)where M is the mobility tensor, which must be positive semi-definite, that is, (cid:80) nα =1 (cid:80) nβ =1 p α · M αβ p β ≥ p . As (10) expresses all thermodynamically consistent choices, we write the terms π α := π ασ , µ α := µ ασ , and ξ α := ξ ασ relative to the Larch´e–Cahn construction given their explicit dependence on theLarch´e–Cahn derivatives as an essential consequence of (5). As a byproduct, we also write the mass flux, α := ασ ( x , t ; grad µ ασ ), and the surface microtraction, ξ α S := ξ α S σ ( x , t ; ξ ασ ) as constructions dependent on theLarch´e–Cahn derivative. Finally, intrinsically in these definitions, we express all quantities relative to the σ -th species.Guided by the original Cahn–Hilliard equation [10], we assume that the evolution of the Ginzburg–Landaufree energy governs the nature of phase separation undergoing spinodal decomposition. In a multicomponentsystem, we determine the constitutive relations in (11) from the Ginzburg-Landau free energy expressed asˆ ψ ( ϕ , grad ϕ ) = N v k B ϑ (cid:32) n (cid:88) α =1 ϕ α ln ϕ α (cid:33) + N v n (cid:88) α =1 n (cid:88) β =1 Ω αβ ϕ α ϕ β + 12 n (cid:88) α =1 n (cid:88) β =1 Γ αβ grad ϕ α · grad ϕ β , (12)where N v is the total number of molecules of the species α per unit volume, k B is the Boltzmann constant,and Ω αβ represents the interaction energy between the mass fraction of the α -th and β -th species, whichis reciprocal; thus, Ω αβ is symmetric. The interaction energy is positive and is related to the critical tem-perature for each pair of species, ϑ αβc , (between the α -th and β -th species). Following standard convention,we adopt that Ω αβ = 0 when α = β and Ω αβ = 2 k B ϑ αβc when α (cid:54) = β [10–12]. Furthermore, Γ αβ = σ αβ (cid:96) αβ [force] (no sum on α and β ) represents the magnitude of the interfacial energy between the α -th and β -thspecies. The parameters σ αβ and (cid:96) αβ are the interfacial tension [force/length] and the interfacial thickness for each pair of species (between the α -th and β -th species) [length], respectively. Cahn & Hilliard [10] definethe force Γ αβ as N v Ω αβ ( (cid:96) αβ ) .We express the relative chemical potential of the α -th species, in the Larch´e–Cahn sense, by combiningthe expressions (11a), (11b), the microforce balance (7), and the constitutive relation for the free energy (12),we arrive at µ ασ = ∂ ( σ ) ˆ ψ∂ϕ α − div ∂ ( σ ) ˆ ψ∂ (grad ϕ α ) − (cid:37) ( γ α − γ σ ) . (13)Therefore, the combination of (13) with (12) specializes to µ ασ = N v k B ϑ ln (cid:18) ϕ α ϕ σ (cid:19) + 2 N v n (cid:88) β =1 (Ω αβ − Ω σβ ) ϕ β − n (cid:88) β =1 (Γ αβ − Γ σβ ) div grad ϕ β − (cid:37) ( γ α − γ σ ) . (14)In the following, we assume an isotropic mobility M αβ = M αβ being M αβ symmetric, but we considerthe off-diagonal terms in the Onsager reciprocal relations. We use the standard assumption that the mobilitycoefficients depend on the phase composition. In particular, we express this dependency in terms of theconcentration of each species. We use the definition M αβ := M αβ ϕ α ( δ αβ − ϕ β ) with no summation on α and β and M αβ is the mobility between the α and β species, with dimension of length per unit force andtime [11]. Thus, (2) implies the following relation n (cid:88) β =1 M αβ = 0 , ∀ α. (15) This expression corresponds to the root mean square effective ”interaction distance”, as suggested by Cahn & Hilliard [10].
XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 5
Chemical reaction.
Let ϕ α be the concentration of a species A α , such that ϕ α := [ A α ]. FollowingKrambeck [13], we express the c -th chemical reaction, in a set of n s chemical reactions, n s ∈ N , as n (cid:88) α =1 υ c,αr A α k c + (cid:10) k c − n (cid:88) α =1 υ c,αp A α , ∀ ≤ c ≤ n s , (16)where υ c,αr and υ c,αp are the α -th stoichiometric coefficient of the c -th chemical reaction of the reactants andproducts, respectively. The number of non-zero stoichiometric coefficients υ c,αr and υ c,αp define the numberof reactants n cr and products n cp in the c -th chemical reaction. For the c -th chemical reaction, zeros populate υ c,αr for α > n cr , whereas zeros populate υ c,αp for α ≤ n cr . k c + ( k c − ) denotes the c -th forward (backward)reaction rate (see, for details Appendix A.3). We focus on ideal materials, then, the c -th rates of both theforward and backward reactions read r c + := k c + n (cid:89) α =1 ( ϕ α ) υ c,αr , (17) r c − := k c − n (cid:89) a =1 ( ϕ α ) υ c,αp . (18)Finally, the internal rate of mass supply term for all n s chemical reactions that enters in (6) is s α := − n s (cid:88) c =1 ( υ c,αr − υ c,αp )( r c + − r c − ) . (19)3. Configurational fields
We describe the interfacial evolution, and its thermodynamics using the configurational forces proposedby Gurtin [14], which relate the integrity of the material and the movement of its defects. The configu-rational forces expend the power associated with the transfer of matter, which allow us to interpret themthermodynamically. Using the configurational balance for a part P by Fried [15], we have (cid:90) S Cn d a + (cid:90) P ( f + e ) d v = , (20)which renders, after localization, div C + f + e = , (21)where C is the configurational stress tensor and f ( e ) is the internal (external) force.Following Appendix A.4, we substitute the constitutive relation (78) in the relation (74), allows us toexpress the configurational stress as C := (cid:37) (cid:32) ψ − n (cid:88) α =1 µ α ϕ α (cid:33) − n (cid:88) α =1 grad ϕ α ⊗ ξ α . (22)We obtain explicit forms for the internal and external configurational forces by combining (11) and (21)with (22), that is f := n (cid:88) α =1 (cid:37) ϕ α grad µ α and e := − n (cid:88) α =1 γ α grad ϕ α . (23)By considering the Larch´e-Cahn derivatives, we express the configurational stress (22) as a configurationalstress relative to the σ -th species as follows C σ := (cid:37) (cid:32) ψ − n (cid:88) α =1 µ ασ ϕ α (cid:33) − n (cid:88) α =1 grad ϕ α ⊗ ξ ασ , (24)while f σ := n (cid:88) α =1 (cid:37) ϕ α grad µ ασ , (25)is the relative internal configurational force. The external configurational force is not determined using aconstitutive relation; thus, it does not depend upon the choice of the reference species. EXTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS
Remark 1 ( Invariance of configurational balance to reference species).
Let ϕ σ be the referencespecies. We establish the following relations for the terms appearing in the configurational stress (24) − n (cid:88) α =1 µ ασ ϕ α = − n (cid:88) α =1 µ α ϕ α + µ σ n (cid:88) α =1 ϕ α = − (cid:32) n (cid:88) α =1 µ α ϕ α (cid:33) + µ σ , (26) and n (cid:88) α =1 grad ϕ α ⊗ ξ ασ = n (cid:88) α =1 grad ϕ α ⊗ ( ξ α − ξ σ ) = n (cid:88) α =1 grad ϕ α ⊗ ξ α − (cid:32) n (cid:88) α =1 grad ϕ α (cid:33) ⊗ ξ σ = n (cid:88) α =1 grad ϕ α ⊗ ξ α . (27) while for the internal configurational force (25) n (cid:88) α =1 ϕ α grad µ ασ = n (cid:88) α =1 ϕ α grad µ α − grad µ σ n (cid:88) α =1 ϕ α = (cid:32) n (cid:88) α =1 ϕ α grad µ α (cid:33) − grad µ σ . (28) Analyzing (26) and (27) , we conclude that only one term in C σ (24) depends on the reference species.Therefore, the relative configurational stress becomes C σ := C + (cid:37) µ σ . (29) Meanwhile, we specialize representation of the relative internal configurational force (25) with (28) yielding f σ := f − (cid:37) grad µ σ , (30) Finally, although both the configurational stress and the internal configurational force explicitly dependon the choice ϕ σ ; nevertheless, their dependencies cancel each other’s contribution to the configurationalbalance (21) , div C σ + f σ = div ( C + (cid:37) µ σ ) + f − (cid:37) grad µ σ , = div C + f . (31) (cid:3) Dimensionless multicomponent Cahn–Hilliard equations
The final system resulting from (6), (11), (14), and (17)-(19) reads ˙ ϕ α = s α − div ασ , ασ = − n (cid:88) β =1 M αβ ϕ α ( δ αβ − ϕ β ) grad µ βσ ,µ βσ = N v k B ϑ ln ϕ β ϕ σ + 2 N v n (cid:88) α =1 (Ω βα − Ω σα ) ϕ α − n (cid:88) α =1 (Γ βα − Γ σα ) div grad ϕ α − ( γ β + γ σ ) ,s α = − n s (cid:88) c =1 (cid:40) ( υ c,α − (cid:36) c,α )( k c + n (cid:89) a =1 ( ϕ a ) υ ca − k c − n (cid:89) a =1 ( ϕ a ) (cid:36) ca ) (cid:41) , (32)in D × (0 , T ) with (cid:40) ϕ α ( x ,
0) = ϕ α , in D , subject to periodic boundary conditions on ∂ D × (0 , T ) , (33)where D is the domain of interest.To make the equations dimensionless, we introduce the reference energy density ψ := 2 N v k B ϑ and definethe set of diffusion coefficients D αβ , D αβ = ψ M αβ ϕ α ( δ αβ − ϕ β ) no sum on α and β. (34) XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 7
The reference energy density relates the species mobilities with the species diffusion as proposed in [16, 17].We also define the following dimensionless variables x = L − x , t = T − t, ϑ αβc = ϑ − ϑ αβc . (35)Conventionally, the definition of the reference time T for the Cahn–Hilliard system relates the diffusioncoefficient, the interfacial thickness, and domain length, that is, T = D (cid:96) L − where L (cid:29) (cid:96) [18, 19].We set D and (cid:96) as the reference diffusion coefficient and interface thickness of a reference species, andintroduce the following dimensionless numbers for the multicomponent system ¯ k + c = k + c D − (cid:96) − L , ¯ k − c = k − c D − (cid:96) − L , ψ = ψ − ψ, ¯ σ αβ = σ αβ ( ψ L ) − , ¯ D αβ = D αβ D − (cid:96) − L , ¯ (cid:96) αβ = L − (cid:96) αβ , ¯ γ α = ψ − γ α . (36)Thus, by inserting the dimensionless quantities in (32), we find the following dimensionless forms ˙ ϕ α = s α − div¯ ασ , ¯ ασ = − n (cid:88) β =1 ¯ D αβ grad ¯ µ βσ , ¯ µ βσ = 12 ln ϕ β ϕ σ + 2 n (cid:88) α =1 ( ¯ ϑ βαc − ¯ ϑ σαc ) ϕ α − n (cid:88) α =1 (¯ σ βα ¯ (cid:96) βα − ¯ σ σα ¯ (cid:96) σα ) div grad ϕ α − (¯ γ β − ¯ γ σ ) , ¯ s α int = − n s (cid:88) c =1 (cid:40) ( υ c,α − (cid:36) c,α )(¯ k c + n (cid:89) a =1 ( ϕ a ) υ ca − ¯ k c − n (cid:89) a =1 ( ϕ a ) (cid:36) ca ) (cid:41) , (37)in D × (0 , T ) , with the initial condition (33).5. Numerical simulation: merging of circular inclusions
We now simulate the interactions between three species where A and A represent the reactants, while A the reaction products. The inclusions (represented by species 1, A ) are embedded in species 2, A . Weexpress the chemical reaction as A + A k + (cid:42) A (38)which takes place at the interface producing the third species, A . Table 1.
Chemical and physical parametersPhysical parameter Value Name ψ [J m − ] 2 × Energy density L [m] 10 − Domain length ϑ [K] 1000.0 Absolute temperature ϑ c [K] 1100.0 Critical temperature between phases 1 and 2 ϑ c [K] 1200.0 Critical temperature between phases 1 and 3 ϑ c [K] 1300.0 Critical temperature between phases 2 and 3 D [m s − ] 10 − Diffusion coefficient (for all phases) k + [m s − ] 10 − Forward reaction rate σ [J m − ] 0.816 Interfacial energy between phases 1 and 2 σ [J m − ] 0.625 Interfacial energy between phases 1 and 3 σ [J m − ] 0.921 Interfacial energy between phases 2 and 3 (cid:96) [m] 1 . × − Interface thickness between phases 1 and 2 (cid:96) [m] 2 × − Interface thickness between phases 1 and 3 (cid:96) [m] 10 − Interface thickness between phases 2 and 3We state the problem as: find ϕ satisfying (37) given (33) subject to periodic boundary conditions up tothe fourth derivative of ϕ with respect to x in a square open region D = (0 , × (0 , EXTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS resulting system of partial differential equations using PetIGA [20], a high-performance isogeometric analysisframework. We solve this system of equations in their primal form using a 256 ×
256 element mesh of apolynomial degree 4 and continuity 3. The initial and boundary conditions are h = 0 . δ ( x ,
0) = − . (cid:18) . (cid:18) ( x − . + ( y − . h ( h + 0 . (cid:19) + 0 . (cid:19) + 0 . δ ( x ,
0) = − . (cid:18) . (cid:18) ( x − . + ( y − . h ( h + 0 . (cid:19) + 0 . (cid:19) + 0 . ϕ ( x ,
0) = 1 + δ + δ , in D ,ϕ ( x ,
0) = 0 . − ϕ , in D ,ϕ ( x ,
0) = 1 − ϕ − ϕ , in D , subject to periodic boundary conditions on ∂ D × (0 , T ) , (39)and the three subfigures on top of Figure 1 depict this initial condition.Table 1 summarizes the dimensional parameters used to obtain the dimensionless parameters in (40)and (41). The diffusion matrix for each entry α and β reads¯ D αβ = 1 × ϕ α ( δ αβ − ϕ β ) ∀ ≤ α, β ≤ n. (40)Next, for clarity, we represent α and β as matrix-columns and -rows indices, which render the remainingdimensionless parameters as follows.¯ σ αβ ¯ (cid:96) αβ = − − .
121 6 . .
121 0 4 . .
250 4 .
605 0 , υ αβ = (cid:2) (cid:3) , (cid:36) αβ = (cid:2) (cid:3) , ¯ k + = 0 . , (41)where we choose D = D and (cid:96) = (cid:96) as the reference diffusion coefficient and interface thickness of areference species, respectively.Here, the configurational tractions drive the interfacial motion in this multicomponent system undergoingreactions. We express the configurational traction along a level curve L α ∗ , upon which ϕ α = ϕ α ∗ . We thenintroduce the normal and tangential coordinates n α and m α on L α ∗ , with unit vectors ν α and τ α definedsuch that grad ϕ α = | grad ϕ α | ν α , | ν α | = 1 , (42)augmented by a sign convention which ensures that rotating τ α clockwise by π/ ν α . In reckoningthe relative configurational stress in a { n α , m α } -frame, we arrive at C σ = ζ − n (cid:88) α =1 | grad ϕ α | ν α ⊗ ξ ασ , (43)with ζ := (cid:37) ( ψ − (cid:80) nα =1 µ α ϕ α ), see (78) in Appendix A.4. We can now specialize (43) with a free-energy ofthe form ˆ ψ ( ϕ , grad ϕ ) = f ( ϕ α ) + 12 n (cid:88) α =1 n (cid:88) β =1 Γ αβ grad ϕ α · grad ϕ β , = f ( ϕ α ) + 12 n (cid:88) α =1 n (cid:88) β =1 Γ αβ | grad ϕ α || grad ϕ β | ν α · ν β , (44)which renders the following relative configurational stress C σ = ζ + n (cid:88) α =1 | grad ϕ α | ν α ⊗ n (cid:88) β =1 (Γ αβ − Γ σβ ) | grad ϕ β | ν β . (45) XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 9
Thus, the configurational tractions C σ ν α are C σ ν α = ζ + n (cid:88) ˆ α =1 n (cid:88) ˆ β =1 (cid:16) (Γ ˆ α ˆ β − Γ σ ˆ β ) | grad ϕ ˆ α || grad ϕ ˆ β | ν ˆ α ⊗ ν ˆ β (cid:17) ν α = ζ ν α + n (cid:88) ˆ α =1 (cid:0) (Γ ˆ αα − Γ σα ) | grad ϕ ˆ α || grad ϕ α | ν ˆ α (cid:1) . (46)In the simulations, we compute the relative physical and chemical quantities, such as the relative chemicalpotential, mass fluxes, microstresses, and byproducts, by setting the reaction product species, that is, A asthe reference phase field. This simulation shows that the configurational fields can describe the behavior ofthe phase evolution. However, this initial work does not exploit this tool exhaustively nor comprehensively.Figure 1 depicts the merging process of two circular inclusions of distinct size into a single one. The figurespans from the early stages until the merged inclusion becomes stationary. From left to right, we depictphases ϕ , ϕ , and ϕ , while from top to bottom, the evolution of the three phases for the dimensionlesstimes t = 0, 6 . × − , 2 . × − , 7 . × − , and 5 . × − .Figures 2 and 3, respectively, present e · C ν ( x = 0 . , x ) and e · C ν ( x = 0 . , x ) on the leftpanel, and e · C ν ( D ) and e · C ν ( D ) on the right panel. That is, the left panels display the profile ofthe relative configurational traction along x , while the right panels display the vertical component of therelative configurational traction on the whole domain. These figures show the x axis in red. From top tobottom, we present these configurational fields at the dimensionless times t = 0, 5 . × − , 6 . × − ,6 . × − , 2 . × − , 7 . × − , and 5 . × − . Figure 2 shows that configurational tractionsbetween the inclusions have opposite directions pushing against one another. As the inclusions approacheach other, the configurational traction profiles become antisymmetric in the region where the merging takesplace (second plot, from top to bottom, in Figure 2). In this region, the ridge and the valley propagatetowards each other until the interfaces merge. Later, the configurational tractions annihilate one another(third plot, from top to bottom, in Figure 2). The third species appears as the chemical reaction takes place.Figure 3 shows how the relative configurational traction C ν pushes apart the boundaries of the doublering, formed by this species. This traction drives the growth of the area encircled by the double ring, whichoccurs at the expense of the other two species through the chemical reaction. Figure 3 (second plot, fromtop to bottom) shows the tractions on each ring as they push against each other, which favours merging. Atlater stages, a single ring-like structure remains, formed by the product species. This ring lies in between theinterface formed by the reactant species. Consequently, the process reaches a semblance of a steady-statewhen the product species obstructs further the chemical reactions.Figure 4 presents a snapshot sequence detailing the merging process from left to right and top to bottom.We use the relative configurational traction C ν , their streamlines (white), L . (black) on top of the phasefield ϕ at the dimensionless times t = 0, 5 . × − , 6 . × − , 6 . × − , and 2 . × − to exemplifythis evolution. In these snapshots, we show the configurational traction C ν with black arrows. Beforethe merging occurs, two node sinks arise, see Figure 4b. These node sinks pull the phase ϕ initiating themerging process. Soon after the node sinks are formed, Figure 4c, phase ϕ migrates and leaves a ‘bridge’between the inclusions. This ‘bridge’ is formed by phase ϕ . After, the merging of the level curve L . (blackline) occurs, see Figure 4d. Figure 1.
Phase-field evolution during merging. From left to right ϕ , ϕ , and ϕ . Fromtop to bottom t = 0, t = 6 . × − , t = 2 . × − , t = 7 . × − , t = 5 . × − . XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 11
Figure 2.
Vertical component of the configurational traction C ν along x = 0 . Figure 3.
Vertical component of the configurational traction C ν along x = 0 . XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 13
Figure 4.
Node sinks triggering the merging C ν . From left to right, the relative con-figurational traction C ν , streamlines of this traction (white), L . (black) on top of thephase field ϕ at the dimensionless times t = 0, 5 . × − , 6 . × − , 6 . × − , and2 . × − . Final remarks
In this work, we present a continuum framework to model phase separation processes such as spinodaldecomposition during cooling as a result of uphill diffusion. These phases are composed of solid solutions ofminerals and diffuse at different rates. In this first attempt to model solid diffusion and chemical reactionsbetween rock minerals, we neglect deformation and heat transfer. To this end, we derive a thermodynamicallyconsistent continuum theory for the multicomponent Cahn–Hilliard equations while accounting for multiplechemical reactions. We consider multiple balances of microforces augmented by multiple constituent contentbalance equations within an extended Larch´e–Cahn framework. Moreover, we derive a configurational bal-ance that includes all the associated configurational fields in agreement with the Larch´e–Cahn framework.In a simple simulation, we depict the role of the configurational tractions during the merging process coupledwith a chemical reaction. Last, in upcoming works, we plan to model the contributions of deformation inthe thermodynamic pressure arising from chemical processes such as mass transport, chemical reactions, andinterfacial effects. 7.
Acknowledgments
We are indebted to Professor Eliot Fried. We had many exhaustive discussions in which he gave usvaluable ideas, constructive comments, and encouragement. This publication was made possible in part bythe CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imag-ing Enterprise Future Science Platforms of the Commonwealth Scientific Industrial Research Organisation,CSIRO, of Australia. The European Union’s Horizon 2020 Research and Innovation Program of the MarieSk(cid:32)lodowska-Curie grant agreement No. 777778, and the Mega-grant of the Russian Federation Government(N 14.Y26.31.0013) provided additional support. Lastly, we acknowledge the support provided at CurtinUniversity by The Institute for Geoscience Research (TIGeR) and by the Curtin Institute for Computation.
Appendix A. Thermodynamically consistent continuum theory for the multicomponentCahn–Hilliard equations
A.1.
Larch´e–Cahn derivatives.
Let ϕ = { ϕ , . . . , ϕ n } (47)be a list of species concentrations and assume that the function F depends on ϕ such that F ( ϕ ) = F ( ϕ , . . . , ϕ n ) . (48)Constraint (5), with (4), implies that the set of concentrations ϕ must be 0 < ϕ α <
1. If we vary oneconcentration ϕ α while holding all others fixed violates the constraint (5). Thus, the conventional partialderivative on functions such as F , on which the constraint (5) is active, is not appropriately defined. Toovercome this shortcoming, Larch´e and Cahn [21] defined the following operation ∂ ( σ ) F ( ϕ ) ∂ϕ α = dd (cid:15) F ( ϕ , . . . , ϕ α + (cid:15), . . . , ϕ σ − (cid:15), . . . , ϕ n ) (cid:12)(cid:12)(cid:12) (cid:15) =0 (49)in which we choose any two concentrations ϕ α and ϕ σ from the set of variables. Then, we introduce an infin-itesimal change (cid:15) in ϕ α , which induces the opposite infinitesimal variation (cid:15) onto ϕ σ , while holding all othervariables unchanged. Thus, this definition satisfies (5) by construction while we express the concentration ϕ σ as ϕ σ = 1 − n (cid:88) α =1 α (cid:54) = σ ϕ α . (50)In multicomponent Cahn–Hilliard systems, we incorporate cross-diffusion gradient energy coefficients Γ αβ into the free-energy definition and obtain the following free-energy densityˆ ψ ( ϕ , grad ϕ ) := f ( ϕ ) + n (cid:88) α =1 n (cid:88) β =1 Γ αβ grad ϕ α · grad ϕ β . (51)Elliott & Garcke in [11] prove that multicomponent systems are well-posed when Γ αβ is positive definite,among other conditions. We show that this condition is sufficient but not necessary. To do so, we extend XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 15 the ideas of Larch´e–Cahn and define a constrained inner product on a constrained space. We consider a setof vectors { p α } subject to the following constraint n (cid:88) α =1 p α = , (52)and use the following inner product n (cid:88) α =1 n (cid:88) β =1 Γ αβ p α · p β . (53)Let each entry of Λ αβ be a single number κ . Thus, due to (52), { p α } is in the null space of Λ αβ , thatis, Null(Λ αβ ) = { p α } . Similarly, if each row of Λ αβ is given by the same entry κ β , we arrive to the sameconclusion. For any of these cases, we have that n (cid:88) α =1 n (cid:88) β =1 Γ αβ p α · p β = n (cid:88) α =1 n (cid:88) β =1 (Γ αβ + Λ αβ ) p α · p β . (54)We impose the constraint (52) with respect to the component σ to the quadratic form (53) to obtain n (cid:88) α =1 n (cid:88) β =1 Γ αβ p α · p β = n (cid:88) α =1 α (cid:54) = σ n (cid:88) β =1 β (cid:54) = σ (Γ αβ + Γ σσ − Γ ασ − Γ σβ (cid:124) (cid:123)(cid:122) (cid:125) Γ αβσ p α · p β . (55)We reinterpret this result as an inner product in an unconstrained space of dimension n − non-invertible mapping Γ αβ (cid:55)→ Γ αβσ defined asΓ αβσ := Γ αβ + Γ σσ − Γ ασ − Γ σβ . (56)Consequently, the problem is well-posed if Γ αβσ is positive definite. Moreover, Γ αβ can be indefinite withoutcompromising the well-posedness of the problem. Now, let Γ αβ be a diagonal matrix such thatΓ αβ = κ δ αβ . (57)From (54), we rewrite Γ αβ as Γ αβ = − κ (1 αβ − δ αβ ) , (58)where 1 αβ is a constant matrix populated by ones and δ αβ is the Kronecker delta, both of dimension n .Although the matrix (58) has a null diagonal, the mapping defined by (56) is identical to the one of thediagonal matrix (57) for all vectors that satisfy the constraint (52).A.2. Thermodynamics.
Here, we establish the first and second law of thermodynamics. First, we augmentthe species balances (6) ˙ (cid:37) α = − div α + s α + s α ext , (59)to consider an external mass supply s α ext as well as an internal one s α arising from chemical reactions.We treat chemical reactions in a similar fashion as Gurtin & Vargas [22]. Moreover, following Gurtin [23]and Cherfils et al. [24], see also [25, 26], we separate conservation statements from constitutive equations.Thus, we introduce the external power expenditure W ext to P done by the external microforces on P andmicrotractions on S to describe the thermodynamics of this system as follows W ext ( P ) := n (cid:88) α =1 (cid:90) P γ α ˙ ϕ α d v + (cid:90) S ξ α S ˙ ϕ α d a , (60)where n is the total number of species and ξ α S = ξ α · n is the α -th microtraction.The first law of thermodynamics states the energy balance between the interleaving of internal energyand the expenditure rate of the chemical (diffusion and reaction) power. The entropy production imbalance,or the second law of thermodynamics in the form of the Clausius-Duhem inequality, states that the rate of growth of the entropy is at least commensurate with the entropy flux and its supply. Thus, we can expressthese two laws as ˙ (cid:90) P (cid:37) ε d v = W ext ( P ) − (cid:90) S q · n d a + (cid:90) P (cid:37) r d v + n (cid:88) α =1 − (cid:90) S µ α α · n d a + (cid:90) P µ α s α ext d v , ˙ (cid:90) P (cid:37) η d v ≥ − (cid:90) S q ϑ · n d a + (cid:90) P (cid:37) rϑ d v, (61)where ε and η represent the internal-energy and entropy densities, respectively, q is the heat flux, r is theheat supply, and ϑ > s α to the energy balance (61).Using the external power expenditure (60), the microforce balance (59), and the constituent contentbalance (59), we can localize the first two laws of thermodynamics (61) to (cid:37) ˙ ε = n (cid:88) α =1 {− π α ˙ ϕ α + ξ α · grad ˙ ϕ α − α · grad µ α + µ α ( (cid:37) ˙ ϕ − s α ) } − div q + (cid:37) r,(cid:37) ˙ η ≥ − div q ϑ + (cid:37) rϑ . (62)Rewriting (62) , we obtain (cid:37) ˙ η ≥ − ϑ div q + 1 ϑ q · grad ϑ + (cid:37) rϑ . (63)We now define the free-energy density as ψ := ε − ϑη, (64)which allow us to rewrite the equation system in terms of ϑ and ψ . To employ this transformation, wemultiply (63) by ϑ and subtract the result from (62) to express the pointwise free-energy imbalance as (cid:37) ( ˙ ψ + ˙ ϑη ) + n (cid:88) α =1 { ( π α − (cid:37) µ α ) ˙ ϕ α − ξ α · grad ˙ ϕ α + α · grad µ α + µ α s α } + 1 ϑ q · grad ϑ ≤ . (65) Remark 2 (Alternative derivation–Principle of virtual power) . The definition of virtual power expenditureencompasses internal and external contributions. Internally to P , the power exerted by internal microforcesand the microstresses; while externally to P , the power effected by the external microforces on P andmicrotractions on S . This definition assumes that these contributions equilibrate each other, that is, V int ( P , χ α ) = V ext ( P ; χ α ) (66) where the definitions of the internal and external virtual powers are V int ( P ; χ α ) := n (cid:88) α =1 (cid:90) P ( − π α χ α + ξ α · grad χ α )d v (67) and V ext ( P ; χ α ) := n (cid:88) α =1 (cid:90) P γ α χ α d v + (cid:90) S ξ α S χ α d a , (68) where { χ α } is a set of n kinematically admissible fields. Finally, we apply the divergence theorem to (66) and use standard variational arguments to localize the balance of microforces (59) to the following. For amore general approach, see [27]. A.3.
Theory of reacting materials.
Theoretically, the total number m of possible independent chemicalreactions, where m ≥ n s ∈ N , is not arbitrary. We seek to fit our framework in the thermochemistry theoryof reacting materials (see, [28] and [29, 30]). Thus, we also postulate the indestructibility of the atomicsubstances n (cid:88) α =1 t αι s α m α = 0 , ≤ ι ≤ n a , (69) XTENDED LARCH´E–CAHN FRAMEWORK FOR REACTIVE CAHN–HILLIARD MULTICOMPONENT SYSTEMS 17 where n a ∈ N is the number of atomic substances making up all the components A , m α is the molecularweight of the α -th component, and t αι is a non-negative integer expressing the number of atoms of the ι -th atomic substance present in the α -th component. This postulate assumes that the atomic substanceare indestructible. Moreover, usually t αι is not a square matrix and rank( t αι ) = min( n, n a ). Finally, themaximum number of possible chemical reactions is m := n − rank( t αι ) . (70)In this setting, forward reactions and their reciprocal backward reaction are not independent. Thus, werepresent them as a single, effective, chemical reaction.A.4. Configurational stress and force.