A Cosserat–phase-field theory of crystal plasticity and grain boundary migration at finite deformation

In metallic polycrystals, an important descriptor of the underlying microstructure is the orientation of the crystal lattice of each grain. During thermomechanical processing, the microstructure can be significantly altered through deformation, nucleation of new subgrains and grain boundary migration. Cosserat crystal plasticity provides orientation as a degree of freedom and is therefore a natural choice for the development of a coupled framework to deal with concurrent viscoplasticity and grain growth. In order to take into account grain boundary motion, the Cosserat theory is adapted with inspiration from orientation phase-field models. This allows for the microstructure at a material point to evolve on the one hand due to deformation-induced lattice reorientation and on the other hand due to a sweeping grain boundary. With a proper separation of plastic evolution in the bulk of the grain and in the grain boundary, the model can successfully capture grain boundary migration due to lattice curvature and due to statistically stored dislocations.


Introduction
Microstructure evolution during thermomechanical processing of metals is largely due to viscoplastic deformation together with grain boundary (GB) migration and nucleation of new subgrains. An important descriptor of the grain structure at the microscale is the crystal lattice orientation. Grains of various orientations are separated by boundaries which consist of a few atom layers that are distinct from the ordered structure inside the bulk of the grains [67]. Viscoplastic deformation may result in significant reorientation of the crystal lattice and in the development of heterogeneous orientation distributions within grains [37,57]. The energy stored in the lattice structure during deformation in the form of statistically stored or geometrically necessary dislocations acts as a driving force for grain boundary migration [33]. GB migration is necessarily associated with lattice reorientation when a grain boundary sweeps a material point. Crystal plasticity at finite deformation has been recently combined with a vertex GB migration model in [52] and with a phase-field recrystallization model in [73]. Both models are used successively, the mechanical model providing the initial conditions of the phase-field part. The present work is dedicated to fully coupled mechanics-phase-field models.
With the view to modeling thermomechanical processing of metallic materials, both reorientation due to viscoplastic deformation and GB migration should be considered. A Cosserat theory [20,22,24] is a natural choice for such a coupled approach as it contains rotation degrees of freedom accounting for microstructure evolution. Each material point in a Cosserat theory is endowed with a triad of directors which can be associated with the lattice crystal directions by means of suitable internal constraints [27,28,50,51]. Plastic deformation inside the grains is due to dislocation glide on preferred slip systems that can be computed by flow rules according to standard crystal plasticity [17,60]. The slip processes result in a relative rotation between the lattice directions and the material lines. Non-homogeneous plastic deformation gives rise to lattice rotation gradients, i.e., lattice curvature that must be accommodated by so-called geometrically necessary dislocations (GNDs) [6]. The Cosserat theory provides lattice curvature as a deformation measure in the form of the torsioncurvature tensor which can be identified as an essential part of the dislocation density tensor [44,45,55,68].
In a recent work [12], a Cosserat crystal plasticity approach together with a level-set method was used to model dynamic recrystallization. The Cosserat directors were associated with the lattice orientations using a penalty method as in [11,27,51]. However, the proposed methodology represents rather a sequential algorithmic procedure and cannot be regarded as a unified field theory of viscoplastic deformation and GB kinematics. A fully coupled theory was recently proposed in [7] for small deformations, combining Cosserat crystal plasticity with a phase-field model due to [42,74]. The Kobayashi-Warren-Carter (KWC) phase-field approach considers the crystal orientation as a phase field with its own relaxational dynamics. Unlike frequently used multiphasefield models [26,65,66], the KWC model is able to accommodate the evolution of heterogeneous lattice orientation inside a grain without the need to introduce new phases (orientations) in the reoriented regions. In addition to the orientation field, the KWC model contains a phase field which is interpreted in the case of solid-to-solid transformations as a measure of crystalline order. In [4], a coupled model for large deformation crystal plasticity was formulated inspired by the KWC approach. The approach is very similar to the one in [7] and the one proposed in this work. Instead of working with the torsion-curvature tensor of the Cosserat theory, the energy density in [4] contains corresponding contributions due to the full dislocation density tensor. Evolution of orientation is assumed to be due to plastic slip processes also in the case of grain boundary migration. The model proposed in [4] was applied to grain boundary sliding as well as curvature driven GB migration and grain rotation. These approaches strongly differ from multiphase-field models of GB migration which assign one phase-field variable to each initial grain orientation, thus requiring the introduction of new phase fields at each material point inside a grain to allow for heterogeneous lattice rotation to occur. That is why the recent multiphase-field models [38,77] rely on a Taylor assumption for the proposed coupling between phase field and crystal plasticity.
The main objective of the present work is, first, to extend the Cosserat-phase-field model proposed in [7] to the finite deformation framework and, second, to provide finite element simulation of dynamic recovery and GB migration, also not presented in the Ref. [7]. The free energy density function of the model contains contributions from geometrically necessary dislocations in the form of the torsion-curvature tensor and comprises terms coming from the KWC model. The Cosserat directors are identified with crystal lattice directions. This results in the description of grain boundaries as diffuse interfaces. In the bulk of the grains, lattice orientation changes are due to plastic slip processes according to classical crystal plasticity. Inside the grain boundary, an additional driving force term for plastic spin is introduced in order to differentiate bulk and GB behavior.
In the original KWC phase-field model, the energy of the system can be lowered by bulk grain rotation in the Allen-Cahn evolution equation for lattice orientation [42]. This is mainly due to the absence of mechanical coupling in the formulation. There is some evidence supporting gross rotation toward preferential misorientations in nanograins (cf [16,70,72] and references therein), but in larger micron-sized grains such gross rotation of the crystal lattice is not observed. It was proposed in the literature to limit the bulk rotation in the KWC model by distinguishing the kinematic coefficients for the orientation evolution in the bulk of the grain and the grain boundary [48]. In the coupled approach proposed here, it is expected that the Cosserat coupling can instead be used to control and even suppress the spurious bulk grain rotation. This is an original contribution of the present work allowing to spare one constitutive evolution law coming from the KWC model whose physical interpretation was still unclear.
The presentation is structured as follows. In Sect. 2, the general framework of the model is established. The large deformation kinematics is introduced, and balance equations and general formats for the constitutive and evolution equations are derived from energy principles. While the general framework is applicable between the nano-and micron scales, the constitutive choices given in Sect. 3 are made with a view to dynamic recrystallization simulations in polycrystalline aggregates on the micron scale. The corresponding linearization of the theory in Sect. 4 allows to recover the previous small strain and small rotation formulation from [7]. It is also compared with the governing equations of the KWC model. In Sect. 5, the proposed theory is applied numerically to grain boundary migration due to curvature and grain boundary migration due to the buildup of statistically stored dislocations (SSDs) during plastic deformation. The influence of the Cosserat coupling on the evolution of lattice orientation and in particular bulk grain rotation is studied. In particular, the essential role of the skew-symmetric part of the stress tensor as a driving force for lattice reorientation in the GB region is highlighted.

Cosserat-phase-field framework
This section is dedicated to the description of the kinematics and order parameters, the derivation of the balance laws and the general formulation of the constitutive laws within the framework of large deformations.

Notation
In what follows, vectors a i are denoted by a and second-order tensors A i j by A ∼ and fourth-order tensors C i jkl are written as C ≈ . In some cases, both the index free and indexed versions of an expression are given. In this case, indices in uppercase indicate material quantities and lowercase indices indicate spatial quantities. The summation convention then only holds for indices of the same case, i.e., A i I is not the trace of the mixed tensor A ∼ . I ∼ is the second-order identity tensor. Quite frequently, skew-symmetric tensors will be represented by their respective pseudo-vectors (denoted by a superposed cross) according to and likewise, the skew-symmetric tensor can be found from the pseudo-vector through where ∼ is the third-order Levi-Civita permutation tensor i jk . Double contraction A i j B i j will be denoted by a colon A ∼ : B ∼ and simple contraction, and a i b i is written as the usual dot product a · b. The tensor product A ∼ = a ⊗ b indicates the construction A i j = a i b j . The (material) gradient ∂a i /∂ X J and (spatial) divergence ∂a i /∂ x i of a vector a are written as a ⊗ ∇ 0 and ∇ · a, respectively. The (material) divergence of a tensor is given by A ∼ · ∇ 0 with differentiation acting on the second index, i.e., ∂ A i J /∂ X J . Differentiation with respect to the current configuration ∂/∂ x i is denoted by ∇.

Cosserat kinematics, deformation measures and order parameter
In the framework of finite deformations, the position of a material point in the reference (material) configuration B 0 at time t = t 0 is described by the vector X and its position in the current (spatial) configuration B t at time t is given by x = ϕ(X, t). The deformation gradient tensor can then be introduced as the material gradient of the map ϕ such that The volumetric part of the deformation is then given by where ρ 0 and ρ are the mass densities in the reference and deformed configurations, respectively. The translational degrees of freedom of each material point are given by the displacement vector u = x − X. In a Cosserat continuum, each material point is endowed with additional rotational degrees of freedom. The Cosserat rota- 3 of a triad of orthogonal lattice vectors attached to a material point between the initial state and a stress released state at time t. The rotation tensor is orthogonal meaning that Furthermore, R ∼ can be described in terms of the vector field Θ (the associated Lie group), containing the three rotational degrees of freedom, through the relation where δ j J is the shifter (c.f. [22,23]). By defining it is possible to derive δΘ as where It can be noted that δΘ is not a total differential and, unlike F ∼ , Γ ∼ is not generally invertible. The Cosserat deformation tensor F ∼ and the wryness or torsion-curvature tensor Γ ∼ are given by where • denotes the frame of the attached directors, which transform with the Cosserat microrotation. These relative measures are frame invariant under Euclidean transformations [39] and therefore make up suitable strain measures for the development of constitutive relations. An equivalent expression for the wryness tensor is [22,28] The velocity and gyration tensors are defined by v =u, The time derivatives of the Cosserat strains can be expressed in terms of the (spatial) gradients of these quantities as where the gyration vector × ω R (note that R is not an index) can be found from The relative velocity gradient v ⊗ ∇ − ω ∼ R in Eq. (13) describes the local motion of a material element with respect to the microstructure.
In order to study moving grain boundaries, the Cosserat model is enhanced with an additional degree of freedom in the form of a phase-field variable φ ∈ [0, 1] and its gradient ∇φ = F ∼ −T · ∇ 0 φ. The variable φ can be considered as a coarse-grained or macroscopic order parameter where φ = 1 in the bulk of a grain of perfect crystalline order and φ < 1 in the diffuse grain boundary where crystal order is lower. In the presence of plastic deformation, the requirement φ = 1 in the bulk of the grain must be relaxed, allowing for the disordering of the crystal lattice through dislocation densities building up and associated stored energy [1]. The study is restricted to single-phase materials. A separate orientation phase-field variable is not introduced. Rather, the Cosserat directors will be associated with the lattice orientations such that the degrees of freedom Θ can be taken to represent the crystalline orientations. The KWC phase-field model takes into consideration the energy of grain boundaries through the lattice curvature (the gradient of the orientation phase-field variable). In the present model, the torsion-curvature tensor is used instead resulting in a diffuse interface representation of the Cosserat orientation field.

Balance laws
The method of virtual power is used to introduce the generalized stresses of the theory and derive the corresponding balance laws. The virtual motions are the velocity v and the gyration rate × ω R introduced in the previous section. In order to fulfill the requirements of frame invariance, only v ⊗∇ −ω ∼ R and × ω R ⊗∇, together with their duals, appear in the formulation (c.f. [28]). The set of virtual field variables, including those due to the phase-field φ, entering the principal of virtual power is then given by The dual quantities for the phase field and its gradient are the microstresses π φ and ξ φ , respectively, according to the formalism in [34]. For the mechanical part, the dual quantities are the stress tensor σ ∼ and the couple stress tensor m ∼ . It should be noted that σ ∼ is in general not symmetric and thus does not represent the usual Cauchy stress tensor. The virtual power density of internal forces is then given by The virtual power densities due to external forces and contact forces are given by where f , c, π ext φ , respectively, are volume forces and couples and a volume source term, and where t is the traction vector, m c is a surface density of couples and π c φ the generalized traction for the phase field. The principle of virtual power states that over any region D with boundary ∂D of B t and for all virtual fields. After integrating by parts and applying the divergence theorem, the principle of virtual power can then be written for the coupled Cosserat-phase-field problem as where n is the outward unit normal to ∂D. By the requirement that (21) should hold over any region D and ∀(φ, v, × ω R ), the local forms of the balance laws and boundary conditions follow directly and are given by and respectively.

Constitutive and evolution equations
The Clausius-Duhem inequality states that 2 where Ψ is the Helmholtz free energy density. Inserting the density of internal power (17) gives Multiplying the above expression with J and taking note of Eqs. (13) and (14) together with (4), the Clausius-Duhem inequality can equivalently be written as where the transformed stresses are given by To extend the treatment to viscoplastic behavior, elastic and plastic strain measures must be defined. A multiplicative split of the Cosserat deformation gradient is adopted [28,49], such that 3 This decomposition is illustrated in Fig. 1. Volumetric measures related to the respective deformations can be introduced according to 2 Isothermal conditions are assumed for a more compact presentation. 3 The multiplicative split of the deformation gradient is assumed here to be of the form F ∼ = F ∼ e · F ∼ p so that the elastic contribution can be related to the elastic contribution of the Cosserat deformation according to F ∼ e = R ∼ T · F ∼ e . The decomposition is made unique (up to a symmetry transformation of the crystal) by assuming an isoclinic intermediate configuration according to [49]. Consider the polar decomposition F ∼ e = R ∼ e · U ∼ e . In the case where the microrotation is required to follow the elastic rotation so that R ∼ = R ∼ e , the intermediate configuration coincides with the frame of the microstructure if the elastic stretches are small, i.e., U ∼ e ≈ I ∼ . The corresponding partition of the wryness tensor is taken to be based on the arguments given in [21,28,61,64]. Adopting decompositions (35) and (37), rather than a purely additive decomposition of the curvature Γ ∼ , enables the definition of a released intermediate configuration where both force and couple stresses are removed. The deformation rates appearing in Eq. (30) can be written in terms of the energetic and dissipative contributions according to Inserted into (30),this gives is the sum of the Mandel-type stress tensors defined as and It is now assumed that the Helmholtz free energy takes the form where r α are internal variables related to the viscoplastic behavior (later they will be associated with the densities of statistically stored dislocations for all N slip systems α = 1, . . . , N ). Following the approach by Gurtin [34] or [35], in order to account for the relaxational behavior of a phase-field model, the stress measure Π φ associated with the phase-field φ is assumed to have equilibrium and non-equilibrium contributions in analogy with a viscous Kelvin element, such that The Clausius-Duhem inequality now takes the form where the thermodynamic force associated with the internal variable r α is given by The terms on the first two lines of (47) are assumed to be purely energetic. Setting each of them in turn to zero yields the constitutive relations The remaining terms on the bottom two lines of (47) are assumed to be dissipative. It is postulated that there exists a dissipation potential with contributions related to each term, such that The intrinsic dissipation inequality can finally be written as

Alternative formulation of constitutive equations
The constitutive relations derived above from (47) contain a mix of mechanical quantities defined in the local frame of the microstructure and Lagrangian variables for the phase-field part, defined in the referential configuration B 0 . An alternative formulation with the phase-field part defined in the intermediate configuration is possible [47]. 4 To this end, a quantity K φ is now introduced in the intermediate configuration through appropriate pull-back/push-forward operations on the gradients of φ, such that A generalized stress ξ p that is work conjugate to K φ can be found according to the requirement that the work generated by ξ φ over ∇φ on a volume element dv in B t should be equal to that generated by ξ p over K φ on a (fictitious) volume element dV p in the released state, i.e., where The Clausius-Duhem inequality (30) can now be written in terms of the newly introduced, intermediate quantities by noting that (62) Assuming now a free energy density of the format (to be compared to 45) the dissipation inequality can be written as whereΠ According to (64), state laws (49) to (51) still hold, but Eq. (52) is replaced by Usual finite deformation phase-field theories mainly rely on a Lagrangian or Eulerian format for the phasefield gradient [13,14]. The framework proposed in this subsection is a new formulation of phase-field models coupled with plasticity compared to existing ones [9,18,62,63]. It is also necessary to define the dissipation potential for the plasticity in terms of the modified Mandel-type stressΠ ∼ M so that

Explicit choices for the free energy and dissipation potentials
The system of equations is closed in this section by explicit expressions for the Helmholtz energy density and the dissipation potentials of the model.

Helmholtz free energy
The following free energy function is adopted The three lines in the previous equation, respectively, contain contributions related to elasticity, GB and GND energy and SSD-based energy. For simplicity, the free energy density is taken as an isotropic function of its arguments. Anisotropy for the linear Cosserat theory has been treated in [40]. The material parameters for the elastic behavior are the Lamé parameters μ and λ together with the Cosserat coupling modulus μ c [11,46,54]. The latter relates the relative rotation to the skew-symmetric part of the stress tensor and acts as a penalty parameter in order to limit the skew-symmetric Cosserat deformation.
The elastic moduli all have unit [Pa], as does the parameter f 0 . The elasticity parameters could be made to depend on the value of φ; however, this option is not explored in this work for mainly two reasons. First, single-phase materials are considered and the phase-field order parameter only distinguishes between grain and grain boundary regions. It is assumed that any difference in elastic behavior in the GB regions compared to the bulk of the grains can be ignored. Second, introducing such a dependence would also introduce additional driving forces for GB migration. This is outside of the scope of the current investigation which focuses on driving forces due to accumulated dislocation densities.
The term in brackets in Eq. (69) related to this latter parameter is nonstandard and is inspired by the KWC phase-field model [42,74]. A more classical formulation with separate contributions from the symmetric and skew-symmetric curvature tensor can be found in, e.g., [54]. The linear term in lattice curvature was shown by [42] to be necessary to localize the grain boundaries, whereas a higher-order term diffuses them and is necessary for GB mobility. In terms of a physical interpretation, the linear term in the norm of the curvature tensor can interpreted as the self-energy of the geometrically necessary dislocations, whereas the quadratic term represents the elastic interaction between them [53,56]. This interpretation relies on the relation between the torsion-curvature tensor and the dislocation density tensor, also called GND tensor. The linearized lattice curvature tensor is often used as an approximation of the full dislocation density tensor following [55]. Potentials involving the full dislocation density tensor were proposed in [19,29,56,75], and it was likewise used in the KWC inspired energy density in [4].
For GB migration in a single-phase material, the single-well potential 1 2 [ 1 − φ ] 2 is adopted in line with [42]. The coupling functions g(φ) and h(φ) must be nonnegative for ensuring positive grain boundary energies. For simplicity, we have considered increasing monotonous functions as in the original KWC model. It is worth noting that, with such a choice, it is still possible to fine tune g(φ) so as to recover the Read-Shockley dependency of the grain boundary energy with respect to the misorientation [1,42]. The parameters a, s and ε all have unit [m]. These parameters together with f 0 must be calibrated to yield the correct grain boundary energy and mobility for the GB migration. In the last term, χ is a constant parameter close to 0.3 [36].
From (69) together with Eqs. (49)- (52), the material behavior is now described by The inclusion of the linear term || Γ ∼ e || in the energy density results in a singular term in Eq. (71) which requires special attention in the numerical treatment [41,42,75].

Dissipation potentials for single crystal plasticity and GB migration
It is assumed for simplicity in this presentation that the wryness tensor is purely energetic, thus setting Γ ∼ p = 0 in Eq. (37). A format for the kinematics for plastic wryness due to edge and screw dislocations was proposed in [28]. A diffuse grain boundary description is adopted in this work, and it will be assumed that the evolution of plastic deformation in the bulk of the grains and inside the diffuse grain boundary take different forms. In the bulk of the grains, plastic deformation is assumed to take place through slip on preferred slip systems according to the classical description of single crystal plasticity. A viscoplastic flow rule taken from [17] is adopted. In order to allow for different behaviors inside the grain boundary, the plastic potential is assumed to have an additional term in the form: where the resolved shear stress τ α on slip system α is given by and α and n α are the slip direction and normal to the slip plane of system α. The second term depends only on the skew-symmetric part ofΠ ∼ M , represented by the pseudo-vector × Π M . In Eq. (74), K v and n are viscosity parameters. The thermodynamical forces R α represent the critical resolved shear stress on each slip system. The viscosity parameter η , that has unit [Pa s], is allowed to depend in a general way on the phase-field parameter φ and its gradient, as well as the curvature tensor Γ ∼ e , to distinguish the behavior in the bulk from the behavior in the grain boundary. Following [32], the thermodynamical variables are treated as parameters when they are included in the dissipation potential. An alternative to (74) which is not pursued here is to adopt a mixture-type format [5,59]. From Eq. (54) together with (74), the evolution for the plastic deformation is derived to beḞ where the slip rateγ α is given byγ where • = Max(•, 0). Two contributions are visible in Eq. (77): the usual plastic flow rule for single crystals, on the one hand, and a Newtonian contribution proportional to the skew-symmetric part of the Mandel stress tensor which is active in the GB only, on the other hand. The internal variables r α are related to the densities of statistically stored dislocations as follows. First, the associated forces R α are given by It is assumed that each r α depends on the statistically stored dislocation densities ρ β according to where b is the norm of the Burgers vector of the considered slip system family and h αβ is the interaction matrix accounting for self and latent hardening. The evolution of dislocations during plastic deformation is due to multiplication and annihilation mechanisms which can be described by the following Kocks-Mecking-Teodosiu evolution equation [43,69]: Here, a term has been added to account for the annihilation of dislocations due to grain boundary migration [1,7]. It is active ifφ > 0 in order to ensure that the recovery only takes place behind the sweeping grain boundary. The parameter C D determines the dynamics of the recovery with a lower value leading to only partial recovery and a higher value leading to full recovery. The function A(|| Γ ∼ e ||) serves to localize the process inside the grain boundary region. The other terms in (80) describe, respectively, the competing dislocation storage and dynamic recovery taking place during the deformation process. The parameter K is a dislocation mobility constant, and d is the critical annihilation distance between opposite sign dislocations. It may be noted that Eq. (80) does not imply transfer of dislocation densities between slip systems. Instead, it is assumed in the model that the lattice rotates to change its orientation due to GB migration and that each slip system is rotated in the same way keeping track of each associated dislocation density. This may not reflect the actual physical mechanism of dislocation sweeping by the GB but is a consequence of the proposed phase-field approach. In general, however, parameter C D is large enough for the dislocations to disappear after the sweeping of the GB [7]. Equation (80) together with (78) and (79) completely describe the kinematics for the internal variables, but they have not been derived from the dissipation potential Ω p (Π ∼ M , R α ) because such a potential for the Kocks-Mecking-Teodosiu model is not currently known [15].
Inserted into (47), the dissipation related to the internal variables becomes Recovery processes and annihilation of dislocations lead to dissipation of energy (sinceρ α < 0 ⇒ D α > 0), whereas energy is stored if the dislocation density increases. Lastly, the dissipation potential related to the phase-field variable is taken to be a simple quadratic format following [1,34] such that which clearly leads to nonnegative dissipation in (58). In (82), η φ is a viscosity parameter which is allowed to depend in a general way on the phase field and curvature. It has unit [Pa s]. The kinematics of the phase-field φ is then given by

Description using intermediate phase-field quantities
For the alternative formulation with the phase-field part formulated in the intermediate configuration, K φ should be used in the energy density rather than ∇φ and constitutive relation (66) replaces Eq. (52). The energy density is taken to be The intermediate microstress measure ξ p is then given by In this alternative formulation, the evolution of viscoplastic deformation is associated with an enhanced Mandeltype stress measure according to (65). The plastic potential in terms ofΠ ∼ M is taken to be The term x α = −J p α · K φ n α · ξ p plays the role of a back-stress for each slip system. A consequence of the introduction of the intermediate quantity K φ is therefore additional kinematic hardening, due to the higher-order terms. This fact has been discussed in a different context of strain gradient plasticity (cf. [47]).

Initial conditions on lattice orientation
In a polycrystal, the grains have different initial orientations even in an undeformed state. In this work, the Cosserat rotations are associated with the lattice directions and it is therefore considered that R ∼ 0 = R ∼ t=0 = 0 ∼ . In the absence of initial residual stresses, there should be zero stress in the undeformed polycrystal which requires that F ∼ e 0 = F ∼ e t=0 = I ∼ . This can be achieved by introducing nonvanishing initial values for the plastic distortion: so that With F ∼ e 0 = R ∼ T 0 · F ∼ e 0 , this implies that F ∼ e 0 = R ∼ 0 . Such a decomposition was also proposed in [3] for diffuse interface-based models for polycrystals. A piecewise constant field R ∼ 0 would correspond to a sharp interface description, whereas a diffuse interface description is obtained by assuming that R ∼ 0 varies smoothly over the grain boundaries. The latter approach will be adopted here.

Linearization and comparison with KWC model
In this section, the linearized model is derived for the case of small displacements, small rotations and small curvature. The resulting equations are then compared with the governing equations of the KWC phase-field model.

Kinematics
The linearized rotation tensor is given by so that the Cosserat deformation can be written approximately as and its rate as where ε ∼ is the usual small strain tensor and ω ∼ is the spin tensor taken as the skew-symmetric part of the velocity gradient. The linearized wryness tensor is given by The deformation measures in the geometrically linear theory are therefore taken to be together with the curvature tensor In the linearized setting, the elastic-viscoplastic multiplicative split of Cosserat deformation (35) becomes an additive decomposition according to No corresponding additive split of the curvature tensor is proposed here because it is assumed to be purely energetic in this work for simplicity. By defining which in turns implies that× This last relation says that the skew-symmetric part of the elastic deformation in the model is the difference between the lattice rotation rate and the Cosserat spin. This is an important result as it provides a link between the lattice and microrotation rates. The lattice and Cosserat rotations coincide if the elastic deformation is symmetric, i.e., if × e e ≡ 0.
The previous kinematic equations are the same as the ones postulated within the small deformation framework in [7].

Constitutive and evolution equations
Helmholtz energy density (69) for the linear case is given by where ε ∼ e is the symmetric part of the linear elastic deformation e ∼ e . The corresponding state laws are given by The evolution of plastic strain in the linearized case is given bẏ where the resolved shear stress on slip system α is τ α = α · σ ∼ · n α and the slip rateγ α is given by (77). Equation (106) can be split into symmetric and antisymmetric parts so that The last term of Eq. (108) will be called ω ∼ , and the evolution law can then be split into two distinct contributions The viscosity parameter η (φ, ∇φ, κ ∼ ) is chosen so that it is very large inside the bulk, with the consequence that × ω is only allowed to change significantly inside the grain boundary. In the interior of the grains, classical crystal plasticity by slip will still be the dominant mechanism for plastic evolution. Based on (109) and (110), the skew-symmetric plastic deformation is decomposed as with× The quantity × e , which was introduced as an eigendeformation in [7], has now been associated explicitly with the viscoplastic spin.

Initial conditions on lattice orientation
In the interest of modeling a polycrystal which is initially not elastically deformed (no residual stresses) but still consists of grains of different orientations, the orientation field at t = 0 is allowed to be nonzero at time it is possible to have an initially oriented configuration without displacements and without elastic deformation if It is assumed that the initial plastic deformation is entirely due to × e , so that

Comparison with KWC model
The proposed model can now be compared with the KWC model [42,74]. The KWC model is formulated in 2D for two fields, φ and the orientation θ . Evolution equations are found by the usual relaxational dynamics of a phase-field method. Assuming an energy density the resulting phase-field evolution equations for the KWC model are given by where η φ and η θ are mobility functions. They may in general depend on the phase fields and their gradients to distinguish the behavior in the bulk of the grains from the behavior inside the grain boundary region.

Inserting (104), (105) and (83) into (24) results in
In the same manner, inserting (102) and (103) into (23) results in The last expression can be rewritten as which together with (110) gives In particular, for a plane problem where rotation takes place only around one axis, there is only one nonzero component θ in Θ and the above expressions (120) and (121) can be written as where × e e is the only nonzero component of × e e . Alternatively, the last expression can be rewritten according to (123) as where × ω is the only nonzero component of × ω . By construction, the proposed coupled model does not change the kinematics for φ compared to the KWC model except for the added driving force due to statistically stored dislocations. This important addition, which was first introduced in [2], accounts for the migration of grain boundaries due to energy stored during plastic deformation. The present coupled model omits the relaxational dynamics for the orientation field, as evidenced by comparison between (119) and (121) where the η θ term is absent. This is also a new feature compared to the coupled approach for small deformations proposed in [7], where the skew-symmetric stress × σ was assumed to contain a dissipative contribution resulting in an evolution equation for the relative rotation × ω e −Θ. A new contribution compared to the KWC model is the coupling between lattice rotation and elastic rotation due to the Cosserat term proportional to μ c in (121). In the bulk of the grains, this term ensures that the lattice directors are properly reoriented to accommodate elastoplastic deformation. Reorientation of the lattice due to grain boundary migration on the other hand is a different process. In the proposed model, a migrating grain boundary gives rise to a (skew-symmetric) stress in the GB region. This stress is relaxed by the viscoplastic spin

Application to GB migration and viscoplastic deformation
In this section, numerical examples are presented to illustrate several important aspects and capabilities of the proposed model. It is shown that the relaxational dynamics for the orientation phase field in the KWC model can be replaced by the additional grain boundary term in plastic evolution law (54) together with the Cosserat penalty term. Both curvature driven and stored energy driven GB migration are considered. Lastly, the evolution of orientation in the grains due to plastic slip in the vicinity of moving GB is demonstrated.

Some elements of the numerical implementation
The linearized coupled phase-field Cosserat theory summarized in Sect. 4 has been implemented in the finite element code Z-set [76]. The implementation combines the approaches used in [27] for Cosserat crystal plasticity and in [1] for the KWC model, as described in [7]. An implicit iterative resolution scheme is used to solve the balance equations based on the Newton-Raphson algorithm. A fourth-order Runge-Kutta method  with automatic time stepping is applied to solve the differential equations driving the internal variables of the model. The finite element discretization uses quadratic elements with reduced integration. The Cosserat pseudo-vector Θ is associated with the lattice orientations with respect to a fixed reference frame, which is here taken to be (x 1 , x 2 , x 3 ) (see Figs. 2,8). With the calculations restricted to the (x 1 , x 2 )plane, the axis of rotation is the x 3 -axis. This gives Θ = [ 0 0 θ 3 ] T . In what follows, the subscript will be dropped when referring to the rotation angle so that θ = θ 3 . In the two-dimensional case, the monolithic weak formulation then involves 4 nodal degrees of freedom: 2 components of displacement, the Cosserat microrotation component θ and the phase-field φ. The angle of rotation θ is used to construct the appropriate orthogonal transformation tensors between the global and the local frame. This is different from a model using a sharp interface description since Θ in the present model varies as a phase-field variable over grain boundaries, i.e., it is continuous everywhere in the computational domain.
Initial conditions on the fields φ and Θ are found from solving the KWC model with the mobility function τ θ in Eq. (119) chosen such that it prevents GB migration and grain rotation (c.f. [7]). The plastic spin × e is also initialized making use of the KWC solution according to Eq. (116).

Coupling functions, mobility functions and model parameters
For the purpose of illustration, material parameters in the reasonable range for pure copper are used. The parameters related to the phase-field contribution were calibrated in [7] using the results found for the sharp interface limit of the KWC model in [48]. They are given in Table 1. The parameters for the elasto-viscoplastic behavior, given in Table 2, are chosen based on reported values in the literature [30,31] with the exception of the parametersη and C which are taken from [7], where it was demonstrated thatη < η φ should hold for mobile grain boundaries.
Calculations are performed on a non-dimensionalized system, where a suitable length scale Λ of unit [m] is introduced to obtain the dimensionless coordinates (x, y, z) = 1 Λ (x, y, z). The dimensionless differential operator ∇ is then given by ∇ = Λ ∇. Furthermore, κ ∼ = Θ ⊗ ∇ = Λ κ ∼ . In the two-dimensional case, there is only one nonzero component θ of Θ and ||κ ∼ || can be replaced by |∇θ | in the energy density. By division of (101) with f 0 , the dimensionless free energy function is then found to be where a = a/Λ, s = s/Λ and ε = ε/Λ, together with A regularization scheme proposed in [74] is adopted to treat the singularity resulting from the linear term |∇θ |, which is replaced by the function A γ (|∇θ |): where the constant γ is a large positive number. Simple formats are chosen for the coupling functions in the energy density: It is assumed that the viscosity-type parameter η φ in (120) is constant. The viscosity-type parameter η in (110) is assumed to depend on the position so that it is small inside the grain boundary and very large inside the bulk of the grain and is therefore taken as where C = C /Λ is dimensionless. Temperature dependence is omitted for the time being but can easily be modeled using an Arrhenius-type law as in [1] to reflect the mobility behavior of real grain boundaries. The viscosity parameters can be non-dimensionalized by choosing an appropriate timescale τ 0 so that t = t/τ 0 is dimensionless and The evolution of plasticity can be non-dimensionalized by . The parameters of the Kocks-Mecking-Teodosiu model, given in Table 3, are chosen based on values found in the literature for copper [8,25]. They are non-dimensionalized according to b = b/Λ and d = d/Λ, and the dimensionless dislocation density is given by ρ α = ρ α Λ 2 . The function A(|∇θ |) in the recovery term due to GB migration is taken to be where C A = C A /Λ is dimensionless. This term localizes the recovery to the grain boundary region. The influence of the (dimensionless) parameter C D on the recovery was studied in [1], where it was shown that full recovery takes place for sufficiently large values. The length scale parameter is chosen to be Λ = 1 μm. The parameter f 0 can be found from the expression whereγ gb is a dimensionless grain boundary energy calculated using the sharp interface limit of the phasefield model [48] and γ gb is the true grain boundary energy which must be measured experimentally. Here, it is assumed that at 15 • misorientation the grain boundary energy is 0.5 J/m 2 , which is in the reasonable range for copper (see, e.g., [71]). This gives f 0 = 8.2 MPa. A suitable timescale can be chosen based on experimental values of GB mobility according to where M is a dimensionless mobility calculated based on the results of the sharp interface limit of the KWC model (c.f. [7,48]).

Circular grain and grain rotation
A small circular grain is embedded in a larger grain according to Fig. 2 (left). Similar calculations were done for a circular grain using the pure KWC model in [42]. The misorientation between the two grains is 15 • ≈ 0.2618 rad. Dirichlet conditions are assumed on the boundaries of the computational domain. In the FE discretization, 100 by 100 elements with quadratic interpolation are used. According to the pure KWC model, the grain is able to both shrink (due to capillary forces) and rotate. At the space scale of interest in this work (micron size), bulk grain rotation is not expected to take place. Suppression of rotation in the KWC model can only be achieved by separating the timescales of evolution of orientation in the grain boundary and in the bulk of the grain. This is done by choosing appropriate formats for the mobility functions and parameters [42,74]. In the coupled model proposed here, it is assumed that there is no relaxational behavior related to the orientation and it is therefore not possible to control the rotation inside the grain in this manner. Instead, the evolution of crystal orientation in the bulk is controlled by the Cosserat coupling term. This term penalizes differences between elastic lattice rotation and Cosserat rotation. In the present case where no deformation is imposed, no rotation inside the grain is expected. Figure 3 shows the initial fields φ and θ (top and bottom, respectively) and the corresponding fields at t = 35. The Cosserat parameter in this simulation was chosen to be μ c = μ c / f 0 = 1000 (so that μ c = 8.2 GPa, to be compared with μ = 46.1 GPa). For this value of the Cosserat coupling parameter, it is clear from Fig. 3 that the desired behavior of shrinking without grain rotation is obtained. On the other hand, μ c = 1 leads to considerable grain rotation. This is demonstrated in Fig. 4 by comparing the progress of the simulation at t = 5, 10, 15, 20, 25, 30, 35 (left to right) for μ c = 1000 (top row) and μ c = 1 (bottom row).
For the lower value of the Cosserat modulus, the grain rotates as it shrinks and will eventually be completely aligned with the surrounding grain. Simulations were performed using several different values of μ c to study its influence. The results are shown in Fig. 5 for four different values of the penalty parameter μ c = μ c / f 0 : 1 (dotted line), 10 (dashed/red line), 100 (dashed-dotted/blue line) and 1000 (solid line). Figure 5 (left) shows the profiles of θ along a radial cut from the center of the grain along x 1 in Fig. 2 at time t = 35. For the lowest value of the Cosserat parameter, μ c = 1, it is clear that the grain in the center has rotated. This tendency is also seen for μ c = 10 although at a slower rate. For μ c = 100 or μ c = 1000, there seems to be little or no rotation.  which is why it seems as if the grain boundary position tends toward a fixed value. In reality, the grain has nearly vanished and the difference in orientation is so small that the grain boundary migration has stopped. In Helmholtz free energy (101), the Cosserat coupling term is in competition with the terms due to the phase-field and grain boundary contribution (including the curvature terms). By rotating the grains to lower the misorientation, the energy due to curvature is also lowered. Inside the bulk however, grain rotation increases the energy due to the Cosserat term. Specifically, in the 2D case (127), the two (non-dimensionalized) contributions are given by In the left-hand side Fig. 6, the evolution over time of the total (dimensionless) energy F due to these densities (integrated over the entire computational domain) is shown, together with the separate contributions from each term for μ c = μ c / f 0 = 1, 10, 1000 (top to bottom). The separate contributions at times t = 0 (black lines) and t = 35 (red lines) taken along a radial cut as in Fig. 5 are shown to the right in Fig. 6. For μ c = 1, it is clear that the gain from lowering the misorientation is larger than the cost due to the Cosserat term ψ C , and therefore, the circular grain rotates which lowers ψ P F but increases ψ C . For μ c = 1000 on the other hand, the misorientation remains constant which indicates that the penalization factor is high enough to prevent rotation. For the case μ c = 10, there is some slight grain rotation toward the very end of the simulation which is evidenced by the fact that the contribution to the energy from the Cosserat term increases slightly. The profiles of φ (left) and θ (right) are plotted in Fig. 7 to illustrate the progress of the simulations. Note that in the absence of displacements and plastic deformations due to slip processes, This example shows that the proposed coupled model captures GB migration due to capillary forces. The influence of the Cosserat coupling parameter μ c was investigated in particular, showing that if this parameter is chosen large enough, there is no net energy gain from lowering the misorientation and no rotation takes place in the bulk of the grain.

Grain boundary migration due to SSDs
The second example is a periodic laminate structure. The computational domain consists of two grains of different orientations which are assumed to be periodically repeated according to Fig. 8. In this case, the grain boundary is flat, and unlike the previous example, there is no GB migration due to curvature. Instead, a singleslip system is assumed and an initial dislocation content of ρ = 10 15 m −2 [37] is assigned to the 15 • oriented grain. For purposes of illustration, the other grain is assumed to be dislocation free. There is no deformation or stress so that no plastic slip is induced. The entire evolution of SSD content is therefore due to the static recovery term in Eq. (80). The stored energy density is given by which in dimensionless quantities together with the material parameters in Tables 1, 2  (139) Figure 9 shows the profile of θ initially and at t = 300. The grain boundary has migrated so that the dislocation free grain expands into the grain with a high stored energy. Since μ c = μ c / f 0 = 1000 in this simulation, there is no bulk grain rotation. Figure 10 shows the profiles of θ (top left) and φ (bottom left) initially and at t = 0, 150, 300 over the left GB. In the 15 • grain, the phase field takes values φ < 1 in the bulk due to the presence of dislocations which changes the solution to the static equilibrium of evolution Eq. (120) (see [1] for the calculation of this equilibrium value). There is complete static recovery behind the sweeping grain boundary as evidenced in Fig. 10 (top right). Finally, the position of the grain boundary over time is shown in Fig. 10 (bottom right). As expected, the velocity is constant. 5 The velocity of the grain boundary is largely independent of the parameter μ c when it is chosen sufficiently large, as evidenced by the fact that the calculations using μ c = 1000 (black line) and μ c = 10 4 (dashed/red line) give the same result. 5 If the true grain boundary velocity v at a given temperature and driving force (corresponding to the energy due to the given dislocation content) is known, then the characteristic time τ 0 = t/t can be calculated as τ 0 = Λv/v, where v is the slope in Fig. 10 (bottom right). Due to the periodic boundary conditions applied in this example, the KWC model would predict bulk grain rotation in addition to grain boundary migration. Such grain rotation in the KWC model can be prevented by choosing the mobility function η θ (φ, ∇φ, ∇θ) in (119) so that it takes much higher values in the bulk of the grain than in the grain boundary [7,74]. To illustrate this, the equilibrium solution to the KWC Eqs. (118) and (119) is calculated using a mobility function [7] on the form with the constants chosen asη θ = 1000 f 0 τ 0 , μ P = 10 5 Λ and β P = 10 4 where τ 0 and Λ are the characteristic time and length scales introduced to obtain a dimensionless system. The remaining parameters are given in Table 1. The shape of the dimensionless mobility functionη θ /( f 0 τ 0 ) is shown in Fig. 11 (right). The driving force for evolution of the orientation is plotted in Fig. 11 (left). In order to be consistent with the formalism of the coupled model, the notation of Eq. (122) is adopted and the driving force (or rather force stress) is formally written as where ξ θ is a couple stress which only has one nonzero component ξ θ in the simulations. Evidently, inside the grain boundary where the mobility is small, the stress is zero as it has relaxed. Outside the grain boundary region, there is still a nonzero stress. In this region, evolution of the orientation is prevented due to the very large value of the mobility function, and therefore, the solution is stable also here on the time scale of the simulation. In the coupled model, there is not a separate evolution law for Θ. Recalling Eqs. (121) and (122) and making use of the fact that in the two-dimensional case κ ∼ can be replaced by ∇θ , the skew-symmetric stress is given by where × e e is the only nonzero component of the skew-symmetric elastic deformation × e e . There are two competing terms. On the one hand, the Cosserat term with the modulus μ c penalizes relative rotations inside the bulk. On the other hand, there is the same term ∇ · ξ θ that produced bulk rotation in the KWC model (ξ θ corresponds to the nonzero contributions to m ∼ in the two-dimensional case). The equilibrium solutions with a static grain boundary are calculated for the coupled model using different values of the Cosserat penalty parameter. The skew-symmetric stress and the couple stress term are shown in Fig. 12 for μ c = 100 (left) and μ c = 1000 (right). For the lower value of the coupling modulus, the skew-symmetric stress contribution is nonzero everywhere, indicating that the grains have rotated.
Finally, the evolution of the skew-symmetric stress contribution is traced in the case of a migrating grain boundary. The result is shown in Fig. 13 for μ c = 1000 (left) or μ c = 10 4 (right). The behavior is similar although for the higher value of the coupling modulus there is a smaller region of nonzero stress around the grain boundary than for the lower value, which is expected. The skew-symmetric stress increases first due to lattice reorientation in the moving GB. It is then relaxed by the viscoplastic spin × ω according to Eq. (110). Inside the grain boundary, the related viscosity parameter η (φ, ∇φ, κ ∼ ) is small and there is complete relaxation of the stress, as evidenced in Fig. 13. Outside of the grain boundary region, this parameter is very large and the skew-symmetric stress is governed only by the Cosserat coupling term. This is why the stress to the left of the initial position of the grain boundary vanishes very slowly as it is outside of the grain boundary region.
This example demonstrates that the model is capable of predicting GB migration due to energy stored in the form of statistically stored dislocations and also the static recovery behind the sweeping grain boundary. The same behavior as in [7] is recovered which further demonstrates that it is not necessary to include a separate relaxational behavior for the relative rotation. This example also highlights the essential role played by the skew-symmetric part of the stress which represents the material resisting force to lattice reorientation. The latter impedes relative material vs. lattice rotation in the bulk and is relaxed in the GB region to allow for further GB motion.

Plastic deformation in simple shear
Periodic simple shear loading is applied to the same geometry as in the previous example according to 0.08 0.18 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 Fig. 14 The orientation field in the sheared bigrain structure at t = 10 (top) when the total deformation has been applied and at t = 10,000 (bottom), where it can be seen that the difference in stored energy between the two grains has resulted in GB migration where r is the periodic displacement fluctuation taking the same value at homologous points of the boundary and the components of the mean strain tensor E ∼ are given by The simple shear contains overall straining and rotation imposed to the unit cell. The deformation is ramped up from t = 0 to t = 10 and then held constant until t = 10 4 (note that the timescale is dimensionless). For illustrative purposes, only one slip system is supposed to be activated, with the slip direction and slip plane normal given by = (0, 1, 0), n = (1, 0, 0).
For this choice, slip is favored in the 0 • grain region (where the lattice is aligned with the slip direction) which results in a net difference of statistically stored dislocations between the two grains. This difference in SSD content drives GB migration as evidenced in Fig. 14. The orientation field is shown at t = 10 (top) and t = 10,000 (bottom). The activation of the slip system in the second grain is due to stress increase induced by SSD hardening. An initial dislocation density of ρ 0 = 10 11 m −2 is assumed in both grains. Due to plastic slip, the dislocation content in the 0 • grain increases to ρ max = 2.5 × 10 13 m −2 , whereas it only goes to ρ min = 5.7 × 10 12 m −2 in the 15 • grain. This can be seen in Fig. 15 (bottom right) which shows the dislocation density across the left grain boundary of the simulation domain at t = 10 (black line) and t = 10,000 (dashed/red line). It can also be seen that the migrating grain boundary has resulted in some, although not complete, static recovery. Figure 15 (top left) shows the profile of the orientation θ at t = 0 (dotted line) and t = 10 (black line) together with the elastic rotation due to the shear deformation (dashed/red line). The evolution of the orientation due to the grain boundary migration is shown in Fig. 15 (bottom left). The evolution of the orientation inside the bulk of the grains from t = 0 to t = 10 is due to plastic slip processes and associated elastic rotation. Due to the relative difference in SSD content in the two grains, the GB then eventually migrates. Figure 15 (top right) shows the stress-strain curve for the shear stress and shear strain in the local coordinate system of the grain at initially 0 • orientation. This is the grain that carries most of the deformation as it is aligned with the slip direction. This example clearly demonstrates that the model is capable of predicting, in a coupled manner, grain reorientation due to both plastic slip processes, which change the overall bulk orientation of the grains, and GB migration due to statistically stored dislocations which accumulate during the plastic deformation.

Conclusions
The original contributions of the present work are the following: 1. The KWC GB migration model and crystal plasticity theory have been coupled within a thermodynamically consistent framework including finite deformation, rotation and curvature. The previous theory proposed in [7] has been shown to result from consistent linearization of the present theory. 2. The Cosserat framework allows to get rid of the orientation relaxation term of the original KWC which is thought to be redundant with the phase-field relaxation contribution. The skew-symmetric stress components prevents bulk rotation and is relaxed in the GB region thus permitting further GB motion. 3. A first example illustrating non-homogeneous plastic deformation in a laminate microstructure and subsequent GB migration has been provided.
The large deformation Cosserat crystal plasticity model extended with a phase-field order parameter allows modeling of concurrent viscoplastic deformation and grain boundary migration. The energy function describing the single and polycrystal material is assumed to contain contributions of lattice curvature inspired by orientation phase-field models and given physical significance through the association between the lattice curvature tensor and the density of geometrically necessary dislocations. The internal thermodynamic variables of the plasticity model are associated with the statistically stored dislocations on each slip system. Their evolution is assumed to follow a classical Kocks-Mecking-Teodosiu law extended with a term for static recovery due to GB migration. Using a linearized version of the general large deformation framework, it was demonstrated by numerical examples that the proposed framework can predict GB migration due to GNDs as well as SSDs.
The large deformation framework proposed in this work can be considered a full 3D generalization of the KWC model to a coupled phase-field-crystal plasticity approach. The general theory is not restricted to the constitutive choices made and can be extended to consider, e.g., elastic anisotropy as well as anisotropic properties of the grain boundaries. Linearization of the model essentially recovers the formulation proposed in [7] with some modifications. Firstly, the eigenrotation introduced in [7] to relax skew-symmetric stress in the grain boundaries has been fully integrated with the plastic deformation flow rule including plastic spin. Secondly, it was shown that it is not necessary to include a separate kinematics for the evolution of relative rotation. Instead, the skew-symmetric part of the Mandel-type stress tensor is the driving force for material reorientation inside grain boundaries, whereas lattice rotation inside the bulk of the grains is driven by plastic slip processes due to deformation.
As a next step, the model will be applied to realistic examples such as the experimental results given in [10] for GB migration in metallic polycrystals undergoing thermomechanical loading. In order to model processes such as cold rolling, the numerical implementation of the full large deformation model is necessary. An extension of the numerical implementation to 3D requires some care regarding the interpolation of orientation across the diffuse grain boundary regions [3,58]. In addition to modeling examples of deformation and GB migration in larger polycrystalline samples, the model can be used to study the behavior under loading of individual grain boundaries with respect to, e.g., dislocation pileup. being the sum of the Mandel-type stress tensors defined as Noting that it appears that the second contribution to (152) cancels out. Thus, the formulation in the local configuration presented above yields formally the same Mandel-type stress as for (41) together with (49) of Sect. 2.3 where the material phase-field quantities were used.