The phase field crystal (PFC) model[1, 2] was proposed as an approach to simulate crystals at the atomic scale but on a coarse-grained diffusive time scale. Many physical processes, such as the formation of ordered structures, phase separation of polynary systems, can be described using this model. The PFC model can also explain elastic and plastic deformations of the lattice, dislocations, grain boundaries, multiple crystal orientations and many other observable phenomena[4, 6].
There are several kinds of PFC models. In general, they can be classified into two classes according to characteristic length scale: one-length-scale and multi-length-scale. One-length-scale PFC models can be used to describe the phase behavior of periodic structures[7, 8, 9]. Accordingly, multi-length-scale PFC models can be employed to explain the formation of quasicrystals[10, 11]. In this work, we focus on the development of numerical methods of one-length-scale PFC model. In particular, the classic Landau-Brazovskii (LB) model[7, 12, 13]
will be used to demonstrate our proposed method. The LB model was built to investigate the character of phase transition. It has been discovered in many scientific fields. For example, the LB model can be derived from more complicated self-consistent field theory of diblock copolymers. Compared with the typical Swift-Hohenberg (SH) model with double-well bulk energy , LB energy functional includes a cubic term which can be used to study the first-order phase transition.
The (Allen-Cahn) or
(Cahn-Hilliard) gradient flow equation is usually adopted to describe the dynamic behavior of the phase-field or PFC model. These dynamic equations are time-dependent nonlinear partial differential equations (PDEs). It is hard to find non-trivial analytical solutions. Therefore, numerically solving these nonlinear PDEs is an efficient approach. To guarantee convergence, numerical schemes of these equations are required to satisfy the energy dissipation property. Meanwhile, an accurate and efficient approach should be designed to deal with nonlinear terms. In terms of time discretization, there have been several effective methods which can preserve energy dissipation law, including the convex splitting method[4, 15, 16, 17, 18, 19, 20, 21], stabilized approach[22, 23, 25, 24], invariant energy quadratization (IEQ) method[26, 19] and recently developed scalar auxiliary variable (SAV) approach. By introduing a scale auxiliary variable to the nonlinear part of energy functional, the SAV approach has a modified energy dissipation property for a large class of gradient flows. The convergent and error analysis of semi-discrete SAV scheme has been given by Shen and Xu. The analysis of energy stability and convergence of fully discretized SAV block-centered finite difference method has been established for gradient flows. More studies about the PFC problem can be found in recent literature [16, 43, 17, 44, 45].
In the study to the PFC model, finite difference methods[2, 3, 4, 5, 46] or spectral methods[27, 28, 29] are limited to regular regions, such as two-dimensional square region or three-dimensional cube region. For complex geometries, finite element method (FEM)[30, 55, 42] is a better choice. Furthermore, the FEM can be further combined with adaptive technologies, which are well suitable for the phase behavior of PFC models, such as the formation of ordered structures, phase transition processes, and coarse-grained processes. The adaptive method can effectively decrease the cost of computing and accurately capture the phase interface.
In this work, we will combine SAV time discretization and FEM spatial discretization to solve the gradient flow equation of LB model. Based on the energy dissipation and the SAV scheme, the derivation process of bounds of the solution is shown in detail. For our fully discrete scheme, we demonstrate its energy stability, and carry out error estimate. Applying our method, we can effectively simulate the mesoscale self-assembly of the diblock copolymer system in two-dimensional convex geometries. In addition, we will consider an adaptive FEM for the PFC model. There are many adaptive finite element methods for phase field equation[49, 50, 51, 52, 53]. To reduce computational cost, we first apply an adaptive method which is effective against phase field equation to the PFC model. Numerical results demonstrate that directly using the gradient as the indicator is more efficient than the posterior error estimator does in solving this problem. Since the gradient obtained from the numerical solution may be discontinuous, therefore, a smooth recovered gradient is employed as the adaptive indicator in our adaptive FEM.
The rest of this paper is organized as follows. In Sec. 2, we introduce the LB free energy functional and take the gradient flow as an example to derive its Allen-Cahn dynamic equation. Sec. 3 details our numerical method, which consists of discretization schemes, energy dissipation and error estimate for and dynamical equations. In Sec. 4, numerical experiments are given to illustrate the accuracy and effectiveness of our scheme. Several standard ordered structures on two-dimensional convex regions are also obtained in this section. Sec. 5 gives a simple but efficient adaptive FEM to PFC model. In Sec. 6, we give conclusions and outlooks.
2 Physical model
The dimensionless free energy functional of LB model is
where is a two dimensional bounded domain with Lipschitz boundary , is the density deviation of a kind of monomer from the disordered phase, and are the parameters of the model, is the Laplace operator.
The gradient flow of LB model is
where is a negative operator. For gradient flow , while for gradient flow, . For brevity, in the following, we take the gradient flow as an example to derive the Allen-Cahn dynamical equation. The gradient flow equation can be derived similarly.
Now we will give the boundary conditions to guarantee the energy dissipation property of the Allen-Cahn dynamical equation. Denote that
therefore, the LB model becomes
The free energy take a derivative with respect to time is
here we introduce two Neumann boundary conditions
Define a function space as
In the sense of the Gateaux differential for all , we have the following equation
where , denote the inner product. According to the variational principle, we have
By introducing a new function , we can write Eqns. (11) as the coupled system
3 Numerical Methods
The main aim of this section is to propose numerical methods to solve the gradient flows of LB model. For brevity, for the Allen-Cahn equation, we first present the discretization scheme, prove the energy stability, and give the error estimate in detail. Then the corresponding results about the Cahn-Hilliard equation are also given.
3.1 SAV discretization
Recently, Shen et al. proposed an efficient time discretization scheme, i.e., the SAV scheme, to a class of gradient flows. Here, we shall apply the idea of the SAV approach to discretize Eqns. (12) in time direction.
Let . Then we introduce the scalar auxiliary variable , where is a constant to ensure , and write Eqns. (12) as
Assume is a fixed time step and is the approximation of at time , we can construct the first-order SAV scheme:
Remark 3.1 The higher-order schemes based on SAV technique can be easily constructed, see Ref. for more details.
3.2 FEM discretization
We discretize Eqns. (14) in space using the FEM. Let denote both the trial and test function spaces
The corresponding Galerkin form of Eqns. (14) can be stated as follows: for , find such that:
Let be a comforming mesh of with , be line segment in 1D or triangle in 2D, and be the linear finite element space over defined by
Thus, Eqns. (16) is transformed as follow: find , such that for :
3.3 Energy stability
We shall prove the energy stability of Eqns. (17). The norm of is denoted by .
Theorem 3.1 If we denote the modified energy
then Eqns. (17) is unconditionally energy stable with the modified energy.
Proof: We take in Eqn. (17a) , and find
According to Eqn. (17b), we have
Setting to Eqn. (20) and using the identity
Multiplying Eqn. (17c) with , we obtain
Then the discretized energy dissipative property is satisfied, i.e.
The energy stability of Eqns. (17) is derived for the modified free energy , not for the original one , owing to the introduction of a function and a scaler variable in the constructing processes of Eqns. (17). The modified free energy plays an important role in the error analysis of the discrete scheme.
3.4 Error estimate
We first give some lemmas below.
Lemma 3.1 Assume . Let be solutions of Eqns. (13). There exists a constant depending only on and such that,
We take the inner product of Eqn. (13a) with , then for all ,
Integrating Eqn. (13b) with , we obtain
Thanks to Eqn. (26), then
By Minkowski’s inequality and Sobolve embedding, and choosing the appropriate such that we have
Next we shall prove using mathematical induction.
When , .
If , , by Minkowski’s inequality and Hölder’s inequality, we have the results as follow:
For Eqns. (17), we can obtain Lemma 3.2, as its proofs essentially the same as for Lemma 3.1.
Lemma 3.2 Assume . Let be solutions of Eqns. (17). Then for all and , we have
where and is the discrete Laplace operator.
Here, we shall derive error estimates of Eqns. (17). Denote , , , where .
Let be the standard elliptic (Ritz) projection operator, satisfying
For , such that
Proof: By Lemma 3.1 and Lemma 3.2, we know that
Note that . Therefore, we can find a constant such that
Using Hölder’s inequality, we get