# A new level set-finite element formulation for anisotropic grain boundary migration

Grain growth in polycrystals is one of the principal mechanisms that take place during heat treatment of metallic components. This work treats an aspect of the anisotropic grain growth problem. By applying the first principles of thermodynamics and mechanics, an expression for the velocity field of a migrating grain boundary with an inclination dependent energy density is expressed. This result is used to generate the first, to the authors' knowledge, analytical solution (for both the form and kinetics) to an anisotropic boundary configuration. This new benchmark is simulated in order to explore the convergence properties of the proposed level set finite element numerical model in an anisotropic setting. Convergence of the method being determined, another configuration, using a more general grain boundary energy density, is investigated in order to show the added value of the new formulation.

## Authors

• 1 publication
• 1 publication
• 1 publication
• 1 publication
• 1 publication
01/19/2020

### A Ritz-based Finite Element Method for a Fractional-Order Boundary Value Problem of Nonlocal Elasticity

We present the analytical formulation and the finite element solution of...
01/02/2020

### Numerical investigation into fracture resistance of bone following adaptation

Bone adapts in response to its mechanical environment. This evolution of...
06/09/2021

### A 2D front-tracking Lagrangian model for the modeling of anisotropic grain growth

Grain growth is a well-known and complex phenomenon occurring during ann...
06/03/2021

### Comparative study and limits of different level-set formulations for the modeling of anisotropic grain growth

Four different finite element level-set (FE-LS) formulations are compare...
04/16/2021

### Coupling of complex function theory and finite element method for crack propagation through energetic formulation: conformal mapping approach and reduction to a Riemann-Hilbert

In this paper we present a theoretical background of a coupled analytica...
09/17/2020

### A new front-tracking Lagrangian model for the modeling of dynamic and post-dynamic recrystallization

A new method for the simulation of evolving multi-domains problems has b...
10/25/2019

### Towards understanding the boundary propagation speeds in tumor growth models

At the continuous level, we consider two types of tumor growth models: t...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

During metal forming operations the microstructures of metallic components are modified by a host of phenomena ranging from solid state transformations to recrystallization [humphreys2012recrystallization]. Perhaps the most critical effect of these mechanisms is the grain boundary motion they induce. Since the in-service properties of metal pieces depend on the material’s microstructural characteristics (grain size, crystal orientation, composition, etc…) [suttonballuf2006], it is important to study how these boundaries evolve under thermo-mechanical loads. Crystalline interfaces migrate differently during the different stages of annealing [humphreys2012recrystallization]: deformation, recovery, recrystallization and grain growth can take place. These dynamics have been widely studied experimentally and numerically. Even so, the long characteristic times associated with grain growth allow investigators to decouple its effects from other processes. This is most likely the reason for which the theory of grain growth is the most established in monographs on the subject.

Grain boundary motion during grain growth is thought to be driven by the reduction of the interfacial free energy [herring1999surface]. Classical models for grain growth in polycrystals use homogenized grain boundary properties to describe crystal interfaces (i.e. constant energy density, constant mobility, etc…) [anderson1984computer, gao1996real, lazar2011more, bernacki2011level, garcke1999multiphase]

. However, at the mesoscopic scale, the grain boundary can be parameterized by five macroscopic crystalline parameters: a boundary plane unit normal vector and a misorientation element

[suttonballuf2006]. The main challenge in the current study of grain boundary motion is the dependence of intrinsic properties such as the grain boundary energy and mobility on these multiple structural parameters. Notwithstanding the difficulty in determining both the energy and mobility of a crystalline interface experimentally [Holm01, Rohrer2010, Adams1997, Morawiec2000, Saylor2003] or numerically [OlmstedI2009, OlmstedII2009], the grain boundary configuration space is itself highly non-Euclidean. Indeed, current work towards elucidating the structure of this space has been oriented towards defining five dimensional identification spaces [olmsted2009new] as well as defining higher dimensional embeddings into the unit octonions for example [francis2019geodesic].

In order to better predict microstructural evolution using numerical models, the intrinsic properties of the crystalline interface must be taken into account. However, many nuances of models taking into account different aspects of boundary variability can be developed. The authors have chosen to differentiate three classes of models and will refer to them using the terminology: isotropic, heterogeneous and anisotropic. In isotropic models, boundary properties are defined as constants for the entire system. While being able to reproduce mean value evolutions, such as mean grain size or even grain size distributions, rather efficiently, local heterogeneities in microstructures, such as the twin boundary, can not be modeled correctly [anderson1984computer, gao1996real, lazar2011more, bernacki2011level, garcke1999multiphase]. Heterogeneous models may employ homogenized intrinsic properties along each grain boundary but differentiate different boundaries between each other [rollett1989simulation, hwang1998simulation, upmanyu2002boundary, Fausty2018, zollner2019texture, kazaryan2002grain, miyoshi2016validation, chang2019effect, miyoshi2019accuracy, miessen2015advanced, Fausty2020]. As such, in a polycrystalline setting, the misorientation dependence of boundary properties can be modeled by these methods but not the inclination dependence. Fully anisotropic approaches attempt to take into account the five parameter dependence of grain boundary properties and as such constitute the most general of the three types.

In the journey towards these fully anisotropic models, able to account for general energy densities and mobilities, this work is constrained to treating the anisotropic grain boundary energy density in a one boundary setting. As such, the energy density function will be inclination dependent and the mobility will remain constant. This constraint on the energy density is rather easily scaled up to the polycrystalline case since a grain boundary, by definition, has a constant misorientation in fully recrystallized microstructures and thus any property variation can be fully attributed to changes in the boundary plane. As such, the developments here can be readily integrated into heterogeneous models, such as [Fausty2018, chang2019effect], in order to model anisotropic grain boundary energy densities. The homogeneity of the mobility is another matter however. Very few investigations deal with fully anisotropic mobilities mostly because the mobility of the boundary is really a derived notion which does not have a clear definition. Indeed, the grain boundary mobility, contrary to the thermodynamic definition of the energy density, is a kinetic parameter which is often fitted to produce observed migration rates. In this work a global grain boundary mobility definition will be proposed in relation to a normalized rate of energy dissipation. However, the question of using fully anisotropic and possibly tensorial values for the grain boundary mobility is still an open one.

As such, this manuscript, starting from a differential geometry description of the grain boundary and thermo-dynamical first principles, proposes an expression for the velocity of a migrating grain boundary in a general anisotropic energy density setting. This velocity is then applied to the transport equation of the level set method. In order to test this new formulation, an analytical benchmark using a collapsing ellipse is developed. Subsequently, a finite element level set numerical model is proposed to simulate this analytical test case. This new numerical model is then used to compute dynamics in a more general boundary energy density configuration.

## 2 The interface model

The model presented here is developed using elements from the field of differential geometry [lee2003graduate, spivak2005comprehensive]. While more rudimentary mathematics can treat problems dealing with surfaces, the language of differential geometry seems like the correct one to treat the anisotropic problem. The necessary mathematical tools are briefly introduced before developing the formalism itself. While these tools might be well-known by experts, the readership specializing in physical metallurgy will, most likely, not be familiar with the terminology. For this reason, a consequential section of the text is devoted to definitions of known mathematical concepts.

### 2.1 Tools and Definitions

###### Definition 2.1.

A smooth -manifold is a triplet comprised of a set , a given topology on the set and a smooth atlas made up of smooth charts.

The manifold description of interfaces is chosen because it is technically the minimal structure with which one must endow a space in order to define derivatives. One could of course weaken the smooth condition to a or perhaps even constraint, however, for the sake of simplicity, smoothness () is considered here.

###### Notation.

Let is the set of all smooth functions that can be defined on the smooth manifold .

For the following definitions let be a smooth manifold.

###### Definition 2.2.

The tangent space at the point is the vector space comprised of elements such that there exists a smooth curve of

 C:R→Mt↦C(t)

with and

 X:C∞(M)→Rf↦Xf:=ddt(f∘C)(0)

The elements of the tangent space, , are also called tangent vectors. Indeed, the elements of the tangent space to a point in the manifold and the classical notion of tangent vectors in Euclidean space are related. As an example using these definitions, if one chooses a chart such that and a function then and element acts on through its equivalent curve

 Xf=ddt(f∘C)=ddt(f∘x−1∘x∘C)

which, using the multidimensional chain rule brings one too

 Xf=ddt(xi∘C)∂i(f∘x−1)

where is the th component function of the chart , is the derivative operator of a multidimensional function with respect to its th component and the Einstein summation convention is in effect, which will be implied from here on unless stated otherwise.

###### Theorem 2.1.

One may construct an orthonormal basis for with the vectors defined as

 ∂∂xif:=∂i(f∘x−1) (1)

As such, for any element of , one may define its components in this basis, using

 Xi=ddt(xi∘C)(0) (2)

where is the curve associated with . As such, its decomposition is written

 X=Xi∂∂xi (3)

Seeing as is a vector space, it admits a dual space.

###### Definition 2.3.

The dual vector space or co-tangent space to the tangent space is the space of linear maps such that

 ω:TpM→RX↦ω(X)

More generally, the local tensor spaces are constructed from the tangent space and its dual.

###### Definition 2.4.

The space of -tensors, , at is defined as

 Tqp,sM:=T∗pM⊗⋯⊗T∗pMs⊗TpM⊗⋯⊗TpMq

where is the tensor product of spaces.

From the tangent spaces at each point of the tangent bundle can be constructed.

###### Definition 2.5.

Let be defined as

such that the tangent bundle is defined as

 TMπ⟶M

where is a continuous surjective map.

Analogously, the -tensor bundles are defined in the same manner.

###### Definition 2.6.

A section of a bundle is a continuous map such that

 σ:B→Eπ(σ(p))=p

Colloquially, the sections of the tangent bundle are called vector fields and in the same manner sections of the tangent bundles are called tensor fields.

###### Notation.

is the space of all smooth sections of the bundle .

###### Definition 2.7.

A Riemannian -manifold is a smooth -manifold equipped with a symmetric -tensor field , called a metric, such that is a positive-definite tensor.

The positive definiteness of means that for any

 g(p)(X,X)>0

.

Riemannian manifolds are of general interest since the metric structure defines inner products on the tangent spaces. As such, a Riemannian manifold is convenient for defining lengths of curves and more general measures of volume. Indeed, this metric structure is what allows one to define the Riemannian integral on the manifold.

###### Definition 2.8.

A differential -form on a smooth manifold is a completely anti-symmetric -tensor field.

###### Corollary 2.8.1.

The volume form of an oriented Riemmannian -manifold is the differential -form such that for a given chart the volume form may be expressed as

 dM=√det(g)dx1∧⋯∧dxn

where is the determinant of the matrix composed by the components of in the chart , is the dual basis of the co-vector space and is the exterior product of differential forms.

Using this machinery, any function can be integrated over the manifold. One is also able to define a relatively straightforward connection on the space called the Levi-Civita connection.

###### Definition 2.9.

A connection over a bundle is a set of linear maps

 ∇: Γ(TqsB)→Γ(TqsB⊗T∗B)

that respect the Leibniz rule,

 ∇(fσ)=σ⊗df+f∇σ (4) ∇(τ⊗σ)=∇τ⊗σ+τ⊗∇σ (5)

where is the classic differential of a smooth function .

###### Corollary 2.9.1.

From a given connection , one may construct a covariant derivative

 ∇⋅: Γ(TB)×Γ(TqsB)→Γ(TqqB) ∇⋅(X,σ)=∇Xσ=(∇σ)(X)

where, when working in a chart, one may use

 (∇Xσ)k…j…=(∇σ)k…j…iXi=∇iσk…j…Xi
###### Definition 2.10.

The Levi-Civita connection on a Riemannian manifold is the unique connection on the tensor bundles which satisfies

 ∇g=0

and has no torsion.

### 2.2 An energetic embedded smooth manifold

Let be a Riemannian -manifold with metric and is a smooth -manifold with . Let be a smooth embedding from to

 φ:S→MS≡homeoφ(S) (6)

where describes a homeomorphism equivalence and Figure 1 provides an illustration. The embedding also provides a map from the tangent bundle of to the tangent bundle of called the push-forward.

###### Definition 2.11.

The push-forward of a map from to , two smooth manifolds, is the linear map such that

for

Much in the same manner, the embedding generates a map from the co-tangent bundles to .

###### Definition 2.12.

The pull-back of a map from to , two smooth manifolds, is the linear map such that

 φ∗: T0qM|φ(S)→T0qS (φ(p),σ)↦(p,φ∗σ) (φ∗σ)(X(1),…,X(q))=σ(φ∗X(1),…,φ∗X(q))

Using the charts and and using the convention by which objects in are indexed by Greek letters and objects in are indexed by Latin letters one can express the components of the pushforward of a vector using its action on a function

 (φ∗X)f =(φ∗X)α∂f∂Zα

and

 (φ∗X)f =Xi∂(Zα∘φ)∂xi∂f∂Zα

which, defining

 φα: S→R p→Zα(φ(p))

 (φ∗X)α=Xi∂φα∂xi (7)

Using the pull-back one may define the metric on and therefore turn into a Riemannian manifold with

 g(p)=(φ∗m)(φ(p)) (8)

which, using the same charts as above and two vectors

 g(X,Y) =(φ∗m)(X,Y) gijXiYj =m(φ∗X,φ∗Y) =mαβ(φ∗X)α(φ∗Y)β =mαβ∂φα∂xi∂φβ∂xjXiYj.

This results defines the components of the induced metric by identification

 gij=mαβ∂φα∂xi∂φβ∂xj (9)

Now, let be an internal property space. For example, , when applied to grain boundaries, would be the five-dimensional space created by the misorientation and inclination parameters . Let

 SB=⋃p∈S(p,B)=S×B (10)

and define the trivial property bundle

 SBπB⟶S (11)

such that a section of the property bundle describes exactly the properties of the -manifold at each point. If one was to define an energy density map

 γ: B→R+

then one could calculate the energy density at any point through the property field as . Given that is now a Riemannian manifold, this energy density can be integrated in order to give the total interface energy of the embedding as

 I=∫S(γ∘b)dS

The model developed here for the interface is thus a triple from which, with an energy density map , the total energy of the interface may be expressed. By design, this model puts no lower bound on . Therefore, this structural model is readily generalized to objects that are not strictly interfaces but can be of lower dimension, such as lines if . This is an important aspect of this model, even if it might be out of the scope of this article, if ever one was to attempt to attribute properties and therefore energies to other defects in the polycrystal microstructure.

### 2.3 Interface thermodynamics

Consider a closed thermodynamic system made up of a Riemannian -manifold of volume , an embedded interface with a boundary energy density map , in a heat bath at absolute temperature , a system entropy and a homogeneous pressure field . One defines an idealized case with the following conditions:

• isothermal heat treatment - and are constant

• the time is parameterized so that

• and the system is closed

the change in the internal energy during the free evolution of the system is defined as

 dUdt=dIdt+Tdηdt−pdVdt (12)

If one looks at the evolution of the Gibbs free energy:

 dGdt =dUdt+d(pV)dt−d(Tη)dt =dIdt+Vdpdt−ηdTdt,

which considering the isobaric and isothermal conditions,

 dGdt=dIdt (13)

such that the change in free energy is exactly equal to the change in interface energy. Additionally, according to the second law of thermodynamics, the closed system must tend to minimize its free energy such that

 dIdt≤0. (14)

The principle of least action affirms that the energy dissipation must be maximal and thus must be minimal .

The flow of the interface, , is defined as

 ψ: S×[0;1]→M (p,t)↦ψ(p,t) ψ(p,0)=φ(p)

thus the embedding is defined.

As the interface evolves only its geometry changes. This means that the misorientation remains constant. Following this statement, the property field should depend, in some manner, on the embedding. If we consider that the boundary is parameterized by the misorientation-inclination pair, , the misorientation is invariant

 dMdt=0. (15)

However, since the inclination of the boundary is a geometrical characteristic of the boundary it does change. The field depends exclusively on the embedding and they are related by means of the push-forward of the tangent vectors to the interface. For any at any time the value of is in completely determined by

 {m(n,n)=1m(n,ψ∗X)=0

or, in component form,

 ⎧⎪⎨⎪⎩mαβnαnβ=1mαβnα∂ψβ∂xiXi=0

 mαβnα∂ψβ∂xi=0,∀i=1,…,s (16)

In the context of the one boundary, and thus one misorientation, the following simplification holds

 γ(M,n(…,∂ψα∂xi,…))=γ(…,∂ψα∂xi,…). (17)

For and with the velocity field is defined as

 (vf)(q) =ddt(f∘ψ(p,⋅))(t) =dψαdt(p,t)∂f∂Zα(q) =vα∂f∂Zα

such that, by identification

 vα(ψ(p,t))=dψαdt(p,t). (18)

Using the statement in corollary 2.8.1, equation (18) and considering the Levi-Civita connection of , knowing that , the energy dissipation may be defined as

 dIdt =∫Sddt(γdS) =∫S1√det(g)∂(γ√det(g))∂∇iψα∇ivαdS.

Expressing

 ∂(γ√det(g))∂∇iψα =√det(g)∂γ∂∇iψα+γ∂√det(g)∂∇iψα,

using Jacobi’s formula to define the derivative of the determinant of matrix and defining the components of an inverse metric tensor as , the energy dissipation may be rewritten as

 dIdt=∫S(∂γ∂∇iψα+γgiqmσα∇qψσ)∇ivαdS (19)

One can define the boundary of as and apply Stokes’ theorem such that

 dIdt= ∫∂Sgikτk(∂γ∂∇iψα+γgiqmσα∇qψσ)vαd∂S −∫S∇i(∂γ∂∇iψα+γgiqmσα∇qψσ)vαdS

where is the outside pointing unitary normal field to .

In order to encapsulate the quantities of interest, we define the following restricted vector fields, with components

 Bα =mαβgikτk(∂γ∂∇iψβ+γgiqmσβ∇qψσ) (20)

and with

 Aα=mαβ∇i(∂γ∂∇iψβ+γgiqmσβ∇qψσ) (21)

such that

 dIdt(t)=∫∂Sm(B,v))|ψ(p,t)d∂S−∫Sm(A,v)|ψ(p,t)dS (22)

Given the one interface restriction used in this work, the boundary of the interface can only be either empty or part of the boundary of the base manifold . Using the case, the boundary term disappears and

 dIdt=−∫Sm(A,v)dS, (23)

As such, the velocity of the boundary is that which minimizes the previous expression. The bilinear form

 ⟨⋅,⋅⟩: Γ(TM|ψ(S,t))2⟶R (24) (X,Y)↦∫Sm(X,Y)dS (25)

defines an inner product on the space, turning it into a Hilbert space. As such, using the Cauchy-Schwarz inequality, one may show that the velocity field that minimizes the energy dissipation has the form

 v=μA,dIdt=−μ⟨A,A⟩ (26)

with being classically the mobility of the boundary.

By replacing the interfacial energy dissipation term in equation (13) with its equivalent expression one determines a definition of the mobility parameter as

 μ=−1⟨A,A⟩dGdt. (27)

The mobility is thus proportional to a normalized value of the energy dissipation of the system. Thus, the mobility of the boundary appears as a kinetic parameter related to the capacity of the system to dissipate energy in the form of heat or work (by a contraction due to the excess volume of the boundaries for example). As such, the mobility of the grain boundaries may have more to do with the boundary conditions imposed on the system then previously imagined.

### 2.4 From embeddings to level-set fields

###### Definition 2.13.

A level-set map or function is a smooth scalar field over the smooth manifold such that, given an embedding

 ϕ(φ(p))=0 (28)

.

Most often, one defines the level-set function as a signed distance function to the interface such that with

 d: M×M→R+ (p,q)↦minC(p,q)∫C(p,q)dC

where is any curve from to , one may then fix

 ϕ(q)=±d(q,φ(S)):=±minp∈Sd(q,φ(p))

where one makes a choice of sign over the domains that the interface separates. The evolution of the interface is simulated by solving the transport equation

 ∂ϕ∂t+v⋅~∇ϕ=0 (29)

everywhere in , where is the Levi-Civita connection on the Riemannian manifold . One may replace the velocity field with the expression developed in the previous paragraphs after a suitable extension of the fields defined on to the entire manifold. As such,

 ∂ϕ∂t+μAα~∇αϕ=0, (30)

Using the following identity

 ∇j∇iφα~∇αϕ=−∇iφα∇jφβ~∇β~∇αϕ, (31)

derived from the orthogonality condition of tangent vectors and the gradient of the level-set, one may express the transport equation in a fully level-set form

 ∂ϕ∂t−μ(γmαβ+∂2γ∂~∇βϕ∂~∇αϕ)~∇α~∇βϕ=0 (32)

where the full derivation is reported in A. Of course, some direct analogies can be made with the derivations proposed in [herring1999surface].

### 2.5 Constraints on the anisotropic grain boundary energy density function

Let be the symmetrized tensor such that

 Dαβ:=γmαβ+12(∂2γ∂~∇βϕ∂~∇αϕ+∂2γ∂~∇αϕ∂~∇βϕ) (33)

where the tensor can be symmetrized because is already symmetric. Now, equation (32) is a purely diffusive equation with as a diffusive coefficient tensor. As such, the well-posedness of the problem depends largely on the positive definiteness of . For solutions to be unique, one must have

 D(ω,ω)>0 (34)

and, therefore,

 Dαβωαωβ>0. (35)

Given the arbitrariness of , applying (35) to the basis vectors of the dual tangent spaces at each point, one quickly obtains (not using the summation convention)

 Dααω2α>0 Dαα>0.

More complex conditions exist for the off-diagonal components. For example, for and , one can show that

 |D12|

 |D12| <√D11D22. (36)

As such, given that the tensor depends entirely on the grain boundary energy function , these constraints are actually directly transferable to the function. Thus, in order to preserve uniqueness of the grain boundary flow, the anisotropy of the function is restricted to maps which satisfy these conditions. While determining the space of functions that satisfy these relations would be a valuable discovery for the community, this kind of development is out of the scope of this article.

## 3 An elliptical benchmark

To the authors’ knowledge, no analytical test case exists for the anisotropic one boundary setting of grain growth. Indeed, while the shrinking sphere is a viable benchmark for the isotropic case and the “Grim Reaper” [Garcke1999] is very useful for testing heterogeneous models, no equivalent configuration has been developed for more general grain boundary energy densities. Theoretical studies have proven that minimal energy surfaces can be constructed for virtually any inclination dependent energy density function [HoffmanI1972] using Wulff shapes, the kinetics with which these shapes should evolve, in a isolated grain undergoing coarsening for example, are completely unknown. Thus, these semi-analytical benchmarks remain incomplete cases for numerical testing. This section is devoted to generating such a completely analytical solution to the problem with constrained kinetics as well as definite morphology.

### 3.1 The setting

Consider a circle as a smooth manifold with the circle topology and smooth structure and the Riemannian manifold equipped with the standard topology and differentiable structures and the flat metric . Using the chart and the Cartesian chart one may construct the following elliptical embedding

 φ: [0;2π]⟶R2 θ↦(acosθ,bsinθ)

where and Figure 2 illustrates this embedding.

All of the relevant geometrical information may thus be extracted from the embedding. The pushforward of the tangent space

 ∂φx∂θ =−asinθ (37) ∂φy∂θ =bcosθ (38)

and the induced metric tensor

 gij =mαβ∇iφα∇jφβ gθθ =(∇θφx)2+(∇θφy)2

such that

 gθθ=a2sin2θ+b2cos2θ (39)

The Levi-Civita connection is thus defined by

 ∇θgθθ =0 ∂gθθ∂θ−2Γθθθgθθ =0 Γθθθ =12gθθ∂gθθ∂θ

where are Christoffel symbols. Therefore,

 Γθθθ=(a2−b2)cosθsinθa2sin2θ+b2cos2θ (40)

### 3.2 A solution

Now consider the boundary energy

 γ(θ)=Gθθgθθ (41)

where is a -tensor field of whose only component is actually a constant in this chart. As such, using equations (26) and (21) the velocity field of the minimizing energy flow is

 vα=μmαβ∇i(∂γ∂∇iψβ+γgiqmσβ∇qφσ)

where, replacing with the expression for in equation (41), one has

 vα =μmαβ∇i(∂Gskgks∂∇iψβ+Gskgksgiqmσβ∇qφσ) =μmαβ∇i(2Gsimβζ∇sφζ+Gskgksgiqmσβ∇qφσ) =μmαβ∇θ(2Gθθmβζ∇θφζ+Gθθgθθgθθmσβ∇θφσ) =3μGθθ∇θ∇θφα

using

 ∇θ∇θφα=∂2φα∂θ2−Γθθθ∂φα∂θ

one arrives at

 (vxvy)=−3μGθθ(acosθbsinθ)−(a2−b2)cosθsinθa2sin2θ+b2cos2θ(−asinθbcosθ)

However, any tangential terms in the velocity, such as the second term in the above equation, have no influence on the flow of the interface such that the flow generated by the velocity field above is equivalent to the flow generated by

 (vxvy)=−3μGθθ(acosθbsinθ) (42)

Thus, turning into a flow , one has

 dφαdt(θ,t)=−3μGθθφα(θ,t)

for which there is only one solution

 φα(θ,t)=φα(θ,0)e−3μGθθt

 (φx(θ,t)φy(θ,t))=e−3μGθθt(acosθbsinθ)

Now given that the minimizing energy flow of the embedding is just the original embedding multiplied by a time dependent function, the flow is actually simply shrinking the ellipse in a homothetic manner to the center point of . Thus, assuming , the eccentricity is a constant of the flow

 e= ⎷1−(φy(π2,t)φx(0,t))2= ⎷1−(e−3μGθθtbe−3μGθθta)2=√1−(ba)2 (43)

and the scalar velocity of any point of the ellipse is

 v(θ,t)=√(vx)2+(vy)2=3μGθθe−3μGθθt√a2cos2θ+b2sin2θ (44)

with, in particular,

 v(0,t)=3μGθθe−3μGθθta (45) v(π2,t)=3μGθθe−3μGθθtb (46)

## 4 The numerical model and test case applications

In order to numerically simulate grain boundary configurations, a numerical model capable of representing boundary dynamics must be developed. As such, the following paragraphs are devoted to first reporting on the level set [Merriman1994, Osher1988] finite element method applied to this kind of boundary transport. Subsequently, the analytical shrinking ellipse test case is simulated and convergence of the numerical model is studied. Finally, a more general grain boundary energy is applied to a circular case in order to compare the classical isotropic formulation for the velocity field with the expression proposed in this work.

### 4.1 The finite element model

In order to solve the minimizing energy flow for the level set function using the Finite Element (FE) method, the problem must first expressed in a weak form, then it can be discretized in both time and space.

Consider the transport equation (32) and the definition in equation (33) where the relevant fields have already been extended from the smooth manifold to the enclosing manifold and is known. With any test function a weak form of the equation can be derived as

 ∫M∂ϕ∂tωdM−∫MμDαβ~∇α~∇βϕωdM =0 ∫M∂ϕ∂tωdM+∫M~∇α(μDαβω)~∇βϕdM−∫∂M~∇α(μDαβω~∇βϕ)d∂M =0

such that

 ∫M∂ϕ∂tωdM+∫MμDαβ~∇αω~∇βϕdM+∫Mμ~∇βDβαω~∇αϕdM =0 (47)

With respectively three distinct terms: the time derivative, a diffusive term and a convective contribution.

In this numerical framework, the Riemannian manifold is meshed using an unstructured simplicial grid generated using Gmsh [Geuzaine2009]. Thus, the smooth Riemannian manifold is approximated by a by parts manifold and any initially smooth field is approximated by a field whose component functions are in (i.e. a P1 field). As such, the level set field is approximated by a linear by parts (inside each cell) field . The details of the algorithm used to compute the distance function can be found in [Shakoor2015amm].

Thus, given a boundary energy map , with the boundary property space, the geometry dependence of can easily be evaluated at each node of the mesh . Considering that is only parameterized by the normal to the boundary for a given boundary, both values for and can be evaluated everywhere on the mesh. As such, the level set field induces a natural discretized extension of both and from to the entire discretized space . Outside the interface the

field has no physical meaning. However, this extension is necessary for solving the problem in a FE setting. The interpolated values of the fields at the interface

are also guaranteed to be the correct values given the linear by parts interpolation of . Figure 3 illustrates the construction for a circle and a particular choice of .

The is computed numerically on the mesh using a Superconvergent Patch Recovery method inspired from [Belhamadia2004] to obtain P1 fields. As such, both the diffusive tensor and the convective velocity are introduced explicitly into the formulation so as to create linearised approximations of the equation (47). Thus, solving the problem is completely linear without need for non-linear solvers or algorithms.

In this work a Generalized Minimal Residual (GMRES) type solver along with an Incomplete LU (ILU) type preconditionner, both linked from the PetsC open source libraries, are used unless specified otherwise. The system is assembled using typical P1 FE elements with a Streamline Upwind Petrov-Galerkin (SUPG) stabilization for the convective term

[Brooks1982]. The boundary conditions used are classical von Neumann conditions which guarantees the orthogonality of the level sets to the boundary of the domain. The discretization of time is obtained using a fully implicit backward Euler method with time step .

Because the resolution of the transport equation does not conserve the distance property of the level set field, the solution is reinitialized using the algorithm developed in [Shakoor2015]. Also, since the geometry of the interface evolves after each time increment, all the other fields must also be recomputed from the reinitialized level set at each step of the simulation. The complete procedure for the minimizing interface energy flow simulation is reported in Algorithm 1.

### 4.2 The shrinking ellipse

One now has an embedding and a way to represent it as a level set field on an unstructured mesh. One also has the FE formulation needed to simulate the dynamics of the minimizing energy flow of the interface. However, the boundary energy is not readily computable on the finite element mesh since it does not explicitly depend on the normal to the interface. Using equation (39) such that

 gθθ =(b2a2b2sin2θ+b2cos2θ) =b2(a2b2sin2θ+cos2θ)

which, if one considers , which should remain constant throughout the simulation if looking at the large and small axes of the ellipse at each instant, then can easily be extended throughout the mesh using

 nynx =rtanθ θ =arctan(1rnynx)

with

 γ(θ,t)=b(t)2(r(t)2sin2θ+cos2θ) (48)

Given the definition of the level set field, takes maximal values at the points within the ellipse furthest away from the interface, i.e. the center of the ellipse. Seeing as is the smallest of both ellipse axes and the level set is minimal distance valued, the value of the level set at the center of the ellipsis should be the value of the small axis. Therefore

 b(t)=maxq∈Mϕ(q,t) (49)

Also, implicit in the calculations in the previous section is the fact that

 ∂2γ∂~∇αϕ<