The Lattice Boltzmann Equation (LBE) Book_LBM2017 is an attractive method to simulate flow and transport phenomena in several areas of science and engineering. Because of its local collision term and its ease implementation of bounce-back method, the LBE has been extensively applied in porous media literature for simulating two-phase flows and transport at pore scale Pan-Luo-Miller_CF2006 ; Genty-Pot_TRT_TiPM2013 ; Pot_etal_AdWR2015 (see Review_LBM_PorousMed_IJHMT2019 for a recent review). When the surface of separation between solid () and liquid () does not depend on time, it is sufficient to identify the nodes located at the interface and to apply the bounce-back method. However, when physico-chemical processes occur on the surface of solid, such as those involved in matrix dissolution or pore clogging, it is necessary of considering the free-boundary problem because the interface position is now a function of time. The general sharp interface model of dissolution and precipitation without fluid flows writes:
Eq. (1a) is the mass conservation of solute in bulk domains (where ), is the composition, is the diffusion coefficient of liquid () and solid (). Although is supposed to be zero in most of the dissolution studies, two diffusion coefficients and are considered here for more generality. Two conditions hold at the interface . The first one (Eq. (1b)) is the balance of advective and diffusive fluxes where is the normal velocity of interface. In that equation the right-hand side is the difference of diffusive fluxes between liquid and solid,
is the unit normal vector of interface pointing intoliquid, and is the composition of solid phase. The second interface equation (Eq. (1c)) is the Gibbs-Thomson condition that relates the driving force (left-hand side) giving birth to the interface motion (right-hand side). In literature, the most common form of is proportional to the difference between the interface composition and the solid composition : . Two terms contributes to the interface motion: the first one is the curvature-driven motion where is the curvature and is a capillary length coefficient. The second term is the normal velocity where is a kinetic coefficient representing the dissipation of energy.
In literature, the lattice Boltzmann methods often simulate an approximation of that sharp interface problem. In Kang_etal_WRR2007 , the equilibrium distribution functions are designed to fulfill the mass conservation at interface (Eq. (1b)). However, the method only simulates an approximation of the Gibbs-Thomson condition because the curvature term is neglected (). For instance in KANG20141049 , Eq. (1c) is replaced by an evolution equation of volume fraction of solid: the time variation of mineral volume is related to the reaction flux by (KANG20141049, , Eq. (5)) where is the molar volume of mineral and is the product of solid area times a kinetic coefficient. The model has been applied recently in Zhang_etal_ChemEngSci2022 for studying the influence of pore space heterogeneity on mineral dissolution. When the surface tension of material can be neglected, then the assumption hold. But in most general cases and the accurate position of interface must be computed while maintaining the two conditions Eqs. (1b)-(1c) at each time-step.
Alternative methods exist for simulating the interface tracking problem. In the “phase-field method”, a phase index is introduced to describe the solid matrix if (solid) and the pore volume if (liquid). The phase index varies continuously between those two extreme values () i.e. the method considers the interface as a diffuse zone. That diffuse interface is characterized by a diffusivity coefficient and an interface width . The interface, initially a surface, becomes a volumic region of transition between liquid and solid. The model is composed of two coupled Partial Derivative Equations (PDEs) defined on whole computational domain. The first equation describes the dynamics of the phase-field and the second one describes the dynamics of composition . Those two PDEs recover the sharp interface problem Eqs. (1a)–(1c) when . The Gibbs-Thomson condition Eq. (1c) is replaced by the phase-field equation which contains implicitly the curvature term . By “implicitly” we mean that the phase-field models always include the curvature-driven motion when they derive from a double-well potential.
Various phase-field models have already been proposed for simulating the processes of precipitation and dissolution Xu-Meakin_PhaseField-Preci_JChemPhys2008 ; Mai_etal_CorroSci2016 ; Bringedal_etal_MMS2020 ; Gao_etal_JCAM2020 . The main feature of those works is the model derivation from a free energy functional . Phase-field models that derive from such a functional have been successfully applied for solid/liquid phase change such as those encountered in crystal growth (e.g. Karma-Rappel_PRE1998 for pure substance and Echebarria_etal_PhysRevE.70.061604 ; Ramirez_etal_BinaryAlloy_PRE2004 for dilute binary mixture). For those applications, the functional depends on the phase-field and the temperature , which is an intensive thermodynamic variable. In spite of those successes for solid/liquid phase change, an issue occurs for models involving composition. The composition is an extensive thermodynamic quantity and the models do not necessarily insure the equality of chemical potentials at equilibrium. In order to fulfill that condition, the Kim-Kim-Suzuki (KKS) model KKS_model_PhysRevE.60.7186 introduces two fictitious compositions and in addition to the global composition . The two PDEs are formulated in and and the source term of depends on and . With a Newton method, those two compositions are explicitly computed inside the interface by imposing the equality of chemical potential (Provatas-Elder_Book_2010, , p. 126). That model has been applied for dissolution in Mai_etal_CorroSci2016 ; Gao_etal_JCAM2020 .
A formulation based on the grand-potential thermodynamic functional avoids that supplementary numerical stage. That approach, proposed in Plapp_PhysRevE.84.031601 , yields to a phase-field model that is totally equivalent to KKS model. That theoretical framework contains the construction of common tangent and insures the equality of chemical potential at equilibrium. As the same way they are derived from , the PDEs are established by minimizing . Hence, we retrieve the same features in the definition of . The density of grand-potential is composed of two terms. The first one, noted , contains the standard double-well potential and the gradient energy term of the interface. The second one, noted
is an interpolation of bulk grand-potentials. Those latter come from the Legendre transform of free energy densities . The main dynamical variables of are the phase-field and the chemical potential . The chemical potential is the conjugate variable to , and like temperature, it is an intensive thermodynamic quantity. Whereas it is inappropriate to make an analogy between and when deriving models, an analogy can be done between and . Thus, the asymptotics are quite similar for establishing the equivalence between sharp interface and phase-field models. That theoretical framework is already extended to study multi-component phase transformation Choudhury-Nestler_PRE2012 . It has been applied for dendritic electro-deposition in Cogswell_GrdPot_PRE2015 . The capability of grand-potential phase-field models to simulate spinodal decomposition is presented in Aagesen_etal_GrdPot_PRE2018 . In reference Simon_etal_GrdPot_CompMatSci2020 effects are presented of introducing elasticity with different interpolation schemes in grand-potential framework.
Contrary to solidification, the curvature term is often neglected in models of dissolution and precipitation. For instance in Xu-Meakin_PhaseField-Preci_JChemPhys2008 , the Gibbs-Thomson condition simply relates the normal velocity proportionally to . In Mai_etal_CorroSci2016 the normal velocity is only equal to the Tafel’s equation (Mai_etal_CorroSci2016, , Eqs. (2)-(3)). In Bringedal_etal_MMS2020 , the curvature term appears in the sharp interface model but the coefficient in front of the curvature is considered very small. As already mentioned, the curvature-driven motion is always contained in the phase-field model. If that motion is undesired in simulations, it is necessary to add a counter term in the phase-field equation as proposed in the pioneer work Folch_etal_PRE1999 . The counter term has been included for interface tracking in Allen-Cahn equation in reference Sun-Beckermann_JCP2007 . For two-phase flows a “conservative Allen-Cahn” equation has been formulated in Chiu-Lin_JCP2011 and coupled with the incompressible Navier-Stokes equations. For dissolution, the same term has been considered in the phase-field equation of Xu-Meakin_PhaseField-Preci_JChemPhys2008 . Here, the effect of the counter term is presented on the dissolution of 2D porous medium. That term has an impact on the shape of porous medium and the heterogeneity of composition inside the solid phase.
|[mol].[L]||Global concentration depending on and|
|Eq. (3)||[E].[L]||Grand-potential density of interface|
|see Tab. 2||[–]||Double-well potential of minima and|
|Eq. (4)||[E].[L]||Interpolation of bulk grand-potential density|
|[–]||Index for bulk phases: solid and liquid|
|[L].[E].[T]||Mobility coefficient of the interface|
|[–]||Hyperbolic tangent solution|
|[E].[L]||Coefficient of gradient energy term|
|[E].[L]||Height of double-well function|
|[E].[L]||Free energy density of bulk phases|
|[–]||Compositions for which is minimum|
|[E].[L]||Grand-potential density of each bulk phase|
|[E].[L]||Curvature of quadratic free energies|
|[E].[L]||Reference volumic energy for dimensionless quantities|
|[E].[L]||Difference of minimum values of free energy densities|
|[–]||Dimensionless grand-potential of bulk phases|
|[–]||Dimensionless free energy densities|
|[–]||values of in bulk phases: and|
|[L]||Interface width of -equation|
|[L].[T]||Diffusivity of -equation|
|[–]||Coupling coefficient of -equation|
|[–]||Unit normal vector of interface|
|[–]||Global composition depending on and|
|[–]||Dimensionless chemical potential|
|[–]||Equilibrium chemical potential of interface Eq. (29a)|
|[–]||Coexistence (or equilibrium) compositions of each phase|
|Eq. (29b)||[–]||Interpolation of coexistence compositions and|
|[L].[T]||Diffusion coefficient of solid () and liquid ()|
|see Tab. 2||[–]||Interpolation function of derivative zero for and|
|see Tab. 2||[–]||Interpolation function for|
|see Tab. 2||[–]||Interpolation function for diffusion coefficients|
|Eq. (25)||[–]||Source term of phase-field equation|
|[T]||Phenomenological counter term|
|Eq. (32)||[L].[T]||Phenomenological anti-trapping current|
|[–]||Coefficient of anti-trapping current|
|[L].[T]||Normal velocity of interface|
|[L]||Capillary length in Gibbs-Thomson condition|
|[L].[T]||Kinetic coefficient in Gibbs-Thomson condition|
|[–]||Ratio of diffusion|
In this paper, we derive in Section 2 a phase-field model based on the grand-potential functional for simulating dissolution and precipitation processes. In Section 2.1, the phase-field equation is presented without counter term for keeping the curvature-driven motion. Next in Section 2.2, the counter term is included in the phase-field equation which is reformulated in conservative form. Although the second main dynamical variable is the chemical potential , we use in this work a mixed formulation between the composition and the chemical potential (Section 2.3). The reason of this choice is explained by a better mass conservation when simulating the model. For sake of simplicity, the link to a thermodynamic database is not considered in this work. The grand-potential densities of each bulk phase derive from two analytical forms of free energy densities. We assume they are quadratic with different curvatures and for each parabola (Section 2.3). Next, Section 2.4 is dedicated to a discussion about the relationships of phase-field parameters , and with the sharp interface parameters, the capillary length and the kinetic coefficient . Those relationships will give indications to set the coupling parameter in simulations.
The model is implemented in LBM_saclay, a numerical code running on various High Performance Computing architectures. With simple modifications of compilation flags, the code can run on CPUs (Central Process Units) or GPUs (Graphics Process Units) Verdier_etal_CMAME2020 . The LBM schemes of phase-field model are presented in Section 3. A special care is taken for canceling diffusion in solid phase and accounting for the anti-trapping current. Validations are carried out in Section 4. LBM results are compared with analytical solutions for precipitation and next for dissolution. The first case is performed for (Section 4.1) to show the discontinuity of composition on each side of interface. The second one presents for (Section 4.2) the impact of anti-trapping current on the profiles of composition. Finally, in Section 5, we present the dissolution of a porous medium characterized by microtomography. Two simulations compare the effect of the counter term on the composition and the shape of porous medium.
2 Phase-field model of dissolution/precipitation
The purpose of this Section is to present the phase-field model of dissolution and precipitation. Its derivation introduces a great quantity of mathematical notations. The reason is inherent to the whole methodology: the diffuse interface method, which originates from out-of-equilibrium thermodynamics, recovers the sharp-interface model through the matched asymptotic expansions. Each keyword introduces its own mathematical notations. All those relative to physical modeling are summarized in Tab. 1.
In Section 2.1, we remind the theoretical framework of grand-potential , and we present the general evolution equations on and . Section 2.2 reminds the equilibrium properties of phase-field equation and introduces the counter term for canceling the curvature-driven motion. Equations on and require the densities of grand-potential for each phase and . In Section 2.3 their expressions are derived from analytical forms of free energies and . The phase-field model will be re-written with a mixed formulation between and with the compositions of coexistence and the equilibrium chemical potential. Finally, in Section 2.4 a discussion will be done regarding the links between phase-field model and free-boundary problem.
2.1 General equations on and in grand-potential theoretical framework
The grand-potential is a thermodynamic functional which depends on the phase-field and the chemical potential , two functions of position and time . In comparison, and the composition are two main dynamical variables of free energy . The functional of grand-potential contains the contribution of two terms:
The first term inside the brackets is the grand-potential density of interface which is defined by the contribution of two terms depending respectively on and :
In Eq. (3), the first term is the double-well potential and is its height. The second term is the gradient energy term which is proportional to the coefficient . A quick dimensional analysis shows that the physical dimension of is an energy per volume unit ([E].[L]) and has the dimension of energy per length unit ([E].[L]). Those two contributions are identical for models that are formulated with a free energy functional . The mathematical form of the double-well used in this work will be specified in Section 2.2.
In Eq. (2), the second term interpolates the grand-potential densities of each bulk phase and by:
where is an interpolation function. It is sufficient to define it (see Section 2.4) as a monotonous function such as and in the bulk phases with null derivatives (w.r.t. ) . In this work we choose
and its derivative w.r.t. is
With that convention, if then and if then .
In this paper, we work with the dimensionless composition describing the local fraction of one chemical species and varying between zero and one. It is related to the concentration (physical dimension [mol].[L]) by where is the molar volume of ([L].[mol]). For both chemical species, the molar volume is assumed to be constant and identical. In the rest of this paper will appear in the equations for reasons of physical dimension, but it will be considered equal to for all numerical simulations.
The concentration is defined by an interpolation of derivatives of and w.r.t. . Each derivative defines the concentration of bulk phase and .
In Eq. (4), the grand-potential densities of each bulk phase and are defined by the Legendre transform of free energy densities and :
where . Finally, the phase-field equations are obtained from the minimization of the grand-potential functional . The most general PDEs write (see (Plapp_PhysRevE.84.031601, , Eq. (43) and Eq. (47))):
Eq. (8a) is the evolution equation on which tracks the interface between solid and liquid. The phase-field equation is derived from where is a coefficient of dimension [L].[E].[T]. The equilibrium properties of that equation are remembered in Section 2.2. The derivative of the double-well function w.r.t. is noted . Compared to the model of reference Plapp_PhysRevE.84.031601 , we notice the opposite sign of the last term because our convention is for solid and for liquid. In the reference, is solid and is liquid and the interpolation function is opposite. In order to make appear the diffusivity coefficient of dimension [L].[T], the coefficient can be put in factor of the right-hand side. In that case, the second term is multiplied by whereas the last term is divided by .
Eq. (8b) is the evolution equation on chemical potential . It is obtained from the conservation equation where the diffusive flux is given by
. The time derivative term has been expressed by the chain rule. The function , called the generalized susceptibility, is defined by the partial derivative of with respect to . For most general cases, the coefficient is the diffusion coefficient which depends on and . Here we assume that the diffusion coefficients and are only interpolated by , i.e. . Actually, in section 2.3, that equation on will be transformed back to an equation on (or ) for reasons of mass conservation in simulations. Eq. (6) will be used to supply a relationship between and .
For simulating Eqs. (8a) and (8b), it is necessary to define the grand-potential densities of each bulk phase and . They both derive from Legendre transforms (Eq. (7)) which require the knowledge of free energy densities and . The free energy densities and depend on the phase diagram of chemical species (or materials), temperature and number of species involved in the process (binary or ternary mixtures). When the model is implemented in a numerical code coupled with a thermodynamic database , those values are updated at each time step of computation. A method for coupling a phase-field model based on the grand-potential with a thermodynamic database is proposed in Choudhury_etal_COSSMS2015 . A coupling of a phase-field model with the “thermodynamics advanced fuel international database” is presented in Introini_etal_JNM2021 with OpenCalphad Sundman_etal_IMMI2015 ; Sundman_etal_CompMatSci2015 . In this work, we assume in Section 2.3 that the densities of free energies and are quadratic.
The variational formulation based on the grand-potential yields to evolution equations on and (Eqs. (8a)-(8b)). Two ingredients miss in those equations: the first one is the counter term and the second one is the anti-trapping current . In our work, both are not contained in the definition of grand-potential and have no variational origin. The counter term has been derived in Folch_etal_PRE1999 . It is used in the phase-field equation (Section 2.2.2) to make vanish the curvature-driven motion. The anti-trapping current has been derived in Karma_AntiTrapping_PRL2001 . It is used in the chemical potential equation (Section 2.3.4) to cancel spurious effects at interface when the diffusion is supposed to be null in the solid. Their use is justified by the matched asymptotic expansions carried out on the phase-field model. The links between the phase-field model and the free-boundary problem will be discussed in Section 2.4.
2.2 Equilibrium properties of phase-field equation
The phase-field equation Eq. (8a) has the same structure as those derived from functionals of free energy. Hence, the equilibrium properties such as the hyperbolic tangent solution , the interface width and the surface tension remain the same. Those equilibrium properties are reminded in Section 2.2.1 with one particular choice of double-well potential . This is done for two reasons. The phase-field equation is written with “thermodynamic” parameters , and . The phase-field equation is re-written with “macroscopic” parameters , and the dimensionless coupling coefficient because they are directly related to the capillary length and kinetic coefficient of sharp interface model. The equilibrium properties are also necessary for introducing in Section 2.2.2 the kernel function and the counter term .
2.2.1 Hyperbolic tangent solution , width and surface tension
When the system is at equilibrium, the construction of common tangent hold and the chemical potential is identical in both phases of value . The construction of common tangent is mathematically equivalent to . When the two phases are at equilibrium, we define the corresponding compositions of coexistence (or equilibrium) by and for solid and liquid respectively. Hence, the last term proportional to in Eq. (8a) vanishes at equilibrium and the time derivative is zero (). We recognize the standard equilibrium equation for interface i.e. in one dimension . After multiplying by , the first term is the derivative of and the second term becomes a derivative of the double-well w.r.t. . After gathering those two terms inside the same brackets, it yields:
In this work we define the double-well by
for which the two minima are and and its derivative w.r.t. is:
For that form of double-well, the solution of Eq. (9) is the usual hyperbolic tangent function
where the interface width and the surface tension are defined by
We can check that the square root of ratio is homogeneous to a length as expected for the physical dimension of the width . Moreover, the square root of the product is homogeneous to an energy per surface unit as expected for the surface tension . The two relationships Eq. (12a) can be easily inverted to yield
As a matter of fact, the double-well function Eq. (10a) is a special case of other popular choices of double-well. For example in two-phase flows of immiscible fluids, the double-well is Lee-Lin_JCP2005 with for which the two minima are and . For that form of double-well, the equilibrium solution is , the surface tension is and the interface width is . Eqs. (11) and (12a) can be recovered by setting and . Another popular choice of double-well is Zheng_etal_LargeDensityRatio_JCP2006 for which the two minima are . Once again, that double-well function is a particular case of the previous one by setting and . The equilibrium solution writes , the surface tension is and the interface width . In this work the choice of Eq. (10a) is done by simplicity.
2.2.2 Removing the curvature-driven motion in Eq. (8a)
Another useful relationship that derives from Eq. (9) is the kernel function . The square root of the term inside the brackets yields where the coefficient was replaced by the interface width with Eq. (12b) (). Thus, with a double-well function defined by Eq. (10a), the kernel function writes:
For canceling the curvature-driven interface motion, a counter term is simply added in the right-hand side of the phase-field equation. The counter term is proportional to the interface diffusivity , the curvature and the kernel function . The curvature is defined by where is the unit normal vector of the interface
In Section 2.4.3, we check that adding such a counter term in the phase-field equation cancels the curvature motion in the Gibbs-Thomson equation. In order to write the phase-field equation in a more compact form, we remark that the second term involving the derivative of the double-well is equivalent to
where the definition of the curvature has been applied for the second term. The right-hand side of Eq. (15b) is and by using the kernel function the phase-field equation writes
where . In simulations of Sections 4 and 5 two versions of the phase-field equation are used: Eq. (8a) when the curvature-driven motion is desired and Eq. (16) when that motion is undesired. When the source term of that equation is null, and when an advective term is considered, Eq. (16) is the conservative Allen-Cahn equation that is applied for interface tracking of two immiscible fluids Chiu-Lin_JCP2011 ; Fakhari_etal_JCP2017 .
2.3 Phase-field model derived from quadratic free energies
The source terms of Eqs. (8a) and (8b) contain the bulk densities of grand-potential and . They need to be specified. Here, we work with analytical expressions which define explicitly and as functions of . The main advantage of that choice is to simplify their expressions by involving several scalar parameters representative of the thermodynamics. The densities of grand-potential are defined by the Legendre transform of free energy densities and . In Plapp_PhysRevE.84.031601 , several choices for are proposed in order to relate the grand-potential framework to the well-known models derived from free energy. The simplest phenomenological approximation is a quadratic free energy for each phase and :
where , of physical dimension [E].[L], are the curvature of each parabola and are two values of composition for which are minimum of values .
In this Section, all terms of Eqs. (8a) and (8b) involving are simplified with that hypothesis of Eq. (17). First, Section 2.3.1 deals with the difference of grand-potential densities which will be written with the dimensionless chemical potential and the thermodynamical parameters , and of Eq. (17). Section 2.3.2 introduces the coexistence compositions of interface and the equilibrium chemical potential . The difference will be re-expressed with , and . In Section 2.3.3 the composition equation is re-written with a mixed formulation between and , and in Section 2.3.4 the anti-trapping current will be formulated as a function of . Finally, the complete model is summarized in Section 2.3.5.
2.3.1 Difference of grand-potential densities in -equation
We start with the difference where the chemical potential is defined by (for ). By inverting those relationships to obtain as a function of , the Legendre transforms of each bulk phase yield the grand-potential densities as function of (see intermediate steps in Plapp_PhysRevE.84.031601 ):
Before going further we set and we define the quantity (dimension [E].[L]) for introducing the dimensionless quantities , and by (with ), and . With those reduced variables, the difference writes
Finally, if we define the dimensionless coefficient of coupling by , the last term of Eq. (8a) writes
where for future use we have set defined by:
2.3.2 Coexistence compositions and chemical potential of equilibrium
In Eq. (21), and are two specific values of for which the quadratic free energies and are minimum. A close link exists between and the coexistence (or equilibrium) compositions . Two relationships allow deriving them: the first one is the equality of chemical potential :
and the second one is the equality of grand-potential densities
where and has been defined in Section 2.3.1. Finally, those two conditions yield two simple relationships between and the parameters and :
In binary case this couple of coexistence compositions is unique, and the mathematical model can be re-defined with and . More precisely in Eq. (21), and are replaced with and by using Eqs. (24a)-(24b). In addition, the ratio is simply replaced by . The source term simplifies to
Here the source term has been formulated with and provided that . If the relationships between and are more complicated because they are solutions of second degree equations. That case will be studied in a future work. Finally, we can relate the compositions to the chemical potential . By using the definition and the dimensionless notation , we obtain and . Those relationships will be useful in Section 4.
2.3.3 Mixed formulation and closure relationship between and in -equation
Even though the equation on chemical potential (Eq. (8b)) could be directly simulated, we prefer using a mixed formulation that involves both variables and . The time derivative is expressed with and the flux is expressed with . The advantage of such a formulation, inspired from (Bayle_PhD2020, , p. 62), is explained by a better mass conservation. With the chain rule, the PDE on (Eq. (8b)) is transformed back to the diffusion equation . Although, the diffusion coefficient is a function of in general cases, here we assume that it is only a function of i.e. . It is relevant to define with ) because the interpolation function
appears naturally during the asymptotic analysis of Section2.4 when switching to a dimensionless timescale. In addition, the coefficient is defined by where is defined by Eq. (27) when the free energies are quadratic. When that coefficient is simply equal to (see Eq. (27)). Finally, with and , the composition equation writes:
The closure equation between and is simply obtained with Eqs. (6) and (18) for expressing the composition . In Eq. (6), the interpolation function can be replaced by another one . The form of is discussed below. The closure equation writes:
Next, by inverting Eq. (27), we find a relationship that relates the dimensionless chemical potential to compositions , and :