Computational investigation of porous media phase field formulations: Microscopic, effective macroscopic, and Langevin equations

Abstract

We consider upscaled/homogenized Cahn–Hilliard/Ginzburg–Landau phase field equations as mesoscopic formulations for interfacial dynamics in strongly heterogeneous domains such as porous media. A recently derived effective macroscopic formulation, which takes systematically the pore geometry into account, is computationally validated. To this end, we compare numerical solutions obtained by fully resolving the microscopic pore-scale with solutions of the upscaled/homogenized porous media formulation. The theoretically derived convergence rate $\mathcal{O}({ϵ}^{1/4})$ is confirmed for circular pore-walls. An even better convergence of $\mathcal{O}({ϵ}^1)$ holds for square shaped pore-walls. We also compute the homogenization error over time for different pore geometries. We find that the quality of the time evolution shows a complex interplay between pore geometry and heterogeneity. Finally, we study the coarsening of interfaces in porous media with computations of the homogenized equation and the microscopic formulation fully resolving the pore space. We recover the experimentally validated and theoretically rigorously derived coarsening rate of $\mathcal{O}(t^{1/3})$ in the periodic porous media setting. In the case of critical quenching and after adding thermal noise to the microscopic porous media formulation, we observe that the influence of thermal fluctuations on the coarsening rate shows after a short, expected phase of universal coarsening, a sharp transition towards a different regime.

Introduction

Over the last decades, phase field modelling has received increasing interest for theoretical and computational investigation of physical, chemical and even experimental systems inspired by the work of Cahn and Hilliard [1]. However, the idea of diffuse interface modelling seems to go back to van der Waals [2]. The variational structure based on free energies allows for thermodynamic modelling of phase transitions [3], [4] and it serves as a predictive tool in engineering of fluid mechanics [5], multiphase flow [6], [7], [8], fuel cells [9], [10], batteries [11], and porous media [12]. Since many systems and applications involve strongly heterogeneous media, we refer to these by the general term Complex Heterogeneous Multiphase Systems. From a numerical point of view, strong heterogeneities lead to computationally high dimensional systems since the mesh size $h>0$ has to be chosen much smaller than the heterogeneity $ϵ>0$, i.e., $0<h\ll ϵ$. The heterogeneity parameter is defined by $ϵ=\frac{\ell }{\mathrm{\Lambda }}$ where ℓ denotes a material specific microscale, e.g., characteristic pore size, and Λ is the macroscopic size of the porous medium. As a consequence, an effective macroscopic phase field equation has been derived in [13], [14] that does not depend on such a restricting mesh constraint. In fact, a first attempt of extending the framework towards fluid flow is [15] albeit requiring specific assumptions such as flows with large Péclet number. Recently, diffuse interface formulations are gaining increasing interest in studying contact lines and droplets on solid substrates [16], [17], [18], [19]. Phase field formulations provide also a convenient computational alternative to sharp interface models. It has been applied in modelling dynamics of multicomponent vesicles [20] or as a computational tool for bubble dynamics [21]. In [22] the authors study the problem of surface diffusion based on a diffuse interface formulation and similarly in [23], the authors look at the problem of phase transition and coarsening on surfaces. Similarly, coarsening is considered for an interacting particle system in [24]. Another more recent and promising extension of Cahn and Hilliard's diffuse interface concept is the phase field crystal method [25] which takes atomistic information into account for modelling crystal growth as proposed in [26], [27]. The Ginzburg–Landau functional leads also to a mathematical theory in superconductivity where the study of minimizers and asymptotic limits is equally important, e.g. [28], [29].

In this article, we computationally investigate the models recently derived in [13], [14]. We introduce numerical methods for validating the rigorously derived error estimates from [30] and for studying dynamics of contactlines [16], [18], [31] in porous media. Finally, we validate the well-accepted physical phenomenon of coarsening of phase separating systems in strongly heterogeneous media [32], [33], [34]. It is well-known that late stages of first-order phase transitions in homogeneous domains [32], [35], [36] show a power law growth with exponent 1/3. This universal behaviour has been put on a rigorous basis in [33]. Our computations based on a full microscopic as well as an upscaled/homogenized phase field formulation recover this universal scaling. The structure of the paper is as follows: in the next paragraph we first recall the basic concepts and results of the Ginzburg–Landau/Cahn–Hilliard phase field theory. In Section 2, we introduce the different phase field formulations of interest, i.e., a microscopic, a novel effective macroscopic, and a microscopic description that accounts for thermal fluctuations. We provide numerical discretizations of these models and study convergence in Section 3. The influence of pore geometry/heterogeneities on the coarsening is the topic of Section 4.

Phase field modelling. Phase ordering/transition is generally described by a coarse grained local order parameter $\varphi (\mathbf{x},t):D\times (0,T)\to \mathbb{R}$, on a bounded domain $D\subset {\mathbb{R}}^d$ where d is the dimension of space and $0<T<\infty$ the maximum time of observation. We choose $D:={[0,1]}^2$ to be the unit square for our computations later on. The phase field variable ϕ is phenomenologically characterized by the Ginzburg–Landau/Cahn–Hilliard free energy in absence of external fields, that is,

$$F(\varphi )=\frac{1}{|D|}\underset{D}{\int }f(\varphi )d\mathbf{x},$$

here normalized by the Lebesgues measure $|D|=1$ of the domain D. The energy density $f(\varphi ):=f_L(\varphi )+\frac{{\lambda }^2}{2}{\left|\mathrm{\nabla }\varphi \right|}^2$ consists of the homogeneous free energy

$$f_L(\varphi )=a(\theta ){\varphi }^2+\frac{b}{2}{\varphi }^4,$$

with $a(\theta ):=a_0(\theta -{\theta }_c)<0$, ${\theta }_c$ a given critical temperature, and $a_0>0$, and $b>0$. The parameter $\lambda >0$ is proportional to the interfacial width. The gradient term in (1) allows for diffuse interfaces. A dynamic description of ϕ is generally obtained by minimizing (1) over time with the help of a gradient descent/flow, i.e., we are looking for solutions ϕ of

$$({\partial }_t\varphi ,v)=-\frac{{\delta }^{V}F(\varphi )}{{\delta }^V\varphi }=-{({\mathrm{\nabla }}_{\varphi }^{V}F(\varphi ),v)}_V,\text{for all }v\in V\text{\hspace{0.17em},}$$

where $(\cdot ,\cdot ):={(\cdot ,\cdot )}_{L^2}$ denotes the $L^2$-scalar product if $V=L^2$, and it is the $H^1$ semi-inner product ${(\cdot ,\cdot )}_V:=(\hat{\mathrm{M}}\mathrm{\nabla }\cdot ,\mathrm{\nabla }\cdot )$ for symmetric, positive definite tensors $\hat{\mathrm{M}}\in {\mathbb{R}}^{d\times d}$ if $V=H^1$. The Fréchet derivative $\frac{{\delta }^{V}F(\varphi )}{{\delta }^V\varphi }={({\mathrm{\nabla }}_{\varphi }^{V}F(\varphi ),v)}_V$ allows us to uniquely identify the functional derivative ${\mathrm{\nabla }}_{\varphi }^V$ by Riesz representation theorem [37, p. 163]. The conservation of mass is obtained for $V=H^1$ which leads to the well-known Cahn–Hilliard equation [1], that is,

$${\partial }_t\varphi =\mathrm{div}\left(\hat{\mathrm{M}}\mathrm{\nabla }(f_L^{\prime }(\varphi )-{\lambda }^2\mathrm{\Delta }\varphi )\right).$$

We note that the $V=L^2$-gradient flow leads to the Allen–Cahn equation, which is not mass conserving in difference to (4). The double-well character of the free energy (1) arises also in the regular solution theory in the form of the Flory–Huggins energy density [38], i.e.,

$$f_L(\varphi )=z{c}_I\varphi (1-\varphi )+k_B\theta (\varphi \ln \varphi -(1-\varphi )\ln (1-\varphi )),$$

where $k_B$ is the Boltzmann constant, z counts the number of bonds with neighbouring species and is called coordination number, and $c_I$ represents an mean field interaction energy. In this article, we will work with the well accepted double-well potential

$$f_L(\varphi )=a{\varphi }^2{(1-\varphi )}^2,a>0,$$

which reliably captures the phenomenological nature of the energy densities (5) and (2) but generally shows more stable behaviour in computations.

Coarsening and coarsening rates. The coarsening, e.g. [24], [32], [33], describes the time evolution of a characteristic length $L(t)$ which is defined as

$$L:=1/F(\varphi ),$$

where F represents the interfacial area per unit volume(/perimeter per unit volume) (1) and hence L has units of length. We note that Kohn and Otto [33] made a step change in the rigorous understanding by proving the time-averaged coarsening rates

$$\frac{1}{T}\underset{0}{\overset{T}{\int }}{F}^{2}dt\ge \frac{C}{T}\underset{0}{\overset{T}{\int }}{(t^{-1/3})}^{2}dt,\text{for }T\gg 1.$$

However, the following classical statement

$$L\le C{t}^{1/3}$$

still lacks a rigorous argument. In Section 4, we will computationally investigate the influence of pore geometries to the coarsening rates computed by the characteristic length L defined in (7).