A Three-Dimensional Continuum Simulation Method for Grain Boundary Motion Incorporating Dislocation Structure
AA Three-Dimensional Continuum Model for Grain BoundaryMotion Incorporating Dislocation Structure
Xiaoxue Qin, Luchan Zhang ∗ , Yang Xiang ∗ Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay,Kowloon, Hong Kong
Abstract
We develop a continuum model for the dynamics of grain boundaries in three dimensionsthat incorporates the motion and reaction of the constituent dislocations. The continuummodel includes evolution equations for both the motion of the grain boundary and theevolution of dislocation structure on the grain boundary. The critical but computationallyexpensive long-range elastic interaction of dislocations is replaced by a projection formulationthat maintains the constraint of the Frank’s formula describing the equilibrium of the stronglong-range interaction. This continuum model is able to describe the grain boundary motionand grain rotation due to both coupling and sliding effects, to which the classical motion bymean curvature model does not apply. Comparisons with atomistic simulation results showthat our continuum model is able to give excellent predictions of evolutions of low anglegrain boundaries and their dislocation structures.
Keywords:
Grain boundary dynamics, Coupling and sliding motions, Dislocationdynamics, Frank’s formula, Projection method
1. Introduction
Grain boundaries are indispensable components in polycrystalline materials. The energyand dynamics of grain boundaries play essential roles in the mechanical and plastic behaviorsof the materials [1]. Most of the available continuum models for the dynamics of grain ∗ Corresponding author
Email addresses: [email protected] (Luchan Zhang), [email protected] (Yang Xiang)
Preprint submitted to Elsevier February 2, 2021 a r X i v : . [ c ond - m a t . m t r l - s c i ] J a n oundaries are mainly based on the motion driven by the capillary force which is proportionalto the local mean curvature of the grain boundary [1, 2, 3]. This motion is a process to reducethe interfacial energy (cid:82) S γdS , where S is the grain boundary and γ is the grain boundaryenergy density. If the energy density γ is fixed, the driving force based on variation ofthe total energy is in the normal direction of the grain boundary and is proportional toits mean curvature. There are many atomistic simulations and continuum models in theliterature for the grain boundary motion driven by the capillary force (by mean curvature),e.g., [4, 5, 6, 7, 8, 9, 10, 11].The decreasing of grain boundary energy density γ ( θ ) can also reduce the total energy.For a low angle grain boundary, this implies the decreasing the misorientation angle θ . Inthis case, the two grains on different sides of the grain boundary will rotate and cause arelatively rigid-body translation of the two grains along the boundary. This process is calledsliding motion of grain boundaries [12, 13, 14, 15, 16, 17, 18].There is a different type of grain boundary motion which is called coupling motion[19, 20, 21], in which the normal motion of the grain boundary induces a tangential motionproportionally. In the coupling motion, the energy density γ ( θ ) can increase although thetotal energy (cid:82) S γdS is decreasing. Cahn and Taylor [21] proposed a unified theory for thecoupling and sliding motions of the grain boundary and demonstrated the theory based ondislocation mechanisms for a circular low angle grain boundary in two dimensions. Espe-cially, the coupling motion of the grain boundary is associated with dislocation conservationduring the motion of the grain boundary. The Cahn-Taylor theory and mechanisms ofmotion and reaction of the constituent dislocations have been examined by atomistic simu-lations and experiments [20, 22, 23, 24, 25, 26, 27, 28, 29]. It has been shown in Ref. [30] bya dislocation model and experimental observations that conservation and annihilation of theconstituent dislocations may lead to cancelation of the coupling and sliding motions of thegrain boundary, leading to the classical motion by curvature. A continuum model has beendeveloped based on the motion and reaction of the constituent dislocations for the dynamicsof low angle grain boundaries in two dimensions [31]. Their model can describe both thecoupling and sliding motions of low angle grain boundaries. Recently, they have derived a2ore efficient formulation [32]. A continuum model that generalizes the Cahn-Taylor theorybased on mass transfer by diffusion confined on the grain boundary has been proposed [33],and numerical simulations based on this generalization were performed using the level setmethod [34]. Crystal plasticity models that include shear-coupled grain boundary motion inthe phase field framework of Ref. [15] have been developed [35, 36], in which the geometricnecessary dislocation (GND) tensor/lattice curvature tensor was used to represent the actualdislocation distributions on the grain boundaries. All these continuum models are for grainboundaries in two dimensions.There are only limited studies in the literature for the three-dimensional coupling andsliding motions of grain boundaries. Grain boundary motion and grain rotation in bccand fcc bicrystals composed of a spherical grain embedded in a single crystal matrix werestudied by using three-dimensional phase field crystal models [28, 29], and properties ofgrain boundaries and their dislocation structures as the grain boundary evolves have beenexamined. Although these atomistic-level phase field crystal simulations are able to providedetailed information associated with the coupling and sliding motions of grain boundaries inthree dimensions, three-dimensional continuum models of the dynamics of grain boundariesincorporating their dislocation structures are still desired for larger scale simulations.In this paper, we generalize the two-dimensional continuum model for grain boundarydynamics in Ref. [31, 32] to three dimensions, where grain boundaries and their constituentdislocations are curved in general. The three-dimensional continuum model for the dynamicsof grain boundaries incorporates the motion and reaction of the constituent dislocations, andis able to describe the grain boundary motion and grain rotation due to both coupling andsliding effects, to which the classical motion by mean curvature model does not apply. Thecontinuum model includes evolution equations for both the motion of the grain boundary andthe evolution of dislocation structure on the grain boundary. The evolution of orientation-dependent continuous distributions of dislocation lines on the grain boundary is based onthe simple representation using dislocation density potential functions [37].The continuum model contains a long-range force in the form of singular integrals, whoseevaluation is time-consuming. With the fact that the long-range force is so strong that an3quilibrium state described by the Frank’s formula [38, 39, 37] is quickly reached duringthe evolution, we replace the long-range force by a constraint of the Frank’s formula. Wefurther solve the constraint evolution problem using the projection method, leading to anew, efficient continuum formulation. Using the obtained continuum model, we performnumerical simulations for the evolution of initially spherical low angle grain boundaries infcc and bcc crystals. Simulation results are compared with atomistic simulation resultsusing phase field crystal models [28, 29]. In particular, we explain the anisotropic motionof the grain boundary in the directions along and normal the rotation axis in terms of theconstraint of the Frank’s formula.This paper is organized as follows. In Sec. 2, we review the description of a grainboundary in three dimensions and the continuum framework for grain boundary dislocationstructure in Ref. [37]. We also briefly review the Frank’s formula for the equilibrium dislo-cation structure of a grain boundary. In Sec. 3, we present our continuum model for grainboundary dynamics in three dimensions incorporating the evolution of the underlying dis-location structure, and develop a new, efficient formulation based on constrained evolution.Simulations for the evolution of initially spherical low angle grain boundaries in fcc andbcc crystals using our continuum model and comparisons with results of phase field crystalsimulations [28, 29] are presented in Sec. 4.
2. Geometric descriptions and continuum framework for grain boundaries andtheir dislocation structures
Consider a grain boundary S in three dimensions parameterized by ( u, v ). That is, apoint on grain boundary S is r = r ( u, v ) = ( x ( u, v ) , y ( u, v ) , z ( u, v )) . (1)With this parametrization, we have the following formulas for two unit tangent vectors T ( ) , T ( ) and the unit normal vector n : T ( ) = r u (cid:107) r u (cid:107) , T ( ) = r v (cid:107) r v (cid:107) , n = T ( ) × T ( ) (cid:107) T ( ) × T ( ) (cid:107) = r u × r v (cid:107) r u × r v (cid:107) , (2)4here r u = (cid:0) ∂x∂u , ∂y∂u , ∂z∂u (cid:1) and r v = (cid:0) ∂x∂v , ∂y∂v , ∂z∂v (cid:1) . The area of the grain boundary is S A = (cid:90) S (cid:107) r u × r v (cid:107) d u d v. (3)Given a function f on the grain boundary S , its surface gradient is ∇ S f = ( ∇ − n ( n · ∇ )) f = (cid:18) (cid:107) r v (cid:107) (cid:107) r u × r v (cid:107) f u − r u · r v (cid:107) r u × r v (cid:107) f v (cid:19) r u + (cid:18) (cid:107) r u (cid:107) (cid:107) r u × r v (cid:107) f v − r u · r v (cid:107) r u × r v (cid:107) f u (cid:19) r v . (4) In this subsection, we review the continuum framework for the dislocation structure oflow angle grain boundaries proposed in Ref. [37]. In this continuum framework, dislocationdensity potential functions (scalar functions) defined on the grain boundary are employed todescribe the dislocation structure on the grain boundary. This provides a simple and accuratedescription of the orientation dependent densities of dislocations which are connected lines.Energy of the grain boundary and driving forces for the dynamics are given based on thedislocation structure.
Figure 1: A dislocation density potential function η defined on a grain boundary S . Its integer-value contourlines represents the array of dislocations with the same Burgers vector. For a dislocation density potential function η defined on a grain boundary S , the con-stituent dislocations of Burgers vector b are given by the contour lines of η : η = j , for5nteger j. See Fig. 1 for an example of dislocation structure on a spherical grain boundaryand η defined on it. From the dislocation density potential function η , the inter-dislocationdistance D can be calculated by D = 1 (cid:107)∇ S η (cid:107) , (5)and the dislocation direction is given by t = ∇ S η × n (cid:107)∇ S η (cid:107) , (6)where ∇ S η is the surface gradient of η on S : ∇ S η = ( ∇ − n ( n · ∇ )) η. (7)Multiple dislocation density potential functions are used for dislocations with different Burg-ers vectors.Suppose there are J arrays of dislocations with Burgers vectors b ( j ) , j = 1 , , · · · , J ,respectively, on the grain boundary S , and they are described by the dislocation densitypotential functions η j , j = 1 , , · · · , J , respectively. The continuum formulation of theelastic energy is E tot = E long + E local , (8)where E long is the long-range interaction energy of dislocations and E local is the local dislo-cation line energy: E long = 12 J (cid:88) i =1 J (cid:88) j =1 (cid:90) S d S i (cid:90) S d S j [ µ π ( ∇ S η i × n i · b ( i ) )( ∇ S η j × n j · b ( j ) ) r ij − µ π ( ∇ S η i × n i ) × ( ∇ S η j × n j ) · ( b ( i ) × b ( j ) ) r ij + µ π (1 − ν ) ( ∇ S η i × n i · b ( i ) ) · ( ∇ ⊗ ∇ r ij ) · ( ∇ S η j × n j · b ( j ) ) (cid:21) , (9) E local = J (cid:88) j =1 (cid:90) S µ ( b ( j ) ) π (1 − ν ) (cid:34) − ν ( ∇ S η j × n · b ( j ) ) ( b ( j ) ) (cid:107)∇ S η j (cid:107) (cid:35) (cid:107)∇ S η j (cid:107) log 1 r g (cid:107)∇ S η j (cid:107) d S j . (10)Here r ij = (cid:107) X i − X j (cid:107) , where X i and X j are the points varying on the grain boundary S andare associated with the surface integral dS i and dS j , respectively. n j is the normal direction6f the surface S associated with the surface integral dS j , notation ⊗ is the tensor productoperator, the gradient in the term ∇ ⊗ ∇ r ij is taken with respect to X i , and b ( j ) = (cid:107) b ( j ) (cid:107) .The elastic constants µ is the shear modulus and ν is the Poisson’s ratio.The driving forces for the dynamics of the grain boundary and the dislocation structureare associated with the variations of the total energy in Eqs. (8)–(10): δE tot δr = − J (cid:88) j =1 (cid:107)∇ S η j (cid:107) ( f ( j )long + f ( j )local ) · n , (11) δE tot δη j = ( f ( j )long + f ( j )local ) · ∇ S η j (cid:107)∇ S η j (cid:107) , (12)where f ( j )long is the continuum long-range force, f ( j )local the local force on the b ( j ) -dislocations(i.e., the array of dislocations with Burgers vector b ( j ) ), and r is the signed distance in thenormal direction of the grain boundary. In fact, the total force f ( j )long + f ( j )local is the continuumlimit [37] from the Peach-Koehler force on dislocations in discrete dislocation dynamics [40].The long-range force f ( j )long and the local force f ( j )local have the following formulations: f ( j )long =( σσσ tot · b ( j ) ) × (cid:18) ∇ S η j (cid:107)∇ S η j (cid:107) × n (cid:19) , (13) f ( j )local = f ( j ) lt + f ( j ) p , (14) f ( j ) lt = µ π (1 − ν ) κ d (cid:104) (1 + ν )( b ( j ) t ) + (1 − ν )( b ( j ) N ) + ( b ( j ) B ) (cid:105) log 1 r g (cid:107)∇ S η j (cid:107) n ( j ) d − µν π (1 − ν ) κ d b ( j ) N b ( j ) B t ( j ) × n ( j ) d , (15) f ( j ) p = µ π (1 − ν ) κ ( j ) p (cid:104) ( b ( j ) ) − ν ( b ( j ) t ) (cid:105) n ( j ) p . (16)Here σσσ tot is the total stress field, which includes the long-range stress field generated by thedislocation arrays on S and other stress fields. The long-range stress field generated by thedislocation arrays is σσσ ( X ) = J (cid:88) j =1 µ π (cid:90) S (cid:20)(cid:18) ∇ r × b ( j ) (cid:19) ⊗ ( ∇ S η j × n ) + ( ∇ S η j × n ) × (cid:18) ∇ r × b ( j ) (cid:19) + 11 − ν (cid:0) b ( j ) × ( ∇ S η j × n ) · ∇ (cid:1) ( ∇ ⊗ ∇ − I ∆) r (cid:21) d S, (17)7here r = (cid:107) X − X S (cid:107) with points X S varying on the grain boundary S and ∇ S η j and n being evaluated at X S . The local force f ( j )local consists of the dislocation line tension force f ( j ) lt and the curvature force f ( j ) p associated the curve on S that is normal to the local dislocation,where κ d is the curvature of dislocation line, n d is the normal direction of dislocation, κ p and n p are the curvature and normal direction of the curve on S that is normal to the locationdislocation, respectively, and κ d n d = ( ∇ S t ) · t , κ p n p = (cid:16) ∇ S ∇ S η j (cid:107)∇ S η j (cid:107) (cid:17) · ∇ S η j (cid:107)∇ S η j (cid:107) , b t = b · t , b N = b · n d , b B = b · ( t × n d ). The equilibrium dislocation structure on a grain boundary is governed by the Frank’sformula [38, 39, 37]: θ ( V × a ) − J (cid:88) j =1 b ( j ) ( ∇ S η j · V ) = , (18)where θ is the misorientation angle, a = ( a , a , a ) is the rotation axis, V is any vectorin the grain boundary’s tangent plane, and there are J arrays of dislocations with Burgersvectors b ( j ) , j = 1 , , · · · , J , respectively, on the grain boundary.The Frank’s formula holds if and only if the long-range elastic fields generated by thegrain boundary cancel out [38, 39, 37]. The classical Frank’s formula is based on planar grainboundaries [38, 39]. It has been shown in Ref. [37] that the equivalence of the Frank’s formulato the cancelation of the long-range elastic fields also holds for curved grain boundaries.Eq. (18) is the Frank’s formula for a curved grain boundary, where the dislocation structureis represented by the dislocation density potential functions { η j } . Eq. (18) holds for anyvector V in the grain boundary’s tangent plane if and only if it holds for the two tangentvectors V = r u and V = r v , where u and v are the two parameters of the grain boundary.
3. Continuum model for grain boundary dynamics in three dimensions
In this section, we present our continuum model for the dynamics of grain boundaries inthree dimensions that incorporates the coupling and sliding motions. The continuum modelconsists of an equation for the motion of the grain boundary and an equation for the evolution8f its dislocation structure. The continuum model is based on the simple representation ofdislocation structure on grain boundaries using dislocation density potential functions andthe continuum energy formulation in Ref. [37]. The continuum model is able to describeboth the coupling and sliding motions of the grain boundary. The misorientation angleassociated with the grain boundary is calculated by the grain boundary profile and itsdislocation structure.This continuum model has a long-range force in the form of a singular integral over theentire grain boundary. Simulations directly using this continuum formulation will be time-consuming. Moreover, when the source term of dislocations (for annihilation and generationof dislocations) is directly expressed in terms of variation of the grain boundary energy, theresulting evolution equation of the dislocation structure is not well-posed due to the factthat the (local) energy of the grain boundary is a concave function of dislocation densities.Modifications of the continuum formulation will be developed to address these limitations.
We first present the three-dimensional continuum model with long range force for thedynamics of a grain boundary and its dislocation structure.
Continuum model with long range force v n = M d J (cid:88) j =1 (cid:107)∇ S η j (cid:107) (cid:80) Jk =1 (cid:107)∇ S η k (cid:107) ( f ( j )long + f ( j )local ) · n , (19) ∂η j ∂t = − M d f ( j )long · ∇ S η j + s j . (20)Eq. (19) gives the normal velocity of the grain boundary. In this equation, f ( j )long (given inEqs. (13) and (17)) and f ( j )local (given in Eqs. (14)–(16)) are the long-range force and the localforce, respectively. The velocity of a point on a grain boundary is the weighted average ofthe velocities of dislocations on the grain boundaries with different Burgers vectors, and M d is the mobility of dislocations.Eq. (20) describes the evolution of the dislocations structure on the grain boundary.The first term on the right-hand side describes the motion of dislocations along the grain9oundary, which is in the form of conservation law driven by the continuum long-rangeelastic force. The continuum long-range force f long is the leading order continuum force,and it maintains a stable dislocation structure [41]. Recall that the cancelation of the long-range force is equivalent to the Frank’s formula for an equilibrium dislocation structure.The second term s j on the right hand side of Eq. (20) is the source term that accounts fordislocation reaction [20, 25, 28], which on the continuum level may be modeled based onthe variation of the local energy of the grain boundary and whose exact formulation will bediscussed in Sec. 3.3.The misorientation angle θ during the evolution can be calculated using Frank’s formulain Eq. (18) based on the grain boundary profile and the dislocation structure, and theformula is θ = 1 S A (cid:90) S J (cid:88) j =1 ( η ju + η jv )( r u + r v ) × a · b ( j ) (cid:107) ( r u + r v ) × a (cid:107) d S. (21)where η ju = ∂η j ∂u , η jv = ∂η j ∂v , and recall that S A is the area of the grain boundary S thatcan be calculated using Eq. (3). Note that although the dislocation structure on the grainboundary is driven into equilibrium by the long-range force generated by the dislocations,at any moment during the evolution, the dislocation structure on the grain boundary is notnecessarily in equilibrium exactly, i.e., there may not be a constant value θ that makes theFrank’s formula in Eq. (18) holds at every point on the grain boundary. The formula of θ inEq. (21) is an estimation for the misorientation angle based on some average over the grainboundary. Its derivation is as follows.Substituting V = r u and V = r v into Frank’s formula in Eq. (18), we have θ ( r u × a ) − J (cid:88) j =1 b ( j ) η ju =0 , (22) θ ( r v × a ) − J (cid:88) j =1 b ( j ) η jv =0 . (23)Here we have used ∇ S η j · r u = η ju and ∇ S η j · r v = η jv . Adding the two equations (22) and1023), multiplying both size of the summation by ( r u + r v ) × a , we have θ (cid:107) ( r u + r v ) × a (cid:107) = J (cid:88) j =1 ( η ju + η jv )( r u + r v ) × a · b ( j ) . (24)Integrating over the entire grain boundary S , we obtain the formula of θ in Eq. (21), whichis exact if the Franks’ formula holds and serves as an estimate formula otherwise. The continuum formulation given by Eqs. (19) and (20) contains the long-range elasticforce f ( j )long (given in Eqs. (13) and (17)), which is a singular integral over the entire grainboundary surface. Numerically, computation of such long-range force with reasonable accu-racy is complicated and time-consuming even in two-dimensional cases [31, 32]. Actually,the long-range interaction between the grain boundary dislocations is so strong that anequilibrium state described by the Frank’s formula [38, 39, 37] is quickly reached during theevolution of the grain boundaries. We simply assume that the Frank’s formula always holdsduring the evolution of the grain boundary, which will not lead to significant change in grainboundary motion. Under on this assumption, the long-range elastic force is replaced by theconstraint of Frank’s formula in the continuum dynamics model. This assumption has beenvalidated in the two-dimensional model [31, 32], and here we apply this assumption to thecontinuum model in three dimensions. The new formulation is: Continuum model with constraint v n = M d J (cid:88) j =1 (cid:107)∇ S η j (cid:107) (cid:80) Jk =1 (cid:107)∇ S η k (cid:107) f ( j )local · n , (25) ∂η j ∂t = s j , (26)subject to h = θ ( V × a ) − J (cid:88) j =1 b ( j ) ( ∇ S η j · V ) = . (27)Here, the constraint in Eq. (27) is the Frank’s formula.The constrained evolution model in Eqs. (25)-(27) can be implemented using the pro-jection method. We evolve the grain boundary and dislocation without the constraint first,11nd then project the grain boundary and its dislocation structure to a new configurationthat satisfies the constraint (i.e. the Frank’s formula). The projection procedure can besolved, leading to the following formulation. Here without loss of generality, suppose thatthe rotation axis is in the + z direction, i.e., a = (0 , , Continuum model using projection method
From t n to t n +1 = t n + δt , v ∗ = (cid:32) M d J (cid:88) j =1 (cid:107)∇ S η j (cid:107) (cid:80) Jk =1 (cid:107)∇ S η k (cid:107) f ( j )local · n (cid:33) n , (28) r ∗ = r n + v ∗ δt, (29) η n +1 j = η nj + s nj δt, (30) δθ = θ ( r ∗ u , r ∗ v , η n +1 ju , η n +1 jv ) − θ ( r nu , r nv , η nju , η njv ) , (31) v = (cid:18) − δθθδt ( x − c ) , − δθθδt ( y − c ) , v ∗ (cid:19) + (cid:32) − θ J (cid:88) j =1 b ( j )2 s j , θ J (cid:88) j =1 b ( j )1 s j , (cid:33) , (32) r n +1 = r n + v δt. (33)where constanta c and c can be determined numerically by the condition that the projectionprocedure alone does not lead to extra rigid translation of the grain boundary. When thetop point of the grain boundary in the + z direction always has a velocity in the z directionduring the evolution due to some symmetry, we simply have c = c = 0 if the z axis is setpassing through that point. The values of misoreintation angle θ in Eq. (31) are calculatedby using Eq. (21).Eq. (29) is the virtual evolution of the grain boundary after the time step δt , using thevelocity v ∗ = ( v ∗ , v ∗ , v ∗ ) due to the local force without the constraint as given in Eq. (28).Eq. (30) is the evolution of dislocation structure due to dislocation reaction without theconstraint. The change of misorientation angle δθ during this time step is determined bythis virtual evolution as given in Eq. (31). Based on this obtained δθ , the grain boundaryvelocity is adjusted by projection of the new configuration of the grain boundary to a statethat satisfies the constraint of Frank’s formula, and the obtained velocity v is given inEq. (32), following which the grain boundary is actually evolved.12ote that in the projection procedure, we essentially adjust the local value of θ deter-mined by the Frank’s formula in Eq. (27) to achieve a uniform misorientation angle θ overthe entire grain boundary. This procedure should not lead to additional rigid translationof the grain boundary. The two constants c and c in the projected velocity formula inEq. (32) can be determined by this condition. For some symmetric configuration of thegrain boundary, we have c = c = 0 as described above. Also note that the Frank’s formulain Eq. (27) does not impose any constraint in the direction parallel to the rotation axis a ,which is in + z direction. As a result, we simply keep the z -component v ∗ of the virtualvelocity in the projected velocity formula in Eq. (32).In the projected velocity formula in Eq. (32), the first term describes the pure couplingmotion of the grain boundary, the second term describes the additional effect of the slidingmotion of the grain boundary due to dislocation reaction. Derivation of this formulation isgiven below. Derivation of the solution of projection method
Suppose that the grain boundary velocity is v , and the Frank’s formula in Eqs. (22) and(23) hold at the time t n . After a small time step δt , if the Frank’s formula still holds, wehave ( δθ r u + θδt v u ) × a − J (cid:88) j =1 b ( j ) δη ju =0 , (34)( δθ r v + θδt v v ) × a − J (cid:88) j =1 b ( j ) δη jv =0 . (35)Here we have used δ r u = δt v u and δ r v = δt v v , where v u = ∂ v ∂u and v v = ∂ v ∂v .Integrating Eq. (34) with respect to u , we have( δθ r ( u, v ) + θδt v ( u, v )) × a − J (cid:88) j =1 b ( j ) δη j ( u, v )= ( δθ r (0 , v ) + θδt v (0 , v ))) × a − J (cid:88) j =1 b ( j ) δη j (0 , v ) . (36)13imilarly, integrating Eq. (35) with respect to v , we have( δθ r ( u, v ) + θδt v ( u, v )) × a − J (cid:88) j =1 b ( j ) δη j ( u, v )= ( δθ r ( u,
0) + θδt v ( u, × a − J (cid:88) j =1 b ( j ) δη j ( u, . (37)Notice that the left-hand sides of Eqs. (36) and (37) are equal, whereas the right-handside of Eq. (36) depends only on v and the right-hand side of Eq. (37) depends only on u .Thus the right-hand sides of Eqs. (36) and (37) must equal to the same constant independentof u and v , denoted by C . That is,( δθ r ( u, v ) + θδt v ( u, v )) × a − J (cid:88) j =1 b ( j ) δη j ( u, v ) = C . (38)As discussed above, the constant C can be determined by the condition that the projectionprocedure alone does not lead to extra rigid translation of the grain boundary, and can becalculated numerically during the evolution. Note that the dislocation structure dependsonly on {∇ S η j } , and does not change by adding a constant.In the case when the top point of the grain boundary in the + z direction always has avelocity in the z direction (e.g., due to some symmetry), we set the z axis passing throughthat point, i.e., that point is r = (0 , , z ) during the evolution. In this case, at that point,we have ( δθ r + θδt v ) × a = . Thus, we have C = , and v = − δθθδt x − θ J (cid:88) j =1 b ( j )2 δη j δt , (39) v = − δθθδt y + 1 θ J (cid:88) j =1 b ( j )1 δη j δt , (40)where v = ( v , v , v ).The Frank’s formula in Eqs. (22) and (23) does not impose any condition on the velocityin the direction of the rotation axis, i.e., the z direction. Thus we keep the projection ofthe velocity due to the local force in the direction of the rotation axis, which is v ∗ where v ∗ = ( v ∗ , v ∗ , v ∗ ) is the velocity due to the local force given in Eq. (28).14n summary, the grain boundary velocity after the projection based on the constraint ofFrank’s formula is v = ( v , v , v ∗ ), where v and v are given in Eqs. (39) and (40), and v ∗ is the z -component of the velocity v ∗ due to the local force given in Eq. (28). This is theprojected velocity formula in Eq. (32). A commonly adopted formulation for the dislocation reaction term s j in the evolutionof dislocation structure in Eq. (20) or Eq. (26) is based on variation of the grain boundaryenergy, i.e., E local in Eq. (10). (Recall that equilibrium with respect to the long-range energyis almost reached during the evolution of the grain boundary.) However, the (local) grainboundary energy density γ gb is not convex as a function of {∇ S η j } , where γ gb = J (cid:88) j =1 µ ( b ( j ) ) π (1 − ν ) (cid:32) − ν ( ∇ S η j × n · b ( j ) ) ( b ( j ) ) (cid:107)∇ S η j (cid:107) (cid:33) (cid:107)∇ S η j (cid:107) log 1 r g (cid:107)∇ S η j (cid:107) d S j . (41)This nonconvexity will lead to an ill-posed formulation when using gradient flow of the grainboundary energy for the evolution of dislocation structure in Eq. (20) or Eq. (26). This canbe understood as follows. Neglecting the orientation dependence factor, the contributionof b ( j ) -dislocations in the energy density γ gb is essentially −(cid:107)∇ S η j (cid:107) log (cid:107)∇ S η j (cid:107) , which is aconcave function of (cid:107)∇ S η j (cid:107) . Recall that from the definition, we always have (cid:107)∇ S η j (cid:107) ≤ η j .In order to obtain a well-posed gradient flow formulation, we use the components of ∇ S η j , i.e., η ju and η jv as independent variables instead of η j itself in the evolution equationof dislocation structure. Variations of the grain boundary energy are taken with respect to η ju and η jv instead of η j . This leads to a new problem that η ju and η jv are not independent,and they are related by ∂η ju ∂v − ∂η jv ∂u = 0. The physical meaning of these relations is thatthe dislocations are connected lines on the grain boundary. We including these relationsas further constraints in our continuum model. Using these treatments, the continuumformulation with constraint in Eqs. (25)–(27) can be rewritten as: Well-posed continuum model with constraints n = M d J (cid:88) j =1 (cid:107)∇ S η j (cid:107) (cid:80) Jk =1 (cid:107)∇ S η k (cid:107) f ( j )local · n , (42) ∂η ju ∂t = − M r ∂γ gb ∂η ju , ∂η jv ∂t = − M r ∂γ gb ∂η jv , (43)subject to h = θ ( V × a ) − J (cid:88) j =1 b ( j ) ( ∇ S η j · V ) = , (44) ∂η ju ∂v − ∂η jv ∂u = 0 . (45)Here M r is the mobility associated with dislocation reaction.Numerically, we implement the additional constraint in Eq. (45) using a projectionmethod similar to that for fluid dynamics problems [42]. Recall that the evolution of η ju and η jv in Eq. (43) is to minimize the energy E local = (cid:82) S γ gb dS by gradient flow. In order toimplement the constraint in Eq. (45), we introduce a Lagrangian function: L = (cid:90) S (cid:32) γ gb + J (cid:88) j =1 λ j (cid:18) ∂η ju ∂v − ∂η jv ∂u (cid:19)(cid:33) d S, (46)where λ j , j = 1 , , · · · , J , are Lagrange multipliers associated with the constraints. Theevolution of dislocation structure in Eq. (43) now becomes ∂η ju ∂t = − M r ∂γ gb ∂η ju + ∂λ j ∂v , (47) ∂η jv ∂t = − M r ∂γ gb ∂η jv − ∂λ j ∂u . (48)Here the coefficients of ∂λ j ∂v and ∂λ j ∂u in these equations are set to be 1.During the evolution in the time step from t n to t n +1 = t n + δt , we separate the contri-butions from γ gb and λ j into two steps: η ∗ ju = η nju − M r ∂γ gb ∂η ju (cid:12)(cid:12)(cid:12)(cid:12) t n · δt, η ∗ jv = η njv − M r ∂γ gb ∂η jv (cid:12)(cid:12)(cid:12)(cid:12) t n · δt, (49) η n +1 ju = η ∗ ju + ∂λ n +1 j ∂v δt, η n +1 jv = η ∗ jv − ∂λ n +1 j ∂u δt. (50)In order to satisfy the constraint ∂η n +1 ju ∂v − ∂η n +1 jv ∂u = 0, using Eq. (50), we have the formula forupdating λ j : (cid:52) λ n +1 j = 1 dt (cid:18) ∂η ∗ jv ∂u − ∂η ∗ ju ∂v (cid:19) , (51)16here (cid:52) is the Laplace operator.Combining this numerical algorithm with that for the implementation of the continuummodel with the constraint of Frank’s formula given in Eqs. (28)–(33), we have the numericalalgorithm for the well-posed continuum model with constraints in Eqs. (42)–(45): Numerical Algorithm
From t n to t n +1 = t n + δt , v ∗ = (cid:32) M d J (cid:88) j =1 (cid:107)∇ S η j (cid:107) (cid:80) Jk =1 (cid:107)∇ S η k (cid:107) f ( j )local · n (cid:33) n , (52) r ∗ = r n + v ∗ δt, (53) η ∗ ju = η nju − M r ∂γ gb ∂η ju (cid:12)(cid:12)(cid:12)(cid:12) t n δt, η ∗ jv = η njv − M r ∂γ gb ∂η jv (cid:12)(cid:12)(cid:12)(cid:12) t n δt, (54) (cid:52) λ n +1 j = 1 δt (cid:18) ∂η ∗ jv ∂u − ∂η ∗ ju ∂v (cid:19) , (55) η n +1 ju = η ∗ ju + ∂λ n +1 j ∂v δt, η n +1 jv = η ∗ jv − ∂λ n +1 j ∂u δt, (56) δθ = θ ( r ∗ u , r ∗ v , η n +1 ju , η n +1 jv ) − θ ( r nu , r nv , η nju , η njv ) , (57) v = (cid:18) − δθθδt ( x − c ) , − δθθδt ( y − c ) , v ∗ (cid:19) + (cid:32) − θ J (cid:88) j =1 b ( j )2 δη j δt , θ J (cid:88) j =1 b ( j )1 δη j δt , (cid:33) , (58) r n +1 = r n + v δ t . (59)Recall that in this numerical formulation, f ( j )local is the local force on the grain boundarygiven in Eqs. (14)–(16), γ gb is the local grain boundary energy given in Eq. (41), the misori-entation angle θ is calculated using Eq. (21), and the constants c and c in the projectedvelocity v is determined numerically by the condition that the projection procedure alonedoes not lead to extra rigid translation of the grain boundary as discussed before (afterEq. (33)). The Poisson equation in (55) can be solved using a finite difference method. Initial dislocation structure
We assume that the initial grain boundary has an equilibrium dislocation structure thatsatisfies the Frank’s formula and has the lowest energy. See Ref. [43] for the method based17n constrained energy minimization to find the equilibrium dislocation structure on a curvedlow angle grain boundary, which is a generalization of the model for planar low angle grainboundaries [44] examined extensively by comparisons with atomistic simulation results.
The classical grain boundary dynamics models are based upon the motion by mean curva-ture (or curvature flow) to reduce the total grain boundary energy, with fixed misorientationangle and energy density of the grain boundary [2, 3, 1, 45]. Using the notation E local forthe total (local) grain boundary energy, a general, anisotropic curvature flow is v n = − M δE local δr , (60)where v n is the velocity in the normal direction of the grain boundary, r is the signed distancein the normal direction, and M > M = M d (cid:80) Jk =1 (cid:107)∇ S η k (cid:107) , using Eq. (11)), and in addition to that, our modelalso has a constraint due to Frank’s formula and evolution of the dislocation structure onthe grain boundary. The misorientation angle θ is fixed during the evolution of the grainboundary by curvature flow; whereas in our model, θ also evolves and is determined bythe grain boundary profile and dislocation structure at any moment (by Eq. (21)), andespecially, θ and the grain boundary energy density γ gb may increase during the evolution.(See numerical examples in the next section.)Our continuum model can be reduced to the classical curvature flow in a special casewhere evolutions of the grain boundary and the dislocation structure are balanced to keep aconstant misorientation angle θ . We use the numerical formulation using projection methodin Eqs. (28)–(33) for this comparison. Now misorientation angle θ is fixed, i.e., δθ = 0.Under this condition, the projected velocity formulation in Eq. (32) becomes v = (cid:32) − θ J (cid:88) j =1 b ( j )2 s j , θ J (cid:88) j =1 b ( j )1 s j , v ∗ (cid:33) . (61)18t reduces to the velocity of the curvature flow, i.e., v ∗ , if the dislocation reaction terms { s j } satisfy − θ (cid:80) Jj =1 b ( j )2 s j = v ∗ and θ (cid:80) Jj =1 b ( j )1 s j = v ∗ .
4. Numerical simulations
In this section, we perform numerical simulations of grain boundary dynamics using ournumerical algorithm in Eqs. (52)-(59), which is a numerical implementation of the continuummodel of constrained evolution in Eqs. (42)–(45).
We first consider grain boundaries in fcc Al. We choose the directions [¯110], [¯1¯12], [111]to be the x , y and z directions, respectively. In this coordinate system, the six Burgersvectors are b (1) = (1 , , b , b (2) = (cid:16) , √ , (cid:17) b , b (3) = (cid:16) , − √ , (cid:17) b , b (4) = (cid:16) , √ , − √ (cid:17) b , b (5) = (cid:16) , √ , √ (cid:17) b , and b (6) = (cid:16) − , √ , √ (cid:17) b , where b is the magnitude of the Burgersvectors. In Al, b = 0 . ν = 0 . a is inthe [111] direction, i.e., + z direction.We study the evolution of an initially spherical grain boundary, whose radius is R = 20 b and misorientation angle is θ = 5 ◦ . There are three sets of dislocations with Burgers vectors b (1) , b (2) , and b (3) , respectively, in the equilibrium dislocation structure on this initial,spherical grain boundary; see the top image in Fig. 2(a).In the dynamics simulation, the grain boundary is parameterized using spherical coor-dinates R = R ( α, β ), for 0 ≤ α < π and 0 ≤ β ≤ π . Here α is the angle between theposition vector of a point on the grain boundary and the x axis, and β is the angle betweenthe position vector of the point and the z axis. Initially, R ( α, β ) = 20 b . The ( α, β ) domainis discretized into 40 ×
20 uniform grids during the evolution. The center of the sphericalgrain boundary is the origin (0 , ,
0) in the coordinate system. Due to symmetry, the twoconstants c = c = 0 in the projected velocity formula in Eq. (58). We first consider the grain boundary motion without dislocation reaction, i.e. the reac-tion mobility M r = 0 in Eq. (54) (and Eq. (43)), and accordingly δη j = 0 in Eq. (58). This19s the pure coupling motion. (a) (b) (c) (d)Figure 2: Shrinkage of an initially spherical grain boundary in fcc under pure coupling motion, i.e., withoutdislocation reaction. The rotation axis is the z direction ([111]), and the initial misorientation angle θ = 5 ◦ .The upper panel of images show the three-dimensional view of the grain boundary during evolution. Themiddle panel of images show the grain boundary during evolution viewed from the + z direction ([111]),and the lower panel of images show the grain boundary during evolution viewed from the + x direction([¯110]). Dislocations with Burgers vectors b (1) , b (2) and b (3) are shown by blue, black and red lines,respectively. Length unit: b . (a) The initial spherical grain boundary. (b), (c), and (d) Configurations attime t = 10 /M d µ, /M d µ, /M d µ , respectively. Fig. 2 shows the shrinkage of the spherical grain boundary under this pure couplingmotion. The grain boundary eventually disappears. In this case, since M r = 0 and δη j = 0,the grain boundary velocity in Eq. (58) becomes v = − δθθδt ( x, y, , , v ∗ ). In the directionnormal the rotation axis, i.e., in the xy plane, the velocity component is in the inwardradial direction, as in the two-dimensional model [31, 32]; this is adjusted from the velocity20omponent due to curvature flow in order to satisfy the Frank’s formula. Whereas in thedirection of the rotation axis, i.e., the z direction, there is no constraint imposed by theFrank’s formula, and the velocity component is the same as that in the curvature flow.As an example, we consider the cross-section of the grain boundary with the z = 0plane (i.e., cross-section normal to the [111] rotation axis), which is the equator of thegrain boundary in the three dimensional view in the upper panel in Fig. 2 and is a circle(the outer circle) as shown in the second panel in Fig. 2 for the view from + z direction.The grain boundary along this circular cross-section is pure tilt, which is similar to thetwo-dimensional grain boundary discussed in Ref. [31, 32]. Along this circle, during theevolution, we have v ∗ = 0, and the grain boundary velocity is v = − δθθδt ( x, y, z = 0 plane. Thus the cross-sectionkeeps the circular shape as it shrinks during the evolution, as shown in the second panel inFig. 2. This shape-preserving evolution agrees with the results of the two-dimensional grainboundary dynamics models [33, 31, 32] and shrinkage of circular grain boundaries in twodimensions by molecular dynamics [20] and phase field crystal [26] simulations. However,here the changing rate of misorientation angle δθδt in the velocity formula is depending onthe entire grain boundary in three dimensions by Eq. (57), and is not just depending on thecircular cross-section itself as in the two dimensional continuum model [32].Next, we consider the cross-section of the grain boundary with the x = 0 plane (i.e., cross-section normal to the [¯110] direction); see the lower panel in Fig. 2 (the outer boundary ofthe projected grain boundary surface). Initially, the cross-section is a circle, and it graduallychanges to an ellipse as it shrinks during the evolution. This shows that the velocity in therotation axis direction, i.e. z direction is larger than that in the x and y directions. Thereason for this is that there is no constraint of Frank’s formula in the z direction which is thedirection of the rotation axis, and the velocity at the two poles on the grain boundary withrespect to the z direction (where the grain boundary is pure twist) is the same as that inthe curvature flow; whereas the velocity components in the x and y directions are adjustedfrom those in the curvature flow by the constraint of the Frank’s formula, and the resultingvelocity in the xy plane are depending on the entire grain boundary through the coefficient21 θ , as discussed above.Evolution of this initially spherical grain boundary and its dislocation structure, espe-cially the property that the shrinkage of the grain boundary is faster in the direction of therotation axis than in other directions, agree with the results of atomistic simulations usingphase field crystal models [28, 29]. Notice that with respect to the rotation axis of the z direction, at the two poles on the grain boundary, the grain boundary is pure twist, and theconstituent dislocations are screw dislocations; along the equator of the grain boundary, thegrain boundary is pure tilt, and the dislocations are edge dislocations. It was suggested inRef. [28] that such an anisotropic motion may be due to the difference in dislocation densitiesor that in the mobilities of screw and edge dislocations. Here, our continuum model providesa further explanation that this anisotropic motion of the grain boundary (and accordingly,the anisotropic motion of the screw and edge portions of the constituent dislocations) isdue to the constraint of Frank’s formula in order to maintain an equilibrium dislocationstructure.Fig. 3(a) shows the change of misorientation angle θ during the evolution, which is con-tinuously increasing. This behavior agrees with Cahn-Taylor theory [21], three dimensionalphase-field crystal simulations [28], and two-dimensional atomistic [20, 25], phase field crys-tal [26], and continuum [31, 32] simulations. Such increasing of misorientation angle cannotbe obtained by the classical motion by mean curvature models or pure sliding models, inwhich the misorientation angle is constant or is decreasing during the evolution.Fig. 3(b) shows evolution of the area of the grain boundary, which reveals the relation: S A ( t ) S A = 1 − At, (62)where A is some constant, and S A ( t ) and S A are the grain boundary area at time t andthat of the initial configuration, respectively. This agrees with the results of nearly lineardecrease of the grain boundary area using phase field crystal model for grain boundaries inboth fcc and bcc crystals [29]. The phase field crystal simulations in Ref. [28] showed thatthe decrease of the volume of the grain enclosed by an initially spherical low angle grainboundary in a bcc crystal approximately follows the relation V / ( t ) V / = 1 − A t , where A is22 Time(/M d ) M i s o r i en t a t i on ang l e ( ° ) (a) Time(/M d ) A r ea c hange S A / S A (b) d )00.20.40.60.8 D i s l o c a t i on den s i t y ( b - ) (c) d )020406080 D i s l o c a t i on l eng t h ( b ) (d)Figure 3: Shrinkage of an initially spherical grain boundary in fcc under pure coupling motion. The rotationaxis is the z direction ([111]), and the initial misorientation angle θ = 5 ◦ . (a) Evolution of misorientationangle θ . (b) Evolution of grain boundary area S A , where S A is the area of the initial grain boundary. (c)Evolution of density of dislocations with Burgers vector b (1) / b (2) / b (3) on the grain boundary. (d) Evolutionof the total length of dislocations with Burgers vector b (1) / b (2) / b (3) on the grain boundary. In (c) and (d),the densities and total lengths of dislocations with these three Burgers vectors are almost identical. V ( t ) and V are the volume of the grain enclosed by the grain boundary attime t and that of the initial configuration, respectively. It was argued in Ref. [28] that theirresults are consistent with the result of classical Von Neumann-Mullins relation [3] for a twodimensional grain boundary driven by curvature with constant energy density, i.e., Eq. (62)if S A denotes the area enclosed by the grain boundary in two dimensions, considering theapproximate relation V / ∼ S A . In this sense, simulation results using our continuum modeland the phase field crystal simulations in [29] are consistent with the results in Ref. [28] aswell as the result of the classical Von Neumann-Mullins relation. The nearly linear decreaseof the grain boundary area in Eq. (62) obtained by our continuum model and the phasefield crystal model in Ref. [29] is also in consistent with the result that the area enclosed bya two dimensional grain boundary is linearly decreasing in the two dimensional phase fieldcrystal simulations for circular grain boundaries [26] and continuum model simulations forcircular [31] and general shape [32] grain boundaries in two dimensions.Evolutions of dislocation densities on the grain boundary and total length of dislocationsare shown in Figs. 3(c) and (d). It can be seen from Fig. 3(c) that the densities of thedislocations with all the three Burgers vectors are increasing during the evolution. This isconsistent with the increase of misorientation angle θ during the evolution. Recall that inthe two dimensional case with dislocation conservation, the inter-dislocation distance is de-creasing as the grain boundary shrinks, which leads to the increases of both the dislocationdensities and the misorientation angle during the evolution [21]. The total length of dislo-cations is decreasing during the evolution as shown in Fig. 3(d). This is in agreement withthe phase field crystal simulation results in Ref. [28]. Unlike in the two dimensional casewith dislocation conservation [20, 21, 25, 26, 31, 32], in three dimensions without dislocationreaction, the grain boundary and its constituent dislocations (which are closed loops insteadof infinite, straight lines in the two dimensional case) are moving freely in the direction ofthe rotation axis, although the motion is subject to Frank’s formula in a direction normal tothe rotation axis. As a result, although along the equator of the grain boundary, the totalnumber of dislocations does not change during the evolution, all the dislocation loops areshrinking and the total length of dislocations is decreasing as the grain boundary shrinks.24 .1.2. Motion with dislocation reaction Now we perform simulations using our continuum model considering dislocation reaction,i.e. M r (cid:54) = 0. Dislocation reaction leads to removal of dislocations, resulting in the couplingmotion of the grain boundary [20, 21, 25, 28, 31, 32]. We use the same initial spherical grainboundary as in Sec. 4.1.1 without dislocation reaction. (a) (b) (c) (d)Figure 4: Shrinkage of an initially spherical grain boundary in fcc with dislocation reaction: M r b /M d =1 . × − . The rotation axis is the z direction ([111]), and the initial misorientation angle θ = 5 ◦ .The upper panel of images show the three-dimensional view of the grain boundary during evolution. Themiddle panel of images show the grain boundary during evolution viewed from the + z direction ([111]),and the lower panel of images show the grain boundary during evolution viewed from the + x direction([¯110]). Dislocations with Burgers vectors b (1) , b (2) and b (3) are shown by blue, black and red lines,respectively. Length unit: b . (a) The initial spherical grain boundary. (b), (c), and (d) Configurations attime t = 10 /M d µ, /M d µ, /M d µ , respectively. Fig. 4 shows the shrinkage of the initially spherical grain boundary with dislocation25eaction, where the reaction mobility M r b /M d = 1 . × − . We consider the cross-sectionof the grain boundary with the z = 0 plane (i.e., cross-section normal to the [111] rotationaxis), which is the equator of the grain boundary in the three dimensional view in the upperpanel in Fig. 4 and the outer curve in the view from the + z axis in the second panel inFig. 4. Along this curve, the grain boundary is pure tilt everywhere, and we have v ∗ = 0,i.e., the velocity is always in the z = 0 plane during the evolution. The evolution of thiscurve is similar to that of the two-dimensional grain boundary discussed in [31, 32]. Theinitial circular cross-section gradually changes to a hexagonal shape as it shrinks. Eachedge in this hexagon is pure tilt that consists of dislocations of only one Burgers vector.This behavior is consistent with the fact that the energy density of the grain boundary isanisotropic and the pure tilt boundary has the minimum energy of all tilt boundaries, andis the same as the evolution of two dimensional grain boundary with dislocation reactionobtained in [31, 32].The lower panel of Fig. 4 shows the evolution of the grain boundary in the view from the+ x direction ([¯110] direction). The cross-section of the grain boundary with the x = 0 planegradually changes to an ellipse as it shrinks, meaning that the grain boundary moves faster inthe rotation axis direction than in a direction normal to the rotation axis. This behavior andthe underlying reason are the same as those in the simulation without dislocation reactionshown and discussed in Sec. 4.1.1.These behaviors of the evolution of the initially spherical grain boundary with dislocationreaction in an fcc crystals are similar to the phase field crystal simulation results of an initiallyspherical grain boundary in a bcc crystal [28].Evolution of the misorietation angle θ with different values of reaction mobility M r isshown in Fig. 5(a). When M r (cid:54) = 0, the evolution of misorientation angle is controlled by boththe coupling effect and sliding effect. As can be seen from Fig. 5(a), the misorientation angle θ is increasing during the evolution except for the case with very high dislocation reactionmobility; as the dislocation reaction mobility M r increases, meaning the sliding effect dueto dislocation reaction is becoming stronger, the increase rate of θ decreases, and when thesliding effect is strong enough, the misorientation angle θ is decreasing. These properties26 Time(/M d ) M i s o r i en t a t i on ang l e ( ° ) (a) Time(/M d ) A r ea c hange S A / S A (b) d )00.20.40.60.8 D i s l o c a t i on den s i t y ( b - ) (c) d )020406080 D i s l o c a t i on l eng t h ( b ) (d)Figure 5: Shrinkage of an initially spherical grain boundary in fcc with different values of reaction mobility M r . The rotation axis is the z direction ([111]), and the initial misorientation angle θ = 5 ◦ . The reactionmobility M r b /M d = 0, 9 . × − , 1 . × − , and 3 . × − from the top curve to the bottom onein (a), (c) and (d), and from the bottom to the top ones in (b). (a) Evolution of misorientation angle θ .(b) Evolution of grain boundary area S A , where S A is the area of the initial grain boundary. (c) Evolutionof density of dislocations with Burgers vector b (1) / b (2) / b (3) on the grain boundary. (d) Evolution of thetotal length of dislocations with Burgers vector b (1) / b (2) / b (3) on the grain boundary. In (c) and (d), thedensities and total lengths of dislocations with these three Burgers vectors are almost identical. θ during the evolution, and the sliding motion generated by dislocation reaction willdecrease θ .Fig. 5(b) shows the evolution of grain boundary area with different values of dislocationreaction mobility M r . It can be seen that except for the case with very high dislocationreaction mobility, the decrease of grain boundary area still follows the linear law in Eq. (62),and is almost unchanged with different values of dislocation reaction mobility. In the casewith very high dislocation reaction mobility M r = 3 . × − M d /b , the decrease of grainboundary area starts to deviate from the linear law with slower deceasing rate, which is dueto the resulting significant decrease in the grain boundary energy density that slows downthe shrinking of the grain boundary. Again, the linear decrease of grain boundary area isconsistent with the available phase field crystal simulation results [28, 29] and results oftwo dimensional simulations using different methods as discussed in Sec. 4.1.1 for the casewithout dislocation reaction.Evolutions of dislocation densities on the grain boundary and total length of dislocationswith different values of reaction mobility M r are shown in Figs. 5(c) and (d). As can be seenfrom Fig. 5(c), the densities of the dislocations with all the three Burgers vectors are increas-ing during the evolution except for the case with very high dislocation reaction mobility;as the dislocation reaction mobility M r increases, the increase rate of dislocation densitiesdecreases, and when the dislocation reaction mobility is high enough, the dislocation den-sities are decreasing. These behaviors are consistent with the increase of misorientationangle θ during the evolution shown in Fig. 5(a). Fig. 5(d) shows that the total length ofdislocations is decreasing as the grain boundary shrinks, and the decrease rate is higher forhigher dislocation reaction mobility M r . The decrease of the total length of dislocations isin agreement with the phase field crystal simulation results in Ref. [28].Note that detailed dislocation reaction mechanisms on a curved grain boundary in threedimensions have been analyzed in Ref. [28] based on their phase field crystal simulations.There are also mechanisms observed/proposed based two dimensional molecular dynamics28imulations [20, 25]. In general, the small dislocation loops and segments of dislocationnetworks on a curved grain boundary in three dimensions are easier to react than the infinite,straight dislocations in the two dimensional cases. Now we consider grain boundaries in bcc Fe. We choose the directions [100], [01¯1] and[011] to be the x , y and z directions, respectively. There seven possible Burgers vectors,and in this coordinate system, they are b (1) = (cid:16) , , √ (cid:17) a , b (2) = (cid:16) , √ , (cid:17) a , b (3) = (cid:16) , − √ , (cid:17) a , b (4) = (cid:16) − , , √ (cid:17) a , b (5) = (1 , , a , b (6) = (cid:16) , √ , √ (cid:17) a , and b (7) = (cid:16) , − √ , √ (cid:17) a , where a is the lattice constant. For bcc Fe, a= 0 . ν = 0 .
29. The rotation axis is a = (0 , , R = 20 b and misorientation angle is θ = 4 ◦ . There are three sets of dislocations withBurgers vectors b (2) , b (3) , b (5) , respectively, in the equilibrium dislocation structure onthis initial, spherical grain boundary; see the top image in Fig. 6(a). Parametrization anddiscretization of the grain boundary are the same as those for the grain boundary in fccdescribed in Sec. 4.1. We first consider the motion of this grain boundary without dislocation reaction, i.e. thereaction mobility M r = 0 in Eq. (54) (and Eq. (43)). This is the pure coupling motion.Fig. 6 shows the shrinkage of this spherical grain boundary in bcc under the pure couplingmotion. Evolutions of the grain boundary and its dislocation structure are similar to thoseof the grain boundary in fcc discussed in Sec. 4.1.1. Especially, the shrinkage of the grainboundary is faster in the direction of the rotation axis than in other directions, which agreeswith the results of atomistic simulations using phase field crystal models [28, 29]. Again,this anisotropic motion can be explained based on our continuum model by the constraintof Frank’s formula in any direction normal to the rotation axis. The shape-preservingevolution of the equator of the grain boundary (with respect to the rotation axis) agreeswith the results of the two-dimensional grain boundary dynamics models [20, 33, 26, 31, 32].29 a) (b) (c) (d)Figure 6: Shrinkage of an initially spherical grain boundary in bcc under pure coupling motion, i.e., withoutdislocation reaction. The rotation axis is the z direction ([011]), and the initial misorientation angle θ = 4 ◦ .The upper panel of images show the three-dimensional view of the grain boundary during evolution. Themiddle panel of images show the grain boundary during evolution viewed from the + z direction ([011]),and the lower panel of images show the grain boundary during evolution viewed from the + x direction([100]). Dislocations with Burgers vectors b (2) , b (3) and b (5) are shown by blue, black and red lines,respectively. Length unit: b . (a) The initial spherical grain boundary. (b), (c), and (d) Configurations attime t = 3 /M d µ, /M d µ, /M d µ , respectively. Time(/M d ) M i s o r i en t a t i on ang l e ( ° ) (a) Time(/M d ) A r ea c hange S A / S A (b) d )00.20.40.60.8 D i s l o c a t i on den s i t y ( b - ) b (2) b (3) b (5) (c) d )020406080 D i s l o c a t i on l eng t h ( b ) b (2) b (3) b (5) (d)Figure 7: Shrinkage of an initially spherical grain boundary in bcc under pure coupling motion. The rotationaxis is the z direction ([011]), and the initial misorientation angle θ = 4 ◦ . (a) Evolution of misorientationangle θ . (b) Evolution of grain boundary area S A , where S A is the area of the initial grain boundary. (c)Evolution of densities of dislocations on the grain boundary. These dislocations have Burgers vectors b (2) , b (3) , and b (5) on the grain boundary. (d) Evolution of the total lengths of dislocations with these Burgersvectors. The densities and total lengths of dislocations with the Burgers vectors b (2) , b (3) , and b (5) areshown by blue lines, black dots, and red lines, respectively, in (c) and (d). θ , linear decrease of the grainboundary area, increase of densities of dislocations, and decrease of the total lengths ofdislocations during the shrinkage of the grain boundary. These results are similar to those ofthe grain boundary in fcc discussed in Sec. 4.1.1 and are consistent with available phase fieldcrystal and two dimensional simulation results; see the discussion there. Note that for thedynamics of this grain boundary in bcc, the densities and total lengths of dislocations withBurgers vectors b (2) and b (3) are almost identical, and are greater than those of dislocationswith Burgers vector b (5) . Recall that the lengths of Burgers vectors b (2) and b (3) are equal,and are smaller than the length of Burgers vector b (5) . Now we perform simulations using our continuum model considering dislocation reaction,i.e. M r (cid:54) = 0, for the motion of a grain boundary in bcc Fe. We use the same initial sphericalgrain boundary as in Sec. 4.2.1 without dislocation reaction.Fig. 8 shows the shrinkage of the initially spherical grain boundary with dislocationreaction, where the reaction mobility M r b /M d = 7 . × − . As in the evolution of thegrain boundary in fcc with dislocation reaction shown in Sec. 4.1.2, the equator of the grainboundary (with respect to the rotation axis) gradually evolves from a circle into a hexagon.However, unlike in fcc, here the grain boundary in bcc evolves into a non-regular hexagon.This is due to the fact that here the lengths of Burgers vectors b (2) and b (3) are smallerthan the length of Burgers vector b (5) , while in the grain boundary in fcc, the three Burgersvectors have the same length. Again, we can see that the shrinkage of the grain boundaryis faster in the direction of the rotation axis than in other directions. These results agreewith the those of phase field crystal simulations [28].Evolution of the misorietation angle θ with different values of reaction mobility M r isshown in Fig. 9(a). When M r (cid:54) = 0, the evolution of misorientation angle is controlled by boththe coupling effect (which is associated with the conservation of dislocations and increases θ ) and sliding effect (which is associated with dislocation reaction and decreases θ ). Whilethe misorientation angle θ is increasing during the evolution, the increase rate of θ decreases32 a) (b) (c) (d)Figure 8: Shrinkage of an initially spherical grain boundary in bcc with dislocation reaction: M r b /M d =7 . × − . The rotation axis is the z direction ([011]), and the initial misorientation angle θ = 4 ◦ . Theupper panel of images show the three-dimensional view of the grain boundary during evolution. The middlepanel of images show the grain boundary during evolution viewed from the + z direction ([011]), and the lowerpanel of images show the grain boundary during evolution viewed from the + x direction ([100]). Dislocationswith Burgers vectors b (2) , b (3) and b (5) are shown by blue, black and red lines, respectively. Length unit: b .(a) The initial spherical grain boundary. (b), (c), and (d) Configurations at time t = 3 /M d µ, /M d µ, /M d µ ,respectively. Time(/M d ) M i s o r i en t a t i on ang l e ( ° ) (a) Time(/M d ) A r ea c hange S A / S A (b) d )00.511.52 D i s l o c a t i on den s i t y ( b - ) (c) d )050100150 D i s l o c a t i on l eng t h ( b ) (d)Figure 9: Shrinkage of an initially spherical grain boundary in bcc with different values of reaction mobility M r . The rotation axis is the z direction ([011]), and the initial misorientation angle θ = 4 ◦ . The reactionmobility M r b /M d = 0, 7 . × − , 1 . × − , and 2 . × − from the top curve to the bottom onein (a), (c) and (d), and from the bottom to the top ones in (b). (a) Evolution of misorientation angle θ .(b) Evolution of grain boundary area S A , where S A is the area of the initial grain boundary. (c) Evolutionof the density of all the dislocations on the grain boundary. (d) Evolution of the total length of all thedislocations on the grain boundary.
34s the dislocation reaction mobility M r increases. Fig. 9(b) shows the evolution of grainboundary area with different values of dislocation reaction mobility M r . The decrease ofgrain boundary area still follows the linear law in Eq. (62) except in the later stage ofthe evolution with high dislocation reaction mobility, and is slower for higher dislocationreaction mobility. Figs. 9(c) and (d) show the evolutions of dislocation densities on thegrain boundary and total length of dislocations with different values of reaction mobility M r . The density of the dislocations on the grain boundary is increasing whereas the totallength of dislocations is decreasing for these values of dislocation reaction mobility M r .As the dislocation reaction mobility M r increases, the increase rate of dislocation densitydecreases, and the decrease rate of the total length of dislocations increases except for thelater stage of evolution without dislocation reaction. These behaviors are similar to those ofthe grain boundary in fcc shown in Fig. 5, and more discussion can be found there. Theseresults also agree with the available phase field crystal simulation results [28, 29].
5. Conclusions
We have developed a continuum model for the dynamics of grain boundaries in threedimensions that incorporates the motion and reaction of the constituent dislocations. Thecontinuum model includes evolution equations for both the motion of the grain boundary andthe evolution of dislocation structure on the grain boundary. The evolution of orientation-dependent continuous distributions of dislocation lines on the grain boundary is based onthe simple representation using dislocation density potential functions. This simple repre-sentation method also guarantees continuity of the dislocation lines on the grain boundariesduring the evolution. The critical but computationally expensive long-range elastic inter-action of dislocations is replaced by a projection formulation that maintains the constraintof the Frank’s formula describing the equilibrium of the strong long-range interaction. Thiscontinuum model is able to describe the grain boundary motion and grain rotation due toboth coupling and sliding effects, to which the classical motion by mean curvature modeldoes not apply.In order to overcome the ill-posedness in formulation that comes from the nonconvexity35f the energy density, we use the components of the surface gradients of the dislocationdensity potential functions instead of these functions directly. Relationship between thecomponents of these surface gradients (i.e. continuity of dislocation lines) is maintained bythe projection method during the evolution.Using the obtained continuum model, simulations are performed for the dynamics ofinitially spherical low angle grain boundaries in fcc Al and bcc Fe, under the conditionswithout dislocation reaction (pure coupling motion) and with dislocation reaction (withsliding motion). The simulations have shown increase of the misorientation angle as thegrain boundary shrinks under the effect of conservation of dislocations, anisotropic motionin the directions along and normal the rotation axis, anisotropic motion in the normalplane with respect to the rotation axis due to dislocation reaction, and linear decrease ofgrain boundary area. These results agree well with those of atomistic simulations (phase-field crystal simulations) [28, 29]. The simulation results are also consistent with previouslyobtained results using continuum model in two dimensions [31, 32]. In particular, we explainthe anisotropic motion in the directions along and normal the rotation axis by the fact thatthe constraint of Frank’s formula only has effect in a direction normal to the rotation axis,and the motion is free in the direction of the rotation axis.
Acknowledgement
This work was supported by the Hong Kong Research Grants Council General ResearchFund 16301720 (LCZ) and 16302818 (YX).
References [1] A. Sutton, R. Balluffi, Interfaces in Crystalline Materials, Clarendon Press, Oxford, 1995.[2] C. Herring, Surface tension as a motivation for sintering, in: W. E. Kingston (Ed.), The Physics ofPowder Metallurgy, McGraw-Hill, New York, 1951, pp. 143–179.[3] W. W. Mullins, Two-dimensional motion of idealized grain boundaries, J. Appl. Phys. 27 (1956) 900–904.
4] L.-Q. Chen, W. Yang, Computer simulation of the domain dynamics of a quenched system with alarge number of nonconserved order parameters: The grain-growth kinetics, Phys. Rev. B 50 (1994)15752–15756.[5] A. Kazaryan, Y. Wang, S. A. Dregia, B. R. Patton, Generalized phase-field model for computer simu-lation of grain growth in anisotropic systems, Phys. Rev. B 61 (2000) 14275–14278.[6] D. Kinderlehrer, C. Liu, Evolution of grain boundaries, Math. Models Methods Appl. Sci. 4 (2001)713–729.[7] M. Upmanyu, G. N. Hassold, A. Kazaryan, E. A. Holm, Y. Wang, B. Patton, D. J. Srolovitz, Boundarymobility and energy anisotropy effects on microstructural evolution during grain growth, Interface Sci.10 (2002) 201–216.[8] H. Zhang, M. Upmanyu, D. J. Srolovitz, Curvature driven grain boundary migration in aluminum:molecular dynamics simulations, Acta Mater. 53 (2005) 79–86.[9] D. M. Kirch, E. Jannot, L. A. Barrales-Mora, D. A. Molodov, G. Gottstein, Inclination dependence ofgrain boundary energy and its impact on the faceting and kinetics of tilt grain boundaries in aluminum,Acta Mater. 56 (2006) 4998–5011.[10] M. Elsey, S. Esedoglu, P. Smereka, Diffusion generated motion for grain growth in two and threedimensions, J. Comput. Phys. 228 (2009) 8015–8033.[11] E. A. Lazar, R. D. MacPherson, D. J. Srolovitz, A more accurate two-dimensional grain growth algo-rithm, Acta Mater. 58 (2010) 364–372.[12] J. C. Li, Possibility of subgrain rotation during recrystallization, J. Appl. Phys. 33 (1962) 2958–2965.[13] P. G. Shewmon, in: H. Margolin (Ed.), Recrystallization, grain growth and textures, American Societyof Metals, Metals Park, 1966, pp. 165–199.[14] K. Harris, V. Singh, A. King, Grain rotation in thin films of gold, Acta Mater. 46 (1998) 2623–2633.[15] R. Kobayashi, J. A. Warren, W. C. Carter, A continuum model of grain boundaries, Phys. D 140 (2000)141–150.[16] M. Upmanyu, D. J. Srolovitz, A. E. Lobkovsky, J. A. Warren, W. C. Carter, Simultaneous grainboundary migration and grain rotation, Acta Mater. 54 (2006) 1707–1719.[17] S. Esedoglu, Grain size distribution under simultaneous grain boundary migration and grain rotationin two dimensions, Comput. Mater. Sci. 121 (2016) 209–216.[18] Y. Epshteyn, C. Liu, M. Mizuno, Motion of grain boundaries with dynamic lattice misorientations andwith triple junctions drag, arXiv preprint arXiv:1903.11512.[19] C. H. Li, E. H. Edwards, J. Washburn, E. R. Parker, Stress-induced movement of crystal boundaries,Acta Metall. 1 (1953) 223–229.[20] S. G. Srinivasan, J. W. Cahn, Challenging some free-energy reduction criteria for grain growth, in: . Ankem, C. S. Pande, I. Ovid’ko, S. Ranganathan (Eds.), Science and Technology of Interfaces,TMS, Seattle, 2002, pp. 3–14.[21] J. W. Cahn, J. E. Taylor, A unified approach to motion of grain boundaries, relative tangential trans-lation along grain boundaries, and grain rotation, Acta Mater. 52 (2004) 4887–4898.[22] J. W. Cahn, Y. Mishin, A. Suzuki, Coupling grain boundary motion to shear deformation, Acta Mater.54 (2006) 4953–4975.[23] D. A. Molodov, V. A. Ivanov, G. Gottstein, Low angle tilt boundary migration coupled to sheardeformation, Acta Mater. 55 (2007) 1843–1848.[24] T. Gorkaya, D. A. Molodov, G. Gottstein, Stress-driven migration of symmetrical < > tilt grainboundaries in al bicrystals, Acta Mater. 57 (2009) 5396–5405.[25] Z. Trautt, Y. Mishin, Grain boundary migration and grain rotation studied by molecular dynamics,Acta Mater. 60 (2012) 2407–2424.[26] K. W. Wu, P. W. Voorhees, Phase field crystal simulations of nanocrystalline grain growth in twodimensions, Acta Mater. 60 (2012) 407–419.[27] K. McReynolds, K.-A. Wu, P. Voorhees, Grain growth and grain translation in crystals, Acta Mater.120 (2016) 264–272.[28] A. Yamanaka, K. McReynolds, P. W. Voorhees, Phase field crystal simulation of grain boundary motion,grain rotation and dislocation reactions in a bcc bicrystal, Acta Mater. 133 (2017) 160–171.[29] M. Salvalaglio, R. Backofen, K. Elder, A. Voigt, Defects at grain boundaries: A coarse-grained, three-dimensional description by the amplitude expansion of the phase-field crystal model, Phys. Rev. Mater.2 (2018) 053804.[30] B. B. Rath, M. Winning, J. C. M. Li, Coupling between grain growth and grain rotation, Appl. Phys.Lett. 90 (2007) 161915.[31] L. C. Zhang, Y. Xiang, Motion of grain boundaries incorporating dislocation structure, J. Mech. Phys.Solids 117 (2018) 157–178.[32] L. Zhang, Y. Xiang, A new formulation of coupling and sliding motions of grian boundaries based ondislocation structure, SIAM J. Appl. Math. 80 (2020) 2365–2387.[33] J. E. Taylor, J. W. Cahn, Shape accommodation of a rotating embedded crystal via a new variationalformulation, Interfaces and Free Boundaries 9 (2007) 493–512.[34] A. Basak, A. Gupta, A two-dimensional study of coupled grain boundary motion using the level setmethod, Modell. Simul. Mater. Sci. Eng. 22 (2014) 055022.[35] N. C. Admal, G. Po, J. Marian, A unified framework for polycrystal plasticity with grain boundaryevolution, Int. J. Plasticity 106 (2018) 1–30.[36] A. Ask, S. Forest, B. Appolaire, K. Ammar, O. U. Salman, A cosserat crystal plasticity and phase field heory for grain boundary migration, J. Mech. Phys. Solids 115 (2018) 167–194.[37] X. H. Zhu, Y. Xiang, Continuum framework for dislocation structure, energy and dynamics of disloca-tion arrays and low angle grain boundaries, J. Mech. Phys. Solids 69 (2014) 175–194.[38] F. C. Frank, The resultant content of dislocations in an arbitrary intercrystalline boundary, Office ofNaval Research, Pittsburgh, 1950, pp. 150–154.[39] B. A. Bilby, Bristol conference report on defects in crystalline materials, Phys. Soc., London (1955)123.[40] J. P. Hirth, J. Lothe, Theory of Dislocations, 2nd Edition, Wiley, New York, 1982.[41] Y. Xiang, X. Yan, Stability of dislocation networks of low angle grain boundaries using a continuumenergy formulation, Discrete Contin. Dyn. Syst. Ser. B 23 (2018) 2989.[42] A. J. Chorin, Numerical solution of the navier-stokes equations, Math. Comp. 22 (1968) 745–762.[43] X. X. Qin, Y. J. Gu, L. C. Zhang, Y. Xiang, Continuum model and numerical method for dislocationstructure and energy of grain boundaries, arXiv (2021) arXiv:2101.02596.[44] L. C. Zhang, Y. J. Gu, Y. Xiang, Energy of low angle grain boundaries based on continuum dislocationstructure, Acta Mater. 126 (2017) 11–24.[45] Y. Giga, Surface Eolution Equations, Birkhauser, Berlin, 2006.heory for grain boundary migration, J. Mech. Phys. Solids 115 (2018) 167–194.[37] X. H. Zhu, Y. Xiang, Continuum framework for dislocation structure, energy and dynamics of disloca-tion arrays and low angle grain boundaries, J. Mech. Phys. Solids 69 (2014) 175–194.[38] F. C. Frank, The resultant content of dislocations in an arbitrary intercrystalline boundary, Office ofNaval Research, Pittsburgh, 1950, pp. 150–154.[39] B. A. Bilby, Bristol conference report on defects in crystalline materials, Phys. Soc., London (1955)123.[40] J. P. Hirth, J. Lothe, Theory of Dislocations, 2nd Edition, Wiley, New York, 1982.[41] Y. Xiang, X. Yan, Stability of dislocation networks of low angle grain boundaries using a continuumenergy formulation, Discrete Contin. Dyn. Syst. Ser. B 23 (2018) 2989.[42] A. J. Chorin, Numerical solution of the navier-stokes equations, Math. Comp. 22 (1968) 745–762.[43] X. X. Qin, Y. J. Gu, L. C. Zhang, Y. Xiang, Continuum model and numerical method for dislocationstructure and energy of grain boundaries, arXiv (2021) arXiv:2101.02596.[44] L. C. Zhang, Y. J. Gu, Y. Xiang, Energy of low angle grain boundaries based on continuum dislocationstructure, Acta Mater. 126 (2017) 11–24.[45] Y. Giga, Surface Eolution Equations, Birkhauser, Berlin, 2006.