## Abstract

This work shows that the Cahn–Hilliard theory of diffusive phase separation is related to an intrinsic *mixed variational principle* that determines the rate of concentration and the chemical potential. The principle characterizes a canonically compact model structure, where the two balances involved for the species content and microforce appear as the Euler equations of a variational statement. The existence of the variational principle underlines an *inherent symmetry* in the two-field representation of the Cahn–Hilliard theory. This can be exploited in the numerical implementation by the construction of time- and space-discrete *incremental potentials*, which fully determine the update problems of typical time-stepping procedures. The mixed variational principles provide the most fundamental approach to the finite-element solution of the Cahn–Hilliard equation based on *low-order* basis functions, leading to monolithic *symmetric algebraic systems* of iterative update procedures based on a linearization of the nonlinear problem. They induce in a natural format the choice of *symmetric solvers* for Newton-type iterative updates, providing a speed-up and reduction of data storage when compared with non-symmetric implementations. In this sense, the potentials developed are believed to be fundamental ingredients to a deeper understanding of the Cahn–Hilliard theory.

## 1. Introduction

The modelling of phase separation is crucial for the understanding of the formation of microstructures in materials. The Cahn–Hilliard equation (2.10), attributed to Cahn & Hilliard [1] and Cahn [2], characterizes important qualitative features of the evolution of two-phase systems in material science. It governs the evolution of a *concentration* of species (order parameter) as a conservation law that describes the transport of the species between unit cells. The phase separation is achieved by a combination of down- and up-hill diffusion, where material moves in and against the direction of the concentration gradient, driven by the gradient of a *chemical potential*. These effects are governed by the competition of two terms in the free energy of the system, a non-convex *configurational energy* function of the concentration and a *surface energy* function of the concentration gradient. Initial stages of the phase separation process minimize the configurational energy by driving the concentration into minima related to the pure phases. Subsequent stages minimize the surface energy by forming larger single-phase domains, which successively absorb smaller phase regions, i.e. minimizing the interface area between the two phases.

### (a) A mixed variational principle for the phase field evolution

The fundamental constitutive structure behind the Cahn–Hilliard equation has been outlined in Gurtin [3] in the context of extended continuum theories based on *microforce balances*. This comprises a rational approach to the derivation of the Cahn–Hilliard equation by combining the *two balances* for the species content and the microforce with *constitutive assumptions* for species flux, microforce and microtraction. Related to the two balances, the concentration field and the chemical potential appear in a natural format as the two primary variables. In this sense, a decomposition of the fourth-order Cahn–Hilliard equation into two second-order equations is consistent with an embedding into standard structures of continuum mechanics. In this context, the Cahn–Hilliard equation may be related to gradient-extend continuum formulations, such as gradient damage, gradient plasticity or phase field theories. These models can be motivated in an elegant format based on microforce balance equations. A rigorous variational framework of these non-standard continuum models based on *rate-type* and *incremental variational principles* was outlined recently by Miehe [4] and applied to a wide range of problems (e.g. [5–7]). However, an additional difficulty appears in the construction of related variational principles for the Cahn–Hilliard theory owing to the *dissipative transport equation*. This calls for a completely new variational framework. With regard to this challenge, the Cahn–Hilliard theory provides the simplest application for the construction of variational principles for *problems which include dissipation effects owing to transport phenomena.* It can be considered as the prototype example for the construction of variational principles for complex coupled problems in continuum mechanics, which include diffusive mass or fluid transport.

This work shows that the decoupled Cahn–Hilliard structure is related to a *mixed variational principle* for the evolution problem, formulated in terms of the rate of the concentration and the chemical potential. Its Euler equations are two partial differential equations, which express the two balances of species content and microforce in terms of two constitutive functions for the energy storage and the dissipation potential, including the Neumann-type boundary conditions. The variational principle *highlights a symmetry* of the evolution problem, characterized by the rate concentration and the chemical potential. It is considered as a canonical ingredient of the Cahn–Hilliard theory that can directly be exploited in the numerical treatment.

### (b) Numerical implementation of the mixed variational principle

A straightforward implementation of the Cahn–Hilliard equation with its fourth-order gradients in the concentration requires higher order interpolations, i.e. basis functions that are globally C1-continuous. With this regard, *finite difference approaches*, such as applied by Choo *et al.* [8] provide straightforward solutions on simple regular and non-moving domains. However, for complex non-regular domains or situations where the Cahn–Hilliard equation is defined on deforming solids, *finite-element schemes* are the appropriate discretization techniques. We refer to the works of Elliott & French [9], Barrett *et al.* [10], Stogner *et al.* [11], Rajagopal *et al.* [12] and Zhang *et al.* [13]. The treatment of C1-continuous finite-element spaces requires an increase in the polynomial degree and favours *isogeometric methods* as outlined in Gomez *et al.* [14] and Cottrell *et al.* [15]. We also refer to the implementation of the Cahn–Hilliard equation by *discontinuous Galerkin methods* outlined in Wells *et al.* [16]. To circumvent the complexities associated with higher order finite-element spaces, the fourth-order Cahn–Hilliard equation can be rewritten as a coupled system of two second-order partial differential equations, reducing the continuity requirement to C0. We refer to the works of Barrett *et al.* [10], Feng & Prohl [17], Elliott & French [9] and Zhang *et al.* [13]. An alternative formulation in this direction has been considered by Ubachs *et al.* [18] and Kuhl & Schmid [19] based on the introduction of local and non-local concentrations as primary unknowns at each node. However, all of these formulations do not exploit a variational principle for the coupled problem. As a consequence, the formulations are not symmetric in the chosen primary fields.

By contrast, this work develops, based on the continuous rate-type variational principle mentioned above, a new time-discrete incremental variational principle that provides a *potential for the incremental update problem*. The subsequent space discretization of the coupled problem by C0 elements with nodal degrees for the concentration and the chemical potential results in a canonical *space–time-discrete incremental potential*, whose first and second derivatives govern the finite-element residual and tangent arrays. The *symmetric structure* of the algebraic system resulting from the nonlinear incremental problem is a consequence of the proposed variational principle. In this sense, we believe that our mixed variational principle provides the most fundamental structure of the Cahn–Hilliard theory.

The work is organized as follows. Section 2 summarizes the classical approach to the Cahn–Hilliard equation based on the notion of a microforce balance. This provides a clear classification of balance equations and constitutive equations. It is shown in §3 to be consistent with variational principles, which are developed first in the continuous rate-type setting, then in a time- and finally in a space-discrete formulation. Finally, §4 considers a specification of the constitutive functions to a standard model problem and investigates numerical simulations of a phase separation processes based on the developed variational approach.

## 2. Generalized Cahn–Hilliard theory of species diffusion

### (a) Primary fields, species flux and microforce traction

The boundary value problem for Cahn–Hilliard-type diffusion is considered as a coupled two-field problem, characterized by the *species concentration* and the *chemical potential fields*
2.1The dimensionless concentration field *c*∈[0,1] represents at a point of the solid the fraction of species referred to a given reference value *r*, such that *cr* is the number of species per unit volume. The chemical potential *μ* drives the flux of the species. Considering a control volume , we have a chemical microtraction *ξ* and the species outflux *h* on the surface element *d*** a** of , which both depend linearly on the outward normal

**2.2through the**

*n**microforce traction vector*and the

*species flux vector*, see figure 1. Both have to be related to the primary fields defined in (2.1) by constitutive functions.

### (b) Balance equations of Cahn–Hilliard-type diffusion

Following the approach of Gurtin [3], the balance equations related to the Cahn–Hilliard theory are formulated for a part , which is considered as a control volume. The *balance of chemical microforce* and the *conservation of species content* read
2.3where *g* is the chemical microforce. Inserting definitions (2.2) for the microforce traction *ξ* and the species out-flux *h*, using the Gauss theorem and a standard localization argument, *two local balance equations* defined on the domain are obtained for the species content and the chemical microforce
2.4

### (c) Principle of irreversibility and constitutive equations

#### (i) Dissipation principle

The formulation of the constitutive equations must be consistent with the principle of irreversibility, i.e. the second axiom of thermodynamics. Let *ψ* be the *free energy density* per unit volume stored at a spatial point . Then, the evolution of the energy in the part must be less than the *power of the chemical external actions* on . This is expressed by the *global dissipation postulate*
2.5The first term is the power of chemical microforces. Note again that given external source terms are not taken into account. The last term characterizes the in-flux of energy owing to the species transport, driven by the chemical potential *μ*. Inserting definitions (2.2) for the microtraction *ξ* and the species out-flux *h*, using the Gauss theorem, the standard localization argument and inserting the balance equations (2.4) gives the *local dissipation postulate*
2.6The dissipation can be split up into a part owing to local actions and a part owing to the species diffusion .

#### (ii) Energy storage function

The constitutive equations are constructed such that dissipation condition (2.6) is *a priori* satisfied for all processes. Assuming a local theory of the grade one, the free energy is assumed to depend on the primary variables (2.1) and its first gradients . Demanding consistency with the dissipation principle in (2.6), cannot be a function of the chemical potential *μ* and its gradient ∇*μ*. We obtain the reduced form . Insertion into (2.6) and application of Coleman’s argument then gives the constitutive equations for the microforce and the microtraction
2.7consistent with the demand for the local dissipation in (2.6).

#### (iii) Dissipation potential function

In order to describe the flux of the species, we define a dissipation potential function that depends on the negative gradient −∇*μ* of the chemical potential at a *given* state {*c*} of concentration. Then the species flux vector is defined by the constitutive equation
2.8The demand for the dissipation owing to diffusion in (2.6) is then satisfied if is a *convex* function with respect to the argument −∇*μ*, positive and zero for −∇*μ*=**0**.

### (d) Summary of equations and boundary conditions

In summary, the initial-boundary value problem of the Cahn–Hilliard-type diffusion is governed by the coupled set of balance and constitutive equations
2.9The combination of these five ingredients finally gives the fourth-order *Cahn–Hilliard equation*
2.10when suitable constitutive functions defined in (4.1) and (4.4) for and are inserted.

In addition, initial and boundary conditions have to be defined. To this end, the surface of the domain is decomposed according to the two primary fields defined in (2.1) by , with and . We consider Dirichlet- and Neumann-type boundary conditions for the species concentration and the microtraction
2.11as well as for the chemical potential and the species flux
2.12with *prescribed* concentration and species outflux , chemical potential and microtraction , respectively. The above pairing of the Dirichlet and Neumann conditions already accounts for the underlying variational principle developed below. The initial condition for the concentration is
2.13

## 3. Variational principles for the evolution problem

### (a) Governing field equations for Cahn–Hilliard-type diffusion

The governing coupled field equations *combine* the balance equations with the constitutive equations in (2.9). First, the insertion of (2.9)_{3,4} into (2.9)_{2} gives the constitutive representation of the *microforce balance*
3.1in terms of the stored energy function . Next, the insertion of (2.9)_{5} into (2.9)_{1} gives the constitutive representation for the *balance of species content*
3.2in terms of the dissipation potential function . Here, we have used the *variational* or *functional derivatives* for a more compact notation.^{1} Two equations (3.1) and (3.2), together with Dirichlet conditions (2.11) and (2.12), provide a modelling framework for the Cahn–Hilliard-type diffusion problem. They are in what follows shown to be related to the Euler equation of a variational principle.

### (b) Stored energy, dissipation potential and load functionals

We introduce two functionals related to the energy storage and the dissipative transport mechanisms, respectively, both depending on the primary fields introduced in (2.1). The *stored energy functional* depends on the concentration field *c*
3.3and characterizes the chemical energy stored in the domain . It is governed by the constitutive energy storage function introduced above.

The suitable definition of a *dissipation potential functional* is the key ingredient to the formulation of a rate-type variational principle for the Cahn–Hilliard-type diffusion problem. It is assumed to be determined by the variational principle
3.4depending on the rate of the concentration, where we introduced the canonical dissipation potential function for later use. It contains an external load functional
3.5of prescribed species flux related to Neumann condition (2.12). Observe that when taking the variation based on the admissible functions , which satisfy the homogeneous form of Dirichlet condition (2.12), the Euler equation of the optimization problem (3.4) yields constitutive form (3.2) for the *balance of species content* including the Neumann condition for the species flux. Definition (3.4) of the dissipation functional , governed by the constitutive function , provides the key approach for incorporating dissipative transport phenomena into a rate-type variational framework. The variational problem (3.4) can be viewed as a *generalized Legendre transformation*, where the chemical potential *μ* plays the role of a *mixed variable* dual to the rate of concentration. The formulation extends the local form to a global definition that accounts for boundary conditions of the dual variable *μ*.

Finally, related to Neumann boundary condition (2.11), the second external load functional is introduced as
3.6that contains a prescribed microtraction . For simplicity, all external actions are assumed to be dead loads, i.e. to be independent of the primary fields *c* and *μ*.

### (c) The canonical minimization principle

Based on the energy and the dissipation potential functionals *E* and *D* introduced and the external work functional above, the *rate-type potential*
3.7is postulated at given state {*c*} and time *t*. We rewrite this potential with its internal and external contributions as
3.8in terms of the *internal rate potential* per unit volume
3.9that contains the evolution of the energy storage function and the canonical dissipation potential function introduced in (3.4). The evolution of the macro- and micro-motion fields at a given state is determined by the *minimization principle*
3.10Note that the minimization structure of this variational principle is governed by the convexity of the dissipation potential function . Taking the variation of the rate potential (3.7) for admissible virtual rates of the concentration, which satisfies the homogeneous form of the Dirichlet boundary conditions, we obtain as the Euler equations of the variational principle (3.10) the *microforce balance* (3.1) along with the associated Neumann-type boundary conditions. They determine for a *given* chemical potential *μ* the concentration field *c*.

### (d) The mixed two-field saddle point principle

Canonical variational principle (3.10) can be combined with the generalized Legendre transformation (3.4) to an extended (mixed) variational principle that contains the chemical potential *μ* as a mixed variable. To this end, the *extended dissipation potential functional*
3.11is defined. Note that this functional now explicitly contains the constitutive dissipation potential function that governs the diffusion process. Based on the energy *E* introduced in (3.3), the above dual dissipation potential functional *D** and the external work functional , we define the *mixed two-field potential*
3.12at a given state {*c*}. Note that the only modification of (3.7) affects the dissipation, which is now expressed in terms of the rate of the concentration *and* the chemical potential *μ*. We rewrite this potential with its internal and external contributions as
3.13in terms of the *mixed internal rate potential* per unit volume
3.14Note that the potential *π** is *linear* with respect to the rates of the constitutive state variables and *concave* with respect to the chemical potential *μ*. The latter property is fully determined by the dissipation potential function , whose convexity is implied by the second axiom of thermodynamics. The evolution of the concentration as well as the chemical potential at a given state is determined by the *mixed two-field variational principle*
3.15which defines at the given state {*c*} and time *t* the rate of the concentration field along with the chemical potential *μ*. Taking the variation of potential (3.13), we find as the Euler equations of the variational principle (3.15) governing equations (3.1) and (3.2), i.e. the *microforce balance equation*
3.16and the *balance of species content*
3.17along with the Neumann boundary conditions. Note that the above two coupled field equations are *governed by variational derivatives of the rate potential π** defined in (3.14).

### (e) Time-discrete incremental variational principle

The rate-type variational structure outlined above is of great importance with regard to the time-discrete incremental setting of Cahn–Hilliard-type diffusion problems. To this end, consider a discrete time interval [*t*_{n},*t*] with step length *τ*:=*t*−*t*_{n}>0 and assume all field variables at time *t*_{n} to be *known*. The goal is then to determine the fields at the current discrete time *t* (variables without subscript) based on variational principles valid for the time increment under consideration. A mixed incremental variational principle of gradient-type dissipative solids is based on the modified functional
3.18in terms of the incremental energy functional *E*^{τ}, the incremental mixed dissipative work functional *D*^{*τ} and the load functional . It can be written in the form
3.19where we call *π*^{*τ} the *mixed incremental potential density* per unit volume. It is related to the rate potential *π** defined in (3.14) by the algorithm
3.20and governed by the energy function and the dissipation potential , respectively. Then, the finite-step-sized incremental variational principle
3.21determines the concentration field *c* as well as the chemical potential *μ at the current time t* as a *saddle point* of the incremental functional *Π*^{*τ}. The algorithm Algo is constructed such that the variation of the incremental potential (3.19) for admissible variations and yields as the Euler equations *consistent algorithmic counterparts* of the continuous equations (3.16) and (3.17). In what follows, we use a compact notation based on the arrays of (mixed) *primary* and *constitutive state variables*
3.22respectively. With these definitions at hand, consider as an example the simple algorithmic form of *mixed incremental potential* (3.20)
3.23based on a backward Euler approximation, which is a function of the constitutive state array defined in (3.22).^{2}

The subsequent finite-element treatment is then governed by the first and the second derivatives of this potential. This defines a *generalized stress array*
3.24and a symmetric *generalized tangent moduli array*
3.25where a decoupled energy function was assumed as considered in (4.1). These two arrays are the constitutive input for the subsequent discretization by the use of the finite-element method that uses Newton-type iterative solver, and make the notation extremely clear and compact. The symmetry of tangent moduli (3.25) for the coupled two-field representation is an important ingredient of the proposed incremental variational principle (3.21) based on algorithmic potential density (3.20). Note the indefinite structure of the moduli resulting from inf-sup problem (3.21). This intrinsic symmetry is not visible in previous treatments of computational splitting methods reported in Elliott *et al.* [20], Kuhl & Schmid [19] and Zhang *et al.* [13]. It is considered as an added value of the proposed variational theory.

### (f) Finite-element potential for Dirichlet problems

Based on time-discrete incremental variational principle (3.21), we now consider a spatial discretization of the coupled problem by a finite-element method. In order to keep the notation compact, we focus on *pure Dirichlet problems* with . The incorporation of Neumann conditions is straightforward. Let denote a finite-element triangulation of the solid domain . The index *h* indicates a typical mesh size based on *E*^{h} finite-element domains and *N*^{h} global nodal points. Associated with the triangulation , we write the finite-element interpolations of the macro- and micro-motion and the constitutive state vector in the compact form
3.26in terms of the *global nodal state vector* ** d**, which contains the primary variable array defined in (3.22) at a typical nodal point of the finite-element mesh. and are symbolic representations of the global matrices of shape functions for the coupled problem. In the model problem (4.7), we use the low-order C0 basis functions for

*both*the concentration

*c*and the chemical potential

*μ*. The spatial discretization of incremental two-field functional (3.19) takes the form 3.27for the Dirichlet problem with . Then, the finite-step-sized discrete minimization principle 3.28determines the nodal state vector of the finite-element mesh at the current time

*t*. The necessary condition of this discrete variational principle reads 3.29in terms of the

*generalized stress array*defined in (3.24) and provides a nonlinear algebraic system for the determination of the nodal state . For smooth problems, a standard Newton-type iteration of this algebraic system updates the state vector by the algorithm 3.30in terms of the monolithic tangent matrix , which is governed by the

*generalized tangent moduli array*defined in (3.25). The iteration update is performed until convergence is achieved. Note that the finite-element residual and tangent of the coupled problem are simply the first and the second derivatives of the potential

*Π*

^{*τh}defined in (3.27). Thereby, the

*symmetry of the tangent matrix*is a consequence of the proposed two-field variational principle in the extended (mixed) space of the concentration

*and*chemical potential fields

*c*and

*μ*, respectively. This property can be exploited by using solvers for the symmetric, indefinite problem (3.30)

_{1}. As pointed out by Toh & Phoon [21], sparse direct solvers taking advantage of symmetries of the underlying formulation are about 50% more efficient than non-symmetric sparse direct solvers.

## 4. Constitutive model and simulations of the two-field problem

### (a) Standard model problem of Cahn–Hilliard-type diffusion

The standard form of the Cahn–Hilliard equation is obtained by the specific choices of the two constitutive functions and . The *energy storage function* decomposes into a configurational contribution and an interface term
4.1The *configurational energy* is assumed to be identical to that considered in Cahn & Hilliard [1]
4.2The first term, governed by the parameter *A*, ensures that the concentration remains in the range *c*∈[0,1]. This is achieved by growth conditions for and , as depicted in figure 2*a*–*c*. The second term, governed by the parameter *B*, induces for *B*>2*A* a non-convexity, as shown in figure 2*d*–*f*, that favours a phase separation. Finally, consistent with the classical Cahn–Hilliard theory, the interface contribution to the free energy
4.3is assumed to be a quadratic function of the concentration gradient. The *dissipation potential* depends on the negative gradient of the chemical potential and is assumed to have the convex quadratic form
4.4at a given concentration *c*. The parameter *M*>0 is the species diffusivity. Table 1 summarizes the material parameters used in the subsequent simulations.

### (b) Numerical benchmark simulations of phase separation

In order to explore the modelling capabilities of the proposed symmetric two-field framework for the Cahn–Hilliard theory based on the time–space-discrete *two-field variational potential* (3.27), we carry out three simulations, which may serve as elementary benchmarks for future comparative studies. To this end, the formulations (3.29) and (3.30) were implemented in a straightforward manner into the finite-element research tool Feap [22], using its interface for parallel computing based on the Petsc library [23,24]. The Newton–Raphson updates (3.30) have been applied together with Jacobian-type preconditioning and the generalized minimal residual (Gmres) method developed by Saad & Schultz [25]. The chosen element type is Q1–Q1, see (4.7). A detailed comparison of alternative space–time discretizations with respect to accuracy and computing time is beyond the scope of the paper and a task for future investigations. We refer in this context to the recent works of Zhang *et al.* [13] on comparative studies with higher order finite elements for solutions of the standard CH equation and Gomez *et al.* [14] on its implementation by isogeometric analysis.

We consider a cubic domain with side length *L*=1/2 given by . On its surface , zero Neumann boundary conditions are prescribed for both the concentration *c* and the chemical potential *μ*, i.e.
4.5according to (2.11) and (2.12). The evolution of the concentration *c* starts from a given initial condition (2.13). Note that boundary condition (4.5)_{2} enforces a constant average concentration in and that hence the choice of *c*_{0}(** x**) determines as
4.6The domain is discretized by

*E*

^{h}=50

^{3}trilinear Q1–Q1 brick elements with parametric shape functions , yielding a mesh with

*N*

^{h}=51

^{3}nodes and a mesh size of

*h*=0.01. Using the

*N*

^{h}global nodal shape functions

*N*

^{I}(

**) that determine the finite-element arrays and in (3.26), both the concentration**

*x**c*and the chemical potential

*μ*can be approximated as 4.7Note from (4.7) that the initial conditions are thus fully determined by the nodal values as . As we want to

*investigate coarsening in mixtures with initially fine microstructure*, we start our simulations from random nodal initial conditions of the form 4.8where

*r*

^{I}are

*N*

^{h}(pseudo-)random numbers between 0 and 1 and where

*p*>0 is an exponent that (for sufficiently large

*N*

^{h}) determines the resulting value of as 4.9During the monolithic solution of the resulting symmetric linear system of equations, a constant time step of Δ

*t*=0.1 is applied.

The material parameters chosen for the benchmark simulations are summarized in table 1. Note that the choice of *B*=2.5 *A* leads to a double-well separation energy with minima at *c*=0.145 and *c*=0.855, see also figure 2*d*–*f*, such that for we *expect the formation and the consequent coarsening of precipitates.* Figures 3–5 show the predicted evolution of *c*(** x**,

*t*) from initial conditions, as specified by (4.8) with

*p*∈{2,3/2,1} and accordingly . As can be seen from the simulation results, the different initial conditions lead to different equilibrium states that are all at least close to local minimizers of the interface surface between the two phases characterized by species contents

*c*=0.145 (blue) and

*c*=0.855 (red), see online version. This is related to analytical considerations: For the limit of

*C*→0 (and hence zero interface thickness), it can be shown that solutions of stationary Cahn–Hilliard type problems indeed minimize the interface surface (e.g. [26]).

The results for *p*=2 and hence (see (4.8)) are depicted in figure 3. Here, the initial random distribution quickly evolves to a very fine distribution of (non-connected) precipitates with *c*=0.855 in a bulk material with *c*=0.145. The further evolution then consists of the typical precipitation coarsening. For the given finite value of *C*, the interface evolves towards an eighth sphere.

Figure 4 shows the evolution from random initial conditions with *p*=3/2, and hence . Again, the initial random distribution quickly evolves to a distribution of (non-connected) precipitates with *c*=0.855 in a bulk material with *c*=0.145. Owing to the higher average species content of , the precipitates have a much higher density than for . Nonetheless, we observe precipitation coarsening, which now leads to an equilibrium distribution in the form of a quarter cylinder.

Finally, figure 5 depicts the numerical experiment resulting from *p*=1, and hence . In contrast to the previous examples, the initial random distribution now evolves to an equal mixture of two strongly interwoven bulk phase networks of phases *c*=0.855 and *c*=0.145. Similar to the precipitates, this microstructure now evolves to a coarser and coarser structure. Here, the final state does indeed minimize the surface for the given volume fraction and leads to a flat interface separating the cube into two equal tetrahedral subdomains.

## 5. Conclusion

In summary, results of this work are: (i) the formulation of *generalized modelling structure* (2.9) for the Cahn–Hilliard-theory in terms of two constitutive functions and , (ii) the subsequent construction of the new *mixed variational principles* (3.15), (3.21) and (3.28), governing the continuous evolution, time-discrete update and time–space-discrete finite-element problems in terms of the *potential densities π** and *π*^{*τ} in (3.14) and (3.23), as well as (iii) their exploitation by a finite-element method with low-order shapes and *solver for symmetric, indefinite problems*. It was shown that the Cahn–Hilliard theory of diffusive phase separation is related to an intrinsic *mixed variational principle* (3.15) that determines the rate of concentration and the chemical potential. The principle characterizes a canonically compact model structure, where the two involved balance for the species content and microforce appear as the Euler equations of a variational statement. The existence of the variational principle underlines an *inherent symmetry* in the multi-field representation of the Cahn–Hilliard theory. This was exploited in the numerical implementation by the construction of time- and space-discrete *incremental potentials* (3.19) and (3.27), which fully determine the update problem of typical time-stepping procedures. The mixed variational principles provide the most fundamental approach to the finite-element solution of the Cahn–Hilliard-theory based on *low-order basis functions* and induce in a natural format the choice of *symmetric solvers* for Newton-type iterative updates, providing a speed-up when compared with non-symmetric implementations. This is a strong argument the first-choice use of the developed variational principles in the computational design. In this sense, the potentials developed are believed to be fundamental ingredients in the continuous and algorithmic setting of the Cahn–Hilliard theory.

## Funding statement

Support for this research was provided by the German Research Foundation (DFG) within the Cluster of Excellence Exc 310 *Simulation Technology* at the University of Stuttgart.

## Footnotes

↵1

*Definition of variational derivatives.*In order to get a compact notation, we use in the subsequent treatment the notation of variational derivatives*δ*_{y}*f*(,∇*y*):=∂*y*_{y}*f*(,∇*y*)−Div[∂*y*_{∇y}*f*(,∇*y*)].*y*↵2

*Variational consistency of incremental potential.*Note that formulation (3.23) of the incremental potential*π*^{*τ}for a dissipation function at*frozen statec*_{n}ensures the variational consistency with the continuous rate type potential*π** in (3.14).

- Received September 26, 2013.
- Accepted January 16, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.