A Cahn-Hilliard-Biot system and its generalized gradient flow structure

In this work, we propose a new model for flow through deformable porous media, where the solid material has two phases with distinct material properties. The two phases of the porous material follow a Cahn-Hilliard type evolution, with additional impact from both elastic and fluid effects, and the coupling between flow and deformation is governed by Biot's theory. This results in a three-way coupled system which can be seen as an extension of the Cahn-Larch\'e equations with the inclusion of a fluid flowing through the medium. The model covers essential coupling terms for several relevant applications, including solid tumor growth, biogrout, and wood growth simulation. Moreover, we show that this coupled set of equations follow a generalized gradient flow framework. This opens a toolbox of analysis and solvers which can be used for further study of the model. Additionally, we provide a numerical example showing the impact of the flow on the solid phase evolution in comparison to the Cahn-Larch\'e system.


Introduction
In this letter, we develop a general model with the ability to capture situations with flow through a deformable porous medium that changes character in terms of stiffness, permeability, compressibility, and poroelastic coupling strength due to Cahn-Hilliard-type phase changes in the solid matrix. There are several applications where this type of behavior exists. One example being solid tumor evolution, where it is argued that stress effects resulting from tumor growth impact the tumor evolution itself [1,2], and that stress can inhibit tumor growth [3,4,5]. Moreover, the elastic properties of the surrounding. matrix and the interstitial fluid pressure are elevated in most solid malignant tumors [6]. One can then consider the two-phase porous medium as cancerous and healthy cells with the surrounding extracellular matrix, and the fluid as the interstitial fluid. Similar models involving Cahn-Hilliard-type evolution of tumor growth can be seen in [7,8,9,10,11,12,13,14,15,16,17]. Additional applications of poroelastic media with solid phase changes range from biogrout to wood growth, where sapwood transforms to heartwood.
The proposed system is an extension of the Cahn-Hilliard model and the quasi-static linear Biot equations, where the Cahn-Hilliard contribution governs the solid phase changes in the system through a smooth phase-field variable, and the Biot equations govern flow and elasticity. The Cahn-Hilliard equation originates from the work of Cahn and Hilliard [18], where the interfacial free energy of a non-uniform composition was introduced to model phase separation. Coupling the Cahn-Hilliard model with elasticity, is often called the Cahn-Larché model due to its origination [19], and several applications have been considered with this model in mind, including li-ion batteries [20], and tumor evolution [9,8]. In this work, we include fluid in the system, which is assumed to flow through the poroelastic medium with Biot-type coupling between flow and elasticity [21].
We show that the resulting model has a generalized gradient flow structure, i.e., a dissipative system where the state of the system evolves with the negative gradient of its free energy. The extension to generalized gradient flows allows for non-quadratic, and partially degenerate, dissipation potentials, and there is currently an increasing interest in the mathematics of generalized gradient flows, both with respect to modeling [22,23], abstract analysis [24,25,26,27] and numerical solution strategies [27,28]. It is long known that the Cahn-Hilliard equation and single-phase flow through porous media can be written as standard gradient flows, and it was showed in [27] that the Biot equations have a generalized gradient flow 2 The derivation of the Cahn-Hilliard-Biot model We consider a saturated porous medium with one fluid phase, and two solid phases with distinct material properties. The solid phases are modeled by a diffuse interface approach of Cahn-Hilliard type, where surface tension, deformation of the solid material, and pore pressure are acting as driving forces.
Let the medium Ω ⊂ R d be a bounded domain, d the spatial dimension, and [0, T ] be a time interval where T denotes the final time. In the matrix, the smooth phase-field, ϕ : Ω × [0, T ] → [−1, 1], tracks the two phases ϕ = −1 and ϕ = 1. We consider linearized elasticity with infinitesimal displacement u, and ∇u 1, the pore pressure is denoted by p, and q is the fluid flux.

Balance laws
Balance laws are imposed for each of the three coupled systems. For the phase-field equation, we assume that the phase-change is conserved through a phase-field flux J and reactions R, The elastic behavior of the material is governed by a quasi-static force balance equation where σ denotes the stress tensor and f external body forces Finally, the fluid is assumed to follow a volume balance law with negligible density gradients, where θ is the volumetric fluid content which changes due to the fluid flux q and source S f .

Free energy
The system is then closed through its free energy together with appropriate constitutive relations. We assume that the energy can be decomposed into three parts; the regularized surface energy, containing chemical energy and interfacial energy between the solid phases, the elastic energy, and the fluid energy The regularized surface energy [18] is given as where deviations from pure phases are penalized through the double-well potential Ψ(ϕ), and transitions between phases are penalized by the second term which is related to the interfacial energy. Here, the parameter γ corresponds to interfacial tension between the phases and will account for adhesive and cohesive forces. The double-well potential takes minimal values in the two phases, ϕ = −1 and ϕ = 1, and is, in this work, given as We assume that the elastic energy takes the form that is typical to the Cahn-Larché equations, where ε(u) = 1 2 ∇u + ∇u is the linearized strain at displacement u. The second term, T (ϕ), is the eigenstrain at ϕ (often called stress-free strain, or intrinsic strain) which corresponds to the state of the strain tensor if the material was uniform and unstressed [29]. Moreover, it can be considered to account for swelling effects [20] and takes different values depending on the solid phase ϕ. Here, we consider the form T (ϕ) = ξϕI, where ξ is a swelling parameter. The elastic stiffness tensor C(ϕ), which can be anisotropic, depends on the phase-field.
Finally, we consider a natural extension of the classical fluid energy which is given as in [27] by where both the compressibility parameter M (ϕ) and the Biot-Willis coupling coefficient α(ϕ) depend on the phase-field ϕ.

Constitutive relations
Assuming that the phase-field follows Fick's law for non-ideal mixtures, the flux J is proportional to the negative gradient of the chemical potential where m(ϕ) is the chemical mobility. The chemical potential µ is defined to be the variational derivative of the free energy with respect to ϕ. Here, we denote the variational derivative of E with respect to y by δ y E, and standard computations yield where zero Neumann or periodic boundary conditions have been applied to ϕ, and According to thermodynamical principles [21], we define the stress tensor to be the rate of change of energy with respect to strain and the pore pressure p to be the rate of change of energy with respect to volumetric fluid content Finally, the flow through the porous medium is assumed to follow Darcy's law where the permeability κ(ϕ) is assumed to depend on the solid phase.
Combining the balance laws with the constitutive relations, and making the identification (14) in (12) and (13), the Cahn-Hilliard-Biot model becomes q + κ(ϕ)∇p = 0, equipped with suitable boundary and initial conditions.
3 The Cahn-Hilliard-Biot model as a generalized gradient flow In this section, we identify the proposed Cahn-Hilliard-Biot model (16)-(20) as a generalized gradient flow, which in contrast to regular gradient flows allow for non-quadratic and even degenerate dissipation potentials. By making this identification for the newly proposed model, a wide toolbox of well-posedness analysis [24,27], numerical error analysis [25,26] and numerical solution algorithms [27,28] are made available, which will be a valuable asset for further study. A generalized gradient flow takes the form where z is a state variable, R is a dissipation potential, E is the energy at state z, D x is the Gateaux gradient with respect to x, and P ext corresponds to external forces. Alternatively, one can reformulate the generalized gradient flow and split between states evolving with (z d ) and without (z df ) dissipation to get the constrained minimization problem subject to s d +∇·l = S, where R(∂ t z d , z d ) =R(F, z d ), ·, · is the canonical inner-product, and the balance law ∂ t z d + ∇ · F = S with flux F, and source S holds.
For the Cahn-Hilliard-Biot system, consider the state variables z = (ϕ, u, θ), the energy E(z) from (4), and the state-dependent dissipation potential with together with the conservation laws ∂ t ϕ + ∇ · J = R and ∂ t θ + ∇ · q = S f .
As the deformation is assumed to be dissipation free, the generalized gradient flow reads: Find ϕ, u, and θ such that (27) subject to η + ∇ · l = R and s + ∇ · v = S f with balance laws (25), P ext,e , w := Ω f · w dx and P ext,f corresponding to external forces related to the fluid (e.g., boundary conditions or gravitational force). Calculating optimality conditions, and substituting the phase-field flux J by the chemical potential µ through Fick's law (9), and the volumetric fluid content θ with the fluid pressure p through the relation (14), one obtains the variational form of the system (16)-(20).

Numerical example
Here, we present a numerical example that emphasizes the need for the Cahn-Hilliard-Biot model. We compare a simulation of the Cahn-Hilliard-Biot model with and without a pressure boundary condition acting as an external force (in order to enforce flow in the domain), with a Cahn-Larché simulation (Cahn-Hilliard coupled with only elasticity). The example clearly shows that when the fluid flow is dominant, it also plays a crucial role in the evolution of the phase-field. However, in regimes with little, to no flow, the phase-field is unaffected compared to the Cahn-Larché model. We consider a unit square domain where four circular shapes of phase ϕ = 1 are surrounded by phase ϕ = −1 initially, see Figure 1a,1e,1i. For both pressure and displacement, we apply zero initial data. The variational system (16)-(20) is discretized in time by a semi-implicit Euler method, where the deviation from fully implicit Euler is an application of the first order convex splitting method of the double-well potential Ψ(ϕ) as proposed in [30]. The three-way coupled nonlinear system is then solved by an iterative decoupling scheme, starting with the Cahn-Hilliard subsystem (16)- (17), then elasticity (18), and finally, flow (19)- (20). The Cahn-Hilliard subsystem (16)- (17) is discretized in space with bilinear rectangular finite elements for both phase-field ϕ and chemical potential µ, and the nonlinear equations are solved by a Newton method in each iterative decoupling-iteration. The flow subsystem (19)- (20) is discretized in space by lowest-order Raviart-Thomas elements, RT0, for the flux and constant elements for pressures, and the elasticity equation (18) is discretized with bilinear finite elements. We have used modules from the DUNE project, specifically dune-functions [31], for the implementation.
The material parameters can be found in Table 1, and the permeability κ(ϕ), compressibility M (ϕ), Biot-Willlis coefficient α(ϕ) and elasticity tensor C(ϕ) are depending on the phase-field through the interpolation function π(ϕ); κ(ϕ) as in [32], with the two elasticity tensors written in Voigt notation in two spatial dimensions. Zero Neumann boundary conditions are applied to both the phase-field and the chemical potential, while the displacement is equipped with zero Dirichlet conditions on the entire boundary. For the flow subsystem, we enforce a pressure drop from p = 0.25 to p = 0 from top to bottom while no-flow conditions are applied on the left and right parts of the boundary. In Figure 1a-1d, the phase-field function ϕ is plotted after a series of time steps for the Cahn-Hilliard-Biot model with a drop in pressure from p = 0.25 to p = 0 from top to bottom. In Figure 1e-1h the solution is plotted at the same time steps, but with zero pressure on the entire boundary, and similarly in Figure 1i-1l the plots are from a simulation of the Cahn-Larché system. We observe that when the flow is prominent in the simulation the phase-field is also significantly affected and takes a directional preference to that of the flow direction. When, on the other hand, the system merely is filled with a fluid that has no driving force in itself, the phase-field evolution is close to unaffected compared to the system without a fluid. We emphasize also that the system energies (including external forces) are decreasing over the scope of the simulation, as is expected from dissipative systems of gradient flow type. This is showed in Figure 1m, where the energy is a combination of the free energy of the system (4), and the external forces applied through the pressure boundary condition, E Tot = E(ϕ, u, p) − Γ Top p Top (q · n) dx, n being the outwards pointing normal vector.

Conclusions
The Cahn-Hilliard-Biot system was derived through balance laws and constitutive relations, i.e., Fick's law for the phase-field, and Darcy's law for the fluid flow. Key quantities are defined, following thermodynamical principles, as rates of change of the free energy. The equations feature a three-way coupling, and the impact from flow to the phase-field was showed to be significant through a numerical example; the phase-field does not only evolve as it would through the Cahn-Larché equations, but its evolution is aligned and magnified in the flow direction. Moreover, we showed that the system follows a generalized gradient flow framework and that the energy dissipates numerically as expected. By this, we lay the groundwork for a general model, showing numerical properties and highlighting important coupling terms, that can be further tailored and  studied depending on the specific application in mind.