Convergence of the finite element method applied to an anisotropic phase-field model

We formulate a finite element method for the computation of solutions to an anisotropic phase-field model for a binary alloy. Convergence is proved in the H1-norm. The convergence result holds for anisotropy below a certain threshold value. We present some numerical experiments verifying the theoretical results. For anisotropy below the threshold value we observe optimal order convergence, whereas in the case where the anisotropy is strong the numerical solution to the phase-field equation does not converge.


Introduction
In this paper we study a finite element method for the numerical computation of an anisotropic phase-field model for a binary alloy.For details on modelling and physical background to this model we refer to Warren and Boettinger [12] and Kessler et al [7].Existence for the anisotropic model was proved by Burman and Rappaz in [1] and convergence of a finite element method in the isotropic case was proved in Kessler and Scheid [8].Other work on convergence of the finite element method for isotropic phase-field models include Chen and Hoffman [3], Feng and Prohl [5], Chen et al. [2].The anisotropy (as introduced in Kobayashi [9]) permits the modelling of branches in models of dendritic growth but makes the second order operator strongly nonlinear.However we show that this operator is strongly monotone, under a certain convexity condition, and Lipschitz continuous.The convergence of the finite element method is proved under regularity assumptions close to the regularity proved for the isotropic model.
We consider a binary alloy of two pure elements in both liquid and solid states inside a lipschitzian domain Ω ⊂ R 2 .The system is characterized by a relative concentration c = c(x, t), where the value c = 1 corresponds to the situation with only one element present and c = 0 with only the other, and by an order parameter φ = φ(x, t) (the phase-field), which takes values between 0 and 1.The value φ = 0 corresponds to a solid region and the value φ = 1 to a liquid region.The nonlinear parabolic system then takes the following form (see Burman and Rappaz [1] for details) ∂φ ∂t − div(A(∇φ)∇φ) − S(c, φ) = 0 in Ω × (0, +∞), (1.1) A(∇φ)∇φ • n = 0 on ∂Ω × (0, +∞), ( φ(0) = φ 0 , c(0) = c 0 in Ω, (1.5) where ∂Ω is the boundary of Ω and n is the unit normal to ∂Ω.We define the anisotropy matrix A(ξ) for ξ ∈ R 2 by where θ ξ denotes the angle between the x-axis and the vector ξ, the function a(θ) is given by a(θ) = 1 + ā cos(kθ) with k > 1 an integer corresponding to the number of branching directions.The functions S, D 1 and D 2 appearing in (1.1)-(1.5)are Lipschitz continuous functions, with first derivatives with respect φ and c uniformly bounded, satisfying S(c, 0) = S(c, 1) = 0, 0 < D s ≤ D 1 (φ) < D l and D 2 (c, φ) = 0 for c = 0 and 1 and φ = 0 and 1.In practice for the numerical computations we choose the following form of the non-linear functions (see Kessler et al [7].) where for 0 ≤ φ, c ≤ c.Outside this interval all these functions are extended continuously by a constant.We remark that by integrating equation (1.2) on Ω and by using (1.4) we obtain conservation of mass In Burman and Rappaz [1], we proved existence of a weak solution for this strongly nonlinear parabolic system under certain assumptions on the parameter ā.To be more specific we denote by V = H 1 (Ω), V the dual space of V , T the final time.The L 2 -scalar product is denoted (•, •) Ω and the corresponding norm • Ω .We have proved for all v, w ∈ H 1 (Ω) and a.e. in (0,T).Furthermore, if where Q T = Ω × (0, T ).Moreover, if we extend the mappings S and D 2 by zero outside the unit square (0, 1) × (0, 1) and if we assume

A finite element method
We discretize the above system of equations using P 1 -lagrangian finite elements in space and a semi-implicit Euler-scheme in time.Let T be a triangulation of Ω.For any triangle K ∈ T , we denote by h K its diameter and set h = max K∈T h K .Let V h be the finite element space defined by where P 1 (K) denotes the set of polynomials of degree 1 on K.For an integer N > 0 we introduce τ = T /N and t n = nτ , n = 0, 1, 2, . . .We consider the following fully discrete scheme.Given for all and n = 0, 1, 2, 3, . . .Note that to compute (φ n+1 h , c n+1 h ) we only need to solve a nonlinear system for φ n+1 h .We can prove in the same way as in Ref. 1 that given (φ n , c n ) both equations have a unique solution [1] (φ n+1 h , c n+1 h ) for timesteps τ sufficiently small.Existence of φ n+1 h is proved using direct methods in the calculus of variations (see Dacorogna [4]) and the existence of c n+1 h follows by a standard application of the Lax-Milgram lemma.

The nonlinear operator
To prove the convergence of the finite element scheme we need to extend the analysis concerning boundedness and continuity of the nonlinear operator.First we recall some fundamental results [1], which we state here without proofs.

Lemma 3.1: When
the Ginzburg-Landau potential

Lemma 3.2:
The Gateaux derivative of the Ginzburg-Landau potential, exists for each φ ∈ V and is given by Lemma 3.3: The anisotropic operator satisfies the following upper and lower bounds: where We also recall some results on Eulerian operators derived from Gateaux-differentiable functionals.

Definition 3.4:
where (•, •) denotes the duality pairing between V and V .In addition to these results we need the following lemma, stating that the nonlinear operator defined by (3.3) is strongly monotone and Lipschitz continuous with respect to the H 1 (Ω)-seminorm, | • | V .Lemma 3.5: If the convexity condition ā < 1 k 2 −1 holds, then we have the following inequalities for all φ, ψ ∈ V : ) where L is a constant independent of φ, ψ and • Ω is the L 2 -norm of vectorial functions.
Proof: The first inequality is a consequence of the fact that the Eulerian operator derived from a convex functional is monotone.We consider the perturbed Ginzburg-Landau functional It is easy to see that for some sufficiently small µ ā > 0 this functional will remain convex.We consider the corresponding Hessian matrix of ξ → a(ξ) 2  2 − µā 2 |ξ| 2 in polar coordinates, for ξ = 0: where O θ denotes the matrix of rotation of the angle θ.It is easy to show [1] that if ā < (k 2 − 1) −1 , then a(θ) + a (θ) > 0 and H(ξ) will remain positive definite for ξ when .
Hence, by the monotonicity of the corresponding Eulerian operator we have where • denotes the scalar product in R 2 , η T is η-transposed, and H is the above hessian matrix with µ ā = 0. Since the spectral norm of H is bounded independently of ξ, we easily obtain inequality (3.5) with The convexity of the functional is of essential importance for the wellposedness of the system.To illustrate how convexity is lost we plot the contourlines of the integrand of the functional (3.2) with k = 4 for the case ā = 0.05 and ā = 0.15 in Fig. 1.Note the non-convex zones appearing around θ = nπ/2, n = 0, 1, 2, 3 when the anisotropy parameter is larger than 1/15.These non-convex zones corresponds to "forbidden" gradient directions and will give rise to corners and rapidly oscillating gradients in the finite element approximation (2.1) of the equation for φ as we will see in the numerical section.

Regularity hypothesis
The strong non-linearity in the anisotropic operator makes a priori estimates on higher order derivatives very difficult to prove, especially considering that A(ξ)ξ is not differentiable at ξ = 0.For the isotropic problem on the other hand quite extensive regularity results were proved in Rappaz and Scheid [10].In fact we have provided that the initial data are sufficiently regular, that is to say φ 0 ∈ H 2 (Ω), ∂φ0 ∂n = 0 on ∂Ω and c 0 ∈ H 1 (Ω).This however does not suffice to show convergence of the finite element method.In the sequel we will assume that there exists a unique solution (φ, c) of the system (1.13) such that both φ and c enjoy the regularity proved for φ and in addition that the gradients are bounded on the space time interval.So we will suppose that (φ, c) ∈ W where This assumption is reasonable as long as the anisotropic functional remains strictly convex such that the strong monotonicity (3.4) holds.To make these assumptions sufficient for the convergence proof to hold we still need to show that this implies sufficient regularity of the time derivatives.For this we need to assume that the non-linear terms S, D 1 , D 2 have bounded first derivatives in φ and c.We show formally in the following lemma that this implies the necessary regularity of the time derivatives.

Lemma 4.1: Under the above regularity hypothesis we have
Proof: We only give the proof for the strongly non-linear equation for φ since the proof for the concentration is similar.Formally differentiating equation (1.1) with respect to t we get We study this equation on weak form: Clearly we have A straightforward calculation using the definition of the non-linear operator and the relation where from which we conclude that Taking the square of this inequality and integrating in time yields (4.1) for φ.

Convergence of the finite element method and error estimate
Proof: In the sequel L A , L S , L D1 and L D2 will denote the Lipschitz constant associated with operators A(∇φ)∇φ, S(φ, c), D 1 (φ) and D 2 (c, φ) furthermore D max and the finite element formulation (2.1), we may write Now using the notation φ n ∆ = φ n h − φ(t n ) and the equality (5.2) we obtain for all It now follows by Lemma 3.5 that Now multiplying by τ and summing over n we obtain, using summation by parts, We proceed by adding and subtracting φ(t n+1 ) in the two last terms in the right hand side and using the Cauchy-Schwarz inequality in combination with Young's inequality We now eliminate the second term on each side of the inequality, we use the Lipschitz continuity of the source terms and a duality argument for the second derivative in time to obtain (using Cauchy-Schwarz and Young's inequality repeatedly) and (5.7) Now using (5.6) and (5.7) in (5.5) and collecting terms we obtain where c n ∆ = c n h − c(t n ).We turn to the equation for the concentration c and obtain, in the same fashion (5.9) We use the formulation (2.1) to replace the c n h in c n ∆ by some arbitrary function Proceeding as above we may write for terms I 1 and I 5 as and using the Young inequality

Convergence of the finite element method
We treat I 2 to I 4 in the same spirit as (5.6) leading to and analogously and Multiplying by τ and summing over n in equation ( 5.10) we have Finally we multiply (5.12) by η = µāDs 64D max 2 , add (5.8) and (5.12) and apply the discrete Gronwall lemma to obtain Since in particular our regularity assumptions include we may chose where π h denotes the interpolation operator.The theorem now follows by a standard interpolation estimate.
Remark 5.2: Note that the exponential factor α is of the order of ∇φ 2 which under the regularity hypothesis should be of the order δ −2 if δ denotes the interface thickness.This is the typical worst case estimate for phase-field equations (see for instance Kessler and Scheid [8] or Chen and Hoffman [3].) However in a recent paper Feng and Prohl [5] show that for the isotropic, thermal phase-field equation, an estimation of the smallest eigenvalue of the linearized operator permits a priori estimates which show growth only in low polynomial order of δ −1 provided that all interior layers are developed in the initial data.

Numerical tests
Implementation of the numerical scheme (2.1) was done using the finite element package ALBERT developed by Schmidt and Siebert [11].We have set up tests to obtain the experimental numerical convergence order of the scheme in the norm L 2 (0, T ; H 1 (Ω)) and compare it with the theoretical result of Theorem 5.1.We have also measured experimental orders of convergence in the L ∞ (0, T ; L 2 (Ω)) norm.Tests have been run using both low and high anisotropy.

Implementation of numerical tests
For the tests, parameters of the nonlinear functions defined in (1.6)-(1.12)are set to We have treated the nonlinearities of numerical scheme (2.1) with just one step of a fixed-point method.Quadrature is exact for polynomials of degree 3, except for the mass matrices for which we used mass lumping.
By adding extra artificial source terms, we have imposed exact solutions to both equations, with which to compare the numerical solutions.We chose solutions that reproduce some features expected of the solutions of system (1.1)-(1.5).Namely, accross the solid-liquid interface, the phase-field is known to have a hyperbolic-tangent-like profile while its values change from 0 to 1, while the concentration goes smoothly from values close to a small constant c s in the solid region, to a large constant c l on the liquid side of the interface, and then down to an intermediate value c 0 in the liquid bulk phase.[6] Also, we assume that propagation velocities of the interface vary depending on the interface's normal direction proportionally to the anisotropy function a(θ).
Since we would like to see all interface directions in our test, we choose exact solutions whose initial conditions represent a circular interface separating bulk solid and liquid phases with different concentrations.The transition is smooth as described above.The system then evolves, and the interface moves outward, with local velocities depending on the interface's normal direction, assimilated in the definition of the test solutions to the angle in local polar coordinates.
Let us define polar coordinates associated to position x by where ρ = x and θ are the cylindrical coordinates of x.In the sequel, we will be always working in the space-time domain [0, 1] 3 .
Let ρ 0 and v be given constants representing the radius of the initial circular interface and the solidification front velocity.We define the auxiliary function and the actual imposed phase-field solution Let c 0 , c s , cl and ρ ∆ be given constants representing the values of the concentration in the liquid and solid bulks, a value proportional to the concentration on the liquid side of the interface, and an δ-rescaled shift from the center of the interface.We then define the two auxiliary functions and the imposed concentration solution For the tests, we fix the previous numerical constants to ρ 0 = 0.2, v = 0.6, c 0 = 0.4, c s = 0.2, cl = 0.8 and ρ ∆ = 1.
As an illustration, radial profiles of the imposed solutions (6.8) and (6.10) for t = 0.5 and θ = 0 are shown in Figure 2, as well as level sets of φ for both low and high anisotropy at the final time t = 1 in Figure 3.
In the implementation, exact solutions (6.8) and (6.10) are used as initial and Dirichlet boundary conditions, and their derivatives are combined to define artificial source terms added to both equations, ensuring that they are then solutions of the differential system.

Results of numerical tests for low anisotropy
We now present numerical results for a low anisotropy ā = 0.05, which is in the scope of the theory presented in this paper.We performed two series of tests: one in which the timestep size was decreased linearly with the space mesh size, and another where the timestep size was decreased quadratically with the mesh size.We are interested in the experimental orders of convergence (OC) for the L 2 (0, T ; H 1 (Ω)) norm of the error, which is in the scope of the theory, and also the L ∞ (0, T ; L 2 (Ω)) norm, for which the convergence rate predicted by our theoretical result is expected to be suboptimal.From the results in Tables 1-4, we verify experimentally that the order of convergence h + τ predicted by Theorem 5.1 for e L 2 (0,T ;H 1 (Ω)) is optimal.However, we can also conclude that this order of convergence is suboptimal for e L ∞ (0,T ;L 2 (Ω)) , which converges faster, at a rate h 2 +τ .The experimental orders of convergence are also illustrated by Figures 4-5, made in log-log scale using the data from Tables 1-4.

Results of numerical tests for high anisotropy
For the high anisotropy, we take ā = 0.10.We observed in numerical tests with imposed solutions (6.8) and (6.10) that the numerical solution is unable to reproduce the features of a presumably H 1 regular solution in regions where the normal direction to the level sets of solutions is pointing at an angle close to 0, π/2, π or 3π/2.This corresponds to nonconvex portions of the Franck diagram, i.e. the graph of a level set of the anisotropy energy  as a function of ∇φ.These are the angles that physicists call "forbidden angles".Qualitatively, it seems that level sets of the numerical solution avoid "forbidden angles" in the imposed solution by zigzagging at "permitted angles", much like a sailboat would zigzag in an effort to sail upwind, using only directions in which it is possible to sail.
We have also tried imposing different solutions, planar front equivalents of (6.8) and (6.10), which contain only one direction.We have chosen a forbidden direction (θ = 0) and a permitted direction (θ = π/4) as examples.The corresponding imposed solutions are the same as (6.8) and (6.10), with a change of definition for ρ and θ, which are no longer the polar coordinates.The angle θ is instead the constant 0 or π/4, whereas ρ is defined as respectively x 0 or (x 0 + x 1 )/ √ 2. The qualitative behavior of these solutions can be seen in Figure 6.We now present numerical evidence of convergence in the L 2 norm even for the "forbidden direction" test, and in the H 1 norm only in the case of a "permitted direction".Apparently, the zigzagging behaviour still allows the solid-liquid front to evolve with a correct average velocity, but with wrong local gradients.For the sake of brevity, in this section we present only results for φ, and decreasing the timestep quadratically with the mesh size.We have observed that in the high anisotropy regime, c always converges better than φ, and its level sets never present the zigzagging behavior.From the results presented in Tables 5-6 and Figure 7, we conjecture that the result of Theorem 5.1 still holds in the high anisotropy case, whenever the level sets of the solution have normals in the region where the anisotropy operator is still convex.

Figure 6 :
Figure 6: Level sets of φ for forced solutions with high anisotropy