Elsevier

Journal of Computational Physics

Volume 305, 15 January 2016, Pages 360-371
Journal of Computational Physics

Isogeometric analysis of the Cahn–Hilliard equation – a convergence study

https://doi.org/10.1016/j.jcp.2015.10.047Get rights and content

Abstract

Herein, we present a numerical convergence study of the Cahn–Hilliard phase-field model within an isogeometric finite element analysis framework. Using a manufactured solution, a mixed formulation of the Cahn–Hilliard equation and the direct discretisation of the weak form, which requires a C1-continuous approximation, are compared in terms of convergence rates. For approximations that are higher than second-order in space, the direct discretisation is found to be superior. Suboptimal convergence rates occur when splines of order p=2 are used. This is validated with a priori error estimates for linear problems. The convergence analysis is completed with an investigation of the temporal discretisation. Second-order accuracy is found for the generalised-α method. This ensures the functionality of an adaptive time stepping scheme which is required for the efficient numerical solution of the Cahn–Hilliard equation. The isogeometric finite element framework is eventually validated by two numerical examples of spinodal decomposition.

Introduction

Phase-field models have become a powerful tool for the modelling of phase transformations and morphological changes in different fields of physics as well as materials and engineering science. Compared to sharp-interface models their advantage is that topological changes are avoided, since interfaces are treated in a diffuse manner which is achieved through a parameter which varies continuously. This phase-field parameter accounts for the different material phases and/or the concentration of the different components. Such an approach allows to fully capture the physics of the individual interfaces without the need to explicitly track them. Typical applications include the modelling and simulation of solidification processes, spinodal decomposition, coarsening of precipitate phases, shape memory effects, re-crystallisation, and dislocation dynamics [1], [2], [3], [4]. Phase-field models have been successfully applied to predict microstructural changes under external fields [5], [6], model tumor growth [7], [8], [9], and image impainting [10]. The phase-field approach has also been used to model crack propagation [11], [12], [13], [14], [15].

Different numerical techniques have been utilised to solve the governing equations of phase-field models [16], including finite difference, finite element, and spectral methods. A difficulty regarding the solution of phase-field models is that they typically involve spatial differential operators that are higher than second-order. Therefore, standard finite elements based on C0-continuous Lagrangian polynomial shape functions do not provide converging solutions when directly applied to the phase-field equation. Instead, mixed formulations [17], [18], where the governing equations are split into a coupled system of lower-order differential equations, can be used. However, this approach introduces additional unknowns and is prone to stability issues [19], [20]. Alternatively, Wells et al. [21] have combined features of continuous and discontinuous Galerkin methods, cf. [22], to use a lower-order continuity basis than is actually required by the weak form. Hence, C0-continuous Lagrangian basis functions can be used, while the first-order inter-element continuity is weakly enforced by Nitsche's method. Although efficient as no additional unknowns are introduced, the approach changes the assembly procedures for the global system of equations and requires a careful choice of the penalty parameter. Starting with the work of Gómez [23] a new approach to the numerical modelling of structural evolution processes has emerged, which combines phase-field models with spline-based approximations. Such an isogeometric analysis framework allows for an accurate and efficient resolution of steep gradients that can occur when higher-order derivatives are introduced. Moreover, the higher-order continuity that is provided by these approximations allows for a direct discretisation without additional degrees of freedom. Differently from earlier approaches based on Hermite elements [24], [25] it can straightforwardly be applied to two- and three-dimensional problems.

Analytic solutions for multidimensional phase-field problems are generally not available. For this reason, the comparison of different discretisations in terms of a quantitative convergence analysis is cumbersome. Approaches commonly use linearised or one-dimensional models in order to carry out accuracy tests. Alternatively, the numerical result that stems from the finest discretisation is considered as the reference solution [26]. Zhang et al. [27] have presented a convergence analysis based on the method of manufactured solutions [28]. By assuming an appropriate solution and extending the underlying differential equation by a source term, different discretisations of the same phase-field model can be compared, see also Aristotelous et al. [29] and Kay et al. [30]. To the authors' knowledge a quantitative convergence study for diffuse interface models has never been carried out within the framework of isogeometric analysis. For this reason, the Cahn–Hilliard equation is used here as a prototypical higher-order equation, and the numerical properties of different discretisations are analysed by a manufactured solution approach.

Originally derived to model spinodal decomposition of binary mixtures [31], [32], the Cahn–Hilliard equation is one of the most commonly used diffuse interface models. Taking the concentration c of one of the mixtures' components as an appropriate phase-field parameter, the governing equation can be stated asct=[M(c)(f(c)λc,kk),l],l. Herein, λ is the interface parameter that governs the width of the interface which separates two adjacent phases, M(c) is the concentration dependant mobility, while () and (),k designate the derivative ()/c and the spatial derivative ()/xk, respectively. Eventually, f(c) represents the configurational or bulk free energy, which can be described by the generalised logarithmic potential [33]f(c)=A[cln(c)+(1c)ln(1c)+Bc(1c)] with constants A and B. The bulk free energy accounts for phase separation which dominates the early stages of the evolution process. A second term, called the gradient or the interfacial free energy, contributes to the total Ginsburg–Landau free energyF=Fbulk+Fint=Ω(f(c)+λ2c,k2)dV and accounts for the influence of the individual interfaces on the system. While the total free energy F is minimised in the course of the structural evolution, different phenomena can occur. Both quantities, Fbulk and Fint, can be used to explain these phenomena from an energetic point of view.

This paper is organised as follows: In Section 2 different finite element formulations are presented for the numerical solution of the Cahn–Hilliard problem, i.e., a mixed form and a direct approach. For the latter, approximations with at least C1 inter-element continuity are required. In Section 3, these discretisations are compared. The numerical properties of both approaches are analysed in terms of the error norms and the convergence rates by making use of a manufactured solution. In Section 4 numerical examples for the spinodal decomposition are presented.

Section snippets

Finite element method

In this contribution, isogeometric analysis is applied to the spatial discretisation of the Cahn–Hilliard equation. The spline-based discretisations offer higher-order continuity properties and are therefore suitable for the direct discretisation of weak forms involving higher-order spatial derivatives. For the temporal discretisation, the generalised-α method [34], [35] is used. Following Gómez et al. [23], the method allows for the development of an efficient time stepping algorithm based on

Convergence analysis based on a manufactured solution

In this section the properties of the presented FE formulations are compared for different approximations in terms of a quantitative analysis of error levels and convergence rates. As no analytical solution is available, a manufactured solution [27], [30] is utilised. The general idea is to use an arbitrary function as analytic reference solution. Since this function will normally not fulfil the governing differential equations exactly, a residual has to be added to the right hand side of the

Numerical examples

In the following examples the predictions of the mixed formulation and the discretisation for the spinodal decomposition of binary systems will be compared. As the different stages of the structural evolution can be identified by the free energy and its individual parts, special attention is paid to them. Starting from an initial concentration distribution in which small perturbations promote the evolution of the system, each problem is considered until a steady state is reached.

The course of

Conclusion

The present work provides a detailed convergence analysis of the Cahn–Hilliard phase-field model. Two different formulations of the model have been implemented in a Matlab-based finite element code that provides Lagrange polynomials and B-splines. Due to the lack of an analytic reference solution, the method of manufactured solutions has been adopted. With respect to the spatial discretisation, convergence rates are obtained that match analytical error estimates available for linear problems.

Acknowledgements

The present study is funded by the German Research Foundation (DFG), Priority Programme (SPP) 1713, grant KA 3309/5-1. This support is gratefully acknowledged. M. Kästner acknowledges financial support of the German Academic Exchange Service (DAAD) during his research visit to the University of Glasgow.

References (44)

  • L. Zhang et al.

    A quantitative comparison between C0 and C1 elements for solving the Cahn–Hilliard equation

    J. Comput. Phys.

    (2013)
  • O. Wodo et al.

    Computationally efficient solution to the Cahn–Hilliard equation: adaptive implicit time schemes, mesh sensitivity analysis and the 3D isoperimetric problem

    J. Comput. Phys.

    (2011)
  • K.E. Jansen et al.

    A generalized-α method for integrating the filtered Navier–Stokes equations with a stabilized finite element method

    Comput. Methods Appl. Mech. Eng.

    (2000)
  • D. Wang et al.

    An improved NURBS-based isogeometric analysis with enhanced treatment of essential boundary conditions

    Comput. Methods Appl. Mech. Eng.

    (2010)
  • J. Kiendl et al.

    Isogeometric shell analysis with Kirchhoff–Love elements

    Comput. Methods Appl. Mech. Eng.

    (2009)
  • A. Tagliabue et al.

    Isogeometric analysis and error estimates for high order partial differential equations in fluid dynamics

    Comput. Fluids

    (2014)
  • J. Liu et al.

    Isogeometric analysis of the advective Cahn–Hilliard equation: spinodal decomposition under shear flow

    J. Comput. Phys.

    (2013)
  • L.-Q. Chen

    Phase-field models for microstructure evolution

    Annu. Rev. Mater. Res.

    (2002)
  • H. Emmerich

    Advances of and by phase-field modelling in condensed-matter physics

    Adv. Phys.

    (2008)
  • I. Steinbach

    Phase-field models in materials science

    Model. Simul. Mater. Sci. Eng.

    (2009)
  • T. Koyama

    Phase-field modeling of microstructure evolutions in magnetic materials

    Sci. Technol. Adv. Mater.

    (2008)
  • A. Hawkins-Daarud et al.

    Numerical simulation of a thermodynamically consistent four-species tumor growth model

    Int. J. Numer. Methods Biomed. Eng.

    (2012)
  • Cited by (59)

    • A local meshless method for transient nonlinear problems: Preliminary investigation and application to phase-field models

      2022, Computers and Mathematics with Applications
      Citation Excerpt :

      The second term in the integrand of the energy functional (1), penalizes the large gradients and hence accounts for the interfacial energy [89–91]. The interfacial parameter λ governs the thickness of the interface separating adjacent phases [46,92]. The Allen-Cahn equation can be applied to describe the evolution of a nonconserved field variable, resulting in formation of a two-phase configuration in a system [93].

    • A spatio-temporal adaptive phase-field fracture method

      2022, Computer Methods in Applied Mechanics and Engineering
    View all citing articles on Scopus
    View full text