A meshless technique based on generalized moving least squares combined with the second-order semi-implicit backward differential formula for numerically solving time-dependent phase field models on the spheres
Introduction
Many phenomena in the applied and natural sciences can be formulated as partial differential equations (PDEs) on surfaces. We can point out texture synthesis in computer graphics [87], flow and solidification of a thin fluid film [72], brain imaging [67], phase separation patterns for diblock copolymers on spherical surfaces [82] and surfactant distribution on a moving interface [93]. One of the most important mathematical models defined on surfaces is the phase field equation. As mentioned in [54], a phase field model is a mathematical formulation to solve interfacial problems [54]. Additionally, in this model, the phase field or order parameter is introduced to identify one phase from the other [54]. This parameter also takes distinct values (for example 0 and 1 in a binary system) in each of the phase with a smooth transition between both values [54].
Some applications of the phase field models can be observed in solidification dynamics [13], viscous fingering [12], fracture dynamics [38], vesicle dynamics [49] and the other applications [10], [35], [47], [84]. Usually, these models are constructed to reproduce a given interfacial dynamics. For example, in solidification problems, the front dynamics are given by a diffusion equation for either concentration or temperature in the bulk and some boundary conditions at the interface (a local equilibrium condition and a conservation law) which constitutes the sharp interface model [51]. Some well-known phase field models are Allen-Cahn (AC) and Allen-Cahn vesicle equations [56], [97], conservative Allen-Cahn (CAC) equation [50], Cahn-Hilliard (CH) equation [14], [15], [26] and the phase field crystal (PFC) equation [33], [34], [62]. As mentioned in [63], the dynamic of these equations derived as a gradient flow of the free energy, admits an energy dissipation law, which justifies its thermodynamic consistency and leads to a mathematically well–posed model [63]. As it is said in [48], the directed self–assembly of block copolymers in thin films is an emerging technology for nanoscale patterning. A diblock copolymer consists of two blocks, each of a different type of monomer, which are joined chemically to each other. When the temperature is lowered below a critical point, the two sequences become incompatible and the copolymer melt undergoes phase separation. This results in the occurrence of periodic structures such as lamellae, spheres, cylinders, and gyroids [48]. This phase ordering or separation may occur on static or dynamic surfaces [82] such as lipid bilayer membranes [8], crystal growth on curved surfaces [9], and phase separation within thin film [78]. The above explanations can be considered as a reason for studying the phase field models on the sphere.
Finding the analytical solution of PDEs on the surfaces specially the phase field equation is somewhat complicated. Thus, simple and accurate numerical techniques play an important role for solving such problems. In recent years, different numerical methods have been developed for the solution of PDEs on surfaces. Generally, these techniques are divided into intrinsic and embedded, narrow-band methods [40]. The first group of these methods uses coordinates intrinsic to the surface, and for discretization the differential operators a surface based mesh will be applied [16], [17], [30], [31], [32], while the second group extends the partial differential equation (PDE) in to a narrow band around the surface, and then modifies the differential operators so that the solution is restricted to the surface [2], [11], [66]. As mentioned in [40], intrinsic methods have the benefit, in which the obtained discretization scheme is consistent with the dimension of the main problem, but in narrow-band methods, the surface differential operators are posed in extrinsic coordinates so that all coordinate singularities can be avoided, and the standard methods such as finite difference (FD) schemes (see also [25], [27], [39], which show finding the numerical solution of some applicable models via FD methods) and finite element methods (FEMs) can be applied on 3D Cartesian grids and 3D unstructured meshes, respectively [30], [40], [86]. As it is said in [40], these methods also lead to degenerate surface differential operators because they allow diffusion occurs only in directions tangential to the surface, and it makes difficulty for using implicit discretization techniques in time-dependent problems [40]. Due to these problems in solving PDEs on surfaces, Fuselier and Wright [40] provided a new numerical scheme based on radial basis functions (RBFs) to solve equations defined on surfaces. Also, they derived the error analysis of the RBFs approximation for target functions in Sobolev spaces on the two-dimensional smooth embedded submanifold of . Besides, authors of [70] have provided a new error analysis for completing the error analysis given in [40], and they have found a new numerical scheme based on RBFs method combined with a first–order time splitting technique for solving AC equation on surfaces.
Authors of [6] studied the arrangement of N particles on the sphere by solving the PFC equation on the spheres via different numerical methods such as an implicit approach, which describes [6] the surface using a phase field description, a parametric finite element approach and a spectral method based on nonequispaced fast Fourier transforms on the spheres [6]. To solve PFC and CH equations on the sphere numerically, a finite element method (FEM) with the block-preconditioner, which is considered and analyzed for the obtained linear system from the proposed discretization is applied [73], [74]. In [75], authors have solved numerically a type of the PFC model that involves a density and a polarization field on the sphere, and they also have found three types of crystals in their simulations. Authors of [83] presented a method based on the application of a diffuse interface framework for solving two-phase problems involving a material quantity on an interface. They [83] also used FD methods with a block-structured adaptive grid, and they applied the nonlinear multigrid approach to solve nonlinear systems. In [71], a new mathematical model is introduced for amoeboid motion, which has application in cell migration, and it is solved via isogeometric analysis combined with spline-based FEM. In [61], the authors applied a narrow band neighborhood of a curved surface using FD scheme for the solution of PFC equation. Finite element approximation is implemented and analyzed for CH equation on surfaces [29]. Another recent research work for finding the numerical solution of PFC equation on surfaces can be found in [58], in which an efficient linear second-order unconditionally stable direct discretization method is employed. An unconditionally energy-stable second-order time-accurate scheme is applied to solve CH equation on surfaces [57]. The CH equation on surfaces is solved using FD scheme based on the closest point method [42]. Also, this equation is simulated via FEM in space and semi-implicit Euler scheme in time [76]. In [24], a local RBFs method is employed to solve numerically CH, Swift–Hohenberg and PFC equations in two- and three-dimensional spaces. The interested reader can also refer to [40] and references therein for more details.
The main concern in solving the phase field mathematical models such as PFC and CH is to show the energy stability during a long time simulation [20], [36], [46], [55], [63], [64], [91], [94], [95]. Therefore, choosing a proper time discretization plays an important role for obtaining the numerical solution of these equations. Different types of the second-order semi-implicit backward time formula (SBDF2) have been considered for discretization the time variable of the phase field models in the previous works [20], [63], [64], [94]. In such methods, the stabilized term, i.e., a second-order Douglas-Dupont regularization is also added to the proposed time discretization, which ensures the energy stability under certain condition. As can be seen in next section, the PFC and CH models are sixth and fourth-orders nonlinear PDEs, respectively. To avoid solving the nonlinear algebraic equations at each time step, the linear and nonlinear terms can be discretized implicity and explicitly, respectively. Of course, in some research works mentioned here, both linear and nonlinear terms are discretized implicitly, which it can be solved via the nonlinear multigrid method or Newton's iteration procedure.
To the best of our knowledge, we use a SBDF2 method by adding a stabilized term for approximating the time variable of the CH, nonlocal CH and PFC equations defined on the spheres, and a meshless technique based on generalized moving least squares (GMLS) is applied to approximate the spatial variables. In this technique, there is no need to use a background mesh or triangulation for approximation [37], [69], [89], [90], which can be applied on the studied equations simply. For the first time, the moving least squares (MLS) approximation on the unit sphere was introduced by Wendland in [90]. The implementation of this approach may not be simple for solving PDEs on spheres and the other manifolds because the PDE operators should be approximated by non-closed-form and complicated shape functions [69]. Therefore, Mirzaei in his research work [69] introduced a direct approximation based on GMLS [68] on the spheres with an error analysis. Another formulation of GMLS approximation based on projected gradient of the shape functions is given in [23] for solving reaction-diffusion equations on surfaces. Recently, authors of [45] have introduced new formulations for GMLS approximation defined on different manifolds, and they have employed this technique for finding the numerical solutions of hydrodynamic flows. For more details about the GMLS approximation on spheres and the other manifolds, the interested reader can refer to [23], [45], [69], [90].
To see more studies about MLS approximation and its improvements, here we mention some important research works as follows. The improved moving least squares (IMLS) approximation is considered for introducing the boundary element–free method (BEFM), which could be used to solve 2D Helmholtz problems efficiently [18]. A combination of the field integral equation (FIE) with the complex variable moving least squares (CVMLS) approximation is done in [19] that is called a complex variable boundary element–free method (CVBEFM), and it is applied for solving the exterior Neumann problem of the Helmholtz equation [19]. The complex variable element–free Galerkin (CVEFG) method is introduced and applied for solving 3D elliptic problems and 3D wave equations in [59]. This technique uses the CVMLS approximation to form shape functions. Also, in [96], the variational multiscale interpolating element–free Galerkin (VMIEFG) method is developed and applied to obtain the numerical solution of the nonlinear Darcy-Forchheimer model, in which the interpolating MLS approximation is used to construct the shape functions. Recently, authors of [1] have applied the interpolating stabilized element–free Galerkin technique for solving the Oldroyd model as a generalized incompressible Navier-Stokes equation. Besides, the uniqueness and existence of the developed numerical method in [1] are investigated. Also, the interested reader can refer to [18], [19], [59], [96] for more details.
The CH equation on the sphere of radius is formulated as follows [42], [57], [73], [76] with suitable initial conditions for continuous functions and , and . is the well-known Laplace-Beltrami operator. Also, denotes the concentration of one component of a binary mixture, and it is the order parameter. is the chemical potential, is a free energy density, is the mobility, and ϵ is a positive constant related to interfacial thickness [42], [57], [73], [76]. By considering , Eq. (1.1) is obtained from a constrained gradient flow in the Hilbert space of the Helmholtz free energy functional as follows [42], [57], [73], [76] This model is used to describe a wide number of phase separation processes from co-polymer systems to lipid membrane [42]. Also, it is applied for spinodal decomposition of two component alloys [76], and it is considered for mathematical description of a quick quench scenario followed by that of a much slower coarsening [76]. As it is said in [60], one of the vital concepts of the CH equation is the interface between two phases such as α and β. It has a finite thickness where the composition u changes gradually. When the binary system approaches the equilibrium state composed of α phase with and β phase with , the regions where and correspond to the α and β phases, respectively whereas the region where varies gradually from to denotes the interface between α and β phases. Some important applications of CH equation can be observed in image processing [28], planet formation [85] and cancer growth [65]. Other applications of CH equation can be found in [60] and references therein.
The crystal structure is a crystal (crystalline solid) based on the microscopic arrangement of atoms inside it. When the atoms find a periodic organization, a crystal is a solid. There are three different types of a crystal, which contain large crystal such as snowflakes, diamonds, and table salt, polycrystals such as most metals, rocks, ceramics, and ice and amorphous solids such as glass, wax and many plastics (see [98], where most materials of this paragraph are taken from [98]).
The PFC equation was first introduced by Elder and his co-workers [33], [34], [81] that is a model for simulation of imperfections in material properties such as vacancies, grain boundaries and dislocations. This model also has some important properties, which can be seen in [7], [33], [34], [81].
The PFC equation on the sphere of radius is formulated as follows [5], [6], [61], [73] with the proper initial conditions for continuous functions , and , and , in which F is known as the free energy density. For this equation, the function f is defined as . is the Laplace-Beltrami operator that is defined on . Also, represents density field, is the chemical potential, where represents free-energy, which it is defined as follows [73] where ϵ is a positive constant [5], [6], [61], [73].
As mentioned in [6], the original formulation of the PFC equation (1.3) returns to one of the Hilbert problems that is the distribution of N points on a sphere. This problem is formulated as an optimization problem, which is difficult to be solved numerically when N increases [6]. Instead of solving this problem, authors of [6] derived a mathematical formulation based on PDE defined on surface, i.e., PFC equation to study the ordering of interacting particles on a sphere. For further information, the interested reader can refer to [6] and references therein.
This equation has many applications, which can be addressed as follows. A simulation of heterogeneous nucleation and growth in an undercooled melt [43], crystal growth in a supercooled liquid [44], the PFC simulations in nanostructure formation for example, in semiconductor nanowires, and it also can be applied for the quantitative simulations: critical wavelength. Moreover, to simulate the solidification microstructures formation, they play important role in various applications such as commercial casting [35], and it is known as eutectic and dendritic solidification. The potential application (epitaxial growth) of the PFC model can be found in the technologically important process of thin film growth [35].
The reminder of this paper is as follows. The mass conservation and energy stable properties of both models are proved in Section 2. In Section 3, a SBDF2 scheme with adding a stabilized term, i.e., (a second-order Douglas-Dupont regularization) is given to the time discretization of CH and PFC equations, where A is a positive constant. In Section 4, we show that the time discretization presented here satisfies the mass conservation property, and also energy stability under the assumption on A. The error analysis of the proposed time discretization is given in Section 5. In Section 6, the GMLS approximation on the spheres is discussed briefly. In next section, i.e., in Section 7, the fully discrete schemes of the equations are derived. Some numerical results and simulations are presented in Section 8 for CH, nonlocal CH and PFC equations defined on the two-dimensional spheres. Finally some concluding remarks are collected at Section 9.
Section snippets
The mass conservation and energy stability properties
In this section, we have shown the mass conservation and the energy stability properties for Eqs. (1.3) on the spheres of radius due to the time variable. In [57], these properties are proved for Eq. (1.1).
We need the following lemma, namely “Green-Beltrami identity”, which will be used to illustrate the mass conservation and energy stability properties for the mathematical models studied here.
Lemma 2.1 [4] For any and , the following relation is concluded.
The time discretization
In this section, a SBDF2 scheme is considered for discretization Eqs. (1.1) and (1.3) in time. By dividing the time interval into M sub-intervals such that , where Δt is the time step and by defining , the SBDF2 time-stepping approximation (by adding a second-order Douglas-Dupont regularization) for Eq. (1.1) can be written as follows where A is a positive constant and . Also, and
Stability analysis
Here, we only show the mass conservation and energy stability of the proposed time discretization for Eq. (3.4). The same procedure can be done for Eq. (3.1), which we omitted it for brevity.
Theorem 4.1 The scheme (3.4) is mass conservative. Proof For , according to Eqs. (3.4) and (3.6) and using the definition of inner-product in , we can write Therefore Using Green-Beltrami identity for terms in the
Error estimate
Our goal of this section is to present the error estimate for the time discretization proposed here. We only give the error estimate of the time discretization for PFC equation. This also can be done in a similar manner for CH equation. We first note that using a careful Taylor expansion gives where is independent of Δt.
Assumption 3 We assume that the initial value of u has sufficient regularity. With this, we can suppose that the exact solution of Eq. (1.3) satisfies the following
The GMLS approximation on the spheres
The GMLS approximation on the d-dimensional spheres has been introduced by Mirzaei [69]. Here, we briefly review this technique for the two-dimensional spheres, i.e., to approximate the Laplace-Beltrami operator. More details on this subject can be found in [69], where most materials of the current section are taken from that.
Suppose that is a function defined on the two-dimensional sphere for some . As it is said in [69], the smoothness of a function u on the spheres (or any
The fully discrete scheme
The GMLS technique presented in Section 6 will be applied to approximate Eqs. (3.4) and (3.6) (i.e., obtaining the fully discrete scheme for PFC equation) with respect to the spatial variables. It also should be noted that the fully discrete scheme for CH equation can be obtained in the same manner, which we omit it for brevity.
Note that we use Eqs. (6.7) and (6.5) without considering sign “⁎”. We now suppose that is a set of N scattered points on of radius . The
Numerical results
In this section, some numerical simulations are provided for demonstrating the ability of the developed numerical scheme in solving the CH, nonlocal CH and PFC equations on the two-dimensional spheres. The two-dimensional sphere of radius is defined as follows The spherical harmonics of at most degree are used in GMLS approximation [4], [69]. Besides, the weight function is chosen to be a compactly supported Wendland's function, and it is defined as follows
Conclusion
The generalized moving least squares approximation in combination with a second-order semi-implicit backward differential formula has been successfully applied to solve numerically the high-order nonlinear time-dependent phase field crystal (PFC) and Cahn-Hilliard (CH) equations on the two-dimensional spheres. In order to ensure the energy stability of the time discretization presented here, we followed the literature, and we added a stabilized term that is the second-order Douglas–Dupont
Acknowledgement
The authors are very grateful to three reviewers for carefully reading this paper and for their comments and suggestions, which have improved the quality of the paper.
References (98)
- et al.
Investigation of the Oldroyd model as a generalized incompressible Navier–Stokes equation via the interpolating stabilized element free Galerkin technique
Appl. Numer. Math.
(2020) - et al.
Transport and diffusion of material quantities on propagating interfaces via level set methods
J. Comput. Phys.
(2003) - et al.
Energy stable and efficient finite–difference nonlinear multigrid schemes for the modified phase field crystal equation
J. Comput. Phys.
(2013) On spinodal decomposition
Acta Metall.
(1961)- et al.
The boundary element–free method for 2D interior and exterior Helmholtz problems
Comput. Math. Appl.
(2019) - et al.
A complex variable boundary element–free method for the Helmholtz equation using regularized combined field integral equations
Appl. Math. Lett.
(2020) - et al.
An energy stable fourth order finite difference scheme for the Cahn-Hilliard equation
J. Comput. Appl. Math.
(2019) - et al.
Polynomial approximation in Sobolev spaces on the unit sphere and the unit ball
J. Approx. Theory
(2011) - et al.
The meshless local collocation method for solving multi–dimensional Cahn-Hilliard, Swift–Hohenberg and phase field crystal equations
Eng. Anal. Bound. Elem.
(2017) Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices
Math. Comput. Simul.
(2006)
A numerical method based on the boundary integral equation and dual reciprocity methods for one–dimensional Cahn–Hilliard equation
Eng. Anal. Bound. Elem.
Finite element approximation of the Cahn-Hilliard equation on surfaces
Comput. Methods Appl. Mech. Engrg.
Cahn–Hilliard on surfaces: a numerical study
Appl. Math. Lett.
An unconditionally energy–stable method for the phase field crystal equation
Comput. Methods Appl. Mech. Engrg.
Meshfree methods on manifolds for hydrodynamic flows on curved surfaces: a generalized moving least-squares (GMLS) approach
J. Comput. Phys.
Phase field modelling of stressed grain growth: analytical study and the effect of microstructural length scale
J. Comput. Phys.
A conservative Allen–Cahn equation with a space–time dependent Lagrange multiplier
Int. J. Eng. Sci.
Continuous and discrete least-squares approximation by radial basis functions on spheres
J. Approx. Theory
Multiphase image segmentation using a phase–field model
Comput. Math. Appl.
An efficient and stable compact fourth-order finite difference scheme for the phase field crystal equation
Comput. Methods Appl. Mech. Engrg.
Fast local image inpainting based on the Allen–Cahn model
Digit. Signal Process.
An unconditionally energy-stable second-order time-accurate scheme for the Cahn–Hilliard equation on surfaces
Commun. Nonlinear Sci. Numer. Simul.
An efficient linear second order unconditionally stable direct discretization method for the phase–field crystal equation on surfaces
Appl. Math. Model.
Three–dimensional complex variable element-free Galerkin method
Appl. Math. Model.
Physical, mathematical, and numerical derivations of the Cahn–Hilliard equation
Comput. Mater. Sci.
A simple and efficient finite difference method for the phase–field crystal equation on curved surfaces
Comput. Methods Appl. Mech. Engrg.
First and second order operator splitting methods for the phase field crystal equation
J. Comput. Phys.
A second-order, uniquely solvable, energy stable BDF numerical scheme for the phase field crystal model
Appl. Numer. Math.
A benchmark for the surface Cahn–Hilliard equation
Appl. Math. Lett.
First and second order numerical methods based on a new convex splitting for phase–field crystal equation
J. Comput. Phys.
An adaptive time-stepping strategy for solving the phase field crystal model
J. Comput. Phys.
Variational multiscale interpolating element-free Galerkin method for the nonlinear Darcy–Forchheimer model
Comput. Math. Appl.
Theoretical Numerical Analysis
Spherical Harmonics and Approximations on the Unit Sphere: An Introduction
Particles on curved surfaces: a dynamic approach by a phase-field-crystal model
Phys. Rev. E
A continuous approach to discrete ordering on
Multiscale Model. Simul.
Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension
Nature
Grain boundary scars and spherical crystallography
Science
Diffusive atomistic dynamics of edge dislocations in two dimensions
Phys. Rev. E
Variational problems and partial differential equations on implicit surfaces
J. Comput. Phys.
Phase-field approach to three-dimensional vesicle dynamics
Phys. Rev. E
Phase-field simulation of solidification
Annu. Rev. Mater. Sci.
Free energy of a nonuniform system I, interfacial free energy
J. Chem. Phys.
Logically rectangular grids and finite volume methods for PDEs in circular and spherical domains
SIAM Rev.
A finite volume method for solving parabolic equations on logically Cartesian curved surface meshes
SIAM J. Sci. Comput.
Approximation Theory and Harmonic Analysis on Spheres and Balls
Approximation of continuous surface differential operators with the generalized moving least–squares (GMLS) method for solving reaction–diffusion equation
Comput. Appl. Math.
On the solution of an initial-boundary value problem that combines Neumann and integral condition for the wave equation
Numer. Methods Partial Differ. Equ.
Area preserving curve-shortening flows: from phase separation to image processing
Interfaces Free Bound.
Cited by (15)
Maximum principle preserving and unconditionally stable scheme for a conservative Allen–Cahn equation
2023, Engineering Analysis with Boundary ElementsHighly efficient, decoupled and unconditionally stable numerical schemes for a modified phase-field crystal model with a strong nonlinear vacancy potential
2023, Computers and Mathematics with ApplicationsCitation Excerpt :For these models, various methods have been proposed for solving them. For example, for the CH model, We can find some numerical methods such as multigrid finite element method [1], adaptive multigrid finite difference methods [2], and meshless methods [3–7]. Among all these models, the phase field crystal (PFC) model is a suitable method to capture phase transition phenomena, and because of this, it is widely used in many researches.
Fast evolution numerical method for the Allen–Cahn equation
2023, Journal of King Saud University - ScienceCitation Excerpt :Recently, various extensions of the AC equation have received increased research attention such as the time-fractional AC equation with volume constraint (Ji et al., 2020). In addition, various numerical studies for other phase-field mathematical models have been researched (Rasoulizadeh and Rashidinia, 2020; Mohammadi et al., 2021; Mohammadi et al., 2022; Yadav et al., 2021; Mohammadi and Dehghan, 2015; Mohammadi and Dehghan, 2020; Mohammadi and Dehghan, 2021; Ghassabzadeh et al., 2021; Dehghan and Taleei, 2010). One of efficient numerical methods for the AC equation is the operator splitting method (OSM) (Li et al., 2010; Xiao et al., 2017; Huang et al., 2019; Lee and Lee, 2015; Weng and Tang, 2016; Li et al., 2020; Sun et al., 2019; Ayub et al., 2019).
Electro-chemo-mechanical induced fracture modeling in proton exchange membrane water electrolysis for sustainable hydrogen production
2022, Computer Methods in Applied Mechanics and EngineeringCitation Excerpt :The main driving force for these developments is the possibility to handle complex fracture phenomena within numerical methods in two and three dimensions. In recent years, several brittle [38–66] and ductile [67–84] phase-field fracture formulations have been proposed in the literature. These studies range from the modeling of 2D/3D small and large strain deformations, variational formulations, multi-scale/physics problems, mathematical analysis, different decompositions, and discretization techniques with many applications in science and engineering.
Probabilistic failure mechanisms via Monte Carlo simulations of complex microstructures
2022, Computer Methods in Applied Mechanics and EngineeringCitation Excerpt :The main driving force for these developments is the possibility to handle complex fracture phenomena within numerical methods in two and three dimensions. In recent years, several brittle [10–42] and ductile [7,43–59] phase-field fracture formulations have been proposed in the literature. These studies range from the modeling of 2D/3D small and large strain deformations, variational formulations, multi-scale/physics problems, mathematical analysis, different decompositions, and discretization techniques with many applications in science and engineering.
Radial basis function-generated finite difference scheme for simulating the brain cancer growth model under radiotherapy in various types of computational domains
2020, Computer Methods and Programs in Biomedicine