Analysis and simulation for an isotropic phase-field model describing grain growth

A phase-field system of coupled Allen–Cahn type PDEs describing grain growth is analyzed and simulated. In the periodic setting, we prove the existence and uniqueness of global weak solutions to the problem. Then we investigate the long-time behavior of the solutions within the theory of infinite-dimensional dissipative dynamical systems. Namely, the problem possesses a global attractor as well as an exponential attractor, which entails that the global attractor has finite fractal dimension. Moreover, we show that each trajectory converges to a single equilibrium. A time-adaptive numerical scheme based on trigonometric interpolation is presented. It allows to track the approximated long-time behavior accurately and leads to a convergence rate. The scheme exhibits a physically aspired discrete free energy dissipation.


Introduction
There exist a vast amount of literature concerned with the modeling and simulation of grain growth and investigations are being continued. Especially for the improvement of modern solarcells, one is interested in realistic models for different kinds of re-crystallization phenomena, as it takes very long to improve costs and efficiency in terms of the usual trial and error approaches. Because of the large quantity of publications we refrain from giving a detailed overview of the methods here, instead we would like to refer to one of several existing review articles, by Moelens et al. [23], and give only a short introduction on the topic. Sharp-interface models, e.g. those treated by the group of Kinderlehrer [20], prescribe evolution laws to the interfaces between neighboring grains. As this leads to a rather complicated treatment of topological transformations or the incorporation of contributions from the bulk, phase-field approaches have become popular in recent years. Here, the interface is smoothed out, and the transition happens on a characteristic, small length-scale. Topological changes, such as the vanishing of a grain and the redistribution of the interface network, happen with no additional effort. These kind of models can be divided into two groups, those that consider an order parameter as indicator for the crystallization in a material and an orientation indicator, which attains for each grain a different value (see e.g. [26]), and those that assign a phase-field variable to each of the grain orientations (see e.g. [5,18]) -both have their pros and cons. In this work we shall concentrate on the isotropic phase-field model from Chen and Yang [5] formulated in terms of the latter approach. We shall perform a complete analysis on the PDE system (well-posedness, long-time behavior) and provide a time-adaptive Fourier collocation scheme to simulate the isotropic grain growth model. We remark that the model [5] has been extended throughout the years, i.e. to anisotropic interfacial energies [18] or recently by incorporating elasticity [3]. Our work lays a fundament for the analysis of those extended models, which we plan to continue in the near future.
We assume that there are n g ∈ N grain orientations. As motivated by Chen and Yang [5] the evolution of corresponding grain boundaries is described by the time-dependent vector u = (u i ) i=1,...,ng of order parameters. Each order-parameter field u i represents grains of a given crystallographic orientation. As also introduced in the originating paper, the corresponding potential energy reads Here, ω and λ are phenomenological (positive) parameters. They are usually assumed to fulfill λ > ω such that the minima of (1.1) are located at the Cartesian basis vectors ±e j ∈ R ng . Like in the Cahn-Hilliard equation, the transition between the grains are smoothed out by introducing suitable gradient terms in the energy functional W . Then the total energy for the isotropic grain growth system reads where ν i > 0 are (constant) gradient coefficients. The kinetics of a quenched system and the relaxation of the n g order-parameter fields are described by Ginzburg-Landau type equations (here we ignore possible noise terms, see [5]): Here M i > 0 are constant mobilities for the grain orientations. Without loss of generality, we simply set M i = 1 in the remaining part of the paper, as different constant mobilities only lead to some more algebraic work without further mathematical difficulties. Our results hold analogously true with these mobility terms. Thus we arrive at the following PDE system that can be viewed as a gradient flow of the energy E, namely, Here, we impose periodic boundary conditions on the phase functions u i in the domain Ω x ∈ Ω, t ≥ 0, i = 1, . . . , n g , j = 1, . . . , n.
(1.5) Besides, the system (1.4) is equipped with the following initial conditions The periodicity assumption (1.5) is reasonable for larger bulk materials, in view of the fact that the evolving patterns repeat statistically on the, in total, much larger material piece. In this setting we are able to carry out numerical simulations with a pseudospectral method similar to that in Chen and Shen [6]. As visualized in Figure 1, the periodic boundary conditions lead to a realistic growth of the grain network. We see how nuclei grow to real grains that eventually fill the whole rectangular domain. These coarsen on a much slower time-scale, which will be discussed to some detail in the numerics Section 4. The isotropy of the surface energy density and the mobility leads to circular growth of the structures, until neighboring crystalline regions touch. This is well visible in the magnifications of two neighboring grains in Figure 1. For improvement of the runtimes one can implement a reassignment routine, e.g. as in [19]. The plan of the paper is as follows. In Section 2, we show the existence of weak solutions to the system (1.4)-(1.5). Then we prove that these solutions are in fact uniquely defined by the initial condition and that the evolution depends continuously on the data. In Section 3, we study the long-time behavior of the problem (1.4)-(1.6). We prove that it generates a strongly continuous semigroup on a suitable phase-space and the associated dynamical system possesses a global attractor and an exponential attractor. Besides, we show that each trajectory converges to a single equilibrium with respect to the phase-space metric. The numerical scheme is presented in the last Section 4. We show the existence of a discrete analogon of the free energy and that it also dissipates in the Fourier collocation setting. Furthermore we show that the derived analytical convergence rate corresponds well to the rate observed in simulations.

Global Weak Solutions
For the domain Ω = Π n i=1 [0, L i ] (n ≤ 3) and arbitrary m ∈ N, we denote by H m p (Ω), m ∈ N, the space of functions which are in H m loc (R n ) and periodic with the period defined by the extents of Ω. H m p (Ω) is a Hilbert space for the scalar product (u, v) H m = |κ|≤m Ω D κ u(x)D κ v(x)dx (κ a multi-index) and its associated norm is given by (Ω) and the inner product as well as the norm on L 2 p (Ω) are simply denoted by (·, ·) and · , respectively. For two suitable Banach spaces V and H, we write for a vector valued function u: u ∈ (V (0, T ; H)) ng , if for each i ∈ {1, . . . , n g }, u i ∈ V (0, T ; H). In the remaining part of the paper, we assume that ω > 0 and λ > 0.
First, we observe that the phase-field system (1.4)-(1.6) obeys a dissipative energy law, which turns out to be important in the study of well-posedness and long-time behavior.

2)
such that for i = 1, . . . , n g and any ϕ ∈ H 1 p (Ω) Based on the energy dissipation law (2.1), existence of weak solutions can be proved by using the classical Galerkin approach (cf. e.g., [29,Chapter 3] for the scalar Allen-Cahn equation). As the only subtlety lies is in treating the coupling nonlinearity properly, we omit the details here. Next, we show the continuous dependence result on the initial data, which also implies the uniqueness of weak solutions to problem (1.4)-(1.6).

4)
where C is a constant depending on u 0 H 1 , v 0 H 1 , λ, ω, ν i , Ω, but not on t. As a consequence, the global weak solutions to problem (1.4)-(1.6) obtained in Theorem 2 are unique.
Proof Denote w i = u i − v i for i = 1, ..., n g . We subtract the equations for the two solutions u, v such that Testing each equation with w i − ∆w i , we obtain that Since ω is a positive constant, it is easy to see that Next, using the fact (2.2), the Sobolev embedding theorem (n ≤ 3) such that H 1 p (Ω) → L 6 p (Ω) and the Hölder inequality as well as Young inequality, we have For the nonlinear coupling terms on the right-hand side of (2.6), we observe that for every Using the above equality, the solution's regularity (2.2), the Sobolev embedding theorem (n ≤ 3) such that H 1 p (Ω) → L 6 p (Ω) and the Young inequality, we calculate for the coupling nonlinearity Summing up (2.6) with respect to i and taking in (2.7) and (2.8) sufficiently small, such that Then our conclusion (2.4) follows from the Gronwall inequality used for (2.9). Moreover, the continuous dependence result (2.4) yields the uniqueness of weak solution to problem (1.4)-(1.6).
Before ending this section, we provide a useful observation of the order parameters u i , i.e., the following weak maximum principle: Proposition 4 Suppose that the assumptions in Theorem 2 are satisfied.
(1) Assume in addition that Proof For any scalar or vector function v(x), we denote v + to be its positive part, i.e, v + = max{0, v(x)}. We now consider the following two cases: .., n g ) and summing up with respect to i, we get We infer from the definition of h(x, t) and the assumption The above estimates yields that d dt First, from the expression of the nonlinearities f i , we easily see that Then using the idea in [15, Chapter 4], we denote A direction computation yields that It follows from the assumption on the initial data that g(0) = 0. Then the above inequality implies that g(t) = 0 for t ≥ 0, which yields the required result.

Long-time Behavior
In the previous section, we have proved the existence and uniqueness of global weak solutions to problem (1.4)-(1.6). Our aim of this section is to investigate its long-time behavior. Therefore denote the solution operator We infer from Theorem 2 and Theorem 3 that problem (1.4)-(1.6) generates a (strongly continuous) semigroup S(t) on the phase space (H 1 p (Ω)) ng .

Global Attractor
First, we show that the semigroup S(t) possesses a global attractor.
For this purpose, the preliminary step is to prove the dissipativity of the semigroup, namely, S(t) possesses an absorbing set in (H 1 p (Ω)) ng . For any R ≥ 0, m ∈ N, we denote by B m (R) the bounded ball in (H m p (Ω)) ng centered at zero with radius R.
Lemma 6 There exists R 1 ≥ 0 such that the ball B 1 (R 1 ) is an absorbing set in (H 1 p (Ω)) ng . That is, given any R ≥ 0, there exists a time t 1 = t 1 (R) such that (3.1) Proof Testing system (1.4) with u i , summing over all i (i = 1, ..., n g ) and adding the resultant with (2.1), we obtain that d dt By the Young inequality, we easily see that where c 1 is a constant depending on ω, λ, ν i , Ω, but not on u i . Moreover, using the Young inequality once more and (3.2), we can find a sufficiently small constant γ ∈ (0, 1] that may depend on ω, λ, ν i , Ω, but not on u i such that γE 1 ≤ D 1 + 1. As a result, we have which easily yields The desired H 1 -absorbing property follows from the above estimate and (3.3).
Next, we show the existence of a compact absorbing set. The calculations performed in hereafter is formal and can be easily justified by the Galerkin approximation.
Lemma 7 With assumptions and terminology as in Lemma 6, there is a positive constant R 2 independent of the initial data and time and, for any Here and in what follows C stands for a positive constant that is independent of t and of the initial data. This constant may vary from line to line. Besides, Q denotes a positive nondecreasing monotone function that is independent of t.
For the sake of simplicity, here we denote f (u i ) = u 3 i − u i , i = 1, ..., n g . Testing the system (1.4) by −(∆u i ) t , after integration by parts, we get The right-hand side of (3.7) can be estimated as follows We deduce from (3.7) and (3.8) that Using the equation (1.4), the dissipative estimates (3.5) and (3.3), we infer that as well as is the time as in Lemma 6) such that Q(R)e −γt 2 ≤ 1. Then using the so-called uniform Gronwall lemma (see, e.g., [25, Chap. III, Lemma 1.1]), we infer from (3.9), (3.10) that The above estimates together with (3.3) and (3.5) yield the existence of constant R 2 and t 2 =t 2 + 1 such that (3.6) holds.

Exponential Attractors
We will show that the fractal dimension of the global attractor A is in fact finite. This property results from the existence of an exponential attractor which is a finer attracting set. The terminology used below can be found in [25].
Here dist H 1 denotes the Hausdorff distance between sets in (H 1 p (Ω)) ng .
As a direct consequence of Theorem 8, we conclude that

Corollary 9
The global attractor A obtained in Theorem 5 has finite fractal dimension.
In order to prove Theorem 8, we first show certain smoothing property of weak solutions to problem (1.4)-(1.6).
Lemma 10 Suppose that B 1 (R 1 ) is an absorbing ball of (H 1 p (Ω)) ng such that S(t) We have the following smoothing estimate where C and L are positive constants that only depend on R 1 , Ω and on the structural parameters ω, λ, ν i .
Proof As in the proof of Theorem 3, we denote w i = u i − v i for i = 1, ..., n g and w = (w i ) = u − v. For i = 1, .., n g , differentiating (2.5) with respective to time and testing the resultant by (w i ) t , we obtain (3.14) For t ≥ t 2 , we can use the estimate (3.6). Applying the Sobolev embedding H 2 → L ∞ (n ≤ 3) and the Hölder inequality, the right hand side of (3.14) can be estimated as follows Then summing (3.14) over i = 1, ..., n g , we obtain that On account of (3.12) and (2.1), we have On the other hand, we infer from the uniform estimate (3.6) that Then it follows from the equation (2.5) and the standard elliptic regularity theory that Using the estimate (3.16) and the continuous dependence result (2.4), we deduce from (3.15) and the uniform Gronwall inequality that for arbitrary but fixed t ≥ t 2 Our conclusion (3.13) follows from (3.18) and the elliptic estimate (3.17).
Next, we state the Hölder continuity with respect to time of the trajectories defined by the semigroup S(t).
Lemma 11 Let B 1 (R) be a bounded ball in (H 1 p (Ω)) ng . Then there is a positive nondecreasing monotone function Q and t * = t * (R) > 0 such that
Proof It follows from (2.1) that On the other hand, Lemma 7 yields that there exist t * = t * (R) > 0 such that (3.6) holds. Using the above estimate, standard interpolation inequalities, and the Hölder inequality, we have for t >t ∈ [t * , +∞) which easily yields (3.19).
In order to prove Theorem 8, we will employ a well-known result from Efendiev et al. [7] that shows under which assumptions one expects an exponential attractor with finite fractal dimension for discrete semigroups generated by the iterations of a proper map S. We state the result and apply it for our needs in a second step for the semigroup S(t) with continuous time.
Lemma 12 (cf. [7]) Let H and H 1 be two Banach spaces such that H 1 is compactly embedded into H. Suppose that B is a bounded closed subset of H. Let us give a map S : B → B such that  [8]). Finally, we infer from the continuity of S(t) (cf. (2.4) and (3.19)) that M satisfies the properties (i)(ii)(iii) stated in Theorem 8.

Convergence to equilibrium
The dissipative energy law (2.1) (see Lemma 1) implies that E(u) is a (global) Lyapunov functional for the semigroup S(t) defined by problem (1.4)-(1.6). As a consequence, the associated dynamical system (S(t), (H 1 p (Ω)) ng ) is a gradient system. We notice that u ∞ = (u ∞ i ) ∈ (H 1 p (Ω)) ng belongs to the set of steady states S if and only if it is a solution to the elliptic system subject to periodic boundary conditions x ∈ Ω, i = 1, . . . , n g , k = 1, . . . , n. (Ω)) ng by using the method of minimizing sequence for the associated energy E. Moreover, by the standard elliptic regularity theory, the stationary solution u ∞ is smooth. Recall that Lemma 7 implies the relative compactness of any trajectory in (H 1 p (Ω)) ng . Then we conclude from the well-known results in dynamical system (cf. e.g., [25,Lemma I.1

.1]) that
Proposition 13 For any u 0 ∈ (H 1 p (Ω)) ng , the ω-limit set ω(u 0 ) ⊂ S is a nonempty compact connected subset of (H 1 p (Ω)) ng . Furthermore, we have We see that the global attractor A obtained in Theorem 5 coincides with the unstable manifold of the set S of the stationary points (cf., e.g., [25,Chapter 7,Section 4]). In addition to the equilibria, A also contains heteroclinic orbits connecting different equilibria. However, the global attractor in general does not provide information on the asymptotic behavior of single trajectories. In particular, the set S can be a continuum (cf., [13]) and the ω-limit set ω(u 0 ) may not consists of a single point. However, since the nonlinearities of problem (1.4)-(1.6) are real analytic (polynomials), using the fact that the dynamical system (S(t), (H 1 p (Ω)) ng ) is a gradient system, we can apply the Lojasiewicz-Simon technique (cf. [15,16] and the references therein) to show that each trajectory does converge to a unique stationary state. This constitutes the main result of this section. The detailed proof is similar to the scalar parabolic equations [14,16,29] (see also generalized result on parabolic systems in [15, Chapter 4, Section 3]) and thus is omitted here.
Theorem 14 For any initial data u 0 ∈ (H 1 p (Ω)) ng , the unique global weak solution to problem (1.4)-(1.6) converges to a single equilibrium u ∞ in the topology of H 1 . Moreover, there exist θ ∈ (0, 1/2) depending on u ∞ and C ≥ 0 depending on u 0 , u ∞ such that Remark The convergence rate (3.23) follows from the argument in [14,29] (for the L 2 -norm) and energy estimates (for higher-order norms, like e.g., in [27]). It is worth noting that, by using the smoothing property of the solutions to parabolic systems, the convergence result as well as the convergence rate estimate (3.23) can be demonstrated with respect to higher-order norms H m (for t ≥ δ > 0).

A time-adaptive Fourier collocation scheme
In the rest of the paper we show how the discussed system can be simulated numerically with a global interpolation method. Besides, due to an adaptive time-stepping procedure we are able to simulate its long-time behavior accurately.

Global interpolation and time-adaptivity
To implement a semi-implicit scheme for the simulation of the isotropic grain growth model (1.4) that is capable to deal with the periodic boundary conditions (1.5), we decompose the system into linear contributions L i and nonlinear parts N i , Systems of this kind have been simulated frequently in a semi-implicit fashion, e.g. [1,17,22], and here we explain how we exploit this structure. By treating the linear derivatives implicitly and all nonlinearities in an explicit fashion we can get sufficient stability and accuracy with help of an adaptive SBDF1/SBDF1-2step method. We apply a trigonometric interpolation method (or Fourier collocation approach, see [4]). For simplicity we treat the two-dimensional case here, i.e., n = 2, the extension to the three-dimensional case is straightforward, but requires more notational inconvenience. The images shown in the introduction, Figure 1, have been calculated by extending the following approach generically. We consider the domain Ω = [0, L 1 ] × [0, L 2 ] and N 1 × N 2 equidistant grid points in x 1 and x 2 directions, respectively. Then we use the two-dimensional discrete Fourier transform (DFT)û and N = N 1 N 2 is the overall number of grid points. The inverse is Note that the weights can be defined differently, as long as the concatenation of both transforms yields the identity. Replacing Now if u is sufficiently smooth, the derivatives of I n u i are spectrally accurate, as shown for example in the book by Canuto et al. [4]. The wave-number pairs k have to be chosen in a suitable set of wave numbers K that has to be ordered correspondingly to the implementation of the DFT/FFT we used. We work with the set The above interpolation operator I n allows for a simple differentiation procedure that gives spectrally accurate derivative approximations. In particular we can calculate the Laplacian as which becomes, as one can see, a multiplication with squared wave-numbers in the transformed space. We use the approximation I n u i (x) instead of the function itself and for the nonlinear contributions we interpolate the polynomial nonlinearities as follows Here we used the transforms for the nonlinear parts such that Inserting these expressions into the grain growth system (1.4), we derive that for each i ∈ {1, . . . , n g }. This yields the following ODE system for the coefficients, Treating the nonlinear parts explicitly and the linear ones implicitly allows easily to obtain a numerical approximation to the coefficients. The matrices of coefficients are typically arranged in vectors for theoretical aspects of the numerical schemes, i.e.û i (k, l) ↔ U i (p), so that the linear operator L becomes a diagonal matrix. With this notation the linear and nonlinear parts are and respectively, for i = 1, . . . , n g . Similarly as the Euler/Euler-2-step method one can define a SBDF1/SBDF1-2-step variant. Therefore we use superscripts for the indication of the time level we are at and calculate two updates, one for the vectors U n+1 Under the assumption that we start with a correct value, for a single PDE for u ↔ U, V, one can calculate the residual approximation With a tolerance one then proceeds as follows: If R ≤ update via U + = 2V n+1 − U n+1 , else, R > , repeat time-step starting from U n with smaller dt. In both cases we use the typical stepsize update dt ← ν R dt .
Here ν < 1 is a safety factor which we set to ν = 0.9. By employing above update the truncation error is reduced by one order. While u(t + dt) = V n+1 + O(dt 2 ), the new update corresponds to u(t + dt) = 2V n+1 − U n+1 + O(dt 3 ). However, we deal with a system of PDEs. We carry the same approach over for this case, only that now i and then we use the higher order update for each grain, Note that one could define U = (U 1 , . . . , U ng ) in one vector and the operators for U t as repetitions of (4.7) and (4.8). The two-step method applied to this system would just yield the presented individual updates. Figure 2 below shows a result calculated with the above approach, here employing an L 1 approximation of the norm for the residual calculation (via Riemann sums), for N 1 = 256, N 2 = 128 grid points, ν i = 1, ω = 1, λ = 2, n g = 60.
To visualize individual grain orientations with different shades of gray, we introduced the following visualization function Ψ(u i ) = i (log(i)u i ) 2 . After a transient time where the nuclei grow (given as spatially randomly distributed Gaussians in the initial data) and where the time steps are small, as topological changes can be fast, they become much larger after full re-crystallization. Then a comparably slow evolution phase, the coarsening of the grains, takes place. Until a first grain collapses more time can pass as for the initial re-crystallization transient. As for the last moments of a collapse the evolution becomes faster, each time the time-step automatically becomes small to capture the rate change. This is well reflected by the many spikes. Using these one could easily count the number of grains that vanished. In the evolution depicted in Figure 2 we are able to calculate for very long times, ending with a state with five grains only.

Discrete energy descent
In the following section we close this work by showing that our global interpolation method yields a discrete energy descent, theoretically in a semi-discrete approach, and practically by plotting a discrete energy in each time point of one particular simulation run. This shows that the numerics conserves the important property from the physics, descent of the free energy.
We consider a semi-discrete approach, similarly as by Ye [28] in a work on the Cahn-Hilliard equation, to show descent of a discrete energy. Ignoring the treatment of the temporal part of the PDE under consideration, our numerical method aims on finding a set of functions Here the discrete spatial grid for Ω = [0, L 1 ] × [0, L 2 ] is denoted by X = (x 1 j 1 , x 2 j 2 ) : x l j l = L l N l j l , j l = 0, 1 . . . , N l − 1, l = 1, 2 .
Algebraically this system is indeed equivalent to (4.5) with the corresponding initial conditions. In terms of the scalar product (·, ·) n defined as the semi-discrete formulation is also equivalent to (∂ t u i , v) n = ν i (∆u i , v) n + (I n φ i (u), v) n , ∀ v ∈ S n u i (x, 0) = u i 0 (x), x ∈ X , i = 1, . . . , n g .
This gives a generic tool for the numerical analysis as one gains the equality (u, v) n = (u, v) L 2 for all u, v ∈ S n . We can replace the scalar products by the L 2 versions. With the corresponding norm · n we define the discrete analogon of the natural free energy (1.2) of the system as follows with the nonlinearity Then we can see that

Proposition 15
The discrete version of the free energy dissipates as its continuous counterpart, namely, d dt E(u) ≤ 0 . (4.10) Proof By a direct calculation, we have which yields our conclusion (4.10).
For the last result we check by example that the descent (4.10) is valid for the adaptive time stepping scheme presented in Section 4.1. Furthermore, also by example, we pick up the convergence rate estimate derived in (3.23) and show that it corresponds well to our numerical findings. In fact, it describes the coarsening behavior even qualitatively. Figure 3 below shows two characteristic quantities during the re-crystallization process and the subsequent coarsening for the same parameter settings as in the previous subsection. Figure 3(a) shows how that the discrete energy (4.9) indeed decreases monotonically with each iteration. After full re-crystallization, the decrease becomes much slower. We are content with this observed quality of our scheme. Figure 3(b) shows the approximated H 1 -norm of the difference between u and someũ ∞ , which is assumed to be the state with only one existing grain (chosen randomly) and all other order parameters being zero. This is not exact, but sufficiently close to the real numerical u ∞ to analyze the convergence rate. The added curve, with a 'good guess' of θ = 1 5 , the function 9000(1 + t) −θ 1−2θ , shows the good quality of result (3.23). Note that the difference grows in the first stages as the order parameters are zero except of small regions with the nuclei. Although initially the H 1 distance to the equilibrium stateũ ∞ is small it grows in the early stage of evolution. In this way the energy is strongly decreasing as the minima of the multiple-well i Ψ i are adopted. Only later the H 1 distance as predicted in the asymptotic limit.