(2017) Computational phase-field modeling. In: Encyclopedia of Computational Mechanics, Second Edition

Phase-ﬁeld modeling is emerging as a promising tool for the treatment of problems with interfaces. The classical description of interface problems requires the numerical solution of partial diﬀerential equations on moving domains in which the domain motions are also unknowns. The computational treatment of these problems requires moving meshes and is very diﬃcult when the moving domains undergo topological changes. Phase-ﬁeld modeling may be understood as a methodology to reformulate interface problems as equations posed on ﬁxed domains. In some cases, the phase-ﬁeld model may be shown to converge to the moving-boundary problem as a regularization parameter tends to zero, which shows the mathematical soundness of the approach. However, this is only part of the story because phase-ﬁeld models do not need to have a moving-boundary problem associated and can be rigorously derived from classical thermomechanics. In this context, the distinguishing feature is that constitutive models depend on the variational derivative of the free energy. In all, phase-ﬁeld models open the opportunity for the eﬃcient treatment of outstanding problems in computational mechanics, such as, the interaction of a large number of cracks in three dimensions, cavitation, ﬁlm and nucleate boiling, tumor growth or fully three-dimensional air-water ﬂows with surface tension. In addition, phase-ﬁeld models bring a new set of challenges for numerical discretization that will excite the computational mechanics community.


Phase-field modeling in computational mechanics
There are many areas of research in computational mechanics that involve moving boundary problems. Prime examples are fluid-structure interaction, fully three-dimensional air-water flows and crack propagation. The mathematical formulation of these problems often involves partialdifferential equations (PDEs) which are posed on moving domains and coupled to another set of PDEs through boundary conditions which hold on a moving and unknown interface. Although in some cases, such as for example, fluid-structure interaction, these problems have been attacked computationally with significant success (Bazilevs et al., 2013;Farhat, 2004), most of them remain as outstanding problems, especially in complicated three-dimensional geometries.
The phase-field method (Emmerich., 2003;Provatas and Elder, 2010;Chen, 2002;Anderson et al., 1998) may be understood as a mathematical theory to reformulate moving boundary problems as PDEs which hold on a known and fixed computational domain. The usual way to achieve this is to introduce a new field, namely, the phase field, which is defined on the entire domain and is a marker of the location of the different phases 1 . The phase field is governed by an equation formation of layers in the solution. Mathematically, the equations remain well posed, because they are regularized by higher-order operators, but usual time integration schemes may be inadequate. Sect. 4 treats in detail time integration algorithms for phase-field models, focusing both on accuracy and stability.
In all, we feel that the theoretical and computational developments of the last few years have elevated phase-field modeling to a point in which it can be used as a predictive model for a number of problems in computational mechanics which were very difficult to treat with conventional methodologies. Among these problems, we include the interaction of a large number of cracks in three-dimensional solids of complicated geometry (Borden et al., 2014), fully three-dimensional air-water flows with surface tension (Ceniceros et al., 2010), liquid-vapor phase transitions and cavitation (Liu et al., 2013), nucleate and film boiling (Liu et al., 2015), phase-change-driven implosion of thin structures (Bueno et al., 2014), cellular migration (Shao et al., 2012) and tumor growth (Hawkins-Daarud et al., 2012;Vilanova et al., 2013Vilanova et al., , 2014Vilanova et al., , 2012Xu et al., 2015;Lorenzo et al., 2015) and others Juanes, 2012, 2014;Gomez et al., 2013). However, even though the last few years witnessed very significant advances in the area, there are still many challenges ahead. Among these challenges, we identify multiphysics and coupled problems, nonlocal phase-field models, and an efficient treatment of randomness and uncertainty in phase-field theories.
We have organized the contents of this encyclopedic chapter as follows. Following a short subsection on notational conventions, we consider in Sect. 2 the heuristic derivation of phasefield models via diffusification of moving-boundary problems in the context of solidification. Then in Sect. 3, we consider phase-field models from a thermomechanical perspective, and we derive a variety of phase-field models that are encountered in computational mechanics. In Sect. 4 we consider various time-integration methods, and in Sect. 5 we summarize methods for the spatial discretization of phase-field models. Finally, Sect. 6 contains numerical examples.

Notational conventions
Let f : I ⊂ R → R be a real-valued function defined over a real, open interval I. The function f transforms x ∈ I into f (x). We denote the derivative of f as f . Alternatively, we will also use the notation df dx . When we work on several space dimensions, we will employ the notation Ω ⊂ R n d to denote an open subset of R n d , where n d is the number of spatial dimensions. In this paper, n d = 3 in most cases, but there are subsections in which n d = 2 or n d = 1. The closure of Ω is denoted by Ω and its boundary by ∂Ω. The set ∂Ω is supposed to have a well-defined unit outward normal. A spatial point x ∈ Ω is represented in component notation as x = (x 1 , x 2 , . . . , x n d ) T . Throughout the manuscript, the variable t ∈ R + ∪ {0} denotes time. We often use functions which depend upon space and time, for example φ : Ω × [0, T ] → R, where [0, T ] is the time interval of interest. The partial derivatives of φ are interchangeably denoted as We also make use of standard notation based on the operator ∇. For example, ∇φ denotes the spatial gradient of φ. We often resort to the use of material time derivatives as defined, e.g., in (Marsden and Hughes, 1994). For example,φ denotes the material derivative of φ which may be interpreted as a derivative with respect to time holding a material particle fixed. For a velocity field u, we haveφ = ∂ t φ+u·∇φ. Throughout the manuscript, we also utilize functions which are defined on sets of material points, rather than spatial points. Material points are also called particles and we denote them by X. For functions defined on sets of particles, material time derivatives are simply regular partial derivatives with respect to time. Sets of particles are denoted as Ω 0 or P 0 . The spatial position of a set of particles at time t is usually denoted P t . Volumetric integrals defined on sets of spatial points are usually denoted Ω dx, while for surface integrals we use the notation ∂Ω da. When the integrals are defined on sets of material points, we use the notation Ω0 dX and ∂Ω0 dA. The reader should also keep in mind that this paper describes a number of phase-field theories using the thermomechanics framework. Many of these theories involve, for example, Cauchy stress tensors, which are usually denoted T . The reader should not understand that T denotes exactly the same quantity throughout the manuscript. Rather, T denotes the Cauchy stress tensor of a particular mechanical theory ans its precise definition may vary from one section to another. This situation is not specific to T and happens also with other quantities such as, for example, the Helmholtz free energy Ψ, the mass flux h, the phase field φ, the reaction term R, the dissipation D and others. The precise definition is made clear in each subsection such that no confusion arises.
2. From classical moving-boundary problems to phase-field models via diffusification In this section we show how classical problems involving moving boundaries and interfaces can be regularized and written as PDEs posed on a fixed domain using the theory of phase-field modeling. The regularization replaces sharp interfaces by diffuse interfaces that are described with the help of the phase-field variable. As indicated before, we shall refer to this particular type of regularization as diffusification.
We consider the classical solidification theory for pure materials, also known as the generalized Stefan problem. The classical formulation consists of two PDEs which hold on different, but adjacent, time-dependent domains. The PDEs are coupled through conditions on the interface, whose location changes with time and is a priori unknown. We demonstrate heuristically how to diffusify this moving-boundary problem and obtain a phase-field model of solidification, which is a PDE system on a fixed domain.

Moving-boundary problem for solidification: The generalized Stefan problem
Let us consider a solid-liquid system that may undergo phase transformations. The system occupies the spatial domain Ω, which is fixed in time. The set Ω can be decomposed as Ω = Ω S ∪ Ω L , where Ω S and Ω L are the regions of Ω that host the solid and liquid phase, respectively. The sets Ω S and Ω L change with time, and their motions are unknowns of the problem. The solid-liquid interface, which is located where Ω S and Ω L meet is denoted Γ LS . Omitting the boundary conditions on ∂Ω for simplicity, the classical sharp-interface theory is described by the following equations The system composed by Eqs.
(2) and (3) correspond to the energy balance for a body with discontinuous material properties, (4) simply reflects that the temperature is continuous across the solid-liquid interface, and (5) is the so-called Gibbs-Thomson condition that determines the motion of the interface. In Eq. (2), e is the internal energy and q is the heat flux. Importantly, both e and q are discontinuous across the solid-fluid interface Γ LS due to the phase change. The internal energy is defined as where C v is the heat capacity, θ is the temperature, λ is the latent heat, and χ S is the characteristic function of Ω S , that is, a function defined on Ω that takes the value 0 in Ω L and 1 in Ω S . The heat flux is assumed to be governed by Fourier's law 4 , that is, q = −k∇θ, where the conductivity k equals the value k S in Ω S and the (possibly distinct) value k L in Ω L . In Eq. (3), n LS is the unit normal to Γ LS pointing towards the solid, V n is the normal velocity of the interface in the direction of n LS and the operator · denotes the jump across Γ LS such that q = q S − q L . In Eq. (5), ω is the coefficient of kinetic undercooling, σ is the surface tension of the liquid-solid interface, θ m is the melting temperature, and κ is the sum of principal curvatures of Γ LS (with the sign convention that κ is positive for liquid spherical bubbles).

The idea behind diffusification
In the above sharp-interface theory, we have three conditions at the interface Γ LS , which may be understood as one boundary condition for Ω S , one boundary condition for Ω L , and another one to determine the interface motion which is unknown a priori. Problem (2)-(5) is thus a movingboundary problem. Note that in the particular case ω = σ = 0, the Gibbs-Thomson condition becomes θ = θ m on Γ LS . In this particular case, the interface can be defined as the θ m level set of the temperature. However, in the general case in which ω and σ are non-zero, additional information is needed to locate the interface. In a phase-field formulation the additional information for the interface location is addressed by introducing a new variable, which is precisely the phase field, and which is nothing else but a marker of the phases' location. Note that this additional variable is not an artifact of the phase-field theory, but it is fundamentally necessary because the temperature is not enough to determine the position of the interface. Let us proceed to diffusify the sharp-interface theory and obtain the corresponding phase-field model. The key idea is to introduce a smooth field, namely, φ, which is defined on the entire domain Ω = Ω S ∪ Ω L and which is time-dependent, i.e., φ = φ(x, t). This field will be designed in such a way that it takes approximately the value −1 in the liquid phase, +1 in the solid phase, and it transitions quickly across the normal direction to the interface Γ LS . Let us use as a fundamental hypothesis a hyperbolic tangent profile, therefore where d t (x) denotes the signed distance from x to Γ LS (negative in the liquid and positive in the solid) and is a length scale related to the thickness of the diffuse interface. Note that Γ LS can be identified by the zero level set of φ. The subscript t in d t emphasizes that Γ LS changes with time, but will be removed from now on for notational simplicity. By using assumption (7), we will heuristically derive a phase-field model. Then, we will show that as → 0, the ansatz (7) is indeed recovered by the solution of the phase-field model, which gives coherence to the framework.

Derivation of a phase-field model for solidification
We proceed by suitably diffusifying the equations in (2), (3) and (5). We start with (5) and compute the geometric terms using the phase field φ. From basic geometry, it is known that |∇d| = 1, n LS = ∇d at Γ LS , and the curvature tensor κ satisfies κ = ∇∇d at Γ LS ; see, e.g., (Deckelnick et al., 2005). From Eq. (7), one can compute the spatial derivatives Using the hyperbolic-tangent identities f = 1 − f 2 , and f = −2f f , we then obtain the curvature tensor in terms of φ: The additive curvature, defined as κ = ∆d, can then be written as Next, using again the relation f = 1 − f 2 , it may be shown that Let us introduce at this point the function referred to as a double-well potential. The key feature of the function W is that it has two local minima which makes possible the coexistence of solid and liquid phases. Using the expression of W , we obtain We will now compute the normal velocity of Γ LS , which is given by the zero level set of φ. Let r LS (t) be a position vector that points to a material particle of the interface Γ LS , and define the function φ(t) = φ(r LS (t), t). The function φ is identically zero by definition of Γ LS . Therefore, its derivative is also identically zero. Using the chain rule, we obtain If we use the notation v LS = dr LS /dt, we can rewrite (15) as which is known as the level set equation (Fedkiw and Osher, 2003;Sethian, 1999). Using Eq. (8), and the definition of n LS , we conclude that V n is the velocity normal to Γ LS , that is, V n = v LS · n LS , thus introducing Eq. (17) into the level set equation, and using that f = 1 − f 2 = 1 − φ 2 , we get Substituting (14) and (18) into (5), we have now derived the phase-field equation: To obtain a regularized energy equation, we shall now diffusify (2)-(3). We consider smooth test functions v defined on the entire domain Ω and for which v = 0 on the external boundary ∂Ω. The test functions v are not necessarily zero on Γ LS , but [[v]] = 0 at Γ LS . We now claim that both (2) and (3) by Reynolds' Theorem, while the right-hand side equals Collecting the interior terms and the interface terms, and noting that [[e]] = λ, reveals the equivalence with (2) and (3). The final ingredients that we need to diffusify (20) are interpolation functions h(φ) and k(φ), which are required to satisfy h(±1) = 1 2 ± 1 2 , k(+1) = k S and k(−1) = k L . The simplest choices fulfilling these conditions are Using these interpolations, we approximate e and q by With these approximations in mind, e and q are smooth across Γ LS , and Eq. (20) is then simply equivalent to which, upon substitution, yields the diffusified energy equation The system of equations (19) and (28) completes the phase-field model, which represents an elementary phase-field solidification theory.

The diffuse-interface transition profile: Asymptotic analysis
Given the phase-field model of solidification, Eqs. (19) and (28), let us verify that it indeed gives rise to diffuse interfaces with the hyperbolic tangent profile (7) as → 0. We give a heuristic outline, but this procedure can be formalized using the theory of matched asymptotic expansions (Caginalp, 1989;Fife, 1988). Let us first verify separation into phases on short time-scales. Let τ = t/ 2 denote a rescaling of time. Then, Eq. (19) can be written as where we have used Landau notation on the right hand side. For small the coefficient in front of the left-hand side is significantly larger than the right-hand side so that which means that φ quickly evolves towards the stable roots φ = ±1 of −W (φ) = −φ 3 + φ. In other words, φ quickly separates into phases and generates a diffuse interface between them. Let us denote the zero level set of φ by Γ LS . It is expected that Γ LS converges to Γ LS as → 0. We next verify the profile of φ in a neighbourhood of Γ LS . We will zoom in to an O( )neighbourhood of Γ LS by considering a change of the independent variables. For simplicity, we consider a two-dimensional problem. We change from the variables (x 1 , x 2 ) to (z, s), where is a re-scaling of the signed distance to Γ LS (t) and s(x 1 , x 2 , t; ) is the arc length of Γ LS (t). Denoting φ(x 1 , x 2 , t) = Φ(z, s, t), then in a neighbourhood of Γ LS : where we have used the identity ∇s · ∇d = 0. Substitution into (19) yields: so that upon neglecting the right-hand side, one obtains, in a neighbourhood of Γ LS , This is a second-order ordinary differential equation subject to the boundary conditions Φ(z) → ±1 as z → ±∞. Its solution is the anticipated hyperbolic tangent profile: which was the hypothesis in the derivation of the phase-field model in Sect. 2.2.

General framework of thermomechanics and energy dissipation for phase-field models
This section illustrates the rigorous thermomechanical framework of phase-field theories. Our approach to thermomechanics follows closely the classical work of (Truesdell and Noll, 1965). Related work in the context of phase-field models can be found in (Gurtin, 1996) and (Oden et al., 2010). The main idea behind this approach is that theories of mechanics are fundamental balance laws of mass, linear momentum, angular momentum, and energy, while constitutive relations are suitably restricted by demanding that the theory be energy dissipative (or entropy productive). One refers to models derived this way as thermomechanically or thermodynamically consistent.
In the remainder of this section, we derive using the thermomechanics framework, the following models: • The classical Allen-Cahn and Cahn-Hilliard equations, which may be considered as the two canonical models of non-conserved and conserved phase dynamics, • The Navier-Stokes-Cahn-Hilliard equations which constitute a model for two-component immiscible flows with surface tension, • A simple phase-field model of brittle fracture, which is the diffuse version of Griffith's theory of fracture, • A phase-field model of tumor growth, which predicts the dynamics of cancerous tumors and nutrients in avascular tissue, • The thermal Navier-Stokes-Korteweg equations, which are a model for liquid-vapor transformations of a single-component fluid.
Another classical problem in which phase-field modeling is widely used is solidification. A thermomechanics framework may be used to derive solidification models (Wang et al., 1993;Penrose and Fife, 1990), but we do not cover this here because solidification was discussed in Sect. 2 using a different perspective. We remark that the diffusification procedure employed in Sect. 2 does not necessarily lead to thermomechanically-consistent models.

The idea behind thermomechanically-consistent phase-field modeling
Fundamental to phase-field thermomechanics is the dependence of the Helmholtz free energy Ψ (or entropy) on values and gradients of the phase field φ, and, additionally, the dependence of other constitutive relations on the values and/or gradients of the variational derivative of the total energy functional with respect to the phase. For example, the canonical free energy, which goes back to (Cahn and Hilliard, 1958), depends on the phase φ and its gradient ∇φ, and has the form for some function G and parameter . The corresponding energy functional is known as a Ginzburg-Landau functional : The variational derivative is then which is used to define the so-called chemical potential µ = δF δφ = G (φ) − 2 ∆φ. In phase-field models, constitutive relations are allowed to depend on µ and/or ∇µ. This will be worked out in detail for the above-mentioned phase-field models in the following sub-sections.
It is possible to interpret the energy functional F as a diffusification of the surface-energy functional where Γ int is an interior surface and ψ is a surface-energy coefficient. Note that F(Γ int ) associates energy to an interface which is proportional to the length of the interface. The simplest connection is obtained by considering F (φ ) and , with d Γint the distance function to Γ int . Note that this particular φ satisfies the required identity; see (12). It can then be shown that by employing the co-area formula with constant α q = ∞ −∞ q(z) dz. The co-area formula is valid for suitably decaying functions q(z), see (Du et al., 2005, Lemma 2.1). By taking the limit in (41), the parameter ψ can be suitably identified. To sum up, phase-field energies collapse to surface energies as the interface thickness approaches zero. 5

Allen-Cahn and Cahn-Hilliard
The two canonical phase-field theories are the Allen-Cahn equation and the Cahn-Hilliard equation (Emmerich., 2003;Provatas and Elder, 2010). They both emanate from the Ginzburg-Landau energy functional in Eq. (38), with G(φ) a double-well potential. The most usual approach would be to take G(φ) = W (φ) = 1 4 (1 − φ 2 ) 2 , but there are other possible choices, such as the logarithmic double well; see Eq. (194). Essentially the Allen-Cahn and Cahn-Hilliard equations are designed in such a way that they are gradient flows for the energy. This means that the energy is time decreasing along solutions to the equations. The Cahn-Hilliard equation incorporates the additional requirement that the mass of the components is conserved, while the Allen-Cahn equation allows for non-conservative dynamics.
The thermomechanical derivation of the Allen-Cahn equation and the Cahn-Hilliard equation starts by defining the constitutive class 6 for Ψ: 5 A rigorous connection between the functionals F and F can be made using the key mathematical theory of Γ-convergence (Ambrosio and Tortorelli, 1990;Braides, 2002;Dal Maso, 1993). 6 Constitutive class is a common nomenclature in classical mechanics; see, e.g., (Truesdell and Noll, 1965;Gurtin et al., 2009). In this particular case, it only means that the functional Ψ is allowed to depend upon φ and ∇φ. Fundamental principles could be used to restrict further the constitutive class Ψ = Ψ(φ, ∇φ). For example, frame invariance mandates that Ψ only depends on ∇φ through |∇φ|; see, for example, (Cahn and Hilliard, 1958).
Let V ⊂ Ω denote an arbitrary part of the spatial domain Ω. The postulated energy-dissipation property may be expressed as where the dissipation D(V) is ≥ 0 for all conceivable processes, while W(V) denotes the working, which is associated to external forces or energy supplies coming through the boundary of V. Basic manipulations on the left hand side of (44) lead to Let us now integrate by parts the last term of Eq. (45) after flipping the time and space derivatives. Doing so, we obtain, where µ is the chemical potential defined as the variational derivative: Now, we proceed to derive the Allen-Cahn and the Cahn-Hilliard equations.

Allen-Cahn equation
The Allen-Cahn equation can be derived by postulating the mass balance where R is designed to achieve free-energy dissipation and assumed to have the dependence: Note the non-standard dependence of R on the variational derivative of the free energy. This is one of the essential features of the phase-field method. Now, to find restrictions on R, we combine Eqs. (46)-(48), which leads to the expression Here, we identify D and W as D(V) = V µRdx and W(V) = ∂V ∂ ∇φ Ψ · n∂ t φda. The nonstandard work term is caused by internal forces ∂ ∇φ Ψ conjugate to changes in φ (Gurtin, 1996). It is clear that the constitutive choice R(φ, ∇φ, µ) = m(φ)µ with m(φ) ≥ 0 achieves the required property of free-energy dissipation. For the classical choice with W (φ) a double well potential, one has and we have thus derived the Allen-Cahn equation:

Cahn-Hilliard equation
We start with a statement of mass conservation 7 , namely where the constitutive class of h is given by Using again Eq. (46), it follows that Integrating by parts the first term on the right-hand side, we obtain where the first term on the right-hand side is identified as −D(V) and the second as W(V). Therefore, the constitutive choice guarantees free-energy dissipation and leads to the Cahn-Hilliard equation: for the above choice of Ψ. We note that µh in Eq. (56) can be interpreted as a free-energy flux.

Navier-Stokes-Cahn-Hilliard
The Navier-Stokes-Cahn-Hilliard equations represent a model for immiscible two-component fluid flow with surface tension. This theory has been widely studied from the physical, the mathematical and the computational points of view Jacqmin, 1999;Gal and Grasselli, 2010;Kay et al., 2008;Boyer et al., 2010;Gal and Grasselli, 2011;Lowengrub and Truskinovsky, 1998;Colli et al., 2012;Kim et al., 2004;Liu and Shen, 2003). To derive the theory, let us assume that we have two immiscible components with volume fractions ϕ 1 and ϕ 2 , respectively. The volume fractions verify the constraint ϕ 1 + ϕ 2 = 1 (59) The velocities of the components are denoted u 1 and u 2 , and their densities ρ 1 and ρ 2 . The density of the mixture, namely, ρ is defined as ρ = ϕ 1 ρ 1 + ϕ 2 ρ 2 . The mass-averaged velocity of the mixture is defined as For simplicity, both components of the mixture are supposed to be incompressible with constant density. In addition, we will assume ρ 1 = ρ 2 . Taken together, these assumptions imply that the mixture density ρ is also a constant equal to ρ 1 and ρ 2 . Therefore, without loss of generality, we can take In this scenario, mass conservation of the components implies where P t is a set of material particles in the current configuration. Using Reynold's theorem we can rewrite Eq. (62) as Summing Eqns. (63) over α, and using Eqns. (59), (60), and (61), we obtain the mixture mass balance div u = 0 Let us define at this point, the phase field φ = ϕ 1 −ϕ 2 , the so-called diffusion velocities w α = u α −u for α = 1, 2 and phase flux h = ϕ 1 w 1 − ϕ 2 w 2 . If we subtract Eqns. (63), we obtain the phase-field equation If we further assume a classical mixture (Romano and Marasco, 2010, Sec. 6.2), that is, a single momentum equation determines the mixture velocity u, then we have the following standard linear momentum balance equationu where T is the Cauchy stress tensor of the mixture which we assume to be symmetric 8 andu denotes the material derivative of u. The fundamental equations of the theory are (64), (65) and (66). Now, we will use the Coleman-Noll procedure (Coleman and Noll, 1963) to find constitutive equations for h and T which ensure that the total energy is dissipated along solutions of the theory. Under the assumption (61), the total energy (sum of free energy and kinetic energy) for an arbitrary part P t may be written as 9 The first term on the right hand side of Eq. (67) is assumed to pertain to the constitutive class where and ∇u T denotes the transpose of ∇u. Let us introduce the chemical potential µ, as before, where δE/δφ denotes the variational derivative of E with respect to φ. We assume that the mixture stress tensor T and and the mass flux h belong to the constitutive classes We postulate now the following energy dissipation law where W(P t ) denotes the working, which is associated to external forces or energy supplies coming through the boundary of P t and D(P t ) denotes the dissipation. As before, we will require that D(P t ) ≥ 0 for all conceivable processes.
Let us compute the left-hand side of (73). Using basic results of continuum mechanics and the constitutive class of Ψ, defined in (68), it follows that where we have also made use of Eq. (64). Using the relation in Eq. (74), and integrating by parts the term that results from ∇φ, we obtain Using Eq. (65) to computeφ, and integrating by parts the resulting term, we get (77) We can now proceed similarly with the kinetic energy. Using Reynold's theorem, Eq. (64) If we use now integration by parts on the term u · div T , we obtain If we sum (77) and (79), and we compare it with (73), we can identify the working W(P t ) and the dissipation D(P t ) as follows Now, D(P t ) ≥ 0 for all P t implies that the integrand has to be ≥ 0. There are several ways to design constitutive laws that satisfy that requirement, but in order to keep the process as simple as possible, we will impose that all the terms in the integrand of (81) are pointwisely positive or zero by themselves. For example, given that the dependence of Ψ onḊ has not been allowed in the constitutive class (68), the only way to make the term −∂ D Ψ :Ḋ pointwisely positive or zero is by taking ∂ D Ψ = 0 which implies Ψ = Ψ(φ, ∇φ). Considering this, and using the identity Without loss of generality, let us consider now the following splittings, where I represents the identity tensor and p is a scalar field that represents the mechanical pressure, which may be interpreted as a Lagrange multiplier of the incompressibility constraint. Using Eq.
(83), it follows that where we have used basic properties of the double-dot product 10 and Eq. (64). Let us mention now that, although ∂ ∇φ Ψ ⊗ ∇φ is not a symmetric tensor for arbitrary Ψ, physically relevant free energies always make the tensor symmetric. Indeed, the symmetry of ∂ ∇φ Ψ ⊗ ∇φ imposes a constraint on how the term ∇φ enters the free energy density Ψ, but this constraint only rules out terms which would be excluded anyway by using frame-invariance arguments, as shown in (Cahn and Hilliard, 1958). Therefore, we will use the fact that ∂ ∇φ Ψ ⊗ ∇φ is a symmetric tensor for physically relevant forms of Ψ, and thus, As a consequence, we can guarantee that the inequality (82) is satisfied by taking where m and ν are positive functions 11 that belong to the same constitutive classes as h and T , respectively. The function ν represents viscosity, and m is a mobility. From Eq. (88), we obtain Finally, considering a usual form for Ψ, that is, Ψ = γ W (φ) + 2 2 |∇φ| 2 , where γ is the surface tension, and W is a double-well potential, we obtain which finalizes the derivation of the theory. For completeness, we rewrite here the complete model Note that the final energy functional for a part P t is given by Eq. (67) It thus consists of two parts: a kinetic energy part and a mixing energy (or Ginzburg-Landau) part. As explained at the end of Sect. 3.1, this energy represents a diffusification of the classical energy of a two-phase moving boundary problem: consisting of kinetic energy and energy associated to surface tension (Gross and Reusken, 2011, Ch. 6), where Γ ∩ P t is the part of the two-phase interface inside of P t . The constant ψ can be show to be proportional to γ with a proportionality factor depending on the double well function W (φ) (Anderson et al., 1998).

Phase-field fracture theory
Classical approaches to treat fracture problems in computational mechanics include the virtual crack closure technique (Krueger, 2004) and the extended finite element method (Moës et al., 1999). These methods represent cracks as discrete discontinuities. From a computational point of view, this has been addressed by introducing discontinuity lines by way of remeshing or enriching the displacement field by means of the partition of unity method (Babuska and Melenk, 1997). However, this procedure has proven particularly difficult for problems in which a large number of cracks interact dynamically in complicated three-dimensional geometries. Phase-field methods, in contrast, permit solving the problem on a fixed mesh, which justifies their recent popularity in the computational mechanics community (Bourdin et al., 2008;Borden et al., 2012;Miehe et al., 2010;Abdollahi and Arias, 2011). The derivation of phase-field fracture theories can be interpreted as a diffuse version of the classical Griffith's theory of fracture.
To introduce the phase-field fracture theory we consider a reference configuration Ω 0 , which may be associated to the initial position of the solid. Points in Ω 0 are called material points or particles. The motion of the solid is defined by a time-dependent mapping ϕ which takes Ω 0 into the current configuration that we denote Ω. The displacement field is defined as d = ϕ − Id, where Id denotes the identity mapping. To develop the crack theory we also need to introduce the quantities With these definitions, we can introduce the free energy of the system, which is essentially obtained by regularizing the sharp-energy functional associated to Griffith's fracture theory: where Ψ e is the stored elastic energy density, ρ 0 is the density in the reference configuration, G c is the fracture toughness and Γ int is a two-dimensional manifold containing the cracked area. The regularized energy functional is obtained by replacing the surface integral on Γ int by a volume integral on Ω 0 which suitably converges to the surface integral as explained in Sect. 3.1. Therefore, the diffuse approximation to the energy functional may be written as where the variable c is a phase field which describes whether the material is fractured or not and which approaches the value 1 away from the crack and equals 0 inside the crack. The chief advantage of the phase-field approach is that the integrals in (99) are posed on a fixed and known domain, in contrast with the surface integral of Eq. (98). From a computational standpoint, this avoids remeshing and tracking the interface algorithmically. We next derive the phase-field fracture model. Consider a subset P 0 ⊂ Ω 0 of the reference configuration. The total energy associated to P 0 is where ρ 0 is the density in the reference configuration and The elastic energy Ψ e is assumed to depend on C and c, while Ψ c is defined as: where F (c) = 1 2 (c − 1) 2 is a single-well potential. For future reference, let us define We start the theory by stating classical balance laws, namely, where P is the first Piola-Kirchhoff stress tensor and b 0 represents a body force per unit volume 12 . We note that both equations are written in Lagrangian form.
Let us define at this point the following constitutive classes We postulate the energy dissipation inequality, where W(P 0 ) is the working and D(P 0 ) the dissipation. By requiring that the dissipation is nonnegative for all conceivable processes we will find constitutive equations for P and R. In fact, (109) Since we are working in Lagrangian coordinates now, (∇c) · = ∇ċ, and Eq. (109) can be rewritten as Introducing the second Piola-Kirchhoff stress tensor S defined as P = F S, and using the relation P :Ḟ = 1 2 S :Ċ, we obtain Therefore, we identify The choices where m is a positive function, guarantee that the dissipation is non-negative for all conceivable processes. The final form of the theory iṡ We note that in some applications, the termċ is neglected in the balance equations. Impressive simulations using this theory and slightly modified versions (including how to address irreversibility of the crack surface evolution) may be found in, for example, (Borden et al., 2012(Borden et al., , 2014Abdollahi and Arias, 2011).

Mechano-biological mixtures: Phase-field tumor growth theory
Phase-field modeling is also important in a number of emerging fields within computational mechanics. In this section, we consider the application of phase fields to mechano-biological continua. In particular, we consider avascular tumor growth, which is the growth of cancerous tumors without blood vessels. Phase-field tumor growth models have recently been proposed; see, e.g., (Cristini et al., 2009), (Lowengrub et al., 2010) and (Oden et al., 2015). Thermomechanical derivations of such models have been considered in (Wise et al., 2008), (Oden et al., 2010) and (Hawkins-Daarud et al., 2012). The starting point is the continuum theory of mixtures, as developed in (Truesdell, 1984, Lecture 5, 6) and (Bowen, 1976). Consider an isothermal mixture with N constituents (or species), α = 1, . . . , N . For tumors, such constituents include typically a tumor-cell phase, healthy-cell phase and a nutrient phase, but may also include extracellular water, extracellular matrix, other cell phases such as a necrotic-cell phase, and various chemicals. Each constituent can be assigned a mass fraction (or concentration) c α which is the mass of constituent α per unit mass. The mass balance for constituent α is: where ρ is the density of the mixture, u α the velocity of constituent α, and γ α the mass supply of α (owing to reactions). We assume saturation, i.e., α c α = 1. Invoking the axiom of mixture mass balance: the global mass balance is recovered by summing (118) over all α: where u = α c α u α is the (mass-averaged) mixture velocity. A diffusion-form of (118) can be obtained by employing the diffusion velocity w α = u α − u, yielding where h α = c α w α is the mass-fraction diffusion flux, and the dot( ) denotes a material derivative with respect to the mixture velocity u.
Instead of mass fractions c α , one can work with volume fractions ϕ α = c α ρ/ρ α where ρ α is the specific density (mass of α per unit volume of α). In this case the diffusion form is As in Sec. 3.3, continuing with the assumption of a classical mixture, we additionally have the global momentum equation where we have taken a fluid-mixture viewpoint, T being the Cauchy stress tensor and b a body force. A homogenized equation such as u = −K∇p is also possible, with K a tensor and p a pressure field. Mixtures with both fluid and solid components are considered in (Oden et al., 2010). We now develop the constitutive theory for a particular phase-field model. Let us assume for the sake of simplicity constant densities, all equal to 1: hence, by (120), In this case the constituent mass balance (121) simplifies tȯ For volume fractions, (122) simplifies to the same equation in (126) but with c α replaced by ϕ α . In fact, ϕ α = c α . As in the Navier-Stokes-Cahn-Hilliard and phase-field fracture theory, the current theory is based on the dissipation law: and requires D(P t ) ≥ 0 for all conceivable processes. Following the phase-field modeling paradigm, we choose a constitutive class with dependence on gradients: where D is defined in (69) and {(·) α } is a short-hand notation for (·) 1 , (·) 2 , . . . , (·) N . Furthermore, we consider the following classes for the stress, mass fluxes and reaction terms: with As before, by invoking the Coleman-Noll argument, suitable restrictions on the constitutive models are obtained. In particular, the stress T decomposes as where p is the pressure field which acts as a Lagrange multiplier for the incompressibility constraint (125), and T visc is the viscous part of the stress (typically chosen as 2νD, cf. (88)). Furthermore, it can be shown that the dissipation D(P t ) can be identified as The requirement D(P t ) ≥ 0 displays three dissipative processes: viscous dissipation, dissipation owing to reactions, and multi-component diffusion.
To illustrate an elementary thermodynamically-consistent tumor-growth model, we neglect the momentum equation (123) (i.e., u = 0) and consider four constituents: c 3 = nutrient-rich extracellular water , (137) c 4 = nutrient-poor extracellular water , We assume that the cell-phases occupy a fixed fractionc cell ∈ (0, 1), i.e., Thus, it is sufficient to work with two variables, e.g., the tumor phase and nutrient concentration: respectively. We select a free energy with W T (φ) a double well function with wells at φ = 0 and φ =c cell . This free energy corresponds to phase separation for the tumor phase φ and ordinary diffusion of nutrients σ. The relevant chemical potentials are To guarantee that D(P t ) ≥ 0, we choose for some growth function g(φ) ≥ 0, mobility m > 0, diffusivity D > 0, and reversibility parameter δ > 0 (which controls the amount of growth versus decline of the tumor phase). The final model is then given by the constituent mass balance equations in (121). With the above constitutive choices, and assuming that m and D are constants, they result in: which is a Cahn-Hilliard/Allen-Cahn equation coupled to an reaction-diffusion equation.

Thermally-coupled Navier-Stokes-Korteweg equations
The Navier-Stokes-Korteweg equations are a phase-field theory for single-component two-phase flows. Here, the two phases refer to liquid and gaseous states of the same fluid and are naturally located by the density field which identifies vapor with the low-density regions and liquid with the high-density areas. This model is somewhat different from those discussed previously because the phase-field plays a dual role. Indeed, the density field is used to locate the phases, but also represents an important physical quantity used in classical compressible gas dynamics. The Navier-Stokes-Korteweg equations allow for liquid-vapor phase transformations induced by temperature or pressure variations. Therefore, this model opens the possibility to solve vexing problems in computational mechanics such as for example, cavitation, film and nucleate boiling or the phasechange-driven implosion of thin structures (Bueno et al., 2014). The Navier-Stokes-Korteweg equations are rather new in the field, but there are interesting works of theoretical (Dunn and Serrin, 1986) and computational nature (Gomez et al., 2010b;Liu et al., 2013Liu et al., , 2015Tian et al., 2015;Giesselmann et al., 2014;Jamet et al., 2001). Our point of departure to derive the Navier-Stokes-Korteweg equations are classical balance laws for mass, linear momentum, angular momentum and energy, which may be written as followṡ ρ + ρ div u = 0 (149) Here, ρ is the density, u is the velocity, T represents the Cauchy stress tensor, b is a body force per unit mass, e is the total energy density, q is the heat flux, r is the radiation and π represents the energy variation due to interfaces. π is the only non-standard term and will be obtained directly from the Coleman-Noll procedure. To use the Coleman-Noll approach we will use the following entropy-production law, d dt where J (P t ) is an entropy flux associated to external heat sources and boundary terms and D(P t ) is required to be positive or zero. Using Reynold's theorem on the left hand side of Eq. (153), and the mass conservation equation, the entropy-production law can be rewritten as Using standard thermodynamic relationships, the energy density e is computed as the sum of the internal energy density ι and kinetic energy density, namely, Taking the material time derivative of Eq. (155), and multiplying with ρ, we obtain Let us use now the standard definition of the Helmholtz free energy Ψ, which is, Taking the time derivative of Eq. (157) and using Eq. (156) we can derive the relation Let us now use the expressions of ρė and ρu from Eqns. (152) and (150). We substitute them in (158), which results in We consider the following constitutive class for the Helmholtz free energy Ψ: Using the chain rule, it follows thaṫ Substituting the value ofΨ given in Eq. (161), and applying to ρ the relation given in Eq. (75), we obtain from (159) Without loss of generality, let us split T and ∇u as follows where and Note that T d and L d are traceless and W is skew-symmetric. This implies that where in the last identity we have used the mass conservation equation (149). In what follows, we will also use the relation Using Eqns. (167) and (168) in expression (162), we obtain Using now the identities and doing some manipulations, we obtain where the entropy flux turns out to be the classical term and D(P t ) is defined implicitly by Eqns. (172)-(173). Let us introduce at this point the following constitutive classes, Now, it may be easily proven that the following constitutive choices guarantee that D(P t ) remains nonnegative for all conceivable processes We notice that the term −ρ 2 µ in Eq. (178) may also be written as −ρ 2 ∂ ρ Ψ − div(∂ ∇ρ Ψ) and it is customary to identify −ρ 2 ∂ ρ Ψ with the thermodynamic pressure p. This concludes the derivation. If we want to use this theory as a predictive model, we need to give a precise definition of the Helmholtz free energy. The usual choice, associated to van der Waals fluids is where R is the specific gas constant, C v is the specific heat capacity, θ ref is a reference temperature, λ is associated to the surface tension of the liquid-vapor interface, and a and b are positive constants.

Energy-dissipative time-integration schemes
Time discretization methods for phase-field models typically have to deal with a nonlinear term W (φ) originating from the nonconvex (double-well in case of a binary system) free energy W (φ). The nonconvexity allows for 'backwards' diffusion, which is only mildly regularized by the higherorder operator coming from the interfacial energy, so that naive time-discrete schemes can lack stability. Because of this difficulty, the past decade has seen an increase in the development of novel time-discretization methods.
To illustrate the main ideas behind the existing methods, let us consider the time discretization of the elementary Cahn-Hilliard equation in split form 13 subject to natural boundary conditions on ∂Ω ∇φ · n = 0 and ∇µ · n = 0 (186) for all t ∈ (0, T ). As mentioned before, this model satisfies the free-energy dissipation where Let t n , n = 0, 1, . . . , denote discrete time-instances at which we wish to approximate the solution. Corresponding approximations are denoted φ n ≈ φ(·, t n ) and µ n ≈ µ(·, t n ) .
For simplicity we assume t n = nτ , where τ is a constant time-step size, but all schemes can be extended to the case of nonuniform time steps. Schemes are presented without discretization in space because we aim at understanding scheme properties which are independent of the spatial discretization. In the consideration of time-discrete schemes (Hairer et al., 1987;Hughes, 2012), we shall pay particular attention to (1) the order of accuracy of the scheme, (2) the stability of the algorithm, (3) the solvability of the time-discrete equations.
In the context of linear problems, stability is associated to the decay of the time-discrete solution with time, a feature also attained by the exact solution of stable linear problems. When the stability of the time discrete solution is achieved independently of τ , the time-integration scheme is said to be unconditionally stable. If stability holds under some constraint (e.g., on the time-step size) then the scheme is said to be conditionally stable. Explicit time integration algorithms which are consistent are always conditionally stable, while implicit schemes might be unconditionally stable.
For nonlinear problems the situation is much more complex than for linear equations. First, several notions of stability may be defined for different problems. Second, designing unconditionally stable schemes requires much more ingenuity than making the algorithm implicit. Because of the fundamental energetic structure behind phase-field models, see e.g. Eq. (187), natural notions of stability are those related to free-energy dissipation. In particular, a numerical scheme is said to be unconditionally energy-(or gradient-, or nonlinearly) stable if independent of τ . In other words, the scheme preserves the energy-dissipative property of the underlying model, in the sense that it dissipates energy at each time step. Analogously to the linear case, if (190) holds under some constraint, then the scheme is said to be conditionally energy stable. In analyzing time-integration schemes for phase-field theories, one needs to be specific about the smoothness of the particular chosen double well function W (φ). Two common assumptions are  (a, b), and there is a constant κ W > 0 such that (A2) W ∈ C 2,1 (R), i.e., W is globally Lipschitz continuous. This means that there is a constant L W > 0 such that Note that if W satisfies (A2) then it also satisfies (A1) (with a = −∞, b = ∞ and κ W = L W ), but not vice versa. In particular, the classical quartic potential satisfies (A1) with a = −∞ and b = ∞, but it does not satisfy (A2). The logarithmic potential satisfies (A1) with a = −1 and b = 1, and it also does not satisfy (A2). On the other hand, the truncated quartic potential 195) and the truncated logarithmic potential, which is defined e.g. in (Copetti and Elliott, 1992), satisfy both (A1) (with a = −∞ and b = ∞) and (A2).

First-order accurate schemes
We now consider first-order schemes for the Cahn-Hilliard equation. All of the schemes are implicit in some sense, since explicit schemes for fourth-order parabolic problems are infeasible.

Backward Euler
The backward Euler method applied to the Cahn-Hilliard equation (184)-(185) leads to the system This method is nonlinearly implicit, because it requires the solution of a nonlinear system for the pair (φ n+1 , µ n+1 ). As anticipated before, implicit schemes which are unconditionally linearly stable, such as the Backward Euler method, are generally conditionally stable for nonlinear problems. To see this, consider the energy difference Given two real numbers u and v that enclose ξ, we can use a Taylor's series expansion to show that 199) where the inequality holds due to assumption (A1), which is valid for all the potentials W of interest. One then obtains where · denotes the L 2 (Ω)-norm. To derive Eq. (200) we have made use of Eqns. (196)-(197) and the natural boundary conditions defined in (186). Next we bound φ n+1 − φ n 2 using (196): where the last inequality holds for for any δ > 0. Finally choosing δ = 2 Therefore, if then (190) holds. In other words, the backward Euler scheme is conditionally energy stable. Furthermore, it can be shown, that under the same time-step constraint, the system (196)-(197) is uniquely solvable, i.e., the nonlinear system (196)-(197) has a unique solution for (φ n+1 , µ n+1 ) if (202) holds. The proof of this can be found in, e.g., (Elliott, 1989) and the idea behind the proof applies to other schemes. To proof the existence of a solution, one shows that the system (196)-(197) is the necessary condition (or Euler-Lagrange equation) corresponding to a minimization problem for a convex functional. Next, to proof the uniqueness of a solution, one carries out steps similar to proving the above energy stability. Note that since is generally very small, (202) implies a severe constraint on the allowed time-step size.
4.1.2. First-order semi-implicit method A popular scheme (Provatas and Elder, 2010) is the following first-order semi-implicit (or implicit/explicit) method: Because it treats W explicitly, it is a linear (or linearly-implicit) method requiring the solution of a linear system at each time step. The system can be written abstractly as where B is the differential operator defined by Here, Id denotes the identity operator. The linear system of differential equations (205) has a unique solution, independent of τ . This follows from the coercivity estimate as well as mass conservation Ω (φ n+1 − φ n )dx = 0 and the condition Ω µ n+1 dx = Ω W (φ n )dx. However, as may be expected, the method is only conditionally energy stable. In fact, we can only show conditional stability if W satisfies (A2). Let us assume that ξ is an undetermined point of the interval (0, 1). Then, where φ n+ξ ≈ φ(·, t n + ξτ ). The inequality (208) is identical to (200) replacing κ W with L W . Therefore, following the same steps as in Sect. 4.1.1, leads to the similar constraint τ < 8 2 /L 2 W for energy stability.

Convex-splitting method
The stability issues for the backward Euler and semi-implicit method originate from the nonconvexity of W (φ). A groundbreaking idea, which goes back to (Elliott and Stuart, 1993) and was popularized by (Eyre, 1998), is to split W into a convex part and concave part, i.e., The, we treat W + implicitly and W − explicitly, as Using the Taylor's formulae where ξ, ζ ∈ (0, 1), the energy difference is now which shows the unconditional stability of the method. If W satisfies (A1), then the above convex splitting is nonunique, but it is always possible, e.g., defining and subsequently setting W + (φ) = W (φ) − W − (φ). Unique solvability of the system (210) follows by equivalence with a strictly-convex minimization problem (cf. (Wise et al., 2009, Theorem 3.4)).
Hence, the convex-splitting scheme is both unconditionally stable and unconditionally solvable. If, in addition, (A2) applies, then a splitting is possible with W + (φ) a quadratic polynomial, i.e., defining and subsequently setting W − (φ) = W (φ) − W + (φ). In this case, the scheme becomes even linear. See also (He et al., 2007;Shen and Yang, 2010), where splitting is seen as stabilization of the semi-implicit scheme.

Second-order accurate schemes
The development of novel second-order time-accurate schemes with attractive stability and solvability properties is still ongoing. The situation is much more difficult than for first-order methods. Most investigated schemes are based on modifications of the Crank-Nicolson method. We start by analyzing the (poor) properties of the Crank-Nicolson method and then we show how the method can be improved.

Crank-Nicolson The Crank-Nicolson method is defined by:
It requires the solution of a nonlinear system at each step. Using the trapezoidal quadrature rule we have where ξ is an unknown point that lays between u and v. Then, it can be shown that The sign of the last term can not be controlled, therefore the method is generally not unconditionally energy stable. Furthermore, since this method corresponds to a nonconvex nonlinear system, it can be shown (by mimicking the proof for the Backward Euler method in (Elliott, 1989)) that its solvability suffers from a similar time-step constraint, τ ≤ C 2 /κ 2 W , as the backward Euler scheme.

Second-order semi-implicit method
The following second-order semi-implicit method (Guillén-González and Tierra, 2013) shares the same stability and solvability properties as the Crank-Nicolson method: However, it is a linear method corresponding to the system where B n is the differential operator defined by Because of nonconvexity of the potential, W can have an arbitrary sign. Therefore, the method is conditionally uniquely solvable under a time-step constraint τ < C 2 /κ 2 W (Guillén-González and Tierra, 2013).

Secant method
The secant method, which goes back to (Du and Nicolaides, 1991), is designed to exactly mimic the energy dissipation. The method is defined as, or, equivalently, which explains the name of the method. Since it follows that the method is unconditionally energy stable, as shown by the relation Because of nonconvexity of W however, the nonlinear system (225)-(226) is only conditionally solvable (Elliott, 1989). The operator DW is usually referred to as discrete variational derivative. The use of discrete variational derivatives is ubiquitous in the design of exact energy preserving schemes and symplectic methods; see e.g., (Furihata and Matsuo, 2011;Simo et al., 1992;Romero, 2009).

Secant variants
Two schemes which avoid constructing the secant, while still maintaining stability, are the implicit Taylor method by (Kim et al., 2004), given by (225)- (226) with and the method by (Gomez and Hughes, 2011), given by (225)- (226) with These two methods result in The last term is negative or zero if W ≥ 0. The inequality W ≥ 0 is true for all potential functions in (193)-(195). In the case in which the potential W does not satisfy the condition W ≥ 0, (Gomez and Hughes, 2011) propose a splitting of W which achieves second-order accuracy and unconditional stability. A different method based on a second-order perturbation of the midpoint rule was proposed in (Liu et al., 2013). The method in (Liu et al., 2013) achieves essentially the same properties.
As the methods discussed in Sect. 4.2.4 correspond to nonconvex nonlinear systems their solvability is expected to be conditional.

Second-order convex splitting
In an effort to have simultaneously unconditional stability and solvability, a second-order time-accurate convex-splitting scheme was proposed by (Hu et al., 2009). The method is given by where which corresponds to a secant treatment of the convex part and a two-step second-order backward difference treatment of the concave part. This method is unconditionally solvable as it is equivalent to a convex minimization problem at each time step. However, the method is not unconditionally stable as defined in (190). A less restrictive statement referred to as weak energy stability (Hu et al., 2009) may be proven if W − is a quadratic polynomial, which applies to many potential functions. An alternative is the multistep scheme proposed in (Guillén-González and Tierra, 2013), which is however restricted to the quartic potential. This scheme is linear, unconditionally energy stable for a modified energy statement, and unconditionally uniquely solvable.

Stabilization
Most second-order linear schemes suffer from conditional stability and/or conditional solvability. As proposed by (Wu et al., 2014) these issues can be stabilized if W satisfies (A2). For example, the stabilized semi-implicit scheme is where β is the stabilization parameter. Such terms are also referred to as artificial viscosity, and are useful in many applications; see e.g., (Labovsky et al., 2009;Jansen et al., 2000). See also (Gomez and Hughes, 2011) for a stabilization of the extended Crank-Nicolson method. Using a Taylor's formula, it can be shown that Then, assuming that W satisfies (A2), (by Young's ineq. with 0 < δ < L W /2) Therefore, for β > L 2 W /4, one has unconditional energy stability. We note that unconditional solvability follows by mimicking the proof for the Backward Euler method in (Elliott, 1989).

Discussion
Based on the Cahn-Hilliard equation, we have tried to illustrate that the design of stable and solvable time-integration schemes for phase-field models is challenging. The key idea for first-order schemes has been convex-splitting, however, this idea is not sufficient for second-order schemes, and one has to resort to particular double-well potentials, modified energy statements and/or stabilization. There are also other popular methods that can be used, such as the generalized-α algorithm (Jansen et al., 2000), for which little is known about their stability properties in case of phase-field models.
For phase-field models that are coupled to additional equations (e.g., the momentum equation), difficulties arise in maintaining energy stability and solvability. The idea of stable decoupled schemes is also challenging, but it is attractive to achieve modularity of computer codes; see for example (Shen and Yang, 2015) in the context of the Cahn-Hilliard-Navier-Stokes model.

Space discretization
The main challenges associated to the spatial discretization of phase-field theories are the approximation of higher-order partial-differential operators, inherent to diffuse-interface models and the presence of internal layers in the solution, which clearly calls for adaptive meshing. Although these are long-standing problems which remain open to some extent, here we show important developments accomplished in the last few years that have elevated phase-fields to a point in which they can be used for practical engineering problems.

Discretization of higher-order partial-differential operators
Most phase-field theories involve higher-order spatial partial-differential operators which give rise to diffuse interfaces in the solution. For example, the Cahn-Hilliard equation is a fourth order PDE in space. The Navier-Stokes-Korteweg equations include third-order spatial derivatives. The fracture theory presented in Sect. 3.4 is a second-order PDE in space, but a higher-order approximation which involves fourth-order derivatives has been recently proposed. Data indicate that the higherorder theory outperforms the model presented in Sect. 3.4 (Borden et al., 2014). Therefore, it seems clear that an important requirement of a numerical algorithm for phase fields is the possibility to discretize higher-order derivatives effectively.
The higher-order operator problem is an old one in computational mechanics. As a matter of fact, classical theories for thin shells and plates pertaining to the Kirchhoff-Love type also involve fourth-order operators. Numerical solutions to higher-order problems in simple geometries can be easily obtained using classical algorithms, such as, for example, the finite difference method, or pseudo-spectral collocation. The problem shows when one needs to solve a higher-order PDE on a complicated three-dimensional geometry. Complex geometries typically call for the use of the finite element method. However, solving fourth-order PDEs directly with finite elements requires globally C 1 -continuous basis functions which are very difficult or even impossible to generate in complex three-dimensional geometries.
Perhaps, the most obvious solution to the higher-order problem is to split the fourth-order equation into a couple of second-order PDEs, but this would lead to a twofold increase in the global number of degrees of freedom. As a consequence, different procedures have been investigated that aim at keeping the computational cost as low as possible. Relevant examples are the continuous/discontinuous Galerkin method (Wells et al., 2006), finite volume methods (Cueto-Felgueroso and Peraire, 2008) or meshless methods (Rajagopal et al., 2010;Zhou and Li, 2006). Here, we will focus on a novel technology, namely isogeometric analysis, which permits simple discretizations of higher-order operators on non-trivial geometries. The key idea is that isogeometric analysis permits generating globally C p−1 -continuous basis functions in physical space by mapping order-p B-Splines or NURBS (Non-uniform rational B-Splines) from a reference configuration. For the reader not familiar with B-Splines and NURBS we recommend the references (Hughes et al., 2005;Cottrell et al., 2009;Aaa, 201xa,x), which are intended for scholars with a background in finite elements. In Sect. 6 we present more details on the spatial discretization of various phase-field models.

Adaptive mesh and time-step refinement
Solutions to phase-field models inherently exhibit a large range of spatial and temporal scales, therefore adaptive mesh refinement and adaptive time-step selection are indispensable for the efficient numerical discretization of phase-field models. Clearly, the phase-field itself is constant everywhere throughout the domain except at interfaces, where it rapidly, but smoothly, changes from one value to the other. This smooth transition happens on a spatial scale of order , and therefore, in principle, its resolution requires a local mesh-size of that same order. In time, rapid phenomena also occur, in particular when topological changes take place. For example, in case of Cahn-Hilliard phenomena, this happens when mixed phases separate, single phases break up into several bubbles, or when separate bubbles coalesce.
Adaptive mesh refinement crucially relies on refinement indicators, which mark elements for refinement that are deemed important for improving the accuracy of a coarse-mesh approximation. The indicators can be based on heuristic criteria or a posteriori error estimates. An example of a heuristic criteria is one for which those elements are marked if the local phase-field approximation deviates from pure phase-field values (Provatas et al., 2005), or if it has a large gradient (Rosam et al., 2007). Typically these indicators result in reasonable refinements, but it may not work well on very coarse meshes with very poorly resolved solutions.
Alternatively, a posteriori error estimates can be employed in a method-of-lines approach. In such an implementation, the estimates are derived for the fully discrete system for advancing a time step, but they estimate errors only due to the spatial part. Any a posteriori error estimator that is suitable for the spatial operator can be employed, e.g., the Zienkiewicz-Zhu gradient-recovery estimator (Provatas et al., 1999), or an explicit residual-based estimator (Stogner et al., 2008).
Adaptive time-step selection is important for temporal accuracy. In a method-of-lines approach one resorts to classical techniques for adaptive time-stepping of ordinary differential equations (Söderlind, 2002;Deuflhard and Bornemann, 2002). An elementary approach computes the approximation for the next time step using a low-order and a higher-order time-stepping scheme. Their difference can then be employed as an estimate for the temporal error and as a guide to decrease or increase the time-step size. For phase-field models such techniques have been employed in, e.g., (Ceniceros and Roma, 2007;Cueto-Felgueroso and Peraire, 2008;Gomez et al., 2008;Wodo and Ganapathysubramanian, 2011).
Finally, space/time a posteriori error estimates can be employed, which can estimate errors due to space and/or time discretization. For residual-based estimates, the difficulty is the derivation of sharp bounds by eliminating the exponential dependence on the interface thickness; see, e.g., (Kessler et al., 2004;Feng and Wu, 2005;Bartels, 2005) for Allen-Cahn estimates and (Feng and Wu, 2008;Bartels and Müller, 2010) for Cahn-Hilliard estimates. Alternatively, estimates based on the computation of a dual problem have recently been developed (Simsek et al., 2015) in the spirit of duality-based estimates for nonlinear parabolic problems (Eriksson et al., 2004). It is also possible to employ goal-oriented error estimates to control the accuracy in output quantities of interest (van der Zee et al., 2011;Mahnken, 2013). These are based on general frameworks as described in, e.g., (Becker and Rannacher, 2001;Stein and Rüter, 2004).

Applications and numerical examples
Here we illustrate how the space discretization of fourth-order problems can be accomplished using smooth basis functions by way of isogeometric analysis (Sect. 6.1) or resorting to classical mixed methods (Sect. 6.2). In Sect. 6.1 we focus on the classical Cahn-Hilliard equation, while in Sect. 6.2 we illustrate the ideas with the tumor-growth model presented before. Finally, Sect. 6.3 presents isogeometric simulations of the Navier-Stokes-Korteweg equations.

Cahn-Hilliard equation
We start by deriving a weak form of the Cahn-Hilliard equation. For simplicity, we assume periodic boundary conditions. Thus, given the Cahn-Hilliard equation, recall from (58), let us consider a smooth function w ∈ V, where V = H 2 (Ω) is the Sobolev space of square integrable functions with square integrable first and second derivatives in Ω. We multiply Eq. (242) with w, integrate over the domain Ω and integrate by parts twice, which leads to the following variational problem: Find φ ∈ V such that for all w ∈ V The key point is that by way of isogeometric analysis we can construct a discrete space V h ⊂ V using B-Splines or NURBS of order p ≥ 2. Therefore, using the Galerkin method, the semi-discrete problem may be stated as: Eq. (244) represents a semidiscrete form of the Cahn-Hilliard equation. Time integration schemes such as discussed in Sect. 4 can be utilized to advance in time.
Our numerical examples have been computed using the generalized-α algorithm (Jansen et al., 2000). We omit here the implementation details which may be found in (Gomez et al., 2008). We use adaptive time stepping because the equation evolves on very different time scales during the simulation (Gomez et al., 2008;Gomez and Hughes, 2011). To completely define the problem we need to choose specific forms of the double-well potential W and the mobility M . We take the physically-relevant logarithmic double well potential, see (194), since the quartic polynomial is a mere qualitative approximation of the logarithmic double well. The parameter θ in (194) represents the quench temperature. In our computations we will take θ = 3/2. Similarly, the mobility is often assumed constant, although the so-called degenerate mobility leads to better agreement with experiments and is our choice for the numerical examples. The key feature of the degenerate mobility is that pure phases have vanishing mobility and the motion is restricted to the (diffuse) interfaces as it happens in so-called Ostwald ripening.
Using dimensional analysis, it may be shown that the solution to the Cahn-Hilliard equation does not depend on all its parameters, but only on the dimensionless number (Gomez et al., 2008) Figure 1. Two-dimensional simulation of the Cahn-Hilliard equation on the computational domain Ω = [0, 1] 2 , with periodic boundary conditions. We started the simulation with a randomly perturbed homogeneous solution (φ = 0.5), which spontaneously separated into two phases (blue and red colors). As time evolves (left to right and top to bottom), the system undergoes coarsening, which manifests itself by larger-scale morphological features. The mesh is uniform and composed of 256 2 C 1 -quadratic elements. We used the parameter C = 10 4 .
where L 0 is a length scale of the problem. Here, we will assume that L 0 is a measure of the computational domain size. For example, for a cubic domain, L 0 is just the length of the cube's edge. Without loss of generality, we will also assume L 0 = 1. The usual setup for numerical simulations of the Cahn-Hilliard equation is as follows: a randomly perturbed homogeneous solution is allowed to evolve in a sufficiently large computational domain. Periodic boundary conditions are often used in an attempt to isolate the phase-separation physics from boundary effects. It is possible to show that solutions of the type φ(x, ·) = φ + η(x), where η(x) is a small perturbation are linearly unstable when φ is in the so-called spinodal region, that is, when W (φ) < 0 . Thus, if we start the simulation with the initial condition φ(x, 0) = φ+η(x), the perturbation η(x) will grow in time, leading to spontaneous phase separation. Phase separation manifests itself as patches where the solution is almost constant and takes the values φ ≈ ±1. The regions where φ ≈ −1 and φ ≈ +1 are connected through a transition layer.
To illustrate this process, we computed a two-dimensional solution using the parameters φ = 0.5 and C = 10 4 . The top left snapshot of Fig. 1 shows the system early after phase separation. The blue areas represent patches where φ ≈ −1, while the red zones correspond to φ ≈ +1. After the system is fully separated, the process of coarsening starts. During the coarsening process, small bubbles merge together giving rise to larger-scale features (see the remaining snapshots of Fig. 1; time increases from left to right and from top to bottom). Fig. 2 shows a similar computation in three dimensions. The parameters are φ = 0.26 and C = 600. The time evolution may be observed in the plot. As time evolves, the mixture coarsens until it is fully separated. We start with a randomly perturbed homogeneous solution (φ = 0.26). The physics is analogous to that of the two-dimensional example . The mixture separates and then coarsens as time evolves (left to right and top to bottom). The mesh is uniform and composed of 128 3 C 1 -quadratic elements.

Phase-field tumor growth theory
Next, we consider a numerical illustration for the tumor-growth model derived in Sect. 3.5, subject to the boundary conditions and initial conditions As an alternative to the spline-based discretization in the previous section, we keep the chemical potential µ 1 as an unknown and consider a mixed formulation suitable for classical C 0 -continuous finite element spaces V h . After semi-discretization in space the problem is then to find for all (w h , η h , ζ h ) ∈ V h × V h × V h . Note that classical C 0 -continuous shape functions are possible since ∇(·) is the highest-order derivative.
To illustrate a fully discrete method, let us employ a simple first-order time-accurate convexsplitting scheme, extending the one for the Cahn-Hilliard equation described in Sect. 4.1.3: Note that we used first-order extrapolation on the nonlinearity caused by g(φ h ). This scheme can be shown to be unconditionally energy stable so that large time-steps can be taken (Wu et al., 2014).
As an example of a numerical simulation which can be performed on a laptop computer, we consider the domain Ω = [−1, 1] 3 , and we set the parameters m = D = 1,c cell = 1 2 , δ = 0.01, = 0.1, and For the convex spliting in the scheme, we decompose W T (φ) by choosing W T + (φ) = 4φ and setting W T − (φ) = W T (φ) − W T + (φ). The nutrient boundary-condition value is σ D = 1 −c cell = 1 2 , the initial nutrient distribution is σ 0 = 1 −c cell = 1 2 , and the initial condition for the tumor φ corresponds to three spheres: with radii R 1 = 1/6, R 2 = 1/4 and R 3 = 1/6, and center coordinates x √ 2, 1 6 ). We choose V h to correspond to linear finite elements on a coarse uniform mesh with 40 3 elements, and set the time-step size to a large value: τ = 0.01. Figure 3 displays several snapshots of the simulation.

Isothermal Navier-Stokes-Korteweg equations
Here, we show how the Navier-Stokes-Korteweg equations can be discretized using globally smooth splines. For simplicity, we consider only the isothermal version of the equations, although similar ideas could be used to discretize the thermally-coupled theory (Liu et al., 2015). Under the assumption of constant temperature, the Navier-Stokes-Korteweg equations reduce to ∂ρ ∂t + div(ρu) = 0 (258) ∂(ρu) ∂t + div(ρu ⊗ u + pI) − div σ = ρb (259) where σ = ν (∇u + ∇ T u) − 2 3 div u I +λ ρ∆ρ + 1 2 |∇ρ| 2 I −λ∇ρ ⊗ ∇ρ and Eq. (261) is referred to as the van der Waals equation, and acts as an equation of state in the theory. In the van der Waals equation, R is simply the specific gas constant, while a and b are positive parameters that may be found in scientific databases, such as for example, that of the National Institute of Standards and Technology. The temperature θ plays the role of a parameter in the isothermal theory. For notational simplicity, we define T visc as the viscous stress, that is, Using the notation (262), it may be shown that div σ = div T visc +λρ∇(∆ρ) Substituting Eq. (263) in (259), we get the non-conservative form of the Navier-Stokes-Korteweg equations, which may be expressed as ∂ρ ∂t + div(ρu) = 0 (264) ∂(ρu) ∂t + div(ρu ⊗ u + pI) − div T visc −λρ∇(∆ρ) = ρb Eqs. (264) and (265) constitute the starting point of our weak formulation, which may be derived by multiplying them with smooth functions q and w, integrating over the computational domain, and then integrating by parts. Following the ideas presented in the previous subsections, we approximate relevant functional spaces by their discrete counterparts. The discrete spaces are subsets of H 2 (Ω) due to the global smoothness of B-Splines and NURBS (for p ≥ 2) and this makes the integrals in the weak form computable. Proceeding this way, we can recast the problem as: find ρ h ∈ V h and u h ∈ V h n d such that for all q h ∈ V h and w h ∈ V h n d

=
where we have assumed periodic boundary conditions for simplicity. We integrate in time Eqs.
(266)-(267) using the generalized-α method. We use this algorithm to perform two-and threedimensional computations. By way of dimensional analysis, it may be shown that the solution to the unforced Navier-Stokes-Korteweg equations depends only on the dimensionless groups where Re is the Reynolds number and Ca is the capillary number. L 0 = 1 represents again the length of the computational domain. We begin by simulating the coalescence of two vapor bubbles in the computational domain Ω = [0, 1] 2 . The centers of the bubbles are located at points x (1) = (0.40, 0.50) and x (2) = (0.78, 0.50). Their radii are R 1 = 0.25 and R 2 = 0.10, respectively. Inside the bubbles, we set a dimensionless density ρ/b = 0.1, while outside we take ρ/b = 0.6. These are approximate values to the so-called Maxwell states, which correspond to mechanical and chemical equilibrium values. The density profiles are regularized using hyperbolic tangent profiles. The initial velocity is zero. The dimensionless groups are given by Re = 512 and Ca = 1/256. The mesh is uniform and composed of 256 2 C 1 -quadratic elements. The time evolution of the system may be observed in Fig. 4. The background color scale represents density, while solid lines represent streamlines. The color scale superimposed to the streamlines corresponds to the velocity magnitude. The simulation shows how mechanical forces drive the two bubbles together, creating a velocity field. After the two bubbles are connected, capillary forces tend to reduce the interfacial length, eventually leading to a stationary state with a circular steady bubble (not shown). Fig. 5 shows the three-dimensional analogue of the last example. The computational domain is Ω = [0, 1] 3 , and the mesh is uniform and composed of 128 3 C 1 -quadratic examples. The centers of the bubbles are located at the points x (1) = (0.40, 0.50, 0.60) and x (2) = (0.75, 0.50, 0.50). Their radii are R 1 = 0.25 and R 2 = 0.10, respectively. The initial conditions are the same as in the last example, while the dimensionless groups are Re = 256 and Ca = 1/128. Fig. 5 presents the time evolution of the system. The bubbles are represented by density isosurfaces, which have a velocity magnitude color scale superimposed. The system evolves towards a steady spherical bubble (not shown). Figure 5. Three-dimensional vapor bubbles right before and after coalescence due to surface tension. The bubbles are represented by density isosurfaces, which are superimposed with a velocity magnitude color scale. The mesh is uniform and composed of 128 3 C 1 -quadratic elements.